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