diff --git a/DDG4/examples/CLICSidSimu.py b/DDG4/examples/CLICSidSimu.py
index 562deed19431a6c76ac88f7c813c04128209c87d..84370551108eb4573207681008b44298368ec14a 100644
--- a/DDG4/examples/CLICSidSimu.py
+++ b/DDG4/examples/CLICSidSimu.py
@@ -45,8 +45,8 @@ def run():
# Configure I/O
evt_root = simple.setupROOTOutput('RootOutput','CLICSiD_'+time.strftime('%Y-%m-%d_%H-%M'))
- #evt_lcio = simple.setupLCIOOutput('LcioOutput','CLICSiD_'+time.strftime('%Y-%m-%d_%H-%M'))
- #evt_lcio.OutputLevel = Output.ERROR
+ evt_lcio = simple.setupLCIOOutput('LcioOutput','CLICSiD_'+time.strftime('%Y-%m-%d_%H-%M'))
+ evt_lcio.OutputLevel = Output.ERROR
generator_output_level = Output.INFO
diff --git a/DDG4/include/DDG4/Geant4MonteCarloTruth.h b/DDG4/include/DDG4/Geant4MonteCarloTruth.h
index ffcb8c25ac2733ff9884127447bb8337f9f6a99c..27507ad073119d04a51411d154e563ee542b26a6 100644
--- a/DDG4/include/DDG4/Geant4MonteCarloTruth.h
+++ b/DDG4/include/DDG4/Geant4MonteCarloTruth.h
@@ -36,24 +36,12 @@ namespace DD4hep {
* \ingroup DD4HEP_SIMULATION
*/
class Geant4MonteCarloTruth {
- public:
- typedef Geant4Particle Particle;
- typedef std::map<int,Particle*> ParticleMap;
- typedef std::map<int,int> TrackEquivalents;
protected:
/// Standard constructor
Geant4MonteCarloTruth();
public:
/// Default destructor
virtual ~Geant4MonteCarloTruth();
-#ifndef __MAKECINT__
- /// Access the particle map
- virtual const ParticleMap& particles() const = 0;
- /// Access the map of track equivalents
- virtual const TrackEquivalents& equivalents() const = 0;
-#endif
- /// Access the equivalent track id (shortcut to the usage of TrackEquivalents)
- virtual int particleID(int track, bool throw_if_not_found=true) 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
@@ -71,24 +59,11 @@ namespace DD4hep {
* \ingroup DD4HEP_SIMULATION
*/
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();
-#ifndef __MAKECINT__
- /// 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; }
-#endif
- /// Access the equivalent track id (shortcut to the usage of TrackEquivalents)
- virtual int particleID(int track, bool) const { return track; }
/// 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
diff --git a/DDG4/include/DDG4/Geant4Particle.h b/DDG4/include/DDG4/Geant4Particle.h
index 278092f64f8b99e857e347dfbc81d57cdcf59a0d..35bbb66b59a473a08d65074f4cd110859e8bc7d9 100644
--- a/DDG4/include/DDG4/Geant4Particle.h
+++ b/DDG4/include/DDG4/Geant4Particle.h
@@ -58,6 +58,7 @@ namespace DD4hep {
G4PARTICLE_CREATED_TRACKER_HIT = 1<<8,
G4PARTICLE_KEEP_USER = 1<<9,
G4PARTICLE_KEEP_ALWAYS = 1<<10,
+ G4PARTICLE_FORCE_KILL = 1<<11,
// Generator status for a given particles: bit 0...3 come from LCIO, rest is internal
G4PARTICLE_GEN_EMPTY = 1<<0, // Empty line
diff --git a/DDG4/include/DDG4/Geant4ParticleHandler.h b/DDG4/include/DDG4/Geant4ParticleHandler.h
index 787b4c35fd5cd6d290f619bcbed8399ecd47995e..4c64b7a20f024e706d33bad09fb2fe2509b14a17 100644
--- a/DDG4/include/DDG4/Geant4ParticleHandler.h
+++ b/DDG4/include/DDG4/Geant4ParticleHandler.h
@@ -27,8 +27,6 @@ namespace DD4hep {
namespace Simulation {
// Forward declarations
- class Particle;
- class Geant4PrimaryMap;
class Geant4UserParticleHandler;
/// Geant4Action to collect the MC particle information.
@@ -40,6 +38,10 @@ namespace DD4hep {
* \ingroup DD4HEP_SIMULATION
*/
class Geant4ParticleHandler : public Geant4GeneratorAction, public Geant4MonteCarloTruth {
+ public:
+ typedef Geant4ParticleMap::Particle Particle;
+ typedef Geant4ParticleMap::ParticleMap ParticleMap;
+ typedef Geant4ParticleMap::TrackEquivalents TrackEquivalents;
#ifdef __MAKECINT__
public:
#else
@@ -52,21 +54,7 @@ namespace DD4hep {
public:
typedef std::vector<std::string> Processes;
- /// Functor to select particles from a integer map by the identifier.
- /**
- * \author M.Frank
- * \version 1.0
- * \ingroup DD4HEP_SIMULATION
- */
- class FindParticleByID {
- protected:
- int pid;
- public:
- FindParticleByID(int p) : pid(p) {}
- inline bool operator()(const std::pair<int,Geant4Particle*>& p) const {
- return p.second->id == pid;
- }
- };
+
protected:
/// Property: Steer printout at tracking action begin
@@ -101,8 +89,6 @@ namespace DD4hep {
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;
/// Rebase the simulated tracks, so that they fit to the generator particles
void rebaseSimulatedTracks(int base);
@@ -128,14 +114,6 @@ namespace DD4hep {
virtual void begin(const G4Track* track);
/// Post-track action callback
virtual void end(const G4Track* track);
-#ifndef __MAKECINT__
- /// 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; }
-#endif
- /// Access the equivalent track id (shortcut to the usage of TrackEquivalents)
- virtual int particleID(int track, bool throw_if_not_found=true) const;
/// Mark a Geant4 track to be kept for later MC truth analysis. Default flag: CREATED_HIT
virtual void mark(const G4Track* track);
@@ -146,6 +124,9 @@ namespace DD4hep {
/// Store a track produced in a step to be kept for later MC truth analysis
virtual void mark(const G4Step* step, int reason);
+ /// Default callback to be answered if the particle should be kept if NO user handler is installed
+ static bool defaultKeepParticle(Particle& particle);
+
};
} // End namespace Simulation
} // End namespace DD4hep
diff --git a/DDG4/include/DDG4/Geant4UserParticleHandler.h b/DDG4/include/DDG4/Geant4UserParticleHandler.h
index 05f0e83993edf574b81029c3d1e9582c9113417a..57d60c5ad080363bdbf91f11b91564f1cd924783 100644
--- a/DDG4/include/DDG4/Geant4UserParticleHandler.h
+++ b/DDG4/include/DDG4/Geant4UserParticleHandler.h
@@ -44,6 +44,11 @@ namespace DD4hep {
class Geant4UserParticleHandler : public Geant4Action {
public:
typedef Geant4Particle Particle;
+
+ protected:
+ /// Property: Energy cut below which particles are not collected, but assigned to the parent
+ double m_kinEnergyCut;
+
public:
/// Standard constructor
Geant4UserParticleHandler(Geant4Context* context, const std::string& nam);
@@ -90,10 +95,18 @@ namespace DD4hep {
* to set the reason mask to NULL in order to drop it.
* The default implementation is empty.
*
+ * If the reason mask entry is set to G4PARTICLE_FORCE_KILL
+ * or is set to NULL, the particle is ALWAYS removed
+ *
+ * The default implementation calls
+ * Geant4ParticleHandler::defaultKeepParticle(particle)
+ * Please have a look therein if it suffices your needs!
+ *
* Note: This may override all other decisions!
* Default implementation is empty.
+ *
*/
- virtual void keepParticle(Particle& particle);
+ virtual bool keepParticle(Particle& particle);
/// Callback when parent should be combined
/** Called before a particle is removed from the final record.
diff --git a/DDG4/lcio/Geant4Output2LCIO.cpp b/DDG4/lcio/Geant4Output2LCIO.cpp
index 26014ee69f2430b892254d660e5acc97bd2ccd99..7c62ff6b6a1d7f23748cdca2cc0a8075a8cf02ce 100644
--- a/DDG4/lcio/Geant4Output2LCIO.cpp
+++ b/DDG4/lcio/Geant4Output2LCIO.cpp
@@ -17,6 +17,7 @@
#include "lcio.h"
#include "IO/LCWriter.h"
#include "IMPL/LCEventImpl.h"
+#include "IMPL/LCCollectionVec.h"
/// Namespace for the AIDA detector description toolkit
namespace DD4hep {
@@ -25,6 +26,8 @@ namespace DD4hep {
/// Namespace for the Geant4 based simulation part of the AIDA detector description toolkit
namespace Simulation {
+
+ class Geant4ParticleMap;
/// Base class to output Geant4 event data to media
/**
@@ -37,6 +40,8 @@ namespace DD4hep {
lcio::LCWriter* m_file;
int m_runNo;
+ /// Data conversion interface for MC particles to LCIO format
+ lcio::LCCollectionVec* saveParticles(Geant4ParticleMap* particles);
public:
/// Standard constructor
Geant4Output2LCIO(Geant4Context* context, const std::string& nam);
@@ -75,6 +80,7 @@ namespace DD4hep {
// Framework include files
#include "DD4hep/InstanceCount.h"
+#include "DD4hep/LCDD.h"
#include "DDG4/Geant4HitCollection.h"
#include "DDG4/Geant4DataConversion.h"
#include "DDG4/Geant4Context.h"
@@ -82,6 +88,8 @@ namespace DD4hep {
#include "DDG4/Geant4Data.h"
//#include "DDG4/Geant4Output2LCIO.h"
+#include "G4ParticleDefinition.hh"
+#include "G4VProcess.hh"
#include "G4Event.hh"
#include "G4Run.hh"
@@ -89,7 +97,11 @@ namespace DD4hep {
#include "IMPL/LCEventImpl.h"
#include "IMPL/LCRunHeaderImpl.h"
#include "IMPL/LCCollectionVec.h"
-
+#include "IMPL/ClusterImpl.h"
+#include "IMPL/SimTrackerHitImpl.h"
+#include "IMPL/SimCalorimeterHitImpl.h"
+#include "IMPL/MCParticleImpl.h"
+#include "UTIL/ILDConf.h"
using namespace DD4hep::Simulation;
using namespace DD4hep;
@@ -98,7 +110,6 @@ using namespace std;
#include "DDG4/Factories.h"
DECLARE_GEANT4ACTION(Geant4Output2LCIO)
-
/// Standard constructor
Geant4Output2LCIO::Geant4Output2LCIO(Geant4Context* ctxt, const string& nam)
: Geant4OutputAction(ctxt,nam), m_file(0), m_runNo(0)
@@ -143,7 +154,6 @@ void Geant4Output2LCIO::saveRun(const G4Run* run) {
m_file->writeRunHeader(rh);
}
-
void Geant4Output2LCIO::begin(const G4Event* /* event */) {
lcio::LCEventImpl* e = new lcio::LCEventImpl;
//fg: here the event context takes ownership and
@@ -151,6 +161,103 @@ void Geant4Output2LCIO::begin(const G4Event* /* event */) {
context()->event().addExtension<lcio::LCEventImpl>( e );
}
+/// Data conversion interface for MC particles to LCIO format
+lcio::LCCollectionVec* Geant4Output2LCIO::saveParticles(Geant4ParticleMap* particles) {
+ typedef ReferenceBitMask<const int> PropertyMask;
+ typedef Geant4ParticleMap::ParticleMap ParticleMap;
+ const ParticleMap& pm = particles->particleMap;
+ size_t nparts = pm.size();
+ lcio::LCCollectionVec* lc_coll = new lcio::LCCollectionVec(lcio::LCIO::MCPARTICLE);
+ lc_coll->reserve(nparts);
+ if ( nparts > 0 ) {
+ size_t cnt = 0;
+ map<int,int> p_ids;
+ vector<const Geant4Particle*> p_part(pm.size(),0);
+ vector<MCParticleImpl*> p_lcio(pm.size(),0);
+ // First create the particles
+ for(ParticleMap::const_iterator i=pm.begin(); i!=pm.end();++i, ++cnt) {
+ int id = (*i).first;
+ const Geant4Particle* p = (*i).second;
+ PropertyMask mask(p->status);
+ const G4ParticleDefinition* def = p->definition;
+ MCParticleImpl* q = new lcio::MCParticleImpl();
+ q->setPDG(p->pdgID);
+
+ float ps_fa[3] = { p->psx/GeV, p->psy/GeV, p->psz/GeV } ;
+ q->setMomentum( ps_fa );
+
+ double vs_fa[3] = { p->vsx/mm, p->vsy/mm, p->vsz/mm } ;
+ q->setVertex( vs_fa );
+
+ double ve_fa[3] = { p->vex/mm, p->vey/mm, p->vez/mm } ;
+ q->setEndpoint( ve_fa );
+
+ q->setTime(p->time/ns);
+ q->setMass(p->mass/GeV);
+ q->setCharge(def ? def->GetPDGCharge()/3.0 : 0); // Charge(e+) = 1 !
+
+ // Set generator status
+ q->setGeneratorStatus(0);
+ if ( mask.isSet(G4PARTICLE_GEN_STABLE) ) q->setGeneratorStatus(1);
+ else if ( mask.isSet(G4PARTICLE_GEN_DECAYED) ) q->setGeneratorStatus(2);
+ else if ( mask.isSet(G4PARTICLE_GEN_DOCUMENTATION) ) q->setGeneratorStatus(3);
+
+ // Set simulation status
+ q->setSimulatorStatus( 0 ) ;
+ q->setCreatedInSimulation( mask.isSet(G4PARTICLE_SIM_CREATED) );
+ q->setBackscatter( mask.isSet(G4PARTICLE_SIM_BACKSCATTER) );
+ q->setVertexIsNotEndpointOfParent( mask.isSet(G4PARTICLE_SIM_PARENT_RADIATED) );
+ q->setDecayedInTracker( mask.isSet(G4PARTICLE_SIM_DECAY_TRACKER) );
+ q->setDecayedInCalorimeter( mask.isSet(G4PARTICLE_SIM_DECAY_CALO) );
+ q->setHasLeftDetector( mask.isSet(G4PARTICLE_SIM_LEFT_DETECTOR) );
+ q->setStopped( mask.isSet(G4PARTICLE_SIM_STOPPED) );
+ q->setOverlay( false );
+
+ //fg: if simstatus !=0 we have to set the generator status to 0:
+ if( q->getSimulatorStatus() != 0 )
+ q->setGeneratorStatus( 0 ) ;
+
+ q->setSpin(p->spin);
+ q->setColorFlow(p->colorFlow);
+
+ lc_coll->addElement(q);
+ p_ids[id] = cnt;
+ p_part[cnt] = p;
+ p_lcio[cnt] = q;
+ }
+
+ // Now establish parent-daughter relationships
+ for(size_t i=0, n=p_ids.size(); i<n; ++i) {
+ map<int,int>::iterator k;
+ const Geant4Particle* p = p_part[i];
+ MCParticleImpl* q = p_lcio[i];
+ const Geant4Particle::Particles& dau = p->daughters;
+ for(Geant4Particle::Particles::const_iterator j=dau.begin(); j!=dau.end(); ++j) {
+ int idau = *j;
+ if ( (k=p_ids.find(idau)) == p_ids.end() ) { // Error!!!
+ printout(FATAL,"Geant4Conversion","+++ Particle %d: FAILED to find daughter with ID:%d",p->id,idau);
+ continue;
+ }
+ int iqdau = (*k).second;
+ MCParticleImpl* qdau = p_lcio[iqdau];
+ qdau->addParent(q);
+ }
+ const Geant4Particle::Particles& par = p->parents;
+ for(Geant4Particle::Particles::const_iterator j=par.begin(); j!=par.end(); ++j) {
+ int ipar = *j; // A parent ID iof -1 means NO parent, because a base of 0 is perfectly leagal!
+ if ( ipar>=0 && (k=p_ids.find(ipar)) == p_ids.end() ) { // Error!!!
+ printout(FATAL,"Geant4Conversion","+++ Particle %d: FAILED to find parent with ID:%d",p->id,ipar);
+ continue;
+ }
+ int iqpar = (*k).second;
+ MCParticleImpl* qpar = p_lcio[iqpar];
+ q->addParent(qpar);
+ }
+ }
+ }
+ return lc_coll;
+}
+
/// Callback to store the Geant4 event
void Geant4Output2LCIO::saveEvent(OutputContext<G4Event>& ctxt) {
lcio::LCEventImpl* e = context()->event().extension<lcio::LCEventImpl>();
@@ -163,10 +270,7 @@ void Geant4Output2LCIO::saveEvent(OutputContext<G4Event>& ctxt) {
if ( part_map ) {
print("+++ Saving %d LCIO particles....",int(part_map->particleMap.size()));
if ( part_map->particleMap.size() > 0 ) {
- typedef pair<const Geant4Context*,const Geant4ParticleMap*> _Args;
- typedef Geant4Conversion<lcio::LCCollectionVec,_Args> _C;
- const _C& cnv = _C::converter(typeid(Geant4ParticleMap));
- lcio::LCCollectionVec* col = cnv(_Args(context(),part_map));
+ lcio::LCCollectionVec* col = saveParticles(part_map);
evt->addCollection(col,lcio::LCIO::MCPARTICLE);
}
}
diff --git a/DDG4/plugins/Geant4TCUserParticleHandler.cpp b/DDG4/plugins/Geant4TCUserParticleHandler.cpp
index 792f5e3886a878f2770de622a6f07241b9443742..28216b48643cafe22f004b7a2b0ffd670d5487f5 100644
--- a/DDG4/plugins/Geant4TCUserParticleHandler.cpp
+++ b/DDG4/plugins/Geant4TCUserParticleHandler.cpp
@@ -79,16 +79,13 @@ Geant4TCUserParticleHandler::Geant4TCUserParticleHandler(Geant4Context* context,
/// Post-track action callback
void Geant4TCUserParticleHandler::end(const G4Track* /* track */, Particle& p) {
-
-
double r_prod = sqrt(p.vsx*p.vsx + p.vsy*p.vsy);
double z_prod = fabs(p.vsz);
bool starts_in_trk_vol = ( r_prod <= m_rTracker && z_prod <= m_zTracker ) ;
// created in tracking volume but below energy cut
- if( starts_in_trk_vol && ! (p.reason&G4PARTICLE_ABOVE_ENERGY_THRESHOLD) ){
-
- p.reason = 0 ;
+ if( starts_in_trk_vol && ! (p.reason&G4PARTICLE_ABOVE_ENERGY_THRESHOLD) ) {
+ p.reason = 0;
return;
}
@@ -98,12 +95,10 @@ void Geant4TCUserParticleHandler::end(const G4Track* /* track */, Particle& p)
// created and ended in calo
if( !starts_in_trk_vol && !ends_in_trk_vol ){
-
p.reason = 0;
return ;
}
//fg: backscatter ??
// if( !starts_in_trk_vol && ends_in_trk_vol ){ keep ? }
-
}
diff --git a/DDG4/src/Geant4InputHandling.cpp b/DDG4/src/Geant4InputHandling.cpp
index 0a2180100431e02634e6ad8dfb0105e20b838216..baf1a422f0ef7d3ed1f5b6f936dd62bb69925d21 100644
--- a/DDG4/src/Geant4InputHandling.cpp
+++ b/DDG4/src/Geant4InputHandling.cpp
@@ -312,6 +312,10 @@ int DD4hep::Simulation::generatePrimaries(const Geant4Action* caller,
v->x/mm,v->y/mm,v->z/mm,v->time/ns);
for(Geant4Vertex::Particles::const_iterator ip=v->out.begin(); ip!=v->out.end(); ++ip) {
Geant4ParticleHandle p = pm[*ip];
+ if ( p->daughters.size() > 0 ) {
+ PropertyMask mask(p->reason);
+ mask.set(G4PARTICLE_HAS_SECONDARIES);
+ }
if ( p->parents.size() == 0 ) {
Primaries relevant = getRelevant(visited,prim,pm,p);
for(Primaries::const_iterator j=relevant.begin(); j!= relevant.end(); ++j) {
diff --git a/DDG4/src/Geant4Particle.cpp b/DDG4/src/Geant4Particle.cpp
index 3dea8117b55050a27a6b1293f25b7a1e2844ff14..a529c210a4773843275f5be13448489cd9ef71d2 100644
--- a/DDG4/src/Geant4Particle.cpp
+++ b/DDG4/src/Geant4Particle.cpp
@@ -339,19 +339,9 @@ void Geant4ParticleMap::adopt(ParticleMap& pm, TrackEquivalents& equiv) {
/// Access the equivalent track id (shortcut to the usage of TrackEquivalents)
int Geant4ParticleMap::particleID(int g4_id, bool) const {
- int equiv_id = g4_id;
- if ( g4_id != 0 ) {
- ParticleMap::const_iterator ipar;
- while( (ipar=particleMap.find(equiv_id)) == particleMap.end() ) {
- TrackEquivalents::const_iterator iequiv = equivalentTracks.find(equiv_id);
- equiv_id = (iequiv == equivalentTracks.end()) ? -1 : (*iequiv).second;
- if ( equiv_id < 0 ) {
- printout(INFO,"Geant4ParticleMap",
- "+++ No Equivalent particle for track:%d last known is:%d",
- g4_id,equiv_id);
- break;
- }
- }
- }
- return equiv_id;
+ TrackEquivalents::const_iterator iequiv = equivalentTracks.find(g4_id);
+ if ( iequiv != equivalentTracks.end() ) return (*iequiv).second;
+ printout(ERROR,"Geant4ParticleMap","+++ No Equivalent particle for track:%d."
+ " Monte Carlo truth record looks broken!",g4_id);
+ return -1;
}
diff --git a/DDG4/src/Geant4ParticleHandler.cpp b/DDG4/src/Geant4ParticleHandler.cpp
index 0c642a7c06b31c5b526a748a8223a13facb466bb..90a700d9cb9d19fc485067a222bbddd64a2e2076 100644
--- a/DDG4/src/Geant4ParticleHandler.cpp
+++ b/DDG4/src/Geant4ParticleHandler.cpp
@@ -202,8 +202,6 @@ void Geant4ParticleHandler::begin(const G4Track* track) {
G4PrimaryParticle* prim = primary(h.id(),context()->event().event());
Particle* prim_part = 0;
-
-
if ( prim ) {
Geant4PrimaryMap::Primaries::const_iterator iprim = m_primaryMap->primaryMap.find(prim);
if ( iprim == m_primaryMap->primaryMap.end() ) {
@@ -437,9 +435,20 @@ void Geant4ParticleHandler::rebaseSimulatedTracks(int ) {
}
if ( ipar != m_particleMap.end() ) {
equivalents[(*i).first] = (*ipar).second->id; // requires (1) !
+ Particle* p = (*ipar).second;
+ const G4ParticleDefinition* def = p->definition;
+ int pdg = int(fabs(def->GetPDGEncoding())+0.1);
+ if ( pdg<36 && !(pdg > 10 && pdg < 17) && pdg != 22 ) {
+ error("+++ ERROR: Geant4 particle for track:%d last known is:%d -- is gluon or quark!",(*i).second,g4_equiv);
+ }
+ pdg = int(fabs(p->pdgID)+0.1);
+ if ( pdg<36 && !(pdg > 10 && pdg < 17) && pdg != 22 ) {
+ error("+++ ERROR(2): Geant4 particle for track:%d last known is:%d -- is gluon or quark!",(*i).second,g4_equiv);
+ }
}
- else
+ else {
error("+++ No Equivalent particle for track:%d last known is:%d",(*i).second,g4_equiv);
+ }
}
// (3) Compute the particle's parents and daughters.
@@ -465,101 +474,101 @@ void Geant4ParticleHandler::rebaseSimulatedTracks(int ) {
m_particleMap = finalParticles;
}
+/// Default callback to be answered if the particle should be kept if NO user handler is installed
+bool Geant4ParticleHandler::defaultKeepParticle(Particle& particle) {
+ PropertyMask mask(particle.reason);
+ 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);
+
+ /// Remove this track if it has not created a hit and the energy is below threshold
+ if ( mask.isNull() || (secondaries && low_energy && !hits_produced) ) {
+ return true;
+ }
+ /// Remove this track if the energy is below threshold. Reassign hits to parent.
+ else if ( !hits_produced && low_energy ) {
+ return true;
+ }
+ /// Remove this track if the origine is in the calorimeter. Reassign hits to parent.
+ else if ( !tracker_track && calo_track && low_energy ) {
+ return true;
+ }
+ else {
+ //printout(INFO,name(),"+++ Track: %d should be kept for no obvious reason....",id);
+ }
+ return false;
+}
+
/// Clean the monte carlo record. Remove all unwanted stuff.
/// This is the core of the object executed at the end of each event action.
int Geant4ParticleHandler::recombineParents() {
- int break_trackID = 38;
set<int> remove;
+
/// Need to start from BACK, to clean first the latest produced stuff.
- /// Otherwise the daughter list of the earlier produced tracks would not be empty!
for(ParticleMap::reverse_iterator i=m_particleMap.rbegin(); i!=m_particleMap.rend(); ++i) {
- int g4_id = (*i).first;
Particle* p = (*i).second;
- set<int>& daughters = p->daughters;
+ PropertyMask mask(p->reason);
+ bool remove_me = false;
// Allow the user to force the particle handling either by
// or the reason mask with G4PARTICLE_KEEP_USER or
// to set the reason mask to NULL in order to drop it.
+ //
+ // If the mask entry is set to G4PARTICLE_FORCE_KILL
+ // or is set to NULL, the particle is ALWAYS removed
+ //
// Note: This may override all other decisions!
- if ( m_userHandler ) {
- m_userHandler->keepParticle(*p);
- }
+ remove_me = m_userHandler ? m_userHandler->keepParticle(*p) : defaultKeepParticle(*p);
- if ( daughters.size() == 0 ) {
- PropertyMask mask(p->reason);
- int id = p->id;
- 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);
- bool remove_me = false;
-
- if ( id == break_trackID ) { // Used for debugging to set break point
- remove_me = false;
- }
-
- if ( mask.isSet(G4PARTICLE_KEEP_USER) ) {
- /// If user decides it must be kept, it MUST be kept!
- mask.set(G4PARTICLE_KEEP_USER);
- continue;
- }
- else if ( mask.isSet(G4PARTICLE_PRIMARY) ) {
- /// Primary particles MUST be kept!
- continue;
- }
- else if ( mask.isSet(G4PARTICLE_KEEP_ALWAYS) ) {
- continue;
- }
- else if ( keep_parent ) {
- //continue;
- }
- else if ( keep_process ) {
- ParticleMap::iterator ip = m_particleMap.find(p->g4Parent);
- if ( ip != m_particleMap.end() ) {
- Particle* parent_part = (*ip).second;
- PropertyMask parent_mask(parent_part->reason);
- if ( parent_mask.isSet(G4PARTICLE_ABOVE_ENERGY_THRESHOLD) ) {
- parent_mask.set(G4PARTICLE_KEEP_PARENT);
- continue;
- }
+ // Now look at the property mask of the particle
+ if ( mask.isNull() || mask.isSet(G4PARTICLE_FORCE_KILL) ) {
+ remove_me = true;
+ }
+ else if ( mask.isSet(G4PARTICLE_KEEP_USER) ) {
+ /// If user decides it must be kept, it MUST be kept!
+ mask.set(G4PARTICLE_KEEP_USER);
+ continue;
+ }
+ else if ( mask.isSet(G4PARTICLE_PRIMARY) ) {
+ /// Primary particles MUST be kept!
+ continue;
+ }
+ else if ( mask.isSet(G4PARTICLE_KEEP_ALWAYS) ) {
+ continue;
+ }
+ else if ( mask.isSet(G4PARTICLE_KEEP_PARENT) ) {
+ //continue;
+ }
+ else if ( mask.isSet(G4PARTICLE_KEEP_PROCESS) ) {
+ ParticleMap::iterator ip = m_particleMap.find(p->g4Parent);
+ if ( ip != m_particleMap.end() ) {
+ Particle* parent_part = (*ip).second;
+ PropertyMask parent_mask(parent_part->reason);
+ if ( parent_mask.isSet(G4PARTICLE_ABOVE_ENERGY_THRESHOLD) ) {
+ parent_mask.set(G4PARTICLE_KEEP_PARENT);
+ continue;
}
- // Low energy stuff. Remove it. Reassign to parent.
- //remove_me = true;
- }
-
- /// Remove this track if it has not created a hit and the energy is below threshold
- if ( mask.isNull() || (secondaries && low_energy && !hits_produced) ) {
- remove_me = true;
- }
- /// Remove this track if the energy is below threshold. Reassign hits to parent.
- else if ( !hits_produced && low_energy ) {
- remove_me = true;
- }
- /// Remove this track if the origine is in the calorimeter. Reassign hits to parent.
- else if ( !tracker_track && calo_track && low_energy ) {
- remove_me = true;
- }
- else {
- //printout(INFO,name(),"+++ Track: %d should be kept for no obvious reason....",id);
}
+ // Low energy stuff. Remove it. Reassign to parent.
+ //remove_me = true;
+ }
- /// Remove this track from the list and also do the cleanup in the parent's children list
- if ( remove_me ) {
- ParticleMap::iterator ip = m_particleMap.find(p->g4Parent);
- remove.insert(g4_id);
- m_equivalentTracks[g4_id] = p->g4Parent;
- if ( ip != m_particleMap.end() ) {
- Particle* parent_part = (*ip).second;
- PropertyMask(parent_part->reason).set(mask.value());
- parent_part->steps += p->steps;
- parent_part->secondaries += p->secondaries;
- /// Update of the particle using the user handler
- if ( m_userHandler ) {
- m_userHandler->combine(*p, *parent_part);
- }
+ /// Remove this track from the list and also do the cleanup in the parent's children list
+ if ( remove_me ) {
+ int g4_id = (*i).first;
+ ParticleMap::iterator ip = m_particleMap.find(p->g4Parent);
+ remove.insert(g4_id);
+ m_equivalentTracks[g4_id] = p->g4Parent;
+ if ( ip != m_particleMap.end() ) {
+ Particle* parent_part = (*ip).second;
+ PropertyMask(parent_part->reason).set(mask.value());
+ parent_part->steps += p->steps;
+ parent_part->secondaries += p->secondaries;
+ /// Update of the particle using the user handler
+ if ( m_userHandler ) {
+ m_userHandler->combine(*p, *parent_part);
}
}
}
@@ -595,12 +604,17 @@ void Geant4ParticleHandler::checkConsistency() const {
// We assume that particles from the generator have consistent parents
// For all other particles except the primaries, the parent must be contained in the record.
if ( !mask.isSet(G4PARTICLE_PRIMARY) && !status.anySet(G4PARTICLE_GEN_GENERATOR) ) {
- int parent_id = equivalentTrack(p->g4Parent);
- bool in_map = (j=m_particleMap.find(parent_id)) != m_particleMap.end();
- bool in_parent_list = p->parents.find(parent_id) != p->parents.end();
- char parent_list[1024];
- parent_list[0] = 0;
+ TrackEquivalents::const_iterator eq_it = m_equivalentTracks.find(p->g4Parent);
+ bool in_map = false, in_parent_list = false;
+ int parent_id = -1;
+ if ( eq_it != m_equivalentTracks.end() ) {
+ parent_id = (*eq_it).second;
+ in_map = (j=m_particleMap.find(parent_id)) != m_particleMap.end();
+ in_parent_list = p->parents.find(parent_id) != p->parents.end();
+ }
if ( !in_map || !in_parent_list ) {
+ char parent_list[1024];
+ parent_list[0] = 0;
++num_errors;
for(set<int>::const_iterator ip=p->parents.begin(); ip!=p->parents.end();++ip)
::snprintf(parent_list+strlen(parent_list),sizeof(parent_list)-strlen(parent_list),"%d ",*ip);
@@ -610,32 +624,7 @@ void Geant4ParticleHandler::checkConsistency() const {
}
}
- /// 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(), iend=m_equivalentTracks.end(); i!=iend; ++i) {
- int g4_id = (*i).first;
- int equiv_id = equivalentTrack(g4_id);
- if ( equiv_id < 0 ) {
- ++num_errors;
- }
- }
-
if ( num_errors > 0 ) {
except("+++ Consistency check failed. Found %d problems.",num_errors);
}
}
-
-/// 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 ) {
- TrackEquivalents::const_iterator iequiv = m_equivalentTracks.find(equiv_id);
- if ( iequiv != m_equivalentTracks.end() ) return (*iequiv).second;
- }
- return -1;
-}
-
-/// Access the equivalent track id (shortcut to the usage of TrackEquivalents)
-int Geant4ParticleHandler::particleID(int track, bool) const {
- return equivalentTrack(track);
-}
diff --git a/DDG4/src/Geant4UserParticleHandler.cpp b/DDG4/src/Geant4UserParticleHandler.cpp
index 221756edf945238c570c9d38641e4390a723ede1..539242a4816fc520ae2b70f0427533bda5461a71 100644
--- a/DDG4/src/Geant4UserParticleHandler.cpp
+++ b/DDG4/src/Geant4UserParticleHandler.cpp
@@ -9,7 +9,9 @@
// Framework include files
#include "DD4hep/Printout.h"
#include "DD4hep/InstanceCount.h"
+#include "DDG4/Geant4ParticleHandler.h"
#include "DDG4/Geant4UserParticleHandler.h"
+#include "CLHEP/Units/SystemOfUnits.h"
using namespace DD4hep::Simulation;
@@ -18,6 +20,7 @@ Geant4UserParticleHandler::Geant4UserParticleHandler(Geant4Context* context, con
: Geant4Action(context,nam)
{
InstanceCount::increment(this);
+ declareProperty("MinimalKineticEnergy",m_kinEnergyCut = 100e0*CLHEP::MeV);
m_needsControl = true;
}
@@ -57,5 +60,6 @@ void Geant4UserParticleHandler::combine(Particle& /* to_be_deleted */, Particle&
}
/// Callback to be answered if the particle MUST be kept during recombination step
-void Geant4UserParticleHandler::keepParticle(Particle& /* particle */) {
+bool Geant4UserParticleHandler::keepParticle(Particle& particle) {
+ return Geant4ParticleHandler::defaultKeepParticle(particle);
}
diff --git a/cmake/Doxyfile.in b/cmake/Doxyfile.in
index 22e0c53e63fc1f5978983e7b467d6d3dc6db1bc2..20d5c22443c9752abc8983e8ba36c799ddb22d14 100644
--- a/cmake/Doxyfile.in
+++ b/cmake/Doxyfile.in
@@ -650,15 +650,19 @@ WARN_LOGFILE =
# with spaces.
INPUT = @CMAKE_CURRENT_SOURCE_DIR@/doc/doxygen/DD4hepMainpage.dox
-INPUT += @CMAKE_CURRENT_SOURCE_DIR@/doc/doxygen/DD4hepFigures.dox
INPUT += @CMAKE_CURRENT_SOURCE_DIR@/DDCore/src
INPUT += @CMAKE_CURRENT_SOURCE_DIR@/DDCore/src/XML
INPUT += @CMAKE_CURRENT_SOURCE_DIR@/DDCore/include/DD4hep
+INPUT += @CMAKE_CURRENT_SOURCE_DIR@/DDCore/include/DD4hep/objects
INPUT += @CMAKE_CURRENT_SOURCE_DIR@/DDCore/include/XML
INPUT += @CMAKE_CURRENT_SOURCE_DIR@/DDG4/include/DDG4
INPUT += @CMAKE_CURRENT_SOURCE_DIR@/DDG4/src
INPUT += @CMAKE_CURRENT_SOURCE_DIR@/DDG4/plugins
INPUT += @CMAKE_CURRENT_SOURCE_DIR@/DDG4/lcio
+INPUT += @CMAKE_CURRENT_SOURCE_DIR@/DDAlign/include/DDAlign
+INPUT += @CMAKE_CURRENT_SOURCE_DIR@/DDAlign/src
+INPUT += @CMAKE_CURRENT_SOURCE_DIR@/DDCond/include/DDCond
+INPUT += @CMAKE_CURRENT_SOURCE_DIR@/DDCond/src
INPUT += @CMAKE_CURRENT_SOURCE_DIR@/DDEve/include/DDEve
INPUT += @CMAKE_CURRENT_SOURCE_DIR@/DDEve/src
INPUT += @CMAKE_CURRENT_SOURCE_DIR@/DDEve/lcio
@@ -669,6 +673,10 @@ INPUT += @CMAKE_CURRENT_SOURCE_DIR@/DDSense/include/DDSense
INPUT += @CMAKE_CURRENT_SOURCE_DIR@/DDRec/src
INPUT += @CMAKE_CURRENT_SOURCE_DIR@/DDRec/include/DDRec
INPUT += @CMAKE_CURRENT_SOURCE_DIR@/DDSurfaces/include/DDSurfaces
+INPUT += @CMAKE_CURRENT_SOURCE_DIR@/doc/doxygen/Geant4Classes.h
+INPUT += @CMAKE_CURRENT_SOURCE_DIR@/doc/doxygen/ROOTClasses.h
+INPUT += @CMAKE_CURRENT_SOURCE_DIR@/doc/doxygen/TiXMLClasses.h
+INPUT += @CMAKE_CURRENT_SOURCE_DIR@/doc/doxygen/DD4hepGroups.h
#INPUT += @CMAKE_CURRENT_SOURCE_DIR@//src
#INPUT += @CMAKE_CURRENT_SOURCE_DIR@//include