From 19f0822af528934d9e6358607a806e1b9c4d1b29 Mon Sep 17 00:00:00 2001
From: Andre Sailer <andre.philippe.sailer@cern.ch>
Date: Tue, 7 May 2024 15:25:50 +0200
Subject: [PATCH] DDG4: allow configuring of the unstable generator status
codes.
Move generator status setting to common place
---
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 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; 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
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 <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;
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; 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
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 @@ class DD4hepSimulation(object):
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 @@ class Physics(ConfigHelper):
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 @@ class Physics(ConfigHelper):
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;
+}
--
GitLab