File indexing completed on 2025-01-18 09:14:23
0001
0002
0003
0004
0005
0006
0007
0008
0009
0010
0011
0012
0013
0014
0015 #include <CLHEP/Units/SystemOfUnits.h>
0016 #include <DD4hep/DD4hepUnits.h>
0017 #include <DD4hep/Printout.h>
0018 #include <DD4hep/Detector.h>
0019 #include <DDG4/Geant4Particle.h>
0020 #include <DDG4/Geant4Data.h>
0021
0022
0023 #include <vector>
0024 #include <cstdio>
0025 #include <memory>
0026 #include <cerrno>
0027
0028
0029 #include <TFile.h>
0030 #include <TTree.h>
0031 #include <TSystem.h>
0032
0033 using namespace std;
0034 using namespace dd4hep;
0035 using namespace dd4hep::detail;
0036
0037 typedef dd4hep::sim::Geant4Tracker Geant4Tracker;
0038 typedef dd4hep::sim::Geant4Calorimeter Geant4Calorimeter;
0039 typedef dd4hep::sim::Geant4Particle Geant4Particle;
0040
0041 namespace {
0042 static const char* line = "+-------------------------------------------------------------+";
0043 static bool have_geometry = false;
0044
0045 int usage() {
0046 printf("\ndumpDDG4 -opt [-opt] \n"
0047 " -compact <compact-geometry> Supply geometry file to check hits with volume manager.\n"
0048 " -input <root-file> File generated with DDG4 \n"
0049 " -event <event-number> Specify event to be dumped. Default: ALL. \n"
0050 "\n\n");
0051 return EINVAL;
0052 }
0053
0054 int printHits(const string& container, vector<Geant4Tracker::Hit*>* hits) {
0055 typedef vector<Geant4Tracker::Hit*> _H;
0056 if ( !hits ) {
0057 ::printf("+ Invalid Hit container '%s'. No printout\n",container.c_str());
0058 }
0059 else if ( hits->empty() ) {
0060 ::printf("+ Invalid Hit container '%s'. No entries. No printout\n",container.c_str());
0061 }
0062 else if ( have_geometry ) {
0063 string det_name = container;
0064 Detector& description = Detector::getInstance();
0065 det_name = det_name.substr(0,det_name.length()-4);
0066 DetElement det(description.detector(det_name));
0067 SensitiveDetector sd(description.sensitiveDetector(det_name));
0068 Segmentation seg = sd.readout().segmentation();
0069 VolumeManager vm = description.volumeManager();
0070 for(_H::const_iterator i=hits->begin(); i!=hits->end(); ++i) {
0071 const Geant4Tracker::Hit* h = *i;
0072 const Position& pos = h->position;
0073 Position pos_cell = seg.position(h->cellID);
0074 PlacedVolume pv = vm.lookupPlacement(h->cellID);
0075 printout(ALWAYS,container,
0076 "+++ Track:%3d PDG:%6d Pos:(%+.2e,%+.2e,%+.2e)[mm] Pixel:(%+.2e,%+.2e,%+.2e)[mm] %s Deposit:%7.3f MeV CellID:%16lX",
0077 h->truth.trackID,h->truth.pdgID,
0078 pos.x()/CLHEP::mm,pos.y()/CLHEP::mm,pos.z()/CLHEP::mm,
0079 pos_cell.x()/dd4hep::mm,pos_cell.y()/dd4hep::mm,pos_cell.z()/dd4hep::mm,
0080 pv.name(),h->truth.deposit/CLHEP::MeV,h->cellID
0081 );
0082 delete h;
0083 }
0084 }
0085 else {
0086 for(_H::const_iterator i=hits->begin(); i!=hits->end(); ++i) {
0087 const Geant4Tracker::Hit* h = *i;
0088 const Position& pos = h->position;
0089 printout(ALWAYS,container,
0090 "+++ Track:%3d PDG:%6d Pos:(%+.2e,%+.2e,%+.2e)[mm] Deposit:%7.3f MeV CellID:%16lX",
0091 h->truth.trackID,h->truth.pdgID,
0092 pos.x()/CLHEP::mm,pos.y()/CLHEP::mm,pos.z()/CLHEP::mm,
0093 h->truth.deposit/CLHEP::MeV,h->cellID
0094 );
0095 delete h;
0096 }
0097 }
0098 return 0;
0099 }
0100 int printHits(const string& container, vector<Geant4Calorimeter::Hit*>* hits) {
0101 typedef vector<Geant4Calorimeter::Hit*> _H;
0102 if ( !hits ) {
0103 ::printf("+ Invalid Hit container '%s'. No printout\n",container.c_str());
0104 }
0105 else if ( hits->empty() ) {
0106 ::printf("+ Invalid Hit container '%s'. No entries. No printout\n",container.c_str());
0107 }
0108 else {
0109 string det_name = container;
0110 Detector& description = Detector::getInstance();
0111 det_name = det_name.substr(0,det_name.length()-4);
0112 DetElement det(description.detector(det_name));
0113 SensitiveDetector sd(description.sensitiveDetector(det_name));
0114 Segmentation seg = sd.readout().segmentation();
0115 VolumeManager vm = description.volumeManager();
0116 for(_H::const_iterator i=hits->begin(); i!=hits->end(); ++i) {
0117 const Geant4Calorimeter::Hit* h = *i;
0118 const Position& pos = h->position;
0119 Position pos_cell = seg.position(h->cellID);
0120 PlacedVolume pv = vm.lookupPlacement(h->cellID);
0121 printout(ALWAYS,container,
0122 "+++ Pos:(%+.2e,%+.2e,%+.2e)[mm] Pixel:(%+.2e,%+.2e,%+.2e)[mm] %s Deposit:%7.3f MeV CellID:%16lX",
0123 pos.x()/CLHEP::mm,pos.y()/CLHEP::mm,pos.z()/CLHEP::mm,
0124 pos_cell.x()/dd4hep::mm,pos_cell.y()/dd4hep::mm,pos_cell.z()/dd4hep::mm,
0125 pv.name(),h->energyDeposit/CLHEP::MeV,h->cellID
0126 );
0127 delete h;
0128 }
0129 }
0130 return 0;
0131 }
0132
0133 int printParticles(const string& container, vector<Geant4Particle*>* particles) {
0134 typedef vector<Geant4Particle*> _P;
0135 if ( !particles ) {
0136 ::printf("+ Invalid particle container '%s'. No printout\n",container.c_str());
0137 }
0138 else if ( particles->empty() ) {
0139 ::printf("+ Invalid particle container '%s'. No entries. No printout\n",container.c_str());
0140 }
0141 else {
0142 for(_P::const_iterator i=particles->begin(); i!=particles->end(); ++i) {
0143 dd4hep::sim::Geant4ParticleHandle p(*i);
0144 char text[256];
0145 text[0]=0;
0146 if ( p->parents.size() == 1 )
0147 ::snprintf(text,sizeof(text),"/%d",*(p->parents.begin()));
0148 else if ( p->parents.size() > 1 ) {
0149 text[0]='/';text[1]=0;
0150 for(set<int>::const_iterator i=p->parents.begin(); i!=p->parents.end(); ++i)
0151 ::snprintf(text+strlen(text),sizeof(text)-strlen(text),"%d ",*i);
0152 }
0153 printout(ALWAYS,container,
0154 "+++%s %3d stat:%08X PDG:%6d Mom:(%+.2e,%+.2e,%+.2e)[MeV] "
0155 "Vtx:(%+.2e,%+.2e,%+.2e)[mm] #Dau:%3d #Par:%1d%-6s",
0156 "",p->id,p->status,p->pdgID,
0157 p->psx/CLHEP::MeV,p->psy/CLHEP::MeV,p->psz/CLHEP::MeV,
0158 p->vsx/CLHEP::mm,p->vsy/CLHEP::mm,p->vsz/CLHEP::mm,
0159 int(p->daughters.size()),
0160 int(p->parents.size()),
0161 text);
0162 delete (*i);
0163 }
0164 }
0165 return 0;
0166 }
0167 }
0168
0169 int dumpDDG4(const char* fname, int event_num) {
0170 TFile* data = TFile::Open(fname);
0171 if ( !data || data->IsZombie() ) {
0172 printf("+ File seems to not exist. Exiting\n");
0173 usage();
0174 return -1;
0175 }
0176 TTree* tree = (TTree*)data->Get("EVENT");
0177 for(int event=0, num=tree->GetEntries(); event<num; ++event) {
0178 TObjArray* arr = tree->GetListOfBranches();
0179 if ( event_num>= 0 ) event = event_num;
0180 for(int j=0, nj=arr->GetEntries(); j<nj; ++j) {
0181 TBranch* b = (TBranch*)arr->At(j);
0182 typedef vector<void*> _E;
0183 _E* e = 0;
0184 b->SetAddress(&e);
0185 int nbytes = b->GetEvent(event);
0186 if ( nbytes > 0 ) {
0187 if ( e->empty() ) {
0188 continue;
0189 }
0190 string br_name = b->GetName();
0191 string cl_name = b->GetClassName();
0192 if ( cl_name.find("dd4hep::sim::Geant4Tracker::Hit") != string::npos ) {
0193 typedef vector<Geant4Tracker::Hit*> _H;
0194 printHits(br_name,(_H*)e);
0195 }
0196 else if ( cl_name.find("dd4hep::sim::Geant4Calorimeter::Hit") != string::npos ) {
0197 typedef vector<Geant4Calorimeter::Hit*> _H;
0198 printHits(br_name,(_H*)e);
0199 }
0200 else if ( cl_name.find("dd4hep::sim::Geant4Particle") != string::npos ) {
0201 typedef vector<Geant4Particle*> _H;
0202 ::printf("%s\n+ Particle Dump of event %8d [%8d bytes] +\n%s\n",
0203 line,event,nbytes,line);
0204 printParticles(br_name,(_H*)e);
0205 }
0206 }
0207 }
0208 if ( event_num >= 0 ) break;
0209 }
0210 delete data;
0211 return 0;
0212 }
0213
0214 int dumpddg4_load_geometry(const char* fname) {
0215 if ( !have_geometry ) {
0216 have_geometry = true;
0217 gSystem->Load("libDDG4Plugins");
0218 Detector& description = Detector::getInstance();
0219 description.fromXML(fname);
0220 VolumeManager::getVolumeManager();
0221 }
0222 return 1;
0223 }
0224
0225 #if !(defined(G__DICTIONARY) || defined(__CINT__) || defined(__MAKECINT__))
0226 int main(int argc, char** argv) {
0227 int event_num = -1;
0228 const char* fname = 0;
0229 for(int i=1; i<argc;++i) {
0230 if ( argv[i][0]=='-' ) {
0231 if ( strncmp(argv[i],"-input",2) == 0 )
0232 fname = argv[++i];
0233 else if ( strncmp(argv[i],"-compact",2) == 0 )
0234 dumpddg4_load_geometry(argv[++i]);
0235 else if ( strncmp(argv[i],"-event",2) == 0 )
0236 event_num = ::atol(argv[++i]);
0237 }
0238 }
0239 if ( !fname ) {
0240 return usage();
0241 }
0242 return dumpDDG4(fname,event_num);
0243 }
0244 #endif
