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

ServiceHandle< IPartPropSvc > m_particlePropSvc
 ParticlePropertyService and ParticleDataTable.
const HepPDT::ParticleDataTable * m_particleDataTable
 PDT used to look up particle masses.
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 48 of file InputConverter.h.

Constructor & Destructor Documentation

◆ InputConverter()

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

Constructor.

Definition at line 51 of file InputConverter.cxx.

52 : base_class(name, svc)
53 , m_particlePropSvc("PartPropSvc",name)
54 , m_particleDataTable(nullptr)
58{
59 // particle mass from particle data table?
60 declareProperty("UseGeneratedParticleMass",
62 "Use particle mass assigned to GenParticle.");
63 // particle filters
64 declareProperty("GenParticleFilters",
66 "Tools for filtering out GenParticles.");
67 // the particle property service
68 declareProperty("ParticlePropertyService",
70 "ParticlePropertyService to retrieve the PDT.");
71 declareProperty("QuasiStableParticlesIncluded", m_quasiStableParticlesIncluded);
72}
ServiceHandle< IPartPropSvc > m_particlePropSvc
ParticlePropertyService and ParticleDataTable.
ToolHandleArray< IGenParticleFilter > m_genParticleFilters
HepMC::GenParticle filters.
const HepPDT::ParticleDataTable * m_particleDataTable
PDT used to look up particle masses.
bool m_useGeneratedParticleMass
use GenParticle::generated_mass() in simulation

◆ ~InputConverter()

ISF::InputConverter::~InputConverter ( )
virtual

Destructor.

Definition at line 76 of file InputConverter.cxx.

77{
78}

Member Function Documentation

◆ addG4PrimaryVertex()

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

Definition at line 938 of file InputConverter.cxx.

940 {
941 /*
942 see conversion from PrimaryParticleInformation to TrackInformation in
943 http://acode-browser.usatlas.bnl.gov/lxr/source/atlas/Simulation/G4Atlas/G4AtlasAlg/src/AthenaStackingAction.cxx#0044
944
945 need to check with
946 http://acode-browser.usatlas.bnl.gov/lxr/source/atlas/Simulation/G4Atlas/G4AtlasAlg/src/TruthHepMCEventConverter.cxx#0151
947
948 that we don't miss something
949 */
950
951 G4PrimaryParticle *g4particle = this->getG4PrimaryParticle( isp, useHepMC, shadowGenEvent );
952 if (!g4particle) {
953 ATH_MSG_ERROR("Failed to create G4PrimaryParticle for ISParticle (" << isp <<")");
954 return;
955 }// Already printed a warning
956
957 // create a new vertex
958 G4PrimaryVertex *g4vertex = new G4PrimaryVertex(isp.position().x(),
959 isp.position().y(),
960 isp.position().z(),
961 isp.timeStamp());
962 g4vertex->SetPrimary( g4particle );
963 ATH_MSG_VERBOSE("Print G4PrimaryVertex: ");
964 if (msgLevel(MSG::VERBOSE)) { g4vertex->Print(); }
965 g4evt.AddPrimaryVertex(g4vertex);
966 return;
967}
#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 118 of file InputConverter.cxx.

120{
121 for ( auto eventPtr : inputGenEvents ) {
122 // skip empty events
123 if (eventPtr == nullptr) { continue; }
124
125 ATH_MSG_DEBUG("Starting conversion of GenEvent with"
126 " signal_process_id=" << HepMC::signal_process_id(eventPtr) <<
127 " and event_number=" << eventPtr->event_number() );
128
129 // new collection containing all gen particles that passed filters
130 bool legacyOrdering = true; // FIXME: this is only to keep the same order of particles
131 // as the prior 'StackFiller' implementation
132 // TODO: remove this functionality from the
133 // getSelectedParticles function once
134 // happy with the results
135 const auto passedGenParticles = getSelectedParticles(*eventPtr, legacyOrdering);
136
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;
143 }
144 // add to collection that will be returned
145 simParticles.push_back(simParticlePtr);
146
147 } // loop over passed gen particles
148
149 } // loop over gen events
150
151 ATH_MSG_DEBUG( "Created initial simulation particle collection with size " << simParticles.size() );
152
153 return StatusCode::SUCCESS;
154}
#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 156 of file InputConverter.cxx.

