ATLAS Offline Software
Loading...
Searching...
No Matches
ISF::InputConverter Class Referencefinal

Convert simulation input collections to ISFParticles for subsequent ISF simulation. More...

#include <InputConverter.h>

Inheritance diagram for ISF::InputConverter:
Collaboration diagram for ISF::InputConverter:

Public Member Functions

 InputConverter (const std::string &name, ISvcLocator *svc)
 Constructor.
virtual ~InputConverter ()
 Destructor.
virtual StatusCode initialize () override final
 Athena algtool Hooks.
virtual StatusCode finalize () override final
 Athena algtool Hook.
virtual StatusCode convert (McEventCollection &inputGenEvents, ISF::ISFParticleContainer &simParticles) const override final
 Convert selected particles from the given McEventCollection into ISFParticles and push them into the given ISFParticleContainer.
virtual StatusCode convertHepMCToG4Event (McEventCollection &inputGenEvents, G4Event &outputG4Event, McEventCollection &shadowGenEvents) const override final
virtual StatusCode convertHepMCToG4EventLegacy (McEventCollection &inputGenEvents, G4Event &outputG4Event) const override final
void ISF_to_G4Event (G4Event &event, const std::vector< ISF::ISFParticle * > &isp, HepMC::GenEvent *genEvent, HepMC::GenEvent *shadowGenEvent=nullptr, bool useHepMC=false) const override final
 Converts vector of ISF::ISFParticles to G4Event.

Private Member Functions

const G4ParticleDefinition * getG4ParticleDefinition (int pdgcode) const
G4PrimaryParticle * getDaughterG4PrimaryParticle (const HepMC::ConstGenParticlePtr &gp) const
G4PrimaryParticle * getDaughterG4PrimaryParticle (const HepMC::GenParticlePtr &gp, bool makeLinkToTruth=true) const
G4PrimaryParticle * getG4PrimaryParticle (ISF::ISFParticle &isp, bool useHepMC, HepMC::GenEvent *shadowGenEvent) const
void addG4PrimaryVertex (G4Event &g4evt, ISF::ISFParticle &isp, bool useHepMC, HepMC::GenEvent *shadowGenEvent) const
void processPredefinedDecays (const HepMC::ConstGenParticlePtr &genpart, ISF::ISFParticle &isp, G4PrimaryParticle *g4particle) const
void processPredefinedDecays (const HepMC::GenParticlePtr &genpart, ISF::ISFParticle &isp, G4PrimaryParticle *g4particle, bool makeLinkToTruth=true) const
bool matchedGenParticles (const HepMC::ConstGenParticlePtr &p1, const HepMC::ConstGenParticlePtr &p2) const
HepMC::GenParticlePtr findShadowParticle (const HepMC::ConstGenParticlePtr &genParticle, HepMC::GenEvent *shadowGenEvent) const
bool isInsideG4WorldVolume (const ISF::ISFParticle &isp, const G4VSolid *worldSolid) const
 Tests whether the given ISFParticle is within the Geant4 world volume.
double getParticleMass (const HepMC::ConstGenParticlePtr &p) const
 get right GenParticle mass
std::vector< HepMC::GenParticlePtrgetSelectedParticles (HepMC::GenEvent &evnt, bool legacyOrdering=false) const
 get all generator particles which pass filters
bool passesFilters (const HepMC::ConstGenParticlePtr &p) const
 check if the given particle passes all filters
ISF::ISFParticleconvertParticle (const HepMC::GenParticlePtr &genPartPtr) const
 convert GenParticle to ISFParticle

Private Attributes

std::shared_ptr< GenDatam_gendata {nullptr}
 ParticlePropertyService and ParticleDataTable.
bool m_useGeneratedParticleMass
 use GenParticle::generated_mass() in simulation
ToolHandleArray< IGenParticleFilterm_genParticleFilters
 HepMC::GenParticle filters.
bool m_quasiStableParticlesIncluded
BooleanProperty m_useShadowEvent {this, "UseShadowEvent", false, "New approach to selecting particles for simulation" }

Detailed Description

Convert simulation input collections to ISFParticles for subsequent ISF simulation.

Author
Elmar.Ritsch -at- cern.ch

Definition at line 45 of file InputConverter.h.

Constructor & Destructor Documentation

◆ InputConverter()

ISF::InputConverter::InputConverter ( const std::string & name,
ISvcLocator * svc )

Constructor.

Definition at line 47 of file InputConverter.cxx.

48 : base_class(name, svc)
52{
53 // particle mass from particle data table?
54 declareProperty("UseGeneratedParticleMass",
56 "Use particle mass assigned to GenParticle.");
57 // particle filters
58 declareProperty("GenParticleFilters",
60 "Tools for filtering out GenParticles.");
61 declareProperty("QuasiStableParticlesIncluded", m_quasiStableParticlesIncluded);
62}
ToolHandleArray< IGenParticleFilter > m_genParticleFilters
HepMC::GenParticle filters.
bool m_useGeneratedParticleMass
use GenParticle::generated_mass() in simulation

◆ ~InputConverter()

ISF::InputConverter::~InputConverter ( )
virtual

Destructor.

Definition at line 66 of file InputConverter.cxx.

67{
68}

Member Function Documentation

◆ addG4PrimaryVertex()

void ISF::InputConverter::addG4PrimaryVertex ( G4Event & g4evt,
ISF::ISFParticle & isp,
bool useHepMC,
HepMC::GenEvent * shadowGenEvent ) const
private

Definition at line 920 of file InputConverter.cxx.

922 {
923 /*
924 see conversion from PrimaryParticleInformation to TrackInformation in
925 http://acode-browser.usatlas.bnl.gov/lxr/source/atlas/Simulation/G4Atlas/G4AtlasAlg/src/AthenaStackingAction.cxx#0044
926
927 need to check with
928 http://acode-browser.usatlas.bnl.gov/lxr/source/atlas/Simulation/G4Atlas/G4AtlasAlg/src/TruthHepMCEventConverter.cxx#0151
929
930 that we don't miss something
931 */
932
933 G4PrimaryParticle *g4particle = this->getG4PrimaryParticle( isp, useHepMC, shadowGenEvent );
934 if (!g4particle) {
935 ATH_MSG_ERROR("Failed to create G4PrimaryParticle for ISParticle (" << isp <<")");
936 return;
937 }// Already printed a warning
938
939 // create a new vertex
940 G4PrimaryVertex *g4vertex = new G4PrimaryVertex(isp.position().x(),
941 isp.position().y(),
942 isp.position().z(),
943 isp.timeStamp());
944 g4vertex->SetPrimary( g4particle );
945 ATH_MSG_VERBOSE("Print G4PrimaryVertex: ");
946 if (msgLevel(MSG::VERBOSE)) { g4vertex->Print(); }
947 g4evt.AddPrimaryVertex(g4vertex);
948 return;
949}
#define ATH_MSG_ERROR(x)
#define ATH_MSG_VERBOSE(x)
double timeStamp() const
Timestamp of the ISFParticle.
const Amg::Vector3D & position() const
The current position of the ISFParticle.
G4PrimaryParticle * getG4PrimaryParticle(ISF::ISFParticle &isp, bool useHepMC, HepMC::GenEvent *shadowGenEvent) const

◆ convert()

StatusCode ISF::InputConverter::convert ( McEventCollection & inputGenEvents,
ISF::ISFParticleContainer & simParticles ) const
finaloverridevirtual

Convert selected particles from the given McEventCollection into ISFParticles and push them into the given ISFParticleContainer.

Definition at line 102 of file InputConverter.cxx.

