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->GetMaterial();
229  G4double density = material->GetDensity();
230 
231  G4String Particle = fParticle->GetParticleName();
232  G4cout << "\n The run is " << numberOfEvent << " "<< Particle << " of "
233  << G4BestUnit(fEkin,"Energy") << " through "
234  << G4BestUnit(fDetector->GetRadius(),"Length") << " of "
235  << material->GetName() << " (density: "
236  << G4BestUnit(density,"Volumic Mass") << ")" << G4endl;
237 
238  if (numberOfEvent == 0) { G4cout.precision(dfprec); return;}
239 
240  //frequency of processes
241  //
242  G4cout << "\n Process calls frequency :" << G4endl;
243  G4int index = 0;
244  std::map<G4String,G4int>::iterator it;
245  for (it = fProcCounter.begin(); it != fProcCounter.end(); it++) {
246  G4String procName = it->first;
247  G4int count = it->second;
248  G4String space = " "; if (++index%3 == 0) space = "\n";
249  G4cout << " " << std::setw(20) << procName << "="<< std::setw(7) << count
250  << space;
251  }
252  G4cout << G4endl;
253 
254  //particles count
255  //
256  G4cout << "\n List of generated particles:" << G4endl;
257 
258  std::map<G4String,ParticleData>::iterator itc;
259  for (itc = fParticleDataMap1.begin(); itc != fParticleDataMap1.end(); itc++) {
260  G4String name = itc->first;
261  ParticleData data = itc->second;
262  G4int count = data.fCount;
263  G4double eMean = data.fEmean/count;
264  G4double eMin = data.fEmin;
265  G4double eMax = data.fEmax;
266 
267  G4cout << " " << std::setw(13) << name << ": " << std::setw(7) << count
268  << " Emean = " << std::setw(wid) << G4BestUnit(eMean, "Energy")
269  << "\t( " << G4BestUnit(eMin, "Energy")
270  << " --> " << G4BestUnit(eMax, "Energy")
271  << ")" << G4endl;
272  }
273 
274  // compute mean Energy deposited and rms
275  //
276  G4int TotNbofEvents = numberOfEvent;
277  fEnergyDeposit /= TotNbofEvents; fEnergyDeposit2 /= TotNbofEvents;
278  G4double rmsEdep = fEnergyDeposit2 - fEnergyDeposit*fEnergyDeposit;
279  if (rmsEdep>0.) rmsEdep = std::sqrt(rmsEdep);
280  else rmsEdep = 0.;
281 
282  G4cout << "\n Mean energy deposit per event = "
283  << G4BestUnit(fEnergyDeposit,"Energy") << "; rms = "
284  << G4BestUnit(rmsEdep, "Energy")
285  << G4endl;
286 
287  // compute mean Energy flow and rms
288  //
289  fEnergyFlow /= TotNbofEvents; fEnergyFlow2 /= TotNbofEvents;
290  G4double rmsEflow = fEnergyFlow2 - fEnergyFlow*fEnergyFlow;
291  if (rmsEflow>0.) rmsEflow = std::sqrt(rmsEflow);
292  else rmsEflow = 0.;
293 
294  G4cout << " Mean energy flow per event = "
295  << G4BestUnit(fEnergyFlow,"Energy") << "; rms = "
296  << G4BestUnit(rmsEflow, "Energy")
297  << G4endl;
298 
299  //particles flux
300  //
301  G4cout << "\n List of particles emerging from the absorber :" << G4endl;
302 
303  std::map<G4String,ParticleData>::iterator itn;
304  for (itn = fParticleDataMap2.begin(); itn != fParticleDataMap2.end(); itn++) {
305  G4String name = itn->first;
306  ParticleData data = itn->second;
307  G4int count = data.fCount;
308  G4double eMean = data.fEmean/count;
309  G4double eMin = data.fEmin;
310  G4double eMax = data.fEmax;
311  G4double Eflow = data.fEmean/TotNbofEvents;
312 
313  G4cout << " " << std::setw(13) << name << ": " << std::setw(7) << count
314  << " Emean = " << std::setw(wid) << G4BestUnit(eMean, "Energy")
315  << "\t( " << G4BestUnit(eMin, "Energy")
316  << " --> " << G4BestUnit(eMax, "Energy")
317  << ") \tEflow/event = " << G4BestUnit(Eflow, "Energy") << G4endl;
318  }
319 
320  //normalize histograms
321  G4AnalysisManager* analysisManager = G4AnalysisManager::Instance();
322 
323  G4int ih = 2;
324  G4double binWidth = analysisManager->GetH1Width(ih);
325  G4double fac = (1./(numberOfEvent*binWidth))*(mm/MeV);
326  analysisManager->ScaleH1(ih,fac);
327 
328  //remove all contents in fProcCounter, fCount
329  fProcCounter.clear();
330  fParticleDataMap2.clear();
331 
332  //restore default format
333  G4cout.precision(dfprec);
334 }
335 
336 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......