diff --git a/DDG4/include/DDG4/Geant4Context.h b/DDG4/include/DDG4/Geant4Context.h index 355a69fbb6291a4e1ab2174fff571861d8d20fe2..cd59f4aae2a7a93e68c7cf6b1b671203069878e2 100644 --- a/DDG4/include/DDG4/Geant4Context.h +++ b/DDG4/include/DDG4/Geant4Context.h @@ -47,9 +47,6 @@ namespace DD4hep { class Geant4StackingActionSequence; class Geant4GeneratorActionSequence; class Geant4SensDetSequences; - class Geant4MonteCarloTruth; - class Geant4MonteCarloRecordManager; - /// User run context for DDG4 /** diff --git a/DDG4/include/DDG4/Geant4Data.h b/DDG4/include/DDG4/Geant4Data.h index 9cefc8c9978d1fc083eefc52a6579560f05e7e34..3c30784e52de63e16ad6a87a710575e78aede5d8 100644 --- a/DDG4/include/DDG4/Geant4Data.h +++ b/DDG4/include/DDG4/Geant4Data.h @@ -99,14 +99,14 @@ namespace DD4hep { /// 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_CREATED_HIT = 1<<1, + G4PARTICLE_PRIMARY = 1<<2, + G4PARTICLE_HAS_SECONDARIES = 1<<3, + G4PARTICLE_ABOVE_ENERGY_THRESHOLD = 1<<4, + G4PARTICLE_KEEP_PROCESS = 1<<5, + G4PARTICLE_KEEP_PARENT = 1<<6, + G4PARTICLE_CREATED_CALORIMETER_HIT = 1<<7, + G4PARTICLE_CREATED_TRACKER_HIT = 1<<8, G4PARTICLE_LAST_NOTHING = 1<<31 }; diff --git a/DDG4/include/DDG4/Geant4Kernel.h b/DDG4/include/DDG4/Geant4Kernel.h index 7831c71f27a98993b37a1b689c5f2d90180851ce..600027e5d4d1257c532a4b88fa5f91dc68e72d12 100644 --- a/DDG4/include/DDG4/Geant4Kernel.h +++ b/DDG4/include/DDG4/Geant4Kernel.h @@ -53,8 +53,6 @@ namespace DD4hep { class Geant4PhysicsListActionSequence; class Geant4SensDetActionSequence; class Geant4SensDetSequences; - class Geant4MonteCarloTruth; - class Geant4MonteCarloRecordManager; /** @class Invoke Geant4Kernel.h DDG4/Geant4Kernel.h * @@ -92,10 +90,6 @@ namespace DD4hep { Geant4SensDetSequences* m_sensDetActions; /// Reference to the geant4 physics list Geant4PhysicsListActionSequence* m_physicsList; - /// Reference to track persistency manager - Geant4MonteCarloTruth* m_mcTruthMgr; - /// Reference to MC record manager - Geant4MonteCarloRecordManager* m_mcRecordMgr; /// Reference to Geant4 track manager G4TrackingManager* m_trackMgr; @@ -297,12 +291,6 @@ namespace DD4hep { Geant4PhysicsListActionSequence& physicsList() { return *physicsList(true); } -#if 0 - /// Access to the Track Manager from the kernel object - Geant4MonteCarloTruth* mcTruthMgr(bool throw_exception=true); - /// Access to the MC record manager from the kernel object (if instantiated!) - Geant4MonteCarloRecordManager* mcRecordMgr(bool throw_exception=true); -#endif /// Construct detector geometry using lcdd plugin void loadGeometry(const std::string& compact_file); /// Run the simulation diff --git a/DDG4/include/DDG4/Geant4MonteCarloRecordManager.h b/DDG4/include/DDG4/Geant4MonteCarloRecordManager.h deleted file mode 100644 index 39e939788c922c9105285e213e60bc6531f30a2d..0000000000000000000000000000000000000000 --- a/DDG4/include/DDG4/Geant4MonteCarloRecordManager.h +++ /dev/null @@ -1,51 +0,0 @@ -// $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 -#ifndef DD4HEP_DDG4_GEANT4MONTECARLORECORDMANAGER_H -#define DD4HEP_DDG4_GEANT4MONTECARLORECORDMANAGER_H - -// Framework include files -#include "DDG4/Geant4TrackPersistency.h" - -/* - * DD4hep namespace declaration - */ -namespace DD4hep { - - /* - * Simulation namespace declaration - */ - namespace Simulation { - - // Forward declarations - - /** @class Invoke Geant4TrackPersistency.h DDG4/Geant4TrackPersistency.h - * - * Default base class for all geant 4 actions and derivates thereof. - * - * @author M.Frank - * @version 1.0 - */ - class Geant4MonteCarloRecordManager : public Geant4GeneratorAction { - public: - /// Flag to indicate if the track information should be collected - bool m_collectInfo; - public: - /// Standard constructor - Geant4MonteCarloRecordManager(Geant4Context* context, const std::string& nam); - /// Default destructor - virtual ~Geant4MonteCarloRecordManager(); - /// Event generation action callback - virtual void operator()(G4Event* event); - /// Save G4Track data - virtual void save(const Geant4TrackPersistency::TrackInfo& track); - }; - } // End namespace Simulation -} // End namespace DD4hep -#endif // DD4HEP_DDG4_GEANT4MONTECARLORECORDMANAGER_H diff --git a/DDG4/include/DDG4/Geant4TrackPersistency.h b/DDG4/include/DDG4/Geant4TrackPersistency.h deleted file mode 100644 index 4d050b60c3aa1912811204b85a8284d4f9c6fbba..0000000000000000000000000000000000000000 --- a/DDG4/include/DDG4/Geant4TrackPersistency.h +++ /dev/null @@ -1,141 +0,0 @@ -// $Id: Geant4Hits.h 513 2013-04-05 14:31:53Z gaede $ -//==================================================================== -// AIDA Detector description implementation -//-------------------------------------------------------------------- -// -// Author : M.Frank -// -//==================================================================== -#ifndef DD4HEP_DDG4_GEANT4TRACKPERSISTENCY_H -#define DD4HEP_DDG4_GEANT4TRACKPERSISTENCY_H - -// Framework include files -#include "DDG4/Geant4Action.h" -#include "DDG4/Geant4EventAction.h" -#include "DDG4/Geant4GeneratorAction.h" -#include "DDG4/Geant4MonteCarloTruth.h" -#include "Math/PxPyPzE4D.h" -#include "G4VUserTrackInformation.hh" - -// C/C++ include files -#include <map> -#include <string> -#include <typeinfo> - -// Forward declarations -class G4Step; -class G4Track; -class G4Event; -class G4SteppingManager; - -/* - * DD4hep namespace declaration - */ -namespace DD4hep { - - /* - * Simulation namespace declaration - */ - namespace Simulation { - - // Forward declarations - - /** @class Geant4TrackPersistency Geant4TrackPersistency.h DDG4/Geant4TrackPersistency.h - * - * Default base class for all geant 4 actions and derivates thereof. - * - * @author M.Frank - * @version 1.0 - */ - class Geant4TrackPersistency - : public Geant4GeneratorAction, - public Geant4MonteCarloTruth - { - public: - typedef ROOT::Math::PxPyPzE4D<double> FourMomentum; - typedef std::map<int,void*> TrackMap; - - /** @class Geant4TrackPersistency::TrackInfo Geant4TrackPersistency.h DDG4/Geant4TrackPersistency.h - * - * - * - * @author M.Frank - * @version 1.0 - */ - class TrackInfo : public G4VUserTrackInformation { - public: - /// Pointer to self - Geant4TrackPersistency* manager; - /// Pointer to the track - const G4Track* track; - /// Pointer to MC tracking history - void* info; - /// Flag to store the track - bool store; - /// Flag to know if this track created a Geant4 hit - bool createdHit; - /// Initial 4-momentum at the beginning of the tracking action - FourMomentum initialMomentum; - double initialEnergy; - - public: - /// Standard constructor - TrackInfo(); - /// Default destructor - virtual ~TrackInfo(); - /// Clear all stored information for this track - void set(const G4Track* trk, void* opt); - }; - 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 - bool m_printStartTracking; - /// Property: Steer printout at tracking action end - bool m_printEndTracking; - - public: - TrackMap m_trackMap; - - public: - /// Standard constructor - Geant4TrackPersistency(Geant4Context* context, const std::string& nam); - /// Default destructor - virtual ~Geant4TrackPersistency(); - /// Access the Geant4 tracking manager. Only use between tracking pre- and post action - G4TrackingManager* trackMgr() const { return m_context.trackMgr(); } - /// 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); - /// 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, 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, int reason); - - }; - } // End namespace Simulation -} // End namespace DD4hep - -#endif // DD4HEP_DDG4_GEANT4TRACKPERSISTENCY_H diff --git a/DDG4/plugins/Geant4Factories.cpp b/DDG4/plugins/Geant4Factories.cpp index 169105fb561890a0ccc3dc72eba53fbb44cf4c87..4858c64fbf3b3aad52416a3c43e465236aec5f71 100644 --- a/DDG4/plugins/Geant4Factories.cpp +++ b/DDG4/plugins/Geant4Factories.cpp @@ -40,12 +40,6 @@ DECLARE_GEANT4ACTION(Geant4UIManager) #include "DDG4/Geant4MonteCarloTruth.h" DECLARE_GEANT4ACTION(Geant4DummyTruthHandler) -#include "DDG4/Geant4TrackPersistency.h" -DECLARE_GEANT4ACTION(Geant4TrackPersistency) - -#include "DDG4/Geant4MonteCarloRecordManager.h" -DECLARE_GEANT4ACTION(Geant4MonteCarloRecordManager) - #include "DDG4/Geant4ParticleHandler.h" DECLARE_GEANT4ACTION(Geant4ParticleHandler) diff --git a/DDG4/src/Geant4Kernel.cpp b/DDG4/src/Geant4Kernel.cpp index 28335656f753f59514c5a4b7909ec1dee2e570df..fe4b9470e9619f19f44315b1987fa169f16767a8 100644 --- a/DDG4/src/Geant4Kernel.cpp +++ b/DDG4/src/Geant4Kernel.cpp @@ -24,8 +24,6 @@ #include "DDG4/Geant4StackingAction.h" #include "DDG4/Geant4GeneratorAction.h" #include "DDG4/Geant4SensDetAction.h" -#include "DDG4/Geant4MonteCarloRecordManager.h" -#include "DDG4/Geant4TrackPersistency.h" // Geant4 include files #include "G4RunManager.hh" @@ -69,8 +67,7 @@ Geant4ActionPhase& Geant4Kernel::PhaseSelector::operator[](const std::string& na Geant4Kernel::Geant4Kernel(LCDD& lcdd) : m_runManager(0), m_generatorAction(0), m_runAction(0), m_eventAction(0), m_trackingAction(0), m_steppingAction(0), m_stackingAction(0), m_sensDetActions(0), - m_physicsList(0), m_mcTruthMgr(0), m_mcRecordMgr(0), - m_lcdd(lcdd), phase(this) { + m_physicsList(0), m_lcdd(lcdd), phase(this) { #if 0 registerSequence(m_runAction, "RunAction"); registerSequence(m_eventAction, "EventAction"); @@ -97,10 +94,6 @@ Geant4Kernel::~Geant4Kernel() { for_each(m_globalFilters.begin(), m_globalFilters.end(), releaseObjects(m_globalFilters)); for_each(m_globalActions.begin(), m_globalActions.end(), releaseObjects(m_globalActions)); deletePtr (m_runManager); - m_mcTruthMgr = 0; - m_mcRecordMgr = 0; // These are already released by the global actions.... - //deletePtr (m_mcTruthMgr); - //releasePtr (m_mcRecordMgr); releasePtr (m_physicsList); releasePtr (m_stackingAction); releasePtr (m_steppingAction); @@ -210,10 +203,6 @@ void Geant4Kernel::terminate() { for_each(m_globalFilters.begin(), m_globalFilters.end(), releaseObjects(m_globalFilters)); for_each(m_globalActions.begin(), m_globalActions.end(), releaseObjects(m_globalActions)); deletePtr (m_runManager); - m_mcTruthMgr = 0; - m_mcRecordMgr = 0; // These are already released by the global actions.... - //deletePtr (m_mcTruthMgr); - //releasePtr (m_mcRecordMgr); releasePtr (m_physicsList); releasePtr (m_stackingAction); releasePtr (m_steppingAction); @@ -412,35 +401,3 @@ Geant4PhysicsListActionSequence* Geant4Kernel::physicsList(bool create) { registerSequence(m_physicsList, "PhysicsList"); return m_physicsList; } - -#if 0 -/// Access to the Track Manager from the kernel object -Geant4MonteCarloTruth* Geant4Kernel::mcTruthMgr(bool throw_exception) { - if ( m_mcTruthMgr ) return m_mcTruthMgr; - // If not present, check if the action is registered. - Geant4Action* a = globalAction("MonteCarloTruthHandler",false); - if ( 0 != a ) { - m_mcTruthMgr = dynamic_cast<Geant4MonteCarloTruth*>(a); - if ( m_mcTruthMgr ) return m_mcTruthMgr; - } - if ( !throw_exception ) return 0; - // No action registered to handle monte carlo truth. This is fatal - throw runtime_error(format("Geant4Kernel", "DDG4: No Geant4MonteCarloTruth defined. " - "Geant4 monte carlo information cannot be saved!")); -} - -/// Access to the MC record manager from the kernel object -Geant4MonteCarloRecordManager* Geant4Kernel::mcRecordMgr(bool throw_exception) { - if ( m_mcRecordMgr ) return m_mcRecordMgr; - // If not present, check if the action is registered. - Geant4Action* a = globalAction("MonteCarloRecordManager",false); - if ( 0 != a ) { - m_mcRecordMgr = dynamic_cast<Geant4MonteCarloRecordManager*>(a); - if ( m_mcRecordMgr ) return m_mcRecordMgr; - } - if ( !throw_exception ) return 0; - // No action registered to save tracks. This is fatal - throw runtime_error(format("Geant4Kernel", "DDG4: No MonteCarloRecordManager defined. " - "Geant4 tracks cannot be saved!")); -} -#endif diff --git a/DDG4/src/Geant4MonteCarloRecordManager.cpp b/DDG4/src/Geant4MonteCarloRecordManager.cpp deleted file mode 100644 index 8ddb6023d294725208a232c3c0cb3a03500765f8..0000000000000000000000000000000000000000 --- a/DDG4/src/Geant4MonteCarloRecordManager.cpp +++ /dev/null @@ -1,60 +0,0 @@ -// $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 "DDG4/Geant4TrackPersistency.h" -#include "DDG4/Geant4TrackingAction.h" -#include "DDG4/Geant4MonteCarloRecordManager.h" -#include "DDG4/Geant4TrackHandler.h" -#include "DDG4/Geant4StepHandler.h" - -#include "G4TrackingManager.hh" -#include "G4Track.hh" - -using namespace DD4hep::Simulation; - -/// Standard constructor -Geant4MonteCarloRecordManager::Geant4MonteCarloRecordManager(Geant4Context* context, const std::string& nam) - : Geant4GeneratorAction(context,nam) -{ - declareProperty("Collect",m_collectInfo=true); -} - -/// Default destructor -Geant4MonteCarloRecordManager::~Geant4MonteCarloRecordManager() { -} - -/// Event generation action callback -void Geant4MonteCarloRecordManager::operator()(G4Event* ) { - printout(INFO,name(),"+++ Add EVENT extension of type Geant4MonteCarloRecordManager....."); - context()->event().addExtension(this, typeid(*this),0); -} - -/// Save G4Track data -void Geant4MonteCarloRecordManager::save(const Geant4TrackPersistency::TrackInfo& trk) { - if ( m_collectInfo ) { - Geant4TrackHandler track(trk.track); -#if 0 - // const G4Step* step = track.step(); - std::string proc = track.creatorName(); - float ene = trk.initialEnergy; - float loss = ene-track.energy(); - G4ThreeVector pos = track.position(); - G4ThreeVector org = track.vertex(); - - info("Track ID=%4d Parent:%4d E:%7.2f MeV Loss:%7.2f MeV " - "%-6s %-8s L:%7.2f Org:(%8.2f %8.2f %8.2f) Pos:(%8.2f %8.2f %8.2f)", - track.id(), track.parent(), ene, loss, - track.name().c_str(), ("["+proc+"]").c_str(), track.length(), - org.x(),org.y(),org.z(), - pos.x(),pos.y(),pos.z()); -#endif - } -} - diff --git a/DDG4/src/Geant4ParticleHandler.cpp b/DDG4/src/Geant4ParticleHandler.cpp index 033be704e85d8aea32ec84f4bd74dbccdd76f02d..e268e2d88e9521e06dc3faad364b809114f19a60 100644 --- a/DDG4/src/Geant4ParticleHandler.cpp +++ b/DDG4/src/Geant4ParticleHandler.cpp @@ -15,17 +15,16 @@ #include "DDG4/Geant4EventAction.h" #include "DDG4/Geant4TrackingAction.h" #include "DDG4/Geant4SteppingAction.h" -#include "DDG4/Geant4MonteCarloRecordManager.h" #include "DDG4/Geant4ParticleHandler.h" -#include "G4TrackingManager.hh" -#include "G4ParticleDefinition.hh" -#include "G4Track.hh" #include "G4Step.hh" +#include "G4Track.hh" +#include "G4Event.hh" #include "G4TrackStatus.hh" -#include "G4PrimaryParticle.hh" #include "G4PrimaryVertex.hh" -#include "G4Event.hh" +#include "G4PrimaryParticle.hh" +#include "G4TrackingManager.hh" +#include "G4ParticleDefinition.hh" #include "CLHEP/Units/SystemOfUnits.h" #include <set> @@ -159,6 +158,7 @@ void Geant4ParticleHandler::begin(const G4Track* track) { Geant4TrackHandler h(track); double kine = h.kineticEnergy(); G4ThreeVector m = track->GetMomentum(); + const G4ThreeVector& v = h.vertex(); // Extract here all the information from the start of the tracking action // which we will need later to create a proper MC particle. @@ -170,9 +170,9 @@ void Geant4ParticleHandler::begin(const G4Track* track) { 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.vsx = v.x(); + m_currTrack.vsy = v.y(); + m_currTrack.vsz = v.z(); m_currTrack.vex = 0.0; m_currTrack.vey = 0.0; m_currTrack.vez = 0.0; @@ -182,6 +182,7 @@ void Geant4ParticleHandler::begin(const G4Track* track) { m_currTrack.pex = 0.0; m_currTrack.pey = 0.0; m_currTrack.pez = 0.0; + m_currTrack.daughters.clear(); // 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()); @@ -217,9 +218,10 @@ void Geant4ParticleHandler::end(const G4Track* track) { // // 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(); + const G4ThreeVector& v = h.vertex(); + m_currTrack.vex = v.x(); + m_currTrack.vey = v.y(); + m_currTrack.vez = v.z(); m_currTrack.pex = m.x(); m_currTrack.pey = m.y(); m_currTrack.pez = m.z(); @@ -286,7 +288,10 @@ void Geant4ParticleHandler::endEvent(const G4Event* ) { checkConsistency(); } +/// 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! @@ -306,12 +311,16 @@ int Geant4ParticleHandler::recombineParents() { int parent_id = par->g4Parent; bool remove_me = false; + if ( id == break_trackID ) { // Used for debugging to set break point + remove_me = false; + } + /// Primary particles MUST be kept! if ( mask.isSet(G4PARTICLE_PRIMARY) ) { continue; } else if ( keep_parent ) { - continue; + //continue; } else if ( keep_process ) { ParticleMap::iterator ip = m_particleMap.find(parent_id); @@ -364,7 +373,6 @@ int Geant4ParticleHandler::recombineParents() { m_particleMap.erase(ir); } } - printout(INFO,name(),"+++ Size of track container:%d",int(m_particleMap.size())); return int(remove.size()); } diff --git a/DDG4/src/Geant4ParticlePrint.cpp b/DDG4/src/Geant4ParticlePrint.cpp index 152e08bec614977214c15b939da0f1a406c58cae..8a7768f078370a527fcba5229efcaaddc8f73d3c 100644 --- a/DDG4/src/Geant4ParticlePrint.cpp +++ b/DDG4/src/Geant4ParticlePrint.cpp @@ -43,11 +43,13 @@ 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); + if ( truth ) { + const ParticleMap& particles = truth->particles(); + if ( (m_outputType&1) != 0 ) printParticles(particles); // Table printout.... + if ( (m_outputType&2) != 0 ) printParticleTree(particles); // Tree printout.... + return; + } + except("No Geant4MonteCarloTruth object present in the event structure to access the particle information!", c_name()); } void Geant4ParticlePrint::printParticle(const std::string& prefix, const Particle* p) const { @@ -63,7 +65,7 @@ void Geant4ParticlePrint::printParticle(const std::string& prefix, const Particl 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]", + 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(), @@ -94,7 +96,7 @@ void Geant4ParticlePrint::printParticles(const ParticleMap& particles) const { int num_tracker_hits = 0; printout(INFO,name(),"+++ MC Particles #Tracks:%6d %-12s Parent%-7s " - "Primary Secondary Energy %-9s Calo Tracker Process/Par Details", + "Primary Secondary Energy %-8s 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; @@ -109,11 +111,11 @@ void Geant4ParticlePrint::printParticles(const ParticleMap& particles) const { 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","", + "Primary Secondary Energy %-8s Calo Tracker Process/Parent", + int(particles.size()),"","",""); + printout(INFO,name(),"+++ MC Particle Summary: %6s %-12s %6s%-7s " + "%7d %9d %6d %-8s %4d %7d %7d %6d", + "","","","", num_primary, num_secondaries, num_energy,"", num_calo_hits,num_tracker_hits,num_process,num_parent); } @@ -144,6 +146,9 @@ void Geant4ParticlePrint::printParticleTree(const ParticleMap& particles, int le /// 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())); + printout(INFO,name(),"+++ MC Particles %12s #Tracks:%7d %-12s Parent%-7s " + "Primary Secondary Energy %-8s 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); diff --git a/DDG4/src/Geant4SensDetAction.cpp b/DDG4/src/Geant4SensDetAction.cpp index b8136abbd2caae74cdd334e2a7325fda83547d22..7f190989deef8257aa99962f4acc9d3e577d54cf 100644 --- a/DDG4/src/Geant4SensDetAction.cpp +++ b/DDG4/src/Geant4SensDetAction.cpp @@ -176,7 +176,7 @@ void Geant4Sensitive::clear(G4HCofThisEvent* /* HCE */) { /// Mark the track to be kept for MC truth propagation during hit processing void Geant4Sensitive::mark(const G4Track* track) const { Geant4MonteCarloTruth* truth = context()->event().extension<Geant4MonteCarloTruth>(false); - if ( truth ) truth->mark(track,true); + if ( truth ) truth->mark(track); } /// Mark the track of this step to be kept for MC truth propagation during hit processing diff --git a/DDG4/src/Geant4TrackPersistency.cpp b/DDG4/src/Geant4TrackPersistency.cpp deleted file mode 100644 index 272542523814d5608c4d62f2a652520d77b63ae8..0000000000000000000000000000000000000000 --- a/DDG4/src/Geant4TrackPersistency.cpp +++ /dev/null @@ -1,549 +0,0 @@ -// $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/Geant4TrackPersistency.h" -#include "DDG4/Geant4EventAction.h" -#include "DDG4/Geant4TrackingAction.h" -#include "DDG4/Geant4SteppingAction.h" -#include "DDG4/Geant4MonteCarloRecordManager.h" -#include "DDG4/Geant4TrackHandler.h" -#include "DDG4/Geant4StepHandler.h" - -#include "G4TrackingManager.hh" -#include "G4ParticleDefinition.hh" -#include "G4Track.hh" -#include "G4Step.hh" -#include "G4TrackStatus.hh" -#include <set> -#include <stdexcept> - -using namespace std; -using DD4hep::InstanceCount; -using namespace DD4hep; -using namespace DD4hep::Simulation; - - - -#include <map> -namespace DD4hep { - namespace Simulation { - - struct Track; - struct Vertex; - typedef std::vector<int> TrackIDs; - typedef std::vector<Track*> Tracks; - typedef std::vector<Vertex*> Vertices; - typedef Geant4TrackPersistency::FourMomentum FourMomentum; - - struct Vertex { - G4ThreeVector position; - Tracks outgoing; - Track* incoming; - int refCount; - Vertex() : incoming(0), refCount(0) { - InstanceCount::increment(this); - } - virtual ~Vertex() { - InstanceCount::decrement(this); - } - void addRef() { - ++refCount; - } - int release() { - int tmp = --refCount; - if ( tmp < 1 ) { - delete this; - } - return tmp; - } - }; - - typedef std::map<int,void*> TrackMap; - typedef std::map<const void*,Track*> TrackMap2; - - struct Track { - enum State { - KEEP=0x1, - CREATED_HIT=0x2, - VALID_G4TRACK=0x4, - TRACKED_G4TRACK=0x8 - }; - - TrackMap2 secondaries; - TrackIDs contained; - FourMomentum momentum; - Vertices daughters; - Vertex* vertex; - Track* parent; - int id; - int refCount; - int flag; - const G4VProcess* process; - const G4ParticleDefinition* def; - - /// Standard constructor - Track() : vertex(0), parent(0), id(0), refCount(0), flag(0), process(0), def(0) - { - InstanceCount::increment(this); - } - virtual ~Track() - { - InstanceCount::decrement(this); - } - bool toBeKept() const { return (flag&KEEP) == KEEP; } - bool createdHit() const { return (flag&CREATED_HIT) == CREATED_HIT; } - bool validG4Track() const { return (flag&VALID_G4TRACK) == VALID_G4TRACK; } - bool wasTracked() const { return (flag&TRACKED_G4TRACK) == TRACKED_G4TRACK; } - }; - } -} - - -void collectContained(const Track* trk, std::map<int,int>& equiv) { - if ( !trk->contained.empty() ) { - for(TrackIDs::const_iterator i=trk->contained.begin(); i!=trk->contained.end(); ++i) - equiv[*i] = trk->id; - } - for(TrackMap2::const_iterator j=trk->secondaries.begin(); j!=trk->secondaries.end(); ++j) - collectContained((*j).second,equiv); -} - -void printTrack(int flg, const char* msg); -void printTree(const Track* trk,const char* msg, int cnt); -void printProperties(const Track* trk); -void printSecondaries(const Track* trk); -void printTrack(const Track* trk, int flg=0, const char* msg="Interesting"); -/// Release track object -void releaseTrack(Geant4TrackPersistency* p, Track* trk); -bool cleanTree(Geant4TrackPersistency* p, const string& proc, Track* trk); - -bool keepTrackTree(const string& proc, Track* trk) { - // We have to check if any of the daughters has the KEEP flag ON. - // If YES, this track may not be deleted, otherwise, we can release the track - string pnam = trk->process ? trk->process->GetProcessName().c_str() : "---"; - if ( pnam == proc ) trk->flag &= ~(Track::KEEP); - if ( (trk->flag&Track::KEEP) == Track::KEEP ) return true; - for(TrackMap2::const_iterator j=trk->secondaries.begin(); j!=trk->secondaries.end(); ++j) - if ( keepTrackTree(proc,(*j).second) ) return true; - return false; -} - -bool cleanTree(Geant4TrackPersistency* p, const string& proc, Track* trk) { - bool remove = false; - for(TrackMap2::iterator j=trk->secondaries.begin(); j!=trk->secondaries.end(); ++j) - remove |= cleanTree(p,proc,(*j).second); - - if ( remove ) { - string pnam = trk->process ? trk->process->GetProcessName().c_str() : "---"; - //if ( pnam == proc ) { - bool keep = keepTrackTree(proc,trk); - if ( !keep ) { - releaseTrack(p,trk); - return true; - } - //} - } - return false; -} - -void printTree(const Track* trk,const char* msg, int cnt) { - char text[64]; - ::snprintf(text,sizeof(text)," %-8d | ",cnt); - printf(" %-6d ",cnt); - ::snprintf(text,sizeof(text),"| "); - std::string m(text); - m += msg; - printTrack(trk, 1, msg); - int c = 0; - for(TrackMap2::const_iterator j=trk->secondaries.begin(); j!=trk->secondaries.end(); ++j) - printTree((*j).second,m.c_str(),++c); -} - -void printProperties(const Track* trk) { - const char* pnam = trk->process ? trk->process->GetProcessName().c_str() : "---"; - int pid = trk->parent ? trk->parent->id : 0; - int sec = int(trk->secondaries.size()); - int vtx = 0; - for(Vertices::const_iterator i=trk->daughters.begin(); i!=trk->daughters.end();++i) - vtx += (*i)->outgoing.size(); - if ( vtx != sec ) { - ::printf(" ++++ ERROR! ++++ "); - } - ::printf("Track %4d Parent:%4d %-8s %-6s Sec:%3d [%3d] Cont:%3d Keep:%s Tracked:%s Hit:%s", - trk->id, pid, pnam, trk->def->GetParticleName().c_str(), sec, vtx, - int(trk->contained.size()), - yes_no(trk->toBeKept()), yes_no(trk->wasTracked()), - yes_no(trk->createdHit())); -} - -void printContained(const Track* trk) { - const TrackIDs& ids = trk->contained; - if ( !ids.empty() ) { - ::printf(" Equiv(%d):",int(ids.size())); - for(TrackIDs::const_iterator i=ids.begin(); i!=ids.end(); ++i) ::printf(" %d",*i); - } -} - -void printSecondaries(const Track* trk) { - char text[256]; - int ndau = 1; - for(TrackMap2::const_iterator j=trk->secondaries.begin(); j!=trk->secondaries.end(); ++j, ++ndau) { - ::snprintf(text,sizeof(text)," ---> Daughter [%3d]",ndau); - ::printTrack((*j).second,2,text); - } -} -void printParents(const Track* trk) { - char text[256]; - int npar = 1; - for(Track* p=trk->parent; p; p=p->parent, ++npar) { - ::snprintf(text,sizeof(text)," ---> Parent [%3d]",npar); - printTrack(p,4,text); - } -} - -void printTrack(const Track* trk, int flg, const char* msg) { - printf("%s> ",msg); - printProperties(trk); - if ( (flg&1)==1 ) printContained(trk); - printf("\n"); - if ( (flg&2)==2 ) printSecondaries(trk); - if ( (flg&4)==4 ) printParents(trk); -} - -void referenceTrack(Track* trk) { - if ( trk ) { - ++trk->refCount; - referenceTrack(trk->parent); - } -} - -/// Release track object -void releaseTrack(Geant4TrackPersistency* p, Track* trk) { - int cnt = --trk->refCount; - Track* par = trk->parent; - if ( cnt < 1 ) { - int id = trk->id; - Vertex* v = trk->vertex; - if ( v ) { - if ( par ) { - Tracks& trks = v->outgoing; - for(Tracks::iterator it=trks.begin(); it != trks.end(); ++it) { - if ( *it == trk ) { trks.erase(it); break; } - } - if ( trks.size() == 0 ) { - Vertices::iterator iv=find(par->daughters.begin(),par->daughters.end(),v); - if ( iv != par->daughters.end() ) { - par->daughters.erase(iv); - } - } - for(TrackMap2::iterator is=par->secondaries.begin(); is!=par->secondaries.end(); ++is) { - if ( (*is).second == trk ) { - par->secondaries.erase(is); - break; - } - } - par->contained.push_back(id); - } - v->release(); - } - std::for_each(trk->daughters.begin(),trk->daughters.end(),releasePtr<Vertex>); - delete trk; - - TrackMap::iterator m = p->m_trackMap.find(id); - if ( m != p->m_trackMap.end() ) p->m_trackMap.erase(m); - //printf(" ---> DT: REMOVE track %p [%s] %d size:%d\n", (void*)trk, pn, id, int(p->m_trackMap.size())); - } - if ( par ) { - releaseTrack(p, par); - } -} - -/// Initialize track data -Track* initTrack(const G4Track* track, Track* trk, Track* par) { - G4ThreeVector mom = track->GetMomentum(); - trk->momentum.SetPxPyPzE(mom.x(),mom.y(),mom.z(),track->GetTotalEnergy()); - trk->process = track->GetCreatorProcess(); - trk->def = track->GetParticleDefinition(); - trk->id = track->GetTrackID(); - trk->parent = par; - if ( par ) { - par->secondaries[track] = trk; - } - referenceTrack(trk); - return trk; -} - -void deleteTrackTree(Track* trk) { - Vertex* v = trk->vertex; - if ( v ) { - //if ( v->refCount == 2 ) --v->refCount; - v->release(); - } - for(TrackMap2::iterator j=trk->secondaries.begin(); j!=trk->secondaries.end(); ++j) - deleteTrackTree((*j).second); - delete trk; -} - -/// Standard constructor -Geant4TrackPersistency::TrackInfo::TrackInfo() : G4VUserTrackInformation() { - this->info = 0; - this->manager = 0; - this->track = 0; - this->store = false; - this->createdHit = false; - this->initialEnergy = 0; - InstanceCount::increment(this); -} - -/// Default destructor -Geant4TrackPersistency::TrackInfo::~TrackInfo() { - InstanceCount::decrement(this); -} - -/// Clear alkl stored information for this track -void Geant4TrackPersistency::TrackInfo::set(const G4Track* trk, void* inf) { - this->track = trk; - if ( trk ) { - G4ThreeVector mom = trk->GetMomentum(); - this->initialMomentum.SetPxPyPzE(mom.x(),mom.y(),mom.z(),trk->GetTotalEnergy()); - this->initialEnergy = trk->GetTotalEnergy(); - this->info = inf; - } - else { - this->initialMomentum.SetPxPyPzE(0e0,0e0,0e0,0e0); - this->initialEnergy = 0; - this->createdHit = false; - this->store = false; - this->info = 0; - } -} - -/// Standard constructor -Geant4TrackPersistency::Geant4TrackPersistency(Geant4Context* context, const std::string& nam) - : Geant4GeneratorAction(context,nam), Geant4MonteCarloTruth(), m_current() -{ - m_current.manager = this; - generatorAction().adopt(this); - eventAction().callAtBegin(this,&Geant4TrackPersistency::beginEvent); - eventAction().callAtFinal(this,&Geant4TrackPersistency::endEvent); - trackingAction().callAtFinal(this,&Geant4TrackPersistency::end,CallbackSequence::FRONT); - trackingAction().callUpFront(this,&Geant4TrackPersistency::begin,CallbackSequence::FRONT); - steppingAction().call(this,&Geant4TrackPersistency::step); - - declareProperty("PrintEndTracking",m_printEndTracking = false); - declareProperty("PrintStartTracking",m_printStartTracking = false); - - InstanceCount::increment(this); -} - -/// Default destructor -Geant4TrackPersistency::~Geant4TrackPersistency() { - InstanceCount::decrement(this); -} - -/// Mark a Geant4 track to be kept for later MC truth analysis -void Geant4TrackPersistency::mark(const G4Track* track, int reason) { - mark(track); // Does all the checking... - if ( reason ) { - Track* trk = (Track*)m_current.info; - trk->flag |= Track::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, int reason) { - if ( step ) { - mark(step->GetTrack(),reason); - return; - } - except("Cannot mark the G4Track if the step-pointer is invalid!", c_name()); -} - -/// Mark a Geant4 track of the step to be kept for later MC truth analysis -void Geant4TrackPersistency::mark(const G4Step* step) { - if ( step ) { - mark(step->GetTrack()); - return; - } - except("Cannot mark the G4Track if the step-pointer is invalid!", c_name()); -} - -/// Mark a Geant4 track of the step to be kept for later MC truth analysis -void Geant4TrackPersistency::mark(const G4Track* track) { - if ( 0 == m_current.track ) { - except("Miconfigured program. The tracking preaction was never called for the Geant4TrackPersistency"); - } - else if ( track == m_current.track ) { - m_current.store = true; - Track* trk = (Track*)m_current.info; - trk->flag |= Track::KEEP; - referenceTrack(trk); - return; - } - except("Call to Geant4TrackPersistency with inconsistent G4Track pointer: %p <> %p", - track,m_current.track); -} - -/// Event generation action callback -void Geant4TrackPersistency::operator()(G4Event* ) { - typedef Geant4MonteCarloTruth _MC; - m_trackMap.clear(); - printout(INFO,name(),"+++ Add EVENT extension of type Geant4TrackPersistency....."); - context()->event().addExtension((_MC*)this, typeid(_MC), 0); -} - -/// Pre-event action callback -void Geant4TrackPersistency::beginEvent(const G4Event* ) { -} - -/// Pre-event action callback -void Geant4TrackPersistency::endEvent(const G4Event* ) { - int cnt = 1; - printf("---------------------------- End of Event: %d tracks kept ----------------------------\n", - int(m_trackMap.size())); - for (TrackMap::iterator m=m_trackMap.begin(); m != m_trackMap.end(); ++m, ++cnt) { - Track* trk = (Track*)(*m).second; - if ( trk->parent == 0 ) { - printTree(trk,"+--",1); - cleanTree(this,"eBrem",trk); - break; - } - } - - printf("---------------------------- End of Event: %d tracks kept ----------------------------\n", - int(m_trackMap.size())); - for(TrackMap::iterator m=m_trackMap.begin(); m != m_trackMap.end(); ++m, ++cnt) { - Track* trk = (Track*)(*m).second; - if ( trk->parent == 0 ) { - printTree(trk,"+--",1); - cleanTree(this,"eIoni",trk); - break; - } - } - - printf("---------------------------- End of Event: %d tracks kept ----------------------------\n", - int(m_trackMap.size())); - for(TrackMap::iterator m=m_trackMap.begin(); m != m_trackMap.end(); ++m, ++cnt) { - Track* trk = (Track*)(*m).second; - if ( trk->parent == 0 ) { - printTree(trk,"+--",1); - deleteTrackTree(trk); - break; - } - - //char text[64]; - //::sprintf(text,"Intersting [%3d]",cnt); - //trk->print(7,text); - } - m_trackMap.clear(); -} - -/// User stepping callback -void Geant4TrackPersistency::step(const G4Step* step, G4SteppingManager* /* mgr */) { - typedef std::vector<const G4Track*> _Sec; - const _Sec* sec=step->GetSecondaryInCurrentStep(); - if ( sec->size() > 0 ) { - Track* t = (Track*)m_current.info; - Vertex* v = new Vertex(); - v->incoming = t; - t->daughters.push_back(v); - - for(_Sec::const_iterator i=sec->begin(); i!=sec->end(); ++i) { - Track* secondary = initTrack(*i,new Track(),t); - secondary->vertex = v; - secondary->vertex->addRef(); - v->outgoing.push_back(secondary); - } - } -} - -/// Pre-track action callback -void Geant4TrackPersistency::begin(const G4Track* track) { - G4VUserTrackInformation* info = &m_current; - int pid = track->GetParentID(); - Track* trk = 0; - - if ( pid == 0 ) { - trk = (Track*)initTrack(track,new Track(),0); - trk->vertex = new Vertex(); - trk->vertex->incoming = 0; - trk->vertex->position = track->GetPosition(); - m_trackMap.insert(make_pair(trk->id,trk)); - } - else { - Track* par = (Track*)m_trackMap[pid]; - trk = par->secondaries[track]; - trk->id = track->GetTrackID(); - trk->vertex->incoming = par; - trk->vertex->position = track->GetPosition(); - m_trackMap.insert(make_pair(trk->id,trk)); - } - - if ( m_printStartTracking ) { - Geant4TrackHandler t(track); - std::string proc = t.creatorName(); - G4ThreeVector pos = t.position(); - - this->info("START Track %p ID=%4d Parent:%4d E:%7.2f MeV " - "%-6s %-8s L:%7.2f Pos:(%8.2f %8.2f %8.2f) Ref:%d",track, - t.id(), t.parent(), t.energy(), - t.name().c_str(), ("["+proc+"]").c_str(), t.length(), - pos.x(),pos.y(),pos.z(),trk->refCount); - } - - trk->flag |= Track::VALID_G4TRACK; - m_current.set(track,trk); - trackMgr()->SetUserTrackInformation(info); - - Geant4TrackHandler t(track); - std::string proc = t.creatorName(); - if ( proc == "conv" ) { - mark(track); - } - if ( proc == "eBrem" ) { - mark(track); - } - else if ( trk->id < 5 ) - mark(track); -} - -/// Post-track action callback -void Geant4TrackPersistency::end(const G4Track* track) { - Track* trk = (Track*)m_current.info; - - trk->flag |= Track::TRACKED_G4TRACK; - trk->flag &= ~(Track::VALID_G4TRACK); - if ( m_printEndTracking ) { - Geant4TrackHandler t(track); - std::string proc = t.creatorName(); - G4ThreeVector pos = t.position(); - - this->info("END Track %p ID=%4d Parent:%4d E:%7.2f MeV " - "%-6s %-8s L:%7.2f Pos:(%8.2f %8.2f %8.2f) Ref:%d Keep:%s Tracked:%s Hit:%s", - track,t.id(), t.parent(), t.energy(), - t.name().c_str(), ("["+proc+"]").c_str(), t.length(), - pos.x(),pos.y(),pos.z(),trk->refCount, - yes_no(trk->toBeKept()), yes_no(trk->wasTracked()), - yes_no(trk->createdHit())); - } - /// If required save Track record... - if ( m_current.store ) { - Geant4MonteCarloRecordManager* mgr = - context()->event().extension<Geant4MonteCarloRecordManager>(false); - if ( mgr ) mgr->save(m_current); - } - m_current.set(0,0); - trackMgr()->SetUserTrackInformation(0); - // All done. Check if the track can be released/deleted - releaseTrack(this,trk); -}