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