104{
105 for ( auto eventPtr : inputGenEvents ) {
106 // skip empty events
107 if (eventPtr == nullptr) { continue; }
108
109 ATH_MSG_DEBUG("Starting conversion of GenEvent with"
110 " signal_process_id=" << HepMC::signal_process_id(eventPtr) <<
111 " and event_number=" << eventPtr->event_number() );
112
113 // new collection containing all gen particles that passed filters
114 bool legacyOrdering = true; // FIXME: this is only to keep the same order of particles
115 // as the prior 'StackFiller' implementation
116 // TODO: remove this functionality from the
117 // getSelectedParticles function once
118 // happy with the results
119 const auto passedGenParticles = getSelectedParticles(*eventPtr, legacyOrdering);
120
121 for ( auto& genPartPtr : passedGenParticles ) {
122 ATH_MSG_VERBOSE("Picking up following GenParticle for conversion to ISFParticle: " << genPartPtr);
123 auto simParticlePtr = convertParticle(genPartPtr);
124 if (!simParticlePtr) {
125 ATH_MSG_ERROR("Error while trying to convert input generator particles. Aborting.");
126 return StatusCode::FAILURE;
127 }
128 // add to collection that will be returned
129 simParticles.push_back(simParticlePtr);
130
131 } // loop over passed gen particles
132
133 } // loop over gen events
134
135 ATH_MSG_DEBUG( "Created initial simulation particle collection with size " << simParticles.size() );
136
137 return StatusCode::SUCCESS;
138}
#define ATH_MSG_DEBUG(x)
std::vector< HepMC::GenParticlePtr > getSelectedParticles(HepMC::GenEvent &evnt, bool legacyOrdering=false) const
get all generator particles which pass filters
ISF::ISFParticle * convertParticle(const HepMC::GenParticlePtr &genPartPtr) const
convert GenParticle to ISFParticle
int signal_process_id(const GenEvent &evt)
Definition GenEvent.h:572

◆ convertHepMCToG4Event()

StatusCode ISF::InputConverter::convertHepMCToG4Event ( McEventCollection & inputGenEvents,
G4Event & outputG4Event,
McEventCollection & shadowGenEvents ) const
finaloverridevirtual

Definition at line 140 of file InputConverter.cxx.

142 {
143 ISF::ISFParticleContainer simParticleList{}; // particles for ISF simulation
144 ATH_CHECK(this->convert(inputGenEvents, simParticleList));
145 //Convert from ISFParticleContainer to ConstISFParticleVector
146 ISF::ISFParticleVector simParticleVector{
147 std::make_move_iterator(std::begin(simParticleList)),
148 std::make_move_iterator(std::end(simParticleList))
149 };
150 HepMC::GenEvent* shadowGenEvent =
151 !shadowGenEvents.empty()
152 ? static_cast<HepMC::GenEvent*>(shadowGenEvents.back())
153 : nullptr;
154 this->ISF_to_G4Event(outputG4Event, simParticleVector, inputGenEvents.back(),
155 shadowGenEvent);
156 return StatusCode::SUCCESS;
157}
#define ATH_CHECK
Evaluate an expression and check for errors.
const T * back() const
Access the last element in the collection as an rvalue.
bool empty() const noexcept
Returns true if the collection is empty.
void ISF_to_G4Event(G4Event &event, const std::vector< ISF::ISFParticle * > &isp, HepMC::GenEvent *genEvent, HepMC::GenEvent *shadowGenEvent=nullptr, bool useHepMC=false) const override final
Converts vector of ISF::ISFParticles to G4Event.
virtual StatusCode convert(McEventCollection &inputGenEvents, ISF::ISFParticleContainer &simParticles) const override final
Convert selected particles from the given McEventCollection into ISFParticles and push them into the ...
HepMC3::GenEvent GenEvent
Definition GenEvent.h:39
std::list< ISF::ISFParticle * > ISFParticleContainer
generic ISFParticle container (not necessarily a std::list!)
std::vector< ISF::ISFParticle * > ISFParticleVector
ISFParticle vector.

◆ convertHepMCToG4EventLegacy()

StatusCode ISF::InputConverter::convertHepMCToG4EventLegacy ( McEventCollection & inputGenEvents,
G4Event & outputG4Event ) const
finaloverridevirtual

Definition at line 159 of file InputConverter.cxx.

160 {
161 ISF::ISFParticleContainer simParticleList{}; // particles for ISF simulation
162 ATH_CHECK(this->convert(inputGenEvents, simParticleList));
163 //Convert from ISFParticleContainer to ConstISFParticleVector
164 ISF::ISFParticleVector simParticleVector{
165 std::make_move_iterator(std::begin(simParticleList)),
166 std::make_move_iterator(std::end(simParticleList))
167 };
168 this->ISF_to_G4Event(outputG4Event, simParticleVector, inputGenEvents.back(),
169 nullptr);
170 return StatusCode::SUCCESS;
171}

◆ convertParticle()

ISF::ISFParticle * ISF::InputConverter::convertParticle ( const HepMC::GenParticlePtr & genPartPtr) const
private

convert GenParticle to ISFParticle

get all generator particles which pass filters

particle origin (TODO: add proper GeoID, collision/cosmics)

Definition at line 229 of file InputConverter.cxx.

229 {
230 if (!genPartPtr) { return nullptr; }
231
232 auto pVertex = genPartPtr->production_vertex();
233 if (!pVertex) {
234 ATH_MSG_ERROR("Unable to convert following generator particle due to missing production vertex for: " << genPartPtr);
235 return nullptr;
236 }
237 auto parentEvent = genPartPtr->parent_event();
238 if (!parentEvent) {
239 ATH_MSG_ERROR("Cannot convert a GenParticle without a parent GenEvent into an ISFParticle!!!");
240 return nullptr;
241 }
242
243 const Amg::Vector3D pos(pVertex->position().x(), pVertex->position().y(), pVertex->position().z());
244 const auto& pMomentum(genPartPtr->momentum());
245 const Amg::Vector3D mom(pMomentum.px(), pMomentum.py(), pMomentum.pz());
246 const double pMass = this->getParticleMass(genPartPtr);
247 double e=pMomentum.e();
248 if (e>1) { //only test for >1 MeV in momentum
249 double px=pMomentum.px();
250 double py=pMomentum.py();
251 double pz=pMomentum.pz();
252 double teste=std::sqrt(px*px + py*py + pz*pz + pMass*pMass);
253 if (std::abs(e-teste)>0.01*e) {
254 ATH_MSG_WARNING("Difference in energy for: " << genPartPtr<<" Morg="<<pMomentum.m()<<" Mmod="<<pMass<<" Eorg="<<e<<" Emod="<<teste);
255 }
256 if (MC::isDecayed(genPartPtr) && pVertex && genPartPtr->end_vertex()) { //check for possible changes of gamma for quasi stable particles
257 const auto& prodVtx = genPartPtr->production_vertex()->position();
258 const auto& endVtx = genPartPtr->end_vertex()->position();
259 CLHEP::Hep3Vector dist3D(endVtx.x()-prodVtx.x(), endVtx.y()-prodVtx.y(), endVtx.z()-prodVtx.z());
260
261 if (dist3D.mag()>1*Gaudi::Units::mm) {
262 CLHEP::HepLorentzVector mom( pMomentum.x(), pMomentum.y(), pMomentum.z(), pMomentum.t() );
263 double gamma_org=mom.gamma();
264 mom.setE(teste);
265 double gamma_new=mom.gamma();
266
267 if (std::abs(gamma_new-gamma_org)/(gamma_new+gamma_org)>0.001) {
268 ATH_MSG_WARNING("Difference in boost gamma for Quasi stable particle "<<genPartPtr);
269 ATH_MSG_WARNING(" gamma(m="<<mom.m()<<")="<<gamma_org<<" gamma(m="<<pMass<<")="<<gamma_new);
270 } else {
271 ATH_MSG_VERBOSE("Quasi stable particle "<<genPartPtr);
272 ATH_MSG_VERBOSE(" gamma(m="<<mom.m()<<")="<<gamma_org<<" gamma(m="<<pMass<<")="<<gamma_new);
273 }
274 }
275 }
276 }
277
278 const int pPdgId = genPartPtr->pdg_id();
279 const double charge = MC::charge(pPdgId);
280 const double pTime = pVertex->position().t() / Gaudi::Units::c_light;
283 const auto pBarcode = HepMC::barcode(genPartPtr);
284 const auto particleID = HepMC::uniqueID(genPartPtr);
285 auto tBinding = std::make_unique<ISF::TruthBinding>(genPartPtr);
286
287 auto hmpl = std::make_unique<HepMcParticleLink>(particleID, parentEvent->event_number(), HepMcParticleLink::IS_EVENTNUM, HepMcParticleLink::IS_ID);
288
289 auto sParticle = std::make_unique<ISF::ISFParticle>( std::move(pos),
290 std::move(mom),
291 pMass,
292 charge,
293 pPdgId,
294 genPartPtr->status(),
295 pTime,
296 origin,
297 particleID,
298 pBarcode,
299 tBinding.release(),
300 hmpl.release() );
301 return sParticle.release();
302}
#define ATH_MSG_WARNING(x)
double charge(const T &p)
Definition AtlasPID.h:997
double getParticleMass(const HepMC::ConstGenParticlePtr &p) const
get right GenParticle mass
Eigen::Matrix< double, 3, 1 > Vector3D
constexpr double pMass
int barcode(const T *p)
Definition Barcode.h:15
int uniqueID(const T &p)
std::pair< AtlasDetDescr::AtlasRegion, ISF::SimSvcID > DetRegionSvcIDPair
the datatype to be used to store each individual particle hop
Definition ISFParticle.h:30
@ fEventGeneratorSimID
Definition SimSvcID.h:34
bool isDecayed(const T &p)
Identify if the particle decayed.
double charge(const T &p)

