Geant4  9.6.p02
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
RunAction.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 //
30 // $Id$
31 //
32 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
33 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
34 
35 #include "RunAction.hh"
36 #include "HistoManager.hh"
37 #include "PrimaryGeneratorAction.hh"
38 
39 #include "G4Run.hh"
40 #include "G4RunManager.hh"
41 #include "G4UnitsTable.hh"
42 #include "G4PhysicalConstants.hh"
43 #include "G4SystemOfUnits.hh"
44 #include <iomanip>
45 
46 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
47 
48 RunAction::RunAction(PrimaryGeneratorAction* kin)
49 :fPrimary(kin)
50 {
51  fHistoManager = new HistoManager();
52 }
53 
54 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
55 
56 RunAction::~RunAction()
57 {
58  delete fHistoManager;
59 }
60 
61 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
62 
63 void RunAction::BeginOfRunAction(const G4Run*)
64 {
65  //initialize arrays
66  //
67  fDecayCount = fTimeCount = 0;
68  for (G4int i=0; i<3; i++) fEkinTot[i] = fPbalance[i] = fEventTime[i] = 0. ;
69  fPrimaryTime = 0.;
70 
71  //histograms
72  //
73  G4AnalysisManager* analysisManager = G4AnalysisManager::Instance();
74  if ( analysisManager->IsActive() ) {
75  analysisManager->OpenFile();
76  }
77 
78  //inform the runManager to save random number seed
79  //
80  G4RunManager::GetRunManager()->SetRandomNumberStore(false);
81 }
82 
83 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
84 
85 void RunAction::ParticleCount(G4String name, G4double Ekin)
86 {
87  fParticleCount[name]++;
88  fEmean[name] += Ekin;
89  //update min max
90  if (fParticleCount[name] == 1) fEmin[name] = fEmax[name] = Ekin;
91  if (Ekin < fEmin[name]) fEmin[name] = Ekin;
92  if (Ekin > fEmax[name]) fEmax[name] = Ekin;
93 }
94 
95 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
96 
97 void RunAction::Balance(G4double Ekin, G4double Pbal)
98 {
99  fDecayCount++;
100  fEkinTot[0] += Ekin;
101  //update min max
102  if (fDecayCount == 1) fEkinTot[1] = fEkinTot[2] = Ekin;
103  if (Ekin < fEkinTot[1]) fEkinTot[1] = Ekin;
104  if (Ekin > fEkinTot[2]) fEkinTot[2] = Ekin;
105 
106  fPbalance[0] += Pbal;
107  //update min max
108  if (fDecayCount == 1) fPbalance[1] = fPbalance[2] = Pbal;
109  if (Pbal < fPbalance[1]) fPbalance[1] = Pbal;
110  if (Pbal > fPbalance[2]) fPbalance[2] = Pbal;
111 }
112 
113 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
114 
115 void RunAction::EventTiming(G4double time)
116 {
117  fTimeCount++;
118  fEventTime[0] += time;
119  if (fTimeCount == 1) fEventTime[1] = fEventTime[2] = time;
120  if (time < fEventTime[1]) fEventTime[1] = time;
121  if (time > fEventTime[2]) fEventTime[2] = time;
122 }
123 
124 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
125 
126 void RunAction::PrimaryTiming(G4double ptime)
127 {
128  fPrimaryTime += ptime;
129 }
130 
131 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
132 
133 void RunAction::EndOfRunAction(const G4Run* run)
134 {
135  G4int nbEvents = run->GetNumberOfEvent();
136  if (nbEvents == 0) { return; }
137 
138  G4ParticleDefinition* particle = fPrimary->GetParticleGun()
139  ->GetParticleDefinition();
140  G4String partName = particle->GetParticleName();
141  G4double eprimary = fPrimary->GetParticleGun()->GetParticleEnergy();
142 
143  G4cout << "\n ======================== run summary ======================";
144  G4cout << "\n The run was " << nbEvents << " " << partName << " of "
145  << G4BestUnit(eprimary,"Energy");
146  G4cout << "\n ===========================================================\n";
147  G4cout << G4endl;
148 
149  G4int prec = 4, wid = prec + 2;
150  G4int dfprec = G4cout.precision(prec);
151 
152  //particle count
153  //
154  G4cout << " Nb of generated particles: \n" << G4endl;
155 
156  std::map<G4String,G4int>::iterator it;
157  for (it = fParticleCount.begin(); it != fParticleCount.end(); it++) {
158  G4String name = it->first;
159  G4int count = it->second;
160  G4double eMean = fEmean[name]/count;
161  G4double eMin = fEmin[name], eMax = fEmax[name];
162 
163  G4cout << " " << std::setw(13) << name << ": " << std::setw(7) << count
164  << " Emean = " << std::setw(wid) << G4BestUnit(eMean, "Energy")
165  << "\t( " << G4BestUnit(eMin, "Energy")
166  << " --> " << G4BestUnit(eMax, "Energy")
167  << ")" << G4endl;
168  }
169 
170  //energy momentum balance
171  //
172 
173  if (fDecayCount > 0) {
174  G4double Ebmean = fEkinTot[0]/fDecayCount;
175  G4double Pbmean = fPbalance[0]/fDecayCount;
176 
177  G4cout << "\n Ekin Total (Q): mean = "
178  << std::setw(wid) << G4BestUnit(Ebmean, "Energy")
179  << "\t( " << G4BestUnit(fEkinTot[1], "Energy")
180  << " --> " << G4BestUnit(fEkinTot[2], "Energy")
181  << ")" << G4endl;
182 
183  G4cout << "\n Momentum balance (excluding gamma desexcitation): mean = "
184  << std::setw(wid) << G4BestUnit(Pbmean, "Energy")
185  << "\t( " << G4BestUnit(fPbalance[1], "Energy")
186  << " --> " << G4BestUnit(fPbalance[2], "Energy")
187  << ")" << G4endl;
188  }
189 
190  //total time of life
191  //
192  if (fTimeCount > 0) {
193  G4double Tmean = fEventTime[0]/fTimeCount;
194  G4double halfLife = Tmean*std::log(2.);
195 
196  G4cout << "\n Total time of life : mean = "
197  << std::setw(wid) << G4BestUnit(Tmean, "Time")
198  << " half-life = "
199  << std::setw(wid) << G4BestUnit(halfLife, "Time")
200  << " ( " << G4BestUnit(fEventTime[1], "Time")
201  << " --> " << G4BestUnit(fEventTime[2], "Time")
202  << ")" << G4endl;
203  }
204 
205  //activity of primary ion
206  //
207  G4double pTimeMean = fPrimaryTime/nbEvents;
208  G4double molMass = particle->GetAtomicMass()*g/mole;
209  G4double nAtoms_perUnitOfMass = Avogadro/molMass;
210  G4double Activity_perUnitOfMass = 0.0;
211  if (pTimeMean > 0.0)
212  { Activity_perUnitOfMass = nAtoms_perUnitOfMass/pTimeMean; }
213 
214  G4cout << "\n Activity of " << partName << " = "
215  << std::setw(wid) << Activity_perUnitOfMass*g/becquerel
216  << " Bq/g (" << Activity_perUnitOfMass*g/curie
217  << " Ci/g) \n"
218  << G4endl;
219 
220  // remove all contents in fParticleCount
221  //
222  fParticleCount.clear();
223  fEmean.clear(); fEmin.clear(); fEmax.clear();
224 
225  // restore default precision
226  //
227  G4cout.precision(dfprec);
228 
229  //normalize and save histograms
230  //
231  G4AnalysisManager* analysisManager = G4AnalysisManager::Instance();
232  G4double factor = 100./nbEvents;
233  analysisManager->ScaleH1(1,factor);
234  analysisManager->ScaleH1(2,factor);
235  analysisManager->ScaleH1(3,factor);
236  analysisManager->ScaleH1(4,factor);
237  analysisManager->ScaleH1(5,factor);
238  if ( analysisManager->IsActive() ) {
239  analysisManager->Write();
240  analysisManager->CloseFile();
241  }
242 }
243 
244 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
Generated on Sat May 25 2013 14:31:59 for Geant4 by   doxygen 1.8.4