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 "G4Threading.hh"
39 #include "G4AutoLock.hh"
40 #include "G4UnitsTable.hh"
41 #include "G4SystemOfUnits.hh"
42 
43 // mutex in a file scope
44 
45 namespace {
46  //Mutex to lock updating the global ion map
47  G4Mutex ionIdMapMutex = G4MUTEX_INITIALIZER;
48 }
49 
50 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
51 
52 std::map<G4String,G4int> Run::fgIonMap;
53 G4int Run::fgIonId = kMaxHisto1;
54 
55 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
56 
57 Run::Run(DetectorConstruction* det)
58 : G4Run(),
59  fDetector(det), fParticle(nullptr), fEkin(0.)
60 {
61  fEnergyDeposit = fEnergyDeposit2 = 0.;
62  fEnergyFlow = fEnergyFlow2 = 0.;
63 }
64 
65 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
66 
67 Run::~Run()
68 { }
69 
70 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
71 
72 void Run::Merge(std::map<G4String, ParticleData>& destinationMap,
73  const std::map<G4String, ParticleData>& sourceMap) const
74 {
75  for ( const auto& particleData : sourceMap ) {
76  G4String name = particleData.first;
77  const ParticleData& localData = particleData.second;
78  if ( destinationMap.find(name) == destinationMap.end()) {
79  destinationMap[name]
80  = ParticleData(localData.fCount,
81  localData.fEmean,
82  localData.fEmin,
83  localData.fEmax,
84  localData.fTmean);
85  }
86  else {
87  ParticleData& data = destinationMap[name];
88  data.fCount += localData.fCount;
89  data.fEmean += localData.fEmean;
90  G4double emin = localData.fEmin;
91  if (emin < data.fEmin) data.fEmin = emin;
92  G4double emax = localData.fEmax;
93  if (emax > data.fEmax) data.fEmax = emax;
94  data.fTmean = localData.fTmean;
95  }
96  }
97 }
98 
99 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
100 
101 void Run::SetPrimary(G4ParticleDefinition* particle, G4double energy)
102 {
103  fParticle = particle;
104  fEkin = energy;
105 }
106 
107 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
108 
109 void Run::CountProcesses(const G4VProcess* process)
110 {
111  G4String procName = process->GetProcessName();
112  std::map<G4String,G4int>::iterator it = fProcCounter.find(procName);
113  if ( it == fProcCounter.end()) {
114  fProcCounter[procName] = 1;
115  }
116  else {
117  fProcCounter[procName]++;
118  }
119 }
120 
121 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
122 
123 void Run::ParticleCount(G4String name, G4double Ekin, G4double meanLife)
124 {
125  std::map<G4String, ParticleData>::iterator it = fParticleDataMap1.find(name);
126  if ( it == fParticleDataMap1.end()) {
127  fParticleDataMap1[name] = ParticleData(1, Ekin, Ekin, Ekin, meanLife);
128  }
129  else {
130  ParticleData& data = it->second;
131  data.fCount++;
132  data.fEmean += Ekin;
133  //update min max
134  G4double emin = data.fEmin;
135  if (Ekin < emin) data.fEmin = Ekin;
136  G4double emax = data.fEmax;
137  if (Ekin > emax) data.fEmax = Ekin;
138  data.fTmean = meanLife;
139  }
140 }
141 
142 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
143 
144 void Run::AddEdep(G4double edep)
145 {
146  fEnergyDeposit += edep;
147  fEnergyDeposit2 += edep*edep;
148 }
149 
150 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
151 
152 void Run::AddEflow(G4double eflow)
153 {
154  fEnergyFlow += eflow;
155  fEnergyFlow2 += eflow*eflow;
156 }
157 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
158 
159 void Run::ParticleFlux(G4String name, G4double Ekin)
160 {
161  std::map<G4String, ParticleData>::iterator it = fParticleDataMap2.find(name);
162  if ( it == fParticleDataMap2.end()) {
163  fParticleDataMap2[name] = ParticleData(1, Ekin, Ekin, Ekin, -1*ns);
164  }
165  else {
166  ParticleData& data = it->second;
167  data.fCount++;
168  data.fEmean += Ekin;
169  //update min max
170  G4double emin = data.fEmin;
171  if (Ekin < emin) data.fEmin = Ekin;
172  G4double emax = data.fEmax;
173  if (Ekin > emax) data.fEmax = Ekin;
174  data.fTmean = -1*ns;
175  }
176 }
177 
178 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
179 
180 G4int Run::GetIonId(G4String ionName)
181 {
182  G4AutoLock lock(&ionIdMapMutex);
183  // updating the global ion map needs to be locked
184 
185  std::map<G4String,G4int>::const_iterator it = fgIonMap.find(ionName);
186  if ( it == fgIonMap.end()) {
187  fgIonMap[ionName] = fgIonId;
188  if (fgIonId < kMaxHisto2) fgIonId++;
189  }
190  return fgIonMap[ionName];
191 }
192 
193 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
194 
195 void Run::Merge(const G4Run* run)
196 {
197  const Run* localRun = static_cast<const Run*>(run);
198 
199  //primary particle info
200  //
201  fParticle = localRun->fParticle;
202  fEkin = localRun->fEkin;
203 
204  // accumulate sums
205  //
206  fEnergyDeposit += localRun->fEnergyDeposit;
207  fEnergyDeposit2 += localRun->fEnergyDeposit2;
208  fEnergyFlow += localRun->fEnergyFlow;
209  fEnergyFlow2 += localRun->fEnergyFlow2;
210 
211  //map: processes count
212  for ( const auto& procCounter : localRun->fProcCounter ) {
213  G4String procName = procCounter.first;
214  G4int localCount = procCounter.second;
215  if ( fProcCounter.find(procName) == fProcCounter.end()) {
216  fProcCounter[procName] = localCount;
217  }
218  else {
219  fProcCounter[procName] += localCount;
220  }
221  }
222 
223  //map: created particles count
224  Merge(fParticleDataMap1, localRun->fParticleDataMap1);
225 
226  //map: particles flux count
227  Merge(fParticleDataMap2, localRun->fParticleDataMap2);
228 
229  G4Run::Merge(run);
230 }
231 
232 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
233 
234 void Run::EndOfRun()
235 {
236  G4int prec = 5, wid = prec + 2;
237  G4int dfprec = G4cout.precision(prec);
238 
239  //run condition
240  //
241  G4Material* material = fDetector->GetAbsorMaterial();
242  G4double density = material->GetDensity();
243 
244  G4String Particle = fParticle->GetParticleName();
245  G4cout << "\n The run is " << numberOfEvent << " "<< Particle << " of "
246  << G4BestUnit(fEkin,"Energy") << " through "
247  << G4BestUnit(fDetector->GetAbsorThickness(),"Length") << " of "
248  << material->GetName() << " (density: "
249  << G4BestUnit(density,"Volumic Mass") << ")" << G4endl;
250 
251  if (numberOfEvent == 0) { G4cout.precision(dfprec); return;}
252 
253  //frequency of processes
254  //
255  G4cout << "\n Process calls frequency :" << G4endl;
256  G4int index = 0;
257  for ( const auto& procCounter : fProcCounter ) {
258  G4String procName = procCounter.first;
259  G4int count = procCounter.second;
260  G4String space = " "; if (++index%3 == 0) space = "\n";
261  G4cout << " " << std::setw(20) << procName << "="<< std::setw(7) << count
262  << space;
263  }
264  G4cout << G4endl;
265 
266  //particles count
267  //
268  G4cout << "\n List of generated particles:" << G4endl;
269 
270  for ( const auto& particleData : fParticleDataMap1 ) {
271  G4String name = particleData.first;
272  ParticleData data = particleData.second;
273  G4int count = data.fCount;
274  G4double eMean = data.fEmean/count;
275  G4double eMin = data.fEmin;
276  G4double eMax = data.fEmax;
277  G4double meanLife = data.fTmean;
278 
279  G4cout << " " << std::setw(13) << name << ": " << std::setw(7) << count
280  << " Emean = " << std::setw(wid) << G4BestUnit(eMean, "Energy")
281  << "\t( " << G4BestUnit(eMin, "Energy")
282  << " --> " << G4BestUnit(eMax, "Energy") << ")";
283  if (meanLife >= 0.)
284  G4cout << "\tmean life = " << G4BestUnit(meanLife, "Time") << G4endl;
285  else G4cout << "\tstable" << G4endl;
286  }
287 
288  // compute mean Energy deposited and rms
289  //
290  G4int TotNbofEvents = numberOfEvent;
291  fEnergyDeposit /= TotNbofEvents; fEnergyDeposit2 /= TotNbofEvents;
292  G4double rmsEdep = fEnergyDeposit2 - fEnergyDeposit*fEnergyDeposit;
293  if (rmsEdep>0.) rmsEdep = std::sqrt(rmsEdep);
294  else rmsEdep = 0.;
295 
296  G4cout << "\n Mean energy deposit per event = "
297  << G4BestUnit(fEnergyDeposit,"Energy") << "; rms = "
298  << G4BestUnit(rmsEdep, "Energy")
299  << G4endl;
300 
301  // compute mean Energy flow and rms
302  //
303  fEnergyFlow /= TotNbofEvents; fEnergyFlow2 /= TotNbofEvents;
304  G4double rmsEflow = fEnergyFlow2 - fEnergyFlow*fEnergyFlow;
305  if (rmsEflow>0.) rmsEflow = std::sqrt(rmsEflow);
306  else rmsEflow = 0.;
307 
308  G4cout << " Mean energy flow per event = "
309  << G4BestUnit(fEnergyFlow,"Energy") << "; rms = "
310  << G4BestUnit(rmsEflow, "Energy")
311  << G4endl;
312 
313  //particles flux
314  //
315  G4cout << "\n List of particles emerging from the target :" << G4endl;
316 
317  for ( const auto& particleData : fParticleDataMap2 ) {
318  G4String name = particleData.first;
319  ParticleData data = particleData.second;
320  G4int count = data.fCount;
321  G4double eMean = data.fEmean/count;
322  G4double eMin = data.fEmin;
323  G4double eMax = data.fEmax;
324  G4double Eflow = data.fEmean/TotNbofEvents;
325 
326  G4cout << " " << std::setw(13) << name << ": " << std::setw(7) << count
327  << " Emean = " << std::setw(wid) << G4BestUnit(eMean, "Energy")
328  << "\t( " << G4BestUnit(eMin, "Energy")
329  << " --> " << G4BestUnit(eMax, "Energy")
330  << ") \tEflow/event = " << G4BestUnit(Eflow, "Energy") << G4endl;
331  }
332 
333  //histogram Id for populations
334  //
335  G4cout << "\n histo Id for populations :" << G4endl;
336 
337  // Update the histogram titles according to the ion map
338  // and print new titles
339  G4AnalysisManager* analysisManager = G4AnalysisManager::Instance();
340  for ( const auto& ionMapElement : fgIonMap ) {
341  G4String ionName = ionMapElement.first;
342  G4int h1Id = ionMapElement.second;
343  // print new titles
344  G4cout << " " << std::setw(20) << ionName << " id = "<< std::setw(3) << h1Id
345  << G4endl;
346 
347  // update histogram ids
348  if ( ! analysisManager->GetH1(h1Id) ) continue;
349  // Skip inactive histograms, this is not necessary
350  // but it makes the code safe wrt modifications in future
351  G4String title = analysisManager->GetH1Title(h1Id);
352  title = ionName + title;
353  analysisManager->SetH1Title(h1Id, title);
354  }
355  G4cout << G4endl;
356 
357  //normalize histograms
358  G4int ih = 2;
359  G4double binWidth = analysisManager->GetH1Width(ih);
360  G4double fac = (1./(numberOfEvent*binWidth))*(mm/MeV);
361  analysisManager->ScaleH1(ih,fac);
362 
363  for (ih=14; ih<24; ih++) {
364  binWidth = analysisManager->GetH1Width(ih);
365  G4double unit = analysisManager->GetH1Unit(ih);
366  fac = (second/(binWidth*unit));
367  analysisManager->ScaleH1(ih,fac);
368  }
369 
370  //remove all contents in fProcCounter, fCount
371  fProcCounter.clear();
372  fParticleDataMap1.clear();
373  fParticleDataMap2.clear();
374  fgIonMap.clear();
375 
376  //restore default format
377  G4cout.precision(dfprec);
378 }
379 
380 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......