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