158 {
159 ISF::ISFParticleContainer simParticleList{}; // particles for ISF simulation
160 ATH_CHECK(this->convert(inputGenEvents, simParticleList));
161 //Convert from ISFParticleContainer to ConstISFParticleVector
162 ISF::ISFParticleVector simParticleVector{
163 std::make_move_iterator(std::begin(simParticleList)),
164 std::make_move_iterator(std::end(simParticleList))
165 };
166 HepMC::GenEvent* shadowGenEvent =
167 !shadowGenEvents.empty()
168 ? static_cast<HepMC::GenEvent*>(shadowGenEvents.back())
169 : nullptr;
170 this->ISF_to_G4Event(outputG4Event, simParticleVector, inputGenEvents.back(),
171 shadowGenEvent);
172 return StatusCode::SUCCESS;
173}
#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 175 of file InputConverter.cxx.

176 {
177 ISF::ISFParticleContainer simParticleList{}; // particles for ISF simulation
178 ATH_CHECK(this->convert(inputGenEvents, simParticleList));
179 //Convert from ISFParticleContainer to ConstISFParticleVector
180 ISF::ISFParticleVector simParticleVector{
181 std::make_move_iterator(std::begin(simParticleList)),
182 std::make_move_iterator(std::end(simParticleList))
183 };
184 this->ISF_to_G4Event(outputG4Event, simParticleVector, inputGenEvents.back(),
185 nullptr);
186 return StatusCode::SUCCESS;
187}

◆ 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 245 of file InputConverter.cxx.

245 {
246 if (!genPartPtr) { return nullptr; }
247
248 auto pVertex = genPartPtr->production_vertex();
249 if (!pVertex) {
250 ATH_MSG_ERROR("Unable to convert following generator particle due to missing production vertex for: " << genPartPtr);
251 return nullptr;
252 }
253 auto parentEvent = genPartPtr->parent_event();
254 if (!parentEvent) {
255 ATH_MSG_ERROR("Cannot convert a GenParticle without a parent GenEvent into an ISFParticle!!!");
256 return nullptr;
257 }
258
259 const Amg::Vector3D pos(pVertex->position().x(), pVertex->position().y(), pVertex->position().z());
260 const auto& pMomentum(genPartPtr->momentum());
261 const Amg::Vector3D mom(pMomentum.px(), pMomentum.py(), pMomentum.pz());
262 const double pMass = this->getParticleMass(genPartPtr);
263 double e=pMomentum.e();
264 if (e>1) { //only test for >1 MeV in momentum
265 double px=pMomentum.px();
266 double py=pMomentum.py();
267 double pz=pMomentum.pz();
268 double teste=std::sqrt(px*px + py*py + pz*pz + pMass*pMass);
269 if (std::abs(e-teste)>0.01*e) {
270 ATH_MSG_WARNING("Difference in energy for: " << genPartPtr<<" Morg="<<pMomentum.m()<<" Mmod="<<pMass<<" Eorg="<<e<<" Emod="<<teste);
271 }
272 if (MC::isDecayed(genPartPtr) && pVertex && genPartPtr->end_vertex()) { //check for possible changes of gamma for quasi stable particles
273 const auto& prodVtx = genPartPtr->production_vertex()->position();
274 const auto& endVtx = genPartPtr->end_vertex()->position();
275 CLHEP::Hep3Vector dist3D(endVtx.x()-prodVtx.x(), endVtx.y()-prodVtx.y(), endVtx.z()-prodVtx.z());
276
277 if (dist3D.mag()>1*Gaudi::Units::mm) {
278 CLHEP::HepLorentzVector mom( pMomentum.x(), pMomentum.y(), pMomentum.z(), pMomentum.t() );
279 double gamma_org=mom.gamma();
280 mom.setE(teste);
281 double gamma_new=mom.gamma();
282
283 if (std::abs(gamma_new-gamma_org)/(gamma_new+gamma_org)>0.001) {
284 ATH_MSG_WARNING("Difference in boost gamma for Quasi stable particle "<<genPartPtr);
285 ATH_MSG_WARNING(" gamma(m="<<mom.m()<<")="<<gamma_org<<" gamma(m="<<pMass<<")="<<gamma_new);
286 } else {
287 ATH_MSG_VERBOSE("Quasi stable particle "<<genPartPtr);
288 ATH_MSG_VERBOSE(" gamma(m="<<mom.m()<<")="<<gamma_org<<" gamma(m="<<pMass<<")="<<gamma_new);
289 }
290 }
291 }
292 }
293
294 const int pPdgId = genPartPtr->pdg_id();
295 const double charge = MC::charge(pPdgId);
296 const double pTime = pVertex->position().t() / Gaudi::Units::c_light;
299 const auto pBarcode = HepMC::barcode(genPartPtr);
300 const auto particleID = HepMC::uniqueID(genPartPtr);
301 auto tBinding = std::make_unique<ISF::TruthBinding>(genPartPtr);
302
303 auto hmpl = std::make_unique<HepMcParticleLink>(particleID, parentEvent->event_number(), HepMcParticleLink::IS_EVENTNUM, HepMcParticleLink::IS_ID);
304
305 auto sParticle = std::make_unique<ISF::ISFParticle>( std::move(pos),
306 std::move(mom),
307 pMass,
308 charge,
309 pPdgId,
310 genPartPtr->status(),
311 pTime,
312 origin,
313 particleID,
314 pBarcode,
315 tBinding.release(),
316 hmpl.release() );
317 return sParticle.release();
318}
#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 108 of file InputConverter.cxx.

