| Geant4 Cross Reference |
1 // 1 //
2 // ******************************************* 2 // ********************************************************************
3 // * License and Disclaimer 3 // * License and Disclaimer *
4 // * 4 // * *
5 // * The Geant4 software is copyright of th 5 // * The Geant4 software is copyright of the Copyright Holders of *
6 // * the Geant4 Collaboration. It is provided 6 // * the Geant4 Collaboration. It is provided under the terms and *
7 // * conditions of the Geant4 Software License 7 // * conditions of the Geant4 Software License, included in the file *
8 // * LICENSE and available at http://cern.ch/ 8 // * LICENSE and available at http://cern.ch/geant4/license . These *
9 // * include a list of copyright holders. 9 // * include a list of copyright holders. *
10 // * 10 // * *
11 // * Neither the authors of this software syst 11 // * Neither the authors of this software system, nor their employing *
12 // * institutes,nor the agencies providing fin 12 // * institutes,nor the agencies providing financial support for this *
13 // * work make any representation or warran 13 // * work make any representation or warranty, express or implied, *
14 // * regarding this software system or assum 14 // * regarding this software system or assume any liability for its *
15 // * use. Please see the license in the file 15 // * use. Please see the license in the file LICENSE and URL above *
16 // * for the full disclaimer and the limitatio 16 // * for the full disclaimer and the limitation of liability. *
17 // * 17 // * *
18 // * This code implementation is the result 18 // * This code implementation is the result of the scientific and *
19 // * technical work of the GEANT4 collaboratio 19 // * technical work of the GEANT4 collaboration. *
20 // * By using, copying, modifying or distri 20 // * By using, copying, modifying or distributing the software (or *
21 // * any work based on the software) you ag 21 // * any work based on the software) you agree to acknowledge its *
22 // * use in resulting scientific publicati 22 // * use in resulting scientific publications, and indicate your *
23 // * acceptance of all terms of the Geant4 Sof 23 // * acceptance of all terms of the Geant4 Software license. *
24 // ******************************************* 24 // ********************************************************************
25 // 25 //
26 /// \file RunAction.cc <<
27 /// \brief Implementation of the RunAction cla <<
28 // 26 //
29 // << 27 // $Id: RunAction.cc,v 1.2 2010-10-11 14:31:39 maire Exp $
30 //....oooOO0OOooo........oooOO0OOooo........oo << 28 // GEANT4 tag $Name: not supported by cvs2svn $
>> 29 //
31 //....oooOO0OOooo........oooOO0OOooo........oo 30 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
>> 31 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
32 32
33 #include "RunAction.hh" 33 #include "RunAction.hh"
34 <<
35 #include "HistoManager.hh" 34 #include "HistoManager.hh"
36 #include "PrimaryGeneratorAction.hh" 35 #include "PrimaryGeneratorAction.hh"
37 #include "Run.hh" <<
38 36
39 #include "G4PhysicalConstants.hh" <<
40 #include "G4Run.hh" 37 #include "G4Run.hh"
41 #include "G4SystemOfUnits.hh" << 38 #include "G4RunManager.hh"
42 #include "G4UnitsTable.hh" 39 #include "G4UnitsTable.hh"
43 <<
44 #include <iomanip> 40 #include <iomanip>
45 41
46 //....oooOO0OOooo........oooOO0OOooo........oo 42 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
47 43
48 RunAction::RunAction(PrimaryGeneratorAction* k << 44 RunAction::RunAction(HistoManager* histo, PrimaryGeneratorAction* kin)
49 { << 45 :histoManager(histo), primary(kin)
50 fHistoManager = new HistoManager(); << 46 { }
51 } <<
52 47
53 //....oooOO0OOooo........oooOO0OOooo........oo 48 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
54 49
55 RunAction::~RunAction() 50 RunAction::~RunAction()
56 { << 51 { }
57 delete fHistoManager; << 52
>> 53 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
>> 54
>> 55 void RunAction::BeginOfRunAction(const G4Run*)
>> 56 {
>> 57 //initialize arrays
>> 58 //
>> 59 decayCount = 0;
>> 60 for (G4int i=0; i<3; i++) Ebalance[i] = Pbalance[i] = EventTime[i] = 0. ;
>> 61
>> 62 //histograms
>> 63 //
>> 64 histoManager->book();
>> 65
>> 66 //inform the runManager to save random number seed
>> 67 //
>> 68 G4RunManager::GetRunManager()->SetRandomNumberStore(false);
58 } 69 }
>> 70
59 //....oooOO0OOooo........oooOO0OOooo........oo 71 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
60 72
61 G4Run* RunAction::GenerateRun() << 73 void RunAction::ParticleCount(G4String name, G4double Ekin)
62 { 74 {
63 fRun = new Run(); << 75 particleCount[name]++;
64 return fRun; << 76 Emean[name] += Ekin;
>> 77 if (particleCount[name] == 1) Emin[name] = Emax[name] = Ekin;
>> 78 if (Ekin < Emin[name]) Emin[name] = Ekin;
>> 79 if (Ekin > Emax[name]) Emax[name] = Ekin;
65 } 80 }
>> 81
66 //....oooOO0OOooo........oooOO0OOooo........oo 82 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
67 83
68 void RunAction::BeginOfRunAction(const G4Run*) << 84 void RunAction::Balance(G4double Ebal, G4double Pbal)
69 { 85 {
70 // keep run condition << 86 decayCount++;
71 if (fPrimary) { << 87 Ebalance[0] += Ebal;
72 G4ParticleDefinition* particle = fPrimary- << 88 if (decayCount == 1) Ebalance[1] = Ebalance[2] = Ebal;
73 G4double energy = fPrimary->GetParticleGun << 89 if (Ebal < Ebalance[1]) Ebalance[1] = Ebal;
74 fRun->SetPrimary(particle, energy); << 90 if (Ebal > Ebalance[2]) Ebalance[2] = Ebal;
75 } << 91
76 << 92 Pbalance[0] += Pbal;
77 // histograms << 93 if (decayCount == 1) Pbalance[1] = Pbalance[2] = Pbal;
78 // << 94 if (Pbal < Pbalance[1]) Pbalance[1] = Pbal;
79 G4AnalysisManager* analysisManager = G4Analy << 95 if (Pbal > Pbalance[2]) Pbalance[2] = Pbal;
80 if (analysisManager->IsActive()) { <<
81 analysisManager->OpenFile(); <<
82 } <<
83 } 96 }
84 97
85 //....oooOO0OOooo........oooOO0OOooo........oo 98 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
86 99
87 void RunAction::EndOfRunAction(const G4Run*) << 100 void RunAction::EventTiming(G4double time)
88 { << 101 {
89 if (isMaster) fRun->EndOfRun(); << 102 EventTime[0] += time;
>> 103 if (decayCount == 1) EventTime[1] = EventTime[2] = time;
>> 104 if (time < EventTime[1]) EventTime[1] = time;
>> 105 if (time > EventTime[2]) EventTime[2] = time;
>> 106 }
>> 107 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
90 108
91 // save histograms << 109 void RunAction::EndOfRunAction(const G4Run* run)
92 // << 110 {
93 G4AnalysisManager* analysisManager = G4Analy << 111 G4int nbEvents = run->GetNumberOfEvent();
94 if (analysisManager->IsActive()) { << 112 if (nbEvents == 0) return;
95 analysisManager->Write(); << 113
96 analysisManager->CloseFile(); << 114 G4ParticleDefinition* particle = primary->GetParticleGun()
97 } << 115 ->GetParticleDefinition();
>> 116 G4String partName = particle->GetParticleName();
>> 117 G4double eprimary = primary->GetParticleGun()->GetParticleEnergy();
>> 118
>> 119 G4cout << "\n ======================== run summary ======================";
>> 120 G4cout << "\n The run was " << nbEvents << " " << partName << " of "
>> 121 << G4BestUnit(eprimary,"Energy");
>> 122 G4cout << "\n ===========================================================\n";
>> 123 G4cout << G4endl;
>> 124
>> 125 G4int prec = 4, wid = prec + 2;
>> 126 G4int dfprec = G4cout.precision(prec);
>> 127
>> 128 //particle count
>> 129 //
>> 130 G4cout << " Nb of generated particles: \n" << G4endl;
>> 131
>> 132 std::map<G4String,G4int>::iterator it;
>> 133 for (it = particleCount.begin(); it != particleCount.end(); it++) {
>> 134 G4String name = it->first;
>> 135 G4int count = it->second;
>> 136 G4double eMean = Emean[name]/count;
>> 137 G4double eMin = Emin[name], eMax = Emax[name];
>> 138
>> 139 G4cout << " " << std::setw(13) << name << ": " << std::setw(7) << count
>> 140 << " Emean = " << std::setw(wid) << G4BestUnit(eMean, "Energy")
>> 141 << "\t( " << G4BestUnit(eMin, "Energy")
>> 142 << " --> " << G4BestUnit(eMax, "Energy")
>> 143 << ")" << G4endl;
>> 144 }
>> 145
>> 146 //energy momentum balance
>> 147 //
>> 148 G4double Ebmean = Ebalance[0]/decayCount;
>> 149 G4double Pbmean = Pbalance[0]/decayCount;
>> 150
>> 151 G4cout << "\n Energy and momentum balance : final state - initial state"
>> 152 << "\n (excluding gamma desexcitation from momentum balance) : \n"
>> 153 << G4endl;
>> 154
>> 155 G4cout
>> 156 << " Energy: mean = " << std::setw(wid) << G4BestUnit(Ebmean, "Energy")
>> 157 << "\t( " << G4BestUnit(Ebalance[1], "Energy")
>> 158 << " --> " << G4BestUnit(Ebalance[2], "Energy")
>> 159 << ")" << G4endl;
>> 160
>> 161 G4cout
>> 162 << " Momentum: mean = " << std::setw(wid) << G4BestUnit(Pbmean, "Energy")
>> 163 << "\t( " << G4BestUnit(Pbalance[1], "Energy")
>> 164 << " --> " << G4BestUnit(Pbalance[2], "Energy")
>> 165 << ")" << G4endl;
>> 166
>> 167 //time of life
>> 168 //
>> 169 G4double Tmean = EventTime[0]/nbEvents;
>> 170 G4double halfLife = Tmean*std::log(2.);
>> 171
>> 172 G4cout << "\n Time of life : mean = "
>> 173 << std::setw(wid) << G4BestUnit(Tmean, "Time")
>> 174 << " half-life = "
>> 175 << std::setw(wid) << G4BestUnit(halfLife, "Time")
>> 176 << " ( " << G4BestUnit(EventTime[1], "Time")
>> 177 << " --> " << G4BestUnit(EventTime[2], "Time")
>> 178 << ")" << G4endl;
>> 179
>> 180 //activity
>> 181 //
>> 182 G4double molMass = particle->GetAtomicMass()*g/mole;
>> 183 G4double nAtoms = Avogadro/molMass;
>> 184 G4double ActivPerAtom = 1./Tmean;
>> 185 G4double ActivPerMass = ActivPerAtom*nAtoms;
>> 186
>> 187 G4cout << "\n Activity = "
>> 188 << std::setw(wid) << ActivPerMass*g/becquerel
>> 189 << " Bq/g (" << ActivPerMass*g/curie
>> 190 << " Ci/g) \n"
>> 191 << G4endl;
>> 192
>> 193 //normalise histo 9
>> 194 G4double binW = histoManager->GetBinWidth(9);
>> 195 G4double factor = (nAtoms*gram)/(binW*nbEvents*becquerel);
>> 196 histoManager->Normalize(9,factor);
>> 197
>> 198 // remove all contents in particleCount
>> 199 //
>> 200 particleCount.clear();
>> 201 Emean.clear(); Emin.clear(); Emax.clear();
>> 202
>> 203 // restore default precision
>> 204 //
>> 205 G4cout.precision(dfprec);
>> 206
>> 207 //save histograms
>> 208 //
>> 209 histoManager->save();
98 } 210 }
99 211
100 //....oooOO0OOooo........oooOO0OOooo........oo 212 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
101 213