| 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 electromagnetic/TestEm18/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 "HistoManager.hh"
37 #include "PrimaryGeneratorAction.hh"
38
39 #include "G4EmCalculator.hh"
40 #include "G4Run.hh"
41 #include "G4UnitsTable.hh"
42 #include "Randomize.hh"
43
44 #include <iomanip>
45
46 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
47
48 RunAction::RunAction(DetectorConstruction* det, PrimaryGeneratorAction* kin)
49 : fDetector(det), 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 // initialisation
66 //
67 fNbSteps = 0;
68 fTrackLength = 0.;
69 fStepMin = DBL_MAX;
70 fStepMax = 0.;
71
72 fEdepPrimary = fEdepSecondary = fEdepTotal = 0.;
73 fEdepPrimMin = fEdepSecMin = fEdepTotMin = DBL_MAX;
74 fEdepPrimMax = fEdepSecMax = fEdepTotMax = 0.;
75
76 fEnergyTransfered = 0.;
77 fEtransfMin = DBL_MAX;
78 fEtransfMax = 0.;
79
80 fEnergyLost = 0.;
81 fElostMin = DBL_MAX;
82 fElostMax = 0.;
83
84 fEnergyBalance = 0.;
85 fEbalMin = DBL_MAX;
86 fEbalMax = 0.;
87
88 // histograms
89 //
90 G4AnalysisManager* analysisManager = G4AnalysisManager::Instance();
91 if (analysisManager->IsActive()) {
92 analysisManager->OpenFile();
93 }
94
95 // show Rndm status
96 CLHEP::HepRandom::showEngineStatus();
97 }
98
99 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
100
101 void RunAction::CountProcesses(G4String procName)
102 {
103 std::map<G4String, G4int>::iterator it = fProcCounter.find(procName);
104 if (it == fProcCounter.end()) {
105 fProcCounter[procName] = 1;
106 }
107 else {
108 fProcCounter[procName]++;
109 }
110 }
111
112 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
113
114 void RunAction::TrackLength(G4double step)
115 {
116 fTrackLength += step;
117 fNbSteps++;
118 if (step < fStepMin) fStepMin = step;
119 if (step > fStepMax) fStepMax = step;
120 }
121
122 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
123
124 void RunAction::EnergyDeposited(G4double edepPrim, G4double edepSecond)
125 {
126 fEdepPrimary += edepPrim;
127 if (edepPrim < fEdepPrimMin) fEdepPrimMin = edepPrim;
128 if (edepPrim > fEdepPrimMax) fEdepPrimMax = edepPrim;
129
130 fEdepSecondary += edepSecond;
131 if (edepSecond < fEdepSecMin) fEdepSecMin = edepSecond;
132 if (edepSecond > fEdepSecMax) fEdepSecMax = edepSecond;
133 }
134
135 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
136
137 void RunAction::EnergyTransferedByProcess(G4String process, G4double energy)
138 {
139 std::map<G4String, MinMaxData>::iterator it = fEtransfByProcess.find(process);
140 if (it == fEtransfByProcess.end()) {
141 fEtransfByProcess[process] = MinMaxData(1, energy, energy, energy);
142 }
143 else {
144 MinMaxData& data = it->second;
145 data.fCount++;
146 data.fVsum += energy;
147 // update min max
148 G4double emin = data.fVmin;
149 if (energy < emin) data.fVmin = energy;
150 G4double emax = data.fVmax;
151 if (energy > emax) data.fVmax = energy;
152 }
153 }
154
155 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
156
157 void RunAction::EnergyTransfered(G4double energy)
158 {
159 fEnergyTransfered += energy;
160 if (energy < fEtransfMin) fEtransfMin = energy;
161 if (energy > fEtransfMax) fEtransfMax = energy;
162 }
163
164 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
165
166 void RunAction::TotalEnergyLost(G4double energy)
167 {
168 fEnergyLost += energy;
169 if (energy < fElostMin) fElostMin = energy;
170 if (energy > fElostMax) fElostMax = energy;
171 }
172
173 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
174
175 void RunAction::EnergyBalance(G4double energy)
176 {
177 fEnergyBalance += energy;
178 if (energy < fEbalMin) fEbalMin = energy;
179 if (energy > fEbalMax) fEbalMax = energy;
180 }
181
182 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
183
184 void RunAction::TotalEnergyDeposit(G4double energy)
185 {
186 fEdepTotal += energy;
187 if (energy < fEdepTotMin) fEdepTotMin = energy;
188 if (energy > fEdepTotMax) fEdepTotMax = energy;
189 }
190
191 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
192
193 void RunAction::EnergySpectrumOfSecondaries(G4String particle, G4double energy)
194 {
195 std::map<G4String, MinMaxData>::iterator it = fEkinOfSecondaries.find(particle);
196 if (it == fEkinOfSecondaries.end()) {
197 fEkinOfSecondaries[particle] = MinMaxData(1, energy, energy, energy);
198 }
199 else {
200 MinMaxData& data = it->second;
201 data.fCount++;
202 data.fVsum += energy;
203 // update min max
204 G4double emin = data.fVmin;
205 if (energy < emin) data.fVmin = energy;
206 G4double emax = data.fVmax;
207 if (energy > emax) data.fVmax = energy;
208 }
209 }
210
211 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
212
213 void RunAction::EndOfRunAction(const G4Run* aRun)
214 {
215 G4int nbEvents = aRun->GetNumberOfEvent();
216 if (nbEvents == 0) return;
217
218 G4Material* material = fDetector->GetMaterial();
219 G4double length = fDetector->GetSize();
220 G4double density = material->GetDensity();
221
222 G4ParticleDefinition* particle = fPrimary->GetParticleGun()->GetParticleDefinition();
223 G4String partName = particle->GetParticleName();
224 G4double ePrimary = fPrimary->GetParticleGun()->GetParticleEnergy();
225
226 G4int prec = G4cout.precision(3);
227 G4cout << "\n ======================== run summary ======================\n";
228 G4cout << "\n The run was " << nbEvents << " " << partName << " of "
229 << G4BestUnit(ePrimary, "Energy") << " through " << G4BestUnit(length, "Length") << " of "
230 << material->GetName() << " (density: " << G4BestUnit(density, "Volumic Mass") << ")";
231 G4cout << G4endl;
232
233 if (particle->GetPDGCharge() == 0.) return;
234
235 G4cout.precision(4);
236
237 // frequency of processes
238 //
239 G4cout << "\n Process defining step :" << G4endl;
240 G4int index = 0;
241 for (const auto& procCounter : fProcCounter) {
242 G4String procName = procCounter.first;
243 G4int count = procCounter.second;
244 G4String space = " ";
245 if (++index % 4 == 0) space = "\n";
246 G4cout << " " << std::setw(15) << procName << "=" << std::setw(7) << count << space;
247 }
248 G4cout << G4endl;
249
250 // track length
251 //
252 G4double trackLPerEvent = fTrackLength / nbEvents;
253 G4double nbStepPerEvent = double(fNbSteps) / nbEvents;
254 G4double stepSize = fTrackLength / fNbSteps;
255
256 G4cout << "\n TrackLength = " << G4BestUnit(trackLPerEvent, "Length")
257 << " nb of steps = " << nbStepPerEvent
258 << " stepSize = " << G4BestUnit(stepSize, "Length") << " ("
259 << G4BestUnit(fStepMin, "Length") << "--> " << G4BestUnit(fStepMax, "Length") << ")"
260 << G4endl;
261
262 // continuous energy deposited by primary track dE1
263 //
264 G4double energyPerEvent = fEdepPrimary / nbEvents;
265
266 G4cout << "\n Energy continuously deposited along primary track"
267 << " (restricted dE/dx) dE1 = " << G4BestUnit(energyPerEvent, "Energy") << " ("
268 << G4BestUnit(fEdepPrimMin, "Energy") << " --> " << G4BestUnit(fEdepPrimMax, "Energy")
269 << ")" << G4endl;
270
271 // eveluation of dE1 from reading restricted Range table
272 //
273 G4EmCalculator emCal;
274
275 G4double r0 = emCal.GetRangeFromRestricteDEDX(ePrimary, particle, material);
276 G4double r1 = r0 - trackLPerEvent;
277 G4double etry = ePrimary - energyPerEvent;
278 G4double efinal = 0.;
279 if (r1 > 0.) efinal = GetEnergyFromRestrictedRange(r1, particle, material, etry);
280 G4double dEtable = ePrimary - efinal;
281 G4double ratio = 0.;
282 if (dEtable > 0.) ratio = energyPerEvent / dEtable;
283
284 G4cout << "\n Evaluation of dE1 from reading restricted Range table : dE1_table = "
285 << G4BestUnit(dEtable, "Energy") << " ---> dE1/dE1_table = " << ratio << G4endl;
286
287 // energy transfered to secondary particles by process : dE2
288 //
289 G4cout << "\n Energy transfered to secondary particles :" << G4endl;
290 std::map<G4String, MinMaxData>::iterator it1;
291 for (it1 = fEtransfByProcess.begin(); it1 != fEtransfByProcess.end(); it1++) {
292 G4String name = it1->first;
293 MinMaxData data = it1->second;
294 energyPerEvent = data.fVsum / nbEvents;
295 G4double eMin = data.fVmin;
296 G4double eMax = data.fVmax;
297
298 G4cout << " " << std::setw(17) << "due to " + name << ": dE2 = " << std::setw(6)
299 << G4BestUnit(energyPerEvent, "Energy") << " (" << G4BestUnit(eMin, "Energy") << " --> "
300 << G4BestUnit(eMax, "Energy") << ")" << G4endl;
301 }
302
303 // total energy tranfered : dE3 = sum of dE2
304 //
305 energyPerEvent = fEnergyTransfered / nbEvents;
306
307 G4cout << "\n Total energy transfered to secondaries : dE3 = sum of dE2 = "
308 << G4BestUnit(energyPerEvent, "Energy") << " (" << G4BestUnit(fEtransfMin, "Energy")
309 << " --> " << G4BestUnit(fEtransfMax, "Energy") << ")" << G4endl;
310
311 // total energy lost by incident particle : dE4 = dE1 + dE3
312 //
313 energyPerEvent = fEnergyLost / nbEvents;
314
315 G4cout << "\n Total energy lost by incident particle : dE4 = dE1 + dE3 = "
316 << G4BestUnit(energyPerEvent, "Energy") << " (" << G4BestUnit(fElostMin, "Energy")
317 << " --> " << G4BestUnit(fElostMax, "Energy") << ")" << G4endl;
318
319 // calcul of energy lost from energy balance : dE4_bal = E_in - E_out
320 //
321 energyPerEvent = fEnergyBalance / nbEvents;
322
323 G4cout << "\n calcul of dE4 from energy balance : dE4_bal = E_in - E_out = "
324 << G4BestUnit(energyPerEvent, "Energy") << " (" << G4BestUnit(fEbalMin, "Energy")
325 << " --> " << G4BestUnit(fEbalMax, "Energy") << ")" << G4endl;
326
327 // eveluation of dE4 from reading full Range table
328 //
329 r0 = emCal.GetCSDARange(ePrimary, particle, material);
330 r1 = r0 - trackLPerEvent;
331 etry = ePrimary - energyPerEvent;
332 efinal = 0.;
333 if (r1 > 0.) efinal = GetEnergyFromCSDARange(r1, particle, material, etry);
334 dEtable = ePrimary - efinal;
335 ratio = 0.;
336 if (dEtable > 0.) ratio = energyPerEvent / dEtable;
337
338 G4cout << "\n Evaluation of dE4 from reading full Range table : dE4_table = "
339 << G4BestUnit(dEtable, "Energy") << " ---> dE4/dE4_table = " << ratio << G4endl;
340
341 // energy spectrum of secondary particles
342 //
343 G4cout << "\n Energy spectrum of secondary particles :" << G4endl;
344 std::map<G4String, MinMaxData>::iterator it2;
345 for (it2 = fEkinOfSecondaries.begin(); it2 != fEkinOfSecondaries.end(); it2++) {
346 G4String name = it2->first;
347 MinMaxData data = it2->second;
348 G4int count = data.fCount;
349 G4double eMean = data.fVsum / count;
350 G4double eMin = data.fVmin;
351 G4double eMax = data.fVmax;
352
353 G4cout << " " << std::setw(13) << name << ": " << std::setw(7) << count
354 << " Emean = " << std::setw(6) << G4BestUnit(eMean, "Energy") << " ("
355 << G4BestUnit(eMin, "Energy") << " --> " << G4BestUnit(eMax, "Energy") << ")" << G4endl;
356 }
357 G4cout << G4endl;
358
359 // continuous energy deposited by secondary tracks dE5
360 // (only if secondary particles are tracked)
361 //
362 if (fEdepSecondary > 0.) {
363 energyPerEvent = fEdepSecondary / nbEvents;
364
365 G4cout << "\n Energy continuously deposited along secondary tracks"
366 << " (restricted dE/dx) dE5 = " << G4BestUnit(energyPerEvent, "Energy") << " ("
367 << G4BestUnit(fEdepSecMin, "Energy") << " --> " << G4BestUnit(fEdepSecMax, "Energy")
368 << ")" << G4endl;
369
370 // total energy deposited : dE6 = dE1 + dE5
371 //
372 energyPerEvent = fEdepTotal / nbEvents;
373
374 G4cout << "\n Total energy deposited : dE6 = dE1 + dE5 = "
375 << G4BestUnit(energyPerEvent, "Energy") << " (" << G4BestUnit(fEdepTotMin, "Energy")
376 << " --> " << G4BestUnit(fEdepTotMax, "Energy") << ") \n"
377 << G4endl;
378 }
379
380 G4cout.precision(prec);
381
382 // clear maps
383 //
384 fProcCounter.clear();
385 fEtransfByProcess.clear();
386 fEkinOfSecondaries.clear();
387
388 // save histograms
389 G4AnalysisManager* analysisManager = G4AnalysisManager::Instance();
390 if (analysisManager->IsActive()) {
391 analysisManager->Write();
392 analysisManager->CloseFile();
393 }
394
395 // show Rndm status
396 CLHEP::HepRandom::showEngineStatus();
397 }
398
399 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
400
401 G4double RunAction::GetEnergyFromRestrictedRange(G4double range, G4ParticleDefinition* particle,
402 G4Material* material, G4double Etry)
403 {
404 G4EmCalculator emCal;
405
406 G4double Energy = Etry, dE = 0., dEdx;
407 G4double r, dr;
408 G4double err = 1., errmax = 0.00001;
409 G4int iter = 0, itermax = 10;
410 while (err > errmax && iter < itermax) {
411 iter++;
412 Energy -= dE;
413 r = emCal.GetRangeFromRestricteDEDX(Energy, particle, material);
414 dr = r - range;
415 dEdx = emCal.GetDEDX(Energy, particle, material);
416 dE = dEdx * dr;
417 err = std::abs(dE) / Energy;
418 }
419 if (iter == itermax) {
420 G4cout << "\n ---> warning: RunAction::GetEnergyFromRestRange() did not converge"
421 << " Etry = " << G4BestUnit(Etry, "Energy")
422 << " Energy = " << G4BestUnit(Energy, "Energy") << " err = " << err
423 << " iter = " << iter << G4endl;
424 }
425
426 return Energy;
427 }
428
429 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
430
431 G4double RunAction::GetEnergyFromCSDARange(G4double range, G4ParticleDefinition* particle,
432 G4Material* material, G4double Etry)
433 {
434 G4EmCalculator emCal;
435
436 G4double Energy = Etry, dE = 0., dEdx;
437 G4double r, dr;
438 G4double err = 1., errmax = 0.00001;
439 G4int iter = 0, itermax = 10;
440 while (err > errmax && iter < itermax) {
441 iter++;
442 Energy -= dE;
443 r = emCal.GetCSDARange(Energy, particle, material);
444 dr = r - range;
445 dEdx = emCal.ComputeTotalDEDX(Energy, particle, material);
446 dE = dEdx * dr;
447 err = std::abs(dE) / Energy;
448 }
449 if (iter == itermax) {
450 G4cout << "\n ---> warning: RunAction::GetEnergyFromCSDARange() did not converge"
451 << " Etry = " << G4BestUnit(Etry, "Energy")
452 << " Energy = " << G4BestUnit(Energy, "Energy") << " err = " << err
453 << " iter = " << iter << G4endl;
454 }
455
456 return Energy;
457 }
458
459 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
460