| Geant4 Cross Reference |
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 //
26 /// \file medical/electronScattering/src/Run.cc
27 /// \brief Implementation of the Run class
28 //
29 //
30 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
31 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
32
33 #include "Run.hh"
34
35 #include "DetectorConstruction.hh"
36 #include "HistoManager.hh"
37 #include "PrimaryGeneratorAction.hh"
38
39 #include "G4PhysicalConstants.hh"
40 #include "G4Run.hh"
41 #include "G4RunManager.hh"
42 #include "G4SystemOfUnits.hh"
43 #include "G4UnitsTable.hh"
44 #include "Randomize.hh"
45
46 #include <iomanip>
47 #include <stdexcept> // std::out_of_range
48
49 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
50
51 Run::Run(DetectorConstruction* det) : G4Run(), fDetector(det), fParticle(0), fEnergy(0)
52 {
53 InitFluence();
54 }
55
56 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
57
58 Run::~Run() {}
59
60 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
61
62 void Run::SetPrimary(G4ParticleDefinition* particle, G4double energy)
63 {
64 fParticle = particle;
65 fEnergy = energy;
66 }
67
68 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
69
70 void Run::Merge(const G4Run* run)
71 {
72 const Run* localRun = static_cast<const Run*>(run);
73
74 // pass information about primary particle
75 fParticle = localRun->fParticle;
76 fEnergy = localRun->fEnergy;
77
78 // In MT all threads have the same fNbBins and fDr value
79 fNbBins = localRun->fNbBins;
80 fDr = localRun->fDr;
81
82 // Some value by value
83
84 for (unsigned int i = 0; i < localRun->fluence.size(); ++i) {
85 try {
86 fluence[i] += localRun->fluence[i];
87 }
88 catch (const std::out_of_range&) {
89 fluence.push_back(localRun->fluence[i]);
90 }
91 }
92
93 for (unsigned int i = 0; i < localRun->fluence1.size(); ++i) {
94 try {
95 fluence1[i] += localRun->fluence1[i];
96 }
97 catch (const std::out_of_range&) {
98 fluence1.push_back(localRun->fluence1[i]);
99 }
100 }
101
102 for (unsigned int i = 0; i < localRun->fluence2.size(); ++i) {
103 try {
104 fluence2[i] += localRun->fluence2[i];
105 }
106 catch (const std::out_of_range&) {
107 fluence2.push_back(localRun->fluence2[i]);
108 }
109 }
110
111 for (unsigned int i = 0; i < localRun->fNbEntries.size(); ++i) {
112 try {
113 fNbEntries[i] += localRun->fNbEntries[i];
114 }
115 catch (const std::out_of_range&) {
116 fNbEntries.push_back(localRun->fNbEntries[i]);
117 }
118 }
119
120 G4Run::Merge(run);
121 }
122
123 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
124
125 void Run::EndOfRun()
126 {
127 if (numberOfEvent == 0) return;
128
129 // Scatter foil
130
131 G4Material* material = fDetector->GetMaterialScatter();
132 G4double length = fDetector->GetThicknessScatter();
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 " << G4BestUnit(length, "Length") << " of "
142 << material->GetName() << " (density: " << G4BestUnit(density, "Volumic Mass") << ")"
143 << G4endl;
144
145 G4cout.precision(prec);
146
147 // normalize histograms
148 //
149 G4AnalysisManager* analysisManager = G4AnalysisManager::Instance();
150 G4double fac = 1. / (double(numberOfEvent));
151 analysisManager->ScaleH1(1, fac);
152 analysisManager->ScaleH1(2, fac);
153 analysisManager->ScaleH1(3, fac);
154 analysisManager->ScaleH1(5, fac);
155 analysisManager->ScaleH1(6, fac);
156
157 ComputeFluenceError();
158 PrintFluence(numberOfEvent);
159 }
160
161 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
162
163 void Run::InitFluence()
164 {
165 // create Fluence histo in any case
166
167 G4AnalysisManager* analysisManager = G4AnalysisManager::Instance();
168 G4int ih = 4;
169 analysisManager->SetH1(ih, 120, 0 * mm, 240 * mm, "mm");
170
171 // construct vectors for fluence distribution
172
173 fNbBins = analysisManager->GetH1Nbins(ih);
174 fDr = analysisManager->GetH1Width(ih);
175 fluence.resize(fNbBins, 0.);
176 fluence1.resize(fNbBins, 0.);
177 fluence2.resize(fNbBins, 0.);
178 fNbEntries.resize(fNbBins, 0);
179 }
180
181 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
182
183 void Run::SumFluence(G4double r, G4double fl)
184 {
185 G4int ibin = (int)(r / fDr);
186 if (ibin >= fNbBins) return;
187 fNbEntries[ibin]++;
188 fluence[ibin] += fl;
189 fluence2[ibin] += fl * fl;
190 }
191
192 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
193
194 void Run::ComputeFluenceError()
195 {
196 // compute rms
197
198 G4double ds, variance, rms;
199 G4double rmean = -0.5 * fDr;
200
201 for (G4int bin = 0; bin < fNbBins; bin++) {
202 rmean += fDr;
203 ds = twopi * rmean * fDr;
204 fluence[bin] /= ds;
205 fluence2[bin] /= (ds * ds);
206 variance = 0.;
207 if (fNbEntries[bin] > 0)
208 variance = fluence2[bin] - (fluence[bin] * fluence[bin]) / fNbEntries[bin];
209 rms = 0.;
210 if (variance > 0.) rms = std::sqrt(variance);
211 fluence2[bin] = rms;
212 }
213
214 // normalize to first bins, compute error and fill histo
215
216 G4double rnorm(4 * mm), radius(0.), fnorm(0.), fnorm2(0.);
217 G4int inorm = -1;
218 do {
219 inorm++;
220 radius += fDr;
221 fnorm += fluence[inorm];
222 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] << "\t "
273 << fluence[bin] / double(TotEvents) << "\t " << fluence1[bin] << "\t " << error << G4endl;
274 }
275
276 // restaure default formats
277 File.setf(mode, std::ios::floatfield);
278 File.precision(prec);
279 }
280
281 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
282