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