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