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