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