12 #include "GaudiKernel/IPartPropSvc.h"
13 #include "GaudiKernel/PhysicalConstants.h"
29 #include "G4PrimaryParticle.hh"
31 #include "G4Geantino.hh"
32 #include "G4ChargedGeantino.hh"
33 #include "G4ParticleTable.hh"
34 #include "G4LorentzVector.hh"
35 #include "G4TransportationManager.hh"
41 #include "CLHEP/Geometry/Point3D.h"
42 #include "CLHEP/Geometry/Vector3D.h"
43 #include "CLHEP/Units/SystemOfUnits.h"
45 #include "HepPDT/ParticleID.hh"
46 #include "HepPDT/DecayData.hh"
47 #include "HepPDT/ParticleDataTable.hh"
53 , m_particlePropSvc(
"PartPropSvc",
name)
54 , m_particleDataTable(nullptr)
55 , m_useGeneratedParticleMass(false)
56 , m_genParticleFilters(this)
57 , m_quasiStableParticlesIncluded(false)
60 declareProperty(
"UseGeneratedParticleMass",
62 "Use particle mass assigned to GenParticle.");
64 declareProperty(
"GenParticleFilters",
66 "Tools for filtering out GenParticles.");
68 declareProperty(
"ParticlePropertyService",
70 "ParticlePropertyService to retrieve the PDT.");
88 if (!m_useGeneratedParticleMass) {
90 m_particleDataTable = m_particlePropSvc->PDT();
91 if (!m_particleDataTable) {
92 ATH_MSG_FATAL(
"Could not get ParticleDataTable from " << m_particlePropSvc <<
". Abort" );
93 return StatusCode::FAILURE;
97 if (!m_genParticleFilters.empty()) {
98 ATH_CHECK(m_genParticleFilters.retrieve());
102 return StatusCode::SUCCESS;
111 return StatusCode::SUCCESS;
121 for (
auto eventPtr : inputGenEvents ) {
123 if (eventPtr ==
nullptr) {
continue; }
127 " and event_number=" << eventPtr->event_number() );
130 bool legacyOrdering =
true;
137 for (
auto& genPartPtr : passedGenParticles ) {
138 ATH_MSG_VERBOSE(
"Picking up following GenParticle for conversion to ISFParticle: " << genPartPtr);
139 auto simParticlePtr = convertParticle(genPartPtr);
140 if (!simParticlePtr) {
141 ATH_MSG_ERROR(
"Error while trying to convert input generator particles. Aborting.");
142 return StatusCode::FAILURE;
145 simParticles.push_back(simParticlePtr);
151 ATH_MSG_DEBUG(
"Created initial simulation particle collection with size " << simParticles.size() );
153 return StatusCode::SUCCESS;
163 std::make_move_iterator(
std::begin(simParticleList)),
164 std::make_move_iterator(
std::end(simParticleList))
166 HepMC::GenEvent* shadowGenEvent =
167 !shadowGenEvents.
empty()
168 ?
static_cast<HepMC::GenEvent*
>(shadowGenEvents.
back())
170 this->ISF_to_G4Event(outputG4Event, simParticleVector, inputGenEvents.
back(),
172 return StatusCode::SUCCESS;
181 std::make_move_iterator(
std::begin(simParticleList)),
182 std::make_move_iterator(
std::end(simParticleList))
184 this->ISF_to_G4Event(outputG4Event, simParticleVector, inputGenEvents.
back(),
186 return StatusCode::SUCCESS;
191 std::vector<HepMC::GenParticlePtr>
193 auto allGenPartBegin = evnt.particles().begin();
194 auto allGenPartEnd = evnt.particles().end();
197 std::vector<HepMC::GenParticlePtr> passedGenParticles{};
201 if (m_useShadowEvent) {
202 if (legacyOrdering) {
205 for (
auto vtx: evnt.vertices() ) {
207 vtx->particles_out().end(),
208 std::back_inserter(passedGenParticles),
215 std::back_inserter(passedGenParticles),
220 if (legacyOrdering) {
223 for (
auto vtx: evnt.vertices() ) {
225 vtx->particles_out().end(),
226 std::back_inserter(passedGenParticles),
227 [
this](
HepMC::GenParticlePtr p){return this->passesFilters(std::const_pointer_cast<const HepMC3::GenParticle>(p));});
233 std::back_inserter(passedGenParticles),
234 [
this](
HepMC::GenParticlePtr p){
return this->passesFilters(std::const_pointer_cast<const HepMC3::GenParticle>(
p));});
238 passedGenParticles.shrink_to_fit();
240 return passedGenParticles;
243 std::vector<HepMC::GenParticlePtr>
245 auto allGenPartBegin = evnt.particles_begin();
246 auto allGenPartEnd = evnt.particles_end();
249 std::vector<HepMC::GenParticlePtr> passedGenParticles{};
253 if (legacyOrdering) {
256 auto vtxIt = evnt.vertices_begin();
257 auto vtxItEnd = evnt.vertices_end();
258 for ( ; vtxIt != vtxItEnd; ++vtxIt ) {
259 const auto vtxPtr = *vtxIt;
262 std::back_inserter(passedGenParticles),
269 std::back_inserter(passedGenParticles),
273 passedGenParticles.shrink_to_fit();
275 return passedGenParticles;
283 if (!genPartPtr) {
return nullptr; }
285 auto pVertex = genPartPtr->production_vertex();
287 ATH_MSG_ERROR(
"Unable to convert following generator particle due to missing production vertex for: " << genPartPtr);
290 auto parentEvent = genPartPtr->parent_event();
292 ATH_MSG_ERROR(
"Cannot convert a GenParticle without a parent GenEvent into an ISFParticle!!!");
296 const Amg::Vector3D pos(pVertex->position().x(), pVertex->position().y(), pVertex->position().z());
297 const auto& pMomentum(genPartPtr->momentum());
300 const double pMass = this->getParticleMass(genPartPtr);
302 const double pMass = this->getParticleMass(*genPartPtr);
304 double e=pMomentum.e();
306 double px=pMomentum.px();
307 double py=pMomentum.py();
308 double pz=pMomentum.pz();
310 if (std::abs(
e-teste)>0.01*
e) {
311 ATH_MSG_WARNING(
"Difference in energy for: " << genPartPtr<<
" Morg="<<pMomentum.m()<<
" Mmod="<<
pMass<<
" Eorg="<<
e<<
" Emod="<<teste);
313 if (
MC::isDecayed(genPartPtr) && pVertex && genPartPtr->end_vertex()) {
314 const auto& prodVtx = genPartPtr->production_vertex()->position();
315 const auto& endVtx = genPartPtr->end_vertex()->position();
316 CLHEP::Hep3Vector dist3D(endVtx.x()-prodVtx.x(), endVtx.y()-prodVtx.y(), endVtx.z()-prodVtx.z());
319 CLHEP::HepLorentzVector
mom( pMomentum.x(), pMomentum.y(), pMomentum.z(), pMomentum.t() );
320 double gamma_org=
mom.gamma();
322 double gamma_new=
mom.gamma();
324 if (std::abs(gamma_new-gamma_org)/(gamma_new+gamma_org)>0.001) {
325 ATH_MSG_WARNING(
"Difference in boost gamma for Quasi stable particle "<<genPartPtr);
335 const int pPdgId = genPartPtr->pdg_id();
342 auto tBinding = std::make_unique<ISF::TruthBinding>(genPartPtr);
346 auto sParticle = std::make_unique<ISF::ISFParticle>( std::move(
pos),
351 genPartPtr->status(),
358 return sParticle.release();
367 double mass =
part->generated_mass();
371 if ( !m_useGeneratedParticleMass ) {
372 const int absPDG = std::abs(
part->pdg_id());
373 HepPDT::ParticleData
const *pData = (m_particleDataTable)
374 ? m_particleDataTable->particle(absPDG)
377 mass = pData->mass();
381 ATH_MSG_WARNING(
"Unable to find mass of particle with PDG ID '" << absPDG <<
"' in ParticleDataTable. Will set mass to generated_mass: " <<
mass);
391 double mass =
part.generated_mass();
395 if ( !m_useGeneratedParticleMass ) {
396 const int absPDG = std::abs(
part.pdg_id());
397 HepPDT::ParticleData
const *pData = (m_particleDataTable)
398 ? m_particleDataTable->particle(absPDG)
401 mass = pData->mass();
405 ATH_MSG_WARNING(
"Unable to find mass of particle with PDG ID '" << absPDG <<
"' in ParticleDataTable. Will set mass to generated_mass: " <<
mass);
419 for (
const auto&
filter : m_genParticleFilters ) {
424 :
"false, will remove particle."));
444 for (
const auto&
filter : m_genParticleFilters ) {
449 :
"false, will remove particle."));
470 HepMC::GenEvent* genEvent, HepMC::GenEvent* shadowGenEvent,
471 bool useHepMC)
const {
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(
event, isp, useHepMC, shadowGenEvent);
497 if (!atlasG4EvtUserInfo) {
499 event.SetUserInformation(atlasG4EvtUserInfo);
511 return G4ChargedGeantino::Definition();
514 return G4Geantino::GeantinoDefinition();
517 G4ParticleTable *ptable = G4ParticleTable::GetParticleTable();
519 return ptable->FindParticle(pdgcode);
521 ATH_MSG_ERROR(
"getG4ParticleDefinition - Failed to retrieve G4ParticleTable!");
531 const double p2=
std::pow(g4particle.GetTotalMomentum(),2);
532 const double m2=
std::pow(g4particle.GetMass(),2);
533 const double l2=
std::pow(GeneratorDecayLength,2);
535 const double tau=std::sqrt(tau2);
536 g4particle.SetProperTime( tau );
545 const G4ParticleDefinition *particleDefinition = this->getG4ParticleDefinition(
genpart->pdg_id());
547 if (particleDefinition==
nullptr) {
548 ATH_MSG_ERROR(
"ISF_to_G4Event particle conversion failed. ISF_Particle PDG code = " <<
genpart->pdg_id() <<
549 "\n This usually indicates a problem with the evgen step.\n" <<
550 "Please report this to the Generators group, mentioning the release and generator used for evgen and the PDG code above." );
556 auto &genpartMomentum =
genpart->momentum();
557 G4double
px = genpartMomentum.x();
558 G4double
py = genpartMomentum.y();
559 G4double
pz = genpartMomentum.z();
561 std::unique_ptr<G4PrimaryParticle> g4particle = std::make_unique<G4PrimaryParticle>(particleDefinition,
px,
py,
pz);
567 const auto& prodVtx =
genpart->production_vertex()->position();
568 const auto& endVtx =
genpart->end_vertex()->position();
574 CLHEP::Hep3Vector dist3D(endVtx.x()-prodVtx.x(), endVtx.y()-prodVtx.y(), endVtx.z()-prodVtx.z());
578 double pmag2=g4particle->GetTotalMomentum();
580 double e2=g4particle->GetTotalEnergy();
582 double beta2=pmag2/
e2;
586 if (m_quasiStableParticlesIncluded) {
599 for (
auto daughter:
genpart->end_vertex()->particles_out() ) {
600 if (m_quasiStableParticlesIncluded) {
601 ATH_MSG_VERBOSE (
"Attempting to add daughter particle : " << daughter );
604 ATH_MSG_WARNING (
"Attempting to add daughter particle: " << daughter );
606 G4PrimaryParticle *daughterG4Particle = this->getDaughterG4PrimaryParticle( daughter, makeLinkToTruth );
607 if (!daughterG4Particle) {
608 ATH_MSG_ERROR(
"Bailing out of loop over daughters of particle: "<<
609 " due to errors - will not return G4Particle.");
612 g4particle->SetDaughter( daughterG4Particle );
616 if (makeLinkToTruth) {
618 std::unique_ptr<PrimaryParticleInformation> primaryPartInfo = std::make_unique<PrimaryParticleInformation>(
genpart);
621 g4particle->SetUserInformation(primaryPartInfo.release());
624 return g4particle.release();
631 const G4ParticleDefinition *particleDefinition = this->getG4ParticleDefinition(
genpart->pdg_id());
633 if (particleDefinition==
nullptr) {
634 ATH_MSG_ERROR(
"ISF_to_G4Event particle conversion failed. ISF_Particle PDG code = " <<
genpart->pdg_id() <<
635 "\n This usually indicates a problem with the evgen step.\n" <<
636 "Please report this to the Generators group, mentioning the release and generator used for evgen and the PDG code above." );
642 auto &genpartMomentum =
genpart->momentum();
643 G4double
px = genpartMomentum.x();
644 G4double
py = genpartMomentum.y();
645 G4double
pz = genpartMomentum.z();
647 std::unique_ptr<G4PrimaryParticle> g4particle = std::make_unique<G4PrimaryParticle>(particleDefinition,
px,
py,
pz);
653 const auto& prodVtx =
genpart->production_vertex()->position();
654 const auto& endVtx =
genpart->end_vertex()->position();
660 CLHEP::Hep3Vector dist3D(endVtx.x()-prodVtx.x(), endVtx.y()-prodVtx.y(), endVtx.z()-prodVtx.z());
664 double pmag2=g4particle->GetTotalMomentum();
666 double e2=g4particle->GetTotalEnergy();
668 double beta2=pmag2/
e2;
672 if (m_quasiStableParticlesIncluded) {
685 for (
auto daughter:
genpart->end_vertex()->particles_out() ) {
686 if (m_quasiStableParticlesIncluded) {
687 ATH_MSG_VERBOSE (
"Attempting to add daughter particle: " << daughter );
690 ATH_MSG_WARNING (
"Attempting to add daughter particle: " << daughter );
692 G4PrimaryParticle *daughterG4Particle = this->getDaughterG4PrimaryParticle( daughter );
693 if (!daughterG4Particle) {
694 ATH_MSG_ERROR(
"Bailing out of loop over daughters of particle due to errors - will not return G4Particle.");
697 g4particle->SetDaughter( daughterG4Particle );
701 return g4particle.release();
708 const G4ParticleDefinition *particleDefinition = this->getG4ParticleDefinition(
genpart.pdg_id());
710 if (particleDefinition==
nullptr) {
711 ATH_MSG_ERROR(
"ISF_to_G4Event particle conversion failed. ISF_Particle PDG code = " <<
genpart.pdg_id() <<
712 "\n This usually indicates a problem with the evgen step.\n" <<
713 "Please report this to the Generators group, mentioning the release and generator used for evgen and the PDG code above." );
719 auto &genpartMomentum =
genpart.momentum();
720 G4double
px = genpartMomentum.x();
721 G4double
py = genpartMomentum.y();
722 G4double
pz = genpartMomentum.z();
724 std::unique_ptr<G4PrimaryParticle> g4particle = std::make_unique<G4PrimaryParticle>(particleDefinition,
px,
py,
pz);
730 const auto& prodVtx =
genpart.production_vertex()->position();
731 const auto& endVtx =
genpart.end_vertex()->position();
737 CLHEP::Hep3Vector dist3D(endVtx.x()-prodVtx.x(), endVtx.y()-prodVtx.y(), endVtx.z()-prodVtx.z());
741 double pmag2=g4particle->GetTotalMomentum();
743 double e2=g4particle->GetTotalEnergy();
745 double beta2=pmag2/
e2;
749 if (m_quasiStableParticlesIncluded) {
762 for (
auto daughterIter=
genpart.end_vertex()->particles_out_const_begin();
763 daughterIter!=
genpart.end_vertex()->particles_out_const_end(); ++daughterIter ) {
764 if (m_quasiStableParticlesIncluded) {
765 ATH_MSG_VERBOSE (
"Attempting to add daughter particle: " << **daughterIter );
768 ATH_MSG_WARNING (
"Attempting to add daughter particle: " << **daughterIter );
770 G4PrimaryParticle *daughterG4Particle = this->getDaughterG4PrimaryParticle( **daughterIter, makeLinkToTruth );
771 if (!daughterG4Particle) {
772 ATH_MSG_ERROR(
"Bailing out of loop over daughters of particle due to errors - will not return G4Particle.");
775 g4particle->SetDaughter( daughterG4Particle );
779 if (makeLinkToTruth) {
781 std::unique_ptr<PrimaryParticleInformation> primaryPartInfo = std::make_unique<PrimaryParticleInformation>(&
genpart);
784 g4particle->SetUserInformation(primaryPartInfo.release());
787 return g4particle.release();
796 && (
p1->status() ==
p2->status())
797 && (
p1->pdg_id() ==
p2->pdg_id())
798 && ((
p1->momentum().px()) == (
p2->momentum().px()))
799 && ((
p1->momentum().py()) == (
p2->momentum().py()))
800 && ((
p1->momentum().pz()) == (
p2->momentum().pz()))
807 if (!shadowGenEvent) {
808 ATH_MSG_FATAL (
"Found status==2 GenParticle with no end vertex and shadow GenEvent is missing - something is wrong here!");
813 const int shadowId = genParticle->attribute<HepMC3::IntAttribute>(
"ShadowParticleId")->
value();
814 for (
auto& shadowParticle : shadowGenEvent->particles()) {
815 if (shadowParticle->id() == shadowId && matchedGenParticles(genParticle, shadowParticle) ) {
return shadowParticle; }
817 return std::make_shared<HepMC::GenParticle>();
819 for (HepMC::GenEvent::particle_iterator pitr=shadowGenEvent->particles_begin(); pitr != shadowGenEvent->particles_end(); ++pitr) {
821 if (matchedGenParticles(genParticle, shadowParticle) ) {
return shadowParticle; }
835 const auto& prodVtx =
genpart->production_vertex()->position();
836 const auto& endVtx =
genpart->end_vertex()->position();
838 CLHEP::Hep3Vector dist3D(endVtx.x()-prodVtx.x(), endVtx.y()-prodVtx.y(), endVtx.z()-prodVtx.z());
842 double pmag2=g4particle->GetTotalMomentum();
844 double e2=g4particle->GetTotalEnergy();
846 double beta2=pmag2/
e2;
847 double mass2=g4particle->GetMass();
852 const G4LorentzVector lv0( prodVtx.x(), prodVtx.y(), prodVtx.z(), prodVtx.t() );
853 const G4LorentzVector lv1( endVtx.x(), endVtx.y(), endVtx.z(), endVtx.t() );
856 G4LorentzVector dist4D(lv1);
859 double dist4Dgamma=std::numeric_limits<double>::infinity();
860 if (dist4D.t()>0 && dist4D.mag2()>0) {
861 dist4Dgamma=dist4D.gamma();
866 G4LorentzVector fourmom(g4particle->GetMomentum(),g4particle->GetTotalEnergy());
867 double fourmomgamma=std::numeric_limits<double>::infinity();
868 if (fourmom.t()>0 && fourmom.mag2()>0) {
869 fourmomgamma=fourmom.gamma();
874 ATH_MSG_VERBOSE(
"gammaVertex="<<dist4Dgamma<<
" gammamom="<<fourmomgamma<<
" gamma(beta)="<<1/std::sqrt(1-beta2)<<
" lifetime tau(beta)="<<std::sqrt(tau2)<<
" lifetime tau="<<tau);
876 if (m_quasiStableParticlesIncluded) {
884 ATH_MSG_WARNING(
"Detected primary particle with end vertex. This should only be the case if" );
885 ATH_MSG_WARNING(
"you are running with quasi-stable particle simulation enabled. This is not" );
886 ATH_MSG_WARNING(
"yet validated - you'd better know what you're doing. Will add the primary" );
893 for (
auto daughter: *(
genpart->end_vertex())) {
894 if (m_quasiStableParticlesIncluded) {
898 ATH_MSG_WARNING (
"Attempting to add daughter particle: " << daughter );
900 G4PrimaryParticle *daughterG4Particle = this->getDaughterG4PrimaryParticle( daughter );
901 if (!daughterG4Particle) {
902 ATH_MSG_FATAL(
"Bailing out of loop over daughters due to errors.");
904 g4particle->SetDaughter( daughterG4Particle );
915 const auto& prodVtx =
genpart->production_vertex()->position();
916 const auto& endVtx =
genpart->end_vertex()->position();
918 CLHEP::Hep3Vector dist3D(endVtx.x()-prodVtx.x(), endVtx.y()-prodVtx.y(), endVtx.z()-prodVtx.z());
922 double pmag2=g4particle->GetTotalMomentum();
924 double e2=g4particle->GetTotalEnergy();
926 double beta2=pmag2/
e2;
927 double mass2=g4particle->GetMass();
932 const G4LorentzVector lv0( prodVtx.x(), prodVtx.y(), prodVtx.z(), prodVtx.t() );
933 const G4LorentzVector lv1( endVtx.x(), endVtx.y(), endVtx.z(), endVtx.t() );
936 G4LorentzVector dist4D(lv1);
939 double dist4Dgamma=std::numeric_limits<double>::infinity();
940 if (dist4D.t()>0 && dist4D.mag2()>0) {
941 dist4Dgamma=dist4D.gamma();
946 G4LorentzVector fourmom(g4particle->GetMomentum(),g4particle->GetTotalEnergy());
947 double fourmomgamma=std::numeric_limits<double>::infinity();
948 if (fourmom.t()>0 && fourmom.mag2()>0) {
949 fourmomgamma=fourmom.gamma();
954 ATH_MSG_VERBOSE(
"gammaVertex="<<dist4Dgamma<<
" gammamom="<<fourmomgamma<<
" gamma(beta)="<<1/std::sqrt(1-beta2)<<
" lifetime tau(beta)="<<std::sqrt(tau2)<<
" lifetime tau="<<tau);
956 if (m_quasiStableParticlesIncluded) {
964 ATH_MSG_WARNING(
"Detected primary particle with end vertex. This should only be the case if" );
965 ATH_MSG_WARNING(
"you are running with quasi-stable particle simulation enabled. This is not" );
966 ATH_MSG_WARNING(
"yet validated - you'd better know what you're doing. Will add the primary" );
973 for (
auto daughter: *(
genpart->end_vertex())) {
974 if (m_quasiStableParticlesIncluded) {
975 ATH_MSG_VERBOSE (
"Attempting to add daughter particle: " << daughter );
978 ATH_MSG_WARNING (
"Attempting to add daughter particle: " << daughter );
981 G4PrimaryParticle *daughterG4Particle = this->getDaughterG4PrimaryParticle( daughter, makeLinkToTruth );
983 G4PrimaryParticle *daughterG4Particle = this->getDaughterG4PrimaryParticle( *daughter, makeLinkToTruth );
985 if (!daughterG4Particle) {
986 ATH_MSG_FATAL(
"Bailing out of loop over daughters due to errors.");
988 g4particle->SetDaughter( daughterG4Particle );
999 description << G4String(
"getG4PrimaryParticle: ") +
"No ISF::TruthBinding associated with ISParticle (" << isp <<
")";
1000 G4Exception(
"iGeant4::TransportTool",
"NoISFTruthBinding", FatalException,
description);
1006 const G4ParticleDefinition *particleDefinition = this->getG4ParticleDefinition(isp.
pdgCode());
1008 if (particleDefinition==
nullptr) {
1009 ATH_MSG_ERROR(
"ISF_to_G4Event particle conversion failed. ISF_Particle PDG code = " << isp.
pdgCode() <<
1010 "\n This usually indicates a problem with the evgen step.\n" <<
1011 "Please report this to the Generators group, mentioning the release and generator used for evgen and the PDG code above." );
1020 if (useHepMC && currentGenPart) {
1021 auto ¤tGenPartMomentum = currentGenPart->momentum();
1022 px = currentGenPartMomentum.x();
1023 py = currentGenPartMomentum.y();
1024 pz = currentGenPartMomentum.z();
1027 auto &ispMomentum = isp.
momentum();
1028 px = ispMomentum.x();
1029 py = ispMomentum.y();
1030 pz = ispMomentum.z();
1033 std::unique_ptr<G4PrimaryParticle> g4particle = std::make_unique<G4PrimaryParticle>(particleDefinition,
px,
py,
pz);
1035 std::unique_ptr<PrimaryParticleInformation> primaryPartInfo = std::make_unique<PrimaryParticleInformation>(primaryGenpart,&isp);
1046 if ( currentGenPart ) {
1047 if (currentGenPart->end_vertex()) {
1050 ATH_MSG_ERROR (
"getG4PrimaryParticle(): GenParticle has a valid end GenVertexPtr!" );
1051 ATH_MSG_ERROR (
"getG4PrimaryParticle(): currentGenPart: " << currentGenPart <<
", barcode: " <<
HepMC::barcode(currentGenPart) );
1052 ATH_MSG_ERROR (
"getG4PrimaryParticle(): currentGenPart->end_vertex(): " << currentGenPart->end_vertex() <<
", barcode: " <<
HepMC::barcode(currentGenPart->end_vertex()) );
1053 ATH_MSG_FATAL (
"getG4PrimaryParticle(): Passing GenParticles with a valid end GenVertexPtr as input is no longer supported." );
1057 && !currentGenPart->end_vertex()) {
1061 auto A_part = currentGenPart->attribute<HepMC::ShadowParticle>(
"ShadowParticle");
1067 ATH_MSG_FATAL (
"Found a GenParticle with no matching GenParticle in the shadowGenEvent - something is wrong here!");
1070 if (!shadowPart->end_vertex()) {
1071 ATH_MSG_FATAL (
"Found status==2 shadow GenParticle with no end vertex - something is wrong here!");
1075 processPredefinedDecays(shadowPart, isp, g4particle.get());
1077 processPredefinedDecays(shadowPart, isp, g4particle.get(),
false);
1082 const double pmass = g4particle->GetMass();
1083 CLHEP::Hep3Vector gpv = g4particle->GetMomentum();
1084 double g4px=g4particle->GetMomentum().x();
1085 double g4py=g4particle->GetMomentum().y();
1086 double g4pz=g4particle->GetMomentum().z();
1094 px=currentGenPart->momentum().px();
1095 py=currentGenPart->momentum().py();
1096 pz=currentGenPart->momentum().pz();
1113 const double pe = std::sqrt(
mag2 + pmass*pmass);
1115 double originalEnergy=currentGenPart->momentum().e();
1116 if (originalEnergy>0.01) {
1117 if ((originalEnergy-
pe)/originalEnergy>0.01) {
1118 double genpx=currentGenPart->momentum().px();
1119 double genpy=currentGenPart->momentum().py();
1120 double genpz=currentGenPart->momentum().pz();
1121 double genp=sqrt(genpx*genpx + genpy*genpy + genpz*genpz);
1122 ATH_MSG_WARNING(
"Truth change in energy for: " << currentGenPart<<
" Morg="<<currentGenPart->momentum().m()<<
" Mmod="<<pmass<<
" Eorg="<<originalEnergy<<
" Emod="<<
pe<<
" porg="<<genp<<
" pmod="<<gpv.mag());
1127 auto& currentGenPart_nc = currentGenPart;
1129 auto* currentGenPart_nc = currentGenPart;
1131 currentGenPart_nc->set_momentum(HepMC::FourVector(
px,
py,
pz,
pe));
1139 g4particle->SetUserInformation(primaryPartInfo.release());
1141 return g4particle.release();
1147 HepMC::GenEvent* shadowGenEvent)
const {
1158 G4PrimaryParticle *g4particle = this->getG4PrimaryParticle( isp, useHepMC, shadowGenEvent );
1160 ATH_MSG_ERROR(
"Failed to create G4PrimaryParticle for ISParticle (" << isp <<
")");
1165 G4PrimaryVertex *g4vertex =
new G4PrimaryVertex(isp.
position().x(),
1169 g4vertex->SetPrimary( g4particle );
1172 g4evt.AddPrimaryVertex(g4vertex);
1181 const G4ThreeVector g4Pos(
pos.x(),
pos.y(),
pos.z() );
1182 EInside insideStatus = worldSolid->Inside( g4Pos );
1184 bool insideWorld = insideStatus != kOutside;