◆ finalize()

StatusCode ISF::InputConverter::finalize ( )
finaloverridevirtual

Athena algtool Hook.

Definition at line 92 of file InputConverter.cxx.

93{
94 ATH_MSG_DEBUG("Finalizing ...");
95 return StatusCode::SUCCESS;
96}

◆ findShadowParticle()

HepMC::GenParticlePtr ISF::InputConverter::findShadowParticle ( const HepMC::ConstGenParticlePtr & genParticle,
HepMC::GenEvent * shadowGenEvent ) const
private

Definition at line 607 of file InputConverter.cxx.

608{
609 if (!shadowGenEvent) {
610 ATH_MSG_FATAL ("Found status==2 GenParticle with no end vertex and shadow GenEvent is missing - something is wrong here!");
611 abort();
612 }
613 // TODO in the future switch to using an Attribute which stores the shadow GenParticlePtr directly.
614 const int shadowId = genParticle->attribute<HepMC3::IntAttribute>(HepMCStr::ShadowParticleId)->value();
615 for (auto& shadowParticle : shadowGenEvent->particles()) {
616 if (shadowParticle->id() == shadowId && matchedGenParticles(genParticle, shadowParticle) ) { return shadowParticle; }
617 }
618 return std::make_shared<HepMC::GenParticle>();
619}
#define ATH_MSG_FATAL(x)
bool matchedGenParticles(const HepMC::ConstGenParticlePtr &p1, const HepMC::ConstGenParticlePtr &p2) const

◆ getDaughterG4PrimaryParticle() [1/2]

G4PrimaryParticle * ISF::InputConverter::getDaughterG4PrimaryParticle ( const HepMC::ConstGenParticlePtr & gp) const
private

Definition at line 517 of file InputConverter.cxx.

517 {
518 ATH_MSG_VERBOSE("Creating G4PrimaryParticle from GenParticle.");
519
520 const G4ParticleDefinition *particleDefinition = this->getG4ParticleDefinition(genpart->pdg_id());
521
522 if (particleDefinition==nullptr) {
523 ATH_MSG_ERROR("ISF_to_G4Event particle conversion failed. ISF_Particle PDG code = " << genpart->pdg_id() <<
524 "\n This usually indicates a problem with the evgen step.\n" <<
525 "Please report this to the Generators group, mentioning the release and generator used for evgen and the PDG code above." );
526 return nullptr;
527 }
528
529 // create new primaries and set them to the vertex
530 // G4double mass = particleDefinition->GetPDGMass();
531 auto &genpartMomentum = genpart->momentum();
532 G4double px = genpartMomentum.x();
533 G4double py = genpartMomentum.y();
534 G4double pz = genpartMomentum.z();
535
536 std::unique_ptr<G4PrimaryParticle> g4particle = std::make_unique<G4PrimaryParticle>(particleDefinition,px,py,pz);
537
538 if (genpart->end_vertex()) {
539 // Set the lifetime appropriately - this is slow but rigorous, and we
540 // don't want to end up with something like vertex time that we have
541 // to validate for every generator on earth...
542 const auto& prodVtx = genpart->production_vertex()->position();
543 const auto& endVtx = genpart->end_vertex()->position();
544 //const G4LorentzVector lv0 ( prodVtx.x(), prodVtx.y(), prodVtx.z(), prodVtx.t() );
545 //const G4LorentzVector lv1 ( endVtx.x(), endVtx.y(), endVtx.z(), endVtx.t() );
546 //Old calculation, not taken because vertex information is not sufficiently precise
547 //g4particle->SetProperTime( (lv1-lv0).mag()/Gaudi::Units::c_light );
548
549 CLHEP::Hep3Vector dist3D(endVtx.x()-prodVtx.x(), endVtx.y()-prodVtx.y(), endVtx.z()-prodVtx.z());
550 double tau=SetProperTimeFromDetectorFrameDecayLength(*g4particle,dist3D.mag());
551
552 if (msgLvl(MSG::VERBOSE)) {
553 double pmag2=g4particle->GetTotalMomentum(); //magnitude of particle momentum
554 pmag2*=pmag2; //magnitude of particle momentum squared
555 double e2=g4particle->GetTotalEnergy(); //energy of particle
556 e2*=e2; //energy of particle squared
557 double beta2=pmag2/e2; //beta^2=v^2/c^2 for particle
558 double tau2=dist3D.mag2()*(1/beta2-1)/Gaudi::Units::c_light/Gaudi::Units::c_light;
559 ATH_MSG_VERBOSE("lifetime tau(beta)="<<std::sqrt(tau2)<<" tau="<<tau);
560 }
562 ATH_MSG_VERBOSE( "Detected primary particle with end vertex." );
563 ATH_MSG_VERBOSE( "Will add the primary particle set on." );
564 ATH_MSG_VERBOSE( "Primary Particle: " << genpart );
565 ATH_MSG_VERBOSE( "Number of daughters : " << genpart->end_vertex()->particles_out().size()<<" at position "<<genpart->end_vertex() );
566 }
567 else {
568 ATH_MSG_WARNING( "Detected primary particle with end vertex." );
569 ATH_MSG_WARNING( "Will add the primary particle set on." );
570 ATH_MSG_WARNING( "Primary Particle: " << genpart );
571 ATH_MSG_WARNING( "Number of daughters: " << genpart->end_vertex()->particles_out().size()<<" at position "<<genpart->end_vertex() );
572 }
573 // Add all necessary daughter particles
574 for ( auto daughter: genpart->end_vertex()->particles_out() ) {
576 ATH_MSG_VERBOSE ( "Attempting to add daughter particle: " << daughter );
577 }
578 else {
579 ATH_MSG_WARNING ( "Attempting to add daughter particle: " << daughter );
580 }
581 G4PrimaryParticle *daughterG4Particle = this->getDaughterG4PrimaryParticle( daughter );
582 if (!daughterG4Particle) {
583 ATH_MSG_ERROR("Bailing out of loop over daughters of particle due to errors - will not return G4Particle.");
584 return nullptr;
585 }
586 g4particle->SetDaughter( daughterG4Particle );
587 }
588 }
589
590 return g4particle.release();
591}
double SetProperTimeFromDetectorFrameDecayLength(G4PrimaryParticle &g4particle, const double GeneratorDecayLength)
const G4ParticleDefinition * getG4ParticleDefinition(int pdgcode) const
G4PrimaryParticle * getDaughterG4PrimaryParticle(const HepMC::ConstGenParticlePtr &gp) const
double e2(const xAOD::CaloCluster &cluster)
return the uncorrected cluster energy in 2nd sampling
tuple genpart
Check that the actual generators, tune, and main PDF are consistent with the JO name.

