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