| Geant4 Cross Reference |
1 // 1 //
2 // ******************************************* 2 // ********************************************************************
3 // * License and Disclaimer 3 // * License and Disclaimer *
4 // * 4 // * *
5 // * The Geant4 software is copyright of th 5 // * The Geant4 software is copyright of the Copyright Holders of *
6 // * the Geant4 Collaboration. It is provided 6 // * the Geant4 Collaboration. It is provided under the terms and *
7 // * conditions of the Geant4 Software License 7 // * conditions of the Geant4 Software License, included in the file *
8 // * LICENSE and available at http://cern.ch/ 8 // * LICENSE and available at http://cern.ch/geant4/license . These *
9 // * include a list of copyright holders. 9 // * include a list of copyright holders. *
10 // * 10 // * *
11 // * Neither the authors of this software syst 11 // * Neither the authors of this software system, nor their employing *
12 // * institutes,nor the agencies providing fin 12 // * institutes,nor the agencies providing financial support for this *
13 // * work make any representation or warran 13 // * work make any representation or warranty, express or implied, *
14 // * regarding this software system or assum 14 // * regarding this software system or assume any liability for its *
15 // * use. Please see the license in the file 15 // * use. Please see the license in the file LICENSE and URL above *
16 // * for the full disclaimer and the limitatio 16 // * for the full disclaimer and the limitation of liability. *
17 // * 17 // * *
18 // * This code implementation is the result 18 // * This code implementation is the result of the scientific and *
19 // * technical work of the GEANT4 collaboratio 19 // * technical work of the GEANT4 collaboration. *
20 // * By using, copying, modifying or distri 20 // * By using, copying, modifying or distributing the software (or *
21 // * any work based on the software) you ag 21 // * any work based on the software) you agree to acknowledge its *
22 // * use in resulting scientific publicati 22 // * use in resulting scientific publications, and indicate your *
23 // * acceptance of all terms of the Geant4 Sof 23 // * acceptance of all terms of the Geant4 Software license. *
24 // ******************************************* 24 // ********************************************************************
25 // 25 //
26 /// \file medical/electronScattering/src/Run.c 26 /// \file medical/electronScattering/src/Run.cc
27 /// \brief Implementation of the Run class 27 /// \brief Implementation of the Run class
28 // 28 //
29 // 29 //
30 //....oooOO0OOooo........oooOO0OOooo........oo 30 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
31 //....oooOO0OOooo........oooOO0OOooo........oo 31 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
32 32
33 #include "Run.hh" 33 #include "Run.hh"
34 <<
35 #include "DetectorConstruction.hh" 34 #include "DetectorConstruction.hh"
36 #include "HistoManager.hh" <<
37 #include "PrimaryGeneratorAction.hh" 35 #include "PrimaryGeneratorAction.hh"
>> 36 #include "HistoManager.hh"
38 37
39 #include "G4PhysicalConstants.hh" <<
40 #include "G4Run.hh" 38 #include "G4Run.hh"
41 #include "G4RunManager.hh" 39 #include "G4RunManager.hh"
42 #include "G4SystemOfUnits.hh" <<
43 #include "G4UnitsTable.hh" 40 #include "G4UnitsTable.hh"
44 #include "Randomize.hh" <<
45 41
>> 42 #include "G4PhysicalConstants.hh"
>> 43 #include "G4SystemOfUnits.hh"
>> 44 #include "Randomize.hh"
46 #include <iomanip> 45 #include <iomanip>
47 #include <stdexcept> // std::out_of_range << 46 #include <stdexcept> // std::out_of_range
48 47
49 //....oooOO0OOooo........oooOO0OOooo........oo 48 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
50 49
51 Run::Run(DetectorConstruction* det) : G4Run(), << 50 Run::Run(DetectorConstruction* det)
>> 51 : G4Run(),
>> 52 fDetector(det),fParticle(0),
>> 53 fEnergy(0)
52 { 54 {
53 InitFluence(); 55 InitFluence();
54 } 56 }
55 57
56 //....oooOO0OOooo........oooOO0OOooo........oo 58 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
57 59
58 Run::~Run() {} << 60 Run::~Run()
>> 61 {
>> 62 }
59 63
60 //....oooOO0OOooo........oooOO0OOooo........oo 64 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
61 65
62 void Run::SetPrimary(G4ParticleDefinition* par << 66 void Run::SetPrimary(G4ParticleDefinition* particle, G4double energy )
63 { 67 {
64 fParticle = particle; 68 fParticle = particle;
65 fEnergy = energy; 69 fEnergy = energy;
66 } 70 }
67 71
>> 72
68 //....oooOO0OOooo........oooOO0OOooo........oo 73 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
69 74
70 void Run::Merge(const G4Run* run) 75 void Run::Merge(const G4Run* run)
71 { 76 {
72 const Run* localRun = static_cast<const Run* << 77 const Run* localRun = static_cast<const Run*>(run);
73 78
74 // pass information about primary particle << 79 // pass information about primary particle
75 fParticle = localRun->fParticle; << 80 fParticle = localRun->fParticle;
76 fEnergy = localRun->fEnergy; << 81 fEnergy = localRun->fEnergy;
77 << 82
78 // In MT all threads have the same fNbBins a << 83 // In MT all threads have the same fNbBins and fDr value
79 fNbBins = localRun->fNbBins; << 84 fNbBins = localRun->fNbBins ;
80 fDr = localRun->fDr; << 85 fDr = localRun->fDr;
81 << 86
82 // Some value by value << 87 // Some value by value
83 << 88
84 for (unsigned int i = 0; i < localRun->fluen << 89 for (unsigned int i = 0; i < localRun->fluence.size(); ++i)
85 try { << 90 {
86 fluence[i] += localRun->fluence[i]; << 91 try { fluence[i]+=localRun->fluence[i]; }
87 } << 92 catch(const std::out_of_range&)
88 catch (const std::out_of_range&) { << 93 { fluence.push_back(localRun->fluence[i]); }
89 fluence.push_back(localRun->fluence[i]); << 94 }
90 } << 95
91 } << 96 for (unsigned int i = 0; i < localRun->fluence1.size(); ++i)
92 << 97 {
93 for (unsigned int i = 0; i < localRun->fluen << 98 try { fluence1[i]+=localRun->fluence1[i]; }
94 try { << 99 catch(const std::out_of_range&)
95 fluence1[i] += localRun->fluence1[i]; << 100 { fluence1.push_back(localRun->fluence1[i]); }
96 } << 101 }
97 catch (const std::out_of_range&) { << 102
98 fluence1.push_back(localRun->fluence1[i] << 103 for (unsigned int i = 0; i < localRun->fluence2.size(); ++i)
99 } << 104 {
100 } << 105 try { fluence2[i]+=localRun->fluence2[i]; }
101 << 106 catch(const std::out_of_range&)
102 for (unsigned int i = 0; i < localRun->fluen << 107 { fluence2.push_back(localRun->fluence2[i]); }
103 try { << 108 }
104 fluence2[i] += localRun->fluence2[i]; << 109
105 } << 110 for (unsigned int i = 0; i < localRun->fNbEntries.size(); ++i)
106 catch (const std::out_of_range&) { << 111 {
107 fluence2.push_back(localRun->fluence2[i] << 112 try { fNbEntries[i]+=localRun->fNbEntries[i]; }
108 } << 113 catch(const std::out_of_range&)
109 } << 114 { fNbEntries.push_back(localRun->fNbEntries[i]); }
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 } 115 }
119 116
120 G4Run::Merge(run); 117 G4Run::Merge(run);
121 } 118 }
122 119
123 //....oooOO0OOooo........oooOO0OOooo........oo 120 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
124 121
>> 122
125 void Run::EndOfRun() 123 void Run::EndOfRun()
126 { 124 {
127 if (numberOfEvent == 0) return; << 125
128 <<
129 // Scatter foil <<
130 126
>> 127 if (numberOfEvent == 0) return;
>> 128
>> 129 //Scatter foil
>> 130
131 G4Material* material = fDetector->GetMateria 131 G4Material* material = fDetector->GetMaterialScatter();
132 G4double length = fDetector->GetThicknessSca << 132 G4double length = fDetector->GetThicknessScatter();
133 G4double density = material->GetDensity(); 133 G4double density = material->GetDensity();
134 G4String partName = fParticle->GetParticleNa 134 G4String partName = fParticle->GetParticleName();
135 135
136 G4cout << "\n ======================== run s 136 G4cout << "\n ======================== run summary ======================\n";
137 137
138 G4int prec = G4cout.precision(3); 138 G4int prec = G4cout.precision(3);
139 << 139
140 G4cout << "\n The run was " << numberOfEvent 140 G4cout << "\n The run was " << numberOfEvent << " " << partName << " of "
141 << G4BestUnit(fEnergy, "Energy") << " << 141 << G4BestUnit(fEnergy,"Energy") << " through "
142 << material->GetName() << " (density: << 142 << G4BestUnit(length,"Length") << " of "
143 << G4endl; << 143 << material->GetName() << " (density: "
>> 144 << G4BestUnit(density,"Volumic Mass") << ")" << G4endl;
144 145
145 G4cout.precision(prec); 146 G4cout.precision(prec);
146 147
147 // normalize histograms 148 // normalize histograms
148 // 149 //
149 G4AnalysisManager* analysisManager = G4Analy << 150 G4AnalysisManager* analysisManager = G4AnalysisManager::Instance();
150 G4double fac = 1. / (double(numberOfEvent)); << 151 G4double fac = 1./(double(numberOfEvent));
151 analysisManager->ScaleH1(1, fac); << 152 analysisManager->ScaleH1(1,fac);
152 analysisManager->ScaleH1(2, fac); << 153 analysisManager->ScaleH1(2,fac);
153 analysisManager->ScaleH1(3, fac); << 154 analysisManager->ScaleH1(3,fac);
154 analysisManager->ScaleH1(5, fac); << 155 analysisManager->ScaleH1(5,fac);
155 analysisManager->ScaleH1(6, fac); << 156 analysisManager->ScaleH1(6,fac);
156 << 157
157 ComputeFluenceError(); 158 ComputeFluenceError();
158 PrintFluence(numberOfEvent); 159 PrintFluence(numberOfEvent);
>> 160
159 } 161 }
160 162
161 //....oooOO0OOooo........oooOO0OOooo........oo 163 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
162 164
163 void Run::InitFluence() 165 void Run::InitFluence()
164 { 166 {
165 // create Fluence histo in any case 167 // create Fluence histo in any case
166 << 168
167 G4AnalysisManager* analysisManager = G4Analy 169 G4AnalysisManager* analysisManager = G4AnalysisManager::Instance();
168 G4int ih = 4; 170 G4int ih = 4;
169 analysisManager->SetH1(ih, 120, 0 * mm, 240 << 171 analysisManager->SetH1(ih, 120, 0*mm, 240*mm, "mm");
170 << 172
171 // construct vectors for fluence distributio << 173 //construct vectors for fluence distribution
172 << 174
173 fNbBins = analysisManager->GetH1Nbins(ih); 175 fNbBins = analysisManager->GetH1Nbins(ih);
174 fDr = analysisManager->GetH1Width(ih); 176 fDr = analysisManager->GetH1Width(ih);
175 fluence.resize(fNbBins, 0.); << 177 fluence.resize(fNbBins, 0.);
176 fluence1.resize(fNbBins, 0.); << 178 fluence1.resize(fNbBins, 0.);
177 fluence2.resize(fNbBins, 0.); 179 fluence2.resize(fNbBins, 0.);
178 fNbEntries.resize(fNbBins, 0); 180 fNbEntries.resize(fNbBins, 0);
179 } 181 }
180 182
181 //....oooOO0OOooo........oooOO0OOooo........oo 183 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
182 184
183 void Run::SumFluence(G4double r, G4double fl) 185 void Run::SumFluence(G4double r, G4double fl)
184 { 186 {
185 G4int ibin = (int)(r / fDr); << 187 G4int ibin = (int)(r/fDr);
186 if (ibin >= fNbBins) return; 188 if (ibin >= fNbBins) return;
187 fNbEntries[ibin]++; 189 fNbEntries[ibin]++;
188 fluence[ibin] += fl; << 190 fluence[ibin] += fl;
189 fluence2[ibin] += fl * fl; << 191 fluence2[ibin] += fl*fl;
190 } 192 }
191 193
192 //....oooOO0OOooo........oooOO0OOooo........oo 194 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
193 195
194 void Run::ComputeFluenceError() 196 void Run::ComputeFluenceError()
195 { 197 {
196 // compute rms << 198 //compute rms
197 << 199
198 G4double ds, variance, rms; << 200 G4double ds,variance,rms;
199 G4double rmean = -0.5 * fDr; << 201 G4double rmean = -0.5*fDr;
200 << 202
201 for (G4int bin = 0; bin < fNbBins; bin++) { << 203 for (G4int bin=0; bin<fNbBins; bin++) {
202 rmean += fDr; << 204 rmean += fDr;
203 ds = twopi * rmean * fDr; << 205 ds = twopi*rmean*fDr;
204 fluence[bin] /= ds; << 206 fluence[bin] /= ds;
205 fluence2[bin] /= (ds * ds); << 207 fluence2[bin] /= (ds*ds);
206 variance = 0.; << 208 variance = 0.;
207 if (fNbEntries[bin] > 0) << 209 if (fNbEntries[bin] > 0)
208 variance = fluence2[bin] - (fluence[bin] << 210 variance = fluence2[bin] - (fluence[bin]*fluence[bin])/fNbEntries[bin];
209 rms = 0.; << 211 rms = 0.;
210 if (variance > 0.) rms = std::sqrt(varianc << 212 if(variance > 0.) rms = std::sqrt(variance);
211 fluence2[bin] = rms; << 213 fluence2[bin] = rms;
212 } << 214 }
213 << 215
214 // normalize to first bins, compute error an << 216 //normalize to first bins, compute error and fill histo
215 << 217
216 G4double rnorm(4 * mm), radius(0.), fnorm(0. << 218 G4double rnorm(4*mm), radius(0.), fnorm(0.), fnorm2(0.);
217 G4int inorm = -1; 219 G4int inorm = -1;
218 do { 220 do {
219 inorm++; << 221 inorm++; radius += fDr; fnorm += fluence[inorm]; fnorm2 += fluence2[inorm];
220 radius += fDr; <<
221 fnorm += fluence[inorm]; <<
222 fnorm2 += fluence2[inorm]; <<
223 } while (radius < rnorm); 222 } while (radius < rnorm);
224 fnorm /= (inorm + 1); << 223 fnorm /= (inorm+1);
225 fnorm2 /= (inorm + 1); << 224 fnorm2 /= (inorm+1);
226 << 225
227 G4double ratio, error; 226 G4double ratio, error;
228 G4double scale = 1. / fnorm; << 227 G4double scale = 1./fnorm;
229 G4double err0 = fnorm2 / fnorm, err1 = 0.; << 228 G4double err0 = fnorm2/fnorm, err1 = 0.;
230 << 229
231 rmean = -0.5 * fDr; << 230 rmean = -0.5*fDr;
232 << 231
233 for (G4int bin = 0; bin < fNbBins; bin++) { << 232 for (G4int bin=0; bin<fNbBins; bin++) {
234 ratio = fluence[bin] * scale; << 233 ratio = fluence[bin]*scale;
235 error = 0.; << 234 error = 0.;
236 if (ratio > 0.) { << 235 if (ratio > 0.) {
237 err1 = fluence2[bin] / fluence[bin]; << 236 err1 = fluence2[bin]/fluence[bin];
238 error = ratio * std::sqrt(err1 * err1 + << 237 error = ratio*std::sqrt(err1*err1 + err0*err0);
239 } << 238 }
240 fluence1[bin] = ratio; << 239 fluence1[bin] = ratio;
241 fluence2[bin] = error; << 240 fluence2[bin] = error;
242 rmean += fDr; << 241 rmean += fDr;
243 G4AnalysisManager* analysisManager = G4Ana << 242 G4AnalysisManager* analysisManager = G4AnalysisManager::Instance();
244 analysisManager->FillH1(4, rmean, ratio); << 243 analysisManager->FillH1(4,rmean,ratio);
245 } 244 }
246 } 245 }
247 246
248 //....oooOO0OOooo........oooOO0OOooo........oo 247 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
249 248
250 #include <fstream> 249 #include <fstream>
251 250
252 void Run::PrintFluence(G4int TotEvents) 251 void Run::PrintFluence(G4int TotEvents)
253 { 252 {
254 G4AnalysisManager* analysisManager = G4Analy << 253 G4AnalysisManager* analysisManager = G4AnalysisManager::Instance();
255 G4String name = analysisManager->GetFileName << 254 G4String name = analysisManager->GetFileName();
256 G4String fileName = name + ".ascii"; 255 G4String fileName = name + ".ascii";
257 std::ofstream File(fileName, std::ios::out); 256 std::ofstream File(fileName, std::ios::out);
258 << 257
259 std::ios::fmtflags mode = File.flags(); << 258 std::ios::fmtflags mode = File.flags();
260 File.setf(std::ios::scientific, std::ios::fl << 259 File.setf( std::ios::scientific, std::ios::floatfield );
261 G4int prec = File.precision(3); 260 G4int prec = File.precision(3);
262 << 261
263 File << " Fluence density distribution \n " << 262 File << " Fluence density distribution \n "
264 << "\n ibin \t radius (mm) \t Nb \t fl 263 << "\n ibin \t radius (mm) \t Nb \t fluence\t norma fl\t rms/nfl (%) \n"
265 << G4endl; 264 << G4endl;
266 265
267 G4double rmean = -0.5 * fDr; << 266 G4double rmean = -0.5*fDr;
268 for (G4int bin = 0; bin < fNbBins; bin++) { << 267 for (G4int bin=0; bin<fNbBins; bin++) {
269 rmean += fDr; << 268 rmean +=fDr;
270 G4double error = 0.; << 269 G4double error = 0.;
271 if (fluence1[bin] > 0.) error = 100 * flue << 270 if (fluence1[bin] > 0.) error = 100*fluence2[bin]/fluence1[bin];
272 File << " " << bin << "\t " << rmean / mm << 271 File << " " << bin << "\t " << rmean/mm << "\t " << fNbEntries[bin]
273 << fluence[bin] / double(TotEvents) < << 272 << "\t " << fluence[bin]/double(TotEvents) << "\t " << fluence1[bin]
>> 273 << "\t " << error
>> 274 << G4endl;
274 } 275 }
275 << 276
276 // restaure default formats 277 // restaure default formats
277 File.setf(mode, std::ios::floatfield); << 278 File.setf(mode,std::ios::floatfield);
278 File.precision(prec); << 279 File.precision(prec);
279 } 280 }
280 281
281 //....oooOO0OOooo........oooOO0OOooo........oo 282 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
>> 283
282 284