Geant4  10.01
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] = 0. ;
50  fEkinTot[1] = fPbalance[1] = fEventTime[1] = DBL_MAX;
51  fEkinTot[2] = fPbalance[2] = fEventTime[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 
124 void Run::Merge(const G4Run* run)
125 {
126  const Run* localRun = static_cast<const Run*>(run);
127 
128  //primary particle info
129  //
130  fParticle = localRun->fParticle;
131  fEkin = localRun->fEkin;
132 
133  // accumulate sums
134  //
135  fDecayCount += localRun->fDecayCount;
136  fTimeCount += localRun->fTimeCount;
137  fPrimaryTime += localRun->fPrimaryTime;
138 
139  fEkinTot[0] += localRun->fEkinTot[0];
140  fPbalance[0] += localRun->fPbalance[0];
141  fEventTime[0] += localRun->fEventTime[0];
142 
143  G4double min,max;
144  min = localRun->fEkinTot[1]; max = localRun->fEkinTot[2];
145  if (fEkinTot[1] > min) fEkinTot[1] = min;
146  if (fEkinTot[2] < max) fEkinTot[2] = max;
147  //
148  min = localRun->fPbalance[1]; max = localRun->fPbalance[2];
149  if (fPbalance[1] > min) fPbalance[1] = min;
150  if (fPbalance[2] < max) fPbalance[2] = max;
151  //
152  min = localRun->fEventTime[1]; max = localRun->fEventTime[2];
153  if (fEventTime[1] > min) fEventTime[1] = min;
154  if (fEventTime[2] < max) fEventTime[2] = max;
155 
156  //maps
157  std::map<G4String,ParticleData>::const_iterator itn;
158  for (itn = localRun->fParticleDataMap.begin();
159  itn != localRun->fParticleDataMap.end(); ++itn) {
160 
161  G4String name = itn->first;
162  const ParticleData& localData = itn->second;
163  if ( fParticleDataMap.find(name) == fParticleDataMap.end()) {
165  = ParticleData(localData.fCount,
166  localData.fEmean,
167  localData.fEmin,
168  localData.fEmax);
169  }
170  else {
171  ParticleData& data = fParticleDataMap[name];
172  data.fCount += localData.fCount;
173  data.fEmean += localData.fEmean;
174  G4double emin = localData.fEmin;
175  if (emin < data.fEmin) data.fEmin = emin;
176  G4double emax = localData.fEmax;
177  if (emax > data.fEmax) data.fEmax = emax;
178  }
179  }
180 
181  G4Run::Merge(run);
182 }
183 
184 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
185 
186 void Run::EndOfRun()
187 {
188  G4int nbEvents = numberOfEvent;
189  G4String partName = fParticle->GetParticleName();
190 
191  G4cout << "\n ======================== run summary ======================";
192  G4cout << "\n The run was " << nbEvents << " " << partName << " of "
193  << G4BestUnit(fEkin,"Energy");
194  G4cout << "\n ===========================================================\n";
195  G4cout << G4endl;
196  if (nbEvents == 0) { return; }
197 
198  G4int prec = 4, wid = prec + 2;
199  G4int dfprec = G4cout.precision(prec);
200 
201  //particle count
202  //
203  G4cout << " Nb of generated particles: \n" << G4endl;
204 
205  std::map<G4String,ParticleData>::iterator it;
206  for (it = fParticleDataMap.begin(); it != fParticleDataMap.end(); it++) {
207  G4String name = it->first;
208  ParticleData data = it->second;
209  G4int count = data.fCount;
210  G4double eMean = data.fEmean/count;
211  G4double eMin = data.fEmin;
212  G4double eMax = data.fEmax;
213 
214  G4cout << " " << std::setw(13) << name << ": " << std::setw(7) << count
215  << " Emean = " << std::setw(wid) << G4BestUnit(eMean, "Energy")
216  << "\t( " << G4BestUnit(eMin, "Energy")
217  << " --> " << G4BestUnit(eMax, "Energy")
218  << ")" << G4endl;
219  }
220 
221  //energy momentum balance
222  //
223 
224  if (fDecayCount > 0) {
225  G4double Ebmean = fEkinTot[0]/fDecayCount;
226  G4double Pbmean = fPbalance[0]/fDecayCount;
227 
228  G4cout << "\n Ekin Total (Q): mean = "
229  << std::setw(wid) << G4BestUnit(Ebmean, "Energy")
230  << "\t( " << G4BestUnit(fEkinTot[1], "Energy")
231  << " --> " << G4BestUnit(fEkinTot[2], "Energy")
232  << ")" << G4endl;
233 
234  G4cout << "\n Momentum balance (excluding gamma desexcitation): mean = "
235  << std::setw(wid) << G4BestUnit(Pbmean, "Energy")
236  << "\t( " << G4BestUnit(fPbalance[1], "Energy")
237  << " --> " << G4BestUnit(fPbalance[2], "Energy")
238  << ")" << G4endl;
239  }
240 
241  //total time of life
242  //
243  if (fTimeCount > 0) {
244  G4double Tmean = fEventTime[0]/fTimeCount;
245  G4double halfLife = Tmean*std::log(2.);
246 
247  G4cout << "\n Total time of life : mean = "
248  << std::setw(wid) << G4BestUnit(Tmean, "Time")
249  << " half-life = "
250  << std::setw(wid) << G4BestUnit(halfLife, "Time")
251  << " ( " << G4BestUnit(fEventTime[1], "Time")
252  << " --> " << G4BestUnit(fEventTime[2], "Time")
253  << ")" << G4endl;
254  }
255 
256  //activity of primary ion
257  //
258  G4double pTimeMean = fPrimaryTime/nbEvents;
259  G4double molMass = fParticle->GetAtomicMass()*g/mole;
260  G4double nAtoms_perUnitOfMass = Avogadro/molMass;
261  G4double Activity_perUnitOfMass = 0.0;
262  if (pTimeMean > 0.0)
263  { Activity_perUnitOfMass = nAtoms_perUnitOfMass/pTimeMean; }
264 
265  G4cout << "\n Activity of " << partName << " = "
266  << std::setw(wid) << Activity_perUnitOfMass*g/becquerel
267  << " Bq/g (" << Activity_perUnitOfMass*g/curie
268  << " Ci/g) \n"
269  << G4endl;
270 
271  //normalize histograms
272  //
273  G4AnalysisManager* analysisManager = G4AnalysisManager::Instance();
274  G4double factor = 100./nbEvents;
275  analysisManager->ScaleH1(1,factor);
276  analysisManager->ScaleH1(2,factor);
277  analysisManager->ScaleH1(3,factor);
278  analysisManager->ScaleH1(4,factor);
279  analysisManager->ScaleH1(5,factor);
280 
281  // remove all contents in fParticleDataMap
282  //
283  fParticleDataMap.clear();
284 
285  // restore default precision
286  //
287  G4cout.precision(dfprec);
288 }
289 
290 //....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:83
G4double fPrimaryTime
Definition: Run.hh:84
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
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:270
Definition: G4Run.hh:46
ExG4HbookAnalysisManager G4AnalysisManager
Definition: g4hbook_defs.hh:61
G4double fEkinTot[3]
Definition: Run.hh:81
G4int GetAtomicMass() const
static const double curie
Definition: G4SIunits.hh:271
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:162
T min(const T t1, const T t2)
brief Return the smallest of the two arguments
static const double mole
Definition: G4SIunits.hh:265
G4int fTimeCount
Definition: Run.hh:80
void EndOfRun()
Definition: Run.cc:147
#define G4endl
Definition: G4ios.hh:61
G4ParticleDefinition * fParticle
Definition: Run.hh:64
G4int fDecayCount
Definition: Run.hh:80
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