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 "G4UnitsTable.hh"
39 #include "G4SystemOfUnits.hh"
40 
41 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
42 
43 Run::Run(DetectorConstruction* det)
44 : G4Run(),
45  fDetector(det), fParticle(0), fEkin(0.)
46 {
47  fEnergyDeposit = fEnergyDeposit2 = 0.;
48  fEnergyFlow = fEnergyFlow2 = 0.;
49 }
50 
51 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
52 
53 Run::~Run()
54 { }
55 
56 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
57 
58 void Run::SetPrimary(G4ParticleDefinition* particle, G4double energy)
59 {
60  fParticle = particle;
61  fEkin = energy;
62 }
63 
64 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
65 
66 void Run::CountProcesses(const G4VProcess* process)
67 {
68  G4String procName = process->GetProcessName();
69  std::map<G4String,G4int>::iterator it = fProcCounter.find(procName);
70  if ( it == fProcCounter.end()) {
71  fProcCounter[procName] = 1;
72  }
73  else {
74  fProcCounter[procName]++;
75  }
76 }
77 
78 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
79 
80 void Run::ParticleCount(G4String name, G4double Ekin)
81 {
82  std::map<G4String, ParticleData>::iterator it = fParticleDataMap1.find(name);
83  if ( it == fParticleDataMap1.end()) {
84  fParticleDataMap1[name] = ParticleData(1, Ekin, Ekin, Ekin);
85  }
86  else {
87  ParticleData& data = it->second;
88  data.fCount++;
89  data.fEmean += Ekin;
90  //update min max
91  G4double emin = data.fEmin;
92  if (Ekin < emin) data.fEmin = Ekin;
93  G4double emax = data.fEmax;
94  if (Ekin > emax) data.fEmax = Ekin;
95  }
96 }
97 
98 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
99 
100 void Run::AddEdep(G4double edep)
101 {
102  fEnergyDeposit += edep;
103  fEnergyDeposit2 += edep*edep;
104 }
105 
106 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
107 
108 void Run::AddEflow(G4double eflow)
109 {
110  fEnergyFlow += eflow;
111  fEnergyFlow2 += eflow*eflow;
112 }
113 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
114 
115 void Run::ParticleFlux(G4String name, G4double Ekin)
116 {
117  std::map<G4String, ParticleData>::iterator it = fParticleDataMap2.find(name);
118  if ( it == fParticleDataMap2.end()) {
119  fParticleDataMap2[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 
133 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
134 
135 void Run::Merge(const G4Run* run)
136 {
137  const Run* localRun = static_cast<const Run*>(run);
138 
139  //primary particle info
140  //
141  fParticle = localRun->fParticle;
142  fEkin = localRun->fEkin;
143 
144  // accumulate sums
145  //
146  fEnergyDeposit += localRun->fEnergyDeposit;
147  fEnergyDeposit2 += localRun->fEnergyDeposit2;
148  fEnergyFlow += localRun->fEnergyFlow;
149  fEnergyFlow2 += localRun->fEnergyFlow2;
150 
151  //map: processes count
152  std::map<G4String,G4int>::const_iterator itp;
153  for ( itp = localRun->fProcCounter.begin();
154  itp != localRun->fProcCounter.end(); ++itp ) {
155 
156  G4String procName = itp->first;
157  G4int localCount = itp->second;
158  if ( fProcCounter.find(procName) == fProcCounter.end()) {
159  fProcCounter[procName] = localCount;
160  }
161  else {
162  fProcCounter[procName] += localCount;
163  }
164  }
165 
166  //map: created particles count
167  std::map<G4String,ParticleData>::const_iterator itc;
168  for (itc = localRun->fParticleDataMap1.begin();
169  itc != localRun->fParticleDataMap1.end(); ++itc) {
170 
171  G4String name = itc->first;
172  const ParticleData& localData = itc->second;
173  if ( fParticleDataMap1.find(name) == fParticleDataMap1.end()) {
174  fParticleDataMap1[name]
175  = ParticleData(localData.fCount,
176  localData.fEmean,
177  localData.fEmin,
178  localData.fEmax);
179  }
180  else {
181  ParticleData& data = fParticleDataMap1[name];
182  data.fCount += localData.fCount;
183  data.fEmean += localData.fEmean;
184  G4double emin = localData.fEmin;
185  if (emin < data.fEmin) data.fEmin = emin;
186  G4double emax = localData.fEmax;
187  if (emax > data.fEmax) data.fEmax = emax;
188  }
189  }
190 
191  //map: particles flux count
192  std::map<G4String,ParticleData>::const_iterator itn;
193  for (itn = localRun->fParticleDataMap2.begin();
194  itn != localRun->fParticleDataMap2.end(); ++itn) {
195 
196  G4String name = itn->first;
197  const ParticleData& localData = itn->second;
198  if ( fParticleDataMap2.find(name) == fParticleDataMap2.end()) {
199  fParticleDataMap2[name]
200  = ParticleData(localData.fCount,
201  localData.fEmean,
202  localData.fEmin,
203  localData.fEmax);
204  }
205  else {
206  ParticleData& data = fParticleDataMap2[name];
207  data.fCount += localData.fCount;
208  data.fEmean += localData.fEmean;
209  G4double emin = localData.fEmin;
210  if (emin < data.fEmin) data.fEmin = emin;
211  G4double emax = localData.fEmax;
212  if (emax > data.fEmax) data.fEmax = emax;
213  }
214  }
215 
216  G4Run::Merge(run);
217 }
218 
219 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
220 
221 void Run::EndOfRun()
222 {
223  G4int prec = 5, wid = prec + 2;
224  G4int dfprec = G4cout.precision(prec);
225 
226  //run condition
227  //
228  G4Material* material = fDetector->GetAbsorMaterial();
229  G4String Particle = fParticle->GetParticleName();
230  G4cout << "\n The run is " << numberOfEvent << " "<< Particle << " of "
231  << G4BestUnit(fEkin,"Energy") << " within "
232  << material->GetName() << " (D = "
233  << G4BestUnit(2*(fDetector->GetAbsorRadius()),"Length") << " L = "
234  << G4BestUnit(fDetector->GetAbsorLength(),"Length") << ")" << G4endl;
235 
236  if (numberOfEvent == 0) { G4cout.precision(dfprec); return;}
237 
238  //frequency of processes
239  //
240  G4cout << "\n Process calls frequency :" << G4endl;
241  G4int index = 0;
242  std::map<G4String,G4int>::iterator it;
243  for (it = fProcCounter.begin(); it != fProcCounter.end(); it++) {
244  G4String procName = it->first;
245  G4int count = it->second;
246  G4String space = " "; if (++index%3 == 0) space = "\n";
247  G4cout << " " << std::setw(20) << procName << "="<< std::setw(7) << count
248  << space;
249  }
250  G4cout << G4endl;
251 
252  //particles count
253  //
254  G4cout << "\n List of generated particles:" << G4endl;
255 
256  std::map<G4String,ParticleData>::iterator itc;
257  for (itc = fParticleDataMap1.begin(); itc != fParticleDataMap1.end(); itc++) {
258  G4String name = itc->first;
259  ParticleData data = itc->second;
260  G4int count = data.fCount;
261  G4double eMean = data.fEmean/count;
262  G4double eMin = data.fEmin;
263  G4double eMax = data.fEmax;
264 
265  G4cout << " " << std::setw(13) << name << ": " << std::setw(7) << count
266  << " Emean = " << std::setw(wid) << G4BestUnit(eMean, "Energy")
267  << "\t( " << G4BestUnit(eMin, "Energy")
268  << " --> " << G4BestUnit(eMax, "Energy")
269  << ")" << G4endl;
270  }
271 
272  // compute mean Energy deposited and rms
273  //
274  G4int TotNbofEvents = numberOfEvent;
275  fEnergyDeposit /= TotNbofEvents; fEnergyDeposit2 /= TotNbofEvents;
276  G4double rmsEdep = fEnergyDeposit2 - fEnergyDeposit*fEnergyDeposit;
277  if (rmsEdep>0.) rmsEdep = std::sqrt(rmsEdep);
278  else rmsEdep = 0.;
279 
280  G4cout << "\n Mean energy deposit per event = "
281  << G4BestUnit(fEnergyDeposit,"Energy") << "; rms = "
282  << G4BestUnit(rmsEdep, "Energy")
283  << G4endl;
284 
285  // compute mean Energy flow and rms
286  //
287  fEnergyFlow /= TotNbofEvents; fEnergyFlow2 /= TotNbofEvents;
288  G4double rmsEflow = fEnergyFlow2 - fEnergyFlow*fEnergyFlow;
289  if (rmsEflow>0.) rmsEflow = std::sqrt(rmsEflow);
290  else rmsEflow = 0.;
291 
292  G4cout << " Mean energy flow per event = "
293  << G4BestUnit(fEnergyFlow,"Energy") << "; rms = "
294  << G4BestUnit(rmsEflow, "Energy")
295  << G4endl;
296 
297  //particles flux
298  //
299  G4cout << "\n List of particles emerging from the container :" << G4endl;
300 
301  std::map<G4String,ParticleData>::iterator itn;
302  for (itn = fParticleDataMap2.begin(); itn != fParticleDataMap2.end(); itn++) {
303  G4String name = itn->first;
304  ParticleData data = itn->second;
305  G4int count = data.fCount;
306  G4double eMean = data.fEmean/count;
307  G4double eMin = data.fEmin;
308  G4double eMax = data.fEmax;
309  G4double Eflow = data.fEmean/TotNbofEvents;
310 
311  G4cout << " " << std::setw(13) << name << ": " << std::setw(7) << count
312  << " Emean = " << std::setw(wid) << G4BestUnit(eMean, "Energy")
313  << "\t( " << G4BestUnit(eMin, "Energy")
314  << " --> " << G4BestUnit(eMax, "Energy")
315  << ") \tEflow/event = " << G4BestUnit(Eflow, "Energy") << G4endl;
316  }
317 
318  //remove all contents in fProcCounter, fCount
319  fProcCounter.clear();
320  fParticleDataMap2.clear();
321 
322  //restore default format
323  G4cout.precision(dfprec);
324 }
325 
326 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......