109{
110 ATH_MSG_DEBUG("Finalizing ...");
111 return StatusCode::SUCCESS;
112}

◆ findShadowParticle()

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

Definition at line 625 of file InputConverter.cxx.

626{
627 if (!shadowGenEvent) {
628 ATH_MSG_FATAL ("Found status==2 GenParticle with no end vertex and shadow GenEvent is missing - something is wrong here!");
629 abort();
630 }
631 // TODO in the future switch to using an Attribute which stores the shadow GenParticlePtr directly.
632 const int shadowId = genParticle->attribute<HepMC3::IntAttribute>(HepMCStr::ShadowParticleId)->value();
633 for (auto& shadowParticle : shadowGenEvent->particles()) {
634 if (shadowParticle->id() == shadowId && matchedGenParticles(genParticle, shadowParticle) ) { return shadowParticle; }
635 }
636 return std::make_shared<HepMC::GenParticle>();
637}
#define ATH_MSG_FATAL(x)
bool matchedGenParticles(const HepMC::ConstGenParticlePtr &p1, const HepMC::ConstGenParticlePtr &p2) const
const std::string ShadowParticleId

◆ getDaughterG4PrimaryParticle() [1/2]

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

Definition at line 535 of file InputConverter.cxx.

535 {
536 ATH_MSG_VERBOSE("Creating G4PrimaryParticle from GenParticle.");
537
538 const G4ParticleDefinition *particleDefinition = this->getG4ParticleDefinition(genpart->pdg_id());
539
540 if (particleDefinition==nullptr) {
541 ATH_MSG_ERROR("ISF_to_G4Event particle conversion failed. ISF_Particle PDG code = " << genpart->pdg_id() <<
542 "\n This usually indicates a problem with the evgen step.\n" <<
543 "Please report this to the Generators group, mentioning the release and generator used for evgen and the PDG code above." );
544 return nullptr;
545 }
546
547 // create new primaries and set them to the vertex
548 // G4double mass = particleDefinition->GetPDGMass();
549 auto &genpartMomentum = genpart->momentum();
550 G4double px = genpartMomentum.x();
551 G4double py = genpartMomentum.y();
552 G4double pz = genpartMomentum.z();
553
554 std::unique_ptr<G4PrimaryParticle> g4particle = std::make_unique<G4PrimaryParticle>(particleDefinition,px,py,pz);
555
556 if (genpart->end_vertex()) {
557 // Set the lifetime appropriately - this is slow but rigorous, and we
558 // don't want to end up with something like vertex time that we have
559 // to validate for every generator on earth...
560 const auto& prodVtx = genpart->production_vertex()->position();
561 const auto& endVtx = genpart->end_vertex()->position();
562 //const G4LorentzVector lv0 ( prodVtx.x(), prodVtx.y(), prodVtx.z(), prodVtx.t() );
563 //const G4LorentzVector lv1 ( endVtx.x(), endVtx.y(), endVtx.z(), endVtx.t() );
564 //Old calculation, not taken because vertex information is not sufficiently precise
565 //g4particle->SetProperTime( (lv1-lv0).mag()/Gaudi::Units::c_light );
566
567 CLHEP::Hep3Vector dist3D(endVtx.x()-prodVtx.x(), endVtx.y()-prodVtx.y(), endVtx.z()-prodVtx.z());
568 double tau=SetProperTimeFromDetectorFrameDecayLength(*g4particle,dist3D.mag());
569
570 if (msgLvl(MSG::VERBOSE)) {
571 double pmag2=g4particle->GetTotalMomentum(); //magnitude of particle momentum
572 pmag2*=pmag2; //magnitude of particle momentum squared
573 double e2=g4particle->GetTotalEnergy(); //energy of particle
574 e2*=e2; //energy of particle squared
575 double beta2=pmag2/e2; //beta^2=v^2/c^2 for particle
576 double tau2=dist3D.mag2()*(1/beta2-1)/Gaudi::Units::c_light/Gaudi::Units::c_light;
577 ATH_MSG_VERBOSE("lifetime tau(beta)="<<std::sqrt(tau2)<<" tau="<<tau);
578 }
580 ATH_MSG_VERBOSE( "Detected primary particle with end vertex." );
581 ATH_MSG_VERBOSE( "Will add the primary particle set on." );
582 ATH_MSG_VERBOSE( "Primary Particle: " << genpart );
583 ATH_MSG_VERBOSE( "Number of daughters : " << genpart->end_vertex()->particles_out().size()<<" at position "<<genpart->end_vertex() );
584 }
585 else {
586 ATH_MSG_WARNING( "Detected primary particle with end vertex." );
587 ATH_MSG_WARNING( "Will add the primary particle set on." );
588 ATH_MSG_WARNING( "Primary Particle: " << genpart );
589 ATH_MSG_WARNING( "Number of daughters: " << genpart->end_vertex()->particles_out().size()<<" at position "<<genpart->end_vertex() );
590 }
591 // Add all necessary daughter particles
592 for ( auto daughter: genpart->end_vertex()->particles_out() ) {
594 ATH_MSG_VERBOSE ( "Attempting to add daughter particle: " << daughter );
595 }
596 else {
597 ATH_MSG_WARNING ( "Attempting to add daughter particle: " << daughter );
598 }
599 G4PrimaryParticle *daughterG4Particle = this->getDaughterG4PrimaryParticle( daughter );
600 if (!daughterG4Particle) {
601 ATH_MSG_ERROR("Bailing out of loop over daughters of particle due to errors - will not return G4Particle.");
602 return nullptr;
603 }
604 g4particle->SetDaughter( daughterG4Particle );
605 }
606 }
607
608 return g4particle.release();
609}
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 449 of file InputConverter.cxx.

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