◆ getDaughterG4PrimaryParticle() [2/2]

G4PrimaryParticle * ISF::InputConverter::getDaughterG4PrimaryParticle ( const HepMC::GenParticlePtr & gp,
bool makeLinkToTruth = true ) const
private

Definition at line 431 of file InputConverter.cxx.

431 {
432 ATH_MSG_VERBOSE("Creating G4PrimaryParticle from GenParticle.");
433
434 const G4ParticleDefinition *particleDefinition = this->getG4ParticleDefinition(genpart->pdg_id());
435
436 if (particleDefinition==nullptr) {
437 ATH_MSG_ERROR("ISF_to_G4Event particle conversion failed. ISF_Particle PDG code = " << genpart->pdg_id() <<
438 "\n This usually indicates a problem with the evgen step.\n" <<
439 "Please report this to the Generators group, mentioning the release and generator used for evgen and the PDG code above." );
440 return nullptr;
441 }
442
443 // create new primaries and set them to the vertex
444 // G4double mass = particleDefinition->GetPDGMass();
445 auto &genpartMomentum = genpart->momentum();
446 G4double px = genpartMomentum.x();
447 G4double py = genpartMomentum.y();
448 G4double pz = genpartMomentum.z();
449
450 std::unique_ptr<G4PrimaryParticle> g4particle = std::make_unique<G4PrimaryParticle>(particleDefinition,px,py,pz);
451
452 if (genpart->end_vertex()) {
453 // Set the lifetime appropriately - this is slow but rigorous, and we
454 // don't want to end up with something like vertex time that we have
455 // to validate for every generator on earth...
456 const auto& prodVtx = genpart->production_vertex()->position();
457 const auto& endVtx = genpart->end_vertex()->position();
458 //const G4LorentzVector lv0 ( prodVtx.x(), prodVtx.y(), prodVtx.z(), prodVtx.t() );
459 //const G4LorentzVector lv1 ( endVtx.x(), endVtx.y(), endVtx.z(), endVtx.t() );
460 //Old calculation, not taken because vertex information is not sufficiently precise
461 //g4particle->SetProperTime( (lv1-lv0).mag()/Gaudi::Units::c_light );
462
463 CLHEP::Hep3Vector dist3D(endVtx.x()-prodVtx.x(), endVtx.y()-prodVtx.y(), endVtx.z()-prodVtx.z());
464 double tau=SetProperTimeFromDetectorFrameDecayLength(*g4particle,dist3D.mag());
465
466 if (msgLvl(MSG::VERBOSE)) {
467 double pmag2=g4particle->GetTotalMomentum(); //magnitude of particle momentum
468 pmag2*=pmag2; //magnitude of particle momentum squared
469 double e2=g4particle->GetTotalEnergy(); //energy of particle
470 e2*=e2; //energy of particle squared
471 double beta2=pmag2/e2; //beta^2=v^2/c^2 for particle
472 double tau2=dist3D.mag2()*(1/beta2-1)/Gaudi::Units::c_light/Gaudi::Units::c_light;
473 ATH_MSG_VERBOSE("lifetime tau(beta)="<<std::sqrt(tau2)<<" tau="<<tau);
474 }
476 ATH_MSG_VERBOSE( "Detected primary particle with end vertex." );
477 ATH_MSG_VERBOSE( "Will add the primary particle set on." );
478 ATH_MSG_VERBOSE( "Primary Particle: " << genpart );
479 ATH_MSG_VERBOSE( "Number of daughters: " << genpart->end_vertex()->particles_out().size()<<" at position "<<genpart->end_vertex() );
480 }
481 else {
482 ATH_MSG_WARNING( "Detected primary particle with end vertex." );
483 ATH_MSG_WARNING( "Will add the primary particle set on." );
484 ATH_MSG_WARNING( "Primary Particle: " << genpart );
485 ATH_MSG_WARNING( "Number of daughters : " << genpart->end_vertex()->particles_out().size()<<" at position "<<genpart->end_vertex() );
486 }
487 // Add all necessary daughter particles
488 for ( auto daughter: genpart->end_vertex()->particles_out() ) {
490 ATH_MSG_VERBOSE ( "Attempting to add daughter particle : " << daughter );
491 }
492 else {
493 ATH_MSG_WARNING ( "Attempting to add daughter particle: " << daughter );
494 }
495 G4PrimaryParticle *daughterG4Particle = this->getDaughterG4PrimaryParticle( daughter, makeLinkToTruth );
496 if (!daughterG4Particle) {
497 ATH_MSG_ERROR("Bailing out of loop over daughters of particle: "<<
498 " due to errors - will not return G4Particle.");
499 return nullptr;
500 }
501 g4particle->SetDaughter( daughterG4Particle );
502 }
503 }
504
505 if (makeLinkToTruth) {
506 // Set the user information for this primary to point to the HepMcParticleLink...
507 std::unique_ptr<PrimaryParticleInformation> primaryPartInfo = std::make_unique<PrimaryParticleInformation>(genpart);
508 primaryPartInfo->SetRegenerationNr(0);
509 ATH_MSG_VERBOSE("Making primary down the line with barcode " << primaryPartInfo->GetParticleUniqueID());
510 g4particle->SetUserInformation(primaryPartInfo.release());
511 }
512
513 return g4particle.release();
514}

◆ getG4ParticleDefinition()

const G4ParticleDefinition * ISF::InputConverter::getG4ParticleDefinition ( int pdgcode) const
private

Special cases for Geantinos

Standard particles

Definition at line 397 of file InputConverter.cxx.

398{
400 if (pdgcode==998) {
401 return G4ChargedGeantino::Definition();
402 }
403 if (pdgcode==999) {
404 return G4Geantino::GeantinoDefinition();
405 }
407 G4ParticleTable *ptable = G4ParticleTable::GetParticleTable();
408 if (ptable) {
409 return ptable->FindParticle(pdgcode);
410 }
411 ATH_MSG_ERROR("getG4ParticleDefinition - Failed to retrieve G4ParticleTable!");
412 return nullptr;
413}

◆ getG4PrimaryParticle()

G4PrimaryParticle * ISF::InputConverter::getG4PrimaryParticle ( ISF::ISFParticle & isp,
bool useHepMC,
HepMC::GenEvent * shadowGenEvent ) const
private

