test
CMS 3D CMS Logo

All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Properties Friends Macros Pages
G4StepStatistics.h
Go to the documentation of this file.
1 #ifndef HelpfulWatchers_G4StepStatistics_h
2 #define HelpfulWatchers_G4StepStatistics_h
3 // -*- C++ -*-
4 //
5 // Package: HelpfulWatchers
6 // Class : SimTracer
7 //
16 //
17 // Original Author:
18 // Created: Tue Nov 22 16:41:33 EST 2005
19 //
20 
21 // system include files
22 #include <iostream>
23 
24 // user include files
28 #include "G4Step.hh"
29 #include "G4VProcess.hh"
30 #include "G4ParticleDefinition.hh"
31 #include <map>
34 
35 #include <TROOT.h>
36 #include <TTree.h>
37 #include <TFile.h>
38 #include <TVector.h>
39 #include <TString.h>
40 #include <TClonesArray.h>
41 //#include<TObjString.h>
42 
43 // forward declarations
44 class DDDWorld;
45 class BeginOfJob;
46 class BeginOfRun;
47 class BeginOfEvent;
48 class BeginOfTrack;
49 class G4Step;
50 
51 class EndOfRun;
52 class EndOfEvent;
53 class EndOfTrack;
54 
55 //Define a class StepID
56 class StepID {
57 
58  private:
59  //G4 Region
60  G4String theG4RegionName;
61  //G4 Physical Process
62  G4String theG4ProcessName;
63  //Particle PDG ID
65 
66  public:
67  //Constructor using G4Step
68  StepID(const G4Step* theG4Step)
69  :
70  theG4RegionName("UNDEFINED"),
71  theG4ProcessName("UNDEFINED"),
72  theParticlePDGID(theG4Step->GetTrack()->GetDefinition()->GetPDGEncoding())
73  {
74  std::cout<<"Start"<<std::endl;
75  if (theG4Step->GetPreStepPoint()->GetPhysicalVolume()) {
76  theG4RegionName = theG4Step->GetPreStepPoint()->GetPhysicalVolume()->GetLogicalVolume()->GetRegion()->GetName();
77  }
78  std::cout<<"Physical Volume"<<std::endl;
79  if (theG4Step->GetPreStepPoint()->GetProcessDefinedStep()){
80  theG4ProcessName = theG4Step->GetPreStepPoint()->GetProcessDefinedStep()->GetProcessName();
81  }
82  std::cout<<"Process Name"<<std::endl;
83  }
84 
85  //Getters
86  G4String GetRegionName () const { return theG4RegionName; }
87  G4String GetProcessName () const { return theG4ProcessName; }
88  G4int GetParticlePDGID () const { return theParticlePDGID; }
89 
90  //Comparison Operators (necessary in order to use StepID as a key in a map)
91  bool operator==(const StepID& id) const{
92  return ( strcmp(theG4RegionName,id.GetRegionName())==0 && strcmp(theG4ProcessName,id.GetProcessName())==0 && theParticlePDGID==id.GetParticlePDGID() ) ? true : false;
93  }
94 
95  bool operator<(const StepID& id) const
96  {
98  return (theParticlePDGID > id.GetParticlePDGID());
99  }
100  else if (strcmp(theG4RegionName,id.GetRegionName())!=0){
101  return strcmp(theG4RegionName,id.GetRegionName())>0 ? true : false;
102  }
103  else if (strcmp(theG4ProcessName,id.GetProcessName())!=0){
104  return strcmp(theG4ProcessName,id.GetProcessName())>0 ? true : false;
105  }
106  else {//The case in which they are all equal!
107  return false;
108  }
109  }
110 
111  bool operator>(const StepID& id) const
112  {
114  return (theParticlePDGID < id.GetParticlePDGID());
115  }
116  else if(strcmp(theG4RegionName,id.GetRegionName())!=0){
117  return strcmp(theG4RegionName,id.GetRegionName())<0 ? true : false;
118  }
119  else if (strcmp(theG4ProcessName,id.GetProcessName())!=0){
120  return strcmp(theG4ProcessName,id.GetProcessName())<0 ? true : false;
121  }
122  else {//The case in which they are all equal!
123  return false;
124  }
125  }
126 };
127 
128 #define OBSERVES(type) public Observer<const type*>
129 #define UPDATE(type) void update(const type*) { std::cout <<"++ signal " #type<<std::endl; }
130 class G4StepStatistics : public SimWatcher,
131 OBSERVES(DDDWorld),
132 OBSERVES(BeginOfJob),
133 OBSERVES(BeginOfRun),
134 OBSERVES(BeginOfEvent),
135 OBSERVES(BeginOfTrack),
136 OBSERVES(G4Step),
137 OBSERVES(EndOfRun),
138 OBSERVES(EndOfEvent),
139 OBSERVES(EndOfTrack)
140 {
141  public:
143  m_verbose(pSet.getUntrackedParameter<bool>("verbose",false)),
144  Event(0)
145  {
146  //Adding TFile Service output
147  G4StepTree = fs->make<TTree>("G4StepTree","G4Step Tree ");
148  G4StepTree->Branch("Event",&Event,"Event/I");
149  G4StepTree->Branch("PDGID",&PDGID,"PDGID[100000]/I");
150  Region = new TClonesArray("TObjString",100000);
151  G4StepTree->Branch("Region",&Region);
152  Process = new TClonesArray("TObjString",100000);
153  G4StepTree->Branch("Process",&Process);
154  G4StepTree->Branch("G4StepFreq",&G4StepFreq,"G4StepFreq[100000]/I");
155  }
159  // void update(const BeginOfRun* iRun) {
160  //std::cout <<"++ signal BeginOfRun " <<std::endl;
161  //}
164  void update(const G4Step* iStep) {
165  std::cout <<"++ signal G4Step " ;
166  //Dump the relevant information from the G4 step in the object mysteptest
167  StepID mysteptest(iStep);
168  //Add the StepID to the map, or increment if it already exists:
169  if ( G4StatsMap.find(mysteptest) == G4StatsMap.end() )
170  {
171  //Allocating new memory for a pointer to associate with the key mysteptest
172  //in our map. Initializing it to 1,will be incremented working on the value of the pointer.
173  unsigned int* MyValue = new unsigned int(1);
174  //Inserting the new key,value pair
175  G4StatsMap.insert(std::make_pair(mysteptest, MyValue));
176  }
177  else
178  {
179  //Incrementing the value of the pointer by 1
180  *G4StatsMap[mysteptest] = *G4StatsMap[mysteptest] + 1;
181  }
182 
183  //If the verbose flag is set, then dump the information
184  if (m_verbose) {
185  std::cout << " StepID RegionName: "<< mysteptest.GetRegionName();
186  std::cout << " StepID ProcessName: "<< mysteptest.GetProcessName();
187  std::cout << " StepID ParticlePDGID: "<< mysteptest.GetParticlePDGID();
188  }
189  std::cout<<std::endl;
190 }
191 //UPDATE(G4Step)
193  // void update(const EndOfRun* iRun) {
194  //std::cout <<"++ signal EndOfRun " <<std::endl;
195  //}
196 
197  //UPDATE(EndOfEvent)
199  std::cout <<"++ signal EndOfEvent " <<std::endl;
200  Event++;
201 
202  //Dumping the map in the log if verbose is chosen:
203  if(m_verbose) {
204  std::cout <<" G4StatsMap size is: "<<G4StatsMap.size()<<std::endl;
205  }
206  int index(0);
207  for (std::map<const StepID,unsigned int*>::const_iterator step = G4StatsMap.begin(); step != G4StatsMap.end(); ++step, ++index){
208  if(m_verbose) {
209  std::cout <<" G4StatsMap step is: "<<step->first.GetRegionName()<<" "<<step->first.GetProcessName()<<" "<<step->first.GetParticlePDGID();//<<" "<<step->first.GetTrackID() ;
210  std::cout <<" Number of such steps: "<< *step->second <<std::endl;
211  }
212  //Rolling the map into 5 "arrays", containing the StepID information and the G4Step statistics
213  PDGID[index]=step->first.GetParticlePDGID();
214  new ((*Region)[index]) TObjString (step->first.GetRegionName());
215  new ((*Process)[index]) TObjString (step->first.GetProcessName());
216  G4StepFreq[index]=*step->second;
217  }
218 
219  G4StepTree->Fill();
220 }
222 
223  private:
224 
225  bool m_verbose;
226 
227 //Adding the G4StatsMap to keep track of statistics in terms of step information...
228  std::map<const StepID,unsigned int*> G4StatsMap;
230  TTree* G4StepTree;
231  unsigned int Event;
232  Int_t PDGID[100000];
233  TClonesArray* Region;
234  TClonesArray* Process;
235  Int_t G4StepFreq[100000];
236 };
237 
238 #endif
std::map< const StepID, unsigned int * > G4StatsMap
Int_t G4StepFreq[100000]
bool operator==(const StepID &id) const
G4StepStatistics(const edm::ParameterSet &pSet)
bool operator<(const StepID &id) const
T * make(const Args &...args) const
make new ROOT object
Definition: TFileService.h:64
G4String theG4RegionName
void update(const DDDWorld *)
This routine will be called when the appropriate signal arrives.
StepID(const G4Step *theG4Step)
G4String theG4ProcessName
Definition: Event.h:16
G4String GetProcessName() const
G4int theParticlePDGID
Int_t PDGID[100000]
bool operator>(const StepID &id) const
TClonesArray * Process
#define OBSERVES(type)
G4int GetParticlePDGID() const
string const
Definition: compareJSON.py:14
TClonesArray * Region
tuple cout
Definition: gather_cfg.py:145
volatile std::atomic< bool > shutdown_flag false
#define UPDATE(type)
edm::Service< TFileService > fs
G4String GetRegionName() const