diff --git a/DDCore/include/DD4hep/Primitives.h b/DDCore/include/DD4hep/Primitives.h
index 1fa133ad28b34da6ac16ce67c0d80319c52ed511..f68088aec49b06b72d7cf83fc3e1bf15170ed0fe 100644
--- a/DDCore/include/DD4hep/Primitives.h
+++ b/DDCore/include/DD4hep/Primitives.h
@@ -217,6 +217,39 @@ namespace DD4hep {
template <typename M> ReleaseObjects<M> releaseObjects(M& m) {
return ReleaseObjects<M>(m);
}
+
+
+ /// Data structure to manipulate a bitmask held by reference and represented by an integer
+ /**
+ * @author M.Frank
+ * @version 1.0
+ */
+ template <typename T> class ReferenceBitMask {
+ public:
+ /// Reference to the data
+ T& mask;
+ /// Standard constructor
+ ReferenceBitMask(T& m);
+ T value() const {
+ return mask;
+ }
+ void set(const T& m) {
+ mask |= m;
+ }
+ bool isSet(const T& m) {
+ return (mask&m) == m;
+ }
+ bool testBit(int bit) const {
+ T m = T(1)<<bit;
+ return isSet(m);
+ }
+ bool isNull() const {
+ return mask == 0;
+ }
+ };
+ /// Standard constructor
+ template <typename T> ReferenceBitMask<T>::ReferenceBitMask(T& m) : mask(m) {}
+
/*
* Geometry namespace declaration
*/
diff --git a/DDG4/examples/CLICSidSimu.py b/DDG4/examples/CLICSidSimu.py
index b16fb9ee56a34c68f3f7287ecc7d13ecd23b960c..b75a606496b9b30c0c6f3454f4f2b66d1672cc67 100644
--- a/DDG4/examples/CLICSidSimu.py
+++ b/DDG4/examples/CLICSidSimu.py
@@ -52,6 +52,10 @@ def run():
kernel.eventAction().add(evt1)
kernel.eventAction().add(evt2)
+ prt = DDG4.EventAction(kernel,'Geant4ParticlePrint/ParticlePrint')
+ prt.OutputType = 3 # Print both: table and tree
+ kernel.eventAction().add(prt)
+
"""
mc = DDG4.Action(kernel,"Geant4MonteCarloRecordManager/MonteCarloRecordManager")
kernel.registerGlobalAction(mc)
@@ -65,14 +69,12 @@ def run():
kernel.generatorAction().add(gen)
# Setup particle gun
- gun = simple.setupGun('Gun','pi-',100*GeV,True)
+ gun = simple.setupGun('Gun','pi-',energy=10*GeV,isotrop=True,multiplicity=3)
trk = DDG4.GeneratorAction(kernel,"Geant4ParticleHandler/ParticleHandler")
kernel.generatorAction().add(trk)
- trk.SaveProcesses = ['conv','Decay']
+ trk.saveProcesses = ['conv','Decay']
trk.enableUI()
- """
- """
"""
rdr = DDG4.GeneratorAction(kernel,"LcioGeneratorAction/Reader")
diff --git a/DDG4/include/DDG4/Geant4Data.h b/DDG4/include/DDG4/Geant4Data.h
index 410358033a97b0b6091f6321001dd26d1e2fe8c1..9cefc8c9978d1fc083eefc52a6579560f05e7e34 100644
--- a/DDG4/include/DDG4/Geant4Data.h
+++ b/DDG4/include/DDG4/Geant4Data.h
@@ -97,6 +97,43 @@ namespace DD4hep {
virtual ~SimpleEvent();
};
+ /// Track properties
+ enum ParticleProperties {
+ G4PARTICLE_CREATED_HIT = 1<<0,
+ G4PARTICLE_PRIMARY = 1<<1,
+ G4PARTICLE_HAS_SECONDARIES = 1<<2,
+ G4PARTICLE_ABOVE_ENERGY_THRESHOLD = 1<<3,
+ G4PARTICLE_KEEP_PROCESS = 1<<4,
+ G4PARTICLE_KEEP_PARENT = 1<<5,
+ G4PARTICLE_CREATED_CALORIMETER_HIT = 1<<6,
+ G4PARTICLE_CREATED_TRACKER_HIT = 1<<7,
+ G4PARTICLE_LAST_NOTHING = 1<<31
+ };
+
+ /// Data structure to store the MC particle information
+ /**
+ * @author M.Frank
+ * @version 1.0
+ */
+ class Particle {
+ public:
+ int id, g4Parent, parent, reason, steps, secondaries;
+ double vsx, vsy, vsz;
+ double vex, vey, vez;
+ double psx, psy, psz, pex, pey, pez, energy, time;
+ const G4VProcess *process;
+ const G4ParticleDefinition *definition;
+ std::set<int> daughters;
+ /// Default constructor
+ Particle();
+ /// Copy constructor
+ Particle(const Particle& c);
+ /// Default destructor
+ virtual ~Particle();
+ /// Remove daughter from set
+ void removeDaughter(int id_daughter);
+ };
+
/** @class SimpleHit Geant4Data.h DDG4/Geant4Data.h
*
* Base class for geant4 hit structures created by the
diff --git a/DDG4/include/DDG4/Geant4MonteCarloTruth.h b/DDG4/include/DDG4/Geant4MonteCarloTruth.h
index bd1d90f2eecf3e8581022c201bc134c93ce6ea11..3db9049faf06287d650803712dc99b22ede1a569 100644
--- a/DDG4/include/DDG4/Geant4MonteCarloTruth.h
+++ b/DDG4/include/DDG4/Geant4MonteCarloTruth.h
@@ -13,6 +13,7 @@
#include "DDG4/Geant4Action.h"
// C/C++ include files
+#include <map>
// Forward declarations
class G4Step;
@@ -29,6 +30,9 @@ namespace DD4hep {
*/
namespace Simulation {
+ // Forward declarations
+ class Particle;
+
/** @class Geant4MonteCarloTruth Geant4MonteCarloTruth.h DDG4/Geant4MonteCarloTruth.h
*
* Default Interface class to handle monte carlo truth records
@@ -37,20 +41,27 @@ namespace DD4hep {
* @version 1.0
*/
class Geant4MonteCarloTruth {
+ public:
+ typedef std::map<int,Particle*> ParticleMap;
+ typedef std::map<int,int> TrackEquivalents;
protected:
/// Standard constructor
Geant4MonteCarloTruth();
public:
/// Default destructor
virtual ~Geant4MonteCarloTruth();
+ /// Access the particle map
+ virtual const ParticleMap& particles() const = 0;
+ /// Access the map of track equivalents
+ virtual const TrackEquivalents& equivalents() const = 0;
/// Mark a Geant4 track to be kept for later MC truth analysis
virtual void mark(const G4Track* track) = 0;
/// Store a track, with a flag
- virtual void mark(const G4Track* track, bool created_hit) = 0;
+ virtual void mark(const G4Track* track, int reason) = 0;
/// Mark a Geant4 track of the step to be kept for later MC truth analysis
virtual void mark(const G4Step* step) = 0;
/// Store a track produced in a step to be kept for later MC truth analysis
- virtual void mark(const G4Step* step, bool created_hit) = 0;
+ virtual void mark(const G4Step* step, int reason) = 0;
};
/** @class Geant4DummyTruthHandler Geant4DummyTruthHandler.h DDG4/Geant4DummyTruthHandler.h
@@ -63,19 +74,28 @@ namespace DD4hep {
class Geant4DummyTruthHandler :
public Geant4Action, public Geant4MonteCarloTruth
{
+ protected:
+ /// Map with stored MC Particles
+ ParticleMap m_particleMap;
+ /// Map associating the G4Track identifiers with identifiers of existing MCParticles
+ TrackEquivalents m_equivalentTracks;
public:
/// Standard constructor
Geant4DummyTruthHandler(Geant4Context* ctxt,const std::string& nam);
/// Default destructor
virtual ~Geant4DummyTruthHandler();
- /// Mark a Geant4 track to be kept for later MC truth analysis
+ /// Access the particle map
+ virtual const ParticleMap& particles() const { return m_particleMap; }
+ /// Access the map of track equivalents
+ virtual const TrackEquivalents& equivalents() const { return m_equivalentTracks; }
+ /// Mark a Geant4 track to be kept for later MC truth analysis. Default flag: CREATED_HIT
virtual void mark(const G4Track* track);
/// Store a track, with a flag
- virtual void mark(const G4Track* track, bool created_hit);
- /// Mark a Geant4 track of the step to be kept for later MC truth analysis
+ virtual void mark(const G4Track* track, int reason);
+ /// Mark a Geant4 track of the step to be kept for later MC truth analysis. Default flag: CREATED_HIT
virtual void mark(const G4Step* step);
/// Store a track produced in a step to be kept for later MC truth analysis
- virtual void mark(const G4Step* step, bool created_hit);
+ virtual void mark(const G4Step* step, int reason);
};
} // End namespace Simulation
diff --git a/DDG4/include/DDG4/Geant4ParticleHandler.h b/DDG4/include/DDG4/Geant4ParticleHandler.h
index 6f384e25abea41808f7ad29f95f6e157a9cd224e..a31632779619165a0facb62be8d36cc5acbb7533 100644
--- a/DDG4/include/DDG4/Geant4ParticleHandler.h
+++ b/DDG4/include/DDG4/Geant4ParticleHandler.h
@@ -10,26 +10,15 @@
#define DD4HEP_DDG4_GEANT4PARTICLEHANDLER_H
// Framework include files
-#include "DDG4/Geant4Action.h"
-#include "DDG4/Geant4EventAction.h"
+#include "DDG4/Geant4Data.h"
#include "DDG4/Geant4GeneratorAction.h"
#include "DDG4/Geant4MonteCarloTruth.h"
-#include "Math/PxPyPzE4D.h"
-#include "G4VUserTrackInformation.hh"
-
-// C/C++ include files
-#include <map>
-#include <set>
-#include <string>
-#include <typeinfo>
// Forward declarations
class G4Step;
class G4Track;
class G4Event;
-class G4VProcess;
class G4SteppingManager;
-class G4ParticleDefinition;
/*
* DD4hep namespace declaration
@@ -42,88 +31,76 @@ namespace DD4hep {
namespace Simulation {
// Forward declarations
+ class Particle;
- /** @class Geant4ParticleHandler Geant4ParticleHandler.h DDG4/Geant4ParticleHandler.h
+ /** Geant4Action to collect the MC particle information.
*
- * Default base class for all geant 4 actions and derivates thereof.
*
* @author M.Frank
* @version 1.0
*/
class Geant4ParticleHandler : public Geant4GeneratorAction, public Geant4MonteCarloTruth
- {
- public:
- class Particle {
- public:
- int id, parent, theParent, reason;
- double vx, vy, vz;
- double px, py, pz, energy, time;
- const G4VProcess *process;
- const G4ParticleDefinition *definition;
- std::set<int> daughters;
- /// Default constructor
- Particle();
- /// Copy constructor
- Particle(const Particle& c);
- /// Default destructor
- virtual ~Particle();
- /// Remove daughter from set
- void removeDaughter(int id_daughter);
- };
-
- typedef std::map<int,Particle*> ParticleMap;
- typedef std::map<int,int> TrackEquivalents;
- typedef std::vector<std::string> Processes;
- protected:
+ {
+ public:
+ typedef std::vector<std::string> Processes;
+ protected:
- /// Property: Steer printout at tracking action begin
- bool m_printStartTracking;
- /// Property: Steer printout at tracking action end
- bool m_printEndTracking;
-
- Particle m_currTrack;
- ParticleMap m_particleMap;
- TrackEquivalents m_equivalentTracks;
- Processes m_processNames;
+ /// Property: Steer printout at tracking action begin
+ bool m_printStartTracking;
+ /// Property: Steer printout at tracking action end
+ bool m_printEndTracking;
+ /// Property: Energy cut below which particles are not collected, but assigned to the parent
+ double m_kinEnergyCut;
+ /// Local buffer about the 'current' G4Track
+ Particle m_currTrack;
+ /// Map with stored MC Particles
+ ParticleMap m_particleMap;
+ /// Map associating the G4Track identifiers with identifiers of existing MCParticles
+ TrackEquivalents m_equivalentTracks;
+ /// All the processes of which the decay products will be explicitly stored
+ Processes m_processNames;
- int recombineParents();
- /// Print record of kept particles
- void printParticles();
- /// Clear particle maps
- void clear();
- /// Check the record consistency
- void checkConsistency() const;
+ /// Recombine particles and associate the to parents with cleanup
+ int recombineParents();
+ /// Clear particle maps
+ void clear();
+ /// Check the record consistency
+ void checkConsistency() const;
+ /// Get proper equivalent track from the particle map according to the given geant4 track ID
+ int equivalentTrack(int g4_id) const;
- public:
- /// Standard constructor
- Geant4ParticleHandler(Geant4Context* context, const std::string& nam);
- /// Default destructor
- virtual ~Geant4ParticleHandler();
- /// Access the Geant4 tracking manager. Only use between tracking pre- and post action
- G4TrackingManager* trackMgr() const { return m_context.trackMgr(); }
- /// Event generation action callback
- virtual void operator()(G4Event* event);
- /// User stepping callback
- virtual void step(const G4Step* step, G4SteppingManager* mgr);
- /// Pre-event action callback
- virtual void beginEvent(const G4Event* event);
- /// Post-event action callback
- virtual void endEvent(const G4Event* event);
- /// Pre-track action callback
- virtual void begin(const G4Track* track);
- /// Post-track action callback
- virtual void end(const G4Track* track);
+ public:
+ /// Standard constructor
+ Geant4ParticleHandler(Geant4Context* context, const std::string& nam);
+ /// Default destructor
+ virtual ~Geant4ParticleHandler();
+ /// Event generation action callback
+ virtual void operator()(G4Event* event);
+ /// User stepping callback
+ virtual void step(const G4Step* step, G4SteppingManager* mgr);
+ /// Pre-event action callback
+ virtual void beginEvent(const G4Event* event);
+ /// Post-event action callback
+ virtual void endEvent(const G4Event* event);
+ /// Pre-track action callback
+ virtual void begin(const G4Track* track);
+ /// Post-track action callback
+ virtual void end(const G4Track* track);
- /// Mark a Geant4 track to be kept for later MC truth analysis
- virtual void mark(const G4Track* track);
- /// Store a track
- virtual void mark(const G4Track* track, bool created_hit);
- /// Mark a Geant4 track of the step to be kept for later MC truth analysis
- virtual void mark(const G4Step* step);
- /// Store a track produced in a step to be kept for later MC truth analysis
- virtual void mark(const G4Step* step, bool created_hit);
+ /// Access the particle map
+ virtual const ParticleMap& particles() const { return m_particleMap; }
+ /// Access the map of track equivalents
+ virtual const TrackEquivalents& equivalents() const { return m_equivalentTracks; }
+ /// Mark a Geant4 track to be kept for later MC truth analysis. Default flag: CREATED_HIT
+ virtual void mark(const G4Track* track);
+ /// Store a track
+ virtual void mark(const G4Track* track, int reason);
+ /// Mark a Geant4 track of the step to be kept for later MC truth analysis. Default flag: CREATED_HIT
+ virtual void mark(const G4Step* step);
+ /// Store a track produced in a step to be kept for later MC truth analysis
+ virtual void mark(const G4Step* step, int reason);
- };
+ };
} // End namespace Simulation
} // End namespace DD4hep
diff --git a/DDG4/include/DDG4/Geant4ParticlePrint.h b/DDG4/include/DDG4/Geant4ParticlePrint.h
new file mode 100644
index 0000000000000000000000000000000000000000..1c5e2786e226a0bbe861e70ea2ee0cca68b587fd
--- /dev/null
+++ b/DDG4/include/DDG4/Geant4ParticlePrint.h
@@ -0,0 +1,62 @@
+// $Id: Geant4Hits.h 513 2013-04-05 14:31:53Z gaede $
+//====================================================================
+// AIDA Detector description implementation
+//--------------------------------------------------------------------
+//
+// Author : M.Frank
+//
+//====================================================================
+#ifndef DD4HEP_DDG4_GEANT4PARTICLEPRINT_H
+#define DD4HEP_DDG4_GEANT4PARTICLEPRINT_H
+
+// Framework include files
+#include "DDG4/Geant4EventAction.h"
+#include "DDG4/Geant4MonteCarloTruth.h"
+
+/*
+ * DD4hep namespace declaration
+ */
+namespace DD4hep {
+
+ /*
+ * Simulation namespace declaration
+ */
+ namespace Simulation {
+
+ /** Geant4Action to collect the MC particle information.
+ *
+ *
+ * @author M.Frank
+ * @version 1.0
+ */
+ class Geant4ParticlePrint : public Geant4EventAction {
+ public:
+ typedef Geant4MonteCarloTruth::ParticleMap ParticleMap;
+ typedef Geant4MonteCarloTruth::TrackEquivalents TrackEquivalents;
+ protected:
+ /// Property: Flag to indicate output type: 1: TABLE, 2:TREE, 3:BOTH (default)
+ int m_outputType;
+
+ void printParticle(const std::string& prefix, const Particle* p) const;
+ /// Print record of kept particles
+ void printParticles(const ParticleMap& particles) const;
+ /// Print tree of kept particles
+ void printParticleTree(const ParticleMap& particles, int level, const Particle* p) const;
+ /// Print tree of kept particles
+ void printParticleTree(const ParticleMap& particles) const;
+
+
+ public:
+ /// Standard constructor
+ Geant4ParticlePrint(Geant4Context* context, const std::string& nam);
+ /// Default destructor
+ virtual ~Geant4ParticlePrint();
+ /// Pre-event action callback
+ virtual void begin(const G4Event* event);
+ /// Post-event action callback
+ virtual void end(const G4Event* event);
+ };
+ } // End namespace Simulation
+} // End namespace DD4hep
+
+#endif // DD4HEP_DDG4_GEANT4PARTICLEPRINT_H
diff --git a/DDG4/include/DDG4/Geant4TrackPersistency.h b/DDG4/include/DDG4/Geant4TrackPersistency.h
index 2aed25b6e4f5d6d435c8d4c838d8903e21361b49..4d050b60c3aa1912811204b85a8284d4f9c6fbba 100644
--- a/DDG4/include/DDG4/Geant4TrackPersistency.h
+++ b/DDG4/include/DDG4/Geant4TrackPersistency.h
@@ -88,6 +88,10 @@ namespace DD4hep {
};
protected:
+ /// Map with stored MC Particles
+ ParticleMap m_particleMap;
+ /// Map associating the G4Track identifiers with identifiers of existing MCParticles
+ TrackEquivalents m_equivalentTracks;
/// Information block for current track
Geant4TrackPersistency::TrackInfo m_current;
/// Property: Steer printout at tracking action begin
@@ -105,7 +109,10 @@ namespace DD4hep {
virtual ~Geant4TrackPersistency();
/// Access the Geant4 tracking manager. Only use between tracking pre- and post action
G4TrackingManager* trackMgr() const { return m_context.trackMgr(); }
- /// Event generation action callback
+ /// Access the particle map
+ virtual const ParticleMap& particles() const { return m_particleMap; }
+ /// Access the map of track equivalents
+ virtual const TrackEquivalents& equivalents() const { return m_equivalentTracks; } /// Event generation action callback
virtual void operator()(G4Event* event);
/// User stepping callback
virtual void step(const G4Step* step, G4SteppingManager* mgr);
@@ -121,11 +128,11 @@ namespace DD4hep {
/// Mark a Geant4 track to be kept for later MC truth analysis
virtual void mark(const G4Track* track);
/// Store a track
- virtual void mark(const G4Track* track, bool created_hit);
+ virtual void mark(const G4Track* track, int reason);
/// Mark a Geant4 track of the step to be kept for later MC truth analysis
virtual void mark(const G4Step* step);
/// Store a track produced in a step to be kept for later MC truth analysis
- virtual void mark(const G4Step* step, bool created_hit);
+ virtual void mark(const G4Step* step, int reason);
};
} // End namespace Simulation
diff --git a/DDG4/plugins/Geant4Factories.cpp b/DDG4/plugins/Geant4Factories.cpp
index 70bfdccd17ae20b14f9468df5105d5bf4d1327df..169105fb561890a0ccc3dc72eba53fbb44cf4c87 100644
--- a/DDG4/plugins/Geant4Factories.cpp
+++ b/DDG4/plugins/Geant4Factories.cpp
@@ -49,6 +49,9 @@ DECLARE_GEANT4ACTION(Geant4MonteCarloRecordManager)
#include "DDG4/Geant4ParticleHandler.h"
DECLARE_GEANT4ACTION(Geant4ParticleHandler)
+#include "DDG4/Geant4ParticlePrint.h"
+DECLARE_GEANT4ACTION(Geant4ParticlePrint)
+
//=============================
#include "DDG4/Geant4TrackingPreAction.h"
DECLARE_GEANT4ACTION(Geant4TrackingPreAction)
diff --git a/DDG4/src/Geant4Data.cpp b/DDG4/src/Geant4Data.cpp
index 5f08fbe82188371dd18ed5b3082ce31db037943a..684d795679478304dffb586bf26dc5e32e7f6d39 100644
--- a/DDG4/src/Geant4Data.cpp
+++ b/DDG4/src/Geant4Data.cpp
@@ -41,6 +41,45 @@ SimpleEvent::~SimpleEvent() {
InstanceCount::decrement(this);
}
+/// Copy constructor
+Particle::Particle(const Particle& c)
+ : id(c.id), g4Parent(c.g4Parent), parent(c.parent), reason(c.reason), steps(c.steps),
+ vsx(c.vsx), vsy(c.vsy), vsz(c.vsz),
+ vex(c.vex), vey(c.vey), vez(c.vez),
+ psx(c.psx), psy(c.psy), psz(c.psz),
+ pex(c.pex), pey(c.pey), pez(c.pez),
+ energy(c.energy), time(c.time),
+ process(c.process), definition(c.definition),
+ daughters(c.daughters)
+{
+ InstanceCount::increment(this);
+}
+
+/// Default constructor
+Particle::Particle()
+ : id(0), g4Parent(0), parent(0), reason(0), steps(0),
+ vsx(0.0), vsy(0.0), vsz(0.0),
+ vex(0.0), vey(0.0), vez(0.0),
+ psx(0.0), psy(0.0), psz(0.0),
+ pex(0.0), pey(0.0), pez(0.0),
+ energy(0.0), time(0.0),
+ process(0), definition(0),
+ daughters()
+{
+ InstanceCount::increment(this);
+}
+
+/// Default destructor
+Particle::~Particle() {
+ InstanceCount::decrement(this);
+}
+
+/// Remove daughter from set
+void Particle::removeDaughter(int id_daughter) {
+ set<int>::iterator j = daughters.find(id_daughter);
+ if ( j != daughters.end() ) daughters.erase(j);
+}
+
/// Default constructor
SimpleHit::SimpleHit()
: cellID(0) {
diff --git a/DDG4/src/Geant4MonteCarloTruth.cpp b/DDG4/src/Geant4MonteCarloTruth.cpp
index 0f4d9d96167e3fae69a8a77ee2f8e4a86af8eeb0..bfa1a452007a28e3bbc43a020a8d8fd8881e79a3 100644
--- a/DDG4/src/Geant4MonteCarloTruth.cpp
+++ b/DDG4/src/Geant4MonteCarloTruth.cpp
@@ -38,7 +38,7 @@ void Geant4DummyTruthHandler::mark(const G4Track*) {
}
/// Store a track, with a flag
-void Geant4DummyTruthHandler::mark(const G4Track*, bool ) {
+void Geant4DummyTruthHandler::mark(const G4Track*, int ) {
}
/// Mark a Geant4 track of the step to be kept for later MC truth analysis
@@ -46,6 +46,6 @@ void Geant4DummyTruthHandler::mark(const G4Step*) {
}
/// Store a track produced in a step to be kept for later MC truth analysis
-void Geant4DummyTruthHandler::mark(const G4Step*, bool ) {
+void Geant4DummyTruthHandler::mark(const G4Step*, int ) {
}
diff --git a/DDG4/src/Geant4ParticleHandler.cpp b/DDG4/src/Geant4ParticleHandler.cpp
index e0bbe8c1b5a4a7d13001159b663f793e935f01d8..033be704e85d8aea32ec84f4bd74dbccdd76f02d 100644
--- a/DDG4/src/Geant4ParticleHandler.cpp
+++ b/DDG4/src/Geant4ParticleHandler.cpp
@@ -10,13 +10,13 @@
#include "DD4hep/Printout.h"
#include "DD4hep/Primitives.h"
#include "DD4hep/InstanceCount.h"
+#include "DDG4/Geant4StepHandler.h"
#include "DDG4/Geant4TrackHandler.h"
#include "DDG4/Geant4EventAction.h"
#include "DDG4/Geant4TrackingAction.h"
#include "DDG4/Geant4SteppingAction.h"
#include "DDG4/Geant4MonteCarloRecordManager.h"
#include "DDG4/Geant4ParticleHandler.h"
-#include "DDG4/Geant4StepHandler.h"
#include "G4TrackingManager.hh"
#include "G4ParticleDefinition.hh"
@@ -26,88 +26,30 @@
#include "G4PrimaryParticle.hh"
#include "G4PrimaryVertex.hh"
#include "G4Event.hh"
+#include "CLHEP/Units/SystemOfUnits.h"
#include <set>
#include <stdexcept>
using namespace std;
-using DD4hep::InstanceCount;
using namespace DD4hep;
using namespace DD4hep::Simulation;
+typedef ReferenceBitMask<int> PropertyMask;
-static double m_kinEnergyCut = 500e0;
-
-enum { CREATED_HIT = 1<<0,
- PRIMARY_PARTICLE = 1<<1,
- SECONDARIES = 1<<2,
- ENERGY = 1<<3,
- KEEP_PROCESS = 1<<4,
- KEEP_PARENT = 1<<5,
- CREATED_CALORIMETER_HIT = 1<<6,
- CREATED_TRACKER_HIT = 1<<7,
-
- LAST_NOTHING = 1<<31
-};
-
-
-template <typename T> class BitMask {
-public:
- T& refMask;
- BitMask(T& m) : refMask(m) {}
- T value() const {
- return refMask;
- }
- void set(const T& m) {
- refMask |= m;
- }
- bool isSet(const T& m) {
- return (refMask&m) == m;
- }
- bool testBit(int bit) const {
- T m = T(1)<<bit;
- return isSet(m);
- }
- bool isNull() const {
- return refMask == 0;
+namespace {
+ G4PrimaryParticle* primary(int id, const G4Event& evt) {
+ for(int i=0, ni=evt.GetNumberOfPrimaryVertex(); i<ni; ++i) {
+ G4PrimaryVertex* v = evt.GetPrimaryVertex(i);
+ for(int j=0, nj=v->GetNumberOfParticle(); j<nj; ++j) {
+ G4PrimaryParticle* p = v->GetPrimary(j);
+ if ( id == p->GetTrackID() ) {
+ return p;
+ }
+ }
+ }
+ return 0;
}
-};
-template <typename T> inline BitMask<T> bitMask(T& m) { return BitMask<T>(m); }
-template <typename T> inline BitMask<const T> bitMask(const T& m) { return BitMask<const T>(m); }
-
-/// Copy constructor
-Geant4ParticleHandler::Particle::Particle(const Particle& c)
- : id(c.id), parent(c.parent), theParent(c.theParent), reason(c.reason),
- vx(c.vx), vy(c.vy), vz(c.vz),
- px(c.px), py(c.py), pz(c.pz),
- energy(c.energy), time(c.time),
- process(c.process), definition(c.definition),
- daughters(c.daughters)
-{
- InstanceCount::increment(this);
-}
-
-/// Default constructor
-Geant4ParticleHandler::Particle::Particle()
- : id(0), parent(0), theParent(0), reason(0),
- vx(0.0), vy(0.0), vz(0.0),
- px(0.0), py(0.0), pz(0.0),
- energy(0.0), time(0.0),
- process(0), definition(0),
- daughters()
-{
- InstanceCount::increment(this);
-}
-
-/// Default destructor
-Geant4ParticleHandler::Particle::~Particle() {
- InstanceCount::decrement(this);
-}
-
-/// Remove daughter from set
-void Geant4ParticleHandler::Particle::removeDaughter(int id_daughter) {
- set<int>::iterator j = daughters.find(id_daughter);
- if ( j != daughters.end() ) daughters.erase(j);
}
/// Standard constructor
@@ -116,14 +58,15 @@ Geant4ParticleHandler::Geant4ParticleHandler(Geant4Context* context, const std::
{
//generatorAction().adopt(this);
eventAction().callAtBegin(this,&Geant4ParticleHandler::beginEvent);
- eventAction().callAtFinal(this,&Geant4ParticleHandler::endEvent);
+ eventAction().callAtEnd(this,&Geant4ParticleHandler::endEvent);
trackingAction().callAtFinal(this,&Geant4ParticleHandler::end,CallbackSequence::FRONT);
trackingAction().callUpFront(this,&Geant4ParticleHandler::begin,CallbackSequence::FRONT);
steppingAction().call(this,&Geant4ParticleHandler::step);
- declareProperty("PrintEndTracking",m_printEndTracking = false);
- declareProperty("PrintStartTracking",m_printStartTracking = false);
- declareProperty("SaveProcesses",m_processNames);
+ declareProperty("printEndTracking",m_printEndTracking = false);
+ declareProperty("printStartTracking",m_printStartTracking = false);
+ declareProperty("saveProcesses",m_processNames);
+ declareProperty("minimalKineticEnergy",m_kinEnergyCut = 100e0*MeV);
InstanceCount::increment(this);
}
@@ -142,16 +85,20 @@ void Geant4ParticleHandler::clear() {
}
/// Mark a Geant4 track to be kept for later MC truth analysis
-void Geant4ParticleHandler::mark(const G4Track* track, bool created_hit) {
- mark(track); // Does all the checking...
- if ( created_hit ) {
+void Geant4ParticleHandler::mark(const G4Track* track, int reason) {
+ if ( track ) {
+ if ( reason != 0 ) {
+ PropertyMask(m_currTrack.reason).set(reason);
+ return;
+ }
}
+ except("Cannot mark the G4Track if the pointer is invalid!", c_name());
}
/// Store a track produced in a step to be kept for later MC truth analysis
-void Geant4ParticleHandler::mark(const G4Step* step, bool created_hit) {
+void Geant4ParticleHandler::mark(const G4Step* step, int reason) {
if ( step ) {
- mark(step->GetTrack(),created_hit);
+ mark(step->GetTrack(),reason);
return;
}
except("Cannot mark the G4Track if the step-pointer is invalid!", c_name());
@@ -168,24 +115,20 @@ void Geant4ParticleHandler::mark(const G4Step* step) {
/// Mark a Geant4 track of the step to be kept for later MC truth analysis
void Geant4ParticleHandler::mark(const G4Track* track) {
- BitMask<int> mask(m_currTrack.reason);
- mask.set(CREATED_HIT);
+ PropertyMask mask(m_currTrack.reason);
+ mask.set(G4PARTICLE_CREATED_HIT);
/// Check if the track origines from the calorimeter.
// If yes, flag it, because it is a candidate fro removal.
G4VPhysicalVolume* vol = track->GetVolume();
if ( strstr(vol->GetName().c_str(),"cal") ) { // just for test!
- mask.set(CREATED_CALORIMETER_HIT);
+ mask.set(G4PARTICLE_CREATED_CALORIMETER_HIT);
}
- else if ( !mask.isSet(CREATED_TRACKER_HIT) ) {
+ else if ( !mask.isSet(G4PARTICLE_CREATED_TRACKER_HIT) ) {
//printout(INFO,name(),"+++ Create TRACKER HIT: id=%d",m_currTrack.id);
- mask.set(CREATED_TRACKER_HIT);
+ mask.set(G4PARTICLE_CREATED_TRACKER_HIT);
}
}
-/// Check the record consistency
-void Geant4ParticleHandler::checkConsistency() const {
-}
-
/// Event generation action callback
void Geant4ParticleHandler::operator()(G4Event*) {
typedef Geant4MonteCarloTruth _MC;
@@ -194,12 +137,153 @@ void Geant4ParticleHandler::operator()(G4Event*) {
clear();
}
+/// User stepping callback
+void Geant4ParticleHandler::step(const G4Step* step, G4SteppingManager* /* mgr */) {
+ typedef std::vector<const G4Track*> _Sec;
+ ++m_currTrack.steps;
+ if ( m_currTrack.energy > m_kinEnergyCut ) {
+ //
+ // Tracks below the energy threshold are NOT stored.
+ // If these tracks produce hits or are selected due to another signature,
+ // this criterium will anyhow take precedence.
+ //
+ const _Sec* sec=step->GetSecondaryInCurrentStep();
+ if ( sec->size() > 0 ) {
+ PropertyMask(m_currTrack.reason).set(G4PARTICLE_HAS_SECONDARIES);
+ }
+ }
+}
+
+/// Pre-track action callback
+void Geant4ParticleHandler::begin(const G4Track* track) {
+ Geant4TrackHandler h(track);
+ double kine = h.kineticEnergy();
+ G4ThreeVector m = track->GetMomentum();
+
+ // Extract here all the information from the start of the tracking action
+ // which we will need later to create a proper MC particle.
+ m_currTrack.id = h.id();
+ m_currTrack.reason = (kine > m_kinEnergyCut) ? G4PARTICLE_ABOVE_ENERGY_THRESHOLD : 0;
+ m_currTrack.energy = kine;
+ m_currTrack.g4Parent = h.parent();
+ m_currTrack.parent = h.parent();
+ m_currTrack.definition = h.trackDef();
+ m_currTrack.process = h.creatorProcess();
+ m_currTrack.time = h.globalTime();
+ m_currTrack.vsx = h.vertex().x();
+ m_currTrack.vsy = h.vertex().y();
+ m_currTrack.vsz = h.vertex().z();
+ m_currTrack.vex = 0.0;
+ m_currTrack.vey = 0.0;
+ m_currTrack.vez = 0.0;
+ m_currTrack.psx = m.x();
+ m_currTrack.psy = m.y();
+ m_currTrack.psz = m.z();
+ m_currTrack.pex = 0.0;
+ m_currTrack.pey = 0.0;
+ m_currTrack.pez = 0.0;
+ // If the creator process of the track is in the list of process products to be kept, set the proper flag
+ if ( m_currTrack.process ) {
+ Processes::iterator i=find(m_processNames.begin(),m_processNames.end(),m_currTrack.process->GetProcessName());
+ if ( i != m_processNames.end() ) {
+ PropertyMask(m_currTrack.reason).set(G4PARTICLE_KEEP_PROCESS);
+ }
+ }
+}
+
+/// Post-track action callback
+void Geant4ParticleHandler::end(const G4Track* track) {
+ Geant4TrackHandler h(track);
+ int id = h.id();
+ PropertyMask mask(m_currTrack.reason);
+ G4PrimaryParticle* prim = primary(id,context()->event().event());
+
+ if ( prim ) {
+ m_currTrack.reason |= (G4PARTICLE_PRIMARY|G4PARTICLE_ABOVE_ENERGY_THRESHOLD);
+ const G4ParticleDefinition* def = m_currTrack.definition;
+ printout(INFO,name(),"+++ Add Primary particle %d def:%p [%s, %s] reason:%d",
+ id, def,
+ def ? def->GetParticleName().c_str() : "UNKNOWN",
+ def ? def->GetParticleType().c_str() : "UNKNOWN", m_currTrack.reason);
+ }
+
+ if ( !mask.isNull() ) {
+ // These are candate tracks with a probability to be stored due to their properties:
+ // - primary particle
+ // - hits created
+ // - secondaries
+ // - above energy threshold
+ // - to be kept due to creator process
+ //
+ // Update vertex end point and final momentum
+ G4ThreeVector m = track->GetMomentum();
+ m_currTrack.vex = h.vertex().x();
+ m_currTrack.vey = h.vertex().y();
+ m_currTrack.vez = h.vertex().z();
+ m_currTrack.pex = m.x();
+ m_currTrack.pey = m.y();
+ m_currTrack.pez = m.z();
+
+ // Create a new MC particle from the current track information saved in the pre-tracking action
+ Particle* p = m_particleMap[id] = new Particle(m_currTrack);
+ // Add this track to it's parents list of daughters
+ ParticleMap::iterator ipar = m_particleMap.find(p->g4Parent);
+ if ( ipar != m_particleMap.end() ) {
+ Particle* p_par = (*ipar).second;
+ p_par->daughters.insert(id);
+ }
+ m_equivalentTracks[id] = id;
+#if 0
+ /// Some printout
+ const char* hit = mask.isSet(G4PARTICLE_CREATED_HIT) ? "CREATED_HIT" : "";
+ const char* sec = mask.isSet(G4PARTICLE_HAS_SECONDARIES) ? "SECONDARIES" : "";
+ const char* prim = mask.isSet(G4PARTICLE_PRIMARY) ? "PRIMARY" : "";
+ printout(INFO,name(),"+++ Tracking postaction: %d/%d par:%d [%s, %s] created by:%s E:%e state:%s %s %s %s",
+ id,int(m_particleMap.size()),
+ h.parent(),h.name().c_str(),h.type().c_str(),h.creatorName().c_str(),
+ p->energy, h.statusName(), prim, hit, sec);
+#endif
+ }
+ else {
+ // These are tracks without any special properties.
+ //
+ // We will not store them on the record, but have to memorise the
+ // track identifier in order to restore the history for the created hits.
+ m_equivalentTracks[id] = h.parent();
+ }
+}
+
/// Pre-event action callback
void Geant4ParticleHandler::beginEvent(const G4Event* ) {
- for(ParticleMap::const_iterator i=m_particleMap.begin(); i!=m_particleMap.end(); ++i) {
- //Particle* part = (*i).second;
- //int equiv_id = part->parent;
+}
+
+/// Post-event action callback
+void Geant4ParticleHandler::endEvent(const G4Event* ) {
+ int count = 0;
+ do {
+ printout(INFO,name(),"+++ Iteration:%d Tracks:%d Equivalents:%d",++count,m_particleMap.size(),m_equivalentTracks.size());
+ } while( recombineParents() > 0 );
+ // Update the particle's parents and replace the original Geant4 track with the
+ // equivalent particle still present in the record.
+ for(ParticleMap::const_iterator ipar, iend=m_particleMap.end(), i=m_particleMap.begin(); i!=iend; ++i) {
+ Particle* p = (*i).second;
+ p->parent = equivalentTrack(p->g4Parent);
+ if ( p->parent > 0 ) {
+ ipar = m_particleMap.find(p->parent);
+ if ( ipar != iend ) {
+ set<int>& daughters = (*ipar).second->daughters;
+ if ( daughters.find(p->id) == daughters.end() ) daughters.insert(p->id);
+ }
+ }
+ }
+ /// No we have to updte the map of equivalent tracks and assign the 'equivalentTrack' entry
+ for(TrackEquivalents::iterator i=m_equivalentTracks.begin(); i!=m_equivalentTracks.end(); ++i) {
+ int& track_id = (*i).second;
+ int equiv_id = equivalentTrack(track_id);
+ track_id = equiv_id;
}
+ // Consistency check....
+ checkConsistency();
}
int Geant4ParticleHandler::recombineParents() {
@@ -210,20 +294,20 @@ int Geant4ParticleHandler::recombineParents() {
Particle* par = (*i).second;
set<int>& daughters = par->daughters;
if ( daughters.size() == 0 ) {
- BitMask<const int> mask(par->reason);
+ PropertyMask mask(par->reason);
int id = par->id;
- bool secondaries = mask.isSet(SECONDARIES);
- bool tracker_track = mask.isSet(CREATED_TRACKER_HIT);
- bool calo_track = mask.isSet(CREATED_CALORIMETER_HIT);
- bool hits_produced = mask.isSet(CREATED_HIT);
- bool low_energy = !mask.isSet(ENERGY);
- bool keep_process = mask.isSet(KEEP_PROCESS);
- bool keep_parent = mask.isSet(KEEP_PARENT);
- int parent_id = par->parent;
+ bool secondaries = mask.isSet(G4PARTICLE_HAS_SECONDARIES);
+ bool tracker_track = mask.isSet(G4PARTICLE_CREATED_TRACKER_HIT);
+ bool calo_track = mask.isSet(G4PARTICLE_CREATED_CALORIMETER_HIT);
+ bool hits_produced = mask.isSet(G4PARTICLE_CREATED_HIT);
+ bool low_energy = !mask.isSet(G4PARTICLE_ABOVE_ENERGY_THRESHOLD);
+ bool keep_process = mask.isSet(G4PARTICLE_KEEP_PROCESS);
+ bool keep_parent = mask.isSet(G4PARTICLE_KEEP_PARENT);
+ int parent_id = par->g4Parent;
bool remove_me = false;
/// Primary particles MUST be kept!
- if ( mask.isSet(PRIMARY_PARTICLE) ) {
+ if ( mask.isSet(G4PARTICLE_PRIMARY) ) {
continue;
}
else if ( keep_parent ) {
@@ -233,9 +317,9 @@ int Geant4ParticleHandler::recombineParents() {
ParticleMap::iterator ip = m_particleMap.find(parent_id);
if ( ip != m_particleMap.end() ) {
Particle* parent_part = (*ip).second;
- BitMask<const int> parent_mask(parent_part->reason);
- if ( parent_mask.isSet(ENERGY) ) {
- parent_part->reason |= KEEP_PARENT;
+ PropertyMask parent_mask(parent_part->reason);
+ if ( parent_mask.isSet(G4PARTICLE_ABOVE_ENERGY_THRESHOLD) ) {
+ parent_mask.set(G4PARTICLE_KEEP_PARENT);
continue;
}
}
@@ -266,7 +350,7 @@ int Geant4ParticleHandler::recombineParents() {
m_equivalentTracks[id] = parent_id;
if ( ip != m_particleMap.end() ) {
Particle* parent_part = (*ip).second;
- bitMask(parent_part->reason).set(mask.value());
+ PropertyMask(parent_part->reason).set(mask.value());
// Remove track from parent's list of daughters
parent_part->removeDaughter(id);
}
@@ -284,173 +368,66 @@ int Geant4ParticleHandler::recombineParents() {
return int(remove.size());
}
-/// Pre-event action callback
-void Geant4ParticleHandler::endEvent(const G4Event* ) {
- int count = 0;
- do {
- printout(INFO,name(),"+++ Iteration:%d Tracks:%d Equivalents:%d",++count,m_particleMap.size(),m_equivalentTracks.size());
- } while( recombineParents() > 0 );
-
- for(ParticleMap::const_iterator ipar, i=m_particleMap.begin(); i!=m_particleMap.end(); ++i) {
- Particle* part = (*i).second;
- int equiv_id = part->parent;
- if ( equiv_id != 0 ) {
- while( (ipar=m_particleMap.find(equiv_id)) == m_particleMap.end() ) {
- TrackEquivalents::const_iterator iequiv = m_equivalentTracks.find(equiv_id);
- equiv_id = (iequiv == m_equivalentTracks.end()) ? -1 : (*iequiv).second;
- if ( equiv_id < 0 ) {
- printout(INFO,name(),"+++ No Equivalent particle for track:%d last known is:%d",part->id,equiv_id);
- break;
- }
+/// Check the record consistency
+void Geant4ParticleHandler::checkConsistency() const {
+ int num_errors = 0;
+
+ /// First check the consistency of the particle map itself
+ for(ParticleMap::const_iterator j, i=m_particleMap.begin(); i!=m_particleMap.end(); ++i) {
+ Particle* p = (*i).second;
+ PropertyMask mask(p->reason);
+ set<int>& daughters = p->daughters;
+ // For all particles, the set of daughters must be contained in the record.
+ for(set<int>::const_iterator id=daughters.begin(); id!=daughters.end(); ++id) {
+ int id_dau = *id;
+ j = m_particleMap.find(id_dau);
+ if ( j == m_particleMap.end() ) {
+ ++num_errors;
+ printout(INFO,name(),"+++ Particle:%d Daughter %d is not in particle map!",p->id,id_dau);
}
- if ( equiv_id>0 && equiv_id != part->parent ) part->theParent = equiv_id;
}
- }
- checkConsistency();
- printParticles();
-}
-
-/// Print record of kept particles
-void Geant4ParticleHandler::printParticles() {
- char equiv[32];
- ParticleMap::const_iterator ipar;
- printout(INFO,name(),"+++ MC Particles #Tracks:%6d %-12s Parent%-7s "
- "Primary Secondaries Energy %9s Calo Tracker Process",
- int(m_particleMap.size()),"ParticleType","","in [MeV]");
- for(ParticleMap::const_iterator i=m_particleMap.begin(); i!=m_particleMap.end(); ++i) {
- Particle* part = (*i).second;
- BitMask<const int> mask(part->reason);
- equiv[0] = 0;
- if ( part->parent != part->theParent ) {
- ::snprintf(equiv,sizeof(equiv),"/%d",part->theParent);
+ // For all particles except the primaries, the parent must be contained in the record.
+ if ( !mask.isSet(G4PARTICLE_PRIMARY) ) {
+ int id_par = p->parent;
+ j = m_particleMap.find(id_par);
+ if ( j == m_particleMap.end() ) {
+ ++num_errors;
+ printout(INFO,name(),"+++ Particle:%d Parent %d is not in particle map!",p->id,id_par);
+ }
}
- printout(INFO,name(),"+++ MC Particle TrackID: %6d %-12s %6d%-7s "
- "%-7s %-11s %-7s %8.3g %-4s %-7s %-7s [%s/%s]",
- part->id,
- part->definition->GetParticleName().c_str(),
- part->parent,equiv,
- yes_no(mask.isSet(PRIMARY_PARTICLE)),
- yes_no(mask.isSet(SECONDARIES)),
- yes_no(mask.isSet(ENERGY)),
- part->energy,
- yes_no(mask.isSet(CREATED_CALORIMETER_HIT)),
- yes_no(mask.isSet(CREATED_TRACKER_HIT)),
- yes_no(mask.isSet(KEEP_PROCESS)),
- part->process ? part->process->GetProcessName().c_str() : "???",
- part->process ? G4VProcess::GetProcessTypeName(part->process->GetProcessType()).c_str() : ""
- );
}
- printout(INFO,"Summary","+++ MC Particles #Tracks:%6d %-12s Parent%-7s "
- "Primary Secondaries Energy %9s Calo Tracker Process",
- int(m_particleMap.size()),"ParticleType","","in [MeV]");
-}
-/// User stepping callback
-void Geant4ParticleHandler::step(const G4Step* step, G4SteppingManager* /* mgr */) {
- typedef std::vector<const G4Track*> _Sec;
- if ( m_currTrack.energy > m_kinEnergyCut ) {
- //
- // Tracks below the energy threshold are NOT stored.
- // If these tracks produce hits or are selected due to another signature,
- // this criterium will anyhow take precedence.
- //
- const _Sec* sec=step->GetSecondaryInCurrentStep();
- if ( sec->size() > 0 ) {
- bitMask(m_currTrack.reason).set(SECONDARIES);
+ /// No we have to check the consistency of the map of equivalent tracks used to assign the
+ /// proper MC particle to the created hits
+ for(TrackEquivalents::const_iterator i=m_equivalentTracks.begin(); i!=m_equivalentTracks.end(); ++i) {
+ int track_id = (*i).second;
+ int equiv_id = equivalentTrack(track_id);
+ ParticleMap::const_iterator j = m_particleMap.find(equiv_id);
+ if ( j == m_particleMap.end() ) {
+ int g4_id = (*i).first;
+ ++num_errors;
+ printout(INFO,name(),"+++ Geant4 Track:%d Track:%d Equivalent track %d is not in particle map!",g4_id,track_id,equiv_id);
}
}
-}
-
-/// Pre-track action callback
-void Geant4ParticleHandler::begin(const G4Track* track) {
- Geant4TrackHandler h(track);
- double kine = h.kineticEnergy();
- G4ThreeVector m = track->GetMomentum();
- // Extract here all the information from the start of the tracking action
- // which we will need later to create a proper MC particle.
- m_currTrack.id = h.id();
- m_currTrack.reason = (kine > m_kinEnergyCut) ? ENERGY : 0;
- m_currTrack.energy = kine;
- m_currTrack.parent = h.parent();
- m_currTrack.theParent = h.parent();
- m_currTrack.definition = h.trackDef();
- m_currTrack.process = h.creatorProcess();
- m_currTrack.time = h.globalTime();
- m_currTrack.vx = h.vertex().x();
- m_currTrack.vy = h.vertex().y();
- m_currTrack.vz = h.vertex().z();
- m_currTrack.px = m.x();
- m_currTrack.py = m.y();
- m_currTrack.pz = m.z();
- // If the creator process of the track is in the list of process products to be kept, set the proper flag
- if ( m_currTrack.process ) {
- Processes::iterator i=find(m_processNames.begin(),m_processNames.end(),m_currTrack.process->GetProcessName());
- if ( i != m_processNames.end() ) {
- bitMask(m_currTrack.reason).set(KEEP_PROCESS);
- }
+ if ( num_errors > 0 ) {
+ printout(ERROR,name(),"+++ Consistency check failed. Found %d problems.",num_errors);
}
}
-static G4PrimaryParticle* primary(int id, const G4Event& evt) {
- for(int i=0, ni=evt.GetNumberOfPrimaryVertex(); i<ni; ++i) {
- G4PrimaryVertex* v = evt.GetPrimaryVertex(i);
- for(int j=0, nj=v->GetNumberOfParticle(); j<nj; ++j) {
- G4PrimaryParticle* p = v->GetPrimary(i);
- if ( id == p->GetTrackID() ) {
- return p;
+/// Get proper equivalent track from the particle map according to the given geant4 track ID
+int Geant4ParticleHandler::equivalentTrack(int g4_id) const {
+ int equiv_id = g4_id;
+ if ( g4_id != 0 ) {
+ ParticleMap::const_iterator ipar;
+ while( (ipar=m_particleMap.find(equiv_id)) == m_particleMap.end() ) {
+ TrackEquivalents::const_iterator iequiv = m_equivalentTracks.find(equiv_id);
+ equiv_id = (iequiv == m_equivalentTracks.end()) ? -1 : (*iequiv).second;
+ if ( equiv_id < 0 ) {
+ printout(INFO,name(),"+++ No Equivalent particle for track:%d last known is:%d",g4_id,equiv_id);
+ break;
}
}
}
- return 0;
-}
-
-/// Post-track action callback
-void Geant4ParticleHandler::end(const G4Track* track) {
- Geant4TrackHandler h(track);
- int id = h.id();
- BitMask<const int> mask(m_currTrack.reason);
- G4PrimaryParticle* prim = primary(id,context()->event().event());
-
- if ( prim ) {
- m_currTrack.reason |= (PRIMARY_PARTICLE|ENERGY);
- const G4ParticleDefinition* def = m_currTrack.definition;
- printout(INFO,name(),"+++ Add Primary particle %d def:%p [%s, %s] reason:%d",
- id, def,
- def ? def->GetParticleName().c_str() : "UNKNOWN",
- def ? def->GetParticleType().c_str() : "UNKNOWN", m_currTrack.reason);
- }
-
- if ( !mask.isNull() ) {
- // These are candate tracks with a probability to be stored due to their properties:
- // - no primary particle
- // - no hits created
- // - no secondaries
- // - not above energy threshold
- //
- // Create a new MC particle from the current track information saved in the pre-tracking action
- Particle* p = m_particleMap[id] = new Particle(m_currTrack);
- // Add this track to it's parents list of daughters
- ParticleMap::iterator ipar = m_particleMap.find(p->parent);
- if ( ipar != m_particleMap.end() ) {
- (*ipar).second->daughters.insert(id);
- }
-#if 0
- /// Some printout
- const char* hit = mask.isSet(CREATED_HIT) ? "CREATED_HIT" : "";
- const char* sec = mask.isSet(SECONDARIES) ? "SECONDARIES" : "";
- const char* prim = mask.isSet(PRIMARY_PARTICLE) ? "PRIMARY" : "";
- printout(INFO,name(),"+++ Tracking postaction: %d/%d par:%d [%s, %s] created by:%s E:%e state:%s %s %s %s",
- id,int(m_particleMap.size()),
- h.parent(),h.name().c_str(),h.type().c_str(),h.creatorName().c_str(),
- p->energy, h.statusName(), prim, hit, sec);
-#endif
- }
- else {
- // These are tracks without any special properties.
- //
- // We will not store them on the record, but have to memorise the track identifier
- // in order to restore the history for the created hits.
- m_equivalentTracks[id] = h.parent();
- }
+ return equiv_id;
}
diff --git a/DDG4/src/Geant4ParticlePrint.cpp b/DDG4/src/Geant4ParticlePrint.cpp
new file mode 100644
index 0000000000000000000000000000000000000000..152e08bec614977214c15b939da0f1a406c58cae
--- /dev/null
+++ b/DDG4/src/Geant4ParticlePrint.cpp
@@ -0,0 +1,152 @@
+// $Id: Geant4Field.cpp 888 2013-11-14 15:54:56Z markus.frank@cern.ch $
+//====================================================================
+// AIDA Detector description implementation for LCD
+//--------------------------------------------------------------------
+//
+// Author : M.Frank
+//
+//====================================================================
+// Framework include files
+#include "DD4hep/Printout.h"
+#include "DD4hep/Primitives.h"
+#include "DD4hep/InstanceCount.h"
+#include "DDG4/Geant4ParticlePrint.h"
+#include "DDG4/Geant4Data.h"
+
+#include "G4ParticleDefinition.hh"
+#include "G4VProcess.hh"
+#include "G4Event.hh"
+
+using namespace std;
+using namespace DD4hep;
+using namespace DD4hep::Simulation;
+
+typedef ReferenceBitMask<const int> PropertyMask;
+
+/// Standard constructor
+Geant4ParticlePrint::Geant4ParticlePrint(Geant4Context* context, const std::string& nam)
+ : Geant4EventAction(context,nam)
+{
+ declareProperty("OutputType",m_outputType=3);
+ InstanceCount::increment(this);
+}
+
+/// Default destructor
+Geant4ParticlePrint::~Geant4ParticlePrint() {
+ InstanceCount::decrement(this);
+}
+
+/// Pre-event action callback
+void Geant4ParticlePrint::begin(const G4Event* ) {
+}
+
+/// Post-event action callback
+void Geant4ParticlePrint::end(const G4Event* ) {
+ Geant4MonteCarloTruth* truth = context()->event().extension<Geant4MonteCarloTruth>();
+ const ParticleMap& particles = truth->particles();
+ // Table printout....
+ if ( (m_outputType&1) != 0 ) printParticles(particles);
+ // Tree printout....
+ if ( (m_outputType&2) != 0 ) printParticleTree(particles);
+}
+
+void Geant4ParticlePrint::printParticle(const std::string& prefix, const Particle* p) const {
+ char equiv[32];
+ PropertyMask mask(p->reason);
+ const char* proc_name = "???";
+ const char* proc_type = p->process ? G4VProcess::GetProcessTypeName(p->process->GetProcessType()).c_str() : "";
+
+ if ( p->process ) proc_name = p->process->GetProcessName().c_str();
+ else if ( mask.isSet(G4PARTICLE_PRIMARY) ) proc_name = "Primary";
+
+ equiv[0] = 0;
+ if ( p->g4Parent != p->parent ) {
+ ::snprintf(equiv,sizeof(equiv),"/%d",p->g4Parent);
+ }
+ printout(INFO,name(),"+++ %sID: %6d %12s %6d%-7s %7s %3s %5d %6s %8.3g %4s %7s %7s %3s [%s%s%s]",
+ prefix.c_str(),
+ p->id,
+ p->definition->GetParticleName().c_str(),
+ p->parent,equiv,
+ yes_no(mask.isSet(G4PARTICLE_PRIMARY)),
+ yes_no(mask.isSet(G4PARTICLE_HAS_SECONDARIES)),
+ int(p->daughters.size()),
+ yes_no(mask.isSet(G4PARTICLE_ABOVE_ENERGY_THRESHOLD)),
+ p->energy,
+ yes_no(mask.isSet(G4PARTICLE_CREATED_CALORIMETER_HIT)),
+ yes_no(mask.isSet(G4PARTICLE_CREATED_TRACKER_HIT)),
+ yes_no(mask.isSet(G4PARTICLE_KEEP_PROCESS)),
+ mask.isSet(G4PARTICLE_KEEP_PARENT) ? "YES" : "",
+ proc_name,
+ p->process ? "/" : "",
+ proc_type
+ );
+}
+
+/// Print record of kept particles
+void Geant4ParticlePrint::printParticles(const ParticleMap& particles) const {
+ int num_energy = 0;
+ int num_parent = 0;
+ int num_process = 0;
+ int num_primary = 0;
+ int num_secondaries = 0;
+ int num_calo_hits = 0;
+ int num_tracker_hits = 0;
+
+ printout(INFO,name(),"+++ MC Particles #Tracks:%6d %-12s Parent%-7s "
+ "Primary Secondary Energy %-9s Calo Tracker Process/Par Details",
+ int(particles.size()),"ParticleType","","in [MeV]");
+ for(ParticleMap::const_iterator i=particles.begin(); i!=particles.end(); ++i) {
+ const Particle* p = (*i).second;
+ PropertyMask mask(p->reason);
+ printParticle("MC Particle Track",p);
+ num_secondaries += int(p->daughters.size());
+ if ( mask.isSet(G4PARTICLE_ABOVE_ENERGY_THRESHOLD) ) ++num_energy;
+ if ( mask.isSet(G4PARTICLE_PRIMARY) ) ++num_primary;
+ if ( mask.isSet(G4PARTICLE_CREATED_TRACKER_HIT) ) ++num_tracker_hits;
+ if ( mask.isSet(G4PARTICLE_CREATED_CALORIMETER_HIT) ) ++num_calo_hits;
+ if ( mask.isSet(G4PARTICLE_KEEP_PARENT) ) ++num_parent;
+ else if ( mask.isSet(G4PARTICLE_KEEP_PROCESS) ) ++num_process;
+ }
+ printout(INFO,name(),"+++ MC Particles #Tracks:%6d %-12s Parent%-7s "
+ "Primary Secondary Energy %-9s Calo Tracker Process/Parent",
+ int(particles.size()),"ParticleType","","");
+ printout(INFO,name(),"+++ MC Particles #Tracks:%6d %-12s Parent%-7s "
+ "%7d %9d %6d %-9s %4d %7d %7d %6d",
+ int(particles.size()),"ParticleType","",
+ num_primary, num_secondaries, num_energy,"",
+ num_calo_hits,num_tracker_hits,num_process,num_parent);
+}
+
+void Geant4ParticlePrint::printParticleTree(const ParticleMap& particles, int level, const Particle* p) const {
+ char txt[32];
+ size_t len = sizeof(txt)-1;
+ // Ensure we do not overwrite the array
+ if ( level>int(len)-3 ) level=len-3;
+
+ ::snprintf(txt,sizeof(txt),"%5d ",level);
+ ::memset(txt+6,' ',len-6);
+ txt[len] = 0;
+ txt[len-2] = '>';
+ txt[level+6]='+';
+ ::memset(txt+level+6+1,'-',len-level-3-6);
+
+ printParticle(txt, p);
+ const set<int>& daughters = p->daughters;
+ // For all particles, the set of daughters must be contained in the record.
+ for(set<int>::const_iterator id=daughters.begin(); id!=daughters.end(); ++id) {
+ int id_dau = *id;
+ const Particle* dau = (*particles.find(id_dau)).second;
+ printParticleTree(particles,level+1,dau);
+ }
+}
+
+/// Print tree of kept particles
+void Geant4ParticlePrint::printParticleTree(const ParticleMap& particles) const {
+ printout(INFO,name(),"+++ MC Particle Parent daughter relationships. [%d particles]",int(particles.size()));
+ for(ParticleMap::const_iterator i=particles.begin(); i!=particles.end(); ++i) {
+ const Particle* p = (*i).second;
+ PropertyMask mask(p->reason);
+ if ( mask.isSet(G4PARTICLE_PRIMARY) ) printParticleTree(particles,0,p);
+ }
+}
diff --git a/DDG4/src/Geant4TrackPersistency.cpp b/DDG4/src/Geant4TrackPersistency.cpp
index 01692ab8be6bdad2677f40f3a60a30a4bb96df7a..272542523814d5608c4d62f2a652520d77b63ae8 100644
--- a/DDG4/src/Geant4TrackPersistency.cpp
+++ b/DDG4/src/Geant4TrackPersistency.cpp
@@ -351,19 +351,19 @@ Geant4TrackPersistency::~Geant4TrackPersistency() {
}
/// Mark a Geant4 track to be kept for later MC truth analysis
-void Geant4TrackPersistency::mark(const G4Track* track, bool created_hit) {
+void Geant4TrackPersistency::mark(const G4Track* track, int reason) {
mark(track); // Does all the checking...
- if ( created_hit ) {
+ if ( reason ) {
Track* trk = (Track*)m_current.info;
trk->flag |= Track::CREATED_HIT;
- m_current.createdHit = created_hit;
+ m_current.createdHit = reason != 0;
}
}
/// Store a track produced in a step to be kept for later MC truth analysis
-void Geant4TrackPersistency::mark(const G4Step* step, bool created_hit) {
+void Geant4TrackPersistency::mark(const G4Step* step, int reason) {
if ( step ) {
- mark(step->GetTrack(),created_hit);
+ mark(step->GetTrack(),reason);
return;
}
except("Cannot mark the G4Track if the step-pointer is invalid!", c_name());