In the case that particles are being passed back to Geant4 then we may have particles which have already interacted, so we should set the regeneration number accordingly.

Definition at line 779 of file InputConverter.cxx.

780{
781 ATH_MSG_VERBOSE("Creating G4PrimaryParticle from ISFParticle.");
782
783 auto* truthBinding = isp.getTruthBinding();
784 if (!truthBinding) {
785 G4ExceptionDescription description;
786 description << G4String("getG4PrimaryParticle: ") + "No ISF::TruthBinding associated with ISParticle (" << isp <<")";
787 G4Exception("iGeant4::TransportTool", "NoISFTruthBinding", FatalException, description);
788 return nullptr; //The G4Exception call above should abort the job, but Coverity does not seem to pick this up.
789 }
790 HepMC::GenParticlePtr currentGenPart = truthBinding->getCurrentGenParticle();
791 HepMC::GenParticlePtr primaryGenpart = truthBinding->getPrimaryGenParticle();
792
793 const G4ParticleDefinition *particleDefinition = this->getG4ParticleDefinition(isp.pdgCode());
794
795 if (particleDefinition==nullptr) {
796 ATH_MSG_ERROR("ISF_to_G4Event particle conversion failed. ISF_Particle PDG code = " << isp.pdgCode() <<
797 "\n This usually indicates a problem with the evgen step.\n" <<
798 "Please report this to the Generators group, mentioning the release and generator used for evgen and the PDG code above." );
799 return nullptr;
800 }
801
802 // create new primaries and set them to the vertex
803 // G4double mass = particleDefinition->GetPDGMass();
804 G4double px(0.0);
805 G4double py(0.0);
806 G4double pz(0.0);
807 if (useHepMC && currentGenPart) {
808 auto &currentGenPartMomentum = currentGenPart->momentum();
809 px = currentGenPartMomentum.x();
810 py = currentGenPartMomentum.y();
811 pz = currentGenPartMomentum.z();
812 }
813 else {
814 auto &ispMomentum = isp.momentum();
815 px = ispMomentum.x();
816 py = ispMomentum.y();
817 pz = ispMomentum.z();
818 }
819
820 std::unique_ptr<G4PrimaryParticle> g4particle = std::make_unique<G4PrimaryParticle>(particleDefinition,px,py,pz);
821 // UserInformation
822 std::unique_ptr<PrimaryParticleInformation> primaryPartInfo = std::make_unique<PrimaryParticleInformation>(primaryGenpart,&isp);
823
827 const int regenerationNr = HepMC::StatusBased::generations(&isp);
828 if (HepMC::BarcodeBased::generations(&isp) != regenerationNr) {
829 ATH_MSG_WARNING ("StatusBased::generations() = " << regenerationNr << ", BarcodeBased::generations() = " << HepMC::BarcodeBased::generations(&isp) << ", isp: " << isp);
830 }
831 primaryPartInfo->SetRegenerationNr(regenerationNr);
832
833 if ( currentGenPart ) {
834 if (currentGenPart->end_vertex()) {
835 // Old approach particle had an end vertex - predefined decays taken from the main GenEvent
836 // No longer supported
837 ATH_MSG_ERROR ( "getG4PrimaryParticle(): GenParticle has a valid end GenVertexPtr!" );
838 ATH_MSG_ERROR ( "getG4PrimaryParticle(): currentGenPart: " << currentGenPart << ", barcode: " << HepMC::barcode(currentGenPart) );
839 ATH_MSG_ERROR ( "getG4PrimaryParticle(): currentGenPart->end_vertex(): " << currentGenPart->end_vertex() << ", barcode: " << HepMC::barcode(currentGenPart->end_vertex()) );
840 ATH_MSG_FATAL ( "getG4PrimaryParticle(): Passing GenParticles with a valid end GenVertexPtr as input is no longer supported." );
841 abort();
842 }
843 else if (MC::isDecayed(currentGenPart) // Some assumptions about main GenEvent here
844 && !currentGenPart->end_vertex()) {
845 // New approach - predefined decays taken from shadow GenEvent
846 // Find the matching particle in the shadowGenEvent
847 auto A_part = currentGenPart->attribute<HepMC::ShadowParticle>(HepMCStr::ShadowParticle);
848 HepMC::ConstGenParticlePtr shadowPart = (A_part) ? A_part->value() : findShadowParticle(currentGenPart, shadowGenEvent);
849 if (!shadowPart) {
850 ATH_MSG_FATAL ("Found a GenParticle with no matching GenParticle in the shadowGenEvent - something is wrong here!");
851 abort();
852 }
853 if (!shadowPart->end_vertex()) {
854 ATH_MSG_FATAL ("Found status==2 shadow GenParticle with no end vertex - something is wrong here!");
855 abort();
856 }
857 processPredefinedDecays(shadowPart, isp, g4particle.get());
858 }
859
860 double px,py,pz;
861 const double pmass = g4particle->GetMass();
862 CLHEP::Hep3Vector gpv = g4particle->GetMomentum();
863 double g4px=g4particle->GetMomentum().x();
864 double g4py=g4particle->GetMomentum().y();
865 double g4pz=g4particle->GetMomentum().z();
866 if (useHepMC) {
867 //Code adapted from TruthHepMCEventConverter::TransformHepMCParticle
868 px=g4px;
869 py=g4py;
870 pz=g4pz;
871 } else {
872 //Take mass from g4particle, put keep momentum as in currentGenPart
873 px=currentGenPart->momentum().px();
874 py=currentGenPart->momentum().py();
875 pz=currentGenPart->momentum().pz();
876 //Now a dirty hack to keep backward compatibility in the truth:
877 //When running AtlasG4 or FullG4 between 21.0.41 and 21.0.111, the currentGenPart 3-momentum and mass was reset to the values from the g4particle
878 //together with the mass of the g4particle after the 1st initialization of the g4particle from the genevent. This is done for a consistent mass
879 //value in the truth record compared to the used g4 mass. Since g4particles don't store the 3-momentum directly, but rather a
880 //unit direction vector, the mass and the kinetic energy, this reduces the numeric accuracy.
881 //For backward compatibility, if all 3-momentum components agree to the g4particle momentum within 1 keV, we keep
882 //this old method. This comparison is needed, since in ISF this code could be rerun after the ID or CALO simulation, where
883 //real energy was lost in previous detectors and hence currentGenPart should NOT be changed to some g4particle values!
884 //TODO: find a way to implement this in a backward compatible way in ISF::InputConverter::convertParticle(HepMC::GenParticlePtr genPartPtr)
885 if (std::abs(px-g4px)<CLHEP::keV && std::abs(py-g4py)<CLHEP::keV && std::abs(pz-g4pz)<CLHEP::keV) {
886 px=g4px;
887 py=g4py;
888 pz=g4pz;
889 }
890 }
891 const double mag2=px*px + py*py + pz*pz;
892 const double pe = std::sqrt(mag2 + pmass*pmass); // this does only change for boosts, etc.
893
894 double originalEnergy=currentGenPart->momentum().e();
895 if (originalEnergy>0.01) { //only test for >1 MeV in momentum
896 if ((originalEnergy-pe)/originalEnergy>0.01) {
897 double genpx=currentGenPart->momentum().px();
898 double genpy=currentGenPart->momentum().py();
899 double genpz=currentGenPart->momentum().pz();
900 double genp=sqrt(genpx*genpx + genpy*genpy + genpz*genpz);
901 ATH_MSG_WARNING("Truth change in energy for: " << currentGenPart<<" Morg="<<currentGenPart->momentum().m()<<" Mmod="<<pmass<<" Eorg="<<originalEnergy<<" Emod="<<pe<<" porg="<<genp<<" pmod="<<gpv.mag());
902 }
903 }
904
905 auto& currentGenPart_nc = currentGenPart;
906 currentGenPart_nc->set_momentum(HepMC::FourVector(px,py,pz,pe));
907 } // Truth was detected
908
909 ATH_MSG_VERBOSE("PrimaryParticleInformation:");
910 ATH_MSG_VERBOSE(" GetParticleUniqueID = " << primaryPartInfo->GetParticleUniqueID());
911 ATH_MSG_VERBOSE(" GetRegenerationNr = " << primaryPartInfo->GetRegenerationNr());
912 ATH_MSG_VERBOSE(" GetHepMCParticle = " << primaryPartInfo->GetHepMCParticle());
913 ATH_MSG_VERBOSE(" GetISFParticle = " << primaryPartInfo->GetISFParticle());
914 g4particle->SetUserInformation(primaryPartInfo.release());
915
916 return g4particle.release();
917}
Scalar mag2() const
mag2 method - forward to squaredNorm()
const TruthBinding * getTruthBinding() const
pointer to the simulation truth - optional, can be 0
const Amg::Vector3D & momentum() const
The current momentum vector of the ISFParticle.
int pdgCode() const
PDG value.
void processPredefinedDecays(const HepMC::ConstGenParticlePtr &genpart, ISF::ISFParticle &isp, G4PrimaryParticle *g4particle) const
HepMC::GenParticlePtr findShadowParticle(const HepMC::ConstGenParticlePtr &genParticle, HepMC::GenEvent *shadowGenEvent) const
HepMC::GenParticlePtr getCurrentGenParticle()
pointer to the particle in the simulation truth
std::string description
glabal timer - how long have I taken so far?
Definition hcg.cxx:93
int generations(const T &p)
Method to return how many interactions a particle has undergone during simulation (only to be used in...
int generations(const T &p)
Method to return how many interactions a particle has undergone during simulation based on the status...
HepMC3::FourVector FourVector
HepMC3::GenParticlePtr GenParticlePtr
Definition GenParticle.h:19
HepMC3::ConstGenParticlePtr ConstGenParticlePtr
Definition GenParticle.h:20

◆ getParticleMass()

double ISF::InputConverter::getParticleMass ( const HepMC::ConstGenParticlePtr & p) const
private

get right GenParticle mass

Definition at line 307 of file InputConverter.cxx.

307 {
308 // default value: generated particle mass
309 double mass = part->generated_mass();
310 ATH_MSG_VERBOSE("part->generated_mass, mass="<<mass);
311
312 // 1. use PDT mass?
314 const int absPDG = std::abs(part->pdg_id());
315 auto pData = m_gendata->particleMass(absPDG);
316 if (pData) {
317 mass = pData.value();
318 ATH_MSG_VERBOSE("using pData mass, mass="<<mass);
319 }
320 else {
321 ATH_MSG_WARNING( "Unable to find mass of particle with PDG ID '" << absPDG << "' in ParticleDataTable. Will set mass to generated_mass: " << mass);
322 }
323 }
324 return mass;
325}
std::shared_ptr< GenData > m_gendata
ParticlePropertyService and ParticleDataTable.

