| 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 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 "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........oooOO0OOooo........oooOO0OOooo......
47
48 Run::Run(DetectorConstruction* det) : fDetector(det)
49 {
50 // initialize energy deposited per absorber
51 //
52 for (G4int k = 0; k < kMaxAbsor; k++) {
53 fSumEAbs[k] = fSum2EAbs[k] = fSumLAbs[k] = fSum2LAbs[k] = 0.;
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()) * (fDetector->GetNbOfAbsor()) + 2;
72 fEnergyFlow.resize(nbPlanes);
73 for (G4int k = 0; k < nbPlanes; k++) {
74 fEnergyFlow[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::CountProcesses(const G4VProcess* process)
89 {
90 if (process == nullptr) return;
91 G4String procName = process->GetProcessName();
92 std::map<G4String, G4int>::iterator it = fProcCounter.find(procName);
93 if (it == fProcCounter.end()) {
94 fProcCounter[procName] = 1;
95 }
96 else {
97 fProcCounter[procName]++;
98 }
99 }
100
101 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
102
103 void Run::SumEdepPerAbsorber(G4int kAbs, G4double EAbs, G4double LAbs)
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........oooOO0OOooo........oooOO0OOooo......
114
115 void Run::SumEnergies(G4double edeptot, G4double eleak0, G4double eleak1)
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........oooOO0OOooo........oooOO0OOooo......
132
133 void Run::SumEnergyFlow(G4int plane, G4double Eflow)
134 {
135 fEnergyFlow[plane] += Eflow;
136 }
137
138 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
139
140 void Run::Merge(const G4Run* run)
141 {
142 const Run* localRun = static_cast<const Run*>(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()) * (fDetector->GetNbOfAbsor()) + 2;
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 itp;
176 for (itp = localRun->fProcCounter.begin(); itp != localRun->fProcCounter.end(); ++itp) {
177 G4String procName = itp->first;
178 G4int localCount = itp->second;
179 if (fProcCounter.find(procName) == fProcCounter.end()) {
180 fProcCounter[procName] = localCount;
181 }
182 else {
183 fProcCounter[procName] += localCount;
184 }
185 }
186
187 G4Run::Merge(run);
188 }
189
190 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
191
192 void Run::EndOfRun()
193 {
194 // run condition
195 //
196 G4String Particle = fParticle->GetParticleName();
197 G4cout << "\n ---> The run is " << numberOfEvent << " " << Particle << " of "
198 << G4BestUnit(fEkin, "Energy") << " through calorimeter" << G4endl;
199
200 // frequency of processes
201 //
202 G4cout << "\n Process calls frequency :" << G4endl;
203 G4int index = 0;
204 std::map<G4String, G4int>::iterator it;
205 for (it = fProcCounter.begin(); it != fProcCounter.end(); it++) {
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 << "=" << std::setw(10) << count << space;
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, resolution, rmsres;
225 G4double MeanLAbs, MeanLAbs2, rmsLAbs;
226
227 std::ios::fmtflags mode = G4cout.flags();
228 G4int prec = G4cout.precision(2);
229 G4cout << "\n------------------------------------------------------------\n";
230 G4cout << std::setw(16) << "material" << std::setw(22) << "Edep rmsE" << std::setw(31)
231 << "sqrt(E0(GeV))*rmsE/Edep" << std::setw(23) << "total tracklen \n \n";
232
233 for (G4int k = 1; k <= fDetector->GetNbOfAbsor(); k++) {
234 MeanEAbs = fSumEAbs[k] * norm;
235 MeanEAbs2 = fSum2EAbs[k] * norm;
236 rmsEAbs = std::sqrt(std::abs(MeanEAbs2 - MeanEAbs * MeanEAbs));
237
238 resolution = 100. * sqbeam * rmsEAbs / MeanEAbs;
239 rmsres = resolution * qnorm;
240
241 MeanLAbs = fSumLAbs[k] * norm;
242 MeanLAbs2 = fSum2LAbs[k] * norm;
243 rmsLAbs = std::sqrt(std::abs(MeanLAbs2 - MeanLAbs * MeanLAbs));
244
245 // print
246 //
247 G4cout << std::setw(2) << k << std::setw(14) << fDetector->GetAbsorMaterial(k)->GetName()
248 << std::setprecision(5) << std::setw(10) << G4BestUnit(MeanEAbs, "Energy")
249 << std::setprecision(4) << std::setw(8) << G4BestUnit(rmsEAbs, "Energy") << std::setw(10)
250 << resolution << " +- " << std::setprecision(3) << std::setw(5) << rmsres << " %"
251 << std::setprecision(4) << std::setw(12) << G4BestUnit(MeanLAbs, "Length") << " +- "
252 << std::setprecision(3) << std::setw(5) << G4BestUnit(rmsLAbs, "Length") << G4endl;
253 }
254
255 // total energy deposited
256 //
257 fEdepTot /= nEvt;
258 fEdepTot2 /= nEvt;
259 G4double rmsEdep = std::sqrt(std::abs(fEdepTot2 - fEdepTot * fEdepTot));
260
261 G4cout << "\n Total energy deposited = " << std::setprecision(4) << G4BestUnit(fEdepTot, "Energy")
262 << " +- " << G4BestUnit(rmsEdep, "Energy") << G4endl;
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(fEleakTot2 - fEleakTot * fEleakTot));
271
272 G4cout << " Leakage : primary = " << G4BestUnit(fEnergyLeak[0], "Energy")
273 << " secondaries = " << G4BestUnit(fEnergyLeak[1], "Energy")
274 << " ---> total = " << G4BestUnit(fEleakTot, "Energy") << " +- "
275 << G4BestUnit(rmsEleak, "Energy") << G4endl;
276
277 // total energy released
278 //
279 fEtotal /= nEvt;
280 fEtotal2 /= nEvt;
281 G4double rmsEtotal = std::sqrt(std::abs(fEtotal2 - fEtotal * fEtotal));
282
283 G4cout << " Total energy released : Edep + Eleak = " << G4BestUnit(fEtotal, "Energy") << " +- "
284 << G4BestUnit(rmsEtotal, "Energy") << G4endl;
285 G4cout << "------------------------------------------------------------\n";
286
287 // Energy flow
288 //
289 G4AnalysisManager* analysis = G4AnalysisManager::Instance();
290 G4int Idmax = (fDetector->GetNbOfLayers()) * (fDetector->GetNbOfAbsor());
291 for (G4int Id = 1; Id <= Idmax + 1; Id++) {
292 analysis->FillH1(2 * kMaxAbsor + 1, (G4double)Id, fEnergyFlow[Id] / nEvt);
293 }
294
295 // normalize histograms
296 //
297 for (G4int ih = kMaxAbsor + 1; ih < 2 * kMaxAbsor + 1; ih++) {
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........oooOO0OOooo........oooOO0OOooo......
309