◆ getG4ParticleDefinition()

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

Special cases for Geantinos

Standard particles

Definition at line 415 of file InputConverter.cxx.

416{
418 if (pdgcode==998) {
419 return G4ChargedGeantino::Definition();
420 }
421 if (pdgcode==999) {
422 return G4Geantino::GeantinoDefinition();
423 }
425 G4ParticleTable *ptable = G4ParticleTable::GetParticleTable();
426 if (ptable) {
427 return ptable->FindParticle(pdgcode);
428 }
429 ATH_MSG_ERROR("getG4ParticleDefinition - Failed to retrieve G4ParticleTable!");
430 return nullptr;
431}

◆ 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 797 of file InputConverter.cxx.

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

323 {
324 // default value: generated particle mass
325 double mass = part->generated_mass();
326 ATH_MSG_VERBOSE("part->generated_mass, mass="<<mass);
327
328 // 1. use PDT mass?
330 const int absPDG = std::abs(part->pdg_id());
331 HepPDT::ParticleData const *pData = (m_particleDataTable)
332 ? m_particleDataTable->particle(absPDG)
333 : nullptr;
334 if (pData) {
335 mass = pData->mass();
336 ATH_MSG_VERBOSE("using pData mass, mass="<<mass);
337 }
338 else {
339 ATH_MSG_WARNING( "Unable to find mass of particle with PDG ID '" << absPDG << "' in ParticleDataTable. Will set mass to generated_mass: " << mass);
340 }
341 }
342 return mass;
343}

