diff --git a/DDG4/hepmc/HepMC3EventReader.cpp b/DDG4/hepmc/HepMC3EventReader.cpp index 88c1124e5..bf9f46094 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 b8a0a5bc4..8da8ad550 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; @@ -178,6 +179,9 @@ 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 = {}; + /// Perform some actions before the run starts, like opening the event inputs void beginRun(const G4Run*); @@ -188,6 +192,10 @@ namespace dd4hep { 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 75f6e4289..bad047fcd 100644 --- a/DDG4/python/DDSim/DD4hepSimulation.py +++ b/DDG4/python/DDSim/DD4hepSimulation.py @@ -440,6 +440,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 ac92bd52d..6a3d93754 100644 --- a/DDG4/python/DDSim/Helper/Physics.py +++ b/DDG4/python/DDSim/Helper/Physics.py @@ -26,6 +26,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 = [] self._closeProperties() @@ -53,6 +54,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 04917c45c..07ed4e0a9 100644 --- a/DDG4/src/Geant4InputAction.cpp +++ b/DDG4/src/Geant4InputAction.cpp @@ -126,6 +126,7 @@ Geant4InputAction::Geant4InputAction(Geant4Context* ctxt, const std::string& nam declareProperty("MomentumScale", m_momScale = 1.0); declareProperty("HaveAbort", m_abort = true); declareProperty("Parameters", m_parameters = {}); + declareProperty("AlternativeDecayStatuses", m_alternativeDecayStatuses = {}); m_needsControl = true; runAction().callAtBegin(this, &Geant4InputAction::beginRun); @@ -285,3 +286,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; +}