| 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 "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........oooOO0OOooo........oooOO0OOooo......
53
54 std::map<G4String, G4int> Run::fgIonMap;
55 G4int Run::fgIonId = kMaxHisto1;
56
57 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
58
59 Run::Run(DetectorConstruction* det) : fDetector(det) {}
60
61 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
62
63 void Run::Merge(std::map<G4String, ParticleData>& destinationMap,
64 const std::map<G4String, ParticleData>& sourceMap) const
65 {
66 for (const auto& particleData : sourceMap) {
67 G4String name = particleData.first;
68 const ParticleData& localData = particleData.second;
69 if (destinationMap.find(name) == destinationMap.end()) {
70 destinationMap[name] = ParticleData(localData.fCount, localData.fEmean, localData.fEmin,
71 localData.fEmax, localData.fTmean);
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........oooOO0OOooo........oooOO0OOooo......
87
88 void Run::SetPrimary(G4ParticleDefinition* particle, G4double energy)
89 {
90 fParticle = particle;
91 fEkin = energy;
92 }
93
94 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
95
96 void Run::CountProcesses(const G4VProcess* process)
97 {
98 if (process == nullptr) return;
99 G4String procName = process->GetProcessName();
100 std::map<G4String, G4int>::iterator it = fProcCounter.find(procName);
101 if (it == fProcCounter.end()) {
102 fProcCounter[procName] = 1;
103 }
104 else {
105 fProcCounter[procName]++;
106 }
107 }
108
109 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
110
111 void Run::ParticleCount(G4String name, G4double Ekin, G4double meanLife)
112 {
113 std::map<G4String, ParticleData>::iterator it = fParticleDataMap1.find(name);
114 if (it == fParticleDataMap1.end()) {
115 fParticleDataMap1[name] = ParticleData(1, Ekin, Ekin, Ekin, meanLife);
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........oooOO0OOooo........oooOO0OOooo......
131
132 void Run::AddEdep(G4double edep)
133 {
134 fEnergyDeposit += edep;
135 fEnergyDeposit2 += edep * edep;
136 }
137
138 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
139
140 void Run::AddEflow(G4double eflow)
141 {
142 fEnergyFlow += eflow;
143 fEnergyFlow2 += eflow * eflow;
144 }
145 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
146
147 void Run::ParticleFlux(G4String name, G4double Ekin)
148 {
149 std::map<G4String, ParticleData>::iterator it = fParticleDataMap2.find(name);
150 if (it == fParticleDataMap2.end()) {
151 fParticleDataMap2[name] = ParticleData(1, Ekin, Ekin, Ekin, -1 * ns);
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........oooOO0OOooo........oooOO0OOooo......
167
168 G4int Run::GetIonId(G4String ionName)
169 {
170 G4AutoLock lock(&ionIdMapMutex);
171 // updating the global ion map needs to be locked
172
173 std::map<G4String, G4int>::const_iterator it = fgIonMap.find(ionName);
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........oooOO0OOooo........oooOO0OOooo......
182
183 void Run::Merge(const G4Run* run)
184 {
185 const Run* localRun = static_cast<const Run*>(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->fProcCounter) {
201 G4String procName = procCounter.first;
202 G4int localCount = procCounter.second;
203 if (fProcCounter.find(procName) == fProcCounter.end()) {
204 fProcCounter[procName] = localCount;
205 }
206 else {
207 fProcCounter[procName] += localCount;
208 }
209 }
210
211 // map: created particles count
212 Merge(fParticleDataMap1, localRun->fParticleDataMap1);
213
214 // map: particles flux count
215 Merge(fParticleDataMap2, localRun->fParticleDataMap2);
216
217 G4Run::Merge(run);
218 }
219
220 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
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->GetAbsorMaterial();
230 G4double density = material->GetDensity();
231
232 G4String Particle = fParticle->GetParticleName();
233 G4cout << "\n The run is " << numberOfEvent << " " << Particle << " of "
234 << G4BestUnit(fEkin, "Energy") << " through "
235 << G4BestUnit(fDetector->GetAbsorThickness(), "Length") << " of " << material->GetName()
236 << " (density: " << G4BestUnit(density, "Volumic Mass") << ")" << G4endl;
237
238 if (numberOfEvent == 0) {
239 G4cout.precision(dfprec);
240 return;
241 }
242
243 // frequency of processes
244 //
245 G4cout << "\n Process calls frequency :" << G4endl;
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 << "=" << std::setw(7) << count << space;
253 }
254 G4cout << G4endl;
255
256 // particles count
257 //
258 G4cout << "\n List of generated particles (with meanLife != 0):" << G4endl;
259
260 for (const auto& particleData : fParticleDataMap1) {
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 << ": " << std::setw(7) << count
270 << " Emean = " << std::setw(wid) << G4BestUnit(eMean, "Energy") << "\t( "
271 << G4BestUnit(eMin, "Energy") << " --> " << G4BestUnit(eMax, "Energy") << ")";
272 if (meanLife >= 0.)
273 G4cout << "\tmean life = " << G4BestUnit(meanLife, "Time") << G4endl;
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 - fEnergyDeposit * fEnergyDeposit;
284 if (rmsEdep > 0.)
285 rmsEdep = std::sqrt(rmsEdep);
286 else
287 rmsEdep = 0.;
288
289 G4cout << "\n Mean energy deposit per event = " << G4BestUnit(fEnergyDeposit, "Energy")
290 << "; rms = " << G4BestUnit(rmsEdep, "Energy") << G4endl;
291
292 // compute mean Energy flow and rms
293 //
294 fEnergyFlow /= TotNbofEvents;
295 fEnergyFlow2 /= TotNbofEvents;
296 G4double rmsEflow = fEnergyFlow2 - fEnergyFlow * fEnergyFlow;
297 if (rmsEflow > 0.)
298 rmsEflow = std::sqrt(rmsEflow);
299 else
300 rmsEflow = 0.;
301
302 G4cout << " Mean energy flow per event = " << G4BestUnit(fEnergyFlow, "Energy")
303 << "; rms = " << G4BestUnit(rmsEflow, "Energy") << G4endl;
304
305 // particles flux
306 //
307 G4cout << "\n List of particles emerging from the target :" << G4endl;
308
309 for (const auto& particleData : fParticleDataMap2) {
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 / TotNbofEvents;
317
318 G4cout << " " << std::setw(13) << name << ": " << std::setw(7) << count
319 << " Emean = " << std::setw(wid) << G4BestUnit(eMean, "Energy") << "\t( "
320 << G4BestUnit(eMin, "Energy") << " --> " << G4BestUnit(eMax, "Energy")
321 << ") \tEflow/event = " << G4BestUnit(Eflow, "Energy") << G4endl;
322 }
323
324 // histogram Id for populations
325 //
326 G4cout << "\n histo Id for populations :" << G4endl;
327
328 // Update the histogram titles according to the ion map
329 // and print new titles
330 G4AnalysisManager* analysisManager = G4AnalysisManager::Instance();
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 << " id = " << std::setw(3) << h1Id << G4endl;
336
337 // update histogram ids
338 if (!analysisManager->GetH1(h1Id)) continue;
339 // Skip inactive histograms, this is not necessary
340 // but it makes the code safe wrt modifications in future
341 G4String title = analysisManager->GetH1Title(h1Id);
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->GetH1Width(ih);
350 G4double fac = (1. / (numberOfEvent * binWidth)) * (mm / MeV);
351 analysisManager->ScaleH1(ih, fac);
352
353 for (ih = 14; ih < 24; ih++) {
354 binWidth = analysisManager->GetH1Width(ih);
355 G4double unit = analysisManager->GetH1Unit(ih);
356 fac = (second / (binWidth * unit));
357 analysisManager->ScaleH1(ih, fac);
358 }
359
360 // remove all contents in fProcCounter, fCount
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........oooOO0OOooo........oooOO0OOooo......
371