◆ getSelectedParticles()

std::vector< HepMC::GenParticlePtr > ISF::InputConverter::getSelectedParticles ( HepMC::GenEvent & evnt,
bool legacyOrdering = false ) const
private

get all generator particles which pass filters

Definition at line 175 of file InputConverter.cxx.

175 {
176 auto allGenPartBegin = evnt.particles().begin();
177 auto allGenPartEnd = evnt.particles().end();
178
179 // reserve destination container with maximum size, i.e. number of particles in input event
180 std::vector<HepMC::GenParticlePtr> passedGenParticles{};
181 size_t maxParticles = std::distance(allGenPartBegin, allGenPartEnd);
182 passedGenParticles.reserve(maxParticles);
183
184 if (m_useShadowEvent) {
185 if (legacyOrdering) {
186 // FIXME: remove this block and the 'legacyOrdering' flag
187 // once we don't need the legacy order any longer
188 for (auto vtx: evnt.vertices() ) {
189 std::copy_if (vtx->particles_out().begin(),
190 vtx->particles_out().end(),
191 std::back_inserter(passedGenParticles),
192 [](HepMC::GenParticlePtr p){return p->attribute<HepMC3::IntAttribute>(HepMCStr::ShadowParticleId);});
193 }
194 }
195 else {
196 std::copy_if (allGenPartBegin,
197 allGenPartEnd,
198 std::back_inserter(passedGenParticles),
199 [](HepMC::GenParticlePtr p){return p->attribute<HepMC3::IntAttribute>(HepMCStr::ShadowParticleId);});
200 }
201 }
202 else {
203 if (legacyOrdering) {
204 // FIXME: remove this block and the 'legacyOrdering' flag
205 // once we don't need the legacy order any longer
206 for (auto vtx: evnt.vertices() ) {
207 std::copy_if (vtx->particles_out().begin(),
208 vtx->particles_out().end(),
209 std::back_inserter(passedGenParticles),
210 [this](HepMC::GenParticlePtr p){return this->passesFilters(std::const_pointer_cast<const HepMC3::GenParticle>(p));});
211 }
212 }
213 else {
214 std::copy_if (allGenPartBegin,
215 allGenPartEnd,
216 std::back_inserter(passedGenParticles),
217 [this](HepMC::GenParticlePtr p){return this->passesFilters(std::const_pointer_cast<const HepMC3::GenParticle>(p));});
218 }
219 }
220
221 passedGenParticles.shrink_to_fit();
222
223 return passedGenParticles;
224}
BooleanProperty m_useShadowEvent
bool passesFilters(const HepMC::ConstGenParticlePtr &p) const
check if the given particle passes all filters
constexpr uint8_t maxParticles()

◆ initialize()

StatusCode ISF::InputConverter::initialize ( )
finaloverridevirtual

Athena algtool Hooks.

Definition at line 73 of file InputConverter.cxx.

74{
75 ATH_MSG_VERBOSE("initialize() begin");
76
78 m_gendata = std::make_shared<GenData>();
79 }
80
81 if (!m_genParticleFilters.empty()) {
83 }
84
85 ATH_MSG_VERBOSE("initialize() successful");
86 return StatusCode::SUCCESS;
87}

◆ ISF_to_G4Event()

void ISF::InputConverter::ISF_to_G4Event ( G4Event & event,
const std::vector< ISF::ISFParticle * > & isp,
HepMC::GenEvent * genEvent,
HepMC::GenEvent * shadowGenEvent = nullptr,
bool useHepMC = false ) const
finaloverride

Converts vector of ISF::ISFParticles to G4Event.

Definition at line 358 of file InputConverter.cxx.

