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 "G4ProcessTable.hh"
41 #include "G4UnitsTable.hh"
42 #include "G4SystemOfUnits.hh"
43 
44 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
45 
47 : G4Run(),
48  fDetector(det), fParticle(0), fEkin(0.),
49  fTotalCount(0), fGammaCount(0),
50  fSumTrack(0.), fSumTrack2(0.),
51  fTargetXXX(false)
52 {
53  for (G4int i=0; i<3; i++) { fPbalance[i] = 0. ; }
54  for (G4int i=0; i<3; i++) { fNbGamma[i] = 0 ; }
55 }
56 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
57 
58 Run::~Run()
59 { }
60 
61 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
62 
64 {
65  fParticle = particle;
66  fEkin = energy;
67 }
68 
69 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
70 
72 {
73  fTargetXXX = flag;
74 }
75 
76 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
77 
79 {
80  G4String procName = process->GetProcessName();
81  std::map<G4String,G4int>::iterator it = fProcCounter.find(procName);
82  if ( it == fProcCounter.end()) {
83  fProcCounter[procName] = 1;
84  }
85  else {
86  fProcCounter[procName]++;
87  }
88 }
89 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
90 
91 void Run::SumTrack(G4double trackl)
92 {
93  fTotalCount++;
94  fSumTrack += trackl; fSumTrack2 += trackl*trackl;
95 }
96 
97 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
98 
100 {
101  std::map<G4String, NuclChannel>::iterator it = fNuclChannelMap.find(name);
102  if ( it == fNuclChannelMap.end()) {
103  fNuclChannelMap[name] = NuclChannel(1, Q);
104  }
105  else {
106  NuclChannel& data = it->second;
107  data.fCount++;
108  data.fQ += Q;
109  }
110 }
111 
112 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
113 
115 {
116  std::map<G4String, ParticleData>::iterator it = fParticleDataMap.find(name);
117  if ( it == fParticleDataMap.end()) {
118  fParticleDataMap[name] = ParticleData(1, Ekin, Ekin, Ekin);
119  }
120  else {
121  ParticleData& data = it->second;
122  data.fCount++;
123  data.fEmean += Ekin;
124  //update min max
125  G4double emin = data.fEmin;
126  if (Ekin < emin) data.fEmin = Ekin;
127  G4double emax = data.fEmax;
128  if (Ekin > emax) data.fEmax = Ekin;
129  }
130 }
131 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
132 
134 {
135  fPbalance[0] += Pbal;
136  //update min max
137  if (fTotalCount == 1) fPbalance[1] = fPbalance[2] = Pbal;
138  if (Pbal < fPbalance[1]) fPbalance[1] = Pbal;
139  if (Pbal > fPbalance[2]) fPbalance[2] = Pbal;
140 }
141 
142 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
143 
144 void Run::CountGamma(G4int nGamma)
145 {
146  fGammaCount++;
147  fNbGamma[0] += nGamma;
148  //update min max
149  if (fGammaCount == 1) fNbGamma[1] = fNbGamma[2] = nGamma;
150  if (nGamma < fNbGamma[1]) fNbGamma[1] = nGamma;
151  if (nGamma > fNbGamma[2]) fNbGamma[2] = nGamma;
152 }
153 
154 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
155 
156 void Run::Merge(const G4Run* run)
157 {
158  const Run* localRun = static_cast<const Run*>(run);
159 
160  //primary particle info
161  //
162  fParticle = localRun->fParticle;
163  fEkin = localRun->fEkin;
164 
165  // accumulate sums
166  //
167  fTotalCount += localRun->fTotalCount;
168  fGammaCount += localRun->fGammaCount;
169  fSumTrack += localRun->fSumTrack;
170  fSumTrack2 += localRun->fSumTrack2;
171 
172 
173  //map: processes count
174  std::map<G4String,G4int>::const_iterator itp;
175  for ( itp = localRun->fProcCounter.begin();
176  itp != localRun->fProcCounter.end(); ++itp ) {
177 
178  G4String procName = itp->first;
179  G4int localCount = itp->second;
180  if ( fProcCounter.find(procName) == fProcCounter.end()) {
181  fProcCounter[procName] = localCount;
182  }
183  else {
184  fProcCounter[procName] += localCount;
185  }
186  }
187 
188  //map: nuclear channels
189  std::map<G4String,NuclChannel>::const_iterator itc;
190  for (itc = localRun->fNuclChannelMap.begin();
191  itc != localRun->fNuclChannelMap.end(); ++itc) {
192 
193  G4String name = itc->first;
194  const NuclChannel& localData = itc->second;
195  if ( fNuclChannelMap.find(name) == fNuclChannelMap.end()) {
196  fNuclChannelMap[name]
197  = NuclChannel(localData.fCount, localData.fQ);
198  }
199  else {
200  NuclChannel& data = fNuclChannelMap[name];
201  data.fCount += localData.fCount;
202  data.fQ += localData.fQ;
203  }
204  }
205 
206  //map: particles count
207  std::map<G4String,ParticleData>::const_iterator itn;
208  for (itn = localRun->fParticleDataMap.begin();
209  itn != localRun->fParticleDataMap.end(); ++itn) {
210 
211  G4String name = itn->first;
212  const ParticleData& localData = itn->second;
213  if ( fParticleDataMap.find(name) == fParticleDataMap.end()) {
214  fParticleDataMap[name]
215  = ParticleData(localData.fCount,
216  localData.fEmean,
217  localData.fEmin,
218  localData.fEmax);
219  }
220  else {
221  ParticleData& data = fParticleDataMap[name];
222  data.fCount += localData.fCount;
223  data.fEmean += localData.fEmean;
224  G4double emin = localData.fEmin;
225  if (emin < data.fEmin) data.fEmin = emin;
226  G4double emax = localData.fEmax;
227  if (emax > data.fEmax) data.fEmax = emax;
228  }
229  }
230 
231  G4Run::Merge(run);
232 }
233 
234 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
235 
237 {
238  G4int prec = 5, wid = prec + 2;
239  G4int dfprec = G4cout.precision(prec);
240 
241  //run condition
242  //
243  G4Material* material = fDetector->GetMaterial();
244  G4double density = material->GetDensity();
245 
246  G4String Particle = fParticle->GetParticleName();
247  G4cout << "\n The run is " << numberOfEvent << " "<< Particle << " of "
248  << G4BestUnit(fEkin,"Energy") << " through "
249  << G4BestUnit(fDetector->GetSize(),"Length") << " of "
250  << material->GetName() << " (density: "
251  << G4BestUnit(density,"Volumic Mass") << ")" << G4endl;
252 
253  if (numberOfEvent == 0) { G4cout.precision(dfprec); return;}
254 
255  //frequency of processes
256  //
257  G4cout << "\n Process calls frequency:" << G4endl;
258  G4int survive = 0;
259  std::map<G4String,G4int>::iterator it;
260  for (it = fProcCounter.begin(); it != fProcCounter.end(); it++) {
261  G4String procName = it->first;
262  G4int count = it->second;
263  G4cout << "\t" << procName << "= " << count;
264  if (procName == "Transportation") survive = count;
265  }
266  G4cout << G4endl;
267 
268  if (survive > 0) {
269  G4cout << "\n Nb of incident particles surviving after "
270  << G4BestUnit(fDetector->GetSize(),"Length") << " of "
271  << material->GetName() << " : " << survive << G4endl;
272  }
273 
274  if (fTotalCount == 0) fTotalCount = 1; //force printing anyway
275 
276  //compute mean free path and related quantities
277  //
278  G4double MeanFreePath = fSumTrack /fTotalCount;
279  G4double MeanTrack2 = fSumTrack2/fTotalCount;
280  G4double rms = std::sqrt(std::fabs(MeanTrack2 - MeanFreePath*MeanFreePath));
281  G4double CrossSection = 0.0;
282  if(MeanFreePath > 0.0) { CrossSection = 1./MeanFreePath; }
283  G4double massicMFP = MeanFreePath*density;
284  G4double massicCS = 0.0;
285  if(massicMFP > 0.0) { massicCS = 1./massicMFP; }
286 
287  G4cout << "\n\n MeanFreePath:\t" << G4BestUnit(MeanFreePath,"Length")
288  << " +- " << G4BestUnit( rms,"Length")
289  << "\tmassic: " << G4BestUnit(massicMFP, "Mass/Surface")
290  << "\n CrossSection:\t" << CrossSection*cm << " cm^-1 "
291  << "\t\tmassic: " << G4BestUnit(massicCS, "Surface/Mass")
292  << G4endl;
293 
294  //cross section per atom (only for single material)
295  //
296  if (material->GetNumberOfElements() == 1) {
297  G4double nbAtoms = material->GetTotNbOfAtomsPerVolume();
298  G4double crossSection = CrossSection/nbAtoms;
299  G4cout << " crossSection per atom:\t"
300  << G4BestUnit(crossSection,"Surface") << G4endl;
301  }
302  //check cross section from G4HadronicProcessStore
303  //
304  G4cout << "\n Verification: "
305  << "crossSections from G4HadronicProcessStore:";
306 
309  G4double sumc1 = 0.0, sumc2 = 0.0;
310  if (material->GetNumberOfElements() == 1) {
311  const G4Element* element = material->GetElement(0);
312  for (it = fProcCounter.begin(); it != fProcCounter.end(); it++) {
313  G4String procName = it->first;
314  G4VProcess* process = processTable->FindProcess(procName, fParticle);
315  G4double xs1 =
316  store->GetCrossSectionPerVolume(fParticle,fEkin,process,material);
317  G4double massSigma = xs1/density;
318  sumc1 += massSigma;
319  G4double xs2 =
320  store->GetCrossSectionPerAtom(fParticle,fEkin,process,element,material);
321  sumc2 += xs2;
322  G4cout << "\n" << std::setw(20) << procName << "= "
323  << G4BestUnit(massSigma, "Surface/Mass") << "\t"
324  << G4BestUnit(xs2, "Surface");
325 
326  }
327  G4cout << "\n" << std::setw(20) << "total" << "= "
328  << G4BestUnit(sumc1, "Surface/Mass") << "\t"
329  << G4BestUnit(sumc2, "Surface") << G4endl;
330  } else {
331  for (it = fProcCounter.begin(); it != fProcCounter.end(); it++) {
332  G4String procName = it->first;
333  G4VProcess* process = processTable->FindProcess(procName, fParticle);
334  G4double xs =
335  store->GetCrossSectionPerVolume(fParticle,fEkin,process,material);
336  G4double massSigma = xs/density;
337  sumc1 += massSigma;
338  G4cout << "\n" << std::setw(20) << procName << "= "
339  << G4BestUnit(massSigma, "Surface/Mass");
340  }
341  G4cout << "\n" << std::setw(20) << "total" << "= "
342  << G4BestUnit(sumc1, "Surface/Mass") << G4endl;
343  }
344 
345  //nuclear channel count
346  //
347  G4cout << "\n List of nuclear reactions: \n" << G4endl;
348  std::map<G4String,NuclChannel>::iterator ic;
349  for (ic = fNuclChannelMap.begin(); ic != fNuclChannelMap.end(); ic++) {
350  G4String name = ic->first;
351  NuclChannel data = ic->second;
352  G4int count = data.fCount;
353  G4double Q = data.fQ/count;
354  if (print)
355  G4cout << " " << std::setw(50) << name << ": " << std::setw(7) << count
356  << " Q = " << std::setw(wid) << G4BestUnit(Q, "Energy")
357  << G4endl;
358  }
359 
360  //Gamma count
361  //
362  if (print && (fGammaCount > 0)) {
363  G4cout << "\n" << std::setw(58) << "Number of gamma: N = "
364  << fNbGamma[1] << " --> " << fNbGamma[2] << G4endl;
365  }
366 
367  if (print && fTargetXXX) {
368  G4cout
369  << "\n --> NOTE: XXXX because neutronHP is unable to return target nucleus"
370  << G4endl;
371  }
372 
373  //particles count
374  //
375  G4cout << "\n List of generated particles:" << G4endl;
376 
377  std::map<G4String,ParticleData>::iterator itn;
378  for (itn = fParticleDataMap.begin(); itn != fParticleDataMap.end(); itn++) {
379  G4String name = itn->first;
380  ParticleData data = itn->second;
381  G4int count = data.fCount;
382  G4double eMean = data.fEmean/count;
383  G4double eMin = data.fEmin;
384  G4double eMax = data.fEmax;
385  if (print)
386  G4cout << " " << std::setw(13) << name << ": " << std::setw(7) << count
387  << " Emean = " << std::setw(wid) << G4BestUnit(eMean, "Energy")
388  << "\t( " << G4BestUnit(eMin, "Energy")
389  << " --> " << G4BestUnit(eMax, "Energy")
390  << ")" << G4endl;
391  }
392 
393  //energy momentum balance
394  //
395  if (fTotalCount > 1) {
396  G4double Pbmean = fPbalance[0]/fTotalCount;
397  G4cout << "\n Momentum balance: Pmean = "
398  << std::setw(wid) << G4BestUnit(Pbmean, "Energy")
399  << "\t( " << G4BestUnit(fPbalance[1], "Energy")
400  << " --> " << G4BestUnit(fPbalance[2], "Energy")
401  << ") \n" << G4endl;
402  }
403 
404  //normalize histograms
408 
409  //remove all contents in fProcCounter, fCount
410  fProcCounter.clear();
411  fNuclChannelMap.clear();
412  fParticleDataMap.clear();
413 
414  //restore default format
415  G4cout.precision(dfprec);
416 }
417 
418 //....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
void CountGamma(G4int)
Definition: Run.cc:144
void SetTargetXXX(G4bool)
Definition: Run.cc:71
static G4HadronicProcessStore * Instance()
Run()
Definition: Run.cc:43
G4int first(char) const
const G4String & GetName() const
Definition: G4Material.hh:178
G4double GetDensity() const
Definition: G4Material.hh:180
static double Q[]
#define G4BestUnit(a, b)
#define G4_USE_G4BESTUNIT_FOR_VERBOSE 1
void Balance(G4double)
Definition: Run.cc:133
const G4Element * GetElement(G4int iel) const
Definition: G4Material.hh:202
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 CountNuclearChannel(G4String, G4double)
Definition: Run.cc:99
static const double prec
Definition: RanecuEngine.cc:58
G4GLOB_DLL std::ostream G4cout
void SumTrack(G4double track)
Definition: Run.cc:81
bool G4bool
Definition: G4Types.hh:79
static constexpr double cm
Definition: G4SIunits.hh:119
Definition: G4Run.hh:46
const G4String & GetProcessName() const
Definition: G4VProcess.hh:408
void print(G4double elem)
static const G4double emax
G4double GetTotNbOfAtomsPerVolume() const
Definition: G4Material.hh:209
G4double energy(const ThreeVector &p, const G4double m)
G4double GetCrossSectionPerVolume(const G4ParticleDefinition *particle, G4double kineticEnergy, const G4VProcess *process, const G4Material *material)
Detector construction class to define materials and geometry.
void EndOfRun()
Definition: Run.cc:147
#define G4endl
Definition: G4ios.hh:61
size_t GetNumberOfElements() const
Definition: G4Material.hh:186
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
static G4ProcessTable * GetProcessTable()
G4VProcess * FindProcess(const G4String &processName, const G4String &particleName) const
G4double GetCrossSectionPerAtom(const G4ParticleDefinition *particle, G4double kineticEnergy, const G4VProcess *process, const G4Element *element, const G4Material *material=0)
~Run()
Definition: Run.cc:72
void ParticleCount(G4String, G4double)
Definition: Run.cc:114