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