◆ 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 191 of file InputConverter.cxx.

191 {
192 auto allGenPartBegin = evnt.particles().begin();
193 auto allGenPartEnd = evnt.particles().end();
194
195 // reserve destination container with maximum size, i.e. number of particles in input event
196 std::vector<HepMC::GenParticlePtr> passedGenParticles{};
197 size_t maxParticles = std::distance(allGenPartBegin, allGenPartEnd);
198 passedGenParticles.reserve(maxParticles);
199
200 if (m_useShadowEvent) {
201 if (legacyOrdering) {
202 // FIXME: remove this block and the 'legacyOrdering' flag
203 // once we don't need the legacy order any longer
204 for (auto vtx: evnt.vertices() ) {
205 std::copy_if (vtx->particles_out().begin(),
206 vtx->particles_out().end(),
207 std::back_inserter(passedGenParticles),
208 [](HepMC::GenParticlePtr p){return p->attribute<HepMC3::IntAttribute>(HepMCStr::ShadowParticleId);});
209 }
210 }
211 else {
212 std::copy_if (allGenPartBegin,
213 allGenPartEnd,
214 std::back_inserter(passedGenParticles),
215 [](HepMC::GenParticlePtr p){return p->attribute<HepMC3::IntAttribute>(HepMCStr::ShadowParticleId);});
216 }
217 }
218 else {
219 if (legacyOrdering) {
220 // FIXME: remove this block and the 'legacyOrdering' flag
221 // once we don't need the legacy order any longer
222 for (auto vtx: evnt.vertices() ) {
223 std::copy_if (vtx->particles_out().begin(),
224 vtx->particles_out().end(),
225 std::back_inserter(passedGenParticles),
226 [this](HepMC::GenParticlePtr p){return this->passesFilters(std::const_pointer_cast<const HepMC3::GenParticle>(p));});
227 }
228 }
229 else {
230 std::copy_if (allGenPartBegin,
231 allGenPartEnd,
232 std::back_inserter(passedGenParticles),
233 [this](HepMC::GenParticlePtr p){return this->passesFilters(std::const_pointer_cast<const HepMC3::GenParticle>(p));});
234 }
235 }
236
237 passedGenParticles.shrink_to_fit();
238
239 return passedGenParticles;
240}
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 83 of file InputConverter.cxx.

84{
85 ATH_MSG_VERBOSE("initialize() begin");
86
87 // setup PDT if requested (to get particle masses later on)
89 ATH_CHECK(m_particlePropSvc.retrieve());
92 ATH_MSG_FATAL( "Could not get ParticleDataTable from " << m_particlePropSvc << ". Abort" );
93 return StatusCode::FAILURE;
94 }
95 }
96
97 if (!m_genParticleFilters.empty()) {
99 }
100
101 ATH_MSG_VERBOSE("initialize() successful");
102 return StatusCode::SUCCESS;
103}

◆ 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 376 of file InputConverter.cxx.

