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 86064 2014-11-07 08:49:32Z gcosmo $
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 "G4Run.hh"
40 #include "G4RunManager.hh"
41 #include "G4UnitsTable.hh"
42 
43 #include "G4PhysicalConstants.hh"
44 #include "G4SystemOfUnits.hh"
45 #include "Randomize.hh"
46 #include <iomanip>
47 #include <stdexcept> // std::out_of_range
48 
49 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
50 
52 : G4Run(),
53  fDetector(det),fParticle(0),
54  fEnergy(0)
55 {
56  InitFluence();
57 }
58 
59 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
60 
61 Run::~Run()
62 {
63 }
64 
65 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
66 
68 {
69  fParticle = particle;
70  fEnergy = energy;
71 }
72 
73 
74 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
75 
76 void Run::Merge(const G4Run* run)
77 {
78  const Run* localRun = static_cast<const Run*>(run);
79 
80  // pass information about primary particle
81  fParticle = localRun->fParticle;
82  fEnergy = localRun->fEnergy;
83 
84  // In MT all threads have the same fNbBins and fDr value
85  fNbBins = localRun->fNbBins ;
86  fDr = localRun->fDr;
87 
88  // Some value by value
89 
90  for (unsigned int i = 0; i < localRun->fluence.size(); i++)
91  {
92  try { fluence[i]+=localRun->fluence[i]; }
93  catch(const std::out_of_range& oor)
94  { fluence.push_back(localRun->fluence[i]); }
95  }
96 
97  for (unsigned int i = 0; i < localRun->fluence1.size(); i++)
98  {
99  try { fluence1[i]+=localRun->fluence1[i]; }
100  catch(const std::out_of_range& oor)
101  { fluence1.push_back(localRun->fluence1[i]); }
102  }
103 
104  for (unsigned int i = 0; i < localRun->fluence2.size(); i++)
105  {
106  try { fluence2[i]+=localRun->fluence2[i]; }
107  catch(const std::out_of_range& oor)
108  { fluence2.push_back(localRun->fluence2[i]); }
109  }
110 
111  for (unsigned int i = 0; i < localRun->fNbEntries.size(); i++)
112  {
113  try { fNbEntries[i]+=localRun->fNbEntries[i]; }
114  catch(const std::out_of_range& oor)
115  { fNbEntries.push_back(localRun->fNbEntries[i]); }
116  }
117 
118  G4Run::Merge(run);
119 }
120 
121 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
122 
123 
124 void Run::EndOfRun()
125 {
126 
127 
128  if (numberOfEvent == 0) return;
129 
130  //Scatter foil
131 
132  G4Material* material = fDetector->GetMaterialScatter();
133  G4double length = fDetector->GetThicknessScatter();
134  G4double density = material->GetDensity();
135  G4String partName = fParticle->GetParticleName();
136 
137  G4cout << "\n ======================== run summary ======================\n";
138 
139  G4int prec = G4cout.precision(3);
140 
141  G4cout << "\n The run was " << numberOfEvent << " " << partName << " of "
142  << G4BestUnit(fEnergy,"Energy") << " through "
143  << G4BestUnit(length,"Length") << " of "
144  << material->GetName() << " (density: "
145  << G4BestUnit(density,"Volumic Mass") << ")" << G4endl;
146 
147  G4cout.precision(prec);
148 
149  // normalize histograms
150  //
151  G4AnalysisManager* analysisManager = G4AnalysisManager::Instance();
152  G4double fac = 1./(double(numberOfEvent));
153  analysisManager->ScaleH1(1,fac);
154  analysisManager->ScaleH1(2,fac);
155  analysisManager->ScaleH1(3,fac);
156  analysisManager->ScaleH1(5,fac);
157  analysisManager->ScaleH1(6,fac);
158 
159  ComputeFluenceError();
160  PrintFluence(numberOfEvent);
161 
162 }
163 
164 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
165 
166 void Run::InitFluence()
167 {
168  // create Fluence histo in any case
169 
170  G4AnalysisManager* analysisManager = G4AnalysisManager::Instance();
171  G4int ih = 4;
172  analysisManager->SetH1(ih, 120, 0*mm, 240*mm, "mm");
173 
174  //construct vectors for fluence distribution
175 
176  fNbBins = analysisManager->GetH1Nbins(ih);
177  fDr = analysisManager->GetH1Width(ih);
178  fluence.resize(fNbBins, 0.);
179  fluence1.resize(fNbBins, 0.);
180  fluence2.resize(fNbBins, 0.);
181  fNbEntries.resize(fNbBins, 0);
182 }
183 
184 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
185 
187 {
188  G4int ibin = (int)(r/fDr);
189  if (ibin >= fNbBins) return;
190  fNbEntries[ibin]++;
191  fluence[ibin] += fl;
192  fluence2[ibin] += fl*fl;
193 }
194 
195 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
196 
197 void Run::ComputeFluenceError()
198 {
199  //compute rms
200 
201  G4double ds,variance,rms;
202  G4double rmean = -0.5*fDr;
203 
204  for (G4int bin=0; bin<fNbBins; bin++) {
205  rmean += fDr;
206  ds = twopi*rmean*fDr;
207  fluence[bin] /= ds;
208  fluence2[bin] /= (ds*ds);
209  variance = 0.;
210  if (fNbEntries[bin] > 0)
211  variance = fluence2[bin] - (fluence[bin]*fluence[bin])/fNbEntries[bin];
212  rms = 0.;
213  if(variance > 0.) rms = std::sqrt(variance);
214  fluence2[bin] = rms;
215  }
216 
217  //normalize to first bins, compute error and fill histo
218 
219  G4double rnorm(4*mm), radius(0.), fnorm(0.), fnorm2(0.);
220  G4int inorm = -1;
221  do {
222  inorm++; radius += fDr; fnorm += fluence[inorm]; fnorm2 += fluence2[inorm];
223  } while (radius < rnorm);
224  fnorm /= (inorm+1);
225  fnorm2 /= (inorm+1);
226 
227  G4double ratio, error;
228  G4double scale = 1./fnorm;
229  G4double err0 = fnorm2/fnorm, err1 = 0.;
230 
231  rmean = -0.5*fDr;
232 
233  for (G4int bin=0; bin<fNbBins; bin++) {
234  ratio = fluence[bin]*scale;
235  error = 0.;
236  if (ratio > 0.) {
237  err1 = fluence2[bin]/fluence[bin];
238  error = ratio*std::sqrt(err1*err1 + err0*err0);
239  }
240  fluence1[bin] = ratio;
241  fluence2[bin] = error;
242  rmean += fDr;
243  G4AnalysisManager* analysisManager = G4AnalysisManager::Instance();
244  analysisManager->FillH1(4,rmean,ratio);
245  }
246 }
247 
248 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
249 
250 #include <fstream>
251 
252 void Run::PrintFluence(G4int TotEvents)
253 {
254  G4AnalysisManager* analysisManager = G4AnalysisManager::Instance();
255  G4String name = analysisManager->GetFileName();
256  G4String fileName = name + ".ascii";
257  std::ofstream File(fileName, std::ios::out);
258 
259  std::ios::fmtflags mode = File.flags();
260  File.setf( std::ios::scientific, std::ios::floatfield );
261  G4int prec = File.precision(3);
262 
263  File << " Fluence density distribution \n "
264  << "\n ibin \t radius (mm) \t Nb \t fluence\t norma fl\t rms/nfl (%) \n"
265  << G4endl;
266 
267  G4double rmean = -0.5*fDr;
268  for (G4int bin=0; bin<fNbBins; bin++) {
269  rmean +=fDr;
270  G4double error = 0.;
271  if (fluence1[bin] > 0.) error = 100*fluence2[bin]/fluence1[bin];
272  File << " " << bin << "\t " << rmean/mm << "\t " << fNbEntries[bin]
273  << "\t " << fluence[bin]/double(TotEvents) << "\t " << fluence1[bin]
274  << "\t " << error
275  << G4endl;
276  }
277 
278  // restaure default formats
279  File.setf(mode,std::ios::floatfield);
280  File.precision(prec);
281 }
282 
283 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
284 
const XML_Char * name
Definition: expat.h:151
virtual void Merge(const G4Run *)
Definition: G4Run.cc:54
G4int numberOfEvent
Definition: G4Run.hh:59
static constexpr double mm
Definition: G4SIunits.hh:115
tuple bin
Definition: plottest35.py:22
G4Material * GetMaterialScatter()
Run()
Definition: Run.cc:43
const G4String & GetName() const
Definition: G4Material.hh:178
G4double GetDensity() const
Definition: G4Material.hh:180
#define G4BestUnit(a, b)
#define G4_USE_G4BESTUNIT_FOR_VERBOSE 1
int G4int
Definition: G4Types.hh:78
const G4String & GetParticleName() const
static constexpr double twopi
Definition: G4SIunits.hh:76
string material
Definition: eplot.py:19
static const double prec
Definition: RanecuEngine.cc:58
G4GLOB_DLL std::ostream G4cout
typedef int(XMLCALL *XML_NotStandaloneHandler)(void *userData)
Definition: G4Run.hh:46
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 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
void SumFluence(G4double, G4double)
Definition: Run.cc:186
virtual void Merge(const G4Run *)
Definition: Run.cc:115
G4CsvAnalysisManager G4AnalysisManager
Definition: g4csv_defs.hh:77
~Run()
Definition: Run.cc:72