HelpfulWatchers/src/G4StepStatistics.h

0001 #ifndef HelpfulWatchers_G4StepStatistics_h
0002 #define HelpfulWatchers_G4StepStatistics_h
0003 // -*- C++ -*-
0004 //
0005 // Package:     HelpfulWatchers
0006 // Class  :     SimTracer
0007 //
0008 /**\class SimTracer SimTracer.h SimG4Core/HelpfulWatchers/interface/SimTracer.h
0009
0010  Description: Prints a message for each Oscar signal
0011
0012  Usage:
0013     <usage>
0014
0015 */
0016 //
0017 // Original Author:
0018 //         Created:  Tue Nov 22 16:41:33 EST 2005
0019 //
0020
0021 // system include files
0022 #include <iostream>
0023
0024 // user include files
0025 #include "FWCore/ParameterSet/interface/ParameterSet.h"
0026 #include "G4ParticleDefinition.hh"
0027 #include "G4Step.hh"
0028 #include "G4VProcess.hh"
0029 #include "SimG4Core/Notification/interface/Observer.h"
0030 #include "SimG4Core/Watcher/interface/SimWatcher.h"
0031 #include <CommonTools/UtilAlgos/interface/TFileService.h>
0032 #include <FWCore/ServiceRegistry/interface/Service.h>
0033 #include <map>
0034
0035 #include <TClonesArray.h>
0036 #include <TFile.h>
0037 #include <TROOT.h>
0038 #include <TString.h>
0039 #include <TTree.h>
0040 #include <TVector.h>
0041 #include <TObjString.h>
0042
0043 // forward declarations
0044 class DDDWorld;
0045 class BeginOfJob;
0046 class BeginOfRun;
0047 class BeginOfEvent;
0048 class BeginOfTrack;
0049 class G4Step;
0050
0051 class EndOfRun;
0052 class EndOfEvent;
0053 class EndOfTrack;
0054
0055 // Define a class StepID
0056 class StepID {
0057 private:
0058   // G4 Region
0059   G4String theG4RegionName;
0060   // G4 Physical Process
0061   G4String theG4ProcessName;
0062   // Particle PDG ID
0063   G4int theParticlePDGID;
0064
0065 public:
0066   // Constructor using G4Step
0067   StepID(const G4Step *theG4Step)
0068       : theG4RegionName("UNDEFINED"),
0069         theG4ProcessName("UNDEFINED"),
0070         theParticlePDGID(theG4Step->GetTrack()->GetDefinition()->GetPDGEncoding()) {
0071     std::cout << "Start" << std::endl;
0072     if (theG4Step->GetPreStepPoint()->GetPhysicalVolume()) {
0073       theG4RegionName = theG4Step->GetPreStepPoint()->GetPhysicalVolume()->GetLogicalVolume()->GetRegion()->GetName();
0074     }
0075     std::cout << "Physical Volume" << std::endl;
0076     if (theG4Step->GetPreStepPoint()->GetProcessDefinedStep()) {
0077       theG4ProcessName = theG4Step->GetPreStepPoint()->GetProcessDefinedStep()->GetProcessName();
0078     }
0079     std::cout << "Process Name" << std::endl;
0080   }
0081
0082   // Getters
0083   G4String GetRegionName() const { return theG4RegionName; }
0084   G4String GetProcessName() const { return theG4ProcessName; }
0085   G4int GetParticlePDGID() const { return theParticlePDGID; }
0086
0087   // Comparison Operators (necessary in order to use StepID as a key in a map)
0088   bool operator==(const StepID &id) const {
0089     return (strcmp(theG4RegionName, id.GetRegionName()) == 0 && strcmp(theG4ProcessName, id.GetProcessName()) == 0 &&
0090             theParticlePDGID == id.GetParticlePDGID())
0091                ? true
0092                : false;
0093   }
0094
0095   bool operator<(const StepID &id) const {
0096     if (theParticlePDGID != id.GetParticlePDGID()) {
0097       return (theParticlePDGID > id.GetParticlePDGID());
0098     } else if (strcmp(theG4RegionName, id.GetRegionName()) != 0) {
0099       return strcmp(theG4RegionName, id.GetRegionName()) > 0 ? true : false;
0100     } else if (strcmp(theG4ProcessName, id.GetProcessName()) != 0) {
0101       return strcmp(theG4ProcessName, id.GetProcessName()) > 0 ? true : false;
0102     } else {  // The case in which they are all equal!
0103       return false;
0104     }
0105   }
0106
0107   bool operator>(const StepID &id) const {
0108     if (theParticlePDGID != id.GetParticlePDGID()) {
0109       return (theParticlePDGID < id.GetParticlePDGID());
0110     } else if (strcmp(theG4RegionName, id.GetRegionName()) != 0) {
0111       return strcmp(theG4RegionName, id.GetRegionName()) < 0 ? true : false;
0112     } else if (strcmp(theG4ProcessName, id.GetProcessName()) != 0) {
0113       return strcmp(theG4ProcessName, id.GetProcessName()) < 0 ? true : false;
0114     } else {  // The case in which they are all equal!
0115       return false;
0116     }
0117   }
0118 };
0119
0120 #define OBSERVES(type) \
0121 public                 \
0122   Observer<const type *>
0123 #define UPDATE(type) \
0124   void update(const type *) override { std::cout << "++ signal " #type << std::endl; }
0125 class G4StepStatistics : public SimWatcher,
0126                          OBSERVES(DDDWorld),
0127                          OBSERVES(BeginOfJob),
0128                          OBSERVES(BeginOfRun),
0129                          OBSERVES(BeginOfEvent),
0130                          OBSERVES(BeginOfTrack),
0131                          OBSERVES(G4Step),
0132                          OBSERVES(EndOfRun),
0133                          OBSERVES(EndOfEvent),
0134                          OBSERVES(EndOfTrack) {
0135 public:
0136   G4StepStatistics(const edm::ParameterSet &pSet)
0137       : m_verbose(pSet.getUntrackedParameter<bool>("verbose", false)), Event(0) {
0138     // Adding TFile Service output
0139     G4StepTree = fs->make<TTree>("G4StepTree", "G4Step Tree ");
0140     G4StepTree->Branch("Event", &Event, "Event/I");
0141     G4StepTree->Branch("PDGID", &PDGID, "PDGID[100000]/I");
0142     Region = new TClonesArray("TObjString", 100000);
0143     G4StepTree->Branch("Region", &Region);
0144     Process = new TClonesArray("TObjString", 100000);
0145     G4StepTree->Branch("Process", &Process);
0146     G4StepTree->Branch("G4StepFreq", &G4StepFreq, "G4StepFreq[100000]/I");
0147   }
0148   UPDATE(DDDWorld)
0149   UPDATE(BeginOfJob)
0150   UPDATE(BeginOfRun)
0151   //    void update(const BeginOfRun* iRun) {
0152   // std::cout <<"++ signal BeginOfRun " <<std::endl;
0153   //}
0154   UPDATE(BeginOfEvent)
0155   UPDATE(BeginOfTrack)
0156   void update(const G4Step *iStep) override {
0157     std::cout << "++ signal G4Step ";
0158     // Dump the relevant information from the G4 step in the object mysteptest
0159     StepID mysteptest(iStep);
0160     // Add the StepID to the map, or increment if it already exists:
0161     if (G4StatsMap.find(mysteptest) == G4StatsMap.end()) {
0162       // Allocating new memory for a pointer to associate with the key
0163       // mysteptest in our map. Initializing it to 1,will be incremented working
0164       // on the value of the pointer.
0165       unsigned int *MyValue = new unsigned int(1);
0166       // Inserting the new key,value pair
0167       G4StatsMap.insert(std::make_pair(mysteptest, MyValue));
0168     } else {
0169       // Incrementing the value of the pointer by 1
0170       *G4StatsMap[mysteptest] = *G4StatsMap[mysteptest] + 1;
0171     }
0172
0173     // If the verbose flag is set, then dump the information
0174     if (m_verbose) {
0175       std::cout << " StepID RegionName: " << mysteptest.GetRegionName();
0176       std::cout << " StepID ProcessName: " << mysteptest.GetProcessName();
0177       std::cout << " StepID ParticlePDGID: " << mysteptest.GetParticlePDGID();
0178     }
0179     std::cout << std::endl;
0180   }
0181   // UPDATE(G4Step)
0182   UPDATE(EndOfRun)
0183   //  void update(const EndOfRun* iRun) {
0184   // std::cout <<"++ signal EndOfRun " <<std::endl;
0185   //}
0186
0187   // UPDATE(EndOfEvent)
0188   void update(const EndOfEvent *iRun) override {
0189     std::cout << "++ signal EndOfEvent " << std::endl;
0190     Event++;
0191
0192     // Dumping the map in the log if verbose is chosen:
0193     if (m_verbose) {
0194       std::cout << " G4StatsMap size is: " << G4StatsMap.size() << std::endl;
0195     }
0196     int index(0);
0197     for (std::map<const StepID, unsigned int *>::const_iterator step = G4StatsMap.begin(); step != G4StatsMap.end();
0198          ++step, ++index) {
0199       if (m_verbose) {
0200         std::cout << " G4StatsMap step is: " << step->first.GetRegionName() << " " << step->first.GetProcessName()
0201                   << " " << step->first.GetParticlePDGID();  //<<" "<<step->first.GetTrackID() ;
0202         std::cout << " Number of such steps: " << *step->second << std::endl;
0203       }
0204       // Rolling the map into 5 "arrays", containing the StepID information and
0205       // the G4Step statistics
0206       PDGID[index] = step->first.GetParticlePDGID();
0207       new ((*Region)[index]) TObjString(step->first.GetRegionName());
0208       new ((*Process)[index]) TObjString(step->first.GetProcessName());
0209       G4StepFreq[index] = *step->second;
0210     }
0211
0212     G4StepTree->Fill();
0213   }
0214   UPDATE(EndOfTrack)
0215
0216 private:
0217   bool m_verbose;
0218
0219   // Adding the G4StatsMap to keep track of statistics in terms of step
0220   // information...
0221   std::map<const StepID, unsigned int *> G4StatsMap;
0222   edm::Service<TFileService> fs;
0223   TTree *G4StepTree;
0224   unsigned int Event;
0225   Int_t PDGID[100000];
0226   TClonesArray *Region;
0227   TClonesArray *Process;
0228   Int_t G4StepFreq[100000];
0229 };
0230
0231 #endif