ECCE @ EIC Software
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
Run.cc
Go to the documentation of this file. Or view the newest version in sPHENIX GitHub for file Run.cc
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 //
30 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
31 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
32 
33 #include "Run.hh"
34 #include "DetectorConstruction.hh"
35 #include "PrimaryGeneratorAction.hh"
36 #include "HistoManager.hh"
37 
38 #include "G4Run.hh"
39 #include "G4RunManager.hh"
40 #include "G4UnitsTable.hh"
41 #include "G4EmCalculator.hh"
42 #include "G4Electron.hh"
43 
44 #include "G4SystemOfUnits.hh"
45 #include "Randomize.hh"
46 #include <iomanip>
47 
48 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
49 
51 :fDetector(det),fKinematic(kin), fProcCounter(0), fMateWall(0),fMateCavity(0)
52 
53 {
54 
55  //geometry
56  //
61 
69 
70  //process counter
71  //
73 
74  //kinetic energy of charged secondary a creation
75  //
77  fNbSec = 0;
78 
79  //charged particles and energy flow in cavity
80  //
81  fPartFlowCavity[0] = fPartFlowCavity[1] = 0;
82  fEnerFlowCavity[0] = fEnerFlowCavity[1] = 0.;
83 
84  //total energy deposit and charged track segment in cavity
85  //
87  fNbEventCavity = 0;
88 
89  //survey convergence
90  //
91  fOldEmean = fOldDose = 0.;
92 
93  //stepLenth of charged particles
94  //
97 
98  //histograms
99  //
100  G4AnalysisManager* analysisManager = G4AnalysisManager::Instance();
101  if ( analysisManager->IsActive() ) {
102  analysisManager->OpenFile();
103  }
104 }
105 
106 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
107 
108 Run::~Run()
109 {
110 }
111 
112 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
113 
114 void Run::EndOfRun()
115 { // Only call by Master thread
116  std::ios::fmtflags mode = G4cout.flags();
117  G4cout.setf(std::ios::fixed,std::ios::floatfield);
118 
119  if (numberOfEvent == 0) return;
120 
121  //run conditions
122  //
125  G4String partName = particle->GetParticleName();
127 
128  G4cout <<"\n ======================== run summary ======================\n";
129 
130  G4int prec = G4cout.precision(3);
131 
132  G4cout <<"\n The run consists of "<<numberOfEvent<<" "<< partName << " of "
133  << G4BestUnit(energy,"Energy") << " through 2*"
134  << G4BestUnit(fWallThickness,"Length") << " of "
135  << fMateWall->GetName() << " (density: "
136  << G4BestUnit(fDensityWall,"Volumic Mass") << ")" << G4endl;
137 
138  G4cout << "\n the cavity is "
139  << G4BestUnit(fCavityThickness,"Length") << " of "
140  << fMateCavity->GetName() << " (density: "
141  << G4BestUnit(fDensityCavity,"Volumic Mass") << "); Mass = "
142  << G4BestUnit(fMassCavity,"Mass") << G4endl;
143 
144  G4cout<<"\n ============================================================\n";
145 
146  //frequency of processes
147  //
148  G4cout << "\n Process calls frequency --->";
149  for (size_t i=0; i< fProcCounter->size();i++) {
150  G4String procName = (*fProcCounter)[i]->GetName();
151  G4int count = (*fProcCounter)[i]->GetCounter();
152  G4cout << " " << procName << "= " << count;
153  }
154  G4cout << G4endl;
155 
156  //extract cross sections with G4EmCalculator
157  //
158  G4EmCalculator emCalculator;
159  G4cout << "\n Gamma crossSections in wall material :";
160  G4double sumc = 0.0;
161  for (size_t i=0; i< fProcCounter->size();i++) {
162  G4String procName = (*fProcCounter)[i]->GetName();
163  G4double massSigma =
164  emCalculator.ComputeCrossSectionPerVolume(energy,particle,
165  procName,fMateWall)/fDensityWall;
166  if (massSigma > 0.) {
167  sumc += massSigma;
168  G4cout << " " << procName << "= "
169  << G4BestUnit(massSigma, "Surface/Mass");
170  }
171  }
172  G4cout << " --> total= "
173  << G4BestUnit(sumc, "Surface/Mass") << G4endl;
174 
175  //mean kinetic energy of secondary electrons
176  //
177  if (fNbSec == 0) return;
178  G4double meanEsecond = fEsecondary/fNbSec,meanEsecond2= fEsecondary2/fNbSec;
179  G4double varianceEsec = meanEsecond2 - meanEsecond*meanEsecond;
180  G4double dToverT = 0.;
181  if (varianceEsec>0.) dToverT = std::sqrt(varianceEsec/fNbSec)/meanEsecond;
182  G4double csdaRange =
183  emCalculator.GetCSDARange(meanEsecond,G4Electron::Electron(),fMateWall);
184 
185  G4cout.precision(4);
186  G4cout
187  << "\n Mean energy of secondary e- = " << G4BestUnit(meanEsecond,"Energy")
188  << " +- " << 100*dToverT << " %"
189  << " (--> range in wall material = " << G4BestUnit(csdaRange,"Length")
190  << ")"
191  << G4endl;
192 
193  //compute mass energy transfer coefficient
194  //
195  G4double massTransfCoef = sumc*meanEsecond/energy;
196 
197  G4cout << " Mass_energy_transfer coef: "
198  << G4BestUnit(massTransfCoef, "Surface/Mass")
199  << G4endl;
200 
201  //stopping power from EmCalculator
202  //
203  G4double dedxWall =
204  emCalculator.GetDEDX(meanEsecond,G4Electron::Electron(),fMateWall);
205  dedxWall /= fDensityWall;
206  G4double dedxCavity =
207  emCalculator.GetDEDX(meanEsecond,G4Electron::Electron(),fMateCavity);
208  dedxCavity /= fDensityCavity;
209 
210  G4cout
211  << "\n StoppingPower in wall = "
212  << G4BestUnit(dedxWall,"Energy*Surface/Mass")
213  << "\n in cavity = "
214  << G4BestUnit(dedxCavity,"Energy*Surface/Mass")
215  << G4endl;
216 
217  //charged particle flow in cavity
218  //
219  G4cout
220  << "\n Charged particle flow in cavity :"
221  << "\n Enter --> nbParticles = " << fPartFlowCavity[0]
222  << "\t Energy = " << G4BestUnit (fEnerFlowCavity[0], "Energy")
223  << "\n Exit --> nbParticles = " << fPartFlowCavity[1]
224  << "\t Energy = " << G4BestUnit (fEnerFlowCavity[1], "Energy")
225  << G4endl;
226 
227  if (fPartFlowCavity[0] == 0) return;
228 
229  //beam energy fluence
230  //
232  G4double surfaceBeam = CLHEP::pi*rBeam*rBeam;
233 
234  //error on Edep in cavity
235  //
236  if (fNbEventCavity == 0) return;
239  G4double varianceEdep = meanEdep2 - meanEdep*meanEdep;
240  G4double dEoverE = 0.;
241  if(varianceEdep>0.) dEoverE=std::sqrt(varianceEdep/fNbEventCavity)/meanEdep;
242 
243  //total dose in cavity
244  //
245  G4double doseCavity = fEdepCavity/fMassCavity;
246  G4double doseOverBeam = doseCavity*surfaceBeam/(numberOfEvent*energy);
247 
248  //track length in cavity
249  G4double meantrack = fTrkSegmCavity/fPartFlowCavity[0];
250 
251  G4cout.precision(4);
252  G4cout
253  << "\n Total edep in cavity = " << G4BestUnit(fEdepCavity,"Energy")
254  << " +- " << 100*dEoverE << " %"
255  << "\t Total charged trackLength = " <<G4BestUnit(fTrkSegmCavity,"Length")
256  << " (mean value = " << G4BestUnit(meantrack,"Length") << ")"
257  << "\n Total dose in cavity = " << doseCavity/(MeV/mg) << " MeV/mg"
258  << "\n Dose/EnergyFluence = " << G4BestUnit(doseOverBeam,"Surface/Mass")
259  << G4endl;
260 
261  //ratio simulation/theory
262  //
263  G4double ratio = doseOverBeam/massTransfCoef;
264  G4double error = ratio*std::sqrt(dEoverE*dEoverE + dToverT*dToverT);
265 
266  G4cout.precision(5);
267  G4cout
268  << "\n (Dose/EnergyFluence)/Mass_energy_transfer = " << ratio
269  << " +- " << error << G4endl;
270 
271  //compute mean step size of charged particles
272  //
275  if (rms>0.) rms = std::sqrt(rms); else rms = 0.;
276 
277  G4cout.precision(4);
278  G4cout
279  << "\n StepSize of ch. tracks in wall = "
280  << G4BestUnit(fStepWall,"Length") << " +- " << G4BestUnit( rms,"Length")
281  << "\t (nbSteps/track = " << double(fNbStepWall)/fNbSec << ")";
282 
285  if (rms>0.) rms = std::sqrt(rms); else rms = 0.;
286 
287  G4cout
288  << "\n StepSize of ch. tracks in cavity = "
289  << G4BestUnit(fStepCavity,"Length") << " +- " << G4BestUnit( rms,"Length")
290  << "\t (nbSteps/track = "<<double(fNbStepCavity)/fPartFlowCavity[0] << ")";
291 
292  G4cout << G4endl;
293 
294  // reset default formats
295  G4cout.setf(mode,std::ios::floatfield);
296  G4cout.precision(prec);
297 
298  // delete and remove all contents in fProcCounter
299  while (fProcCounter->size()>0){
300  OneProcessCount* aProcCount=fProcCounter->back();
301  fProcCounter->pop_back();
302  delete aProcCount;
303  }
304  delete fProcCounter;
305 
306 }
307 
308 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
309 
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  //
326  G4double surfaceBeam = CLHEP::pi*rBeam*rBeam;
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 << " ---> 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 << " %) \n"
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......