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 "G4ProcessTable.hh"
39 #include "G4HadronicProcessStore.hh"
40 #include "G4UnitsTable.hh"
41 #include "G4SystemOfUnits.hh"
42 
43 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
44 
45 Run::Run(DetectorConstruction* det)
46 : G4Run(),
47  fDetector(det), fParticle(0), fEkin(0.),
48  fTotalCount(0), fGammaCount(0),
49  fSumTrack(0.), fSumTrack2(0.),
50  fTargetXXX(false)
51 {
52  for (G4int i=0; i<3; i++) { fPbalance[i] = 0. ; }
53  for (G4int i=0; i<3; i++) { fNbGamma[i] = 0 ; }
54  fPbalance[1] = DBL_MAX;
55  fNbGamma[1] = 10000;
56 }
57 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
58 
59 Run::~Run()
60 { }
61 
62 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
63 
64 void Run::SetPrimary(G4ParticleDefinition* particle, G4double energy)
65 {
66  fParticle = particle;
67  fEkin = energy;
68 }
69 
70 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
71 
72 void Run::SetTargetXXX(G4bool flag)
73 {
74  fTargetXXX = flag;
75 }
76 
77 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
78 
79 void Run::CountProcesses(G4VProcess* process)
80 {
81  G4String procName = process->GetProcessName();
82  std::map<G4String,G4int>::iterator it = fProcCounter.find(procName);
83  if ( it == fProcCounter.end()) {
84  fProcCounter[procName] = 1;
85  }
86  else {
87  fProcCounter[procName]++;
88  }
89 }
90 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
91 
92 void Run::SumTrack(G4double trackl)
93 {
94  fTotalCount++;
95  fSumTrack += trackl; fSumTrack2 += trackl*trackl;
96 }
97 
98 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
99 
100 void Run::CountNuclearChannel(G4String name, G4double Q)
101 {
102  std::map<G4String, NuclChannel>::iterator it = fNuclChannelMap.find(name);
103  if ( it == fNuclChannelMap.end()) {
104  fNuclChannelMap[name] = NuclChannel(1, Q);
105  }
106  else {
107  NuclChannel& data = it->second;
108  data.fCount++;
109  data.fQ += Q;
110  }
111 }
112 
113 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
114 
115 void Run::ParticleCount(G4String name, G4double Ekin)
116 {
117  std::map<G4String, ParticleData>::iterator it = fParticleDataMap.find(name);
118  if ( it == fParticleDataMap.end()) {
119  fParticleDataMap[name] = ParticleData(1, Ekin, Ekin, Ekin);
120  }
121  else {
122  ParticleData& data = it->second;
123  data.fCount++;
124  data.fEmean += Ekin;
125  //update min max
126  G4double emin = data.fEmin;
127  if (Ekin < emin) data.fEmin = Ekin;
128  G4double emax = data.fEmax;
129  if (Ekin > emax) data.fEmax = Ekin;
130  }
131 }
132 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
133 
134 void Run::Balance(G4double Pbal)
135 {
136  fPbalance[0] += Pbal;
137  //update min max
138  if (fTotalCount == 1) fPbalance[1] = fPbalance[2] = Pbal;
139  if (Pbal < fPbalance[1]) fPbalance[1] = Pbal;
140  if (Pbal > fPbalance[2]) fPbalance[2] = Pbal;
141 }
142 
143 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
144 
145 void Run::CountGamma(G4int nGamma)
146 {
147  fGammaCount++;
148  fNbGamma[0] += nGamma;
149  //update min max
150  if (fGammaCount == 1) fNbGamma[1] = fNbGamma[2] = nGamma;
151  if (nGamma < fNbGamma[1]) fNbGamma[1] = nGamma;
152  if (nGamma > fNbGamma[2]) fNbGamma[2] = nGamma;
153 }
154 
155 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
156 
157 void Run::Merge(const G4Run* run)
158 {
159  const Run* localRun = static_cast<const Run*>(run);
160 
161  //primary particle info
162  //
163  fParticle = localRun->fParticle;
164  fEkin = localRun->fEkin;
165 
166  // accumulate sums
167  //
168  fTotalCount += localRun->fTotalCount;
169  fGammaCount += localRun->fGammaCount;
170  fSumTrack += localRun->fSumTrack;
171  fSumTrack2 += localRun->fSumTrack2;
172 
173  fPbalance[0] += localRun->fPbalance[0];
174  G4double min,max;
175  min = localRun->fPbalance[1]; max = localRun->fPbalance[2];
176  if (fPbalance[1] > min) fPbalance[1] = min;
177  if (fPbalance[2] < max) fPbalance[2] = max;
178 
179  fNbGamma[0] += localRun->fNbGamma[0];
180  G4int nbmin, nbmax;
181  nbmin = localRun->fNbGamma[1]; nbmax = localRun->fNbGamma[2];
182  if (fNbGamma[1] > nbmin) fNbGamma[1] = nbmin;
183  if (fNbGamma[2] < nbmax) fNbGamma[2] = nbmax;
184 
185  //map: processes count
186  std::map<G4String,G4int>::const_iterator itp;
187  for ( itp = localRun->fProcCounter.begin();
188  itp != localRun->fProcCounter.end(); ++itp ) {
189 
190  G4String procName = itp->first;
191  G4int localCount = itp->second;
192  if ( fProcCounter.find(procName) == fProcCounter.end()) {
193  fProcCounter[procName] = localCount;
194  }
195  else {
196  fProcCounter[procName] += localCount;
197  }
198  }
199 
200  //map: nuclear channels
201  std::map<G4String,NuclChannel>::const_iterator itc;
202  for (itc = localRun->fNuclChannelMap.begin();
203  itc != localRun->fNuclChannelMap.end(); ++itc) {
204 
205  G4String name = itc->first;
206  const NuclChannel& localData = itc->second;
207  if ( fNuclChannelMap.find(name) == fNuclChannelMap.end()) {
208  fNuclChannelMap[name]
209  = NuclChannel(localData.fCount, localData.fQ);
210  }
211  else {
212  NuclChannel& data = fNuclChannelMap[name];
213  data.fCount += localData.fCount;
214  data.fQ += localData.fQ;
215  }
216  }
217 
218  //map: particles count
219  std::map<G4String,ParticleData>::const_iterator itn;
220  for (itn = localRun->fParticleDataMap.begin();
221  itn != localRun->fParticleDataMap.end(); ++itn) {
222 
223  G4String name = itn->first;
224  const ParticleData& localData = itn->second;
225  if ( fParticleDataMap.find(name) == fParticleDataMap.end()) {
226  fParticleDataMap[name]
227  = ParticleData(localData.fCount,
228  localData.fEmean,
229  localData.fEmin,
230  localData.fEmax);
231  }
232  else {
233  ParticleData& data = fParticleDataMap[name];
234  data.fCount += localData.fCount;
235  data.fEmean += localData.fEmean;
236  G4double emin = localData.fEmin;
237  if (emin < data.fEmin) data.fEmin = emin;
238  G4double emax = localData.fEmax;
239  if (emax > data.fEmax) data.fEmax = emax;
240  }
241  }
242 
243  G4Run::Merge(run);
244 }
245 
246 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
247 
248 void Run::EndOfRun(G4bool print)
249 {
250  G4int prec = 5, wid = prec + 2;
251  G4int dfprec = G4cout.precision(prec);
252 
253  //run condition
254  //
255  G4Material* material = fDetector->GetMaterial();
256  G4double density = material->GetDensity();
257 
258  G4String Particle = fParticle->GetParticleName();
259  G4cout << "\n The run is " << numberOfEvent << " "<< Particle << " of "
260  << G4BestUnit(fEkin,"Energy") << " through "
261  << G4BestUnit(fDetector->GetSize(),"Length") << " of "
262  << material->GetName() << " (density: "
263  << G4BestUnit(density,"Volumic Mass") << ")" << G4endl;
264 
265  if (numberOfEvent == 0) { G4cout.precision(dfprec); return;}
266 
267  //frequency of processes
268  //
269  G4cout << "\n Process calls frequency:" << G4endl;
270  G4int survive = 0;
271  std::map<G4String,G4int>::iterator it;
272  for (it = fProcCounter.begin(); it != fProcCounter.end(); it++) {
273  G4String procName = it->first;
274  G4int count = it->second;
275  G4cout << "\t" << procName << "= " << count;
276  if (procName == "Transportation") survive = count;
277  }
278  G4cout << G4endl;
279 
280  if (survive > 0) {
281  G4cout << "\n Nb of incident particles surviving after "
282  << G4BestUnit(fDetector->GetSize(),"Length") << " of "
283  << material->GetName() << " : " << survive << G4endl;
284  }
285 
286  if (fTotalCount == 0) fTotalCount = 1; //force printing anyway
287 
288  //compute mean free path and related quantities
289  //
290  G4double MeanFreePath = fSumTrack /fTotalCount;
291  G4double MeanTrack2 = fSumTrack2/fTotalCount;
292  G4double rms = std::sqrt(std::fabs(MeanTrack2 - MeanFreePath*MeanFreePath));
293  G4double CrossSection = 0.0;
294  if(MeanFreePath > 0.0) { CrossSection = 1./MeanFreePath; }
295  G4double massicMFP = MeanFreePath*density;
296  G4double massicCS = 0.0;
297  if(massicMFP > 0.0) { massicCS = 1./massicMFP; }
298 
299  G4cout << "\n\n MeanFreePath:\t" << G4BestUnit(MeanFreePath,"Length")
300  << " +- " << G4BestUnit( rms,"Length")
301  << "\tmassic: " << G4BestUnit(massicMFP, "Mass/Surface")
302  << "\n CrossSection:\t" << CrossSection*cm << " cm^-1 "
303  << "\t\tmassic: " << G4BestUnit(massicCS, "Surface/Mass")
304  << G4endl;
305 
306  //cross section per atom (only for single material)
307  //
308  if (material->GetNumberOfElements() == 1) {
309  G4double nbAtoms = material->GetTotNbOfAtomsPerVolume();
310  G4double crossSection = CrossSection/nbAtoms;
311  G4cout << " crossSection per atom:\t"
312  << G4BestUnit(crossSection,"Surface") << G4endl;
313  }
314  //check cross section from G4HadronicProcessStore
315  //
316  G4cout << "\n Verification: "
317  << "crossSections from G4HadronicProcessStore:";
318 
319  G4ProcessTable* processTable = G4ProcessTable::GetProcessTable();
320  G4HadronicProcessStore* store = G4HadronicProcessStore::Instance();
321  G4double sumc1 = 0.0, sumc2 = 0.0;
322  if (material->GetNumberOfElements() == 1) {
323  const G4Element* element = material->GetElement(0);
324  for (it = fProcCounter.begin(); it != fProcCounter.end(); it++) {
325  G4String procName = it->first;
326  G4VProcess* process = processTable->FindProcess(procName, fParticle);
327  G4double xs1 =
328  store->GetCrossSectionPerVolume(fParticle,fEkin,process,material);
329  G4double massSigma = xs1/density;
330  sumc1 += massSigma;
331  G4double xs2 =
332  store->GetCrossSectionPerAtom(fParticle,fEkin,process,element,material);
333  sumc2 += xs2;
334  G4cout << "\n" << std::setw(20) << procName << "= "
335  << G4BestUnit(massSigma, "Surface/Mass") << "\t"
336  << G4BestUnit(xs2, "Surface");
337 
338  }
339  G4cout << "\n" << std::setw(20) << "total" << "= "
340  << G4BestUnit(sumc1, "Surface/Mass") << "\t"
341  << G4BestUnit(sumc2, "Surface") << G4endl;
342  } else {
343  for (it = fProcCounter.begin(); it != fProcCounter.end(); it++) {
344  G4String procName = it->first;
345  G4VProcess* process = processTable->FindProcess(procName, fParticle);
346  G4double xs =
347  store->GetCrossSectionPerVolume(fParticle,fEkin,process,material);
348  G4double massSigma = xs/density;
349  sumc1 += massSigma;
350  G4cout << "\n" << std::setw(20) << procName << "= "
351  << G4BestUnit(massSigma, "Surface/Mass");
352  }
353  G4cout << "\n" << std::setw(20) << "total" << "= "
354  << G4BestUnit(sumc1, "Surface/Mass") << G4endl;
355  }
356 
357  //nuclear channel count
358  //
359  G4cout << "\n List of nuclear reactions: \n" << G4endl;
360  std::map<G4String,NuclChannel>::iterator ic;
361  for (ic = fNuclChannelMap.begin(); ic != fNuclChannelMap.end(); ic++) {
362  G4String name = ic->first;
363  NuclChannel data = ic->second;
364  G4int count = data.fCount;
365  G4double Q = data.fQ/count;
366  if (print)
367  G4cout << " " << std::setw(60) << name << ": " << std::setw(7) << count
368  << " Q = " << std::setw(wid) << G4BestUnit(Q, "Energy")
369  << G4endl;
370  }
371 
372  //Gamma count
373  //
374  if (print && (fGammaCount > 0)) {
375  G4cout << "\n" << std::setw(58) << "number of gamma or e- (ic): N = "
376  << fNbGamma[1] << " --> " << fNbGamma[2] << G4endl;
377  }
378 
379  if (print && fTargetXXX) {
380  G4cout
381  << "\n --> NOTE: XXXX because neutronHP is unable to return target nucleus"
382  << G4endl;
383  }
384 
385  //particles count
386  //
387  G4cout << "\n List of generated particles:" << G4endl;
388 
389  std::map<G4String,ParticleData>::iterator itn;
390  for (itn = fParticleDataMap.begin(); itn != fParticleDataMap.end(); itn++) {
391  G4String name = itn->first;
392  ParticleData data = itn->second;
393  G4int count = data.fCount;
394  G4double eMean = data.fEmean/count;
395  G4double eMin = data.fEmin;
396  G4double eMax = data.fEmax;
397  if (print)
398  G4cout << " " << std::setw(13) << name << ": " << std::setw(7) << count
399  << " Emean = " << std::setw(wid) << G4BestUnit(eMean, "Energy")
400  << "\t( " << G4BestUnit(eMin, "Energy")
401  << " --> " << G4BestUnit(eMax, "Energy")
402  << ")" << G4endl;
403  }
404 
405  //energy momentum balance
406  //
407  if (fTotalCount > 1) {
408  G4double Pbmean = fPbalance[0]/fTotalCount;
409  G4cout << "\n Momentum balance: Pmean = "
410  << std::setw(wid) << G4BestUnit(Pbmean, "Energy")
411  << "\t( " << G4BestUnit(fPbalance[1], "Energy")
412  << " --> " << G4BestUnit(fPbalance[2], "Energy")
413  << ") \n" << G4endl;
414  }
415 
416  //normalize histograms
420 
421  //remove all contents in fProcCounter, fCount
422  fProcCounter.clear();
423  fNuclChannelMap.clear();
424  fParticleDataMap.clear();
425 
426  //restore default format
427  G4cout.precision(dfprec);
428 }
429 
430 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......