ECCE @ EIC Software
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
RunAction.cc
Go to the documentation of this file. Or view the newest version in sPHENIX GitHub for file RunAction.cc
1 //
2 // ********************************************************************
3 // * License and Disclaimer *
4 // * *
5 // * The Geant4 software is copyright of the Copyright Holders of *
6 // * the Geant4 Collaboration. It is provided under the terms and *
7 // * conditions of the Geant4 Software License, included in the file *
8 // * LICENSE and available at http://cern.ch/geant4/license . These *
9 // * include a list of copyright holders. *
10 // * *
11 // * Neither the authors of this software system, nor their employing *
12 // * institutes,nor the agencies providing financial support for this *
13 // * work make any representation or warranty, express or implied, *
14 // * regarding this software system or assume any liability for its *
15 // * use. Please see the license in the file LICENSE and URL above *
16 // * for the full disclaimer and the limitation of liability. *
17 // * *
18 // * This code implementation is the result of the scientific and *
19 // * technical work of the GEANT4 collaboration. *
20 // * By using, copying, modifying or distributing the software (or *
21 // * any work based on the software) you agree to acknowledge its *
22 // * use in resulting scientific publications, and indicate your *
23 // * acceptance of all terms of the Geant4 Software license. *
24 // ********************************************************************
25 //
28 //
29 //
30 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
31 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
32 
33 #include "RunAction.hh"
34 
35 #include "DetectorConstruction.hh"
36 #include "PrimaryGeneratorAction.hh"
37 #include "G4ParticleDefinition.hh"
38 
39 #include "G4Run.hh"
40 #include "G4RunManager.hh"
41 #include "G4UnitsTable.hh"
42 #include "G4EmCalculator.hh"
43 
44 #include "G4Gamma.hh"
45 #include "G4Electron.hh"
46 #include "G4Positron.hh"
47 
48 #include "G4SystemOfUnits.hh"
49 #include "G4PhysicalConstants.hh"
50 #include <iomanip>
51 
52 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
53 
54 RunAction::RunAction(DetectorConstruction* det, PrimaryGeneratorAction* prim)
55 : G4UserRunAction(),
56  fGamma(0), fElectron(0), fPositron(0),
57  fDetector(det), fPrimary(prim), fProcCounter(0), fAnalysisManager(0),
58  fTotalEventCount(0),
59  fPhotonStats(), fElectronStats(), fPositronStats()
60 {
61  fGamma = G4Gamma::Gamma();
62  fElectron = G4Electron::Electron();
63  fPositron = G4Positron::Positron();
64 
65  BookHisto();
66 }
67 
68 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
69 
70 RunAction::~RunAction()
71 {}
72 
73 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
74 
75 void RunAction::BeginOfRunAction(const G4Run* aRun)
76 {
77  G4cout << "### Run " << aRun->GetRunID() << " start." << G4endl;
78 
79  // save Rndm status
80  // G4RunManager::GetRunManager()->SetRandomNumberStore(false);
81  // CLHEP::HepRandom::showEngineStatus();
82 
83  if (fProcCounter) delete fProcCounter;
84  fProcCounter = new ProcessesCount;
85  fTotalEventCount = 0;
86  fPhotonStats.Clear();
87  fElectronStats.Clear();
88  fPositronStats.Clear();
89 }
90 
91 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
92 
93 void RunAction::FillData(const G4ParticleDefinition* particle,
94  G4double kinEnergy, G4double costheta,
95  G4double phi,
96  G4double longitudinalPolarization)
97 {
98  G4int id = -1;
99  if (particle == fGamma) {
100  fPhotonStats.FillData(kinEnergy, costheta, longitudinalPolarization);
101  if(fAnalysisManager) { id = 1; }
102  } else if (particle == fElectron) {
103  fElectronStats.FillData(kinEnergy, costheta, longitudinalPolarization);
104  if(fAnalysisManager) { id = 5; }
105  } else if (particle == fPositron) {
106  fPositronStats.FillData(kinEnergy, costheta, longitudinalPolarization);
107  if(fAnalysisManager) { id = 9; }
108  }
109  if(id > 0) {
110  fAnalysisManager->FillH1(id,kinEnergy,1.0);
111  fAnalysisManager->FillH1(id+1,costheta,1.0);
112  fAnalysisManager->FillH1(id+2,phi,1.0);
113  fAnalysisManager->FillH1(id+3,longitudinalPolarization,1.0);
114  }
115 }
116 
117 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
118 
119 void RunAction::BookHisto()
120 {
121  // Always creating analysis manager
122  fAnalysisManager = G4AnalysisManager::Instance();
123  fAnalysisManager->SetActivation(true);
124  fAnalysisManager->SetVerboseLevel(1);
125 
126  // Open file histogram file
127  fAnalysisManager->OpenFile("pol01");
128  fAnalysisManager->SetFirstHistoId(1);
129 
130  // Creating an 1-dimensional histograms in the root directory of the tree
131  const G4String id[] = { "h1", "h2", "h3", "h4", "h5",
132  "h6", "h7", "h8", "h9", "h10", "h11", "h12"};
133  const G4String title[] =
134  { "Gamma Energy distribution", //1
135  "Gamma Cos(Theta) distribution", //2
136  "Gamma Phi angular distribution", //3
137  "Gamma longitudinal Polarization", //4
138  "Electron Energy distribution", //5
139  "Electron Cos(Theta) distribution", //6
140  "Electron Phi angular distribution", //7
141  "Electron longitudinal Polarization", //8
142  "Positron Energy distribution", //9
143  "Positron Cos(Theta) distribution", //10
144  "Positron Phi angular distribution", //11
145  "Positron longitudinal Polarization" //12
146  };
147  G4double vmin, vmax;
148  G4int nbins = 120;
149  for(int i=0; i<12; ++i) {
150  G4int j = i - i/4*4;
151  if(0==j) { vmin = 0.; vmax = 12.*MeV; }
152  else if(1==j) { vmin = -1.; vmax = 1.; }
153  else if(2==j) { vmin = 0.; vmax = pi; }
154  else { vmin = -1.5; vmax = 1.5; }
155  G4int ih = fAnalysisManager->CreateH1(id[i],title[i],nbins,vmin,vmax);
156  fAnalysisManager->SetH1Activation(ih, false);
157  }
158 }
159 
160 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
161 
162 void RunAction::SaveHisto(G4int nevents)
163 {
164  if(fAnalysisManager) {
165  G4double norm = 1.0/G4double(nevents);
166  for(int i=0; i<12; ++i) {
167  fAnalysisManager->ScaleH1(i, norm);
168  }
169  fAnalysisManager->Write();
170  fAnalysisManager->CloseFile();
171  delete fAnalysisManager;
172  fAnalysisManager = 0;
173  }
174 }
175 
176 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
177 
178 void RunAction::CountProcesses(G4String procName)
179 {
180  // is the process already counted ?
181  // *AS* change to std::map?!
182  size_t nbProc = fProcCounter->size();
183  size_t i = 0;
184  while ((i<nbProc)&&((*fProcCounter)[i]->GetName()!=procName)) i++;
185  if (i == nbProc) fProcCounter->push_back( new ProcessCount(procName));
186 
187  (*fProcCounter)[i]->Count();
188 }
189 
190 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
191 
192 void RunAction::EndOfRunAction(const G4Run* aRun)
193 {
194  G4int NbOfEvents = aRun->GetNumberOfEvent();
195  if (NbOfEvents == 0) return;
196 
197  G4int prec = G4cout.precision(5);
198 
199  G4Material* material = fDetector->GetMaterial();
200  G4double density = material->GetDensity();
201 
202  G4ParticleDefinition* particle =
203  fPrimary->GetParticleGun()->GetParticleDefinition();
204  G4String Particle = particle->GetParticleName();
205  G4double energy = fPrimary->GetParticleGun()->GetParticleEnergy();
206  G4cout << "\n The run consists of " << NbOfEvents << " "<< Particle << " of "
207  << G4BestUnit(energy,"Energy") << " through "
208  << G4BestUnit(fDetector->GetBoxSizeZ(),"Length") << " of "
209  << material->GetName() << " (density: "
210  << G4BestUnit(density,"Volumic Mass") << ")" << G4endl;
211 
212  //frequency of processes
213  G4cout << "\n Process calls frequency --->\n";
214  for (size_t i=0; i< fProcCounter->size();i++) {
215  G4String procName = (*fProcCounter)[i]->GetName();
216  G4int count = (*fProcCounter)[i]->GetCounter();
217  G4cout << "\t" << procName << " = " << count<<"\n";
218  }
219 
220  if (fTotalEventCount == 0) return;
221 
222  G4cout<<" Gamma: \n";
223  fPhotonStats.PrintResults(fTotalEventCount);
224  G4cout<<" Electron: \n";
225  fElectronStats.PrintResults(fTotalEventCount);
226  G4cout<<" Positron: \n";
227  fPositronStats.PrintResults(fTotalEventCount);
228 
229  //cross check from G4EmCalculator
230  // G4cout << "\n Verification from G4EmCalculator. \n";
231  // G4EmCalculator emCal;
232 
233  //restore default format
234  G4cout.precision(prec);
235 
236  // write out histograms
237  SaveHisto(NbOfEvents);
238 
239  // show Rndm status
240  CLHEP::HepRandom::showEngineStatus();
241 }
242 
243 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
244 
245 void RunAction::EventFinished()
246 {
247  ++fTotalEventCount;
248  fPhotonStats.EventFinished();
249  fElectronStats.EventFinished();
250  fPositronStats.EventFinished();
251 }
252 
253 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
254 
255 RunAction::ParticleStatistics::ParticleStatistics()
256  : fCurrentNumber(0),
257  fTotalNumber(0), fTotalNumber2(0),
258  fSumEnergy(0), fSumEnergy2(0),
259  fSumPolarization(0), fSumPolarization2(0),
260  fSumCosTheta(0), fSumCosTheta2(0)
261 {}
262 
263 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
264 
265 RunAction::ParticleStatistics::~ParticleStatistics()
266 {}
267 
268 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
269 
270 void RunAction::ParticleStatistics::EventFinished()
271 {
272  fTotalNumber+=fCurrentNumber;
273  fTotalNumber2+=fCurrentNumber*fCurrentNumber;
274  fCurrentNumber=0;
275 }
276 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
277 
278 void RunAction::ParticleStatistics:: FillData(G4double kinEnergy,
279  G4double costheta,
280  G4double longitudinalPolarization)
281 {
282  ++fCurrentNumber;
283  fSumEnergy+=kinEnergy;
284  fSumEnergy2+=kinEnergy*kinEnergy;
285  fSumPolarization+=longitudinalPolarization;
286  fSumPolarization2+=longitudinalPolarization*longitudinalPolarization;
287  fSumCosTheta+=costheta;
288  fSumCosTheta2+=costheta*costheta;
289 }
290 
291 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
292 
293 void RunAction::ParticleStatistics::PrintResults(G4int totalNumberOfEvents)
294 {
295  G4cout<<"Mean Number per Event :"
296  <<G4double(fTotalNumber)/G4double(totalNumberOfEvents)<<"\n";
297  if (fTotalNumber==0) fTotalNumber=1;
298  G4double energyMean=fSumEnergy/fTotalNumber;
299  G4double energyRms=std::sqrt(fSumEnergy2/fTotalNumber-energyMean*energyMean);
300  G4cout<<"Mean Energy :"<< G4BestUnit(energyMean,"Energy")
301  <<" +- "<<G4BestUnit(energyRms,"Energy")<<"\n";
302  G4double polarizationMean=fSumPolarization/fTotalNumber;
303  G4double polarizationRms=
304  std::sqrt(fSumPolarization2/fTotalNumber-polarizationMean*polarizationMean);
305  G4cout<<"Mean Polarization :"<< polarizationMean
306  <<" +- "<<polarizationRms<<"\n";
307 }
308 
309 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
310 
311 void RunAction::ParticleStatistics::Clear()
312 {
313  fCurrentNumber=0;
314  fTotalNumber=fTotalNumber2=0;
315  fSumEnergy=fSumEnergy2=0;
316  fSumPolarization=fSumPolarization2=0;
317  fSumCosTheta=fSumCosTheta2=0;
318 }
319 
320 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......