source: trunk/examples/extended/medical/fanoCavity/src/RunAction.cc

Last change on this file was 1337, checked in by garnier, 14 years ago

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