// $Id: $
//==========================================================================
// AIDA Detector description implementation for LCD
//--------------------------------------------------------------------------
// Copyright (C) Organisation europeenne pour la Recherche nucleaire (CERN)
// All rights reserved.
//
// For the licensing terms see $DD4hepINSTALL/LICENSE.
// For the list of contributors see $DD4hepINSTALL/doc/CREDITS.
//
// 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 "DDG4/Geant4HitCollection.h"
// Geant4 include files
#include "G4Event.hh"
// C/C++ include files
#include <cstring>
using namespace std;
using namespace DD4hep;
using namespace DD4hep::Simulation;
typedef ReferenceBitMask<const int> PropertyMask;
/// Standard constructor
Geant4ParticlePrint::Geant4ParticlePrint(Geant4Context* ctxt, const std::string& nam)
: Geant4EventAction(ctxt,nam)
{
declareProperty("OutputType",m_outputType=3);
declareProperty("PrintBegin",m_printBegin=false);
declareProperty("PrintEnd", m_printEnd=true);
declareProperty("PrintGen", m_printGeneration=true);
declareProperty("PrintHits", m_printHits=false);
InstanceCount::increment(this);
}
/// Default destructor
Geant4ParticlePrint::~Geant4ParticlePrint() {
InstanceCount::decrement(this);
}
/// Print particle table
void Geant4ParticlePrint::makePrintout(const G4Event* e) const {
Geant4ParticleMap* parts = context()->event().extension<Geant4ParticleMap>();
if ( parts ) {
const ParticleMap& particles = parts->particles();
print("+++ ******** MC Particle Printout for event ID:%d ********",e->GetEventID());
if ( (m_outputType&2) != 0 ) printParticleTree(e,particles); // Tree printout....
if ( (m_outputType&1) != 0 ) printParticles(0,particles); // Table printout....
return;
}
except("No Geant4MonteCarloTruth object present in the event structure to access the particle information!", c_name());
}
/// Generation action callback
void Geant4ParticlePrint::operator()(G4Event* e) {
if ( m_printGeneration ) makePrintout(e);
}
/// Pre-event action callback
void Geant4ParticlePrint::begin(const G4Event* e) {
if ( m_printBegin ) makePrintout(e);
}
/// Post-event action callback
void Geant4ParticlePrint::end(const G4Event* e) {
if ( m_printEnd ) makePrintout(e);
}
void Geant4ParticlePrint::printParticle(const std::string& prefix, const G4Event* e, Geant4ParticleHandle p) const {
char equiv[32];
PropertyMask mask(p->reason);
PropertyMask status(p->status);
string proc_name = p.processName();
string proc_type = p.processTypeName();
int parent_id = p->parents.empty() ? -1 : *(p->parents.begin());
equiv[0] = 0;
if ( p->parents.end() == p->parents.find(p->g4Parent) ) {
::snprintf(equiv,sizeof(equiv),"/%d",p->g4Parent);
}
print("+++ %s ID:%7d %12s %6d%-7s %7s %3s %5d %3s %+.3e %-4s %-7s %-3s %-3s %2d [%s%s%s] %c%c%c%c",
prefix.c_str(),
p->id,
p.particleName().c_str(),
parent_id,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" : "",
p.numParent(),
proc_name.c_str(),
p->process ? "/" : "",
proc_type.c_str(),
status.isSet(G4PARTICLE_GEN_EMPTY) ? 'E' : '.',
status.isSet(G4PARTICLE_GEN_STABLE) ? 'S' : '.',
status.isSet(G4PARTICLE_GEN_DECAYED) ? 'D' : '.',
status.isSet(G4PARTICLE_GEN_DOCUMENTATION) ? 'd' : '.'
);
if ( e && m_printHits ) {
Geant4ParticleMap* truth = context()->event().extension<Geant4ParticleMap>();
G4HCofThisEvent* hc = e->GetHCofThisEvent();
for (int ihc=0, nhc=hc->GetNumberOfCollections(); ihc<nhc; ++ihc) {
G4VHitsCollection* c = hc->GetHC(ihc);
Geant4HitCollection* coll = dynamic_cast<Geant4HitCollection*>(c);
size_t nhits = coll->GetSize();
for(size_t i=0; i<nhits; ++i) {
Geant4HitData* h = coll->hit(i);
Geant4Tracker::Hit* trk_hit = dynamic_cast<Geant4Tracker::Hit*>(h);
if ( 0 != trk_hit ) {
Geant4HitData::Contribution& t = trk_hit->truth;
int trackID = t.trackID;
int trueID = truth->particleID(trackID);
if ( trueID == p->id ) {
print("+++ %20s %s[%d] (%+.2e,%+.2e,%+.2e)[mm]","",c->GetName().c_str(),i,
trk_hit->position.x(),trk_hit->position.y(),trk_hit->position.z());
}
}
Geant4Calorimeter::Hit* cal_hit = dynamic_cast<Geant4Calorimeter::Hit*>(h);
if ( 0 != cal_hit ) {
Geant4HitData::Contributions& contrib = cal_hit->truth;
for(Geant4HitData::Contributions::iterator j=contrib.begin(); j!=contrib.end(); ++j) {
Geant4HitData::Contribution& t = *j;
int trackID = t.trackID;
int trueID = truth->particleID(trackID);
if ( trueID == p->id ) {
print("+++ %20s %s[%d] (%+.2e,%+.2e,%+.2e)[mm]","",c->GetName().c_str(),i,
cal_hit->position.x(),cal_hit->position.y(),cal_hit->position.z());
}
}
}
}
}
}
}
/// Print record of kept particles
void Geant4ParticlePrint::printParticles(const G4Event* e, 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;
print("+++ MC Particles #Tracks:%7d ParticleType Parent/Geant4 "
"Primary Secondary Energy in [MeV] Calo Tracker Process/Par Details",
int(particles.size()));
for(ParticleMap::const_iterator i=particles.begin(); i!=particles.end(); ++i) {
Geant4ParticleHandle p = (*i).second;
PropertyMask mask(p->reason);
printParticle("MC Particle Track",e, 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;
}
print("+++ MC Particles #Tracks:%7d ParticleType Parent/Geant4 "
"Primary Secondary Energy Calo Tracker Process/Par",
int(particles.size()));
print("+++ MC Particle Summary: %7d %10d %7d %7d %9d %5d %6d",
num_primary, num_secondaries, num_energy,
num_calo_hits,num_tracker_hits,num_process,num_parent);
}
void Geant4ParticlePrint::printParticleTree(const G4Event* e, const ParticleMap& particles, int level, Geant4ParticleHandle 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-1] = 0;
txt[len-2] = '>';
txt[level+6]='+';
::memset(txt+level+6+1,'-',len-level-3-6);
printParticle(txt, e, 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;
Geant4ParticleHandle dau = (*particles.find(id_dau)).second;
printParticleTree(e, particles, level+1, dau);
}
}
/// Print tree of kept particles
void Geant4ParticlePrint::printParticleTree(const G4Event* e, const ParticleMap& particles) const {
print("+++ MC Particle Parent daughter relationships. [%d particles]",int(particles.size()));
print("+++ 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) {
Geant4ParticleHandle p = (*i).second;
PropertyMask mask(p->reason);
if ( mask.isSet(G4PARTICLE_PRIMARY) ) printParticleTree(e, particles, 0, p);
}
}