From 4e4e23de0a9fce3642de3e7cc9d62db08f767d02 Mon Sep 17 00:00:00 2001 From: Ben Rosser Date: Wed, 22 May 2024 01:14:37 +0200 Subject: [PATCH] Backport #1260 (allow configuring of the unstable generator status codes) See this merge request from Andre: https://patch-diff.githubusercontent.com/raw/AIDASoft/DD4hep/pull/1260.patch The diff didn't apply perfectly cleanly. I applied it manually using `patch -p1` rather than through git. But this is just a check to see if we can bring this functionality back to an older version. --- DDG4/hepmc/HepMC3EventReader.cpp | 10 ++-------- DDG4/include/DDG4/Geant4InputAction.h | 10 +++++++++- DDG4/lcio/LCIOEventReader.cpp | 10 ++-------- DDG4/python/DDSim/DD4hepSimulation.py | 1 + DDG4/python/DDSim/Helper/Physics.py | 11 +++++++++++ DDG4/src/Geant4InputAction.cpp | 14 ++++++++++++++ 6 files changed, 39 insertions(+), 17 deletions(-) diff --git a/DDG4/hepmc/HepMC3EventReader.cpp b/DDG4/hepmc/HepMC3EventReader.cpp index 74425aff3..5277c4b0a 100644 --- a/DDG4/hepmc/HepMC3EventReader.cpp +++ b/DDG4/hepmc/HepMC3EventReader.cpp @@ -95,7 +95,6 @@ HEPMC3EventReader::readParticles(int event_number, Vertices& vertices, Particles const float spin[] = {0.0, 0.0, 0.0}; //mcp->getSpin(); // FIXME const int color[] = {colorFlow[0][mcp->id()], colorFlow[1][mcp->id()]}; const int pdg = mcp->pid(); - PropertyMask status(p->status); p->pdgID = pdg; p->charge = 0; // int(mcp->getCharge()*3.0); // FIXME p->psx = mom.get_component(0) * mom_unit; @@ -122,16 +121,11 @@ HEPMC3EventReader::readParticles(int event_number, Vertices& vertices, Particles for(int num=par.size(), k=0; kparents.insert(GET_ENTRY(mcparts,par[k])); + PropertyMask status(p->status); int genStatus = mcp->status(); - if ( genStatus == 0 ) status.set(G4PARTICLE_GEN_EMPTY); - else if ( genStatus == 1 ) status.set(G4PARTICLE_GEN_STABLE); - else if ( genStatus == 2 ) status.set(G4PARTICLE_GEN_DECAYED); - else if ( genStatus == 3 ) status.set(G4PARTICLE_GEN_DOCUMENTATION); - else if ( genStatus == 4 ) status.set(G4PARTICLE_GEN_BEAM); - else - status.set(G4PARTICLE_GEN_OTHER); // Copy raw generator status p->genStatus = genStatus&G4PARTICLE_GEN_STATUS_MASK; + m_inputAction->setGeneratorStatus(genStatus, status); if ( p->parents.size() == 0 ) { // A particle without a parent in HepMC3 can only be (something like) a beam particle, and it is attached to the diff --git a/DDG4/include/DDG4/Geant4InputAction.h b/DDG4/include/DDG4/Geant4InputAction.h index e6b07a045..72cd9364b 100644 --- a/DDG4/include/DDG4/Geant4InputAction.h +++ b/DDG4/include/DDG4/Geant4InputAction.h @@ -31,8 +31,9 @@ #include "Parsers/Parsers.h" // C/C++ include files -#include #include +#include +#include // Forward declarations class G4Event; @@ -174,11 +175,18 @@ namespace dd4hep { /// Property: named parameters to configure file readers or input actions std::map< std::string, std::string> m_parameters; + /// Property: set of alternative decay statuses that MC generators might use for unstable particles + std::set m_alternativeDecayStatuses = {}; + public: /// Read an event and return a LCCollectionVec of MCParticles. int readParticles(int event_number, Vertices& vertices, Particles& particles); + using PropertyMask = dd4hep::detail::ReferenceBitMask; + /// Convert the generator status into a common set of generator status bits + void setGeneratorStatus(int generatorStatus, PropertyMask& status); + /// helper to report Geant4 exceptions std::string issue(int i) const; diff --git a/DDG4/lcio/LCIOEventReader.cpp b/DDG4/lcio/LCIOEventReader.cpp index 7f6383409..4903cd263 100644 --- a/DDG4/lcio/LCIOEventReader.cpp +++ b/DDG4/lcio/LCIOEventReader.cpp @@ -91,7 +91,6 @@ LCIOEventReader::readParticles(int event_number, const float *spin = mcp->getSpin(); const int *color = mcp->getColorFlow(); const int pdg = mcp->getPDG(); - PropertyMask status(p->status); p->pdgID = pdg; p->charge = int(mcp->getCharge()*3.0); p->psx = mom[0]*CLHEP::GeV; @@ -118,16 +117,11 @@ LCIOEventReader::readParticles(int event_number, for(int num=par.size(),k=0; kparents.insert(GET_ENTRY(mcparts,par[k])); + PropertyMask status(p->status); int genStatus = mcp->getGeneratorStatus(); - if ( genStatus == 0 ) status.set(G4PARTICLE_GEN_EMPTY); - else if ( genStatus == 1 ) status.set(G4PARTICLE_GEN_STABLE); - else if ( genStatus == 2 ) status.set(G4PARTICLE_GEN_DECAYED); - else if ( genStatus == 3 ) status.set(G4PARTICLE_GEN_DOCUMENTATION); - else if ( genStatus == 4 ) status.set(G4PARTICLE_GEN_BEAM); - else - status.set(G4PARTICLE_GEN_OTHER); // Copy raw generator status p->genStatus = genStatus&G4PARTICLE_GEN_STATUS_MASK; + m_inputAction->setGeneratorStatus(genStatus, status); //fg: we simply add all particles without parents as with their own vertex. // This might include the incoming beam particles, e.g. in diff --git a/DDG4/python/DDSim/DD4hepSimulation.py b/DDG4/python/DDSim/DD4hepSimulation.py index da83bd6a8..95d1e4405 100644 --- a/DDG4/python/DDSim/DD4hepSimulation.py +++ b/DDG4/python/DDSim/DD4hepSimulation.py @@ -418,6 +418,7 @@ def run(self): else: # this should never happen because we already check at the top, but in case of some LogicError... raise RuntimeError("Unknown input file type: %s" % inputFile) + gen.AlternativeDecayStatuses = self.physics.alternativeDecayStatuses gen.Sync = self.skipNEvents gen.Mask = index actionList.append(gen) diff --git a/DDG4/python/DDSim/Helper/Physics.py b/DDG4/python/DDSim/Helper/Physics.py index db1e82e33..bc703a8dd 100644 --- a/DDG4/python/DDSim/Helper/Physics.py +++ b/DDG4/python/DDSim/Helper/Physics.py @@ -28,6 +28,7 @@ def __init__(self): 4101, 4103, 4201, 4203, 4301, 4303, 4403, # c? diquarks 5101, 5103, 5201, 5203, 5301, 5303, 5401, 5403, 5503} # b? diquarks self._zeroTimePDGs = {11, 13, 15, 17} + self._alternativeDecayStatuses = set() self._userFunctions = [] @property @@ -54,6 +55,16 @@ def zeroTimePDGs(self): def zeroTimePDGs(self, val): self._zeroTimePDGs = self.makeSet(val) + @property + def alternativeDecayStatuses(self): + """Set of Generator Statuses that are used to mark unstable particles that should decay inside of Geant4. + """ + return self._alternativeDecayStatuses + + @alternativeDecayStatuses.setter + def alternativeDecayStatuses(self, val): + self._alternativeDecayStatuses = self.makeSet(val) + @property def rangecut(self): """ The global geant4 rangecut for secondary production diff --git a/DDG4/src/Geant4InputAction.cpp b/DDG4/src/Geant4InputAction.cpp index 999962518..0c41934f2 100644 --- a/DDG4/src/Geant4InputAction.cpp +++ b/DDG4/src/Geant4InputAction.cpp @@ -126,6 +126,7 @@ Geant4InputAction::Geant4InputAction(Geant4Context* ctxt, const string& nam) declareProperty("MomentumScale", m_momScale = 1.0); declareProperty("HaveAbort", m_abort = true); declareProperty("Parameters", m_parameters = {}); + declareProperty("AlternativeDecayStatuses", m_alternativeDecayStatuses = {}); m_needsControl = true; } @@ -274,3 +275,16 @@ void Geant4InputAction::operator()(G4Event* event) { p.dumpWithMomentumAndVertex(outputLevel()-1,name(),"->"); } } + +void Geant4InputAction::setGeneratorStatus(int genStatus, PropertyMask& status) { + if ( genStatus == 0 ) status.set(G4PARTICLE_GEN_EMPTY); + else if ( genStatus == 1 ) status.set(G4PARTICLE_GEN_STABLE); + else if ( genStatus == 2 ) status.set(G4PARTICLE_GEN_DECAYED); + else if ( genStatus == 3 ) status.set(G4PARTICLE_GEN_DOCUMENTATION); + else if ( genStatus == 4 ) status.set(G4PARTICLE_GEN_BEAM); + else if ( m_alternativeDecayStatuses.count(genStatus) ) status.set(G4PARTICLE_GEN_DECAYED); + else + status.set(G4PARTICLE_GEN_OTHER); + + return; +}