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