| Geant4 Cross Reference |
1 // 1
2 // *******************************************
3 // * License and Disclaimer
4 // *
5 // * The Geant4 software is copyright of th
6 // * the Geant4 Collaboration. It is provided
7 // * conditions of the Geant4 Software License
8 // * LICENSE and available at http://cern.ch/
9 // * include a list of copyright holders.
10 // *
11 // * Neither the authors of this software syst
12 // * institutes,nor the agencies providing fin
13 // * work make any representation or warran
14 // * regarding this software system or assum
15 // * use. Please see the license in the file
16 // * for the full disclaimer and the limitatio
17 // *
18 // * This code implementation is the result
19 // * technical work of the GEANT4 collaboratio
20 // * By using, copying, modifying or distri
21 // * any work based on the software) you ag
22 // * use in resulting scientific publicati
23 // * acceptance of all terms of the Geant4 Sof
24 // *******************************************
25 //
26 /// \file medical/electronScattering/src/Run.c
27 /// \brief Implementation of the Run class
28 //
29 //
30 //....oooOO0OOooo........oooOO0OOooo........oo
31 //....oooOO0OOooo........oooOO0OOooo........oo
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........oo
50
51 Run::Run(DetectorConstruction* det) : G4Run(),
52 {
53 InitFluence();
54 }
55
56 //....oooOO0OOooo........oooOO0OOooo........oo
57
58 Run::~Run() {}
59
60 //....oooOO0OOooo........oooOO0OOooo........oo
61
62 void Run::SetPrimary(G4ParticleDefinition* par
63 {
64 fParticle = particle;
65 fEnergy = energy;
66 }
67
68 //....oooOO0OOooo........oooOO0OOooo........oo
69
70 void Run::Merge(const G4Run* run)
71 {
72 const Run* localRun = static_cast<const 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 a
79 fNbBins = localRun->fNbBins;
80 fDr = localRun->fDr;
81
82 // Some value by value
83
84 for (unsigned int i = 0; i < localRun->fluen
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->fluen
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->fluen
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->fNbEn
112 try {
113 fNbEntries[i] += localRun->fNbEntries[i]
114 }
115 catch (const std::out_of_range&) {
116 fNbEntries.push_back(localRun->fNbEntrie
117 }
118 }
119
120 G4Run::Merge(run);
121 }
122
123 //....oooOO0OOooo........oooOO0OOooo........oo
124
125 void Run::EndOfRun()
126 {
127 if (numberOfEvent == 0) return;
128
129 // Scatter foil
130
131 G4Material* material = fDetector->GetMateria
132 G4double length = fDetector->GetThicknessSca
133 G4double density = material->GetDensity();
134 G4String partName = fParticle->GetParticleNa
135
136 G4cout << "\n ======================== run s
137
138 G4int prec = G4cout.precision(3);
139
140 G4cout << "\n The run was " << numberOfEvent
141 << G4BestUnit(fEnergy, "Energy") << "
142 << material->GetName() << " (density:
143 << G4endl;
144
145 G4cout.precision(prec);
146
147 // normalize histograms
148 //
149 G4AnalysisManager* analysisManager = G4Analy
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........oo
162
163 void Run::InitFluence()
164 {
165 // create Fluence histo in any case
166
167 G4AnalysisManager* analysisManager = G4Analy
168 G4int ih = 4;
169 analysisManager->SetH1(ih, 120, 0 * mm, 240
170
171 // construct vectors for fluence distributio
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........oo
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........oo
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]
209 rms = 0.;
210 if (variance > 0.) rms = std::sqrt(varianc
211 fluence2[bin] = rms;
212 }
213
214 // normalize to first bins, compute error an
215
216 G4double rnorm(4 * mm), radius(0.), fnorm(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 +
239 }
240 fluence1[bin] = ratio;
241 fluence2[bin] = error;
242 rmean += fDr;
243 G4AnalysisManager* analysisManager = G4Ana
244 analysisManager->FillH1(4, rmean, ratio);
245 }
246 }
247
248 //....oooOO0OOooo........oooOO0OOooo........oo
249
250 #include <fstream>
251
252 void Run::PrintFluence(G4int TotEvents)
253 {
254 G4AnalysisManager* analysisManager = G4Analy
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::fl
261 G4int prec = File.precision(3);
262
263 File << " Fluence density distribution \n "
264 << "\n ibin \t radius (mm) \t Nb \t fl
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 * flue
272 File << " " << bin << "\t " << rmean / mm
273 << fluence[bin] / double(TotEvents) <
274 }
275
276 // restaure default formats
277 File.setf(mode, std::ios::floatfield);
278 File.precision(prec);
279 }
280
281 //....oooOO0OOooo........oooOO0OOooo........oo
282