361 {
362 // G4Event *g4evt = new G4Event(ctx.eventID().event_number());
363
364 // retrieve world solid (volume)
365 const G4VSolid *worldSolid = G4TransportationManager::GetTransportationManager()->GetNavigatorForTracking()->GetWorldVolume()->GetLogicalVolume()->GetSolid();
366
367 for ( ISF::ISFParticle *ispPtr: ispVector ) {
368 ISF::ISFParticle &isp = *ispPtr;
369 if ( !isInsideG4WorldVolume(isp, worldSolid) ) {
370 ATH_MSG_WARNING("Unable to convert ISFParticle to G4PrimaryParticle!");
371 ATH_MSG_WARNING(" ISFParticle: " << isp );
372 if (worldSolid) {
373 ATH_MSG_WARNING(" is outside Geant4 world volume: ");
374 worldSolid->DumpInfo();
375 G4cout << std::flush;
376 }
377 else {
378 ATH_MSG_WARNING(" is outside Geant4 world volume.");
379 }
380 continue;
381 }
382 this->addG4PrimaryVertex(event, isp, useHepMC, shadowGenEvent);
383 }
384
385 AtlasG4EventUserInfo* atlasG4EvtUserInfo =
386 dynamic_cast<AtlasG4EventUserInfo*>(event.GetUserInformation());
387 if (!atlasG4EvtUserInfo) {
388 atlasG4EvtUserInfo = new AtlasG4EventUserInfo(Gaudi::Hive::currentContext());
389 event.SetUserInformation(atlasG4EvtUserInfo);
390 }
391 atlasG4EvtUserInfo->SetLastProcessedTrackID(0); // TODO Check if it is better to set this to -1 initially
392 atlasG4EvtUserInfo->SetLastProcessedStep(0); // TODO Check if it is better to set this to -1 initially
393 atlasG4EvtUserInfo->SetHepMCEvent(genEvent);
394}
void SetHepMCEvent(HepMC::GenEvent *)
set m_theEvent, the pointer to the HepMC::GenEvent used to create the G4Event.
void SetLastProcessedStep(int stepNumber)
record value of the G4Track::GetCurrentStepNumber() for the current G4Step.
void SetLastProcessedTrackID(int trackID)
record the value of G4Track::GetTrackID() for the current G4Step.
void addG4PrimaryVertex(G4Event &g4evt, ISF::ISFParticle &isp, bool useHepMC, HepMC::GenEvent *shadowGenEvent) const
bool isInsideG4WorldVolume(const ISF::ISFParticle &isp, const G4VSolid *worldSolid) const
Tests whether the given ISFParticle is within the Geant4 world volume.

◆ isInsideG4WorldVolume()

bool ISF::InputConverter::isInsideG4WorldVolume ( const ISF::ISFParticle & isp,
const G4VSolid * worldSolid ) const
private

Tests whether the given ISFParticle is within the Geant4 world volume.

Definition at line 952 of file InputConverter.cxx.

953{
954
955 const Amg::Vector3D &pos = isp.position();
956 const G4ThreeVector g4Pos( pos.x(), pos.y(), pos.z() );
957 EInside insideStatus = worldSolid->Inside( g4Pos );
958
959 bool insideWorld = insideStatus != kOutside;
960 return insideWorld;
961}

◆ matchedGenParticles()

bool ISF::InputConverter::matchedGenParticles ( const HepMC::ConstGenParticlePtr & p1,
const HepMC::ConstGenParticlePtr & p2 ) const
private

Definition at line 594 of file InputConverter.cxx.

596{
597 return (HepMC::barcode(p1) == HepMC::barcode(p2))
598 && (p1->status() == p2->status())
599 && (p1->pdg_id() == p2->pdg_id())
600 && ((p1->momentum().px()) == (p2->momentum().px()))
601 && ((p1->momentum().py()) == (p2->momentum().py()))
602 && ((p1->momentum().pz()) == (p2->momentum().pz()))
603 && (float(p1->momentum().m()) == float(p2->momentum().m()));
604}

◆ passesFilters()

bool ISF::InputConverter::passesFilters ( const HepMC::ConstGenParticlePtr & p) const
private

check if the given particle passes all filters

Definition at line 331 of file InputConverter.cxx.

332{
333 // TODO: implement this as a std::find_if with a lambda function
334 for ( const auto& filter : m_genParticleFilters ) {
335 // determine if the particle passes current filter
336 bool passFilter = filter->pass(part);
337 ATH_MSG_VERBOSE("GenParticleFilter '" << filter.typeAndName() << "' returned: "
338 << (passFilter ? "true, will keep particle."
339 : "false, will remove particle."));
340 const auto& momentum = part->momentum();
341 ATH_MSG_VERBOSE("Particle: ("
342 <<momentum.px()<<", "
343 <<momentum.py()<<", "
344 <<momentum.pz()<<"), pdgCode: "
345 <<part->pdg_id() );
346
347 if (!passFilter) {
348 return false;
349 }
350 }
351
352 return true;
353}

◆ processPredefinedDecays() [1/2]

void ISF::InputConverter::processPredefinedDecays ( const HepMC::ConstGenParticlePtr & genpart,
ISF::ISFParticle & isp,
G4PrimaryParticle * g4particle ) const
private

Set the lifetime appropriately - this is slow but rigorous, and we don't want to end up with something like vertex time that we have to validate for every generator on earth...

Definition at line 622 of file InputConverter.cxx.

623{
627 const auto& prodVtx = genpart->production_vertex()->position();
628 const auto& endVtx = genpart->end_vertex()->position();
629
630 CLHEP::Hep3Vector dist3D(endVtx.x()-prodVtx.x(), endVtx.y()-prodVtx.y(), endVtx.z()-prodVtx.z());
631 double tau=SetProperTimeFromDetectorFrameDecayLength(*g4particle,dist3D.mag());
632
633 if (msgLvl(MSG::VERBOSE)) {
634 double pmag2=g4particle->GetTotalMomentum(); //magnitude of particle momentum
635 pmag2*=pmag2; //magnitude of particle momentum squared
636 double e2=g4particle->GetTotalEnergy(); //energy of particle
637 e2*=e2; //energy^2 of particle
638 double beta2=pmag2/e2; //beta^2=v^2/c^2 for particle
639 double mass2=g4particle->GetMass(); //mass of particle
640 mass2*=mass2; //mass^2 of particle
641
642 double tau2=dist3D.mag2()*(1/beta2-1)/Gaudi::Units::c_squared;
643
644 const G4LorentzVector lv0( prodVtx.x(), prodVtx.y(), prodVtx.z(), prodVtx.t() );
645 const G4LorentzVector lv1( endVtx.x(), endVtx.y(), endVtx.z(), endVtx.t() );
646 //Old calculation, not taken because vertex information is not sufficiently precise
647 //g4particle->SetProperTime( (lv1-lv0).mag()/Gaudi::Units::c_light );
648 G4LorentzVector dist4D(lv1);
649 dist4D-=lv0;
650
651 double dist4Dgamma=std::numeric_limits<double>::infinity();
652 if (dist4D.t()>0 && dist4D.mag2()>0) {
653 dist4Dgamma=dist4D.gamma();
654 } else {
655 ATH_MSG_VERBOSE( "dist4D t="<<dist4D.t()<<" mag2="<<dist4D.mag2());
656 }
657
658 G4LorentzVector fourmom(g4particle->GetMomentum(),g4particle->GetTotalEnergy());
659 double fourmomgamma=std::numeric_limits<double>::infinity();
660 if (fourmom.t()>0 && fourmom.mag2()>0) {
661 fourmomgamma=fourmom.gamma();
662 } else {
663 ATH_MSG_VERBOSE( "fourmom t="<<fourmom.t()<<" mag2="<<fourmom.mag2());
664 }
665
666 ATH_MSG_VERBOSE( "gammaVertex="<<dist4Dgamma<<" gammamom="<<fourmomgamma<<" gamma(beta)="<<1/std::sqrt(1-beta2)<<" lifetime tau(beta)="<<std::sqrt(tau2)<<" lifetime tau="<<tau);
667 }
669 ATH_MSG_VERBOSE( "Detected primary particle with end vertex." );
670 ATH_MSG_VERBOSE( "Will add the primary particle set on." );
671 ATH_MSG_VERBOSE( "ISF Particle: " << isp );
672 ATH_MSG_VERBOSE( "Primary Particle: " << genpart );
673 ATH_MSG_VERBOSE( "Number of daughters: " << genpart->end_vertex()->particles_out().size() << " at position "<< genpart->end_vertex() );
674 }
675 else {
676 ATH_MSG_WARNING( "Detected primary particle with end vertex. This should only be the case if" );
677 ATH_MSG_WARNING( "you are running with quasi-stable particle simulation enabled. This is not" );
678 ATH_MSG_WARNING( "yet validated - you'd better know what you're doing. Will add the primary" );
679 ATH_MSG_WARNING( "particle set on." );
680 ATH_MSG_WARNING( "ISF Particle: " << isp );
681 ATH_MSG_WARNING( "Primary Particle: " << genpart );
682 ATH_MSG_VERBOSE( "Number of daughters: " << genpart->end_vertex()->particles_out().size() );
683 }
684 // Add all necessary daughter particles
685 for ( auto daughter: *(genpart->end_vertex())) {
687 ATH_MSG_VERBOSE ( "Attempting to add daughter particle" << daughter );
688 }
689 else {
690 ATH_MSG_WARNING ( "Attempting to add daughter particle: " << daughter );
691 }
692 G4PrimaryParticle *daughterG4Particle = this->getDaughterG4PrimaryParticle( daughter );
693 if (!daughterG4Particle) {
694 ATH_MSG_FATAL("Bailing out of loop over daughters due to errors.");
695 }
696 g4particle->SetDaughter( daughterG4Particle );
697 }
698}

