| Geant4 Cross Reference |
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 /// \file electromagnetic/TestEm3/src/Run.cc
27 /// \brief Implementation of the Run class
28 //
29 //
30 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
31 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
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........oooOO0OOooo........oooOO0OOooo......
52
53 Run::Run(DetectorConstruction* det) : fDetector(det)
54 {
55 // initialize cumulative quantities
56 //
57 for (G4int k = 0; k < kMaxAbsor; k++) {
58 fSumEAbs[k] = fSum2EAbs[k] = fSumLAbs[k] = fSum2LAbs[k] = 0.;
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()) * (fDetector->GetNbOfAbsor()) + 2;
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........oooOO0OOooo........oooOO0OOooo......
79
80 void Run::SetPrimary(G4ParticleDefinition* particle, G4double energy)
81 {
82 fParticle = particle;
83 fEkin = energy;
84 }
85
86 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
87
88 void Run::FillPerEvent(G4int kAbs, G4double EAbs, G4double LAbs)
89 {
90 // accumulate statistic with restriction
91 //
92 if (fApplyLimit) fEnergyDeposit[kAbs].push_back(EAbs);
93 fSumEAbs[kAbs] += EAbs;
94 fSum2EAbs[kAbs] += EAbs * EAbs;
95 fSumLAbs[kAbs] += LAbs;
96 fSum2LAbs[kAbs] += LAbs * LAbs;
97 }
98
99 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
100
101 void Run::SumEnergies(G4double edeptot, G4double eleak)
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........oooOO0OOooo........oooOO0OOooo......
114
115 void Run::SumEnergyFlow(G4int plane, G4double Eflow)
116 {
117 fEnergyFlow[plane] += Eflow;
118 }
119
120 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
121
122 void Run::SumLateralEleak(G4int cell, G4double Eflow)
123 {
124 fLateralEleak[cell] += Eflow;
125 }
126
127 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
128
129 void Run::AddChargedStep()
130 {
131 fChargedStep += 1.0;
132 }
133
134 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
135
136 void Run::AddNeutralStep()
137 {
138 fNeutralStep += 1.0;
139 }
140
141 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
142
143 void Run::AddSecondaryTrack(const G4Track* track)
144 {
145 const G4ParticleDefinition* d = track->GetDefinition();
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........oooOO0OOooo........oooOO0OOooo......
158
159 void Run::Merge(const G4Run* run)
160 {
161 const Run* localRun = static_cast<const Run*>(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()) * (fDetector->GetNbOfAbsor()) + 2;
186 for (G4int k = 0; k < nbPlanes; k++) {
187 fEnergyFlow[k] += localRun->fEnergyFlow[k];
188 fLateralEleak[k] += localRun->fLateralEleak[k];
189 }
190
191 for (G4int k=0; k<kMaxAbsor; k++) {
192 fEnergyDeposit[k].insert(fEnergyDeposit[k].end(),
193 localRun->fEnergyDeposit[k].begin(), localRun->fEnergyDeposit[k].end());
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........oooOO0OOooo........oooOO0OOooo......
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, resolution, rmsres;
232 G4double MeanLAbs, MeanLAbs2, rmsLAbs;
233
234 std::ios::fmtflags mode = G4cout.flags();
235 G4int prec = G4cout.precision(2);
236 G4cout << "\n------------------------------------------------------------\n";
237 G4cout << std::setw(14) << "material" << std::setw(17) << "Edep RMS" << std::setw(33)
238 << "sqrt(E0(GeV))*rmsE/Emean" << std::setw(23) << "total tracklen \n \n";
239
240 for (G4int k = 1; k <= fDetector->GetNbOfAbsor(); k++) {
241 MeanEAbs = fSumEAbs[k] * norm;
242 MeanEAbs2 = fSum2EAbs[k] * norm;
243 rmsEAbs = std::sqrt(std::abs(MeanEAbs2 - MeanEAbs * MeanEAbs));
244 // G4cout << "k= " << k << " RMS= " << rmsEAbs
245 // << " fApplyLimit: " << fApplyLimit << G4endl;
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 - MeanEAbs * MeanEAbs));
265 }
266
267 resolution = (MeanEAbs > 0.) ? 100. * sqbeam * rmsEAbs / MeanEAbs : 0.0;
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 - MeanLAbs * MeanLAbs));
277
278 // print
279 //
280 G4cout << std::setw(14) << fDetector->GetAbsorMaterial(k)->GetName() << ": "
281 << std::setprecision(5) << std::setw(6) << G4BestUnit(MeanEAbs, "Energy") << " : "
282 << std::setprecision(4) << std::setw(5) << G4BestUnit(rmsEAbs, "Energy") << std::setw(10)
283 << resolution << " +- " << std::setw(5) << rmsres << " %" << std::setprecision(3)
284 << std::setw(10) << G4BestUnit(MeanLAbs, "Length") << " +- " << std::setw(4)
285 << G4BestUnit(rmsLAbs, "Length") << G4endl;
286 }
287
288 // total energy deposited
289 //
290 fEdepTot *= norm;
291 fEdepTot2 *= norm;
292 G4double rmsEdep = std::sqrt(std::abs(fEdepTot2 - fEdepTot * fEdepTot));
293
294 G4cout << "\n Total energy deposited = " << std::setprecision(4) << G4BestUnit(fEdepTot, "Energy")
295 << " +- " << G4BestUnit(rmsEdep, "Energy") << G4endl;
296
297 // Energy leakage
298 //
299 fEleakTot *= norm;
300 fEleakTot2 *= norm;
301 G4double rmsEleak = std::sqrt(std::abs(fEleakTot2 - fEleakTot * fEleakTot));
302
303 G4cout << " Energy leakage = " << G4BestUnit(fEleakTot, "Energy") << " +- "
304 << G4BestUnit(rmsEleak, "Energy") << G4endl;
305
306 // total energy
307 //
308 fEtotal *= norm;
309 fEtotal2 *= norm;
310 G4double rmsEtotal = std::sqrt(std::abs(fEtotal2 - fEtotal * fEtotal));
311
312 G4cout << " Total energy : Edep + Eleak = " << G4BestUnit(fEtotal, "Energy") << " +- "
313 << G4BestUnit(rmsEtotal, "Energy") << G4endl;
314
315 G4cout << "------------------------------------------------------------\n";
316
317 G4cout << " Beam particle " << fParticle->GetParticleName()
318 << " E = " << G4BestUnit(beamEnergy, "Energy") << G4endl;
319 G4cout << " Mean number of gamma " << (G4double)fN_gamma * norm << G4endl;
320 G4cout << " Mean number of e- " << (G4double)fN_elec * norm << G4endl;
321 G4cout << " Mean number of e+ " << (G4double)fN_pos * norm << G4endl;
322 G4cout << std::setprecision(6) << " Mean number of charged steps " << fChargedStep << G4endl;
323 G4cout << " Mean number of neutral steps " << fNeutralStep << G4endl;
324 G4cout << "------------------------------------------------------------\n";
325
326 // Energy flow
327 //
328 G4AnalysisManager* analysis = G4AnalysisManager::Instance();
329 G4int Idmax = (fDetector->GetNbOfLayers()) * (fDetector->GetNbOfAbsor());
330 for (G4int Id = 1; Id <= Idmax + 1; Id++) {
331 analysis->FillH1(2 * kMaxAbsor + 1, (G4double)Id, fEnergyFlow[Id]);
332 analysis->FillH1(2 * kMaxAbsor + 2, (G4double)Id, fLateralEleak[Id]);
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] - fEnergyFlow[Id + 1] - fLateralEleak[Id]);
346 }
347
348 G4cout << std::setprecision(3) << "\n Energy deposition from Energy flow balance : \n"
349 << std::setw(10) << " material \t Total Edep \n \n";
350 G4cout.precision(6);
351
352 for (G4int k = 1; k <= nbOfAbsor; k++) {
353 EdepTot[k] *= norm;
354 G4cout << std::setw(10) << fDetector->GetAbsorMaterial(k)->GetName() << ":"
355 << "\t " << G4BestUnit(EdepTot[k], "Energy") << "\n";
356 }
357
358 G4cout << "\n------------------------------------------------------------\n" << G4endl;
359
360 // Acceptance
361 EmAcceptance acc;
362 G4bool isStarted = false;
363 for (G4int j = 1; j <= fDetector->GetNbOfAbsor(); j++) {
364 if (fLimittrue[j] < DBL_MAX) {
365 if (!isStarted) {
366 acc.BeginOfAcceptance("Sampling Calorimeter", nEvt);
367 isStarted = true;
368 }
369 MeanEAbs = fSumEAbs[j];
370 rmsEAbs = fSum2EAbs[j];
371 G4String mat = fDetector->GetAbsorMaterial(j)->GetName();
372 acc.EmAcceptanceGauss("Edep" + mat, nEvt, MeanEAbs, fEdeptrue[j], fRmstrue[j], fLimittrue[j]);
373 acc.EmAcceptanceGauss("Erms" + mat, nEvt, rmsEAbs, fRmstrue[j], fRmstrue[j],
374 2.0 * fLimittrue[j]);
375 }
376 }
377 if (isStarted) acc.EndOfAcceptance();
378
379 // normalize histograms
380 //
381 for (G4int ih = kMaxAbsor + 1; ih < kMaxHisto - 2; ih++) {
382 analysis->ScaleH1(ih, norm / MeV);
383 }
384
385 G4cout.setf(mode, std::ios::floatfield);
386 G4cout.precision(prec);
387 }
388
389 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
390
391 void Run::SetEdepAndRMS(G4int i, G4double edep, G4double rms, G4double lim)
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........oooOO0OOooo........oooOO0OOooo......
401
402 void Run::SetApplyLimit(G4bool val)
403 {
404 fApplyLimit = val;
405 }
406
407 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
408