| Geant4 Cross Reference |
1 // 1
2 // *******************************************
3 // * License and Disclaimer
4 // *
5 // * The Geant4 software is copyright of th
6 // * the Geant4 Collaboration. It is provided
7 // * conditions of the Geant4 Software License
8 // * LICENSE and available at http://cern.ch/
9 // * include a list of copyright holders.
10 // *
11 // * Neither the authors of this software syst
12 // * institutes,nor the agencies providing fin
13 // * work make any representation or warran
14 // * regarding this software system or assum
15 // * use. Please see the license in the file
16 // * for the full disclaimer and the limitatio
17 // *
18 // * This code implementation is the result
19 // * technical work of the GEANT4 collaboratio
20 // * By using, copying, modifying or distri
21 // * any work based on the software) you ag
22 // * use in resulting scientific publicati
23 // * acceptance of all terms of the Geant4 Sof
24 // *******************************************
25 //
26 /// \file Run.cc
27 /// \brief Implementation of the Run class
28 //
29 //
30 //....oooOO0OOooo........oooOO0OOooo........oo
31 //....oooOO0OOooo........oooOO0OOooo........oo
32
33 #include "Run.hh"
34
35 #include "DetectorConstruction.hh"
36 #include "HistoManager.hh"
37 #include "PrimaryGeneratorAction.hh"
38
39 #include "G4ParticleDefinition.hh"
40 #include "G4SystemOfUnits.hh"
41 #include "G4Track.hh"
42 #include "G4UnitsTable.hh"
43
44 #include <iomanip>
45
46 //....oooOO0OOooo........oooOO0OOooo........oo
47
48 Run::Run(DetectorConstruction* det) : fDetecto
49 {
50 // initialize energy deposited per absorber
51 //
52 for (G4int k = 0; k < kMaxAbsor; k++) {
53 fSumEAbs[k] = fSum2EAbs[k] = fSumLAbs[k] =
54 }
55
56 // initialize total energy deposited
57 //
58 fEdepTot = fEdepTot2 = 0.;
59
60 // initialize leakage
61 //
62 fEnergyLeak[0] = fEnergyLeak[1] = 0.;
63 fEleakTot = fEleakTot2 = 0.;
64
65 // initialize total energy released
66 //
67 fEtotal = fEtotal2 = 0.;
68
69 // initialize Eflow
70 //
71 G4int nbPlanes = (fDetector->GetNbOfLayers()
72 fEnergyFlow.resize(nbPlanes);
73 for (G4int k = 0; k < nbPlanes; k++) {
74 fEnergyFlow[k] = 0.;
75 }
76 }
77
78 //....oooOO0OOooo........oooOO0OOooo........oo
79
80 void Run::SetPrimary(G4ParticleDefinition* par
81 {
82 fParticle = particle;
83 fEkin = energy;
84 }
85
86 //....oooOO0OOooo........oooOO0OOooo........oo
87
88 void Run::CountProcesses(const G4VProcess* pro
89 {
90 if (process == nullptr) return;
91 G4String procName = process->GetProcessName(
92 std::map<G4String, G4int>::iterator it = fPr
93 if (it == fProcCounter.end()) {
94 fProcCounter[procName] = 1;
95 }
96 else {
97 fProcCounter[procName]++;
98 }
99 }
100
101 //....oooOO0OOooo........oooOO0OOooo........oo
102
103 void Run::SumEdepPerAbsorber(G4int kAbs, G4dou
104 {
105 // accumulate statistic with restriction
106 //
107 fSumEAbs[kAbs] += EAbs;
108 fSum2EAbs[kAbs] += EAbs * EAbs;
109 fSumLAbs[kAbs] += LAbs;
110 fSum2LAbs[kAbs] += LAbs * LAbs;
111 }
112
113 //....oooOO0OOooo........oooOO0OOooo........oo
114
115 void Run::SumEnergies(G4double edeptot, G4doub
116 {
117 fEdepTot += edeptot;
118 fEdepTot2 += edeptot * edeptot;
119
120 fEnergyLeak[0] += eleak0;
121 fEnergyLeak[1] += eleak1;
122 G4double eleaktot = eleak0 + eleak1;
123 fEleakTot += eleaktot;
124 fEleakTot2 += eleaktot * eleaktot;
125
126 G4double etotal = edeptot + eleaktot;
127 fEtotal += etotal;
128 fEtotal2 += etotal * etotal;
129 }
130
131 //....oooOO0OOooo........oooOO0OOooo........oo
132
133 void Run::SumEnergyFlow(G4int plane, G4double
134 {
135 fEnergyFlow[plane] += Eflow;
136 }
137
138 //....oooOO0OOooo........oooOO0OOooo........oo
139
140 void Run::Merge(const G4Run* run)
141 {
142 const Run* localRun = static_cast<const Run*
143
144 // pass information about primary particle
145 fParticle = localRun->fParticle;
146 fEkin = localRun->fEkin;
147
148 // accumulate sums
149 //
150 for (G4int k = 0; k < kMaxAbsor; k++) {
151 fSumEAbs[k] += localRun->fSumEAbs[k];
152 fSum2EAbs[k] += localRun->fSum2EAbs[k];
153 fSumLAbs[k] += localRun->fSumLAbs[k];
154 fSum2LAbs[k] += localRun->fSum2LAbs[k];
155 }
156
157 fEdepTot += localRun->fEdepTot;
158 fEdepTot2 += localRun->fEdepTot2;
159
160 fEnergyLeak[0] += localRun->fEnergyLeak[0];
161 fEnergyLeak[1] += localRun->fEnergyLeak[1];
162
163 fEleakTot += localRun->fEleakTot;
164 fEleakTot2 += localRun->fEleakTot2;
165
166 fEtotal += localRun->fEtotal;
167 fEtotal2 += localRun->fEtotal2;
168
169 G4int nbPlanes = (fDetector->GetNbOfLayers()
170 for (G4int k = 0; k < nbPlanes; k++) {
171 fEnergyFlow[k] += localRun->fEnergyFlow[k]
172 }
173
174 // map: processes count
175 std::map<G4String, G4int>::const_iterator it
176 for (itp = localRun->fProcCounter.begin(); i
177 G4String procName = itp->first;
178 G4int localCount = itp->second;
179 if (fProcCounter.find(procName) == fProcCo
180 fProcCounter[procName] = localCount;
181 }
182 else {
183 fProcCounter[procName] += localCount;
184 }
185 }
186
187 G4Run::Merge(run);
188 }
189
190 //....oooOO0OOooo........oooOO0OOooo........oo
191
192 void Run::EndOfRun()
193 {
194 // run condition
195 //
196 G4String Particle = fParticle->GetParticleNa
197 G4cout << "\n ---> The run is " << numberOfE
198 << G4BestUnit(fEkin, "Energy") << " t
199
200 // frequency of processes
201 //
202 G4cout << "\n Process calls frequency :" <<
203 G4int index = 0;
204 std::map<G4String, G4int>::iterator it;
205 for (it = fProcCounter.begin(); it != fProcC
206 G4String procName = it->first;
207 G4int count = it->second;
208 G4String space = " ";
209 if (++index % 3 == 0) space = "\n";
210 G4cout << " " << std::setw(22) << procName
211 }
212
213 G4cout << G4endl;
214 G4int nEvt = numberOfEvent;
215 G4double norm = G4double(nEvt);
216 if (norm > 0) norm = 1. / norm;
217 G4double qnorm = std::sqrt(norm);
218
219 // energy deposit per absorber
220 //
221 G4double beamEnergy = fEkin;
222 G4double sqbeam = std::sqrt(beamEnergy / GeV
223
224 G4double MeanEAbs, MeanEAbs2, rmsEAbs, resol
225 G4double MeanLAbs, MeanLAbs2, rmsLAbs;
226
227 std::ios::fmtflags mode = G4cout.flags();
228 G4int prec = G4cout.precision(2);
229 G4cout << "\n-------------------------------
230 G4cout << std::setw(16) << "material" << std
231 << "sqrt(E0(GeV))*rmsE/Edep" << std::
232
233 for (G4int k = 1; k <= fDetector->GetNbOfAbs
234 MeanEAbs = fSumEAbs[k] * norm;
235 MeanEAbs2 = fSum2EAbs[k] * norm;
236 rmsEAbs = std::sqrt(std::abs(MeanEAbs2 - M
237
238 resolution = 100. * sqbeam * rmsEAbs / Mea
239 rmsres = resolution * qnorm;
240
241 MeanLAbs = fSumLAbs[k] * norm;
242 MeanLAbs2 = fSum2LAbs[k] * norm;
243 rmsLAbs = std::sqrt(std::abs(MeanLAbs2 - M
244
245 // print
246 //
247 G4cout << std::setw(2) << k << std::setw(1
248 << std::setprecision(5) << std::set
249 << std::setprecision(4) << std::set
250 << resolution << " +- " << std::set
251 << std::setprecision(4) << std::set
252 << std::setprecision(3) << std::set
253 }
254
255 // total energy deposited
256 //
257 fEdepTot /= nEvt;
258 fEdepTot2 /= nEvt;
259 G4double rmsEdep = std::sqrt(std::abs(fEdepT
260
261 G4cout << "\n Total energy deposited = " <<
262 << " +- " << G4BestUnit(rmsEdep, "Ene
263
264 // Energy leakage
265 //
266 fEnergyLeak[0] /= nEvt;
267 fEnergyLeak[1] /= nEvt;
268 fEleakTot /= nEvt;
269 fEleakTot2 /= nEvt;
270 G4double rmsEleak = std::sqrt(std::abs(fElea
271
272 G4cout << " Leakage : primary = " << G4Best
273 << " secondaries = " << G4BestUnit(
274 << " ---> total = " << G4BestUnit(fE
275 << G4BestUnit(rmsEleak, "Energy") <<
276
277 // total energy released
278 //
279 fEtotal /= nEvt;
280 fEtotal2 /= nEvt;
281 G4double rmsEtotal = std::sqrt(std::abs(fEto
282
283 G4cout << " Total energy released : Edep +
284 << G4BestUnit(rmsEtotal, "Energy") <<
285 G4cout << "---------------------------------
286
287 // Energy flow
288 //
289 G4AnalysisManager* analysis = G4AnalysisMana
290 G4int Idmax = (fDetector->GetNbOfLayers()) *
291 for (G4int Id = 1; Id <= Idmax + 1; Id++) {
292 analysis->FillH1(2 * kMaxAbsor + 1, (G4dou
293 }
294
295 // normalize histograms
296 //
297 for (G4int ih = kMaxAbsor + 1; ih < 2 * kMax
298 analysis->ScaleH1(ih, norm / MeV);
299 }
300
301 // remove all contents in fProcCounter
302 fProcCounter.clear();
303
304 G4cout.setf(mode, std::ios::floatfield);
305 G4cout.precision(prec);
306 }
307
308 //....oooOO0OOooo........oooOO0OOooo........oo
309