◆ processPredefinedDecays() [2/2]

void ISF::InputConverter::processPredefinedDecays ( const HepMC::GenParticlePtr & genpart,
ISF::ISFParticle & isp,
G4PrimaryParticle * g4particle,
bool makeLinkToTruth = true ) const
private

Set the lifetime appropriately - this is slow but rigorous, and we don't want to end up with something like vertex time that we have to validate for every generator on earth...

Definition at line 701 of file InputConverter.cxx.

702{
706 const auto& prodVtx = genpart->production_vertex()->position();
707 const auto& endVtx = genpart->end_vertex()->position();
708
709 CLHEP::Hep3Vector dist3D(endVtx.x()-prodVtx.x(), endVtx.y()-prodVtx.y(), endVtx.z()-prodVtx.z());
710 double tau=SetProperTimeFromDetectorFrameDecayLength(*g4particle,dist3D.mag());
711
712 if (msgLvl(MSG::VERBOSE)) {
713 double pmag2=g4particle->GetTotalMomentum(); //magnitude of particle momentum
714 pmag2*=pmag2; //magnitude of particle momentum squared
715 double e2=g4particle->GetTotalEnergy(); //energy of particle
716 e2*=e2; //energy^2 of particle
717 double beta2=pmag2/e2; //beta^2=v^2/c^2 for particle
718 double mass2=g4particle->GetMass(); //mass of particle
719 mass2*=mass2; //mass^2 of particle
720
721 double tau2=dist3D.mag2()*(1/beta2-1)/Gaudi::Units::c_squared;
722
723 const G4LorentzVector lv0( prodVtx.x(), prodVtx.y(), prodVtx.z(), prodVtx.t() );
724 const G4LorentzVector lv1( endVtx.x(), endVtx.y(), endVtx.z(), endVtx.t() );
725 //Old calculation, not taken because vertex information is not sufficiently precise
726 //g4particle->SetProperTime( (lv1-lv0).mag()/Gaudi::Units::c_light );
727 G4LorentzVector dist4D(lv1);
728 dist4D-=lv0;
729
730 double dist4Dgamma=std::numeric_limits<double>::infinity();
731 if (dist4D.t()>0 && dist4D.mag2()>0) {
732 dist4Dgamma=dist4D.gamma();
733 } else {
734 ATH_MSG_VERBOSE( "dist4D t="<<dist4D.t()<<" mag2="<<dist4D.mag2());
735 }
736
737 G4LorentzVector fourmom(g4particle->GetMomentum(),g4particle->GetTotalEnergy());
738 double fourmomgamma=std::numeric_limits<double>::infinity();
739 if (fourmom.t()>0 && fourmom.mag2()>0) {
740 fourmomgamma=fourmom.gamma();
741 } else {
742 ATH_MSG_VERBOSE( "fourmom t="<<fourmom.t()<<" mag2="<<fourmom.mag2());
743 }
744
745 ATH_MSG_VERBOSE( "gammaVertex="<<dist4Dgamma<<" gammamom="<<fourmomgamma<<" gamma(beta)="<<1/std::sqrt(1-beta2)<<" lifetime tau(beta)="<<std::sqrt(tau2)<<" lifetime tau="<<tau);
746 }
748 ATH_MSG_VERBOSE( "Detected primary particle with end vertex." );
749 ATH_MSG_VERBOSE( "Will add the primary particle set on." );
750 ATH_MSG_VERBOSE( "ISF Particle: " << isp );
751 ATH_MSG_VERBOSE( "Primary Particle: " << genpart );
752 ATH_MSG_VERBOSE( "Number of daughters: " << genpart->end_vertex()->particles_out_size() << " at position "<< genpart->end_vertex() );
753 }
754 else {
755 ATH_MSG_WARNING( "Detected primary particle with end vertex. This should only be the case if" );
756 ATH_MSG_WARNING( "you are running with quasi-stable particle simulation enabled. This is not" );
757 ATH_MSG_WARNING( "yet validated - you'd better know what you're doing. Will add the primary" );
758 ATH_MSG_WARNING( "particle set on." );
759 ATH_MSG_WARNING( "ISF Particle: " << isp );
760 ATH_MSG_WARNING( "Primary Particle: " << genpart );
761 ATH_MSG_WARNING( "Number of daughters: " << genpart->end_vertex()->particles_out_size() );
762 }
763 // Add all necessary daughter particles
764 for ( auto daughter: *(genpart->end_vertex())) {
766 ATH_MSG_VERBOSE ( "Attempting to add daughter particle: " << daughter );
767 }
768 else {
769 ATH_MSG_WARNING ( "Attempting to add daughter particle: " << daughter );
770 }
771 G4PrimaryParticle *daughterG4Particle = this->getDaughterG4PrimaryParticle( daughter, makeLinkToTruth );
772 if (!daughterG4Particle) {
773 ATH_MSG_FATAL("Bailing out of loop over daughters due to errors.");
774 }
775 g4particle->SetDaughter( daughterG4Particle );
776 }
777}

Member Data Documentation

◆ m_gendata

std::shared_ptr<GenData> ISF::InputConverter::m_gendata {nullptr}
private

ParticlePropertyService and ParticleDataTable.

Definition at line 115 of file InputConverter.h.

115{nullptr};

◆ m_genParticleFilters

ToolHandleArray<IGenParticleFilter> ISF::InputConverter::m_genParticleFilters
private

HepMC::GenParticle filters.

Definition at line 119 of file InputConverter.h.

◆ m_quasiStableParticlesIncluded

bool ISF::InputConverter::m_quasiStableParticlesIncluded
private

Definition at line 121 of file InputConverter.h.

◆ m_useGeneratedParticleMass

bool ISF::InputConverter::m_useGeneratedParticleMass
private

use GenParticle::generated_mass() in simulation

Definition at line 117 of file InputConverter.h.

◆ m_useShadowEvent

BooleanProperty ISF::InputConverter::m_useShadowEvent {this, "UseShadowEvent", false, "New approach to selecting particles for simulation" }
private

Definition at line 122 of file InputConverter.h.

122{this, "UseShadowEvent", false, "New approach to selecting particles for simulation" };

The documentation for this class was generated from the following files: