Geant4  9.6.p02
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
RunAction.cc
Go to the documentation of this file.
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 //
26 // $Id$
27 //
28 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
29 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
30 
31 #include "RunAction.hh"
32 
33 #include "PrimaryGeneratorAction.hh"
34 #include "HistoManager.hh"
35 
36 #include "G4Run.hh"
37 #include "G4RunManager.hh"
38 #include "G4UnitsTable.hh"
39 #include "G4Geantino.hh"
40 
41 #include "Randomize.hh"
42 
43 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
44 
45 RunAction::RunAction(DetectorConstruction* det, PrimaryGeneratorAction* prim,
46  HistoManager* hist)
47 :detector(det), primary(prim), histoManager(hist)
48 {
49  writeFile = false;
50 }
51 
52 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
53 
54 RunAction::~RunAction()
55 { }
56 
57 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
58 
59 void RunAction::BeginOfRunAction(const G4Run* aRun)
60 {
61  G4cout << "### Run " << aRun->GetRunID() << " start." << G4endl;
62 
63  // save Rndm status
64  //
65  G4RunManager::GetRunManager()->SetRandomNumberStore(true);
66  CLHEP::HepRandom::showEngineStatus();
67 
68  //initialize cumulative quantities
69  //
70  G4int nbPixels = detector->GetSizeVectorPixels();
71  G4int size = totalEnergy.size();
72  if (size < nbPixels) {
73  visibleEnergy.resize(nbPixels);
74  totalEnergy.resize(nbPixels);
75 
76  visibleEnergy2.resize(nbPixels);
77  totalEnergy2.resize(nbPixels);
78  }
79 
80  for (G4int k=0; k<nbPixels; k++) {
81  visibleEnergy[k] = visibleEnergy2[k] = totalEnergy[k]= totalEnergy2[k] = 0.0;
82  }
83 
84  G4int n1pxl = detector->GetN1Pixels();
85  size = layerEtot.size();
86  if (size < n1pxl) {
87  layerEvis.resize(n1pxl);
88  layerEtot.resize(n1pxl);
89  layerEvis2.resize(n1pxl);
90  layerEtot2.resize(n1pxl);
91  }
92 
93  for (G4int k=0; k<n1pxl; k++) {
94  layerEvis[k] = layerEvis2[k] = layerEtot[k]= layerEtot2[k] = 0.0;
95  }
96 
97  nbEvents = 0;
98  calorEvis = calorEvis2 = calorEtot = calorEtot2 = Eleak = Eleak2 = 0.;
99  EdLeak[0] = EdLeak[1] = EdLeak[2] = 0.;
100  nbRadLen = nbRadLen2 = 0.;
101 
102  //histograms
103  //
104  histoManager->book();
105 
106  //create ascii file for pixels
107  //
108  if (writeFile) CreateFilePixels();
109 }
110 
111 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
112 
113 void RunAction::fillPerEvent_1(G4int pixel, G4double Evis, G4double Etot)
114 {
115  //accumulate statistic per pixel
116  //
117  visibleEnergy[pixel] += Evis; visibleEnergy2[pixel] += Evis*Evis;
118  totalEnergy[pixel] += Etot; totalEnergy2[pixel] += Etot*Etot;
119 }
120 
121 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
122 
123 void RunAction::fillPerEvent_2(G4int layer, G4double Evis, G4double Etot)
124 {
125  //accumulate statistic per layer
126  //
127  layerEvis[layer] += Evis; layerEvis2[layer] += Evis*Evis;
128  layerEtot[layer] += Etot; layerEtot2[layer] += Etot*Etot;
129 }
130 
131 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
132 
133 void RunAction::fillPerEvent_3(G4double calEvis, G4double calEtot,
134  G4double eleak)
135 {
136  //accumulate statistic
137  //
138  nbEvents++;
139  calorEvis += calEvis; calorEvis2 += calEvis*calEvis;
140  calorEtot += calEtot; calorEtot2 += calEtot*calEtot;
141  Eleak += eleak; Eleak2 += eleak*eleak;
142 }
143 
144 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
145 
146 void RunAction::fillDetailedLeakage(G4int icase, G4double energy)
147 {
148  //forward, backward, lateral leakage
149  //
150  EdLeak[icase] += energy;
151 }
152 
153 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
154 
155 void RunAction::fillNbRadLen(G4double dn)
156 {
157  //total number of radiation length
158  //
159  nbRadLen += dn; nbRadLen2 += dn*dn;
160 }
161 
162 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
163 
164 
165 void RunAction::EndOfRunAction(const G4Run*)
166 {
167  //calorimeter
168  //
169  detector->PrintCalorParameters();
170 
171  //run conditions
172  //
173  G4ParticleDefinition* particle = primary->GetParticleGun()
174  ->GetParticleDefinition();
175  G4String partName = particle->GetParticleName();
176  G4double energy = primary->GetParticleGun()->GetParticleEnergy();
177 
178  G4int prec = G4cout.precision(3);
179 
180  G4cout << " The run was " << nbEvents << " " << partName << " of "
181  << G4BestUnit(energy,"Energy") << " through the calorimeter" << G4endl;
182 
183  G4cout << "------------------------------------------------------------"
184  << G4endl;
185 
186  //if no events, return
187  //
188  if (nbEvents == 0) return;
189 
190  //compute and print statistic
191  //
192  std::ios::fmtflags mode = G4cout.flags();
193 
194  //number of radiation length
195  //
196  if (particle == G4Geantino::Geantino() ) {
197  G4double meanNbRadL = nbRadLen/ nbEvents;
198  G4double meanNbRadL2 = nbRadLen2/nbEvents;
199  G4double varNbRadL = meanNbRadL2 - meanNbRadL*meanNbRadL;
200  G4double rmsNbRadL = 0.;
201  if (varNbRadL > 0.) rmsNbRadL = std::sqrt(varNbRadL);
202  G4double effRadL = (detector->GetCalorThickness())/meanNbRadL;
203  G4cout.precision(5);
204  G4cout
205  << "\n Calor : mean number of Rad Length = "
206  << meanNbRadL << " +- "<< rmsNbRadL
207  << " --> Effective Rad Length = "
208  << G4BestUnit( effRadL,"Length") << G4endl;
209 
210  G4cout << "------------------------------------------------------------"
211  << G4endl;
212  }
213 
214  G4cout.precision(prec);
215  G4cout << "\n "
216  << "visible Energy (rms/mean) "
217  << "total Energy (rms/mean)" << G4endl;
218 
219  G4double meanEvis,meanEvis2,varianceEvis,rmsEvis,resEvis;
220  G4double meanEtot,meanEtot2,varianceEtot,rmsEtot,resEtot;
221 
222  G4int n1pxl = detector->GetN1Pixels();
223 
224  for (G4int i1=0; i1<n1pxl; i1++) {
225  //visible energy
226  meanEvis = layerEvis[i1] /nbEvents;
227  meanEvis2 = layerEvis2[i1]/nbEvents;
228  varianceEvis = meanEvis2 - meanEvis*meanEvis;
229  rmsEvis = 0.;
230  if (varianceEvis > 0.) rmsEvis = std::sqrt(varianceEvis);
231  resEvis = 100*rmsEvis/meanEvis;
232  histoManager->FillHisto(3, i1+0.5, meanEvis);
233 
234  //total energy
235  meanEtot = layerEtot[i1] /nbEvents;
236  meanEtot2 = layerEtot2[i1]/nbEvents;
237  varianceEtot = meanEtot2 - meanEtot*meanEtot;
238  rmsEtot = 0.;
239  if (varianceEtot > 0.) rmsEtot = std::sqrt(varianceEtot);
240  resEtot = 100*rmsEtot/meanEtot;
241  histoManager->FillHisto(4, i1+0.5, meanEtot);
242 
243  //print
244  //
245  G4cout
246  << "\n layer " << i1 << ": "
247  << std::setprecision(5)
248  << std::setw(6) << G4BestUnit(meanEvis,"Energy") << " +- "
249  << std::setprecision(4)
250  << std::setw(5) << G4BestUnit( rmsEvis,"Energy") << " ("
251  << std::setprecision(2)
252  << std::setw(3) << resEvis << " %)"
253  << " "
254  << std::setprecision(5)
255  << std::setw(6) << G4BestUnit(meanEtot,"Energy") << " +- "
256  << std::setprecision(4)
257  << std::setw(5) << G4BestUnit( rmsEtot,"Energy") << " ("
258  << std::setprecision(2)
259  << std::setw(3) << resEtot << " %)";
260  }
261  G4cout << G4endl;
262 
263  //calorimeter: visible energy
264  meanEvis = calorEvis /nbEvents;
265  meanEvis2 = calorEvis2/nbEvents;
266  varianceEvis = meanEvis2 - meanEvis*meanEvis;
267  rmsEvis = 0.;
268  if (varianceEvis > 0.) rmsEvis = std::sqrt(varianceEvis);
269  resEvis = 100*rmsEvis/meanEvis;
270 
271  //calorimeter: total energy
272  meanEtot = calorEtot /nbEvents;
273  meanEtot2 = calorEtot2/nbEvents;
274  varianceEtot = meanEtot2 - meanEtot*meanEtot;
275  rmsEtot = 0.;
276  if (varianceEtot > 0.) rmsEtot = std::sqrt(varianceEtot);
277  resEtot = 100*rmsEtot/meanEtot;
278 
279  //print
280  //
281  G4cout
282  << "\n total calor : "
283  << std::setprecision(5)
284  << std::setw(6) << G4BestUnit(meanEvis,"Energy") << " +- "
285  << std::setprecision(4)
286  << std::setw(5) << G4BestUnit( rmsEvis,"Energy") << " ("
287  << std::setprecision(2)
288  << std::setw(3) << resEvis << " %)"
289  << " "
290  << std::setprecision(5)
291  << std::setw(6) << G4BestUnit(meanEtot,"Energy") << " +- "
292  << std::setprecision(4)
293  << std::setw(5) << G4BestUnit( rmsEtot,"Energy") << " ("
294  << std::setprecision(2)
295  << std::setw(3) << resEtot << " %)";
296 
297  G4cout << "\n------------------------------------------------------------"
298  << G4endl;
299 
300  //leakage
301  G4double meanEleak,meanEleak2,varianceEleak,rmsEleak,ratio;
302  meanEleak = Eleak /nbEvents;
303  meanEleak2 = Eleak2/nbEvents;
304  varianceEleak = meanEleak2 - meanEleak*meanEleak;
305  rmsEleak = 0.;
306  if (varianceEleak > 0.) rmsEleak = std::sqrt(varianceEleak);
307  ratio = 100*meanEleak/energy;
308 
309  G4double forward = 100*EdLeak[0]/(nbEvents*energy);
310  G4double bakward = 100*EdLeak[1]/(nbEvents*energy);
311  G4double lateral = 100*EdLeak[2]/(nbEvents*energy);
312  //print
313  //
314  G4cout
315  << "\n Leakage : "
316  << std::setprecision(5)
317  << std::setw(6) << G4BestUnit(meanEleak,"Energy") << " +- "
318  << std::setprecision(4)
319  << std::setw(5) << G4BestUnit( rmsEleak,"Energy")
320  << "\n Eleak/Ebeam ="
321  << std::setprecision(3)
322  << std::setw(4) << ratio << " % ( forward ="
323  << std::setw(4) << forward << " %; backward ="
324  << std::setw(4) << bakward << " %; lateral ="
325  << std::setw(4) << lateral << " %)"
326  << G4endl;
327 
328  G4cout.setf(mode,std::ios::floatfield);
329  G4cout.precision(prec);
330 
331  //save histograms
332  histoManager->save();
333 
334  // show Rndm status
335  CLHEP::HepRandom::showEngineStatus();
336 }
337 
338 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
339 
340 #include <fstream>
341 
342 void RunAction::CreateFilePixels()
343 {
344  //create file and write run header
345  //
346  G4String name = histoManager->GetFileName();
347  G4String fileName = name + ".pixels.ascii";
348 
349  std::ofstream File(fileName, std::ios::out);
350 
351  G4int n1pxl = detector->GetN1Pixels();
352  G4int n2pxl = detector->GetN2Pixels();
353  G4int n1shift = detector->GetN1Shift();
354  G4int noEvents = G4RunManager::GetRunManager()->GetCurrentRun()
355  ->GetNumberOfEventToBeProcessed();
356  File << noEvents << " " << n1pxl << " " << n2pxl << " " << n1shift
357  << G4endl;
358 }
359 
360 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
361 
362 
Generated on Sat May 25 2013 14:31:59 for Geant4 by   doxygen 1.8.4