| Geant4 Cross Reference |
1 // 1
2 // *******************************************
3 // * License and Disclaimer
4 // *
5 // * The Geant4 software is copyright of th
6 // * the Geant4 Collaboration. It is provided
7 // * conditions of the Geant4 Software License
8 // * LICENSE and available at http://cern.ch/
9 // * include a list of copyright holders.
10 // *
11 // * Neither the authors of this software syst
12 // * institutes,nor the agencies providing fin
13 // * work make any representation or warran
14 // * regarding this software system or assum
15 // * use. Please see the license in the file
16 // * for the full disclaimer and the limitatio
17 // *
18 // * This code implementation is the result
19 // * technical work of the GEANT4 collaboratio
20 // * By using, copying, modifying or distri
21 // * any work based on the software) you ag
22 // * use in resulting scientific publicati
23 // * acceptance of all terms of the Geant4 Sof
24 // *******************************************
25 //
26 /// \file polarisation/Pol01/src/RunAction.cc
27 /// \brief Implementation of the RunAction cla
28 //
29 //
30 //....oooOO0OOooo........oooOO0OOooo........oo
31 //....oooOO0OOooo........oooOO0OOooo........oo
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........oo
54
55 RunAction::RunAction(DetectorConstruction* det
56 : fDetector(det), fPrimary(prim), fAnalysisM
57 {
58 fGamma = G4Gamma::Gamma();
59 fElectron = G4Electron::Electron();
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
77 //....oooOO0OOooo........oooOO0OOooo........oo
78
79 RunAction::~RunAction() {}
80
81 //....oooOO0OOooo........oooOO0OOooo........oo
82
83 void RunAction::BeginOfRunAction(const G4Run*
84 {
85 G4cout << "### Run " << aRun->GetRunID() <<
86
87 auto accumulableManager = G4AccumulableManag
88 accumulableManager->Reset();
89
90 // save Rndm status
91 // G4RunManager::GetRunManager()->SetRandom
92 // CLHEP::HepRandom::showEngineStatus();
93
94 fTotalEventCount = 0;
95
96 // Open file histogram file
97 fAnalysisManager->OpenFile();
98 }
99
100 //....oooOO0OOooo........oooOO0OOooo........oo
101
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 }
135
136 //....oooOO0OOooo........oooOO0OOooo........oo
137
138 void RunAction::BookHisto()
139 {
140 // Always creating analysis manager
141 fAnalysisManager = G4AnalysisManager::Instan
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 t
152 const G4String id[] = {"h1", "h2", "h3", "h4
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", //
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 }
215
216 //....oooOO0OOooo........oooOO0OOooo........oo
217
218 void RunAction::EndOfRunAction(const G4Run* aR
219 {
220 // Total number of events in run (all threa
221 G4int NbOfEvents = aRun->GetNumberOfEventToB
222 // G4int NbOfEvents = aRun->GetNumberOfEven
223
224 if (NbOfEvents == 0) return;
225
226 G4int prec = G4cout.precision(5);
227
228 G4Material* material = fDetector->GetMateria
229 G4double density = material->GetDensity();
230
231 if (fPrimary != nullptr) {
232 G4ParticleDefinition* particle = fPrimary-
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
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........oo
276
277 void RunAction::EventFinished()
278 {
279 auto accManager = G4AccumulableManager::Inst
280
281 ++fTotalEventCount;
282
283 dynamic_cast<ParticleStatistics*>(accManager
284 dynamic_cast<ParticleStatistics*>(accManager
285 dynamic_cast<ParticleStatistics*>(accManager
286 }
287
288 //....oooOO0OOooo........oooOO0OOooo........oo
289
290 RunAction::ParticleStatistics::ParticleStatist
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........oo
304
305 RunAction::ParticleStatistics::~ParticleStatis
306
307 //....oooOO0OOooo........oooOO0OOooo........oo
308
309 void RunAction::ParticleStatistics::EventFinis
310 {
311 fTotalNumber += fCurrentNumber;
312 fTotalNumber2 += fCurrentNumber * fCurrentNu
313 fCurrentNumber = 0;
314 }
315 //....oooOO0OOooo........oooOO0OOooo........oo
316
317 void RunAction::ParticleStatistics::FillData(G
318 G
319 {
320 ++fCurrentNumber;
321 fSumEnergy += kinEnergy;
322 fSumEnergy2 += kinEnergy * kinEnergy;
323 fSumPolarization += longitudinalPolarization
324 fSumPolarization2 += longitudinalPolarizatio
325 fSumCosTheta += costheta;
326 fSumCosTheta2 += costheta * costheta;
327 }
328
329 //....oooOO0OOooo........oooOO0OOooo........oo
330
331 void RunAction::ParticleStatistics::PrintResul
332 {
333 G4cout << "Mean Number per Event :" << G4dou
334 << "\n";
335 if (fTotalNumber == 0) fTotalNumber = 1;
336 G4double energyMean = fSumEnergy / fTotalNum
337 G4double energyRms = std::sqrt(fSumEnergy2 /
338 G4cout << "Mean Energy :" << G4BestUnit(ener
339 << G4BestUnit(energyRms, "Energy") <<
340 G4double polarizationMean = fSumPolarization
341 G4double polarizationRms =
342 std::sqrt(fSumPolarization2 / fTotalNumber
343 G4cout << "Mean Polarization :" << polarizat
344 }
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 }
356
357 //....oooOO0OOooo........oooOO0OOooo........oo
358
359 void RunAction::ParticleStatistics::Merge(cons
360 {
361 auto rstat = dynamic_cast<const RunAction::P
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........oo
375