379 {
380 // G4Event *g4evt = new G4Event(ctx.eventID().event_number());
381
382 // retrieve world solid (volume)
383 const G4VSolid *worldSolid = G4TransportationManager::GetTransportationManager()->GetNavigatorForTracking()->GetWorldVolume()->GetLogicalVolume()->GetSolid();
384
385 for ( ISF::ISFParticle *ispPtr: ispVector ) {
386 ISF::ISFParticle &isp = *ispPtr;
387 if ( !isInsideG4WorldVolume(isp, worldSolid) ) {
388 ATH_MSG_WARNING("Unable to convert ISFParticle to G4PrimaryParticle!");
389 ATH_MSG_WARNING(" ISFParticle: " << isp );
390 if (worldSolid) {
391 ATH_MSG_WARNING(" is outside Geant4 world volume: ");
392 worldSolid->DumpInfo();
393 G4cout << std::flush;
394 }
395 else {
396 ATH_MSG_WARNING(" is outside Geant4 world volume.");
397 }
398 continue;
399 }
400 this->addG4PrimaryVertex(event, isp, useHepMC, shadowGenEvent);
401 }
402
403 AtlasG4EventUserInfo* atlasG4EvtUserInfo =
404 dynamic_cast<AtlasG4EventUserInfo*>(event.GetUserInformation());
405 if (!atlasG4EvtUserInfo) {
406 atlasG4EvtUserInfo = new AtlasG4EventUserInfo(Gaudi::Hive::currentContext());
407 event.SetUserInformation(atlasG4EvtUserInfo);
408 }
409 atlasG4EvtUserInfo->SetLastProcessedTrackID(0); // TODO Check if it is better to set this to -1 initially
410 atlasG4EvtUserInfo->SetLastProcessedStep(0); // TODO Check if it is better to set this to -1 initially
411 atlasG4EvtUserInfo->SetHepMCEvent(genEvent);
412}
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 970 of file InputConverter.cxx.

971{
972
973 const Amg::Vector3D &pos = isp.position();
974 const G4ThreeVector g4Pos( pos.x(), pos.y(), pos.z() );
975 EInside insideStatus = worldSolid->Inside( g4Pos );
976
977 bool insideWorld = insideStatus != kOutside;
978 return insideWorld;
979}

◆ matchedGenParticles()

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

Definition at line 612 of file InputConverter.cxx.

614{
615 return (HepMC::barcode(p1) == HepMC::barcode(p2))
616 && (p1->status() == p2->status())
617 && (p1->pdg_id() == p2->pdg_id())
618 && ((p1->momentum().px()) == (p2->momentum().px()))
619 && ((p1->momentum().py()) == (p2->momentum().py()))
620 && ((p1->momentum().pz()) == (p2->momentum().pz()))
621 && (float(p1->momentum().m()) == float(p2->momentum().m()));
622}

◆ passesFilters()

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

check if the given particle passes all filters

Definition at line 349 of file InputConverter.cxx.

350{
351 // TODO: implement this as a std::find_if with a lambda function
352 for ( const auto& filter : m_genParticleFilters ) {
353 // determine if the particle passes current filter
354 bool passFilter = filter->pass(part);
355 ATH_MSG_VERBOSE("GenParticleFilter '" << filter.typeAndName() << "' returned: "
356 << (passFilter ? "true, will keep particle."
357 : "false, will remove particle."));
358 const auto& momentum = part->momentum();
359 ATH_MSG_VERBOSE("Particle: ("
360 <<momentum.px()<<", "
361 <<momentum.py()<<", "
362 <<momentum.pz()<<"), pdgCode: "
363 <<part->pdg_id() );
364
365 if (!passFilter) {
366 return false;
367 }
368 }
369
370 return true;
371}

◆ 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 640 of file InputConverter.cxx.

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

◆ 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 719 of file InputConverter.cxx.

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

Member Data Documentation

◆ m_genParticleFilters

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

HepMC::GenParticle filters.

Definition at line 123 of file InputConverter.h.

◆ m_particleDataTable

const HepPDT::ParticleDataTable* ISF::InputConverter::m_particleDataTable
private

PDT used to look up particle masses.

Definition at line 119 of file InputConverter.h.

◆ m_particlePropSvc

ServiceHandle<IPartPropSvc> ISF::InputConverter::m_particlePropSvc
private

ParticlePropertyService and ParticleDataTable.

particle properties svc to retrieve PDT

Definition at line 118 of file InputConverter.h.

◆ m_quasiStableParticlesIncluded

bool ISF::InputConverter::m_quasiStableParticlesIncluded
private

Definition at line 125 of file InputConverter.h.

◆ m_useGeneratedParticleMass

bool ISF::InputConverter::m_useGeneratedParticleMass
private

use GenParticle::generated_mass() in simulation

Definition at line 121 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 126 of file InputConverter.h.

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

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