Geant4  10.03.p02
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
Run.cc
Go to the documentation of this file.
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: Run.cc 71035 2013-06-10 09:17:35Z gcosmo $
30 //
31 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
32 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
33 
34 #include "Run.hh"
35 #include "DetectorConstruction.hh"
36 #include "PrimaryGeneratorAction.hh"
37 #include "HistoManager.hh"
38 
39 #include "G4Run.hh"
40 #include "G4RunManager.hh"
41 #include "G4UnitsTable.hh"
42 #include "G4EmCalculator.hh"
43 #include "G4Electron.hh"
44 
45 #include "G4SystemOfUnits.hh"
46 #include "Randomize.hh"
47 #include <iomanip>
48 
49 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
50 
52 :fDetector(det),fKinematic(kin), fProcCounter(0), fMateWall(0),fMateCavity(0)
53 
54 {
55 
56  //geometry
57  //
58  fWallThickness = fDetector->GetWallThickness();
59  fWallRadius = fDetector->GetWallRadius();
60  fMateWall = fDetector->GetWallMaterial();
61  fDensityWall = fMateWall->GetDensity();
62 
63  fCavityThickness = fDetector->GetCavityThickness();
64  fCavityRadius = fDetector->GetCavityRadius();
65  fSurfaceCavity = CLHEP::pi*fCavityRadius*fCavityRadius;
66  fVolumeCavity = fSurfaceCavity*fCavityThickness;
67  fMateCavity = fDetector->GetCavityMaterial();
68  fDensityCavity = fMateCavity->GetDensity();
69  fMassCavity = fVolumeCavity*fDensityCavity;
70 
71  //process counter
72  //
73  fProcCounter = new ProcessesCount;
74 
75  //kinetic energy of charged secondary a creation
76  //
77  fEsecondary = fEsecondary2 = 0.;
78  fNbSec = 0;
79 
80  //charged particles and energy flow in cavity
81  //
82  fPartFlowCavity[0] = fPartFlowCavity[1] = 0;
83  fEnerFlowCavity[0] = fEnerFlowCavity[1] = 0.;
84 
85  //total energy deposit and charged track segment in cavity
86  //
87  fEdepCavity = fEdepCavity2 = fTrkSegmCavity = 0.;
88  fNbEventCavity = 0;
89 
90  //survey convergence
91  //
92  fOldEmean = fOldDose = 0.;
93 
94  //stepLenth of charged particles
95  //
96  fStepWall = fStepWall2 = fStepCavity = fStepCavity2 =0.;
97  fNbStepWall = fNbStepCavity = 0;
98 
99  //histograms
100  //
101  G4AnalysisManager* analysisManager = G4AnalysisManager::Instance();
102  if ( analysisManager->IsActive() ) {
103  analysisManager->OpenFile();
104  }
105 }
106 
107 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
108 
109 Run::~Run()
110 {
111 }
112 
113 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
114 
115 void Run::EndOfRun()
116 { // Only call by Master thread
117  std::ios::fmtflags mode = G4cout.flags();
118  G4cout.setf(std::ios::fixed,std::ios::floatfield);
119 
120  if (numberOfEvent == 0) return;
121 
122  //run conditions
123  //
124  G4ParticleDefinition* particle = fKinematic->GetParticleGun()
126  G4String partName = particle->GetParticleName();
127  G4double energy = fKinematic->GetParticleGun()->GetParticleEnergy();
128 
129  G4cout <<"\n ======================== run summary ======================\n";
130 
131  G4int prec = G4cout.precision(3);
132 
133  G4cout <<"\n The run consists of "<<numberOfEvent<<" "<< partName << " of "
134  << G4BestUnit(energy,"Energy") << " through 2*"
135  << G4BestUnit(fWallThickness,"Length") << " of "
136  << fMateWall->GetName() << " (density: "
137  << G4BestUnit(fDensityWall,"Volumic Mass") << ")" << G4endl;
138 
139  G4cout << "\n the cavity is "
140  << G4BestUnit(fCavityThickness,"Length") << " of "
141  << fMateCavity->GetName() << " (density: "
142  << G4BestUnit(fDensityCavity,"Volumic Mass") << "); Mass = "
143  << G4BestUnit(fMassCavity,"Mass") << G4endl;
144 
145  G4cout<<"\n ============================================================\n";
146 
147  //frequency of processes
148  //
149  G4cout << "\n Process calls frequency --->";
150  for (size_t i=0; i< fProcCounter->size();i++) {
151  G4String procName = (*fProcCounter)[i]->GetName();
152  G4int count = (*fProcCounter)[i]->GetCounter();
153  G4cout << " " << procName << "= " << count;
154  }
155  G4cout << G4endl;
156 
157  //extract cross sections with G4EmCalculator
158  //
159  G4EmCalculator emCalculator;
160  G4cout << "\n Gamma crossSections in wall material :";
161  G4double sumc = 0.0;
162  for (size_t i=0; i< fProcCounter->size();i++) {
163  G4String procName = (*fProcCounter)[i]->GetName();
164  G4double massSigma =
165  emCalculator.ComputeCrossSectionPerVolume(energy,particle,
166  procName,fMateWall)/fDensityWall;
167  if (massSigma > 0.) {
168  sumc += massSigma;
169  G4cout << " " << procName << "= "
170  << G4BestUnit(massSigma, "Surface/Mass");
171  }
172  }
173  G4cout << " --> total= "
174  << G4BestUnit(sumc, "Surface/Mass") << G4endl;
175 
176  //mean kinetic energy of secondary electrons
177  //
178  if (fNbSec == 0) return;
179  G4double meanEsecond = fEsecondary/fNbSec,meanEsecond2= fEsecondary2/fNbSec;
180  G4double varianceEsec = meanEsecond2 - meanEsecond*meanEsecond;
181  G4double dToverT = 0.;
182  if (varianceEsec>0.) dToverT = std::sqrt(varianceEsec/fNbSec)/meanEsecond;
183  G4double csdaRange =
184  emCalculator.GetCSDARange(meanEsecond,G4Electron::Electron(),fMateWall);
185 
186  G4cout.precision(4);
187  G4cout
188  << "\n Mean energy of secondary e- = " << G4BestUnit(meanEsecond,"Energy")
189  << " +- " << 100*dToverT << " %"
190  << " (--> range in wall material = " << G4BestUnit(csdaRange,"Length")
191  << ")"
192  << G4endl;
193 
194  //compute mass energy transfer coefficient
195  //
196  G4double massTransfCoef = sumc*meanEsecond/energy;
197 
198  G4cout << " Mass_energy_transfer coef: "
199  << G4BestUnit(massTransfCoef, "Surface/Mass")
200  << G4endl;
201 
202  //stopping power from EmCalculator
203  //
204  G4double dedxWall =
205  emCalculator.GetDEDX(meanEsecond,G4Electron::Electron(),fMateWall);
206  dedxWall /= fDensityWall;
207  G4double dedxCavity =
208  emCalculator.GetDEDX(meanEsecond,G4Electron::Electron(),fMateCavity);
209  dedxCavity /= fDensityCavity;
210 
211  G4cout
212  << "\n StoppingPower in wall = "
213  << G4BestUnit(dedxWall,"Energy*Surface/Mass")
214  << "\n in cavity = "
215  << G4BestUnit(dedxCavity,"Energy*Surface/Mass")
216  << G4endl;
217 
218  //charged particle flow in cavity
219  //
220  G4cout
221  << "\n Charged particle flow in cavity :"
222  << "\n Enter --> nbParticles = " << fPartFlowCavity[0]
223  << "\t Energy = " << G4BestUnit (fEnerFlowCavity[0], "Energy")
224  << "\n Exit --> nbParticles = " << fPartFlowCavity[1]
225  << "\t Energy = " << G4BestUnit (fEnerFlowCavity[1], "Energy")
226  << G4endl;
227 
228  if (fPartFlowCavity[0] == 0) return;
229 
230  //beam energy fluence
231  //
232  G4double rBeam = fWallRadius*(fKinematic->GetBeamRadius());
233  G4double surfaceBeam = CLHEP::pi*rBeam*rBeam;
234 
235  //error on Edep in cavity
236  //
237  if (fNbEventCavity == 0) return;
238  G4double meanEdep = fEdepCavity/fNbEventCavity;
239  G4double meanEdep2 = fEdepCavity2/fNbEventCavity;
240  G4double varianceEdep = meanEdep2 - meanEdep*meanEdep;
241  G4double dEoverE = 0.;
242  if(varianceEdep>0.) dEoverE=std::sqrt(varianceEdep/fNbEventCavity)/meanEdep;
243 
244  //total dose in cavity
245  //
246  G4double doseCavity = fEdepCavity/fMassCavity;
247  G4double doseOverBeam = doseCavity*surfaceBeam/(numberOfEvent*energy);
248 
249  //track length in cavity
250  G4double meantrack = fTrkSegmCavity/fPartFlowCavity[0];
251 
252  G4cout.precision(4);
253  G4cout
254  << "\n Total edep in cavity = " << G4BestUnit(fEdepCavity,"Energy")
255  << " +- " << 100*dEoverE << " %"
256  << "\t Total charged trackLength = " <<G4BestUnit(fTrkSegmCavity,"Length")
257  << " (mean value = " << G4BestUnit(meantrack,"Length") << ")"
258  << "\n Total dose in cavity = " << doseCavity/(MeV/mg) << " MeV/mg"
259  << "\n Dose/EnergyFluence = " << G4BestUnit(doseOverBeam,"Surface/Mass")
260  << G4endl;
261 
262  //ratio simulation/theory
263  //
264  G4double ratio = doseOverBeam/massTransfCoef;
265  G4double error = ratio*std::sqrt(dEoverE*dEoverE + dToverT*dToverT);
266 
267  G4cout.precision(5);
268  G4cout
269  << "\n (Dose/EnergyFluence)/Mass_energy_transfer = " << ratio
270  << " +- " << error << G4endl;
271 
272  //compute mean step size of charged particles
273  //
274  fStepWall /= fNbStepWall; fStepWall2 /= fNbStepWall;
275  G4double rms = fStepWall2 - fStepWall*fStepWall;
276  if (rms>0.) rms = std::sqrt(rms); else rms = 0.;
277 
278  G4cout.precision(4);
279  G4cout
280  << "\n StepSize of ch. tracks in wall = "
281  << G4BestUnit(fStepWall,"Length") << " +- " << G4BestUnit( rms,"Length")
282  << "\t (nbSteps/track = " << double(fNbStepWall)/fNbSec << ")";
283 
284  fStepCavity /= fNbStepCavity; fStepCavity2 /= fNbStepCavity;
285  rms = fStepCavity2 - fStepCavity*fStepCavity;
286  if (rms>0.) rms = std::sqrt(rms); else rms = 0.;
287 
288  G4cout
289  << "\n StepSize of ch. tracks in cavity = "
290  << G4BestUnit(fStepCavity,"Length") << " +- " << G4BestUnit( rms,"Length")
291  << "\t (nbSteps/track = "<<double(fNbStepCavity)/fPartFlowCavity[0] << ")";
292 
293  G4cout << G4endl;
294 
295  // reset default formats
296  G4cout.setf(mode,std::ios::floatfield);
297  G4cout.precision(prec);
298 
299  // delete and remove all contents in fProcCounter
300  while (fProcCounter->size()>0){
301  OneProcessCount* aProcCount=fProcCounter->back();
302  fProcCounter->pop_back();
303  delete aProcCount;
304  }
305  delete fProcCounter;
306 
307 }
308 
309 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
311 {
312  if (NbofEvents == 0) return;
313 
314  //mean kinetic energy of secondary electrons
315  //
316  G4double meanEsecond = 0.;
317  if (fNbSec > 0) meanEsecond = fEsecondary/fNbSec;
318  G4double rateEmean = 0.;
319  // compute variation rate (%), iteration to iteration
320  if (fOldEmean > 0.) rateEmean = 100*(meanEsecond/fOldEmean - 1.);
321  fOldEmean = meanEsecond;
322 
323  //beam energy fluence
324  //
325  G4double rBeam = fWallRadius*(fKinematic->GetBeamRadius());
326  G4double surfaceBeam = CLHEP::pi*rBeam*rBeam;
327  G4double beamEnergy = fKinematic->GetParticleGun()->GetParticleEnergy();
328 
329  //total dose in cavity
330  //
331  G4double doseCavity = fEdepCavity/fMassCavity;
332  G4double doseOverBeam = doseCavity*surfaceBeam/(NbofEvents*beamEnergy);
333  G4double rateDose = 0.;
334  // compute variation rate (%), iteration to iteration
335  if (fOldDose > 0.) rateDose = 100*(doseOverBeam/fOldDose - 1.);
336  fOldDose = doseOverBeam;
337 
338  std::ios::fmtflags mode = G4cout.flags();
339  G4cout.setf(std::ios::fixed,std::ios::floatfield);
340  G4int prec = G4cout.precision(3);
341 
342  G4cout << "\n ---> NbofEvents= " << NbofEvents
343  << " NbOfelectr= " << fNbSec
344  << " Tkin= " << G4BestUnit(meanEsecond,"Energy")
345  << " (" << rateEmean << " %)"
346  << " NbOfelec in cav= " << fPartFlowCavity[0]
347  << " Dose/EnFluence= " << G4BestUnit(doseOverBeam,"Surface/Mass")
348  << " (" << rateDose << " %)"
349  << G4endl;
350 
351  // reset default formats
352  G4cout.setf(mode,std::ios::floatfield);
353  G4cout.precision(prec);
354 }
355 
356 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
357 
358 void Run::Merge(const G4Run* run) {
359 
360  const Run* localRun = static_cast<const Run*>(run);
361 
362  fPartFlowCavity[0]+= localRun->fPartFlowCavity[0];
363  fPartFlowCavity[1]+= localRun->fPartFlowCavity[1];
364  fEnerFlowCavity[0]+= localRun->fEnerFlowCavity[0];
365  fEnerFlowCavity[1]+= localRun->fEnerFlowCavity[1];
366  fEdepCavity += localRun->fEdepCavity;
367  fEdepCavity2 += localRun->fEdepCavity2;
368  fTrkSegmCavity += localRun->fTrkSegmCavity;
369  fNbEventCavity += localRun->fNbEventCavity;
370 
371  fStepWall += localRun->fStepWall;
372  fStepWall2 += localRun->fStepWall2;
373  fStepCavity += localRun->fStepCavity;
374  fStepCavity2 += localRun->fStepCavity2;
375  fNbStepWall += localRun->fNbStepWall;
376  fNbStepCavity += localRun->fNbStepCavity;
377 
378  fEsecondary += localRun->fEsecondary;
379  fEsecondary2 += localRun->fEsecondary2;
380 
381  fNbSec += localRun->fNbSec;
382 
383  // ??? G4double fOldEmean
384  // ??? G4Double fOldDose;
385 
386  // Merge ProcessCount varaibles
387  std::vector<OneProcessCount*>::iterator it;
388  for ( it = localRun->fProcCounter->begin();it !=localRun->fProcCounter->end();
389  it++ )
390  {
391  OneProcessCount* process = *it;
392  for ( G4int i = 0; i < process->GetCounter() ; i++)
393  this->CountProcesses(process->GetName());
394  }
395 
396  G4Run::Merge(run);
397 
398 }
399 
400 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
401 
402 void Run::CountProcesses(G4String procName)
403 {
404  //does the process already encounted ?
405  size_t nbProc = fProcCounter->size();
406  size_t i = 0;
407  while ((i<nbProc)&&((*fProcCounter)[i]->GetName()!=procName)) i++;
408  if (i == nbProc) fProcCounter->push_back( new OneProcessCount(procName));
409 
410  (*fProcCounter)[i]->Count();
411 }
412 
413 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
std::vector< OneProcessCount * > ProcessesCount
The primary generator action class with particle gun.
virtual void Merge(const G4Run *)
Definition: G4Run.cc:54
G4int numberOfEvent
Definition: G4Run.hh:59
void CountProcesses(G4String procName)
Definition: Run.cc:72
G4Material * GetCavityMaterial()
static constexpr double mg
Definition: G4SIunits.hh:184
Run()
Definition: Run.cc:43
G4Material * GetWallMaterial()
const G4String & GetName() const
Definition: G4Material.hh:178
G4double GetDensity() const
Definition: G4Material.hh:180
G4double GetCSDARange(G4double kinEnergy, const G4ParticleDefinition *, const G4Material *, const G4Region *r=nullptr)
#define G4BestUnit(a, b)
#define G4_USE_G4BESTUNIT_FOR_VERBOSE 1
int G4int
Definition: G4Types.hh:78
const G4String & GetParticleName() const
void SurveyConvergence(G4int)
Definition: Run.cc:310
static const double prec
Definition: RanecuEngine.cc:58
G4GLOB_DLL std::ostream G4cout
G4double GetDEDX(G4double kinEnergy, const G4ParticleDefinition *, const G4Material *, const G4Region *r=nullptr)
Definition: G4Run.hh:46
G4String GetName()
G4ParticleGun * GetParticleGun()
G4double energy(const ThreeVector &p, const G4double m)
G4double ComputeCrossSectionPerVolume(G4double kinEnergy, const G4ParticleDefinition *, const G4String &processName, const G4Material *, G4double cut=0.0)
G4ParticleDefinition * GetParticleDefinition() const
static G4Electron * Electron()
Definition: G4Electron.cc:94
Detector construction class to define materials and geometry.
void EndOfRun()
Definition: Run.cc:147
#define G4endl
Definition: G4ios.hh:61
static constexpr double MeV
Definition: G4SIunits.hh:214
static PROLOG_HANDLER error
Definition: xmlrole.cc:112
double G4double
Definition: G4Types.hh:76
Definition: Run.hh:46
virtual void Merge(const G4Run *)
Definition: Run.cc:115
G4CsvAnalysisManager G4AnalysisManager
Definition: g4csv_defs.hh:77
static constexpr double pi
Definition: SystemOfUnits.h:54
G4double GetParticleEnergy() const
~Run()
Definition: Run.cc:72