DDG4: allow configuring of the unstable generator status codes.

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; k<num; ++k)
       p->parents.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

DDG4/include/DDG4/Geant4InputAction.h

@@ -31,8 +31,9 @@
 #include "Parsers/Parsers.h"
 
 // C/C++ include files
-#include <vector>
 #include <memory>
+#include <set>
+#include <vector>
 
 // 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<int> 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<int>;
+      /// 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;
 

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; k<num; ++k)
       p->parents.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

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)

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

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;
+}