DDG4/src/Geant4ParticlePrint.cpp

0001 //==========================================================================
0002 //  AIDA Detector description implementation
0003 //--------------------------------------------------------------------------
0004 // Copyright (C) Organisation europeenne pour la Recherche nucleaire (CERN)
0005 // All rights reserved.
0006 //
0007 // For the licensing terms see $DD4hepINSTALL/LICENSE.
0008 // For the list of contributors see $DD4hepINSTALL/doc/CREDITS.
0009 //
0010 // Author     : M.Frank
0011 //
0012 //==========================================================================
0013
0014 // Framework include files
0015 #include <DD4hep/Printout.h>
0016 #include <DD4hep/Primitives.h>
0017 #include <DD4hep/InstanceCount.h>
0018 #include <DDG4/Geant4ParticlePrint.h>
0019 #include <DDG4/Geant4Data.h>
0020 #include <DDG4/Geant4HitCollection.h>
0021
0022 // Geant4 include files
0023 #include <G4Event.hh>
0024
0025 // C/C++ include files
0026 #include <cstring>
0027
0028 using namespace dd4hep::sim;
0029 using PropertyMask = dd4hep::detail::ReferenceBitMask<const int>;
0030
0031 /// Standard constructor
0032 Geant4ParticlePrint::Geant4ParticlePrint(Geant4Context* ctxt, const std::string& nam)
0033   : Geant4EventAction(ctxt,nam)
0034 {
0035   declareProperty("OutputType",m_outputType=3);
0036   declareProperty("PrintBegin",m_printBegin=false);
0037   declareProperty("PrintEnd",  m_printEnd=true);
0038   declareProperty("PrintGen",  m_printGeneration=true);
0039   declareProperty("PrintHits", m_printHits=false);
0040   InstanceCount::increment(this);
0041 }
0042
0043 /// Default destructor
0044 Geant4ParticlePrint::~Geant4ParticlePrint()  {
0045   InstanceCount::decrement(this);
0046 }
0047
0048 /// Print particle table
0049 void Geant4ParticlePrint::makePrintout(const G4Event* e) const  {
0050   Geant4ParticleMap* parts = context()->event().extension<Geant4ParticleMap>();
0051   if ( parts )   {
0052     const ParticleMap& particles = parts->particles();
0053     print("+++ ******** MC Particle Printout for event ID:%d ********",e->GetEventID());
0054     if ( (m_outputType&2) != 0 ) printParticleTree(e,particles);  // Tree printout....
0055     if ( (m_outputType&1) != 0 ) printParticles(0,particles);     // Table printout....
0056     return;
0057   }
0058   except("No Geant4MonteCarloTruth object present in the event structure to access the particle information!", c_name());
0059 }
0060
0061 /// Generation action callback
0062 void Geant4ParticlePrint::operator()(G4Event* e)  {
0063   if ( m_printGeneration ) makePrintout(e);
0064 }
0065
0066 /// Pre-event action callback
0067 void Geant4ParticlePrint::begin(const G4Event* e)  {
0068   if ( m_printBegin ) makePrintout(e);
0069 }
0070
0071 /// Post-event action callback
0072 void Geant4ParticlePrint::end(const G4Event* e)  {
0073   if ( m_printEnd ) makePrintout(e);
0074 }
0075
0076 void Geant4ParticlePrint::printParticle(const std::string& prefix, const G4Event* e, Geant4ParticleHandle p) const   {
0077   char equiv[32];
0078   PropertyMask mask(p->reason);
0079   PropertyMask status(p->status);
0080   std::string proc_name = p.processName();
0081   std::string proc_type = p.processTypeName();
0082   int parent_id = p->parents.empty() ? -1 : *(p->parents.begin());
0083
0084   equiv[0] = 0;
0085   if ( p->parents.end() == p->parents.find(p->g4Parent) )  {
0086     ::snprintf(equiv,sizeof(equiv),"/%d",p->g4Parent);
0087   }
0088   print("+++ %s ID:%7d %12s %6d%-7s %7s %3s %5d %3s %+.3e  %-4s %-7s %-3s %-3s %2d  [%s%s%s] %c%c%c%c",
0089         prefix.c_str(),
0090         p->id,
0091         p.particleName().c_str(),
0092         parent_id,equiv,
0093         yes_no(mask.isSet(G4PARTICLE_PRIMARY)),
0094         yes_no(mask.isSet(G4PARTICLE_HAS_SECONDARIES)),
0095         int(p->daughters.size()),
0096         yes_no(mask.isSet(G4PARTICLE_ABOVE_ENERGY_THRESHOLD)),
0097         p.energy(),
0098         yes_no(mask.isSet(G4PARTICLE_CREATED_CALORIMETER_HIT)),
0099         yes_no(mask.isSet(G4PARTICLE_CREATED_TRACKER_HIT)),
0100         yes_no(mask.isSet(G4PARTICLE_KEEP_PROCESS)),
0101         mask.isSet(G4PARTICLE_KEEP_PARENT) ? "YES" : "",
0102         p.numParent(),
0103         proc_name.c_str(),
0104         p->process ? "/" : "",
0105         proc_type.c_str(),
0106         status.isSet(G4PARTICLE_GEN_EMPTY) ? 'E' : '.',
0107         status.isSet(G4PARTICLE_GEN_STABLE) ? 'S' : '.',
0108         status.isSet(G4PARTICLE_GEN_DECAYED) ? 'D' : '.',
0109         status.isSet(G4PARTICLE_GEN_DOCUMENTATION) ? 'd' : '.'
0110         );
0111   if ( e && m_printHits )  {
0112     Geant4ParticleMap* truth = context()->event().extension<Geant4ParticleMap>();
0113     G4HCofThisEvent* hc = e->GetHCofThisEvent();
0114     for (int ihc=0, nhc=hc->GetNumberOfCollections(); ihc<nhc; ++ihc)   {
0115       G4VHitsCollection* c = hc->GetHC(ihc);
0116       Geant4HitCollection* coll = dynamic_cast<Geant4HitCollection*>(c);
0117       if ( coll )  {
0118         size_t nhits = coll->GetSize();
0119         for(size_t i=0; i<nhits; ++i)   {
0120           Geant4HitData* h = coll->hit(i);
0121           Geant4Tracker::Hit* trk_hit = dynamic_cast<Geant4Tracker::Hit*>(h);
0122           if ( 0 != trk_hit )   {
0123             Geant4HitData::Contribution& t = trk_hit->truth;
0124             int trackID = t.trackID;
0125             int trueID  = truth->particleID(trackID);
0126             if ( trueID == p->id )   {
0127               print("+++ %20s           %s[%d]  (%+.2e,%+.2e,%+.2e)[mm]","",c->GetName().c_str(),i,
0128                     trk_hit->position.x(),trk_hit->position.y(),trk_hit->position.z());
0129             }
0130           }
0131           Geant4Calorimeter::Hit* cal_hit = dynamic_cast<Geant4Calorimeter::Hit*>(h);
0132           if ( 0 != cal_hit )   {
0133             Geant4HitData::Contributions& contrib = cal_hit->truth;
0134             for(Geant4HitData::Contributions::iterator j=contrib.begin(); j!=contrib.end(); ++j)  {
0135               Geant4HitData::Contribution& t = *j;
0136               int trackID = t.trackID;
0137               int trueID  = truth->particleID(trackID);
0138               if ( trueID == p->id )   {
0139                 print("+++ %20s           %s[%d]  (%+.2e,%+.2e,%+.2e)[mm]","",c->GetName().c_str(),i,
0140                       cal_hit->position.x(),cal_hit->position.y(),cal_hit->position.z());
0141               }
0142             }
0143           }
0144         }
0145       }
0146       else  {
0147         print("+++ Hit unknown hit collection type: %s --> %s",
0148               c->GetName().c_str(),typeName(typeid(*c)).c_str());
0149       }
0150     }
0151   }
0152 }
0153
0154 /// Print record of kept particles
0155 void Geant4ParticlePrint::printParticles(const G4Event* e, const ParticleMap& particles) const  {
0156   int num_energy = 0;
0157   int num_parent = 0;
0158   int num_process = 0;
0159   int num_primary = 0;
0160   int num_secondaries = 0;
0161   int num_calo_hits = 0;
0162   int num_tracker_hits = 0;
0163
0164   print("+++ MC Particles #Tracks:%7d ParticleType Parent/Geant4 "
0165         "Primary Secondary Energy in [MeV] Calo Tracker Process/Par Details",
0166         int(particles.size()));
0167   for(ParticleMap::const_iterator i=particles.begin(); i!=particles.end(); ++i)  {
0168     Geant4ParticleHandle p = (*i).second;
0169     PropertyMask mask(p->reason);
0170     printParticle("MC Particle Track",e, p);
0171     num_secondaries += int(p->daughters.size());
0172     if ( mask.isSet(G4PARTICLE_ABOVE_ENERGY_THRESHOLD) ) ++num_energy;
0173     if ( mask.isSet(G4PARTICLE_PRIMARY) ) ++num_primary;
0174     if ( mask.isSet(G4PARTICLE_CREATED_TRACKER_HIT) ) ++num_tracker_hits;
0175     if ( mask.isSet(G4PARTICLE_CREATED_CALORIMETER_HIT) ) ++num_calo_hits;
0176     if ( mask.isSet(G4PARTICLE_KEEP_PARENT) ) ++num_parent;
0177     else if ( mask.isSet(G4PARTICLE_KEEP_PROCESS) ) ++num_process;
0178   }
0179   print("+++ MC Particles #Tracks:%7d ParticleType Parent/Geant4 "
0180         "Primary Secondary Energy          Calo Tracker Process/Par",
0181         int(particles.size()));
0182   print("+++ MC Particle Summary:                       %7d %10d %7d  %7d      %9d %5d %6d",
0183         num_primary, num_secondaries, num_energy,
0184         num_calo_hits,num_tracker_hits,num_process,num_parent);
0185 }
0186
0187 void Geant4ParticlePrint::printParticleTree(const G4Event* e,
0188                                             const ParticleMap& particles,
0189                                             int level,
0190                                             Geant4ParticleHandle p)  const
0191 {
0192   char txt[64];
0193   size_t len = sizeof(txt)-33; // Careful about overruns...
0194   // Ensure we do not overwrite the array
0195   if ( level>int(len)-3 ) level = len-3;
0196
0197   ::snprintf(txt,sizeof(txt),"%5d ",level);
0198   ::memset(txt+6,' ',len-6);
0199   txt[len-1] = 0;
0200   txt[len-2] = '>';
0201   if ( size_t(level + 6) < len )    {
0202     txt[level+6]='+';
0203     ::memset(txt+level+6+1,'-',len-level-3-6);
0204   }
0205
0206   printParticle(txt, e, p);
0207   // For all particles, the set of daughters must be contained in the record.
0208   for( int id_dau : p->daughters )  {
0209     auto iter = particles.find(id_dau);
0210     if ( iter != particles.end() )  {
0211       Geant4ParticleHandle dau = iter->second;
0212       printParticleTree(e, particles, level+1, dau);
0213     }
0214   }
0215 }
0216
0217 /// Print tree of kept particles
0218 void Geant4ParticlePrint::printParticleTree(const G4Event* e, const ParticleMap& particles)  const  {
0219   print("+++ MC Particle Parent daughter relationships. [%d particles]",int(particles.size()));
0220   print("+++ MC Particles %12s #Tracks:%7d %-12s Parent%-7s "
0221         "Primary Secondary Energy %-8s Calo Tracker Process/Par  Details",
0222         "",int(particles.size()),"ParticleType","","in [MeV]");
0223   for(ParticleMap::const_iterator i=particles.begin(); i!=particles.end(); ++i)  {
0224     Geant4ParticleHandle p = (*i).second;
0225     PropertyMask mask(p->reason);
0226     if ( mask.isSet(G4PARTICLE_PRIMARY) ) printParticleTree(e, particles, 0, p);
0227   }
0228 }