| 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 polarisation/Pol01/src/RunAction.cc << 26 // $Id: RunAction.cc,v 1.3 2006/11/17 11:44:46 vnivanch Exp $
27 /// \brief Implementation of the RunAction cla << 27 // GEANT4 tag $Name: geant4-08-02 $
28 // << 28 //
29 // <<
30 //....oooOO0OOooo........oooOO0OOooo........oo 29 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
31 //....oooOO0OOooo........oooOO0OOooo........oo 30 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
32 31
33 #include "RunAction.hh" 32 #include "RunAction.hh"
34 33
35 #include "DetectorConstruction.hh" 34 #include "DetectorConstruction.hh"
36 #include "PrimaryGeneratorAction.hh" 35 #include "PrimaryGeneratorAction.hh"
37 #include "ProcessesCount.hh" << 36 #include "HistoManager.hh"
38 37
39 #include "G4AccumulableManager.hh" <<
40 #include "G4Electron.hh" <<
41 #include "G4EmCalculator.hh" <<
42 #include "G4Gamma.hh" <<
43 #include "G4ParticleDefinition.hh" <<
44 #include "G4PhysicalConstants.hh" <<
45 #include "G4Positron.hh" <<
46 #include "G4Run.hh" 38 #include "G4Run.hh"
47 #include "G4RunManager.hh" 39 #include "G4RunManager.hh"
48 #include "G4SystemOfUnits.hh" <<
49 #include "G4UnitsTable.hh" 40 #include "G4UnitsTable.hh"
>> 41 #include "G4EmCalculator.hh"
50 42
>> 43 #include "Randomize.hh"
51 #include <iomanip> 44 #include <iomanip>
52 45
53 //....oooOO0OOooo........oooOO0OOooo........oo 46 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
54 47
55 RunAction::RunAction(DetectorConstruction* det << 48 RunAction::RunAction(DetectorConstruction* det, PrimaryGeneratorAction* prim,
56 : fDetector(det), fPrimary(prim), fAnalysisM << 49 HistoManager* histo)
>> 50 : detector(det), primary(prim), ProcCounter(0), histoManager(histo)
57 { 51 {
58 fGamma = G4Gamma::Gamma(); << 52 totalEventCount=0;
59 fElectron = G4Electron::Electron(); << 53 }
60 fPositron = G4Positron::Positron(); <<
61 <<
62 auto accumulableManager = G4AccumulableManag <<
63 auto fPhotonStats = new ParticleStatistics(" <<
64 auto fElectronStats = new ParticleStatistics <<
65 auto fPositronStats = new ParticleStatistics <<
66 auto fProcCounter = new ProcessesCount("Proc <<
67 <<
68 accumulableManager->Register(fPhotonStats); <<
69 accumulableManager->Register(fElectronStats) <<
70 accumulableManager->Register(fPositronStats) <<
71 <<
72 accumulableManager->Register(fProcCounter); <<
73 <<
74 BookHisto(); <<
75 } <<
76 54
77 //....oooOO0OOooo........oooOO0OOooo........oo 55 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
78 56
79 RunAction::~RunAction() {} << 57 RunAction::~RunAction()
>> 58 { }
80 59
81 //....oooOO0OOooo........oooOO0OOooo........oo 60 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
82 61
83 void RunAction::BeginOfRunAction(const G4Run* 62 void RunAction::BeginOfRunAction(const G4Run* aRun)
84 { << 63 {
85 G4cout << "### Run " << aRun->GetRunID() << 64 G4cout << "### Run " << aRun->GetRunID() << " start." << G4endl;
86 << 65
87 auto accumulableManager = G4AccumulableManag <<
88 accumulableManager->Reset(); <<
89 <<
90 // save Rndm status 66 // save Rndm status
91 // G4RunManager::GetRunManager()->SetRandom 67 // G4RunManager::GetRunManager()->SetRandomNumberStore(false);
92 // CLHEP::HepRandom::showEngineStatus(); 68 // CLHEP::HepRandom::showEngineStatus();
93 << 69
94 fTotalEventCount = 0; << 70 if (ProcCounter) delete ProcCounter;
95 << 71 ProcCounter = new ProcessesCount;
96 // Open file histogram file << 72 totalEventCount = 0;
97 fAnalysisManager->OpenFile(); << 73 photonStats.Clear();
98 } << 74 electronStats.Clear();
99 << 75 positronStats.Clear();
100 //....oooOO0OOooo........oooOO0OOooo........oo << 76
101 << 77 histoManager->book();
102 void RunAction::FillData(const G4ParticleDefin <<
103 G4double costheta, G4 <<
104 { <<
105 auto accManager = G4AccumulableManager::Inst <<
106 G4int id = -1; <<
107 if (particle == fGamma) { <<
108 dynamic_cast<ParticleStatistics*>(accManag <<
109 ->FillData(kinEnergy, costheta, longitud <<
110 if (fAnalysisManager) { <<
111 id = 1; <<
112 } <<
113 } <<
114 else if (particle == fElectron) { <<
115 dynamic_cast<ParticleStatistics*>(accManag <<
116 ->FillData(kinEnergy, costheta, longitud <<
117 if (fAnalysisManager) { <<
118 id = 5; <<
119 } <<
120 } <<
121 else if (particle == fPositron) { <<
122 dynamic_cast<ParticleStatistics*>(accManag <<
123 ->FillData(kinEnergy, costheta, longitud <<
124 if (fAnalysisManager) { <<
125 id = 9; <<
126 } <<
127 } <<
128 if (id > 0) { <<
129 fAnalysisManager->FillH1(id, kinEnergy, 1. <<
130 fAnalysisManager->FillH1(id + 1, costheta, <<
131 fAnalysisManager->FillH1(id + 2, phi, 1.0) <<
132 fAnalysisManager->FillH1(id + 3, longitudi <<
133 } <<
134 } 78 }
135 79
136 //....oooOO0OOooo........oooOO0OOooo........oo 80 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
137 << 81 void RunAction::FillData(const G4String & particleName,
138 void RunAction::BookHisto() << 82 G4double kinEnergy, G4double costheta,
139 { << 83 G4double /* phi*/,
140 // Always creating analysis manager << 84 G4double longitudinalPolarization)
141 fAnalysisManager = G4AnalysisManager::Instan << 85 {
142 fAnalysisManager->SetDefaultFileType("root") << 86 if (particleName=="gamma")
143 fAnalysisManager->SetActivation(true); << 87 photonStats.FillData(kinEnergy, costheta, longitudinalPolarization);
144 fAnalysisManager->SetVerboseLevel(1); << 88 else if (particleName=="e-")
145 << 89 electronStats.FillData(kinEnergy, costheta, longitudinalPolarization);
146 // Open file histogram file << 90 else if (particleName=="e+")
147 fAnalysisManager->SetFileName("pol01"); << 91 positronStats.FillData(kinEnergy, costheta, longitudinalPolarization);
148 << 92 }
149 fAnalysisManager->SetFirstHistoId(1); << 93
150 << 94 void RunAction::CountProcesses(G4String procName)
151 // Creating an 1-dimensional histograms in t << 95 {
152 const G4String id[] = {"h1", "h2", "h3", "h4 << 96 // is the process already counted ?
153 const G4String title[] = { << 97 // *AS* change to std::map?!
154 "Gamma Energy distribution", // 1 << 98 size_t nbProc = ProcCounter->size();
155 "Gamma Cos(Theta) distribution", // 2 << 99 size_t i = 0;
156 "Gamma Phi angular distribution", // 3 << 100 while ((i<nbProc)&&((*ProcCounter)[i]->GetName()!=procName)) i++;
157 "Gamma longitudinal Polarization", // 4 << 101 if (i == nbProc) ProcCounter->push_back( new OneProcessCount(procName));
158 "Electron Energy distribution", // 5 << 102
159 "Electron Cos(Theta) distribution", // 6 << 103 (*ProcCounter)[i]->Count();
160 "Electron Phi angular distribution", // 7 <<
161 "Electron longitudinal Polarization", // <<
162 "Positron Energy distribution", // 9 <<
163 "Positron Cos(Theta) distribution", // 10 <<
164 "Positron Phi angular distribution", // 1 <<
165 "Positron longitudinal Polarization" // 1 <<
166 }; <<
167 G4double vmin, vmax; <<
168 G4int nbins = 120; <<
169 for (int i = 0; i < 12; ++i) { <<
170 G4int j = i - i / 4 * 4; <<
171 if (0 == j) { <<
172 vmin = 0.; <<
173 vmax = 12. * MeV; <<
174 } <<
175 else if (1 == j) { <<
176 vmin = -1.; <<
177 vmax = 1.; <<
178 } <<
179 else if (2 == j) { <<
180 vmin = 0.; <<
181 vmax = pi; <<
182 } <<
183 else { <<
184 vmin = -1.5; <<
185 vmax = 1.5; <<
186 } <<
187 G4int ih = fAnalysisManager->CreateH1(id[i <<
188 fAnalysisManager->SetH1Activation(ih, fals <<
189 } <<
190 } <<
191 <<
192 //....oooOO0OOooo........oooOO0OOooo........oo <<
193 <<
194 void RunAction::SaveHisto(G4int nevents) <<
195 { <<
196 if (fAnalysisManager) { <<
197 if (IsMaster()) { <<
198 G4double norm = 1.0 / G4double(nevents); <<
199 for (int i = 0; i < 12; ++i) { <<
200 fAnalysisManager->ScaleH1(i, norm); <<
201 } <<
202 } <<
203 fAnalysisManager->Write(); <<
204 fAnalysisManager->CloseFile(); <<
205 } <<
206 } <<
207 <<
208 //....oooOO0OOooo........oooOO0OOooo........oo <<
209 <<
210 void RunAction::CountProcesses(G4String& procN <<
211 { <<
212 auto accManager = G4AccumulableManager::Inst <<
213 dynamic_cast<ProcessesCount*>(accManager->Ge <<
214 } 104 }
215 105
216 //....oooOO0OOooo........oooOO0OOooo........oo 106 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
217 107
218 void RunAction::EndOfRunAction(const G4Run* aR 108 void RunAction::EndOfRunAction(const G4Run* aRun)
219 { 109 {
220 // Total number of events in run (all threa << 110 G4int NbOfEvents = aRun->GetNumberOfEvent();
221 G4int NbOfEvents = aRun->GetNumberOfEventToB <<
222 // G4int NbOfEvents = aRun->GetNumberOfEven <<
223 <<
224 if (NbOfEvents == 0) return; 111 if (NbOfEvents == 0) return;
225 << 112
226 G4int prec = G4cout.precision(5); << 113 G4int prec = G4cout.precision(5);
227 << 114
228 G4Material* material = fDetector->GetMateria << 115 G4Material* material = detector->GetMaterial();
229 G4double density = material->GetDensity(); 116 G4double density = material->GetDensity();
>> 117
>> 118 G4ParticleDefinition* particle =
>> 119 primary->GetParticleGun()->GetParticleDefinition();
>> 120 G4String Particle = particle->GetParticleName();
>> 121 G4double energy = primary->GetParticleGun()->GetParticleEnergy();
>> 122 G4cout << "\n The run consists of " << NbOfEvents << " "<< Particle << " of "
>> 123 << G4BestUnit(energy,"Energy") << " through "
>> 124 << G4BestUnit(detector->GetBoxSizeZ(),"Length") << " of "
>> 125 << material->GetName() << " (density: "
>> 126 << G4BestUnit(density,"Volumic Mass") << ")" << G4endl;
>> 127
>> 128 //frequency of processes
>> 129 G4cout << "\n Process calls frequency --->\n";
>> 130 for (size_t i=0; i< ProcCounter->size();i++) {
>> 131 G4String procName = (*ProcCounter)[i]->GetName();
>> 132 G4int count = (*ProcCounter)[i]->GetCounter();
>> 133 G4cout << "\t" << procName << " = " << count<<"\n";
>> 134 }
>> 135
>> 136 if (totalEventCount == 0) return;
>> 137
>> 138 G4cout<<" Gamma: \n";
>> 139 photonStats.PrintResults(totalEventCount);
>> 140 G4cout<<" Electron: \n";
>> 141 electronStats.PrintResults(totalEventCount);
>> 142 G4cout<<" Positron: \n";
>> 143 positronStats.PrintResults(totalEventCount);
>> 144
>> 145 //cross check from G4EmCalculator
>> 146 // G4cout << "\n Verification from G4EmCalculator. \n";
>> 147 // G4EmCalculator emCal;
230 148
231 if (fPrimary != nullptr) { << 149 //restore default format
232 G4ParticleDefinition* particle = fPrimary- << 150 G4cout.precision(prec);
233 G4String Particle = particle->GetParticleN <<
234 G4double energy = fPrimary->GetParticleGun <<
235 G4cout << "\n The run consists of " << fTo <<
236 << G4BestUnit(energy, "Energy") << <<
237 << G4BestUnit(fDetector->GetBoxSize <<
238 << " (density: " << G4BestUnit(dens <<
239 } <<
240 // cross check from G4EmCalculator <<
241 // G4cout << "\n Verification from G4EmCal <<
242 // G4EmCalculator emCal; <<
243 <<
244 auto accManager = G4AccumulableManager::Inst <<
245 accManager->Merge(); <<
246 <<
247 if (IsMaster()) { <<
248 // frequency of processes <<
249 G4cout << "\n Process calls frequency ---> <<
250 dynamic_cast<ProcessesCount*>(accManager-> <<
251 G4cout << " Gamma: \n"; <<
252 dynamic_cast<ParticleStatistics*>(accManag <<
253 ->PrintResults(NbOfEvents); <<
254 G4cout << " Electron: \n"; <<
255 dynamic_cast<ParticleStatistics*>(accManag <<
256 ->PrintResults(NbOfEvents); <<
257 G4cout << " Positron: \n"; <<
258 dynamic_cast<ParticleStatistics*>(accManag <<
259 ->PrintResults(NbOfEvents); <<
260 G4cout << G4endl; <<
261 } <<
262 <<
263 // restore default format <<
264 G4cout.precision(prec); <<
265 151
266 // write out histograms << 152 // write out histograms
267 SaveHisto(NbOfEvents); << 153 histoManager->save();
268 154
269 if (IsMaster()) { << 155 // show Rndm status
270 // show Rndm status << 156 CLHEP::HepRandom::showEngineStatus();
271 CLHEP::HepRandom::showEngineStatus(); <<
272 } <<
273 } 157 }
274 158
275 //....oooOO0OOooo........oooOO0OOooo........oo 159 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
276 160
277 void RunAction::EventFinished() 161 void RunAction::EventFinished()
278 { 162 {
279 auto accManager = G4AccumulableManager::Inst << 163 ++totalEventCount;
280 << 164 photonStats.EventFinished();
281 ++fTotalEventCount; << 165 electronStats.EventFinished();
282 << 166 positronStats.EventFinished();
283 dynamic_cast<ParticleStatistics*>(accManager <<
284 dynamic_cast<ParticleStatistics*>(accManager <<
285 dynamic_cast<ParticleStatistics*>(accManager <<
286 } 167 }
287 168
288 //....oooOO0OOooo........oooOO0OOooo........oo 169 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
289 170
290 RunAction::ParticleStatistics::ParticleStatist << 171 RunAction::ParticleStatistics::ParticleStatistics()
291 : G4VAccumulable(name), << 172 : currentNumber(0),
292 fCurrentNumber(0), << 173 totalNumber(0), totalNumber2(0),
293 fTotalNumber(0), << 174 sumEnergy(0), sumEnergy2(0),
294 fTotalNumber2(0), << 175 sumPolarization(0), sumPolarization2(0),
295 fSumEnergy(0), << 176 sumCosTheta(0), sumCosTheta2(0)
296 fSumEnergy2(0), <<
297 fSumPolarization(0), <<
298 fSumPolarization2(0), <<
299 fSumCosTheta(0), <<
300 fSumCosTheta2(0) <<
301 {} 177 {}
302 178
303 //....oooOO0OOooo........oooOO0OOooo........oo 179 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
304 180
305 RunAction::ParticleStatistics::~ParticleStatis << 181 RunAction::ParticleStatistics::~ParticleStatistics()
>> 182 {}
306 183
307 //....oooOO0OOooo........oooOO0OOooo........oo 184 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
308 185
309 void RunAction::ParticleStatistics::EventFinis 186 void RunAction::ParticleStatistics::EventFinished()
310 { 187 {
311 fTotalNumber += fCurrentNumber; << 188 totalNumber+=currentNumber;
312 fTotalNumber2 += fCurrentNumber * fCurrentNu << 189 totalNumber2+=currentNumber*currentNumber;
313 fCurrentNumber = 0; << 190 currentNumber=0;
314 } 191 }
315 //....oooOO0OOooo........oooOO0OOooo........oo 192 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
316 193
317 void RunAction::ParticleStatistics::FillData(G << 194 void RunAction::ParticleStatistics:: FillData(G4double kinEnergy,
318 G << 195 G4double costheta,
319 { << 196 G4double longitudinalPolarization)
320 ++fCurrentNumber; << 197 {
321 fSumEnergy += kinEnergy; << 198 ++currentNumber;
322 fSumEnergy2 += kinEnergy * kinEnergy; << 199 sumEnergy+=kinEnergy;
323 fSumPolarization += longitudinalPolarization << 200 sumEnergy2+=kinEnergy*kinEnergy;
324 fSumPolarization2 += longitudinalPolarizatio << 201 sumPolarization+=longitudinalPolarization;
325 fSumCosTheta += costheta; << 202 sumPolarization2+=longitudinalPolarization*longitudinalPolarization;
326 fSumCosTheta2 += costheta * costheta; << 203 sumCosTheta+=costheta;
>> 204 sumCosTheta2+=costheta*costheta;
327 } 205 }
328 206
329 //....oooOO0OOooo........oooOO0OOooo........oo 207 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
330 208
331 void RunAction::ParticleStatistics::PrintResul 209 void RunAction::ParticleStatistics::PrintResults(G4int totalNumberOfEvents)
332 { 210 {
333 G4cout << "Mean Number per Event :" << G4dou << 211 G4cout<<"Mean Number per Event :"
334 << "\n"; << 212 <<G4double(totalNumber)/G4double(totalNumberOfEvents)<<"\n";
335 if (fTotalNumber == 0) fTotalNumber = 1; << 213 if (totalNumber==0) totalNumber=1;
336 G4double energyMean = fSumEnergy / fTotalNum << 214 G4double energyMean=sumEnergy/totalNumber;
337 G4double energyRms = std::sqrt(fSumEnergy2 / << 215 G4double energyRms=std::sqrt(sumEnergy2/totalNumber-energyMean*energyMean);
338 G4cout << "Mean Energy :" << G4BestUnit(ener << 216 G4cout<<"Mean Energy :"<< G4BestUnit(energyMean,"Energy")
339 << G4BestUnit(energyRms, "Energy") << << 217 <<" +- "<<G4BestUnit(energyRms,"Energy")<<"\n";
340 G4double polarizationMean = fSumPolarization << 218 G4double polarizationMean=sumPolarization/totalNumber;
341 G4double polarizationRms = << 219 G4double polarizationRms=
342 std::sqrt(fSumPolarization2 / fTotalNumber << 220 std::sqrt(sumPolarization2/totalNumber-polarizationMean*polarizationMean);
343 G4cout << "Mean Polarization :" << polarizat << 221 G4cout<<"Mean Polarization :"<< polarizationMean
344 } << 222 <<" +- "<<polarizationRms<<"\n";
345 <<
346 //....oooOO0OOooo........oooOO0OOooo........oo <<
347 <<
348 void RunAction::ParticleStatistics::Reset() <<
349 { <<
350 fCurrentNumber = 0; <<
351 fTotalNumber = fTotalNumber2 = 0; <<
352 fSumEnergy = fSumEnergy2 = 0; <<
353 fSumPolarization = fSumPolarization2 = 0; <<
354 fSumCosTheta = fSumCosTheta2 = 0; <<
355 } 223 }
356 224
357 //....oooOO0OOooo........oooOO0OOooo........oo 225 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
358 226
359 void RunAction::ParticleStatistics::Merge(cons << 227 void RunAction::ParticleStatistics::Clear()
360 { << 228 {
361 auto rstat = dynamic_cast<const RunAction::P << 229 currentNumber=0;
362 << 230 totalNumber=totalNumber2=0;
363 fCurrentNumber += rstat.fCurrentNumber; << 231 sumEnergy=sumEnergy2=0;
364 fTotalNumber += rstat.fTotalNumber; << 232 sumPolarization=sumPolarization2=0;
365 fTotalNumber2 += rstat.fTotalNumber2; << 233 sumCosTheta=sumCosTheta2=0;
366 fSumEnergy += rstat.fSumEnergy; <<
367 fSumEnergy2 += rstat.fSumEnergy2; <<
368 fSumPolarization += rstat.fSumPolarization; <<
369 fSumPolarization2 += rstat.fSumPolarization2 <<
370 fSumCosTheta += rstat.fSumCosTheta; <<
371 fSumCosTheta2 += rstat.fSumCosTheta2; <<
372 } 234 }
373 235
374 //....oooOO0OOooo........oooOO0OOooo........oo 236 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
375 237