| 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 "HistoManager.hh"
36 #include "PrimaryGeneratorAction.hh"
37
38 #include "G4PhysicalConstants.hh"
39 #include "G4SystemOfUnits.hh"
40 #include "G4UnitsTable.hh"
41
42 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
43
44 void Run::SetPrimary(G4ParticleDefinition* particle, G4double energy)
45 {
46 fParticle = particle;
47 fEkin = energy;
48 }
49
50 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
51
52 void Run::ParticleCount(G4String name, G4double Ekin, G4double meanLife)
53 {
54 std::map<G4String, ParticleData>::iterator it = fParticleDataMap.find(name);
55 if (it == fParticleDataMap.end()) {
56 fParticleDataMap[name] = ParticleData(1, Ekin, Ekin, Ekin, meanLife);
57 }
58 else {
59 ParticleData& data = it->second;
60 data.fCount++;
61 data.fEmean += Ekin;
62 // update min max
63 G4double emin = data.fEmin;
64 if (Ekin < emin) data.fEmin = Ekin;
65 G4double emax = data.fEmax;
66 if (Ekin > emax) data.fEmax = Ekin;
67 data.fTmean = meanLife;
68 }
69 }
70
71 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
72
73 void Run::SetTimeWindow(G4double t1, G4double t2)
74 {
75 fTimeWindow1 = t1;
76 fTimeWindow2 = t2;
77 }
78
79 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
80
81 void Run::CountInTimeWindow(G4String name, G4bool life1, G4bool life2, G4bool decay)
82 {
83 std::map<G4String, ActivityData>::iterator it = fActivityMap.find(name);
84 if (it == fActivityMap.end()) {
85 G4int n1(0), n2(0), nd(0);
86 if (life1) n1 = 1;
87 if (life2) n2 = 1;
88 if (decay) nd = 1;
89 fActivityMap[name] = ActivityData(n1, n2, nd);
90 }
91 else {
92 ActivityData& data = it->second;
93 if (life1) data.fNlife_t1++;
94 if (life2) data.fNlife_t2++;
95 if (decay) data.fNdecay_t1t2++;
96 }
97 }
98
99 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
100
101 void Run::Balance(G4double Ekin, G4double Pbal)
102 {
103 fDecayCount++;
104 fEkinTot[0] += Ekin;
105 // update min max
106 if (fDecayCount == 1) fEkinTot[1] = fEkinTot[2] = Ekin;
107 if (Ekin < fEkinTot[1]) fEkinTot[1] = Ekin;
108 if (Ekin > fEkinTot[2]) fEkinTot[2] = Ekin;
109
110 fPbalance[0] += Pbal;
111 // update min max
112 if (fDecayCount == 1) fPbalance[1] = fPbalance[2] = Pbal;
113 if (Pbal < fPbalance[1]) fPbalance[1] = Pbal;
114 if (Pbal > fPbalance[2]) fPbalance[2] = Pbal;
115 }
116
117 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
118
119 void Run::EventTiming(G4double time)
120 {
121 fTimeCount++;
122 fEventTime[0] += time;
123 if (fTimeCount == 1) fEventTime[1] = fEventTime[2] = time;
124 if (time < fEventTime[1]) fEventTime[1] = time;
125 if (time > fEventTime[2]) fEventTime[2] = time;
126 }
127
128 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
129
130 void Run::PrimaryTiming(G4double ptime)
131 {
132 fPrimaryTime += ptime;
133 }
134
135 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
136
137 void Run::EvisEvent(G4double Evis)
138 {
139 fEvisEvent[0] += Evis;
140 if (fTimeCount == 1) fEvisEvent[1] = fEvisEvent[2] = Evis;
141 if (Evis < fEvisEvent[1]) fEvisEvent[1] = Evis;
142 if (Evis > fEvisEvent[2]) fEvisEvent[2] = Evis;
143 }
144
145 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
146
147 void Run::Merge(const G4Run* run)
148 {
149 const Run* localRun = static_cast<const Run*>(run);
150
151 // primary particle info
152 //
153 fParticle = localRun->fParticle;
154 fEkin = localRun->fEkin;
155
156 // accumulate sums
157 //
158 fDecayCount += localRun->fDecayCount;
159 fTimeCount += localRun->fTimeCount;
160 fPrimaryTime += localRun->fPrimaryTime;
161
162 fEkinTot[0] += localRun->fEkinTot[0];
163 fPbalance[0] += localRun->fPbalance[0];
164 fEventTime[0] += localRun->fEventTime[0];
165 fEvisEvent[0] += localRun->fEvisEvent[0];
166
167 G4double min, max;
168 min = localRun->fEkinTot[1];
169 max = localRun->fEkinTot[2];
170 if (fEkinTot[1] > min) fEkinTot[1] = min;
171 if (fEkinTot[2] < max) fEkinTot[2] = max;
172 //
173 min = localRun->fPbalance[1];
174 max = localRun->fPbalance[2];
175 if (fPbalance[1] > min) fPbalance[1] = min;
176 if (fPbalance[2] < max) fPbalance[2] = max;
177 //
178 min = localRun->fEventTime[1];
179 max = localRun->fEventTime[2];
180 if (fEventTime[1] > min) fEventTime[1] = min;
181 if (fEventTime[2] < max) fEventTime[2] = max;
182 //
183 min = localRun->fEvisEvent[1];
184 max = localRun->fEvisEvent[2];
185 if (fEvisEvent[1] > min) fEvisEvent[1] = min;
186 if (fEvisEvent[2] < max) fEvisEvent[2] = max;
187
188 // maps
189 std::map<G4String, ParticleData>::const_iterator itn;
190 for (itn = localRun->fParticleDataMap.begin(); itn != localRun->fParticleDataMap.end(); ++itn) {
191 G4String name = itn->first;
192 const ParticleData& localData = itn->second;
193 if (fParticleDataMap.find(name) == fParticleDataMap.end()) {
194 fParticleDataMap[name] = ParticleData(localData.fCount, localData.fEmean, localData.fEmin,
195 localData.fEmax, localData.fTmean);
196 }
197 else {
198 ParticleData& data = fParticleDataMap[name];
199 data.fCount += localData.fCount;
200 data.fEmean += localData.fEmean;
201 G4double emin = localData.fEmin;
202 if (emin < data.fEmin) data.fEmin = emin;
203 G4double emax = localData.fEmax;
204 if (emax > data.fEmax) data.fEmax = emax;
205 data.fTmean = localData.fTmean;
206 }
207 }
208
209 // activity
210 fTimeWindow1 = localRun->fTimeWindow1;
211 fTimeWindow2 = localRun->fTimeWindow2;
212
213 std::map<G4String, ActivityData>::const_iterator ita;
214 for (ita = localRun->fActivityMap.begin(); ita != localRun->fActivityMap.end(); ++ita) {
215 G4String name = ita->first;
216 const ActivityData& localData = ita->second;
217 if (fActivityMap.find(name) == fActivityMap.end()) {
218 fActivityMap[name] =
219 ActivityData(localData.fNlife_t1, localData.fNlife_t2, localData.fNdecay_t1t2);
220 }
221 else {
222 ActivityData& data = fActivityMap[name];
223 data.fNlife_t1 += localData.fNlife_t1;
224 data.fNlife_t2 += localData.fNlife_t2;
225 data.fNdecay_t1t2 += localData.fNdecay_t1t2;
226 }
227 }
228
229 G4Run::Merge(run);
230 }
231
232 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
233
234 void Run::EndOfRun()
235 {
236 G4int nbEvents = numberOfEvent;
237 G4String partName = fParticle->GetParticleName();
238
239 G4cout << "\n ======================== run summary ======================";
240 G4cout << "\n The run was " << nbEvents << " " << partName << " of "
241 << G4BestUnit(fEkin, "Energy");
242 G4cout << "\n ===========================================================\n";
243 G4cout << G4endl;
244 if (nbEvents == 0) {
245 return;
246 }
247
248 G4int prec = 4, wid = prec + 2;
249 G4int dfprec = G4cout.precision(prec);
250
251 // particle count
252 //
253 G4cout << " Nb of generated particles: \n" << G4endl;
254
255 std::map<G4String, ParticleData>::iterator it;
256 for (it = fParticleDataMap.begin(); it != fParticleDataMap.end(); it++) {
257 G4String name = it->first;
258 ParticleData data = it->second;
259 G4int count = data.fCount;
260 G4double eMean = data.fEmean / count;
261 G4double eMin = data.fEmin;
262 G4double eMax = data.fEmax;
263 G4double meanLife = data.fTmean;
264
265 G4cout << " " << std::setw(15) << name << ": " << std::setw(7) << count
266 << " Emean = " << std::setw(wid) << G4BestUnit(eMean, "Energy") << "\t( "
267 << G4BestUnit(eMin, "Energy") << " --> " << G4BestUnit(eMax, "Energy") << ")";
268 if (meanLife > 0.)
269 G4cout << "\tmean life = " << G4BestUnit(meanLife, "Time") << G4endl;
270 else if (meanLife < 0.)
271 G4cout << "\tstable" << G4endl;
272 else
273 G4cout << G4endl;
274 }
275
276 // energy momentum balance
277 //
278
279 if (fDecayCount > 0) {
280 G4double Ebmean = fEkinTot[0] / fDecayCount;
281 G4double Pbmean = fPbalance[0] / fDecayCount;
282
283 G4cout << "\n Ekin Total (Q single decay): mean = " << std::setw(wid)
284 << G4BestUnit(Ebmean, "Energy") << "\t( " << G4BestUnit(fEkinTot[1], "Energy") << " --> "
285 << G4BestUnit(fEkinTot[2], "Energy") << ")" << G4endl;
286
287 G4cout << "\n Momentum balance (excluding gamma desexcitation): mean = " << std::setw(wid)
288 << G4BestUnit(Pbmean, "Energy") << "\t( " << G4BestUnit(fPbalance[1], "Energy")
289 << " --> " << G4BestUnit(fPbalance[2], "Energy") << ")" << G4endl;
290 }
291
292 // total time of life
293 //
294 if (fTimeCount > 0) {
295 G4double Tmean = fEventTime[0] / fTimeCount;
296 G4double halfLife = Tmean * std::log(2.);
297
298 G4cout << "\n Total time of life (full chain): mean = " << std::setw(wid)
299 << G4BestUnit(Tmean, "Time") << " half-life = " << std::setw(wid)
300 << G4BestUnit(halfLife, "Time") << " ( " << G4BestUnit(fEventTime[1], "Time")
301 << " --> " << G4BestUnit(fEventTime[2], "Time") << ")" << G4endl;
302 }
303
304 // total visible Energy
305 //
306 if (fTimeCount > 0) {
307 G4double Evmean = fEvisEvent[0] / fTimeCount;
308
309 G4cout << "\n Total visible energy (full chain) : mean = " << std::setw(wid)
310 << G4BestUnit(Evmean, "Energy") << " ( " << G4BestUnit(fEvisEvent[1], "Energy")
311 << " --> " << G4BestUnit(fEvisEvent[2], "Energy") << ")" << G4endl;
312 }
313
314 // activity of primary ion
315 //
316 G4double pTimeMean = fPrimaryTime / nbEvents;
317 G4double molMass = fParticle->GetAtomicMass() * g / mole;
318 G4double nAtoms_perUnitOfMass = Avogadro / molMass;
319 G4double Activity_perUnitOfMass = 0.0;
320 if (pTimeMean > 0.0) {
321 Activity_perUnitOfMass = nAtoms_perUnitOfMass / pTimeMean;
322 }
323
324 G4cout << "\n Activity of " << partName << " = " << std::setw(wid)
325 << Activity_perUnitOfMass * g / becquerel << " Bq/g ("
326 << Activity_perUnitOfMass * g / curie << " Ci/g) \n"
327 << G4endl;
328
329 // activities in time window
330 //
331 if (fTimeWindow2 > 0.) {
332 G4double dt = fTimeWindow2 - fTimeWindow1;
333 G4cout << " Activities in time window [t1, t2] = [" << G4BestUnit(fTimeWindow1, "Time")
334 << ", " << G4BestUnit(fTimeWindow2, "Time")
335 << "] (delta time = " << G4BestUnit(dt, "Time") << ") : \n"
336 << G4endl;
337
338 std::map<G4String, ActivityData>::iterator ita;
339 for (ita = fActivityMap.begin(); ita != fActivityMap.end(); ita++) {
340 G4String name = ita->first;
341 ActivityData data = ita->second;
342 G4int n1 = data.fNlife_t1;
343 G4int n2 = data.fNlife_t2;
344 G4int ndecay = data.fNdecay_t1t2;
345 G4double actv = ndecay / dt;
346
347 G4cout << " " << std::setw(15) << name << ": "
348 << " n(t1) = " << std::setw(7) << n1 << "\tn(t2) = " << std::setw(7) << n2
349 << "\t decays = " << std::setw(7) << ndecay
350 << " ---> <actv> = " << G4BestUnit(actv, "Activity") << "\n";
351 }
352 }
353 G4cout << G4endl;
354
355 // normalize histograms
356 //
357 G4AnalysisManager* analysisManager = G4AnalysisManager::Instance();
358 G4double factor = 100. / nbEvents;
359 analysisManager->ScaleH1(1, factor);
360 analysisManager->ScaleH1(2, factor);
361 analysisManager->ScaleH1(3, factor);
362 analysisManager->ScaleH1(4, factor);
363 analysisManager->ScaleH1(5, factor);
364
365 // remove all contents in fParticleDataMap
366 //
367 fParticleDataMap.clear();
368 fActivityMap.clear();
369
370 // restore default precision
371 //
372 G4cout.precision(dfprec);
373 }
374
375 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
376