Geant4  10.02.p02
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"
40 #include "G4RadioactiveDecay.hh"
41 #include "G4TwoVector.hh"
42 #include "G4UnitsTable.hh"
43 #include "G4SystemOfUnits.hh"
44 
45 #include <fstream>
46 
47 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
48 
50 : G4Run(),
51  fDetector(det), fParticle(0), fEkin(0.)
52 {
54  fEdepDetect = fEdepDetect2 = 0.;
55 }
56 
57 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
58 
59 Run::~Run()
60 { }
61 
62 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
63 
65 {
66  fParticle = particle;
67  fEkin = energy;
68 }
69 
70 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
71 
72 void Run::CountProcesses(const G4VProcess* process, G4int iVol)
73 {
74  G4String procName = process->GetProcessName();
75 
76  if (iVol == 1) {
77  std::map<G4String,G4int>::iterator it1 = fProcCounter1.find(procName);
78  if ( it1 == fProcCounter1.end()) {
79  fProcCounter1[procName] = 1;
80  }
81  else {
82  fProcCounter1[procName]++;
83  }
84  }
85 
86  if (iVol == 2) {
87  std::map<G4String,G4int>::iterator it2 = fProcCounter2.find(procName);
88  if ( it2 == fProcCounter2.end()) {
89  fProcCounter2[procName] = 1;
90  }
91  else {
92  fProcCounter2[procName]++;
93  }
94  }
95 }
96 
97 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
98 
100 {
101  if (iVol == 1) {
102  std::map<G4String, ParticleData>::iterator it = fParticleDataMap1.find(name);
103  if ( it == fParticleDataMap1.end()) {
104  fParticleDataMap1[name] = ParticleData(1, Ekin, Ekin, Ekin);
105  }
106  else {
107  ParticleData& data = it->second;
108  data.fCount++;
109  data.fEmean += Ekin;
110  //update min max
111  G4double emin = data.fEmin;
112  if (Ekin < emin) data.fEmin = Ekin;
113  G4double emax = data.fEmax;
114  if (Ekin > emax) data.fEmax = Ekin;
115  }
116  }
117 
118  if (iVol == 2) {
119  std::map<G4String, ParticleData>::iterator it = fParticleDataMap2.find(name);
120  if ( it == fParticleDataMap2.end()) {
121  fParticleDataMap2[name] = ParticleData(1, Ekin, Ekin, Ekin);
122  }
123  else {
124  ParticleData& data = it->second;
125  data.fCount++;
126  data.fEmean += Ekin;
127  //update min max
128  G4double emin = data.fEmin;
129  if (Ekin < emin) data.fEmin = Ekin;
130  G4double emax = data.fEmax;
131  if (Ekin > emax) data.fEmax = Ekin;
132  }
133  }
134 }
135 
136 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
137 
138 void Run::AddEdep(G4double edep1, G4double edep2)
139 {
140  fEdepTarget += edep1;
141  fEdepTarget2 += edep1*edep1;
142  fEdepDetect += edep2;
143  fEdepDetect2 += edep2*edep2;
144 }
145 
146 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
147 
148 void Run::Merge(const G4Run* run)
149 {
150  const Run* localRun = static_cast<const Run*>(run);
151 
152  //primary particle info
153  //
154  fParticle = localRun->fParticle;
155  fEkin = localRun->fEkin;
156 
157  // accumulate sums
158  //
159  fEdepTarget += localRun->fEdepTarget;
160  fEdepTarget2 += localRun->fEdepTarget2;
161  fEdepDetect += localRun->fEdepDetect;
162  fEdepDetect2 += localRun->fEdepDetect2;
163 
164  //map: processes count in target
165 
166  std::map<G4String,G4int>::const_iterator itp1;
167  for ( itp1 = localRun->fProcCounter1.begin();
168  itp1 != localRun->fProcCounter1.end(); ++itp1 ) {
169 
170  G4String procName = itp1->first;
171  G4int localCount = itp1->second;
172  if ( fProcCounter1.find(procName) == fProcCounter1.end()) {
173  fProcCounter1[procName] = localCount;
174  }
175  else {
176  fProcCounter1[procName] += localCount;
177  }
178  }
179 
180  //map: processes count in detector
181 
182  std::map<G4String,G4int>::const_iterator itp2;
183  for ( itp2 = localRun->fProcCounter2.begin();
184  itp2 != localRun->fProcCounter2.end(); ++itp2 ) {
185 
186  G4String procName = itp2->first;
187  G4int localCount = itp2->second;
188  if ( fProcCounter2.find(procName) == fProcCounter2.end()) {
189  fProcCounter2[procName] = localCount;
190  }
191  else {
192  fProcCounter2[procName] += localCount;
193  }
194  }
195 
196  //map: created particles in target
197  std::map<G4String,ParticleData>::const_iterator itc;
198  for (itc = localRun->fParticleDataMap1.begin();
199  itc != localRun->fParticleDataMap1.end(); ++itc) {
200 
201  G4String name = itc->first;
202  const ParticleData& localData = itc->second;
203  if ( fParticleDataMap1.find(name) == fParticleDataMap1.end()) {
205  = ParticleData(localData.fCount,
206  localData.fEmean,
207  localData.fEmin,
208  localData.fEmax);
209  }
210  else {
211  ParticleData& data = fParticleDataMap1[name];
212  data.fCount += localData.fCount;
213  data.fEmean += localData.fEmean;
214  G4double emin = localData.fEmin;
215  if (emin < data.fEmin) data.fEmin = emin;
216  G4double emax = localData.fEmax;
217  if (emax > data.fEmax) data.fEmax = emax;
218  }
219  }
220 
221  //map: created particle in detector
222  std::map<G4String,ParticleData>::const_iterator itn;
223  for (itn = localRun->fParticleDataMap2.begin();
224  itn != localRun->fParticleDataMap2.end(); ++itn) {
225 
226  G4String name = itn->first;
227  const ParticleData& localData = itn->second;
228  if ( fParticleDataMap2.find(name) == fParticleDataMap2.end()) {
230  = ParticleData(localData.fCount,
231  localData.fEmean,
232  localData.fEmin,
233  localData.fEmax);
234  }
235  else {
236  ParticleData& data = fParticleDataMap2[name];
237  data.fCount += localData.fCount;
238  data.fEmean += localData.fEmean;
239  G4double emin = localData.fEmin;
240  if (emin < data.fEmin) data.fEmin = emin;
241  G4double emax = localData.fEmax;
242  if (emax > data.fEmax) data.fEmax = emax;
243  }
244  }
245 
246  G4Run::Merge(run);
247 }
248 
249 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
250 
251 void Run::EndOfRun()
252 {
253  G4int prec = 5, wid = prec + 2;
254  G4int dfprec = G4cout.precision(prec);
255 
256  // run condition
257  //
258  G4String Particle = fParticle->GetParticleName();
259  G4cout << "\n The run is " << numberOfEvent << " "<< Particle << " of "
260  << G4BestUnit(fEkin,"Energy") << " through : ";
261 
262  G4cout << "\n Target : Length = "
263  << G4BestUnit(fDetector->GetTargetLength(),"Length")
264  << " Radius = "
265  << G4BestUnit(fDetector->GetTargetRadius(),"Length")
266  << " Material = "
268  G4cout << "\n Detector : Length = "
269  << G4BestUnit(fDetector->GetDetectorLength(),"Length")
270  << " Thickness = "
271  << G4BestUnit(fDetector->GetDetectorThickness(),"Length")
272  << " Material = "
274 
275  if (numberOfEvent == 0) { G4cout.precision(dfprec); return;}
276 
277  // compute mean Energy deposited and rms in target
278  //
279  G4int TotNbofEvents = numberOfEvent;
280  fEdepTarget /= TotNbofEvents; fEdepTarget2 /= TotNbofEvents;
282  if (rmsEdep>0.) rmsEdep = std::sqrt(rmsEdep);
283  else rmsEdep = 0.;
284 
285  G4cout << "\n Mean energy deposit in target, in time window = "
286  << G4BestUnit(fEdepTarget,"Energy") << "; rms = "
287  << G4BestUnit(rmsEdep, "Energy")
288  << G4endl;
289 
290  // compute mean Energy deposited and rms in detector
291  //
292  fEdepDetect /= TotNbofEvents; fEdepDetect2 /= TotNbofEvents;
294  if (rmsEdep>0.) rmsEdep = std::sqrt(rmsEdep);
295  else rmsEdep = 0.;
296 
297  G4cout << " Mean energy deposit in detector, in time window = "
298  << G4BestUnit(fEdepDetect,"Energy") << "; rms = "
299  << G4BestUnit(rmsEdep, "Energy")
300  << G4endl;
301 
302  // frequency of processes in target
303  //
304  G4cout << "\n Process calls frequency in target :" << G4endl;
305  G4int index = 0;
306  std::map<G4String,G4int>::iterator it1;
307  for (it1 = fProcCounter1.begin(); it1 != fProcCounter1.end(); it1++) {
308  G4String procName = it1->first;
309  G4int count = it1->second;
310  G4String space = " "; if (++index%3 == 0) space = "\n";
311  G4cout << " " << std::setw(20) << procName << "="<< std::setw(7) << count
312  << space;
313  }
314  G4cout << G4endl;
315 
316  // frequency of processes in detector
317  //
318  G4cout << "\n Process calls frequency in detector:" << G4endl;
319  index = 0;
320  std::map<G4String,G4int>::iterator it2;
321  for (it2 = fProcCounter2.begin(); it2 != fProcCounter2.end(); it2++) {
322  G4String procName = it2->first;
323  G4int count = it2->second;
324  G4String space = " "; if (++index%3 == 0) space = "\n";
325  G4cout << " " << std::setw(20) << procName << "="<< std::setw(7) << count
326  << space;
327  }
328  G4cout << G4endl;
329 
330  // particles count in target
331  //
332  G4cout << "\n List of generated particles in target:" << G4endl;
333 
334  std::map<G4String,ParticleData>::iterator itc;
335  for (itc = fParticleDataMap1.begin(); itc != fParticleDataMap1.end(); itc++) {
336  G4String name = itc->first;
337  ParticleData data = itc->second;
338  G4int count = data.fCount;
339  G4double eMean = data.fEmean/count;
340  G4double eMin = data.fEmin;
341  G4double eMax = data.fEmax;
342 
343  G4cout << " " << std::setw(13) << name << ": " << std::setw(7) << count
344  << " Emean = " << std::setw(wid) << G4BestUnit(eMean, "Energy")
345  << "\t( " << G4BestUnit(eMin, "Energy")
346  << " --> " << G4BestUnit(eMax, "Energy")
347  << ")" << G4endl;
348  }
349 
350  // particles count in detector
351  //
352  G4cout << "\n List of generated particles in detector:" << G4endl;
353 
354  std::map<G4String,ParticleData>::iterator itn;
355  for (itn = fParticleDataMap2.begin(); itn != fParticleDataMap2.end(); itn++) {
356  G4String name = itn->first;
357  ParticleData data = itn->second;
358  G4int count = data.fCount;
359  G4double eMean = data.fEmean/count;
360  G4double eMin = data.fEmin;
361  G4double eMax = data.fEmax;
362 
363  G4cout << " " << std::setw(13) << name << ": " << std::setw(7) << count
364  << " Emean = " << std::setw(wid) << G4BestUnit(eMean, "Energy")
365  << "\t( " << G4BestUnit(eMin, "Energy")
366  << " --> " << G4BestUnit(eMax, "Energy") << ")" << G4endl;
367  }
368  G4cout << G4endl;
369 
370  // activities in VR mode
371  //
373 
374  //remove all contents in fProcCounter, fCount
375  fProcCounter1.clear();
376  fProcCounter2.clear();
377  fParticleDataMap1.clear();
378  fParticleDataMap2.clear();
379 
380  //restore default format
381  G4cout.precision(dfprec);
382 }
383 
384 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
385 
387 {
390  pTable->FindProcess("RadioactiveDecay", "GenericIon");
391 
392  // output the induced radioactivities (in VR mode only)
393  //
394  if ((rDecay == 0) || (rDecay->IsAnalogueMonteCarlo())) return;
395 
396  G4String fileName = G4AnalysisManager::Instance()->GetFileName() + ".activity";
397  std::ofstream outfile (fileName, std::ios::out );
398 
399  std::vector<G4RadioactivityTable*> theTables =
400  rDecay->GetTheRadioactivityTables();
401 
402  for (size_t i = 0 ; i < theTables.size(); i++) {
403  G4double rate, error;
404  outfile << "Radioactivities in decay window no. " << i << G4endl;
405  outfile << "Z \tA \tE \tActivity (decays/window) \tError (decays/window) "
406  << G4endl;
407 
408  map<G4ThreeVector,G4TwoVector> *aMap = theTables[i]->GetTheMap();
409  map<G4ThreeVector,G4TwoVector>::iterator iter;
410  for (iter=aMap->begin(); iter != aMap->end(); iter++) {
411  rate = iter->second.x()/nevent;
412  error = std::sqrt(iter->second.y())/nevent;
413  if (rate < 0.) rate = 0.; // statically it can be < 0.
414  outfile << iter->first.x() <<"\t"<< iter->first.y() <<"\t"
415  << iter->first.z() << "\t" << rate <<"\t" << error << G4endl;
416  }
417  outfile << G4endl;
418  }
419  outfile.close();
420 }
421 
422 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
virtual void Merge(const G4Run *)
Definition: G4Run.cc:54
G4int numberOfEvent
Definition: G4Run.hh:59
G4double fEkin
Definition: Run.hh:65
void CountProcesses(G4String procName)
Definition: Run.cc:72
Run()
Definition: Run.cc:43
G4int first(char) const
G4String name
Definition: TRTMaterials.hh:40
const G4String & GetName() const
Definition: G4Material.hh:178
G4double fEmean
Definition: Run.hh:73
G4double fEdepDetect2
Definition: Run.hh:81
void WriteActivity(G4int)
Definition: Run.cc:386
#define G4BestUnit(a, b)
#define G4_USE_G4BESTUNIT_FOR_VERBOSE 1
G4double fEdepTarget2
Definition: Run.hh:80
int G4int
Definition: G4Types.hh:78
const G4String & GetParticleName() const
void AddEdep(G4double val)
Definition: Run.hh:61
static const double prec
Definition: RanecuEngine.cc:58
G4GLOB_DLL std::ostream G4cout
std::map< G4String, G4int > fProcCounter2
Definition: Run.hh:84
G4Material * GetDetectorMaterial()
std::vector< G4RadioactivityTable * > GetTheRadioactivityTables()
std::map< G4String, G4int > fProcCounter1
Definition: Run.hh:83
G4double fEmax
Definition: Run.hh:75
G4int fCount
Definition: Run.hh:72
Definition: G4Run.hh:46
G4double fEmin
Definition: Run.hh:74
G4double fEdepDetect
Definition: Run.hh:81
const G4String & GetProcessName() const
Definition: G4VProcess.hh:408
static const G4double emax
G4double fEdepTarget
Definition: Run.hh:80
G4double energy(const ThreeVector &p, const G4double m)
std::map< G4String, ParticleData > fParticleDataMap1
Definition: Run.hh:84
std::map< G4String, ParticleData > fParticleDataMap2
Definition: Run.hh:85
DetectorConstruction * fDetector
Definition: Run.hh:63
Detector construction class to demonstrate various ways of placement.
void EndOfRun()
Definition: Run.cc:147
#define G4endl
Definition: G4ios.hh:61
G4ParticleDefinition * fParticle
Definition: Run.hh:64
static PROLOG_HANDLER error
Definition: xmlrole.cc:112
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
~Run()
Definition: Run.cc:72
const G4Material * GetTargetMaterial() const
void ParticleCount(G4String, G4double)
Definition: Run.cc:114