Geant4  10.02.p01
Run.cc
Go to the documentation of this file.
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 // $Id: Run.cc 71376 2013-06-14 07:44:50Z maire $
30 //
31 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
32 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
33 
34 #include "Run.hh"
35 #include "PrimaryGeneratorAction.hh"
36 #include "HistoManager.hh"
37 
38 #include "G4SystemOfUnits.hh"
39 #include "G4UnitsTable.hh"
40 #include "G4PhysicalConstants.hh"
41 
42 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
43 
44 Run::Run()
45 : G4Run(),
46  fParticle(0), fEkin(0.),
47  fDecayCount(0), fTimeCount(0), fPrimaryTime(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 
61 {
62  fParticle = particle;
63  fEkin = energy;
64 }
65 
66 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
67 
69 {
70  std::map<G4String, ParticleData>::iterator it = fParticleDataMap.find(name);
71  if ( it == fParticleDataMap.end()) {
72  fParticleDataMap[name] = ParticleData(1, Ekin, Ekin, Ekin);
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  }
84 }
85 
86 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
87 
88 void Run::Balance(G4double Ekin, G4double Pbal)
89 {
90  fDecayCount++;
91  fEkinTot[0] += Ekin;
92  //update min max
93  if (fDecayCount == 1) fEkinTot[1] = fEkinTot[2] = Ekin;
94  if (Ekin < fEkinTot[1]) fEkinTot[1] = Ekin;
95  if (Ekin > fEkinTot[2]) fEkinTot[2] = Ekin;
96 
97  fPbalance[0] += Pbal;
98  //update min max
99  if (fDecayCount == 1) fPbalance[1] = fPbalance[2] = Pbal;
100  if (Pbal < fPbalance[1]) fPbalance[1] = Pbal;
101  if (Pbal > fPbalance[2]) fPbalance[2] = Pbal;
102 }
103 
104 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
105 
107 {
108  fTimeCount++;
109  fEventTime[0] += time;
110  if (fTimeCount == 1) fEventTime[1] = fEventTime[2] = time;
111  if (time < fEventTime[1]) fEventTime[1] = time;
112  if (time > fEventTime[2]) fEventTime[2] = time;
113 }
114 
115 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
116 
118 {
119  fPrimaryTime += ptime;
120 }
121 
122 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
123 
125 {
126  fEvisEvent[0] += Evis;
127  if (fTimeCount == 1) fEvisEvent[1] = fEvisEvent[2] = Evis;
128  if (Evis < fEvisEvent[1]) fEvisEvent[1] = Evis;
129  if (Evis > fEvisEvent[2]) fEvisEvent[2] = Evis;
130 }
131 
132 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
133 
134 void Run::Merge(const G4Run* run)
135 {
136  const Run* localRun = static_cast<const Run*>(run);
137 
138  //primary particle info
139  //
140  fParticle = localRun->fParticle;
141  fEkin = localRun->fEkin;
142 
143  // accumulate sums
144  //
145  fDecayCount += localRun->fDecayCount;
146  fTimeCount += localRun->fTimeCount;
147  fPrimaryTime += localRun->fPrimaryTime;
148 
149  fEkinTot[0] += localRun->fEkinTot[0];
150  fPbalance[0] += localRun->fPbalance[0];
151  fEventTime[0] += localRun->fEventTime[0];
152  fEvisEvent[0] += localRun->fEvisEvent[0];
153 
154  G4double min,max;
155  min = localRun->fEkinTot[1]; max = localRun->fEkinTot[2];
156  if (fEkinTot[1] > min) fEkinTot[1] = min;
157  if (fEkinTot[2] < max) fEkinTot[2] = max;
158  //
159  min = localRun->fPbalance[1]; max = localRun->fPbalance[2];
160  if (fPbalance[1] > min) fPbalance[1] = min;
161  if (fPbalance[2] < max) fPbalance[2] = max;
162  //
163  min = localRun->fEventTime[1]; max = localRun->fEventTime[2];
164  if (fEventTime[1] > min) fEventTime[1] = min;
165  if (fEventTime[2] < max) fEventTime[2] = max;
166  //
167  min = localRun->fEvisEvent[1]; max = localRun->fEvisEvent[2];
168  if (fEvisEvent[1] > min) fEvisEvent[1] = min;
169  if (fEvisEvent[2] < max) fEvisEvent[2] = max;
170 
171  //maps
172  std::map<G4String,ParticleData>::const_iterator itn;
173  for (itn = localRun->fParticleDataMap.begin();
174  itn != localRun->fParticleDataMap.end(); ++itn) {
175 
176  G4String name = itn->first;
177  const ParticleData& localData = itn->second;
178  if ( fParticleDataMap.find(name) == fParticleDataMap.end()) {
180  = ParticleData(localData.fCount,
181  localData.fEmean,
182  localData.fEmin,
183  localData.fEmax);
184  }
185  else {
186  ParticleData& data = fParticleDataMap[name];
187  data.fCount += localData.fCount;
188  data.fEmean += localData.fEmean;
189  G4double emin = localData.fEmin;
190  if (emin < data.fEmin) data.fEmin = emin;
191  G4double emax = localData.fEmax;
192  if (emax > data.fEmax) data.fEmax = emax;
193  }
194  }
195 
196  G4Run::Merge(run);
197 }
198 
199 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
200 
201 void Run::EndOfRun()
202 {
203  G4int nbEvents = numberOfEvent;
204  G4String partName = fParticle->GetParticleName();
205 
206  G4cout << "\n ======================== run summary ======================";
207  G4cout << "\n The run was " << nbEvents << " " << partName << " of "
208  << G4BestUnit(fEkin,"Energy");
209  G4cout << "\n ===========================================================\n";
210  G4cout << G4endl;
211  if (nbEvents == 0) { return; }
212 
213  G4int prec = 4, wid = prec + 2;
214  G4int dfprec = G4cout.precision(prec);
215 
216  //particle count
217  //
218  G4cout << " Nb of generated particles: \n" << G4endl;
219 
220  std::map<G4String,ParticleData>::iterator it;
221  for (it = fParticleDataMap.begin(); it != fParticleDataMap.end(); it++) {
222  G4String name = it->first;
223  ParticleData data = it->second;
224  G4int count = data.fCount;
225  G4double eMean = data.fEmean/count;
226  G4double eMin = data.fEmin;
227  G4double eMax = data.fEmax;
228 
229  G4cout << " " << std::setw(13) << name << ": " << std::setw(7) << count
230  << " Emean = " << std::setw(wid) << G4BestUnit(eMean, "Energy")
231  << "\t( " << G4BestUnit(eMin, "Energy")
232  << " --> " << G4BestUnit(eMax, "Energy")
233  << ")" << G4endl;
234  }
235 
236  //energy momentum balance
237  //
238 
239  if (fDecayCount > 0) {
240  G4double Ebmean = fEkinTot[0]/fDecayCount;
241  G4double Pbmean = fPbalance[0]/fDecayCount;
242 
243  G4cout << "\n Ekin Total (Q single decay): mean = "
244  << std::setw(wid) << G4BestUnit(Ebmean, "Energy")
245  << "\t( " << G4BestUnit(fEkinTot[1], "Energy")
246  << " --> " << G4BestUnit(fEkinTot[2], "Energy")
247  << ")" << G4endl;
248 
249  G4cout << "\n Momentum balance (excluding gamma desexcitation): mean = "
250  << std::setw(wid) << G4BestUnit(Pbmean, "Energy")
251  << "\t( " << G4BestUnit(fPbalance[1], "Energy")
252  << " --> " << G4BestUnit(fPbalance[2], "Energy")
253  << ")" << G4endl;
254  }
255 
256  //total time of life
257  //
258  if (fTimeCount > 0) {
259  G4double Tmean = fEventTime[0]/fTimeCount;
260  G4double halfLife = Tmean*std::log(2.);
261 
262  G4cout << "\n Total time of life (full chain): mean = "
263  << std::setw(wid) << G4BestUnit(Tmean, "Time")
264  << " half-life = "
265  << std::setw(wid) << G4BestUnit(halfLife, "Time")
266  << " ( " << G4BestUnit(fEventTime[1], "Time")
267  << " --> " << G4BestUnit(fEventTime[2], "Time")
268  << ")" << G4endl;
269  }
270 
271  //total visible Energy
272  //
273  if (fTimeCount > 0) {
274  G4double Evmean = fEvisEvent[0]/fTimeCount;
275 
276  G4cout << "\n Total visible energy (full chain) : mean = "
277  << std::setw(wid) << G4BestUnit(Evmean, "Energy")
278  << " ( " << G4BestUnit(fEvisEvent[1], "Energy")
279  << " --> " << G4BestUnit(fEvisEvent[2], "Energy")
280  << ")" << G4endl;
281  }
282 
283  //activity of primary ion
284  //
285  G4double pTimeMean = fPrimaryTime/nbEvents;
286  G4double molMass = fParticle->GetAtomicMass()*g/mole;
287  G4double nAtoms_perUnitOfMass = Avogadro/molMass;
288  G4double Activity_perUnitOfMass = 0.0;
289  if (pTimeMean > 0.0)
290  { Activity_perUnitOfMass = nAtoms_perUnitOfMass/pTimeMean; }
291 
292  G4cout << "\n Activity of " << partName << " = "
293  << std::setw(wid) << Activity_perUnitOfMass*g/becquerel
294  << " Bq/g (" << Activity_perUnitOfMass*g/curie
295  << " Ci/g) \n"
296  << G4endl;
297 
298  //normalize histograms
299  //
300  G4AnalysisManager* analysisManager = G4AnalysisManager::Instance();
301  G4double factor = 100./nbEvents;
302  analysisManager->ScaleH1(1,factor);
303  analysisManager->ScaleH1(2,factor);
304  analysisManager->ScaleH1(3,factor);
305  analysisManager->ScaleH1(4,factor);
306  analysisManager->ScaleH1(5,factor);
307 
308  // remove all contents in fParticleDataMap
309  //
310  fParticleDataMap.clear();
311 
312  // restore default precision
313  //
314  G4cout.precision(dfprec);
315 }
316 
317 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
std::map< G4String, ParticleData > fParticleDataMap
Definition: Run.hh:100
virtual void Merge(const G4Run *)
Definition: G4Run.cc:54
G4int numberOfEvent
Definition: G4Run.hh:59
G4double fEkin
Definition: Run.hh:65
void PrimaryTiming(G4double)
Definition: Run.cc:117
G4double fEventTime[3]
Definition: Run.hh:84
G4double fPrimaryTime
Definition: Run.hh:85
Run()
Definition: Run.cc:43
G4int first(char) const
G4String name
Definition: TRTMaterials.hh:40
void EventTiming(G4double)
Definition: Run.cc:106
#define G4BestUnit(a, b)
#define G4_USE_G4BESTUNIT_FOR_VERBOSE 1
void Balance(G4double)
Definition: Run.cc:133
int G4int
Definition: G4Types.hh:78
const G4String & GetParticleName() const
void EvisEvent(G4double)
Definition: Run.cc:124
static const double prec
Definition: RanecuEngine.cc:58
G4GLOB_DLL std::ostream G4cout
G4double fPbalance[3]
Definition: Run.hh:103
static const double becquerel
Definition: G4SIunits.hh:288
Definition: G4Run.hh:46
ExG4HbookAnalysisManager G4AnalysisManager
Definition: g4hbook_defs.hh:66
G4double fEkinTot[3]
Definition: Run.hh:82
G4int GetAtomicMass() const
static const double curie
Definition: G4SIunits.hh:289
static const G4double emax
static const G4double factor
T max(const T t1, const T t2)
brief Return the largest of the two arguments
G4double energy(const ThreeVector &p, const G4double m)
static const double g
Definition: G4SIunits.hh:180
G4double fEvisEvent[3]
Definition: Run.hh:86
T min(const T t1, const T t2)
brief Return the smallest of the two arguments
static const double mole
Definition: G4SIunits.hh:283
G4int fTimeCount
Definition: Run.hh:81
void EndOfRun()
Definition: Run.cc:147
#define G4endl
Definition: G4ios.hh:61
G4ParticleDefinition * fParticle
Definition: Run.hh:64
G4int fDecayCount
Definition: Run.hh:81
void SetPrimary(G4ParticleDefinition *particle, G4double energy)
Definition: Run.cc:77
double G4double
Definition: G4Types.hh:76
Definition: Run.hh:46
#define DBL_MAX
Definition: templates.hh:83
virtual void Merge(const G4Run *)
Definition: Run.cc:115
~Run()
Definition: Run.cc:72
void ParticleCount(G4String, G4double)
Definition: Run.cc:114