File indexing completed on 2024-04-06 12:30:25
0001 #ifndef HelpfulWatchers_G4StepStatistics_h
0002 #define HelpfulWatchers_G4StepStatistics_h
0003
0004
0005
0006
0007
0008
0009
0010
0011
0012
0013
0014
0015
0016
0017
0018
0019
0020
0021
0022 #include <iostream>
0023
0024
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
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
0056 class StepID {
0057 private:
0058
0059 G4String theG4RegionName;
0060
0061 G4String theG4ProcessName;
0062
0063 G4int theParticlePDGID;
0064
0065 public:
0066
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
0083 G4String GetRegionName() const { return theG4RegionName; }
0084 G4String GetProcessName() const { return theG4ProcessName; }
0085 G4int GetParticlePDGID() const { return theParticlePDGID; }
0086
0087
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 {
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 {
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
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
0152
0153
0154 UPDATE(BeginOfEvent)
0155 UPDATE(BeginOfTrack)
0156 void update(const G4Step *iStep) override {
0157 std::cout << "++ signal G4Step ";
0158
0159 StepID mysteptest(iStep);
0160
0161 if (G4StatsMap.find(mysteptest) == G4StatsMap.end()) {
0162
0163
0164
0165 unsigned int *MyValue = new unsigned int(1);
0166
0167 G4StatsMap.insert(std::make_pair(mysteptest, MyValue));
0168 } else {
0169
0170 *G4StatsMap[mysteptest] = *G4StatsMap[mysteptest] + 1;
0171 }
0172
0173
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
0182 UPDATE(EndOfRun)
0183
0184
0185
0186
0187
0188 void update(const EndOfEvent *iRun) override {
0189 std::cout << "++ signal EndOfEvent " << std::endl;
0190 Event++;
0191
0192
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();
0202 std::cout << " Number of such steps: " << *step->second << std::endl;
0203 }
0204
0205
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
0220
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