Geant4  10.01.p02
Run.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 //
28 //
29 // $Id: Run.cc 71376 2013-06-14 07:44:50Z maire $
30 //
31 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
32 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
33 
34 #include "Run.hh"
35 #include "DetectorConstruction.hh"
36 #include "HistoManager.hh"
37 
38 #include "G4ParticleDefinition.hh"
39 #include "G4SystemOfUnits.hh"
40 #include "G4UnitsTable.hh"
41 
42 #include <iomanip>
43 
44 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
45 
47 : G4Run(),
48  fDetector(det),
49  fParticle(0), fEkin(0.),
50  nbOfModules(0), nbOfLayers(0), kLayerMax(0),
51  EtotCalor(0.), Etot2Calor(0.), EvisCalor(0.), Evis2Calor(0.),
52  Eleak(0.), Eleak2(0.)
53 {
57 
58  //initialize vectors
59  //
60  EtotLayer.resize(kLayerMax); Etot2Layer.resize(kLayerMax);
61  EvisLayer.resize(kLayerMax); Evis2Layer.resize(kLayerMax);
62  for (G4int k=0; k<kLayerMax; k++) {
63  EtotLayer[k] = Etot2Layer[k] = EvisLayer[k] = Evis2Layer[k] = 0.0;
64  }
65 
67  EdLeak[0] = EdLeak[1] = EdLeak[2] = 0.;
68 }
69 
70 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
71 
73 { }
74 
75 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
76 
78 {
79  fParticle = particle;
80  fEkin = energy;
81 }
82 
83 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
84 
85 void Run::SumEvents_1(G4int layer, G4double Etot, G4double Evis)
86 {
87  //accumulate statistic per layer
88  //
89  EtotLayer[layer] += Etot; Etot2Layer[layer] += Etot*Etot;
90  EvisLayer[layer] += Evis; Evis2Layer[layer] += Evis*Evis;
91 }
92 
93 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
94 
95 void Run::SumEvents_2(G4double etot, G4double evis, G4double eleak)
96 {
97  //accumulate statistic for full calorimeter
98  //
99  EtotCalor += etot; Etot2Calor += etot*etot;
100  EvisCalor += evis; Evis2Calor += evis*evis;
101  Eleak += eleak; Eleak2 += eleak*eleak;
102 }
103 
104 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
105 
107 {
108  //forward, backward, lateral leakage
109  //
110  EdLeak[icase] += energy;
111 }
112 
113 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
114 
115 void Run::Merge(const G4Run* run)
116 {
117  const Run* localRun = static_cast<const Run*>(run);
118 
119  // pass information about primary particle
120  fParticle = localRun->fParticle;
121  fEkin = localRun->fEkin;
122 
123  // accumulate sums
124  //
125  for (G4int k=0; k<kLayerMax; k++) {
126  EtotLayer[k] += localRun->EtotLayer[k];
127  Etot2Layer[k] += localRun->Etot2Layer[k];
128  EvisLayer[k] += localRun->EvisLayer[k];
129  Evis2Layer[k] += localRun->Evis2Layer[k];
130  }
131 
132  EtotCalor += localRun->EtotCalor;
133  Etot2Calor += localRun->Etot2Calor;
134  EvisCalor += localRun->EvisCalor;
135  Evis2Calor += localRun->Evis2Calor;
136  Eleak += localRun->Eleak;
137  Eleak2 += localRun->Eleak2;
138  EdLeak[0] += localRun->EdLeak[0];
139  EdLeak[1] += localRun->EdLeak[1];
140  EdLeak[2] += localRun->EdLeak[2];
141 
142  G4Run::Merge(run);
143 }
144 
145 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
146 
148 {
149  //calorimeter
150  //
152 
153  //run conditions
154  //
155  G4String partName = fParticle->GetParticleName();
156  G4int nbEvents = numberOfEvent;
157 
158  G4int prec = G4cout.precision(3);
159 
160  G4cout << " The run was " << nbEvents << " " << partName << " of "
161  << G4BestUnit(fEkin,"Energy") << " through the calorimeter" << G4endl;
162 
163  G4cout << "------------------------------------------------------------"
164  << G4endl;
165 
166  //if no events, return
167  //
168  if (nbEvents == 0) return;
169 
170  //compute and print statistic
171  //
172  std::ios::fmtflags mode = G4cout.flags();
173 
174  // energy in layers
175  //
176  G4cout.precision(prec);
177  G4cout << "\n "
178  << "total Energy (rms/mean) "
179  << "visible Energy (rms/mean)" << G4endl;
180 
181  G4AnalysisManager* analysisManager = G4AnalysisManager::Instance();
182 
183  G4double meanEtot,meanEtot2,varianceEtot,rmsEtot,resEtot;
184  G4double meanEvis,meanEvis2,varianceEvis,rmsEvis,resEvis;
185 
186  for (G4int i1=1; i1<kLayerMax; i1++) {
187  //total energy
188  meanEtot = EtotLayer[i1] /nbEvents;
189  meanEtot2 = Etot2Layer[i1]/nbEvents;
190  varianceEtot = meanEtot2 - meanEtot*meanEtot;
191  resEtot = rmsEtot = 0.;
192  if (varianceEtot > 0.) rmsEtot = std::sqrt(varianceEtot);
193  if (meanEtot > 0.) resEtot = 100*rmsEtot/meanEtot;
194  analysisManager->FillH1(3, i1+0.5, meanEtot);
195 
196  //visible energy
197  meanEvis = EvisLayer[i1] /nbEvents;
198  meanEvis2 = Evis2Layer[i1]/nbEvents;
199  varianceEvis = meanEvis2 - meanEvis*meanEvis;
200  resEvis = rmsEvis = 0.;
201  if (varianceEvis > 0.) rmsEvis = std::sqrt(varianceEvis);
202  if (meanEvis > 0.) resEvis = 100*rmsEvis/meanEvis;
203  analysisManager->FillH1(4, i1+0.5, meanEvis);
204 
205  //print
206  //
207  G4cout
208  << "\n layer " << i1 << ": "
209  << std::setprecision(5)
210  << std::setw(6) << G4BestUnit(meanEtot,"Energy") << " +- "
211  << std::setprecision(4)
212  << std::setw(5) << G4BestUnit( rmsEtot,"Energy") << " ("
213  << std::setprecision(2)
214  << std::setw(3) << resEtot << " %)"
215  << " "
216  << std::setprecision(5)
217  << std::setw(6) << G4BestUnit(meanEvis,"Energy") << " +- "
218  << std::setprecision(4)
219  << std::setw(5) << G4BestUnit( rmsEvis,"Energy") << " ("
220  << std::setprecision(2)
221  << std::setw(3) << resEvis << " %)";
222  }
223  G4cout << G4endl;
224 
225  //calorimeter: total energy
226  meanEtot = EtotCalor /nbEvents;
227  meanEtot2 = Etot2Calor/nbEvents;
228  varianceEtot = meanEtot2 - meanEtot*meanEtot;
229  resEtot = rmsEtot = 0.;
230  if (varianceEtot > 0.) rmsEtot = std::sqrt(varianceEtot);
231  if (meanEtot > 0.) resEtot = 100*rmsEtot/meanEtot;
232 
233  //calorimeter: visible energy
234  meanEvis = EvisCalor /nbEvents;
235  meanEvis2 = Evis2Calor/nbEvents;
236  varianceEvis = meanEvis2 - meanEvis*meanEvis;
237  resEvis = rmsEvis = 0.;
238  if (varianceEvis > 0.) rmsEvis = std::sqrt(varianceEvis);
239  if (meanEvis > 0.) resEvis = 100*rmsEvis/meanEvis;
240 
241  //print
242  //
243  G4cout
244  << "\n total calor : "
245  << std::setprecision(5)
246  << std::setw(6) << G4BestUnit(meanEtot,"Energy") << " +- "
247  << std::setprecision(4)
248  << std::setw(5) << G4BestUnit( rmsEtot,"Energy") << " ("
249  << std::setprecision(2)
250  << std::setw(3) << resEtot << " %)"
251  << " "
252  << std::setprecision(5)
253  << std::setw(6) << G4BestUnit(meanEvis,"Energy") << " +- "
254  << std::setprecision(4)
255  << std::setw(5) << G4BestUnit( rmsEvis,"Energy") << " ("
256  << std::setprecision(2)
257  << std::setw(3) << resEvis << " %)";
258 
259  G4cout << "\n------------------------------------------------------------"
260  << G4endl;
261 
262  //leakage
263  G4double meanEleak,meanEleak2,varianceEleak,rmsEleak,ratio;
264  meanEleak = Eleak /nbEvents;
265  meanEleak2 = Eleak2/nbEvents;
266  varianceEleak = meanEleak2 - meanEleak*meanEleak;
267  rmsEleak = 0.;
268  if (varianceEleak > 0.) rmsEleak = std::sqrt(varianceEleak);
269  ratio = 100*meanEleak/fEkin;
270 
271  G4double forward = 100*EdLeak[0]/(nbEvents*fEkin);
272  G4double bakward = 100*EdLeak[1]/(nbEvents*fEkin);
273  G4double lateral = 100*EdLeak[2]/(nbEvents*fEkin);
274 
275  //print
276  //
277  G4cout
278  << "\n Leakage : "
279  << std::setprecision(5)
280  << std::setw(6) << G4BestUnit(meanEleak,"Energy") << " +- "
281  << std::setprecision(4)
282  << std::setw(5) << G4BestUnit( rmsEleak,"Energy")
283  << "\n Eleak/Ebeam ="
284  << std::setprecision(3)
285  << std::setw(4) << ratio << " % ( forward ="
286  << std::setw(4) << forward << " %; backward ="
287  << std::setw(4) << bakward << " %; lateral ="
288  << std::setw(4) << lateral << " %)"
289  << G4endl;
290 
291  G4cout.setf(mode,std::ios::floatfield);
292  G4cout.precision(prec);
293 
294  //normalize histograms
295  G4double factor = 1./nbEvents;
296  analysisManager->ScaleH1(5,factor);
297 }
298 
299 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
virtual void Merge(const G4Run *)
Definition: G4Run.cc:54
G4int numberOfEvent
Definition: G4Run.hh:59
G4double fEkin
Definition: Run.hh:65
void SumEvents_2(G4double, G4double, G4double)
Definition: Run.cc:95
Run()
Definition: Run.cc:43
#define G4BestUnit(a, b)
#define G4_USE_G4BESTUNIT_FOR_VERBOSE 1
std::vector< G4double > Evis2Layer
Definition: Run.hh:69
int G4int
Definition: G4Types.hh:78
const G4String & GetParticleName() const
static const double prec
Definition: RanecuEngine.cc:58
G4GLOB_DLL std::ostream G4cout
G4double Eleak
Definition: Run.hh:73
G4double EdLeak[3]
Definition: Run.hh:74
G4double Etot2Calor
Definition: Run.hh:71
G4double Evis2Calor
Definition: Run.hh:72
Definition: G4Run.hh:46
G4double EvisCalor
Definition: Run.hh:72
std::vector< G4double > Etot2Layer
Definition: Run.hh:68
ExG4HbookAnalysisManager G4AnalysisManager
Definition: g4hbook_defs.hh:61
G4double Eleak2
Definition: Run.hh:73
void DetailedLeakage(G4int, G4double)
Definition: Run.cc:106
static const G4double factor
void SumEvents_1(G4int, G4double, G4double)
Definition: Run.cc:85
G4double energy(const ThreeVector &p, const G4double m)
std::vector< G4double > EvisLayer
Definition: Run.hh:69
G4double EtotCalor
Definition: Run.hh:71
G4int nbOfLayers
Definition: Run.hh:67
DetectorConstruction * fDetector
Definition: Run.hh:63
Detector construction class to demonstrate various ways of placement.
void EndOfRun()
Definition: Run.cc:147
#define G4endl
Definition: G4ios.hh:61
G4ParticleDefinition * fParticle
Definition: Run.hh:64
G4int nbOfModules
Definition: Run.hh:67
void SetPrimary(G4ParticleDefinition *particle, G4double energy)
Definition: Run.cc:77
double G4double
Definition: G4Types.hh:76
Definition: Run.hh:46
std::vector< G4double > EtotLayer
Definition: Run.hh:68
virtual void Merge(const G4Run *)
Definition: Run.cc:115
~Run()
Definition: Run.cc:72
G4int kLayerMax
Definition: Run.hh:67