| Geant4 Cross Reference |
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 //
26 /// \file polarisation/Pol01/src/RunAction.cc
27 /// \brief Implementation of the RunAction class
28 //
29 //
30 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
31 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
32
33 #include "RunAction.hh"
34
35 #include "DetectorConstruction.hh"
36 #include "PrimaryGeneratorAction.hh"
37 #include "ProcessesCount.hh"
38
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"
47 #include "G4RunManager.hh"
48 #include "G4SystemOfUnits.hh"
49 #include "G4UnitsTable.hh"
50
51 #include <iomanip>
52
53 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
54
55 RunAction::RunAction(DetectorConstruction* det, PrimaryGeneratorAction* prim)
56 : fDetector(det), fPrimary(prim), fAnalysisManager(0), fTotalEventCount(0)
57 {
58 fGamma = G4Gamma::Gamma();
59 fElectron = G4Electron::Electron();
60 fPositron = G4Positron::Positron();
61
62 auto accumulableManager = G4AccumulableManager::Instance();
63 auto fPhotonStats = new ParticleStatistics("PhotonStats");
64 auto fElectronStats = new ParticleStatistics("ElectronStats");
65 auto fPositronStats = new ParticleStatistics("PositronStats");
66 auto fProcCounter = new ProcessesCount("ProcCounter");
67
68 accumulableManager->Register(fPhotonStats);
69 accumulableManager->Register(fElectronStats);
70 accumulableManager->Register(fPositronStats);
71
72 accumulableManager->Register(fProcCounter);
73
74 BookHisto();
75 }
76
77 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
78
79 RunAction::~RunAction() {}
80
81 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
82
83 void RunAction::BeginOfRunAction(const G4Run* aRun)
84 {
85 G4cout << "### Run " << aRun->GetRunID() << " start." << G4endl;
86
87 auto accumulableManager = G4AccumulableManager::Instance();
88 accumulableManager->Reset();
89
90 // save Rndm status
91 // G4RunManager::GetRunManager()->SetRandomNumberStore(false);
92 // CLHEP::HepRandom::showEngineStatus();
93
94 fTotalEventCount = 0;
95
96 // Open file histogram file
97 fAnalysisManager->OpenFile();
98 }
99
100 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
101
102 void RunAction::FillData(const G4ParticleDefinition* particle, G4double kinEnergy,
103 G4double costheta, G4double phi, G4double longitudinalPolarization)
104 {
105 auto accManager = G4AccumulableManager::Instance();
106 G4int id = -1;
107 if (particle == fGamma) {
108 dynamic_cast<ParticleStatistics*>(accManager->GetAccumulable("PhotonStats"))
109 ->FillData(kinEnergy, costheta, longitudinalPolarization);
110 if (fAnalysisManager) {
111 id = 1;
112 }
113 }
114 else if (particle == fElectron) {
115 dynamic_cast<ParticleStatistics*>(accManager->GetAccumulable("ElectronStats"))
116 ->FillData(kinEnergy, costheta, longitudinalPolarization);
117 if (fAnalysisManager) {
118 id = 5;
119 }
120 }
121 else if (particle == fPositron) {
122 dynamic_cast<ParticleStatistics*>(accManager->GetAccumulable("PositronStats"))
123 ->FillData(kinEnergy, costheta, longitudinalPolarization);
124 if (fAnalysisManager) {
125 id = 9;
126 }
127 }
128 if (id > 0) {
129 fAnalysisManager->FillH1(id, kinEnergy, 1.0);
130 fAnalysisManager->FillH1(id + 1, costheta, 1.0);
131 fAnalysisManager->FillH1(id + 2, phi, 1.0);
132 fAnalysisManager->FillH1(id + 3, longitudinalPolarization, 1.0);
133 }
134 }
135
136 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
137
138 void RunAction::BookHisto()
139 {
140 // Always creating analysis manager
141 fAnalysisManager = G4AnalysisManager::Instance();
142 fAnalysisManager->SetDefaultFileType("root");
143 fAnalysisManager->SetActivation(true);
144 fAnalysisManager->SetVerboseLevel(1);
145
146 // Open file histogram file
147 fAnalysisManager->SetFileName("pol01");
148
149 fAnalysisManager->SetFirstHistoId(1);
150
151 // Creating an 1-dimensional histograms in the root directory of the tree
152 const G4String id[] = {"h1", "h2", "h3", "h4", "h5", "h6", "h7", "h8", "h9", "h10", "h11", "h12"};
153 const G4String title[] = {
154 "Gamma Energy distribution", // 1
155 "Gamma Cos(Theta) distribution", // 2
156 "Gamma Phi angular distribution", // 3
157 "Gamma longitudinal Polarization", // 4
158 "Electron Energy distribution", // 5
159 "Electron Cos(Theta) distribution", // 6
160 "Electron Phi angular distribution", // 7
161 "Electron longitudinal Polarization", // 8
162 "Positron Energy distribution", // 9
163 "Positron Cos(Theta) distribution", // 10
164 "Positron Phi angular distribution", // 11
165 "Positron longitudinal Polarization" // 12
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], title[i], nbins, vmin, vmax);
188 fAnalysisManager->SetH1Activation(ih, false);
189 }
190 }
191
192 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
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........oooOO0OOooo........oooOO0OOooo......
209
210 void RunAction::CountProcesses(G4String& procName)
211 {
212 auto accManager = G4AccumulableManager::Instance();
213 dynamic_cast<ProcessesCount*>(accManager->GetAccumulable("ProcCounter"))->Count(procName);
214 }
215
216 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
217
218 void RunAction::EndOfRunAction(const G4Run* aRun)
219 {
220 // Total number of events in run (all threads)
221 G4int NbOfEvents = aRun->GetNumberOfEventToBeProcessed();
222 // G4int NbOfEvents = aRun->GetNumberOfEvent();
223
224 if (NbOfEvents == 0) return;
225
226 G4int prec = G4cout.precision(5);
227
228 G4Material* material = fDetector->GetMaterial();
229 G4double density = material->GetDensity();
230
231 if (fPrimary != nullptr) {
232 G4ParticleDefinition* particle = fPrimary->GetParticleGun()->GetParticleDefinition();
233 G4String Particle = particle->GetParticleName();
234 G4double energy = fPrimary->GetParticleGun()->GetParticleEnergy();
235 G4cout << "\n The run consists of " << fTotalEventCount << " " << Particle << " of "
236 << G4BestUnit(energy, "Energy") << " through "
237 << G4BestUnit(fDetector->GetBoxSizeZ(), "Length") << " of " << material->GetName()
238 << " (density: " << G4BestUnit(density, "Volumic Mass") << ")" << G4endl;
239 }
240 // cross check from G4EmCalculator
241 // G4cout << "\n Verification from G4EmCalculator. \n";
242 // G4EmCalculator emCal;
243
244 auto accManager = G4AccumulableManager::Instance();
245 accManager->Merge();
246
247 if (IsMaster()) {
248 // frequency of processes
249 G4cout << "\n Process calls frequency --->\n";
250 dynamic_cast<ProcessesCount*>(accManager->GetAccumulable("ProcCounter"))->Print();
251 G4cout << " Gamma: \n";
252 dynamic_cast<ParticleStatistics*>(accManager->GetAccumulable("PhotonStats"))
253 ->PrintResults(NbOfEvents);
254 G4cout << " Electron: \n";
255 dynamic_cast<ParticleStatistics*>(accManager->GetAccumulable("ElectronStats"))
256 ->PrintResults(NbOfEvents);
257 G4cout << " Positron: \n";
258 dynamic_cast<ParticleStatistics*>(accManager->GetAccumulable("PositronStats"))
259 ->PrintResults(NbOfEvents);
260 G4cout << G4endl;
261 }
262
263 // restore default format
264 G4cout.precision(prec);
265
266 // write out histograms
267 SaveHisto(NbOfEvents);
268
269 if (IsMaster()) {
270 // show Rndm status
271 CLHEP::HepRandom::showEngineStatus();
272 }
273 }
274
275 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
276
277 void RunAction::EventFinished()
278 {
279 auto accManager = G4AccumulableManager::Instance();
280
281 ++fTotalEventCount;
282
283 dynamic_cast<ParticleStatistics*>(accManager->GetAccumulable("PhotonStats"))->EventFinished();
284 dynamic_cast<ParticleStatistics*>(accManager->GetAccumulable("ElectronStats"))->EventFinished();
285 dynamic_cast<ParticleStatistics*>(accManager->GetAccumulable("PositronStats"))->EventFinished();
286 }
287
288 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
289
290 RunAction::ParticleStatistics::ParticleStatistics(const G4String& name)
291 : G4VAccumulable(name),
292 fCurrentNumber(0),
293 fTotalNumber(0),
294 fTotalNumber2(0),
295 fSumEnergy(0),
296 fSumEnergy2(0),
297 fSumPolarization(0),
298 fSumPolarization2(0),
299 fSumCosTheta(0),
300 fSumCosTheta2(0)
301 {}
302
303 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
304
305 RunAction::ParticleStatistics::~ParticleStatistics() {}
306
307 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
308
309 void RunAction::ParticleStatistics::EventFinished()
310 {
311 fTotalNumber += fCurrentNumber;
312 fTotalNumber2 += fCurrentNumber * fCurrentNumber;
313 fCurrentNumber = 0;
314 }
315 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
316
317 void RunAction::ParticleStatistics::FillData(G4double kinEnergy, G4double costheta,
318 G4double longitudinalPolarization)
319 {
320 ++fCurrentNumber;
321 fSumEnergy += kinEnergy;
322 fSumEnergy2 += kinEnergy * kinEnergy;
323 fSumPolarization += longitudinalPolarization;
324 fSumPolarization2 += longitudinalPolarization * longitudinalPolarization;
325 fSumCosTheta += costheta;
326 fSumCosTheta2 += costheta * costheta;
327 }
328
329 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
330
331 void RunAction::ParticleStatistics::PrintResults(G4int totalNumberOfEvents)
332 {
333 G4cout << "Mean Number per Event :" << G4double(fTotalNumber) / G4double(totalNumberOfEvents)
334 << "\n";
335 if (fTotalNumber == 0) fTotalNumber = 1;
336 G4double energyMean = fSumEnergy / fTotalNumber;
337 G4double energyRms = std::sqrt(fSumEnergy2 / fTotalNumber - energyMean * energyMean);
338 G4cout << "Mean Energy :" << G4BestUnit(energyMean, "Energy") << " +- "
339 << G4BestUnit(energyRms, "Energy") << "\n";
340 G4double polarizationMean = fSumPolarization / fTotalNumber;
341 G4double polarizationRms =
342 std::sqrt(fSumPolarization2 / fTotalNumber - polarizationMean * polarizationMean);
343 G4cout << "Mean Polarization :" << polarizationMean << " +- " << polarizationRms << "\n";
344 }
345
346 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
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 }
356
357 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
358
359 void RunAction::ParticleStatistics::Merge(const G4VAccumulable& other)
360 {
361 auto rstat = dynamic_cast<const RunAction::ParticleStatistics&>(other);
362
363 fCurrentNumber += rstat.fCurrentNumber;
364 fTotalNumber += rstat.fTotalNumber;
365 fTotalNumber2 += rstat.fTotalNumber2;
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 }
373
374 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
375