11 #include "GaudiKernel/IPartPropSvc.h"
12 #include "GaudiKernel/PhysicalConstants.h"
28 #include "G4PrimaryParticle.hh"
30 #include "G4Geantino.hh"
31 #include "G4ChargedGeantino.hh"
32 #include "G4ParticleTable.hh"
33 #include "G4LorentzVector.hh"
34 #include "G4TransportationManager.hh"
40 #include "CLHEP/Geometry/Point3D.h"
41 #include "CLHEP/Geometry/Vector3D.h"
42 #include "CLHEP/Units/SystemOfUnits.h"
44 #include "HepPDT/ParticleID.hh"
45 #include "HepPDT/DecayData.hh"
46 #include "HepPDT/ParticleDataTable.hh"
52 , m_particlePropSvc(
"PartPropSvc",
name)
53 , m_particleDataTable(nullptr)
54 , m_useGeneratedParticleMass(false)
55 , m_genParticleFilters(this)
56 , m_quasiStableParticlesIncluded(false)
61 "Use particle mass assigned to GenParticle.");
65 "Tools for filtering out GenParticles.");
69 "ParticlePropertyService to retrieve the PDT.");
87 if (!m_useGeneratedParticleMass) {
89 m_particleDataTable = m_particlePropSvc->PDT();
90 if (!m_particleDataTable) {
91 ATH_MSG_FATAL(
"Could not get ParticleDataTable from " << m_particlePropSvc <<
". Abort" );
92 return StatusCode::FAILURE;
96 if (!m_genParticleFilters.empty()) {
97 ATH_CHECK(m_genParticleFilters.retrieve());
101 return StatusCode::SUCCESS;
110 return StatusCode::SUCCESS;
120 for (
auto eventPtr : inputGenEvents ) {
122 if (eventPtr ==
nullptr) {
continue; }
126 " and event_number=" << eventPtr->event_number() );
129 bool legacyOrdering =
true;
136 for (
auto& genPartPtr : passedGenParticles ) {
137 ATH_MSG_VERBOSE(
"Picking up following GenParticle for conversion to ISFParticle: " << genPartPtr);
138 auto simParticlePtr = convertParticle(genPartPtr);
139 if (!simParticlePtr) {
140 ATH_MSG_ERROR(
"Error while trying to convert input generator particles. Aborting.");
141 return StatusCode::FAILURE;
144 simParticles.push_back(simParticlePtr);
150 ATH_MSG_DEBUG(
"Created initial simulation particle collection with size " << simParticles.size() );
152 return StatusCode::SUCCESS;
162 std::make_move_iterator(
std::begin(simParticleList)),
163 std::make_move_iterator(
std::end(simParticleList))
165 if (!shadowGenEvents.
empty()) {
166 outputG4Event = this->ISF_to_G4Event(simParticleVector, inputGenEvents.
back(), shadowGenEvents.
back());
169 outputG4Event = this->ISF_to_G4Event(simParticleVector, inputGenEvents.
back(),
nullptr);
171 return StatusCode::SUCCESS;
176 G4Event*& outputG4Event)
const
182 std::make_move_iterator(
std::begin(simParticleList)),
183 std::make_move_iterator(
std::end(simParticleList))
185 outputG4Event = this->ISF_to_G4Event(simParticleVector, inputGenEvents.
back(),
nullptr);
186 return StatusCode::SUCCESS;
192 std::vector<HepMC::GenParticlePtr>
194 auto allGenPartBegin = evnt.particles().begin();
195 auto allGenPartEnd = evnt.particles().end();
198 std::vector<HepMC::GenParticlePtr> passedGenParticles{};
202 if (m_useShadowEvent) {
203 if (legacyOrdering) {
206 for (
auto vtx: evnt.vertices() ) {
208 vtx->particles_out().end(),
209 std::back_inserter(passedGenParticles),
216 std::back_inserter(passedGenParticles),
221 if (legacyOrdering) {
224 for (
auto vtx: evnt.vertices() ) {
226 vtx->particles_out().end(),
227 std::back_inserter(passedGenParticles),
228 [
this](
HepMC::GenParticlePtr p){return this->passesFilters(std::const_pointer_cast<const HepMC3::GenParticle>(p));});
234 std::back_inserter(passedGenParticles),
235 [
this](
HepMC::GenParticlePtr p){
return this->passesFilters(std::const_pointer_cast<const HepMC3::GenParticle>(
p));});
239 passedGenParticles.shrink_to_fit();
241 return passedGenParticles;
244 std::vector<HepMC::GenParticlePtr>
246 auto allGenPartBegin = evnt.particles_begin();
247 auto allGenPartEnd = evnt.particles_end();
250 std::vector<HepMC::GenParticlePtr> passedGenParticles{};
254 if (legacyOrdering) {
257 auto vtxIt = evnt.vertices_begin();
258 auto vtxItEnd = evnt.vertices_end();
259 for ( ; vtxIt != vtxItEnd; ++vtxIt ) {
260 const auto vtxPtr = *vtxIt;
263 std::back_inserter(passedGenParticles),
270 std::back_inserter(passedGenParticles),
274 passedGenParticles.shrink_to_fit();
276 return passedGenParticles;
284 if (!genPartPtr) {
return nullptr; }
286 auto pVertex = genPartPtr->production_vertex();
288 ATH_MSG_ERROR(
"Unable to convert following generator particle due to missing production vertex for: " << genPartPtr);
291 auto parentEvent = genPartPtr->parent_event();
293 ATH_MSG_ERROR(
"Cannot convert a GenParticle without a parent GenEvent into an ISFParticle!!!");
297 const Amg::Vector3D pos(pVertex->position().x(), pVertex->position().y(), pVertex->position().z());
298 const auto& pMomentum(genPartPtr->momentum());
301 const double pMass = this->getParticleMass(genPartPtr);
303 const double pMass = this->getParticleMass(*genPartPtr);
305 double e=pMomentum.e();
307 double px=pMomentum.px();
308 double py=pMomentum.py();
309 double pz=pMomentum.pz();
311 if (std::abs(
e-teste)>0.01*
e) {
312 ATH_MSG_WARNING(
"Difference in energy for: " << genPartPtr<<
" Morg="<<pMomentum.m()<<
" Mmod="<<
pMass<<
" Eorg="<<
e<<
" Emod="<<teste);
314 if (
MC::isDecayed(genPartPtr) && pVertex && genPartPtr->end_vertex()) {
315 const auto& prodVtx = genPartPtr->production_vertex()->position();
316 const auto& endVtx = genPartPtr->end_vertex()->position();
317 CLHEP::Hep3Vector dist3D(endVtx.x()-prodVtx.x(), endVtx.y()-prodVtx.y(), endVtx.z()-prodVtx.z());
320 CLHEP::HepLorentzVector
mom( pMomentum.x(), pMomentum.y(), pMomentum.z(), pMomentum.t() );
321 double gamma_org=
mom.gamma();
323 double gamma_new=
mom.gamma();
325 if (std::abs(gamma_new-gamma_org)/(gamma_new+gamma_org)>0.001) {
326 ATH_MSG_WARNING(
"Difference in boost gamma for Quasi stable particle "<<genPartPtr);
336 const int pPdgId = genPartPtr->pdg_id();
343 auto tBinding = std::make_unique<ISF::TruthBinding>(genPartPtr);
347 auto sParticle = std::make_unique<ISF::ISFParticle>( std::move(
pos),
352 genPartPtr->status(),
359 return sParticle.release();
368 double mass =
part->generated_mass();
372 if ( !m_useGeneratedParticleMass ) {
373 const int absPDG = std::abs(
part->pdg_id());
374 HepPDT::ParticleData
const *pData = (m_particleDataTable)
375 ? m_particleDataTable->particle(absPDG)
378 mass = pData->mass();
382 ATH_MSG_WARNING(
"Unable to find mass of particle with PDG ID '" << absPDG <<
"' in ParticleDataTable. Will set mass to generated_mass: " <<
mass);
392 double mass =
part.generated_mass();
396 if ( !m_useGeneratedParticleMass ) {
397 const int absPDG = std::abs(
part.pdg_id());
398 HepPDT::ParticleData
const *pData = (m_particleDataTable)
399 ? m_particleDataTable->particle(absPDG)
402 mass = pData->mass();
406 ATH_MSG_WARNING(
"Unable to find mass of particle with PDG ID '" << absPDG <<
"' in ParticleDataTable. Will set mass to generated_mass: " <<
mass);
420 for (
const auto&
filter : m_genParticleFilters ) {
425 :
"false, will remove particle."));
445 for (
const auto&
filter : m_genParticleFilters ) {
450 :
"false, will remove particle."));
471 const int eventID(1);
472 G4Event *g4evt =
new G4Event(eventID);
475 const G4VSolid *worldSolid = G4TransportationManager::GetTransportationManager()->GetNavigatorForTracking()->GetWorldVolume()->GetLogicalVolume()->GetSolid();
479 if ( !isInsideG4WorldVolume(isp, worldSolid) ) {
480 ATH_MSG_WARNING(
"Unable to convert ISFParticle to G4PrimaryParticle!");
484 worldSolid->DumpInfo();
492 this->addG4PrimaryVertex(g4evt,isp,useHepMC,shadowGenEvent);
499 g4evt->SetUserInformation(atlasG4EvtUserInfo);
509 return G4ChargedGeantino::Definition();
512 return G4Geantino::GeantinoDefinition();
515 G4ParticleTable *ptable = G4ParticleTable::GetParticleTable();
517 return ptable->FindParticle(pdgcode);
519 ATH_MSG_ERROR(
"getG4ParticleDefinition - Failed to retrieve G4ParticleTable!");
529 const double p2=
std::pow(g4particle.GetTotalMomentum(),2);
530 const double m2=
std::pow(g4particle.GetMass(),2);
531 const double l2=
std::pow(GeneratorDecayLength,2);
533 const double tau=std::sqrt(tau2);
534 g4particle.SetProperTime( tau );
543 const G4ParticleDefinition *particleDefinition = this->getG4ParticleDefinition(
genpart->pdg_id());
545 if (particleDefinition==
nullptr) {
546 ATH_MSG_ERROR(
"ISF_to_G4Event particle conversion failed. ISF_Particle PDG code = " <<
genpart->pdg_id() <<
547 "\n This usually indicates a problem with the evgen step.\n" <<
548 "Please report this to the Generators group, mentioning the release and generator used for evgen and the PDG code above." );
554 auto &genpartMomentum =
genpart->momentum();
555 G4double
px = genpartMomentum.x();
556 G4double
py = genpartMomentum.y();
557 G4double
pz = genpartMomentum.z();
559 std::unique_ptr<G4PrimaryParticle> g4particle = std::make_unique<G4PrimaryParticle>(particleDefinition,
px,
py,
pz);
565 const auto& prodVtx =
genpart->production_vertex()->position();
566 const auto& endVtx =
genpart->end_vertex()->position();
572 CLHEP::Hep3Vector dist3D(endVtx.x()-prodVtx.x(), endVtx.y()-prodVtx.y(), endVtx.z()-prodVtx.z());
576 double pmag2=g4particle->GetTotalMomentum();
578 double e2=g4particle->GetTotalEnergy();
580 double beta2=pmag2/
e2;
584 if (m_quasiStableParticlesIncluded) {
597 for (
auto daughter:
genpart->end_vertex()->particles_out() ) {
598 if (m_quasiStableParticlesIncluded) {
599 ATH_MSG_VERBOSE (
"Attempting to add daughter particle : " << daughter );
602 ATH_MSG_WARNING (
"Attempting to add daughter particle: " << daughter );
604 G4PrimaryParticle *daughterG4Particle = this->getDaughterG4PrimaryParticle( daughter, makeLinkToTruth );
605 if (!daughterG4Particle) {
606 ATH_MSG_ERROR(
"Bailing out of loop over daughters of particle: "<<
607 " due to errors - will not return G4Particle.");
610 g4particle->SetDaughter( daughterG4Particle );
614 if (makeLinkToTruth) {
616 std::unique_ptr<PrimaryParticleInformation> primaryPartInfo = std::make_unique<PrimaryParticleInformation>(
genpart);
618 g4particle->SetUserInformation(primaryPartInfo.release());
622 return g4particle.release();
629 const G4ParticleDefinition *particleDefinition = this->getG4ParticleDefinition(
genpart->pdg_id());
631 if (particleDefinition==
nullptr) {
632 ATH_MSG_ERROR(
"ISF_to_G4Event particle conversion failed. ISF_Particle PDG code = " <<
genpart->pdg_id() <<
633 "\n This usually indicates a problem with the evgen step.\n" <<
634 "Please report this to the Generators group, mentioning the release and generator used for evgen and the PDG code above." );
640 auto &genpartMomentum =
genpart->momentum();
641 G4double
px = genpartMomentum.x();
642 G4double
py = genpartMomentum.y();
643 G4double
pz = genpartMomentum.z();
645 std::unique_ptr<G4PrimaryParticle> g4particle = std::make_unique<G4PrimaryParticle>(particleDefinition,
px,
py,
pz);
651 const auto& prodVtx =
genpart->production_vertex()->position();
652 const auto& endVtx =
genpart->end_vertex()->position();
658 CLHEP::Hep3Vector dist3D(endVtx.x()-prodVtx.x(), endVtx.y()-prodVtx.y(), endVtx.z()-prodVtx.z());
662 double pmag2=g4particle->GetTotalMomentum();
664 double e2=g4particle->GetTotalEnergy();
666 double beta2=pmag2/
e2;
670 if (m_quasiStableParticlesIncluded) {
683 for (
auto daughter:
genpart->end_vertex()->particles_out() ) {
684 if (m_quasiStableParticlesIncluded) {
685 ATH_MSG_VERBOSE (
"Attempting to add daughter particle: " << daughter );
688 ATH_MSG_WARNING (
"Attempting to add daughter particle: " << daughter );
690 G4PrimaryParticle *daughterG4Particle = this->getDaughterG4PrimaryParticle( daughter );
691 if (!daughterG4Particle) {
692 ATH_MSG_ERROR(
"Bailing out of loop over daughters of particle due to errors - will not return G4Particle.");
695 g4particle->SetDaughter( daughterG4Particle );
699 return g4particle.release();
706 const G4ParticleDefinition *particleDefinition = this->getG4ParticleDefinition(
genpart.pdg_id());
708 if (particleDefinition==
nullptr) {
709 ATH_MSG_ERROR(
"ISF_to_G4Event particle conversion failed. ISF_Particle PDG code = " <<
genpart.pdg_id() <<
710 "\n This usually indicates a problem with the evgen step.\n" <<
711 "Please report this to the Generators group, mentioning the release and generator used for evgen and the PDG code above." );
717 auto &genpartMomentum =
genpart.momentum();
718 G4double
px = genpartMomentum.x();
719 G4double
py = genpartMomentum.y();
720 G4double
pz = genpartMomentum.z();
722 std::unique_ptr<G4PrimaryParticle> g4particle = std::make_unique<G4PrimaryParticle>(particleDefinition,
px,
py,
pz);
728 const auto& prodVtx =
genpart.production_vertex()->position();
729 const auto& endVtx =
genpart.end_vertex()->position();
735 CLHEP::Hep3Vector dist3D(endVtx.x()-prodVtx.x(), endVtx.y()-prodVtx.y(), endVtx.z()-prodVtx.z());
739 double pmag2=g4particle->GetTotalMomentum();
741 double e2=g4particle->GetTotalEnergy();
743 double beta2=pmag2/
e2;
747 if (m_quasiStableParticlesIncluded) {
760 for (
auto daughterIter=
genpart.end_vertex()->particles_out_const_begin();
761 daughterIter!=
genpart.end_vertex()->particles_out_const_end(); ++daughterIter ) {
762 if (m_quasiStableParticlesIncluded) {
763 ATH_MSG_VERBOSE (
"Attempting to add daughter particle: " << **daughterIter );
766 ATH_MSG_WARNING (
"Attempting to add daughter particle: " << **daughterIter );
768 G4PrimaryParticle *daughterG4Particle = this->getDaughterG4PrimaryParticle( **daughterIter, makeLinkToTruth );
769 if (!daughterG4Particle) {
770 ATH_MSG_ERROR(
"Bailing out of loop over daughters of particle due to errors - will not return G4Particle.");
773 g4particle->SetDaughter( daughterG4Particle );
777 if (makeLinkToTruth) {
779 std::unique_ptr<PrimaryParticleInformation> primaryPartInfo = std::make_unique<PrimaryParticleInformation>(&
genpart);
781 g4particle->SetUserInformation(primaryPartInfo.release());
785 return g4particle.release();
794 && (p1->status() == p2->status())
795 && (p1->pdg_id() == p2->pdg_id())
796 && ((p1->momentum().px()) == (p2->momentum().px()))
797 && ((p1->momentum().py()) == (p2->momentum().py()))
798 && ((p1->momentum().pz()) == (p2->momentum().pz()))
799 && (
float(p1->momentum().m()) ==
float(p2->momentum().m()));
805 if (!shadowGenEvent) {
806 ATH_MSG_FATAL (
"Found status==2 GenParticle with no end vertex and shadow GenEvent is missing - something is wrong here!");
811 const int shadowId = genParticle->attribute<HepMC3::IntAttribute>(
"ShadowParticleId")->
value();
812 for (
auto& shadowParticle : shadowGenEvent->particles()) {
813 if (shadowParticle->id() == shadowId && matchedGenParticles(genParticle, shadowParticle) ) {
return shadowParticle; }
815 return std::make_shared<HepMC::GenParticle>();
817 for (HepMC::GenEvent::particle_iterator pitr=shadowGenEvent->particles_begin(); pitr != shadowGenEvent->particles_end(); ++pitr) {
819 if (matchedGenParticles(genParticle, shadowParticle) ) {
return shadowParticle; }
833 const auto& prodVtx =
genpart->production_vertex()->position();
834 const auto& endVtx =
genpart->end_vertex()->position();
836 CLHEP::Hep3Vector dist3D(endVtx.x()-prodVtx.x(), endVtx.y()-prodVtx.y(), endVtx.z()-prodVtx.z());
840 double pmag2=g4particle->GetTotalMomentum();
842 double e2=g4particle->GetTotalEnergy();
844 double beta2=pmag2/
e2;
845 double mass2=g4particle->GetMass();
850 const G4LorentzVector lv0( prodVtx.x(), prodVtx.y(), prodVtx.z(), prodVtx.t() );
851 const G4LorentzVector lv1( endVtx.x(), endVtx.y(), endVtx.z(), endVtx.t() );
854 G4LorentzVector dist4D(lv1);
857 double dist4Dgamma=std::numeric_limits<double>::infinity();
858 if (dist4D.t()>0 && dist4D.mag2()>0) {
859 dist4Dgamma=dist4D.gamma();
864 G4LorentzVector fourmom(g4particle->GetMomentum(),g4particle->GetTotalEnergy());
865 double fourmomgamma=std::numeric_limits<double>::infinity();
866 if (fourmom.t()>0 && fourmom.mag2()>0) {
867 fourmomgamma=fourmom.gamma();
872 ATH_MSG_VERBOSE(
"gammaVertex="<<dist4Dgamma<<
" gammamom="<<fourmomgamma<<
" gamma(beta)="<<1/std::sqrt(1-beta2)<<
" lifetime tau(beta)="<<std::sqrt(tau2)<<
" lifetime tau="<<tau);
874 if (m_quasiStableParticlesIncluded) {
882 ATH_MSG_WARNING(
"Detected primary particle with end vertex. This should only be the case if" );
883 ATH_MSG_WARNING(
"you are running with quasi-stable particle simulation enabled. This is not" );
884 ATH_MSG_WARNING(
"yet validated - you'd better know what you're doing. Will add the primary" );
891 for (
auto daughter: *(
genpart->end_vertex())) {
892 if (m_quasiStableParticlesIncluded) {
896 ATH_MSG_WARNING (
"Attempting to add daughter particle: " << daughter );
898 G4PrimaryParticle *daughterG4Particle = this->getDaughterG4PrimaryParticle( daughter );
899 if (!daughterG4Particle) {
900 ATH_MSG_FATAL(
"Bailing out of loop over daughters due to errors.");
902 g4particle->SetDaughter( daughterG4Particle );
913 const auto& prodVtx =
genpart->production_vertex()->position();
914 const auto& endVtx =
genpart->end_vertex()->position();
916 CLHEP::Hep3Vector dist3D(endVtx.x()-prodVtx.x(), endVtx.y()-prodVtx.y(), endVtx.z()-prodVtx.z());
920 double pmag2=g4particle->GetTotalMomentum();
922 double e2=g4particle->GetTotalEnergy();
924 double beta2=pmag2/
e2;
925 double mass2=g4particle->GetMass();
930 const G4LorentzVector lv0( prodVtx.x(), prodVtx.y(), prodVtx.z(), prodVtx.t() );
931 const G4LorentzVector lv1( endVtx.x(), endVtx.y(), endVtx.z(), endVtx.t() );
934 G4LorentzVector dist4D(lv1);
937 double dist4Dgamma=std::numeric_limits<double>::infinity();
938 if (dist4D.t()>0 && dist4D.mag2()>0) {
939 dist4Dgamma=dist4D.gamma();
944 G4LorentzVector fourmom(g4particle->GetMomentum(),g4particle->GetTotalEnergy());
945 double fourmomgamma=std::numeric_limits<double>::infinity();
946 if (fourmom.t()>0 && fourmom.mag2()>0) {
947 fourmomgamma=fourmom.gamma();
952 ATH_MSG_VERBOSE(
"gammaVertex="<<dist4Dgamma<<
" gammamom="<<fourmomgamma<<
" gamma(beta)="<<1/std::sqrt(1-beta2)<<
" lifetime tau(beta)="<<std::sqrt(tau2)<<
" lifetime tau="<<tau);
954 if (m_quasiStableParticlesIncluded) {
962 ATH_MSG_WARNING(
"Detected primary particle with end vertex. This should only be the case if" );
963 ATH_MSG_WARNING(
"you are running with quasi-stable particle simulation enabled. This is not" );
964 ATH_MSG_WARNING(
"yet validated - you'd better know what you're doing. Will add the primary" );
971 for (
auto daughter: *(
genpart->end_vertex())) {
972 if (m_quasiStableParticlesIncluded) {
973 ATH_MSG_VERBOSE (
"Attempting to add daughter particle: " << daughter );
976 ATH_MSG_WARNING (
"Attempting to add daughter particle: " << daughter );
979 G4PrimaryParticle *daughterG4Particle = this->getDaughterG4PrimaryParticle( daughter, makeLinkToTruth );
981 G4PrimaryParticle *daughterG4Particle = this->getDaughterG4PrimaryParticle( *daughter, makeLinkToTruth );
983 if (!daughterG4Particle) {
984 ATH_MSG_FATAL(
"Bailing out of loop over daughters due to errors.");
986 g4particle->SetDaughter( daughterG4Particle );
997 description << G4String(
"getG4PrimaryParticle: ") +
"No ISF::TruthBinding associated with ISParticle (" << isp <<
")";
998 G4Exception(
"iGeant4::TransportTool",
"NoISFTruthBinding", FatalException,
description);
1004 const G4ParticleDefinition *particleDefinition = this->getG4ParticleDefinition(isp.
pdgCode());
1006 if (particleDefinition==
nullptr) {
1007 ATH_MSG_ERROR(
"ISF_to_G4Event particle conversion failed. ISF_Particle PDG code = " << isp.
pdgCode() <<
1008 "\n This usually indicates a problem with the evgen step.\n" <<
1009 "Please report this to the Generators group, mentioning the release and generator used for evgen and the PDG code above." );
1018 if (useHepMC && currentGenPart) {
1019 auto ¤tGenPartMomentum = currentGenPart->momentum();
1020 px = currentGenPartMomentum.x();
1021 py = currentGenPartMomentum.y();
1022 pz = currentGenPartMomentum.z();
1025 auto &ispMomentum = isp.
momentum();
1026 px = ispMomentum.x();
1027 py = ispMomentum.y();
1028 pz = ispMomentum.z();
1031 std::unique_ptr<G4PrimaryParticle> g4particle = std::make_unique<G4PrimaryParticle>(particleDefinition,
px,
py,
pz);
1033 std::unique_ptr<PrimaryParticleInformation> primaryPartInfo = std::make_unique<PrimaryParticleInformation>(primaryGenpart,&isp);
1044 if ( currentGenPart ) {
1045 if (currentGenPart->end_vertex()) {
1048 ATH_MSG_ERROR (
"getG4PrimaryParticle(): GenParticle has a valid end GenVertexPtr!" );
1049 ATH_MSG_ERROR (
"getG4PrimaryParticle(): currentGenPart: " << currentGenPart <<
", barcode: " <<
HepMC::barcode(currentGenPart) );
1050 ATH_MSG_ERROR (
"getG4PrimaryParticle(): currentGenPart->end_vertex(): " << currentGenPart->end_vertex() <<
", barcode: " <<
HepMC::barcode(currentGenPart->end_vertex()) );
1051 ATH_MSG_FATAL (
"getG4PrimaryParticle(): Passing GenParticles with a valid end GenVertexPtr as input is no longer supported." );
1055 && !currentGenPart->end_vertex()) {
1059 auto A_part = currentGenPart->attribute<HepMC::ShadowParticle>(
"ShadowParticle");
1065 ATH_MSG_FATAL (
"Found a GenParticle with no matching GenParticle in the shadowGenEvent - something is wrong here!");
1068 if (!shadowPart->end_vertex()) {
1069 ATH_MSG_FATAL (
"Found status==2 shadow GenParticle with no end vertex - something is wrong here!");
1073 processPredefinedDecays(shadowPart, isp, g4particle.get());
1075 processPredefinedDecays(shadowPart, isp, g4particle.get(),
false);
1080 const double pmass = g4particle->GetMass();
1081 CLHEP::Hep3Vector gpv = g4particle->GetMomentum();
1082 double g4px=g4particle->GetMomentum().x();
1083 double g4py=g4particle->GetMomentum().y();
1084 double g4pz=g4particle->GetMomentum().z();
1092 px=currentGenPart->momentum().px();
1093 py=currentGenPart->momentum().py();
1094 pz=currentGenPart->momentum().pz();
1111 const double pe = std::sqrt(
mag2 + pmass*pmass);
1113 double originalEnergy=currentGenPart->momentum().e();
1114 if (originalEnergy>0.01) {
1115 if ((originalEnergy-
pe)/originalEnergy>0.01) {
1116 double genpx=currentGenPart->momentum().px();
1117 double genpy=currentGenPart->momentum().py();
1118 double genpz=currentGenPart->momentum().pz();
1119 double genp=sqrt(genpx*genpx + genpy*genpy + genpz*genpz);
1120 ATH_MSG_WARNING(
"Truth change in energy for: " << currentGenPart<<
" Morg="<<currentGenPart->momentum().m()<<
" Mmod="<<pmass<<
" Eorg="<<originalEnergy<<
" Emod="<<
pe<<
" porg="<<genp<<
" pmod="<<gpv.mag());
1125 auto& currentGenPart_nc = currentGenPart;
1127 auto* currentGenPart_nc = currentGenPart;
1129 currentGenPart_nc->set_momentum(HepMC::FourVector(
px,
py,
pz,
pe));
1137 g4particle->SetUserInformation(primaryPartInfo.release());
1139 return g4particle.release();
1155 G4PrimaryParticle *g4particle = this->getG4PrimaryParticle( isp, useHepMC, shadowGenEvent );
1157 ATH_MSG_ERROR(
"Failed to create G4PrimaryParticle for ISParticle (" << isp <<
")");
1162 G4PrimaryVertex *g4vertex =
new G4PrimaryVertex(isp.
position().x(),
1166 g4vertex->SetPrimary( g4particle );
1169 g4evt->AddPrimaryVertex( g4vertex );
1178 const G4ThreeVector g4Pos(
pos.x(),
pos.y(),
pos.z() );
1179 EInside insideStatus = worldSolid->Inside( g4Pos );
1181 bool insideWorld = insideStatus != kOutside;