Geant4  v4-10.4-release
 모두 클래스 네임스페이스들 파일들 함수 변수 타입정의 열거형 타입 열거형 멤버 Friends 매크로 그룹들 페이지들
examples/extended/electromagnetic/TestEm18/src/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 // $Id: RunAction.cc 105930 2017-08-31 08:39:49Z gcosmo $
30 //
31 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
32 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
33 
34 #include "RunAction.hh"
35 #include "DetectorConstruction.hh"
36 #include "PrimaryGeneratorAction.hh"
37 #include "HistoManager.hh"
38 
39 #include "G4Run.hh"
40 #include "G4UnitsTable.hh"
41 #include "G4EmCalculator.hh"
42 
43 #include "Randomize.hh"
44 #include <iomanip>
45 
46 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
47 
49 :G4UserRunAction(),fDetector(det), fPrimary(kin), fHistoManager(0)
50 {
51  fHistoManager = new HistoManager();
52 }
53 
54 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
55 
57 {
58  delete fHistoManager;
59 }
60 
61 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
62 
64 {
65  //initialisation
66  //
67  fNbSteps = 0;
68  fTrackLength = 0.;
69  fStepMin = DBL_MAX;
70  fStepMax = 0.;
71 
75 
76  fEnergyTransfered = 0.;
78  fEtransfMax = 0.;
79 
80  fEnergyLost = 0.;
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
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 
115 {
116  fTrackLength += step; fNbSteps++;
117  if (step<fStepMin) fStepMin = step;
118  if (step>fStepMax) fStepMax = step;
119 }
120 
121 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
122 
123 void RunAction::EnergyDeposited (G4double edepPrim, G4double edepSecond)
124 {
125  fEdepPrimary += edepPrim;
126  if (edepPrim<fEdepPrimMin) fEdepPrimMin = edepPrim;
127  if (edepPrim>fEdepPrimMax) fEdepPrimMax = edepPrim;
128 
129  fEdepSecondary += edepSecond;
130  if (edepSecond<fEdepSecMin) fEdepSecMin = edepSecond;
131  if (edepSecond>fEdepSecMax) fEdepSecMax = edepSecond;
132 }
133 
134 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
135 
137 {
138  std::map<G4String, MinMaxData>::iterator it = fEtransfByProcess.find(process);
139  if ( it == fEtransfByProcess.end()) {
140  fEtransfByProcess[process] = MinMaxData(1, energy, energy, energy);
141  }
142  else {
143  MinMaxData& data = it->second;
144  data.fCount++;
145  data.fVsum += energy;
146  //update min max
147  G4double emin = data.fVmin;
148  if (energy < emin) data.fVmin = energy;
149  G4double emax = data.fVmax;
150  if (energy > emax) data.fVmax = energy;
151  }
152 }
153 
154 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
155 
157 {
159  if (energy<fEtransfMin) fEtransfMin = energy;
160  if (energy>fEtransfMax) fEtransfMax = energy;
161 }
162 
163 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
164 
166 {
167  fEnergyLost += energy;
168  if (energy<fElostMin) fElostMin = energy;
169  if (energy>fElostMax) fElostMax = energy;
170 }
171 
172 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
173 
175 {
177  if (energy<fEbalMin) fEbalMin = energy;
178  if (energy>fEbalMax) fEbalMax = energy;
179 }
180 
181 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
182 
184 {
185  fEdepTotal += energy;
186  if (energy<fEdepTotMin) fEdepTotMin = energy;
187  if (energy>fEdepTotMax) fEdepTotMax = energy;
188 }
189 
190 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
191 
193 {
194  std::map<G4String,MinMaxData>::iterator it = fEkinOfSecondaries.find(particle);
195  if ( it == fEkinOfSecondaries.end()) {
196  fEkinOfSecondaries[particle] = MinMaxData(1, energy, energy, energy);
197  }
198  else {
199  MinMaxData& data = it->second;
200  data.fCount++;
201  data.fVsum += energy;
202  //update min max
203  G4double emin = data.fVmin;
204  if (energy < emin) data.fVmin = energy;
205  G4double emax = data.fVmax;
206  if (energy > emax) data.fVmax = energy;
207  }
208 }
209 
210 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
211 
212 void RunAction::EndOfRunAction(const G4Run* aRun)
213 {
214  G4int nbEvents = aRun->GetNumberOfEvent();
215  if (nbEvents == 0) return;
216 
217  G4Material* material = fDetector->GetMaterial();
218  G4double length = fDetector->GetSize();
219  G4double density = material->GetDensity();
220 
223  G4String partName = particle->GetParticleName();
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 "
230  << G4BestUnit(length,"Length") << " of "
231  << material->GetName() << " (density: "
232  << G4BestUnit(density,"Volumic Mass") << ")";
233  G4cout << G4endl;
234 
235  if (particle->GetPDGCharge() == 0.) return;
236 
237  G4cout.precision(4);
238 
239  //frequency of processes
240  //
241  G4cout << "\n Process defining step :" << G4endl;
242  G4int index = 0;
243  for ( const auto& procCounter : fProcCounter ) {
244  G4String procName = procCounter.first;
245  G4int count = procCounter.second;
246  G4String space = " "; if (++index%4 == 0) space = "\n";
247  G4cout << " " << std::setw(15) << procName << "="<< std::setw(7) << count
248  << space;
249  }
250  G4cout << G4endl;
251 
252  //track length
253  //
254  G4double trackLPerEvent = fTrackLength/nbEvents;
255  G4double nbStepPerEvent = double(fNbSteps)/nbEvents;
256  G4double stepSize = fTrackLength/fNbSteps;
257 
258  G4cout
259  << "\n TrackLength = "
260  << G4BestUnit(trackLPerEvent, "Length")
261  << " nb of steps = " << nbStepPerEvent
262  << " stepSize = " << G4BestUnit(stepSize, "Length")
263  << " (" << G4BestUnit(fStepMin, "Length")
264  << "--> " << G4BestUnit(fStepMax, "Length") << ")"
265  << G4endl;
266 
267  //continuous energy deposited by primary track dE1
268  //
269  G4double energyPerEvent = fEdepPrimary/nbEvents;
270 
271  G4cout
272  << "\n Energy continuously deposited along primary track"
273  << " (restricted dE/dx) dE1 = "
274  << G4BestUnit(energyPerEvent, "Energy")
275  << " (" << G4BestUnit(fEdepPrimMin, "Energy")
276  << " --> " << G4BestUnit(fEdepPrimMax, "Energy") << ")"
277  << G4endl;
278 
279  //eveluation of dE1 from reading restricted Range table
280  //
281  G4EmCalculator emCal;
282 
283  G4double r0 = emCal.GetRangeFromRestricteDEDX(ePrimary,particle,material);
284  G4double r1 = r0 - trackLPerEvent;
285  G4double etry = ePrimary - energyPerEvent;
286  G4double efinal = 0.;
287  if (r1 > 0.) efinal = GetEnergyFromRestrictedRange(r1,particle,material,etry);
288  G4double dEtable = ePrimary - efinal;
289  G4double ratio = 0.;
290  if (dEtable > 0.) ratio = energyPerEvent/dEtable;
291 
292  G4cout
293  << "\n Evaluation of dE1 from reading restricted Range table : dE1_table = "
294  << G4BestUnit(dEtable, "Energy")
295  << " ---> dE1/dE1_table = " << ratio
296  << G4endl;
297 
298  // energy transfered to secondary particles by process : dE2
299  //
300  G4cout << "\n Energy transfered to secondary particles :" << G4endl;
301  std::map<G4String,MinMaxData>::iterator it1;
302  for (it1 = fEtransfByProcess.begin(); it1 != fEtransfByProcess.end(); it1++) {
303  G4String name = it1->first;
304  MinMaxData data = it1->second;
305  energyPerEvent = data.fVsum/nbEvents;
306  G4double eMin = data.fVmin;
307  G4double eMax = data.fVmax;
308 
309  G4cout << " " << std::setw(17) << "due to " + name << ": dE2 = "
310  << std::setw(6) << G4BestUnit(energyPerEvent, "Energy")
311  << " (" << G4BestUnit(eMin, "Energy")
312  << " --> " << G4BestUnit(eMax, "Energy")
313  << ")" << G4endl;
314  }
315 
316  // total energy tranfered : dE3 = sum of dE2
317  //
318  energyPerEvent = fEnergyTransfered/nbEvents;
319 
320  G4cout
321  << "\n Total energy transfered to secondaries : dE3 = sum of dE2 = "
322  << G4BestUnit(energyPerEvent, "Energy")
323  << " (" << G4BestUnit(fEtransfMin, "Energy")
324  << " --> " << G4BestUnit(fEtransfMax, "Energy") << ")"
325  << G4endl;
326 
327  // total energy lost by incident particle : dE4 = dE1 + dE3
328  //
329  energyPerEvent = fEnergyLost/nbEvents;
330 
331  G4cout
332  << "\n Total energy lost by incident particle : dE4 = dE1 + dE3 = "
333  << G4BestUnit(energyPerEvent, "Energy")
334  << " (" << G4BestUnit(fElostMin, "Energy")
335  << " --> " << G4BestUnit(fElostMax, "Energy") << ")"
336  << G4endl;
337 
338  // calcul of energy lost from energy balance : dE4_bal = E_in - E_out
339  //
340  energyPerEvent = fEnergyBalance/nbEvents;
341 
342  G4cout
343  << "\n calcul of dE4 from energy balance : dE4_bal = E_in - E_out = "
344  << G4BestUnit(energyPerEvent, "Energy")
345  << " (" << G4BestUnit(fEbalMin, "Energy")
346  << " --> " << G4BestUnit(fEbalMax, "Energy") << ")"
347  << G4endl;
348 
349  //eveluation of dE4 from reading full Range table
350  //
351  r0 = emCal.GetCSDARange(ePrimary,particle,material);
352  r1 = r0 - trackLPerEvent;
353  etry = ePrimary - energyPerEvent;
354  efinal = 0.;
355  if (r1 > 0.) efinal = GetEnergyFromCSDARange(r1,particle,material,etry);
356  dEtable = ePrimary - efinal;
357  ratio = 0.;
358  if (dEtable > 0.) ratio = energyPerEvent/dEtable;
359 
360  G4cout
361  << "\n Evaluation of dE4 from reading full Range table : dE4_table = "
362  << G4BestUnit(dEtable, "Energy")
363  << " ---> dE4/dE4_table = " << ratio
364  << G4endl;
365 
366  //energy spectrum of secondary particles
367  //
368  G4cout << "\n Energy spectrum of secondary particles :" << G4endl;
369  std::map<G4String,MinMaxData>::iterator it2;
370  for (it2 = fEkinOfSecondaries.begin();it2 != fEkinOfSecondaries.end(); it2++){
371  G4String name = it2->first;
372  MinMaxData data = it2->second;
373  G4int count = data.fCount;
374  G4double eMean = data.fVsum/count;
375  G4double eMin = data.fVmin;
376  G4double eMax = data.fVmax;
377 
378  G4cout << " " << std::setw(13) << name << ": " << std::setw(7) << count
379  << " Emean = " << std::setw(6) << G4BestUnit(eMean, "Energy")
380  << " (" << G4BestUnit(eMin, "Energy")
381  << " --> " << G4BestUnit(eMax, "Energy")
382  << ")" << G4endl;
383  }
384  G4cout << G4endl;
385 
386  //continuous energy deposited by secondary tracks dE5
387  // (only if secondary particles are tracked)
388  //
389  if (fEdepSecondary > 0.) {
390  energyPerEvent = fEdepSecondary/nbEvents;
391 
392  G4cout
393  << "\n Energy continuously deposited along secondary tracks"
394  << " (restricted dE/dx) dE5 = "
395  << G4BestUnit(energyPerEvent, "Energy")
396  << " (" << G4BestUnit(fEdepSecMin, "Energy")
397  << " --> " << G4BestUnit(fEdepSecMax, "Energy") << ")"
398  << G4endl;
399 
400  // total energy deposited : dE6 = dE1 + dE5
401  //
402  energyPerEvent = fEdepTotal/nbEvents;
403 
404  G4cout
405  << "\n Total energy deposited : dE6 = dE1 + dE5 = "
406  << G4BestUnit(energyPerEvent, "Energy")
407  << " (" << G4BestUnit(fEdepTotMin, "Energy")
408  << " --> " << G4BestUnit(fEdepTotMax, "Energy") << ") \n"
409  << G4endl;
410 }
411 
412  G4cout.precision(prec);
413 
414  //clear maps
415  //
416  fProcCounter.clear();
417  fEtransfByProcess.clear();
418  fEkinOfSecondaries.clear();
419 
420  //save histograms
421  G4AnalysisManager* analysisManager = G4AnalysisManager::Instance();
422  if ( analysisManager->IsActive() ) {
423  analysisManager->Write();
424  analysisManager->CloseFile();
425  }
426 
427  // show Rndm status
429 }
430 
431 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
432 
434  G4ParticleDefinition* particle, G4Material* material, G4double Etry)
435 {
436  G4EmCalculator emCal;
437 
438  G4double Energy = Etry, dE = 0., dEdx;
439  G4double r, dr;
440  G4double err = 1., errmax = 0.00001;
441  G4int iter = 0 , itermax = 10;
442  while (err > errmax && iter < itermax) {
443  iter++;
444  Energy -= dE;
445  r = emCal.GetRangeFromRestricteDEDX(Energy,particle,material);
446  dr = r - range;
447  dEdx = emCal.GetDEDX(Energy,particle,material);
448  dE = dEdx*dr;
449  err = std::abs(dE)/Energy;
450  }
451  if (iter == itermax) {
452  G4cout
453  << "\n ---> warning: RunAction::GetEnergyFromRestRange() did not converge"
454  << " Etry = " << G4BestUnit(Etry,"Energy")
455  << " Energy = " << G4BestUnit(Energy,"Energy")
456  << " err = " << err
457  << " iter = " << iter << G4endl;
458  }
459 
460  return Energy;
461 }
462 
463 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
464 
466  G4ParticleDefinition* particle, G4Material* material, G4double Etry)
467 {
468  G4EmCalculator emCal;
469 
470  G4double Energy = Etry, dE = 0., dEdx;
471  G4double r, dr;
472  G4double err = 1., errmax = 0.00001;
473  G4int iter = 0 , itermax = 10;
474  while (err > errmax && iter < itermax) {
475  iter++;
476  Energy -= dE;
477  r = emCal.GetCSDARange(Energy,particle,material);
478  dr = r - range;
479  dEdx = emCal.ComputeTotalDEDX(Energy,particle,material);
480  dE = dEdx*dr;
481  err = std::abs(dE)/Energy;
482  }
483  if (iter == itermax) {
484  G4cout
485  << "\n ---> warning: RunAction::GetEnergyFromCSDARange() did not converge"
486  << " Etry = " << G4BestUnit(Etry,"Energy")
487  << " Energy = " << G4BestUnit(Energy,"Energy")
488  << " err = " << err
489  << " iter = " << iter << G4endl;
490  }
491 
492  return Energy;
493 }
494 
495 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
static const G4double dE
G4double GetRangeFromRestricteDEDX(G4double kinEnergy, const G4ParticleDefinition *, const G4Material *, const G4Region *r=nullptr)
const XML_Char * name
Definition: expat.h:151
static const double prec
Definition: RanecuEngine.cc:58
void EnergyTransferedByProcess(G4String procName, G4double energy)
#define G4endl
Definition: G4ios.hh:61
G4double GetDEDX(G4double kinEnergy, const G4ParticleDefinition *, const G4Material *, const G4Region *r=nullptr)
const G4String & GetParticleName() const
G4double GetPDGCharge() const
static const G4double emax
G4double GetEnergyFromCSDARange(G4double, G4ParticleDefinition *, G4Material *, G4double)
G4double GetEnergyFromRestrictedRange(G4double, G4ParticleDefinition *, G4Material *, G4double)
G4int GetNumberOfEvent() const
Definition: G4Run.hh:79
const XML_Char const XML_Char * data
Definition: expat.h:268
const G4String & GetName() const
Definition: G4Material.hh:179
double G4double
Definition: G4Types.hh:76
const G4ParticleDefinition const G4Material *G4double range
double energy
Definition: plottest35.C:25
G4double ComputeTotalDEDX(G4double kinEnergy, const G4ParticleDefinition *, const G4Material *, G4double cut=DBL_MAX)
The primary generator action class with particle gun.
G4CsvAnalysisManager G4AnalysisManager
Definition: g4csv_defs.hh:77
Definition: G4Run.hh:46
#define G4BestUnit(a, b)
#define G4_USE_G4BESTUNIT_FOR_VERBOSE 1
void EnergyDeposited(G4double edepPrim, G4double edepSecond)
int G4int
Definition: G4Types.hh:78
G4double GetParticleEnergy() const
G4double GetCSDARange(G4double kinEnergy, const G4ParticleDefinition *, const G4Material *, const G4Region *r=nullptr)
G4GLOB_DLL std::ostream G4cout
G4int first(char) const
G4ParticleDefinition * GetParticleDefinition() const
#define DBL_MAX
Definition: templates.hh:83
static void showEngineStatus()
Definition: Random.cc:303
Simple detector construction with a box volume placed in a world.
void EnergySpectrumOfSecondaries(G4String particleName, G4double ekin)
G4double GetDensity() const
Definition: G4Material.hh:181