| 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 electromagnetic/TestEm3/src/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 "EmAcceptance.hh"
37 #include "HistoManager.hh"
38 #include "PrimaryGeneratorAction.hh"
39
40 #include "G4Electron.hh"
41 #include "G4Gamma.hh"
42 #include "G4ParticleDefinition.hh"
43 #include "G4ParticleTable.hh"
44 #include "G4Positron.hh"
45 #include "G4SystemOfUnits.hh"
46 #include "G4Track.hh"
47 #include "G4UnitsTable.hh"
48
49 #include <iomanip>
50
51 //....oooOO0OOooo........oooOO0OOooo........oo
52
53 Run::Run(DetectorConstruction* det) : fDetecto
54 {
55 // initialize cumulative quantities
56 //
57 for (G4int k = 0; k < kMaxAbsor; k++) {
58 fSumEAbs[k] = fSum2EAbs[k] = fSumLAbs[k] =
59 fEnergyDeposit[k].clear();
60 fEdeptrue[k] = fRmstrue[k] = 1.;
61 fLimittrue[k] = DBL_MAX;
62 }
63
64 fEdepTot = fEdepTot2 = 0.;
65 fEleakTot = fEleakTot2 = 0.;
66 fEtotal = fEtotal2 = 0.;
67
68 // initialize Eflow
69 //
70 G4int nbPlanes = (fDetector->GetNbOfLayers()
71 fEnergyFlow.resize(nbPlanes);
72 fLateralEleak.resize(nbPlanes);
73 for (G4int k = 0; k < nbPlanes; k++) {
74 fEnergyFlow[k] = fLateralEleak[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::FillPerEvent(G4int kAbs, G4double EA
89 {
90 // accumulate statistic with restriction
91 //
92 if (fApplyLimit) fEnergyDeposit[kAbs].push_b
93 fSumEAbs[kAbs] += EAbs;
94 fSum2EAbs[kAbs] += EAbs * EAbs;
95 fSumLAbs[kAbs] += LAbs;
96 fSum2LAbs[kAbs] += LAbs * LAbs;
97 }
98
99 //....oooOO0OOooo........oooOO0OOooo........oo
100
101 void Run::SumEnergies(G4double edeptot, G4doub
102 {
103 fEdepTot += edeptot;
104 fEdepTot2 += edeptot * edeptot;
105 fEleakTot += eleak;
106 fEleakTot2 += eleak * eleak;
107
108 G4double etotal = edeptot + eleak;
109 fEtotal += etotal;
110 fEtotal2 += etotal * etotal;
111 }
112
113 //....oooOO0OOooo........oooOO0OOooo........oo
114
115 void Run::SumEnergyFlow(G4int plane, G4double
116 {
117 fEnergyFlow[plane] += Eflow;
118 }
119
120 //....oooOO0OOooo........oooOO0OOooo........oo
121
122 void Run::SumLateralEleak(G4int cell, G4double
123 {
124 fLateralEleak[cell] += Eflow;
125 }
126
127 //....oooOO0OOooo........oooOO0OOooo........oo
128
129 void Run::AddChargedStep()
130 {
131 fChargedStep += 1.0;
132 }
133
134 //....oooOO0OOooo........oooOO0OOooo........oo
135
136 void Run::AddNeutralStep()
137 {
138 fNeutralStep += 1.0;
139 }
140
141 //....oooOO0OOooo........oooOO0OOooo........oo
142
143 void Run::AddSecondaryTrack(const G4Track* tra
144 {
145 const G4ParticleDefinition* d = track->GetDe
146 if (d == G4Gamma::Gamma()) {
147 ++fN_gamma;
148 }
149 else if (d == G4Electron::Electron()) {
150 ++fN_elec;
151 }
152 else if (d == G4Positron::Positron()) {
153 ++fN_pos;
154 }
155 }
156
157 //....oooOO0OOooo........oooOO0OOooo........oo
158
159 void Run::Merge(const G4Run* run)
160 {
161 const Run* localRun = static_cast<const Run*
162
163 // pass information about primary particle
164 fParticle = localRun->fParticle;
165 fEkin = localRun->fEkin;
166
167 // accumulate sums
168 //
169 for (G4int k = 0; k < kMaxAbsor; k++) {
170 fSumEAbs[k] += localRun->fSumEAbs[k];
171 fSum2EAbs[k] += localRun->fSum2EAbs[k];
172 fSumLAbs[k] += localRun->fSumLAbs[k];
173 fSum2LAbs[k] += localRun->fSum2LAbs[k];
174 }
175
176 fEdepTot += localRun->fEdepTot;
177 fEdepTot2 += localRun->fEdepTot2;
178
179 fEleakTot += localRun->fEleakTot;
180 fEleakTot2 += localRun->fEleakTot2;
181
182 fEtotal += localRun->fEtotal;
183 fEtotal2 += localRun->fEtotal2;
184
185 G4int nbPlanes = (fDetector->GetNbOfLayers()
186 for (G4int k = 0; k < nbPlanes; k++) {
187 fEnergyFlow[k] += localRun->fEnergyFlow[k]
188 fLateralEleak[k] += localRun->fLateralElea
189 }
190
191 for (G4int k=0; k<kMaxAbsor; k++) {
192 fEnergyDeposit[k].insert(fEnergyDeposit[k]
193 localRun->fEnergyDeposit[k].begin(), local
194 }
195
196 fChargedStep += localRun->fChargedStep;
197 fNeutralStep += localRun->fNeutralStep;
198
199 fN_gamma += localRun->fN_gamma;
200 fN_elec += localRun->fN_elec;
201 fN_pos += localRun->fN_pos;
202
203 fApplyLimit = localRun->fApplyLimit;
204
205 for (G4int k = 0; k < kMaxAbsor; k++) {
206 fEdeptrue[k] = localRun->fEdeptrue[k];
207 fRmstrue[k] = localRun->fRmstrue[k];
208 fLimittrue[k] = localRun->fLimittrue[k];
209 }
210
211 G4Run::Merge(run);
212 }
213
214 //....oooOO0OOooo........oooOO0OOooo........oo
215
216 void Run::EndOfRun()
217 {
218 G4int nEvt = numberOfEvent;
219 G4double norm = G4double(nEvt);
220 if (norm > 0) norm = 1. / norm;
221 G4double qnorm = std::sqrt(norm);
222
223 fChargedStep *= norm;
224 fNeutralStep *= norm;
225
226 // compute and print statistic
227 //
228 G4double beamEnergy = fEkin;
229 G4double sqbeam = std::sqrt(beamEnergy / GeV
230
231 G4double MeanEAbs, MeanEAbs2, rmsEAbs, resol
232 G4double MeanLAbs, MeanLAbs2, rmsLAbs;
233
234 std::ios::fmtflags mode = G4cout.flags();
235 G4int prec = G4cout.precision(2);
236 G4cout << "\n-------------------------------
237 G4cout << std::setw(14) << "material" << std
238 << "sqrt(E0(GeV))*rmsE/Emean" << std:
239
240 for (G4int k = 1; k <= fDetector->GetNbOfAbs
241 MeanEAbs = fSumEAbs[k] * norm;
242 MeanEAbs2 = fSum2EAbs[k] * norm;
243 rmsEAbs = std::sqrt(std::abs(MeanEAbs2 - M
244 // G4cout << "k= " << k << " RMS= " << r
245 // << " fApplyLimit: " << fApplyLimi
246 if (fApplyLimit) {
247 G4int nn = 0;
248 G4double sume = 0.0;
249 G4double sume2 = 0.0;
250 // compute trancated means
251 G4double lim = rmsEAbs * 2.5;
252 for (G4int i = 0; i < nEvt; i++) {
253 G4double e = (fEnergyDeposit[k])[i];
254 if (std::abs(e - MeanEAbs) < lim) {
255 sume += e;
256 sume2 += e * e;
257 nn++;
258 }
259 }
260 G4double norm1 = G4double(nn);
261 if (norm1 > 0.0) norm1 = 1.0 / norm1;
262 MeanEAbs = sume * norm1;
263 MeanEAbs2 = sume2 * norm1;
264 rmsEAbs = std::sqrt(std::abs(MeanEAbs2 -
265 }
266
267 resolution = (MeanEAbs > 0.) ? 100. * sqbe
268 rmsres = resolution * qnorm;
269
270 // Save mean and RMS
271 fSumEAbs[k] = MeanEAbs;
272 fSum2EAbs[k] = rmsEAbs;
273
274 MeanLAbs = fSumLAbs[k] * norm;
275 MeanLAbs2 = fSum2LAbs[k] * norm;
276 rmsLAbs = std::sqrt(std::abs(MeanLAbs2 - M
277
278 // print
279 //
280 G4cout << std::setw(14) << fDetector->GetA
281 << std::setprecision(5) << std::set
282 << std::setprecision(4) << std::set
283 << resolution << " +- " << std::set
284 << std::setw(10) << G4BestUnit(Mean
285 << G4BestUnit(rmsLAbs, "Length") <<
286 }
287
288 // total energy deposited
289 //
290 fEdepTot *= norm;
291 fEdepTot2 *= norm;
292 G4double rmsEdep = std::sqrt(std::abs(fEdepT
293
294 G4cout << "\n Total energy deposited = " <<
295 << " +- " << G4BestUnit(rmsEdep, "Ene
296
297 // Energy leakage
298 //
299 fEleakTot *= norm;
300 fEleakTot2 *= norm;
301 G4double rmsEleak = std::sqrt(std::abs(fElea
302
303 G4cout << " Energy leakage = " << G4BestUnit
304 << G4BestUnit(rmsEleak, "Energy") <<
305
306 // total energy
307 //
308 fEtotal *= norm;
309 fEtotal2 *= norm;
310 G4double rmsEtotal = std::sqrt(std::abs(fEto
311
312 G4cout << " Total energy : Edep + Eleak = "
313 << G4BestUnit(rmsEtotal, "Energy") <<
314
315 G4cout << "---------------------------------
316
317 G4cout << " Beam particle " << fParticle->Ge
318 << " E = " << G4BestUnit(beamEnergy,
319 G4cout << " Mean number of gamma " <<
320 G4cout << " Mean number of e- " <<
321 G4cout << " Mean number of e+ " <<
322 G4cout << std::setprecision(6) << " Mean num
323 G4cout << " Mean number of neutral steps "
324 G4cout << "---------------------------------
325
326 // Energy flow
327 //
328 G4AnalysisManager* analysis = G4AnalysisMana
329 G4int Idmax = (fDetector->GetNbOfLayers()) *
330 for (G4int Id = 1; Id <= Idmax + 1; Id++) {
331 analysis->FillH1(2 * kMaxAbsor + 1, (G4dou
332 analysis->FillH1(2 * kMaxAbsor + 2, (G4dou
333 }
334
335 // Energy deposit from energy flow balance
336 //
337 G4double EdepTot[kMaxAbsor];
338 for (G4int k = 0; k < kMaxAbsor; k++)
339 EdepTot[k] = 0.;
340
341 G4int nbOfAbsor = fDetector->GetNbOfAbsor();
342 for (G4int Id = 1; Id <= Idmax; Id++) {
343 G4int iAbsor = Id % nbOfAbsor;
344 if (iAbsor == 0) iAbsor = nbOfAbsor;
345 EdepTot[iAbsor] += (fEnergyFlow[Id] - fEne
346 }
347
348 G4cout << std::setprecision(3) << "\n Energy
349 << std::setw(10) << " material \t To
350 G4cout.precision(6);
351
352 for (G4int k = 1; k <= nbOfAbsor; k++) {
353 EdepTot[k] *= norm;
354 G4cout << std::setw(10) << fDetector->GetA
355 << "\t " << G4BestUnit(EdepTot[k],
356 }
357
358 G4cout << "\n-------------------------------
359
360 // Acceptance
361 EmAcceptance acc;
362 G4bool isStarted = false;
363 for (G4int j = 1; j <= fDetector->GetNbOfAbs
364 if (fLimittrue[j] < DBL_MAX) {
365 if (!isStarted) {
366 acc.BeginOfAcceptance("Sampling Calori
367 isStarted = true;
368 }
369 MeanEAbs = fSumEAbs[j];
370 rmsEAbs = fSum2EAbs[j];
371 G4String mat = fDetector->GetAbsorMateri
372 acc.EmAcceptanceGauss("Edep" + mat, nEvt
373 acc.EmAcceptanceGauss("Erms" + mat, nEvt
374 2.0 * fLimittrue[j
375 }
376 }
377 if (isStarted) acc.EndOfAcceptance();
378
379 // normalize histograms
380 //
381 for (G4int ih = kMaxAbsor + 1; ih < kMaxHist
382 analysis->ScaleH1(ih, norm / MeV);
383 }
384
385 G4cout.setf(mode, std::ios::floatfield);
386 G4cout.precision(prec);
387 }
388
389 //....oooOO0OOooo........oooOO0OOooo........oo
390
391 void Run::SetEdepAndRMS(G4int i, G4double edep
392 {
393 if (i >= 0 && i < kMaxAbsor) {
394 fEdeptrue[i] = edep;
395 fRmstrue[i] = rms;
396 fLimittrue[i] = lim;
397 }
398 }
399
400 //....oooOO0OOooo........oooOO0OOooo........oo
401
402 void Run::SetApplyLimit(G4bool val)
403 {
404 fApplyLimit = val;
405 }
406
407 //....oooOO0OOooo........oooOO0OOooo........oo
408