source: trunk/examples/extended/medical/fanoCavity2/src/RunAction.cc @ 1337

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

tag geant4.9.4 beta 1 + modifs locales
File size: 12.5 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.3 2007/11/05 13:19:16 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  // do not save Rndm status
63  G4RunManager::GetRunManager()->SetRandomNumberStore(false);
64  CLHEP::HepRandom::showEngineStatus();
65 
66  G4cout << "### Run " << aRun->GetRunID() << " start." << G4endl;
67   
68  G4int NbofEvents = aRun->GetNumberOfEventToBeProcessed();
69  if (NbofEvents == 0) return;
70 
71  //run conditions
72  //     
73  G4ParticleDefinition* particleGun
74                    = kinematic->GetParticleGun()->GetParticleDefinition();
75  G4String partName = particleGun->GetParticleName();                           
76  energyGun = kinematic->GetParticleGun()->GetParticleEnergy();
77 
78  //geometry : effective wall volume
79  // 
80  G4double cavityThickness = detector->GetCavityThickness();         
81  G4Material* mateCavity   = detector->GetCavityMaterial();
82  G4double densityCavity   = mateCavity->GetDensity();
83  massCavity = cavityThickness*densityCavity;
84     
85  G4double wallThickness = detector->GetWallThickness();
86  G4Material* mateWall   = detector->GetWallMaterial();
87  G4double densityWall   = mateWall->GetDensity();
88 
89  G4EmCalculator emCal;
90  G4double RangeWall = emCal.GetCSDARange(energyGun,particleGun,mateWall);
91  G4double factor = 1.2;
92  G4double effWallThick = factor*RangeWall;
93  if ((effWallThick > wallThickness)||(effWallThick <= 0.))
94    effWallThick = wallThickness;
95  massWall = 2*effWallThick*densityWall; 
96 
97  G4double massTotal     = massWall + massCavity;
98  G4double massWallRatio = massWall/massTotal; 
99  kinematic->RunInitialisation(effWallThick, massWallRatio ); 
100     
101  G4double massRatio = massCavity/massWall;
102 
103  //check radius
104  //
105  G4double worldRadius = detector->GetWorldRadius();
106  G4double RangeCavity = emCal.GetCSDARange(energyGun,particleGun,mateCavity);
107   
108  std::ios::fmtflags mode = G4cout.flags();
109  G4cout.setf(std::ios::fixed,std::ios::floatfield);
110  G4int prec = G4cout.precision(3);
111 
112  G4cout << "\n ======================== run conditions =====================\n";
113 
114  G4cout << "\n The run will be " << NbofEvents << " "<< partName << " of "
115         << G4BestUnit(energyGun,"Energy") << " through 2*" 
116         << G4BestUnit(effWallThick,"Length") << " of "
117         << mateWall->GetName() << " (density: " 
118         << G4BestUnit(densityWall,"Volumic Mass") << "); Mass/cm2 = "
119         << G4BestUnit(massWall*cm2,"Mass") 
120         << "\n csdaRange: " << G4BestUnit(RangeWall,"Length") << G4endl;
121         
122  G4cout << "\n the cavity is "
123         << G4BestUnit(cavityThickness,"Length") << " of "
124         << mateCavity->GetName() << " (density: " 
125         << G4BestUnit(densityCavity,"Volumic Mass") << "); Mass/cm2 = " 
126         << G4BestUnit(massCavity*cm2,"Mass") 
127         << " --> massRatio = " << std::setprecision(6) << massRatio << G4endl;
128         
129  G4cout.precision(3);   
130  G4cout << " World radius: " << G4BestUnit(worldRadius,"Length")
131         << "; range in cavity: " << G4BestUnit(RangeCavity,"Length")
132         << G4endl;     
133                 
134  G4cout << "\n ============================================================\n";
135         
136  //stopping power from EmCalculator
137  //
138  G4double dedxWall = 
139      emCal.GetDEDX(energyGun,G4Electron::Electron(),mateWall);
140  dedxWall /= densityWall;
141  G4double dedxCavity = 
142      emCal.GetDEDX(energyGun,G4Electron::Electron(),mateCavity);
143  dedxCavity /= densityCavity;
144 
145  G4cout << std::setprecision(4)
146         << "\n StoppingPower in wall   = " 
147         << G4BestUnit(dedxWall,"Energy*Surface/Mass")
148         << "\n               in cavity = " 
149         << G4BestUnit(dedxCavity,"Energy*Surface/Mass")
150         << G4endl;
151                 
152  //process counter
153  //
154  ProcCounter = new ProcessesCount;
155 
156  //charged particles and energy flow in cavity
157  //
158  PartFlowCavity[0] = PartFlowCavity[1] = 0;
159  EnerFlowCavity[0] = EnerFlowCavity[1] = 0.;
160       
161  //total energy deposit and charged track segment in cavity
162  //
163  EdepCavity = EdepCavity2 = trkSegmCavity = 0.;
164  nbEventCavity = 0;
165     
166  //stepLenth of charged particles
167  //
168  stepWall = stepWall2 = stepCavity = stepCavity2 =0.;
169  nbStepWall = nbStepCavity = 0;
170   
171  //histograms
172  //
173  histoManager->book();
174 
175  // reset default formats
176  G4cout.setf(mode,std::ios::floatfield);
177  G4cout.precision(prec);     
178}
179
180//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
181
182void RunAction::CountProcesses(G4String procName)
183{
184   //does the process  already encounted ?
185   size_t nbProc = ProcCounter->size();
186   size_t i = 0;
187   while ((i<nbProc)&&((*ProcCounter)[i]->GetName()!=procName)) i++;
188   if (i == nbProc) ProcCounter->push_back( new OneProcessCount(procName));
189
190   (*ProcCounter)[i]->Count();
191}
192
193//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
194
195void RunAction::SurveyConvergence(G4int NbofEvents)
196{ 
197  if (NbofEvents == 0) return;
198   
199               
200  //beam fluence
201  //
202  G4int Nwall   = kinematic->GetWallCount();
203  G4int Ncavity = kinematic->GetCavityCount();
204  G4double Iwall   = Nwall/massWall;   
205  G4double Icavity = Ncavity/massCavity;
206  G4double Iratio  = Icavity/Iwall;
207  G4double Itot    = NbofEvents/(massWall+massCavity);
208  G4double energyFluence = energyGun*Itot;
209           
210  //total dose in cavity
211  //               
212  G4double doseCavity = EdepCavity/massCavity;
213  G4double ratio = doseCavity/energyFluence;
214  G4double err = 100*(ratio-1.);
215
216  std::ios::fmtflags mode = G4cout.flags();
217  G4cout.setf(std::ios::fixed,std::ios::floatfield);
218  G4int prec = G4cout.precision(5);
219 
220  G4cout << "\n--->evntNb= " << NbofEvents
221         << " Nwall= " << Nwall
222         << " Ncav= "  << Ncavity
223         << " Ic/Iw= " << Iratio       
224         << " Ne-_cav= " << PartFlowCavity[0]
225         << " doseCavity/Ebeam= " << ratio
226         << "  (100*(ratio-1) = " << err << " %)"
227         << G4endl;
228         
229  // reset default formats
230  G4cout.setf(mode,std::ios::floatfield);
231  G4cout.precision(prec); 
232}
233
234//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
235
236void RunAction::EndOfRunAction(const G4Run* aRun)
237{
238  std::ios::fmtflags mode = G4cout.flags();
239  G4cout.setf(std::ios::fixed,std::ios::floatfield);
240  G4int prec = G4cout.precision(3);
241 
242  G4int NbofEvents = aRun->GetNumberOfEvent();
243  if (NbofEvents == 0) return;
244
245  //frequency of processes
246  //
247  G4cout << "\n Process calls frequency --->";
248  for (size_t i=0; i< ProcCounter->size();i++) {
249     G4String procName = (*ProcCounter)[i]->GetName();
250     G4int    count    = (*ProcCounter)[i]->GetCounter(); 
251     G4cout << "  " << procName << "= " << count;
252  }
253  G4cout << G4endl;
254       
255  //charged particle flow in cavity
256  //
257  G4cout
258    << "\n Charged particle flow in cavity :"
259    << "\n      Enter --> nbParticles = " << PartFlowCavity[0]
260    << "\t Energy = " << G4BestUnit (EnerFlowCavity[0], "Energy")
261    << "\n      Exit  --> nbParticles = " << PartFlowCavity[1]
262    << "\t Energy = " << G4BestUnit (EnerFlowCavity[1], "Energy")
263    << G4endl;
264             
265  if (PartFlowCavity[0] == 0) return;
266               
267  //beam fluence
268  //
269  G4int Nwall   = kinematic->GetWallCount();
270  G4int Ncavity = kinematic->GetCavityCount(); 
271  G4double Iwall   = Nwall/massWall;
272  G4double Icavity = Ncavity/massCavity;
273  G4double Iratio  = Icavity/Iwall;
274  G4double Itot    = NbofEvents/(massWall+massCavity);
275  G4double energyFluence = energyGun*Itot; 
276 
277  G4cout.precision(5);       
278  G4cout
279    << "\n beamFluence in wall = " << Nwall
280    << "\t in cavity = " << Ncavity
281    << "\t Icav/Iwall = " << Iratio       
282    << "\t energyFluence = " << energyFluence/(MeV*cm2/mg) << " MeV*cm2/mg"
283    << G4endl;
284 
285  //error on Edep in cavity
286  //
287  if (nbEventCavity == 0) return;
288  G4double meanEdep  = EdepCavity/nbEventCavity;
289  G4double meanEdep2 = EdepCavity2/nbEventCavity;
290  G4double varianceEdep = meanEdep2 - meanEdep*meanEdep;
291  G4double dEoverE = 0.;
292  if(varianceEdep>0.) dEoverE = std::sqrt(varianceEdep/nbEventCavity)/meanEdep;
293               
294  //total dose in cavity
295  //               
296  G4double doseCavity = EdepCavity/massCavity;
297  G4double ratio = doseCavity/energyFluence, error = ratio*dEoverE;
298                 
299  G4cout
300    << "\n Total edep in cavity = "      << G4BestUnit(EdepCavity,"Energy")
301    << " +- " << 100*dEoverE << " %"       
302    << "\n Total dose in cavity = " << doseCavity/(MeV*cm2/mg) << " MeV*cm2/mg"
303    << " +- " << 100*dEoverE << " %"         
304    << "\n\n DoseCavity/EnergyFluence = " << ratio
305    << " +- " << error << G4endl;
306   
307
308  //track length in cavity
309  G4double meantrack = trkSegmCavity/PartFlowCavity[0];
310 
311  G4cout.precision(4); 
312  G4cout 
313    << "\n Total charged trackLength in cavity = " 
314    << G4BestUnit(trkSegmCavity,"Length")
315    << "   (mean value = " << G4BestUnit(meantrack,"Length") << ")"       
316    << G4endl;
317                 
318  //compute mean step size of charged particles
319  //
320  stepWall /= nbStepWall; stepWall2 /= nbStepWall;
321  G4double rms = stepWall2 - stepWall*stepWall;       
322  if (rms>0.) rms = std::sqrt(rms); else rms = 0.;
323  G4double nbTrackWall = kinematic->GetWallCount();
324
325  G4cout
326    << "\n StepSize of ch. tracks in wall   = " 
327    << G4BestUnit(stepWall,"Length") << " +- " << G4BestUnit( rms,"Length")
328    << "\t (nbSteps/track = " << double(nbStepWall)/nbTrackWall << ")";
329   
330  stepCavity /= nbStepCavity; stepCavity2 /= nbStepCavity;
331  rms = stepCavity2 - stepCavity*stepCavity;       
332  if (rms>0.) rms = std::sqrt(rms); else rms = 0.;
333
334  G4cout
335    << "\n StepSize of ch. tracks in cavity = " 
336    << G4BestUnit(stepCavity,"Length") << " +- " << G4BestUnit( rms,"Length")
337    << "\t (nbSteps/track = " << double(nbStepCavity)/PartFlowCavity[0] << ")";
338       
339  G4cout << G4endl;
340 
341   // reset default formats
342  G4cout.setf(mode,std::ios::floatfield);
343  G4cout.precision(prec);
344 
345  // delete and remove all contents in ProcCounter
346  while (ProcCounter->size()>0){
347    OneProcessCount* aProcCount=ProcCounter->back();
348    ProcCounter->pop_back();
349    delete aProcCount;
350  }
351  delete ProcCounter;
352 
353  // save histograms
354  histoManager->save();
355 
356  // show Rndm status
357  CLHEP::HepRandom::showEngineStatus();
358}
359
360//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
Note: See TracBrowser for help on using the repository browser.