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 
36 #include "EventAction.hh"
37 #include "HistoManager.hh"
38 #include "PrimaryGeneratorAction.hh"
39 
40 #include "G4Material.hh"
41 #include "G4Event.hh"
42 #include "G4SystemOfUnits.hh"
43 #include "G4UnitsTable.hh"
44 #include <iomanip>
45 
46 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
47 
48 Run::Run(DetectorConstruction* detector)
49 : G4Run(),
50  fDetector(detector),
51  fParticle(0), fEkin(0.)
52 {
53  for (G4int i=0; i<3; ++i) { fStatus[i] = 0; fTotEdep[i] = 0.; }
54  fTotEdep[1] = joule;
55  for (G4int i=0; i<kMaxAbsor; ++i) {
56  fEdeposit[i] = 0.; fEmin[i] = joule; fEmax[i] = 0.;
57  }
58 }
59 
60 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
61 
62 Run::~Run()
63 { }
64 
65 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
66 
67 void Run::SetPrimary (G4ParticleDefinition* particle, G4double energy)
68 {
69  fParticle = particle;
70  fEkin = energy;
71 }
72 
73 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
74 
75 void Run::CountProcesses(const G4VProcess* process)
76 {
77  G4String procName = process->GetProcessName();
78  std::map<G4String,G4int>::iterator it = fProcCounter.find(procName);
79  if ( it == fProcCounter.end()) {
80  fProcCounter[procName] = 1;
81  }
82  else {
83  fProcCounter[procName]++;
84  }
85 }
86 
87 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
88 
89 void Run::ParticleCount(G4int k, G4String name, G4double Ekin)
90 {
91  std::map<G4String, ParticleData>::iterator it = fParticleDataMap[k].find(name);
92  if ( it == fParticleDataMap[k].end()) {
93  (fParticleDataMap[k])[name] = ParticleData(1, Ekin, Ekin, Ekin);
94  }
95  else {
96  ParticleData& data = it->second;
97  data.fCount++;
98  data.fEmean += Ekin;
99  //update min max
100  G4double emin = data.fEmin;
101  if (Ekin < emin) data.fEmin = Ekin;
102  G4double emax = data.fEmax;
103  if (Ekin > emax) data.fEmax = Ekin;
104  }
105 }
106 
107 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
108 
109 void Run::AddEdep (G4int i, G4double e)
110 {
111  if (e > 0.) {
112  fEdeposit[i] += e;
113  if (e < fEmin[i]) fEmin[i] = e;
114  if (e > fEmax[i]) fEmax[i] = e;
115  }
116 }
117 
118 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
119 
120 void Run::AddTotEdep (G4double e)
121 {
122  if (e > 0.) {
123  fTotEdep[0] += e;
124  if (e < fTotEdep[1]) fTotEdep[1] = e;
125  if (e > fTotEdep[2]) fTotEdep[2] = e;
126  }
127 }
128 
129 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
130 
131 void Run::AddTrackStatus (G4int i)
132 {
133  fStatus[i]++ ;
134 }
135 
136 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
137 
138 void Run::Merge(const G4Run* run)
139 {
140  const Run* localRun = static_cast<const Run*>(run);
141 
142  // pass information about primary particle
143  fParticle = localRun->fParticle;
144  fEkin = localRun->fEkin;
145 
146  // Edep in absorbers
147  //
148  G4int nbOfAbsor = fDetector->GetNbOfAbsor();
149  for (G4int i=1; i<=nbOfAbsor; ++i) {
150  fEdeposit[i] += localRun->fEdeposit[i];
151  // min, max
152  G4double min,max;
153  min = localRun->fEmin[i]; max = localRun->fEmax[i];
154  if (fEmin[i] > min) fEmin[i] = min;
155  if (fEmax[i] < max) fEmax[i] = max;
156  }
157 
158  for (G4int i=0; i<3; ++i) fStatus[i] += localRun->fStatus[i];
159 
160  // total Edep
161  fTotEdep[0] += localRun->fTotEdep[0];
162  G4double min,max;
163  min = localRun->fTotEdep[1]; max = localRun->fTotEdep[2];
164  if (fTotEdep[1] > min) fTotEdep[1] = min;
165  if (fTotEdep[2] < max) fTotEdep[2] = max;
166 
167  //map: processes count
168  std::map<G4String,G4int>::const_iterator itp;
169  for ( itp = localRun->fProcCounter.begin();
170  itp != localRun->fProcCounter.end(); ++itp ) {
171 
172  G4String procName = itp->first;
173  G4int localCount = itp->second;
174  if ( fProcCounter.find(procName) == fProcCounter.end()) {
175  fProcCounter[procName] = localCount;
176  }
177  else {
178  fProcCounter[procName] += localCount;
179  }
180  }
181 
182  //map: created particles in absorbers count
183  for (G4int k=0; k<=nbOfAbsor; ++k) {
184  std::map<G4String,ParticleData>::const_iterator itc;
185  for (itc = localRun->fParticleDataMap[k].begin();
186  itc != localRun->fParticleDataMap[k].end(); ++itc) {
187 
188  G4String name = itc->first;
189  const ParticleData& localData = itc->second;
190  if ( fParticleDataMap[k].find(name) == fParticleDataMap[k].end()) {
191  (fParticleDataMap[k])[name]
192  = ParticleData(localData.fCount,
193  localData.fEmean,
194  localData.fEmin,
195  localData.fEmax);
196  }
197  else {
198  ParticleData& data = (fParticleDataMap[k])[name];
199  data.fCount += localData.fCount;
200  data.fEmean += localData.fEmean;
201  G4double emin = localData.fEmin;
202  if (emin < data.fEmin) data.fEmin = emin;
203  G4double emax = localData.fEmax;
204  if (emax > data.fEmax) data.fEmax = emax;
205  }
206  }
207  }
208 
209  G4Run::Merge(run);
210 }
211 
212 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
213 
214 void Run::EndOfRun()
215 {
216  G4int prec = 5, wid = prec + 2;
217  G4int dfprec = G4cout.precision(prec);
218 
219  //run conditions
220  //
221  G4String partName = fParticle->GetParticleName();
222  G4int nbOfAbsor = fDetector->GetNbOfAbsor();
223 
224  G4cout << "\n ======================== run summary =====================\n";
225  G4cout
226  << "\n The run is " << numberOfEvent << " "<< partName << " of "
227  << G4BestUnit(fEkin,"Energy")
228  << " through " << nbOfAbsor << " absorbers: \n";
229  for (G4int i=1; i<= nbOfAbsor; i++) {
230  G4Material* material = fDetector->GetAbsorMaterial(i);
231  G4double thickness = fDetector->GetAbsorThickness(i);
232  G4double density = material->GetDensity();
233  G4cout << std::setw(5) << i
234  << std::setw(10) << G4BestUnit(thickness,"Length") << " of "
235  << material->GetName() << " (density: "
236  << G4BestUnit(density,"Volumic Mass") << ")" << G4endl;
237  }
238 
239  if (numberOfEvent == 0) { G4cout.precision(dfprec); return;}
240 
241  G4cout.precision(3);
242 
243  //frequency of processes
244  //
245  G4cout << "\n Process calls frequency :" << G4endl;
246  G4int index = 0;
247  std::map<G4String,G4int>::iterator it;
248  for (it = fProcCounter.begin(); it != fProcCounter.end(); it++) {
249  G4String procName = it->first;
250  G4int count = it->second;
251  G4String space = " "; if (++index%3 == 0) space = "\n";
252  G4cout << " " << std::setw(20) << procName << "="<< std::setw(7) << count
253  << space;
254  }
255  G4cout << G4endl;
256 
257  //Edep in absorbers
258  //
259  for (G4int i=1; i<= nbOfAbsor; i++) {
260  fEdeposit[i] /= numberOfEvent;
261 
262  G4cout
263  << "\n Edep in absorber " << i << " = "
264  << G4BestUnit(fEdeposit[i],"Energy")
265  << "\t(" << G4BestUnit(fEmin[i], "Energy")
266  << "-->" << G4BestUnit(fEmax[i], "Energy")
267  << ")";
268  }
269  G4cout << G4endl;
270 
271  if (nbOfAbsor > 1) {
272  fTotEdep[0] /= numberOfEvent;
273  G4cout
274  << "\n Edep in all absorbers = " << G4BestUnit(fTotEdep[0],"Energy")
275  << "\t(" << G4BestUnit(fTotEdep[1], "Energy")
276  << "-->" << G4BestUnit(fTotEdep[2], "Energy")
277  << ")" << G4endl;
278  }
279 
280  //particles count in absorbers
281  //
282  for (G4int k=1; k<= nbOfAbsor; k++) {
283  G4cout << "\n List of generated particles in absorber " << k << ":" << G4endl;
284 
285  std::map<G4String,ParticleData>::iterator itc;
286  for (itc = fParticleDataMap[k].begin();
287  itc != fParticleDataMap[k].end(); itc++) {
288  G4String name = itc->first;
289  ParticleData data = itc->second;
290  G4int count = data.fCount;
291  G4double eMean = data.fEmean/count;
292  G4double eMin = data.fEmin;
293  G4double eMax = data.fEmax;
294 
295  G4cout << " " << std::setw(13) << name << ": " << std::setw(7) << count
296  << " Emean = " << std::setw(wid) << G4BestUnit(eMean, "Energy")
297  << "\t( " << G4BestUnit(eMin, "Energy")
298  << " --> " << G4BestUnit(eMax, "Energy")
299  << ")" << G4endl;
300  }
301  }
302  //particles emerging from absorbers
303  //
304  G4cout << "\n List of particles emerging from absorbers :" << G4endl;
305 
306  std::map<G4String,ParticleData>::iterator itc;
307  for (itc = fParticleDataMap[0].begin();
308  itc != fParticleDataMap[0].end(); itc++) {
309  G4String name = itc->first;
310  ParticleData data = itc->second;
311  G4int count = data.fCount;
312  G4double eMean = data.fEmean/count;
313  G4double eMin = data.fEmin;
314  G4double eMax = data.fEmax;
315 
316  G4cout << " " << std::setw(13) << name << ": " << std::setw(7) << count
317  << " Emean = " << std::setw(wid) << G4BestUnit(eMean, "Energy")
318  << "\t( " << G4BestUnit(eMin, "Energy")
319  << " --> " << G4BestUnit(eMax, "Energy")
320  << ")" << G4endl;
321  }
322 
323  //transmission coefficients
324  //
325  G4double dNofEvents = double(numberOfEvent);
326  G4double absorbed = 100.*fStatus[0]/dNofEvents;
327  G4double transmit = 100.*fStatus[1]/dNofEvents;
328  G4double reflected = 100.*fStatus[2]/dNofEvents;
329 
330  G4cout.precision(2);
331  G4cout
332  << "\n Nb of events with primary absorbed = " << absorbed << " %,"
333  << " transmit = " << transmit << " %,"
334  << " reflected = " << reflected << " %" << G4endl;
335 
336  // normalize histograms of longitudinal energy profile
337  //
338  G4AnalysisManager* analysisManager = G4AnalysisManager::Instance();
339  G4int ih = 10;
340  G4double binWidth = analysisManager->GetH1Width(ih)
341  *analysisManager->GetH1Unit(ih);
342  G4double fac = (1./(numberOfEvent*binWidth))*(mm/MeV);
343  analysisManager->ScaleH1(ih,fac);
344 
345  //remove all contents in fProcCounter, fCount
346  fProcCounter.clear();
347  for (G4int k=0; k<= nbOfAbsor; k++) fParticleDataMap[k].clear();
348 
349  // reset default formats
350  G4cout.precision(dfprec);
351 }
352 
353 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......