| 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 26 // G4StatAnalysis inline methods implementation
27 // 27 //
28 // Author: J.Madsen, 25.10.2018 28 // Author: J.Madsen, 25.10.2018
29 // ------------------------------------------- 29 // --------------------------------------------------------------------
30 30
31 #include <cmath> 31 #include <cmath>
32 #include <fstream> 32 #include <fstream>
33 #include <iomanip> 33 #include <iomanip>
34 #include <iostream> 34 #include <iostream>
35 #include <limits> 35 #include <limits>
36 36
37 #include "G4Timer.hh" 37 #include "G4Timer.hh"
38 #include "G4ios.hh" 38 #include "G4ios.hh"
39 #include "globals.hh" 39 #include "globals.hh"
40 #include "tls.hh" 40 #include "tls.hh"
41 41
42 #include "G4Allocator.hh" 42 #include "G4Allocator.hh"
43 #include "G4Types.hh" 43 #include "G4Types.hh"
44 44
45 //-------------------------------------------- 45 //----------------------------------------------------------------------------//
46 46
47 G4StatAnalysis::G4StatAnalysis() {} 47 G4StatAnalysis::G4StatAnalysis() {}
48 48
49 //-------------------------------------------- 49 //----------------------------------------------------------------------------//
50 50
51 G4double G4StatAnalysis::GetMean() const 51 G4double G4StatAnalysis::GetMean() const
52 { 52 {
53 return (fHits > 0) ? fSum1 / ((G4double) fHi 53 return (fHits > 0) ? fSum1 / ((G4double) fHits) : 0.;
54 } 54 }
55 55
56 //-------------------------------------------- 56 //----------------------------------------------------------------------------//
57 57
58 const G4double& G4StatAnalysis::GetSum() const 58 const G4double& G4StatAnalysis::GetSum() const { return fSum1; }
59 59
60 //-------------------------------------------- 60 //----------------------------------------------------------------------------//
61 61
62 const G4double& G4StatAnalysis::GetSumSquared( 62 const G4double& G4StatAnalysis::GetSumSquared() const { return fSum2; }
63 63
64 //-------------------------------------------- 64 //----------------------------------------------------------------------------//
65 65
66 const G4double& G4StatAnalysis::GetSum1() cons 66 const G4double& G4StatAnalysis::GetSum1() const { return fSum1; }
67 67
68 //-------------------------------------------- 68 //----------------------------------------------------------------------------//
69 69
70 const G4double& G4StatAnalysis::GetSum2() cons 70 const G4double& G4StatAnalysis::GetSum2() const { return fSum2; }
71 71
72 //-------------------------------------------- 72 //----------------------------------------------------------------------------//
73 73
74 const G4int& G4StatAnalysis::GetHits() const { 74 const G4int& G4StatAnalysis::GetHits() const { return fHits; }
75 75
76 //-------------------------------------------- 76 //----------------------------------------------------------------------------//
77 77
78 G4int G4StatAnalysis::GetNumNonZero() const { 78 G4int G4StatAnalysis::GetNumNonZero() const { return fHits - fZero; }
79 79
80 //-------------------------------------------- 80 //----------------------------------------------------------------------------//
81 81
82 G4int G4StatAnalysis::GetNumZero() const { ret 82 G4int G4StatAnalysis::GetNumZero() const { return fZero; }
83 83
84 //-------------------------------------------- 84 //----------------------------------------------------------------------------//
85 85
86 void G4StatAnalysis::SetSum(const G4double& va 86 void G4StatAnalysis::SetSum(const G4double& val) { fSum1 = val; }
87 87
88 //-------------------------------------------- 88 //----------------------------------------------------------------------------//
89 89
90 void G4StatAnalysis::SetSumSquared(const G4dou 90 void G4StatAnalysis::SetSumSquared(const G4double& val) { fSum2 = val; }
91 91
92 //-------------------------------------------- 92 //----------------------------------------------------------------------------//
93 93
94 void G4StatAnalysis::SetSum1(const G4double& v 94 void G4StatAnalysis::SetSum1(const G4double& val) { fSum1 = val; }
95 95
96 //-------------------------------------------- 96 //----------------------------------------------------------------------------//
97 97
98 void G4StatAnalysis::SetSum2(const G4double& v 98 void G4StatAnalysis::SetSum2(const G4double& val) { fSum2 = val; }
99 99
100 //-------------------------------------------- 100 //----------------------------------------------------------------------------//
101 101
102 void G4StatAnalysis::SetHits(const G4int& val) 102 void G4StatAnalysis::SetHits(const G4int& val) { fHits = val; }
103 103
104 //-------------------------------------------- 104 //----------------------------------------------------------------------------//
105 105
106 void G4StatAnalysis::SetZero(const G4int& val) 106 void G4StatAnalysis::SetZero(const G4int& val) { fZero = val; }
107 107
108 //-------------------------------------------- 108 //----------------------------------------------------------------------------//
109 109
110 G4double G4StatAnalysis::GetFOM() const 110 G4double G4StatAnalysis::GetFOM() const
111 { 111 {
112 G4double elapsed_time = this->GetCpuTime(); 112 G4double elapsed_time = this->GetCpuTime();
113 G4double relative_err = this->GetRelativeErr 113 G4double relative_err = this->GetRelativeError();
114 // lambda for equation clarity (will be inli 114 // lambda for equation clarity (will be inlined)
115 auto compute_figure_of_merit = [&]() { 115 auto compute_figure_of_merit = [&]() {
116 return (1.0 / (relative_err * relative_err 116 return (1.0 / (relative_err * relative_err) / elapsed_time);
117 }; 117 };
118 return (std::fabs(relative_err) > 0.0 && ela 118 return (std::fabs(relative_err) > 0.0 && elapsed_time > 0.0)
119 ? compute_figure_of_merit() 119 ? compute_figure_of_merit()
120 : ((fHits > 0) ? 1.0 : 0.0); 120 : ((fHits > 0) ? 1.0 : 0.0);
121 } 121 }
122 122
123 //-------------------------------------------- 123 //----------------------------------------------------------------------------//
124 124
125 G4double G4StatAnalysis::GetRelativeError() co 125 G4double G4StatAnalysis::GetRelativeError() const
126 { 126 {
127 // lambda for equation clarity (will be inli 127 // lambda for equation clarity (will be inlined)
128 auto compute_relative_error = [&]() { 128 auto compute_relative_error = [&]() {
129 return (GetStdDev() / GetMean() / std::sqr 129 return (GetStdDev() / GetMean() / std::sqrt((G4double) fHits));
130 }; 130 };
131 return (std::fabs(GetMean()) > 0 && fHits > 131 return (std::fabs(GetMean()) > 0 && fHits > 0) ? compute_relative_error()
132 132 : ((fHits > 0) ? 1.0 : 0.0);
133 } 133 }
134 134
135 //-------------------------------------------- 135 //----------------------------------------------------------------------------//
136 136
137 G4double G4StatAnalysis::GetStdDev() const 137 G4double G4StatAnalysis::GetStdDev() const
138 { 138 {
139 return ::sqrt(std::fabs(GetVariance())); 139 return ::sqrt(std::fabs(GetVariance()));
140 } 140 }
141 141
142 //-------------------------------------------- 142 //----------------------------------------------------------------------------//
143 143
144 G4double G4StatAnalysis::GetVariance() const 144 G4double G4StatAnalysis::GetVariance() const
145 { 145 {
146 // lambda for equation clarity (will be inli 146 // lambda for equation clarity (will be inlined)
147 auto compute_variance = [&]() { 147 auto compute_variance = [&]() {
148 return ((fSum2 - (std::pow(fSum1, 2.0) / f 148 return ((fSum2 - (std::pow(fSum1, 2.0) / fHits)) /
149 (((G4double) fHits) - 1.0)); 149 (((G4double) fHits) - 1.0));
150 }; 150 };
151 return (fHits > 1) ? compute_variance() : 0. 151 return (fHits > 1) ? compute_variance() : 0.0;
152 } 152 }
153 153
154 //-------------------------------------------- 154 //----------------------------------------------------------------------------//
155 155
156 G4double G4StatAnalysis::GetCoeffVariation() c 156 G4double G4StatAnalysis::GetCoeffVariation() const
157 { 157 {
158 // lambda for equation clarity (will be inli 158 // lambda for equation clarity (will be inlined)
159 auto coefficient_of_variation = [&]() { 159 auto coefficient_of_variation = [&]() {
160 G4double hits = fHits; 160 G4double hits = fHits;
161 return ::sqrt((hits / (hits - 1.0)) * 161 return ::sqrt((hits / (hits - 1.0)) *
162 ((fSum2 / (fSum1 * fSum1)) - 162 ((fSum2 / (fSum1 * fSum1)) - (1.0 / hits)));
163 }; 163 };
164 return (fHits > 1) ? coefficient_of_variatio 164 return (fHits > 1) ? coefficient_of_variation() : 0.0;
165 // return (fHits > 0 && fabs(fSum1) > 0.0) 165 // return (fHits > 0 && fabs(fSum1) > 0.0)
166 // ? (100.0*GetStdDev()/GetMean()) : 166 // ? (100.0*GetStdDev()/GetMean()) : 0.0;
167 } 167 }
168 168
169 //-------------------------------------------- 169 //----------------------------------------------------------------------------//
170 170
171 G4double G4StatAnalysis::GetEfficiency() const 171 G4double G4StatAnalysis::GetEfficiency() const
172 { 172 {
173 G4double hits = fHits; 173 G4double hits = fHits;
174 G4double nzero = fHits - fZero; 174 G4double nzero = fHits - fZero;
175 return (fHits > 0) ? (nzero / hits) : 0.0; 175 return (fHits > 0) ? (nzero / hits) : 0.0;
176 } 176 }
177 177
178 //-------------------------------------------- 178 //----------------------------------------------------------------------------//
179 179
180 G4double G4StatAnalysis::GetR2Int() const 180 G4double G4StatAnalysis::GetR2Int() const
181 { 181 {
182 G4double hits = fHits; 182 G4double hits = fHits;
183 return (fHits > 0) 183 return (fHits > 0)
184 ? (fSum2 / (fSum1 * fSum1)) - 1.0 / 184 ? (fSum2 / (fSum1 * fSum1)) - 1.0 / (GetEfficiency() * hits)
185 : 0.0; 185 : 0.0;
186 } 186 }
187 187
188 //-------------------------------------------- 188 //----------------------------------------------------------------------------//
189 189
190 G4double G4StatAnalysis::GetR2Eff() const 190 G4double G4StatAnalysis::GetR2Eff() const
191 { 191 {
192 G4double hits = fHits; 192 G4double hits = fHits;
193 return (fHits > 0) ? (1.0 - GetEfficiency()) 193 return (fHits > 0) ? (1.0 - GetEfficiency()) / (GetEfficiency() * hits) : 0.0;
194 } 194 }
195 195
196 //-------------------------------------------- 196 //----------------------------------------------------------------------------//
197 197
198 G4StatAnalysis::operator G4double() const { re 198 G4StatAnalysis::operator G4double() const { return this->GetSum(); }
199 199
200 //-------------------------------------------- 200 //----------------------------------------------------------------------------//
201 201
202 void G4StatAnalysis::Reset() 202 void G4StatAnalysis::Reset()
203 { 203 {
204 fHits = 0; 204 fHits = 0;
205 fZero = 0; 205 fZero = 0;
206 fSum1 = 0.0; 206 fSum1 = 0.0;
207 fSum2 = 0.0; 207 fSum2 = 0.0;
208 } 208 }
209 209
210 //-------------------------------------------- 210 //----------------------------------------------------------------------------//
211 211
212 void G4StatAnalysis::Add(const G4double& val, 212 void G4StatAnalysis::Add(const G4double& val, const G4double& weight)
213 { 213 {
214 fHits += 1; 214 fHits += 1;
215 fSum1 += val * weight; 215 fSum1 += val * weight;
216 fSum2 += val * val * weight; 216 fSum2 += val * val * weight;
217 if(std::fabs(val * weight) < 217 if(std::fabs(val * weight) <
218 std::fabs(GetMean() * std::numeric_limits 218 std::fabs(GetMean() * std::numeric_limits<double>::epsilon()))
219 fZero += 1; 219 fZero += 1;
220 } 220 }
221 221
222 //-------------------------------------------- 222 //----------------------------------------------------------------------------//
223 223
224 void G4StatAnalysis::Rescale(const G4double& f 224 void G4StatAnalysis::Rescale(const G4double& factor)
225 { 225 {
226 fSum1 *= factor; 226 fSum1 *= factor;
227 fSum2 *= factor * factor; 227 fSum2 *= factor * factor;
228 } 228 }
229 229
230 //-------------------------------------------- 230 //----------------------------------------------------------------------------//
231 231
232 void G4StatAnalysis::PrintInfo(std::ostream& o 232 void G4StatAnalysis::PrintInfo(std::ostream& os, const std::string& tab) const
233 { 233 {
234 G4int _hits = this->GetHits(); 234 G4int _hits = this->GetHits();
235 G4double _sum = this->GetSum(); 235 G4double _sum = this->GetSum();
236 G4double _sigma = this->GetStdDev(); 236 G4double _sigma = this->GetStdDev();
237 G4double _coeff = this->GetCoeffVariation(); 237 G4double _coeff = this->GetCoeffVariation();
238 G4double _error = this->GetRelativeError(); 238 G4double _error = this->GetRelativeError();
239 G4double _eff = this->GetEfficiency(); 239 G4double _eff = this->GetEfficiency();
240 G4double _fom = this->GetFOM(); 240 G4double _fom = this->GetFOM();
241 G4double _r2int = this->GetR2Int(); 241 G4double _r2int = this->GetR2Int();
242 G4double _r2eff = this->GetR2Eff(); 242 G4double _r2eff = this->GetR2Eff();
243 243
244 using std::fixed; 244 using std::fixed;
245 using std::ios; 245 using std::ios;
246 using std::left; 246 using std::left;
247 using std::right; 247 using std::right;
248 using std::scientific; 248 using std::scientific;
249 using std::setprecision; 249 using std::setprecision;
250 using std::setw; 250 using std::setw;
251 251
252 std::stringstream ss; 252 std::stringstream ss;
253 ss << tab //<< scientific 253 ss << tab //<< scientific
254 << setprecision((G4int)os.precision()) << 254 << setprecision((G4int)os.precision()) << right << _sum << left
255 << " [sigma: " << right << _sigma << left 255 << " [sigma: " << right << _sigma << left << " | error: " << right
256 << _error << left << " | coeff: " << righ 256 << _error << left << " | coeff: " << right << _coeff << left
257 << " | eff: " << right << _eff << left << 257 << " | eff: " << right << _eff << left << " | fom: " << right << _fom
258 << left << " | r2int: " << right << _r2in 258 << left << " | r2int: " << right << _r2int << left << " | r2eff: " << right
259 << _r2eff << left << " | hits: " << right 259 << _r2eff << left << " | hits: " << right << _hits << left << " ]";
260 260
261 os << ss.str(); 261 os << ss.str();
262 } 262 }
263 263
264 //-------------------------------------------- 264 //----------------------------------------------------------------------------//
265 265
266 G4double G4StatAnalysis::GetCpuTime() const 266 G4double G4StatAnalysis::GetCpuTime() const
267 { 267 {
268 tms* startTime = GetCpuClock(); 268 tms* startTime = GetCpuClock();
269 tms endTime; 269 tms endTime;
270 times(&endTime); 270 times(&endTime);
271 return ((endTime.tms_stime - startTime->tms_ 271 return ((endTime.tms_stime - startTime->tms_stime) +
272 (endTime.tms_utime - startTime->tms_ 272 (endTime.tms_utime - startTime->tms_utime)) /
273 sysconf(_SC_CLK_TCK); 273 sysconf(_SC_CLK_TCK);
274 } 274 }
275 275
276 //-------------------------------------------- 276 //----------------------------------------------------------------------------//
277 277
278 G4StatAnalysis& G4StatAnalysis::operator+=(con 278 G4StatAnalysis& G4StatAnalysis::operator+=(const G4double& _val)
279 { 279 {
280 this->Add(_val); 280 this->Add(_val);
281 return *this; 281 return *this;
282 } 282 }
283 283
284 //-------------------------------------------- 284 //----------------------------------------------------------------------------//
285 285
286 G4StatAnalysis& G4StatAnalysis::operator/=(con 286 G4StatAnalysis& G4StatAnalysis::operator/=(const G4double& _val)
287 { 287 {
288 fSum1 /= _val; 288 fSum1 /= _val;
289 fSum2 /= (_val * _val); 289 fSum2 /= (_val * _val);
290 return *this; 290 return *this;
291 } 291 }
292 292
293 //-------------------------------------------- 293 //----------------------------------------------------------------------------//
294 294
295 G4StatAnalysis& G4StatAnalysis::operator+=(con 295 G4StatAnalysis& G4StatAnalysis::operator+=(const G4StatAnalysis& rhs)
296 { 296 {
297 fHits += rhs.fHits; 297 fHits += rhs.fHits;
298 fSum1 += rhs.fSum1; 298 fSum1 += rhs.fSum1;
299 fSum2 += rhs.fSum2; 299 fSum2 += rhs.fSum2;
300 fZero += rhs.fZero; 300 fZero += rhs.fZero;
301 return *this; 301 return *this;
302 } 302 }
303 303
304 //-------------------------------------------- 304 //----------------------------------------------------------------------------//
305 305
306 G4StatAnalysis& G4StatAnalysis::operator-=(con 306 G4StatAnalysis& G4StatAnalysis::operator-=(const G4StatAnalysis& rhs)
307 { 307 {
308 fHits -= rhs.fHits; 308 fHits -= rhs.fHits;
309 fSum1 -= rhs.fSum1; 309 fSum1 -= rhs.fSum1;
310 fSum2 -= rhs.fSum2; 310 fSum2 -= rhs.fSum2;
311 fZero -= rhs.fZero; 311 fZero -= rhs.fZero;
312 return *this; 312 return *this;
313 } 313 }
314 314
315 //-------------------------------------------- 315 //----------------------------------------------------------------------------//
316 316
317 inline G4Allocator<G4StatAnalysis>*& _aStatAna 317 inline G4Allocator<G4StatAnalysis>*& _aStatAnalysisAllocator_G4MT_TLS_()
318 { 318 {
319 G4ThreadLocalStatic G4Allocator<G4StatAnalys 319 G4ThreadLocalStatic G4Allocator<G4StatAnalysis>* _instance =
320 new G4Allocator<G4StatAnalysis>(); 320 new G4Allocator<G4StatAnalysis>();
321 return _instance; 321 return _instance;
322 } 322 }
323 323
324 //-------------------------------------------- 324 //----------------------------------------------------------------------------//
325 325
326 void* G4StatAnalysis::operator new(std::size_t 326 void* G4StatAnalysis::operator new(std::size_t)
327 { 327 {
328 G4Allocator<G4StatAnalysis>& _allocator = 328 G4Allocator<G4StatAnalysis>& _allocator =
329 *_aStatAnalysisAllocator_G4MT_TLS_(); 329 *_aStatAnalysisAllocator_G4MT_TLS_();
330 return (void*) _allocator.MallocSingle(); 330 return (void*) _allocator.MallocSingle();
331 } 331 }
332 332
333 //-------------------------------------------- 333 //----------------------------------------------------------------------------//
334 334
335 void G4StatAnalysis::operator delete(void* _pt 335 void G4StatAnalysis::operator delete(void* _ptr)
336 { 336 {
337 G4Allocator<G4StatAnalysis>& _allocator = 337 G4Allocator<G4StatAnalysis>& _allocator =
338 *_aStatAnalysisAllocator_G4MT_TLS_(); 338 *_aStatAnalysisAllocator_G4MT_TLS_();
339 _allocator.FreeSingle((G4StatAnalysis*) _ptr 339 _allocator.FreeSingle((G4StatAnalysis*) _ptr);
340 } 340 }
341 341
342 //-------------------------------------------- 342 //----------------------------------------------------------------------------//
343 343