| 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 // G4StatAnalysis inline methods implementatio <<
27 // 26 //
28 // Author: J.Madsen, 25.10.2018 << 27 //
29 // ------------------------------------------- << 28 //
>> 29 // ----------------------------------------------------------------------
>> 30 // Class typename G4StatAnalysis
>> 31 //
>> 32 // Class description:
>> 33 //
>> 34 // Class for statistical analysis of random variable
>> 35 //
>> 36 // Adapted
>> 37 // Lux, I.
>> 38 // Monte Carlo particle transport methods: neutron and photon
>> 39 // calculations/authors, Ivan Lux and Laszlo Koblinger.
>> 40 // ISBN 0-8493-6074-9
>> 41 // 1. Neutron transport theory. 2. Photon transport theory.
>> 42 // 3. Monte Carlo method. I. Koblinger, Laszlo. II. Title.
>> 43 // QC793.5.N4628L88 1990
>> 44 // 530.1 '38—dc20
>> 45 //
>> 46 // https://gnssn.iaea.org/NSNI/Shared%20Documents/OPEN%20Shared%20Files/MonteCarloParticleTransportMethodsNeutronAndPhotonCalculations.pdf
>> 47 //
>> 48 //
30 49
31 #include <cmath> <<
32 #include <fstream> <<
33 #include <iomanip> <<
34 #include <iostream> 50 #include <iostream>
>> 51 #include <iomanip>
35 #include <limits> 52 #include <limits>
>> 53 #include <fstream>
>> 54 #include <cmath>
36 55
37 #include "G4Timer.hh" <<
38 #include "G4ios.hh" <<
39 #include "globals.hh" 56 #include "globals.hh"
40 #include "tls.hh" 57 #include "tls.hh"
>> 58 #include "G4Timer.hh"
>> 59 #include "G4ios.hh"
41 60
42 #include "G4Allocator.hh" <<
43 #include "G4Types.hh" 61 #include "G4Types.hh"
>> 62 #include "G4Allocator.hh"
44 63
45 //-------------------------------------------- 64 //----------------------------------------------------------------------------//
46 65
47 G4StatAnalysis::G4StatAnalysis() {} << 66 G4StatAnalysis::G4StatAnalysis()
>> 67 : fSum1(0.0),
>> 68 fSum2(0.0),
>> 69 fHits(0),
>> 70 fZero(0)
>> 71 { }
48 72
49 //-------------------------------------------- 73 //----------------------------------------------------------------------------//
50 74
51 G4double G4StatAnalysis::GetMean() const 75 G4double G4StatAnalysis::GetMean() const
52 { 76 {
53 return (fHits > 0) ? fSum1 / ((G4double) fHi << 77 return (fHits > 0) ? fSum1/((G4double) fHits) : 0.;
54 } 78 }
55 79
56 //-------------------------------------------- 80 //----------------------------------------------------------------------------//
57 81
58 const G4double& G4StatAnalysis::GetSum() const << 82 const G4double& G4StatAnalysis::GetSum() const
>> 83 {
>> 84 return fSum1;
>> 85 }
59 86
60 //-------------------------------------------- 87 //----------------------------------------------------------------------------//
61 88
62 const G4double& G4StatAnalysis::GetSumSquared( << 89 const G4double& G4StatAnalysis::GetSumSquared() const
>> 90 {
>> 91 return fSum2;
>> 92 }
63 93
64 //-------------------------------------------- 94 //----------------------------------------------------------------------------//
65 95
66 const G4double& G4StatAnalysis::GetSum1() cons << 96 const G4double& G4StatAnalysis::GetSum1() const
>> 97 {
>> 98 return fSum1;
>> 99 }
67 100
68 //-------------------------------------------- 101 //----------------------------------------------------------------------------//
69 102
70 const G4double& G4StatAnalysis::GetSum2() cons << 103 const G4double& G4StatAnalysis::GetSum2() const
>> 104 {
>> 105 return fSum2;
>> 106 }
71 107
72 //-------------------------------------------- 108 //----------------------------------------------------------------------------//
73 109
74 const G4int& G4StatAnalysis::GetHits() const { << 110 const G4int& G4StatAnalysis::GetHits() const
>> 111 {
>> 112 return fHits;
>> 113 }
75 114
76 //-------------------------------------------- 115 //----------------------------------------------------------------------------//
77 116
78 G4int G4StatAnalysis::GetNumNonZero() const { << 117 G4int G4StatAnalysis::GetNumNonZero() const
>> 118 {
>> 119 return fHits - fZero;
>> 120 }
79 121
80 //-------------------------------------------- 122 //----------------------------------------------------------------------------//
81 123
82 G4int G4StatAnalysis::GetNumZero() const { ret << 124 G4int G4StatAnalysis::GetNumZero() const
>> 125 {
>> 126 return fZero;
>> 127 }
83 128
84 //-------------------------------------------- 129 //----------------------------------------------------------------------------//
85 130
86 void G4StatAnalysis::SetSum(const G4double& va << 131 void G4StatAnalysis::SetSum(const G4double& val)
>> 132 {
>> 133 fSum1 = val;
>> 134 }
87 135
88 //-------------------------------------------- 136 //----------------------------------------------------------------------------//
89 137
90 void G4StatAnalysis::SetSumSquared(const G4dou << 138 void G4StatAnalysis::SetSumSquared(const G4double& val)
>> 139 {
>> 140 fSum2 = val;
>> 141 }
91 142
92 //-------------------------------------------- 143 //----------------------------------------------------------------------------//
93 144
94 void G4StatAnalysis::SetSum1(const G4double& v << 145 void G4StatAnalysis::SetSum1(const G4double& val)
>> 146 {
>> 147 fSum1 = val;
>> 148 }
95 149
96 //-------------------------------------------- 150 //----------------------------------------------------------------------------//
97 151
98 void G4StatAnalysis::SetSum2(const G4double& v << 152 void G4StatAnalysis::SetSum2(const G4double& val)
>> 153 {
>> 154 fSum2 = val;
>> 155 }
99 156
100 //-------------------------------------------- 157 //----------------------------------------------------------------------------//
101 158
102 void G4StatAnalysis::SetHits(const G4int& val) << 159 void G4StatAnalysis::SetHits(const G4int& val)
>> 160 {
>> 161 fHits = val;
>> 162 }
103 163
104 //-------------------------------------------- 164 //----------------------------------------------------------------------------//
105 165
106 void G4StatAnalysis::SetZero(const G4int& val) << 166 void G4StatAnalysis::SetZero(const G4int& val)
>> 167 {
>> 168 fZero = val;
>> 169 }
107 170
108 //-------------------------------------------- 171 //----------------------------------------------------------------------------//
109 172
110 G4double G4StatAnalysis::GetFOM() const 173 G4double G4StatAnalysis::GetFOM() const
111 { 174 {
112 G4double elapsed_time = this->GetCpuTime(); << 175 G4double elapsed_time = this->GetCpuTime();
113 G4double relative_err = this->GetRelativeErr << 176 G4double relative_err = this->GetRelativeError();
114 // lambda for equation clarity (will be inli << 177 // lambda for equation clarity (will be inlined)
115 auto compute_figure_of_merit = [&]() { << 178 auto compute_figure_of_merit = [&] ()
116 return (1.0 / (relative_err * relative_err << 179 {
117 }; << 180 return ( 1.0 / ( relative_err * relative_err ) / elapsed_time );
118 return (std::fabs(relative_err) > 0.0 && ela << 181 };
119 ? compute_figure_of_merit() << 182 return (std::fabs(relative_err) > 0.0 && elapsed_time > 0.0)
120 : ((fHits > 0) ? 1.0 : 0.0); << 183 ? compute_figure_of_merit() : ((fHits > 0) ? 1.0 : 0.0);
121 } 184 }
122 185
123 //-------------------------------------------- 186 //----------------------------------------------------------------------------//
124 187
125 G4double G4StatAnalysis::GetRelativeError() co 188 G4double G4StatAnalysis::GetRelativeError() const
126 { 189 {
127 // lambda for equation clarity (will be inli << 190 // lambda for equation clarity (will be inlined)
128 auto compute_relative_error = [&]() { << 191 auto compute_relative_error = [&] ()
129 return (GetStdDev() / GetMean() / std::sqr << 192 {
130 }; << 193 return ( GetStdDev() / GetMean() / std::sqrt((G4double) fHits) );
131 return (std::fabs(GetMean()) > 0 && fHits > << 194 };
132 << 195 return (std::fabs(GetMean()) > 0 && fHits > 0)
>> 196 ? compute_relative_error() : ((fHits > 0) ? 1.0 : 0.0);
133 } 197 }
134 198
135 //-------------------------------------------- 199 //----------------------------------------------------------------------------//
136 200
137 G4double G4StatAnalysis::GetStdDev() const 201 G4double G4StatAnalysis::GetStdDev() const
138 { 202 {
139 return ::sqrt(std::fabs(GetVariance())); << 203 return ::sqrt(std::fabs(GetVariance()));
140 } 204 }
141 205
142 //-------------------------------------------- 206 //----------------------------------------------------------------------------//
143 207
144 G4double G4StatAnalysis::GetVariance() const 208 G4double G4StatAnalysis::GetVariance() const
145 { 209 {
146 // lambda for equation clarity (will be inli << 210 // lambda for equation clarity (will be inlined)
147 auto compute_variance = [&]() { << 211 auto compute_variance = [&] ()
148 return ((fSum2 - (std::pow(fSum1, 2.0) / f << 212 {
149 (((G4double) fHits) - 1.0)); << 213 return ((fSum2 - (std::pow(fSum1, 2.0)/fHits))/(((G4double) fHits) - 1.0));
150 }; << 214 };
151 return (fHits > 1) ? compute_variance() : 0. << 215 return (fHits > 1) ? compute_variance() : 0.0;
152 } 216 }
153 217
154 //-------------------------------------------- 218 //----------------------------------------------------------------------------//
155 219
156 G4double G4StatAnalysis::GetCoeffVariation() c 220 G4double G4StatAnalysis::GetCoeffVariation() const
157 { 221 {
158 // lambda for equation clarity (will be inli << 222 // lambda for equation clarity (will be inlined)
159 auto coefficient_of_variation = [&]() { << 223 auto coefficient_of_variation = [&] ()
160 G4double hits = fHits; << 224 {
161 return ::sqrt((hits / (hits - 1.0)) * << 225 G4double hits = fHits;
162 ((fSum2 / (fSum1 * fSum1)) - << 226 return ::sqrt(
163 }; << 227 (hits / (hits-1.0)) *
164 return (fHits > 1) ? coefficient_of_variatio << 228 ( (fSum2/(fSum1*fSum1)) - (1.0/hits))
165 // return (fHits > 0 && fabs(fSum1) > 0.0) << 229 );
166 // ? (100.0*GetStdDev()/GetMean()) : << 230 };
>> 231 return (fHits > 1) ? coefficient_of_variation() : 0.0;
>> 232 //return (fHits > 0 && fabs(fSum1) > 0.0)
>> 233 // ? (100.0*GetStdDev()/GetMean()) : 0.0;
>> 234
167 } 235 }
168 236
169 //-------------------------------------------- 237 //----------------------------------------------------------------------------//
170 238
171 G4double G4StatAnalysis::GetEfficiency() const 239 G4double G4StatAnalysis::GetEfficiency() const
172 { 240 {
173 G4double hits = fHits; << 241 G4double hits = fHits;
174 G4double nzero = fHits - fZero; << 242 G4double nzero = fHits - fZero;
175 return (fHits > 0) ? (nzero / hits) : 0.0; << 243 return (fHits > 0) ? (nzero/hits) : 0.0;
176 } 244 }
177 245
178 //-------------------------------------------- 246 //----------------------------------------------------------------------------//
179 247
180 G4double G4StatAnalysis::GetR2Int() const 248 G4double G4StatAnalysis::GetR2Int() const
181 { 249 {
182 G4double hits = fHits; << 250 G4double hits = fHits;
183 return (fHits > 0) << 251 return (fHits > 0)
184 ? (fSum2 / (fSum1 * fSum1)) - 1.0 / << 252 ? (fSum2 / (fSum1 * fSum1)) - 1.0/(GetEfficiency() * hits)
185 : 0.0; << 253 : 0.0;
186 } 254 }
187 255
188 //-------------------------------------------- 256 //----------------------------------------------------------------------------//
189 257
190 G4double G4StatAnalysis::GetR2Eff() const 258 G4double G4StatAnalysis::GetR2Eff() const
191 { 259 {
192 G4double hits = fHits; << 260 G4double hits = fHits;
193 return (fHits > 0) ? (1.0 - GetEfficiency()) << 261 return (fHits > 0)
>> 262 ? (1.0 - GetEfficiency()) / (GetEfficiency() * hits)
>> 263 : 0.0;
194 } 264 }
195 265
196 //-------------------------------------------- 266 //----------------------------------------------------------------------------//
197 267
198 G4StatAnalysis::operator G4double() const { re << 268 G4StatAnalysis::operator G4double() const
>> 269 {
>> 270 return this->GetSum();
>> 271 }
199 272
200 //-------------------------------------------- 273 //----------------------------------------------------------------------------//
201 274
202 void G4StatAnalysis::Reset() 275 void G4StatAnalysis::Reset()
203 { 276 {
204 fHits = 0; << 277 fHits = 0;
205 fZero = 0; << 278 fZero = 0;
206 fSum1 = 0.0; << 279 fSum1 = 0.0;
207 fSum2 = 0.0; << 280 fSum2 = 0.0;
208 } 281 }
209 282
210 //-------------------------------------------- 283 //----------------------------------------------------------------------------//
211 284
212 void G4StatAnalysis::Add(const G4double& val, 285 void G4StatAnalysis::Add(const G4double& val, const G4double& weight)
213 { 286 {
214 fHits += 1; << 287 fHits += 1;
215 fSum1 += val * weight; << 288 fSum1 += val * weight;
216 fSum2 += val * val * weight; << 289 fSum2 += val * val * weight;
217 if(std::fabs(val * weight) < << 290 if(std::fabs(val*weight) < std::fabs(GetMean() * std::numeric_limits<double>::epsilon()))
218 std::fabs(GetMean() * std::numeric_limits << 291 fZero += 1;
219 fZero += 1; <<
220 } 292 }
221 293
222 //-------------------------------------------- 294 //----------------------------------------------------------------------------//
223 295
224 void G4StatAnalysis::Rescale(const G4double& f 296 void G4StatAnalysis::Rescale(const G4double& factor)
225 { 297 {
226 fSum1 *= factor; << 298 fSum1 *= factor;
227 fSum2 *= factor * factor; << 299 fSum2 *= factor * factor;
228 } 300 }
229 301
230 //-------------------------------------------- 302 //----------------------------------------------------------------------------//
231 303
232 void G4StatAnalysis::PrintInfo(std::ostream& o 304 void G4StatAnalysis::PrintInfo(std::ostream& os, const std::string& tab) const
233 { 305 {
234 G4int _hits = this->GetHits(); << 306 G4int _hits = this->GetHits();
235 G4double _sum = this->GetSum(); << 307 G4double _sum = this->GetSum();
236 G4double _sigma = this->GetStdDev(); << 308 G4double _sigma = this->GetStdDev();
237 G4double _coeff = this->GetCoeffVariation(); << 309 G4double _coeff = this->GetCoeffVariation();
238 G4double _error = this->GetRelativeError(); << 310 G4double _error = this->GetRelativeError();
239 G4double _eff = this->GetEfficiency(); << 311 G4double _eff = this->GetEfficiency();
240 G4double _fom = this->GetFOM(); << 312 G4double _fom = this->GetFOM();
241 G4double _r2int = this->GetR2Int(); << 313 G4double _r2int = this->GetR2Int();
242 G4double _r2eff = this->GetR2Eff(); << 314 G4double _r2eff = this->GetR2Eff();
243 << 315
244 using std::fixed; << 316 using std::setprecision;
245 using std::ios; << 317 using std::setw;
246 using std::left; << 318 using std::scientific;
247 using std::right; << 319 using std::fixed;
248 using std::scientific; << 320 using std::left;
249 using std::setprecision; << 321 using std::right;
250 using std::setw; << 322 using std::ios;
251 << 323
252 std::stringstream ss; << 324 std::stringstream ss;
253 ss << tab //<< scientific << 325 ss << tab //<< scientific
254 << setprecision((G4int)os.precision()) << << 326 << setprecision(os.precision()) << right << _sum
255 << " [sigma: " << right << _sigma << left << 327 << left << " [sigma: " << right << _sigma
256 << _error << left << " | coeff: " << righ << 328 << left << " | error: " << right << _error
257 << " | eff: " << right << _eff << left << << 329 << left << " | coeff: " << right << _coeff
258 << left << " | r2int: " << right << _r2in << 330 << left << " | eff: " << right << _eff
259 << _r2eff << left << " | hits: " << right << 331 << left << " | fom: " << right << _fom
>> 332 << left << " | r2int: " << right << _r2int
>> 333 << left << " | r2eff: " << right << _r2eff
>> 334 << left << " | hits: " << right << _hits
>> 335 << left << " ]";
260 336
261 os << ss.str(); << 337 os << ss.str();
262 } 338 }
263 339
264 //-------------------------------------------- 340 //----------------------------------------------------------------------------//
265 341
266 G4double G4StatAnalysis::GetCpuTime() const 342 G4double G4StatAnalysis::GetCpuTime() const
267 { 343 {
268 tms* startTime = GetCpuClock(); << 344 tms* startTime = GetCpuClock();
269 tms endTime; << 345 tms endTime;
270 times(&endTime); << 346 times(&endTime);
271 return ((endTime.tms_stime - startTime->tms_ << 347 return ((endTime.tms_stime - startTime->tms_stime) +
272 (endTime.tms_utime - startTime->tms_ << 348 (endTime.tms_utime - startTime->tms_utime)) / sysconf(_SC_CLK_TCK);
273 sysconf(_SC_CLK_TCK); <<
274 } 349 }
275 350
276 //-------------------------------------------- 351 //----------------------------------------------------------------------------//
277 352
278 G4StatAnalysis& G4StatAnalysis::operator+=(con << 353 G4StatAnalysis&
>> 354 G4StatAnalysis::operator+=(const G4double& _val)
279 { 355 {
280 this->Add(_val); << 356 this->Add(_val);
281 return *this; << 357 return *this;
282 } 358 }
283 359
284 //-------------------------------------------- 360 //----------------------------------------------------------------------------//
285 361
286 G4StatAnalysis& G4StatAnalysis::operator/=(con << 362 G4StatAnalysis&
>> 363 G4StatAnalysis::operator/=(const G4double& _val)
287 { 364 {
288 fSum1 /= _val; << 365 fSum1 /= _val;
289 fSum2 /= (_val * _val); << 366 fSum2 /= (_val*_val);
290 return *this; << 367 return *this;
291 } 368 }
292 369
293 //-------------------------------------------- 370 //----------------------------------------------------------------------------//
294 371
295 G4StatAnalysis& G4StatAnalysis::operator+=(con << 372 G4StatAnalysis&
>> 373 G4StatAnalysis::operator+=(const G4StatAnalysis& rhs)
296 { 374 {
297 fHits += rhs.fHits; << 375 fHits += rhs.fHits;
298 fSum1 += rhs.fSum1; << 376 fSum1 += rhs.fSum1;
299 fSum2 += rhs.fSum2; << 377 fSum2 += rhs.fSum2;
300 fZero += rhs.fZero; << 378 fZero += rhs.fZero;
301 return *this; << 379 return *this;
302 } 380 }
303 381
304 //-------------------------------------------- 382 //----------------------------------------------------------------------------//
305 383
306 G4StatAnalysis& G4StatAnalysis::operator-=(con << 384 G4StatAnalysis&
>> 385 G4StatAnalysis::operator-=(const G4StatAnalysis& rhs)
307 { 386 {
308 fHits -= rhs.fHits; << 387 fHits -= rhs.fHits;
309 fSum1 -= rhs.fSum1; << 388 fSum1 -= rhs.fSum1;
310 fSum2 -= rhs.fSum2; << 389 fSum2 -= rhs.fSum2;
311 fZero -= rhs.fZero; << 390 fZero -= rhs.fZero;
312 return *this; << 391 return *this;
313 } 392 }
314 393
315 //-------------------------------------------- 394 //----------------------------------------------------------------------------//
316 395
317 inline G4Allocator<G4StatAnalysis>*& _aStatAna 396 inline G4Allocator<G4StatAnalysis>*& _aStatAnalysisAllocator_G4MT_TLS_()
318 { 397 {
319 G4ThreadLocalStatic G4Allocator<G4StatAnalys << 398 G4ThreadLocalStatic G4Allocator<G4StatAnalysis>* _instance
320 new G4Allocator<G4StatAnalysis>(); << 399 = new G4Allocator<G4StatAnalysis>();
321 return _instance; << 400 return _instance;
322 } 401 }
323 402
324 //-------------------------------------------- 403 //----------------------------------------------------------------------------//
325 404
326 void* G4StatAnalysis::operator new(std::size_t << 405 void* G4StatAnalysis::operator new(size_t)
327 { 406 {
328 G4Allocator<G4StatAnalysis>& _allocator = << 407 G4Allocator<G4StatAnalysis>& _allocator = *_aStatAnalysisAllocator_G4MT_TLS_();
329 *_aStatAnalysisAllocator_G4MT_TLS_(); << 408 return (void*) _allocator.MallocSingle();
330 return (void*) _allocator.MallocSingle(); <<
331 } 409 }
332 410
333 //-------------------------------------------- 411 //----------------------------------------------------------------------------//
334 412
335 void G4StatAnalysis::operator delete(void* _pt 413 void G4StatAnalysis::operator delete(void* _ptr)
336 { 414 {
337 G4Allocator<G4StatAnalysis>& _allocator = << 415 G4Allocator<G4StatAnalysis>& _allocator = *_aStatAnalysisAllocator_G4MT_TLS_();
338 *_aStatAnalysisAllocator_G4MT_TLS_(); << 416 _allocator.FreeSingle( (G4StatAnalysis*) _ptr);
339 _allocator.FreeSingle((G4StatAnalysis*) _ptr <<
340 } 417 }
341 418
342 //-------------------------------------------- 419 //----------------------------------------------------------------------------//
343 420