| 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 "G4AutoLock.hh"
40 #include "G4SystemOfUnits.hh"
41 #include "G4Threading.hh"
42 #include "G4UnitsTable.hh"
43
44 // mutex in a file scope
45
46 namespace
47 {
48 // Mutex to lock updating the global ion map
49 G4Mutex ionIdMapMutex = G4MUTEX_INITIALIZER;
50 } // namespace
51
52 //....oooOO0OOooo........oooOO0OOooo........oo
53
54 std::map<G4String, G4int> Run::fgIonMap;
55 G4int Run::fgIonId = kMaxHisto1;
56
57 //....oooOO0OOooo........oooOO0OOooo........oo
58
59 Run::Run(DetectorConstruction* det) : fDetecto
60
61 //....oooOO0OOooo........oooOO0OOooo........oo
62
63 void Run::Merge(std::map<G4String, ParticleDat
64 const std::map<G4String, Parti
65 {
66 for (const auto& particleData : sourceMap) {
67 G4String name = particleData.first;
68 const ParticleData& localData = particleDa
69 if (destinationMap.find(name) == destinati
70 destinationMap[name] = ParticleData(loca
71 loca
72 }
73 else {
74 ParticleData& data = destinationMap[name
75 data.fCount += localData.fCount;
76 data.fEmean += localData.fEmean;
77 G4double emin = localData.fEmin;
78 if (emin < data.fEmin) data.fEmin = emin
79 G4double emax = localData.fEmax;
80 if (emax > data.fEmax) data.fEmax = emax
81 data.fTmean = localData.fTmean;
82 }
83 }
84 }
85
86 //....oooOO0OOooo........oooOO0OOooo........oo
87
88 void Run::SetPrimary(G4ParticleDefinition* par
89 {
90 fParticle = particle;
91 fEkin = energy;
92 }
93
94 //....oooOO0OOooo........oooOO0OOooo........oo
95
96 void Run::CountProcesses(const G4VProcess* pro
97 {
98 if (process == nullptr) return;
99 G4String procName = process->GetProcessName(
100 std::map<G4String, G4int>::iterator it = fPr
101 if (it == fProcCounter.end()) {
102 fProcCounter[procName] = 1;
103 }
104 else {
105 fProcCounter[procName]++;
106 }
107 }
108
109 //....oooOO0OOooo........oooOO0OOooo........oo
110
111 void Run::ParticleCount(G4String name, G4doubl
112 {
113 std::map<G4String, ParticleData>::iterator i
114 if (it == fParticleDataMap1.end()) {
115 fParticleDataMap1[name] = ParticleData(1,
116 }
117 else {
118 ParticleData& data = it->second;
119 data.fCount++;
120 data.fEmean += Ekin;
121 // update min max
122 G4double emin = data.fEmin;
123 if (Ekin < emin) data.fEmin = Ekin;
124 G4double emax = data.fEmax;
125 if (Ekin > emax) data.fEmax = Ekin;
126 data.fTmean = meanLife;
127 }
128 }
129
130 //....oooOO0OOooo........oooOO0OOooo........oo
131
132 void Run::AddEdep(G4double edep)
133 {
134 fEnergyDeposit += edep;
135 fEnergyDeposit2 += edep * edep;
136 }
137
138 //....oooOO0OOooo........oooOO0OOooo........oo
139
140 void Run::AddEflow(G4double eflow)
141 {
142 fEnergyFlow += eflow;
143 fEnergyFlow2 += eflow * eflow;
144 }
145 //....oooOO0OOooo........oooOO0OOooo........oo
146
147 void Run::ParticleFlux(G4String name, G4double
148 {
149 std::map<G4String, ParticleData>::iterator i
150 if (it == fParticleDataMap2.end()) {
151 fParticleDataMap2[name] = ParticleData(1,
152 }
153 else {
154 ParticleData& data = it->second;
155 data.fCount++;
156 data.fEmean += Ekin;
157 // update min max
158 G4double emin = data.fEmin;
159 if (Ekin < emin) data.fEmin = Ekin;
160 G4double emax = data.fEmax;
161 if (Ekin > emax) data.fEmax = Ekin;
162 data.fTmean = -1 * ns;
163 }
164 }
165
166 //....oooOO0OOooo........oooOO0OOooo........oo
167
168 G4int Run::GetIonId(G4String ionName)
169 {
170 G4AutoLock lock(&ionIdMapMutex);
171 // updating the global ion map needs to be l
172
173 std::map<G4String, G4int>::const_iterator it
174 if (it == fgIonMap.end()) {
175 fgIonMap[ionName] = fgIonId;
176 if (fgIonId < (kMaxHisto2 - 1)) fgIonId++;
177 }
178 return fgIonMap[ionName];
179 }
180
181 //....oooOO0OOooo........oooOO0OOooo........oo
182
183 void Run::Merge(const G4Run* run)
184 {
185 const Run* localRun = static_cast<const Run*
186
187 // primary particle info
188 //
189 fParticle = localRun->fParticle;
190 fEkin = localRun->fEkin;
191
192 // accumulate sums
193 //
194 fEnergyDeposit += localRun->fEnergyDeposit;
195 fEnergyDeposit2 += localRun->fEnergyDeposit2
196 fEnergyFlow += localRun->fEnergyFlow;
197 fEnergyFlow2 += localRun->fEnergyFlow2;
198
199 // map: processes count
200 for (const auto& procCounter : localRun->fPr
201 G4String procName = procCounter.first;
202 G4int localCount = procCounter.second;
203 if (fProcCounter.find(procName) == fProcCo
204 fProcCounter[procName] = localCount;
205 }
206 else {
207 fProcCounter[procName] += localCount;
208 }
209 }
210
211 // map: created particles count
212 Merge(fParticleDataMap1, localRun->fParticle
213
214 // map: particles flux count
215 Merge(fParticleDataMap2, localRun->fParticle
216
217 G4Run::Merge(run);
218 }
219
220 //....oooOO0OOooo........oooOO0OOooo........oo
221
222 void Run::EndOfRun()
223 {
224 G4int prec = 5, wid = prec + 2;
225 G4int dfprec = G4cout.precision(prec);
226
227 // run condition
228 //
229 G4Material* material = fDetector->GetAbsorMa
230 G4double density = material->GetDensity();
231
232 G4String Particle = fParticle->GetParticleNa
233 G4cout << "\n The run is " << numberOfEvent
234 << G4BestUnit(fEkin, "Energy") << " t
235 << G4BestUnit(fDetector->GetAbsorThic
236 << " (density: " << G4BestUnit(densit
237
238 if (numberOfEvent == 0) {
239 G4cout.precision(dfprec);
240 return;
241 }
242
243 // frequency of processes
244 //
245 G4cout << "\n Process calls frequency :" <<
246 G4int index = 0;
247 for (const auto& procCounter : fProcCounter)
248 G4String procName = procCounter.first;
249 G4int count = procCounter.second;
250 G4String space = " ";
251 if (++index % 3 == 0) space = "\n";
252 G4cout << " " << std::setw(20) << procName
253 }
254 G4cout << G4endl;
255
256 // particles count
257 //
258 G4cout << "\n List of generated particles (w
259
260 for (const auto& particleData : fParticleDat
261 G4String name = particleData.first;
262 ParticleData data = particleData.second;
263 G4int count = data.fCount;
264 G4double eMean = data.fEmean / count;
265 G4double eMin = data.fEmin;
266 G4double eMax = data.fEmax;
267 G4double meanLife = data.fTmean;
268
269 G4cout << " " << std::setw(13) << name <<
270 << " Emean = " << std::setw(wid) <
271 << G4BestUnit(eMin, "Energy") << "
272 if (meanLife >= 0.)
273 G4cout << "\tmean life = " << G4BestUnit
274 else
275 G4cout << "\tstable" << G4endl;
276 }
277
278 // compute mean Energy deposited and rms
279 //
280 G4int TotNbofEvents = numberOfEvent;
281 fEnergyDeposit /= TotNbofEvents;
282 fEnergyDeposit2 /= TotNbofEvents;
283 G4double rmsEdep = fEnergyDeposit2 - fEnergy
284 if (rmsEdep > 0.)
285 rmsEdep = std::sqrt(rmsEdep);
286 else
287 rmsEdep = 0.;
288
289 G4cout << "\n Mean energy deposit per event
290 << "; rms = " << G4BestUnit(rmsEdep,
291
292 // compute mean Energy flow and rms
293 //
294 fEnergyFlow /= TotNbofEvents;
295 fEnergyFlow2 /= TotNbofEvents;
296 G4double rmsEflow = fEnergyFlow2 - fEnergyFl
297 if (rmsEflow > 0.)
298 rmsEflow = std::sqrt(rmsEflow);
299 else
300 rmsEflow = 0.;
301
302 G4cout << " Mean energy flow per event =
303 << "; rms = " << G4BestUnit(rmsEflow
304
305 // particles flux
306 //
307 G4cout << "\n List of particles emerging fro
308
309 for (const auto& particleData : fParticleDat
310 G4String name = particleData.first;
311 ParticleData data = particleData.second;
312 G4int count = data.fCount;
313 G4double eMean = data.fEmean / count;
314 G4double eMin = data.fEmin;
315 G4double eMax = data.fEmax;
316 G4double Eflow = data.fEmean / TotNbofEven
317
318 G4cout << " " << std::setw(13) << name <<
319 << " Emean = " << std::setw(wid) <
320 << G4BestUnit(eMin, "Energy") << "
321 << ") \tEflow/event = " << G4BestUn
322 }
323
324 // histogram Id for populations
325 //
326 G4cout << "\n histo Id for populations :" <<
327
328 // Update the histogram titles according to
329 // and print new titles
330 G4AnalysisManager* analysisManager = G4Analy
331 for (const auto& ionMapElement : fgIonMap) {
332 G4String ionName = ionMapElement.first;
333 G4int h1Id = ionMapElement.second;
334 // print new titles
335 G4cout << " " << std::setw(20) << ionName
336
337 // update histogram ids
338 if (!analysisManager->GetH1(h1Id)) continu
339 // Skip inactive histograms, this is not n
340 // but it makes the code safe wrt modific
341 G4String title = analysisManager->GetH1Tit
342 title = ionName + title;
343 analysisManager->SetH1Title(h1Id, title);
344 }
345 G4cout << G4endl;
346
347 // normalize histograms
348 G4int ih = 2;
349 G4double binWidth = analysisManager->GetH1Wi
350 G4double fac = (1. / (numberOfEvent * binWid
351 analysisManager->ScaleH1(ih, fac);
352
353 for (ih = 14; ih < 24; ih++) {
354 binWidth = analysisManager->GetH1Width(ih)
355 G4double unit = analysisManager->GetH1Unit
356 fac = (second / (binWidth * unit));
357 analysisManager->ScaleH1(ih, fac);
358 }
359
360 // remove all contents in fProcCounter, fCou
361 fProcCounter.clear();
362 fParticleDataMap1.clear();
363 fParticleDataMap2.clear();
364 fgIonMap.clear();
365
366 // restore default format
367 G4cout.precision(dfprec);
368 }
369
370 //....oooOO0OOooo........oooOO0OOooo........oo
371