| Geant4 Cross Reference |
1 //
2 // ********************************************************************
3 // * License and Disclaimer *
4 // * *
5 // * The Geant4 software is copyright of the Copyright Holders of *
6 // * the Geant4 Collaboration. It is provided under the terms and *
7 // * conditions of the Geant4 Software License, included in the file *
8 // * LICENSE and available at http://cern.ch/geant4/license . These *
9 // * include a list of copyright holders. *
10 // * *
11 // * Neither the authors of this software system, nor their employing *
12 // * institutes,nor the agencies providing financial support for this *
13 // * work make any representation or warranty, express or implied, *
14 // * regarding this software system or assume any liability for its *
15 // * use. Please see the license in the file LICENSE and URL above *
16 // * for the full disclaimer and the limitation of liability. *
17 // * *
18 // * This code implementation is the result of the scientific and *
19 // * technical work of the GEANT4 collaboration. *
20 // * By using, copying, modifying or distributing the software (or *
21 // * any work based on the software) you agree to acknowledge its *
22 // * use in resulting scientific publications, and indicate your *
23 // * acceptance of all terms of the Geant4 Software license. *
24 // ********************************************************************
25 //
26 // G4StatAnalysis inline methods implementation
27 //
28 // Author: J.Madsen, 25.10.2018
29 // --------------------------------------------------------------------
30
31 #include <cmath>
32 #include <fstream>
33 #include <iomanip>
34 #include <iostream>
35 #include <limits>
36
37 #include "G4Timer.hh"
38 #include "G4ios.hh"
39 #include "globals.hh"
40 #include "tls.hh"
41
42 #include "G4Allocator.hh"
43 #include "G4Types.hh"
44
45 //----------------------------------------------------------------------------//
46
47 G4StatAnalysis::G4StatAnalysis() {}
48
49 //----------------------------------------------------------------------------//
50
51 G4double G4StatAnalysis::GetMean() const
52 {
53 return (fHits > 0) ? fSum1 / ((G4double) fHits) : 0.;
54 }
55
56 //----------------------------------------------------------------------------//
57
58 const G4double& G4StatAnalysis::GetSum() const { return fSum1; }
59
60 //----------------------------------------------------------------------------//
61
62 const G4double& G4StatAnalysis::GetSumSquared() const { return fSum2; }
63
64 //----------------------------------------------------------------------------//
65
66 const G4double& G4StatAnalysis::GetSum1() const { return fSum1; }
67
68 //----------------------------------------------------------------------------//
69
70 const G4double& G4StatAnalysis::GetSum2() const { return fSum2; }
71
72 //----------------------------------------------------------------------------//
73
74 const G4int& G4StatAnalysis::GetHits() const { return fHits; }
75
76 //----------------------------------------------------------------------------//
77
78 G4int G4StatAnalysis::GetNumNonZero() const { return fHits - fZero; }
79
80 //----------------------------------------------------------------------------//
81
82 G4int G4StatAnalysis::GetNumZero() const { return fZero; }
83
84 //----------------------------------------------------------------------------//
85
86 void G4StatAnalysis::SetSum(const G4double& val) { fSum1 = val; }
87
88 //----------------------------------------------------------------------------//
89
90 void G4StatAnalysis::SetSumSquared(const G4double& val) { fSum2 = val; }
91
92 //----------------------------------------------------------------------------//
93
94 void G4StatAnalysis::SetSum1(const G4double& val) { fSum1 = val; }
95
96 //----------------------------------------------------------------------------//
97
98 void G4StatAnalysis::SetSum2(const G4double& val) { fSum2 = val; }
99
100 //----------------------------------------------------------------------------//
101
102 void G4StatAnalysis::SetHits(const G4int& val) { fHits = val; }
103
104 //----------------------------------------------------------------------------//
105
106 void G4StatAnalysis::SetZero(const G4int& val) { fZero = val; }
107
108 //----------------------------------------------------------------------------//
109
110 G4double G4StatAnalysis::GetFOM() const
111 {
112 G4double elapsed_time = this->GetCpuTime();
113 G4double relative_err = this->GetRelativeError();
114 // lambda for equation clarity (will be inlined)
115 auto compute_figure_of_merit = [&]() {
116 return (1.0 / (relative_err * relative_err) / elapsed_time);
117 };
118 return (std::fabs(relative_err) > 0.0 && elapsed_time > 0.0)
119 ? compute_figure_of_merit()
120 : ((fHits > 0) ? 1.0 : 0.0);
121 }
122
123 //----------------------------------------------------------------------------//
124
125 G4double G4StatAnalysis::GetRelativeError() const
126 {
127 // lambda for equation clarity (will be inlined)
128 auto compute_relative_error = [&]() {
129 return (GetStdDev() / GetMean() / std::sqrt((G4double) fHits));
130 };
131 return (std::fabs(GetMean()) > 0 && fHits > 0) ? compute_relative_error()
132 : ((fHits > 0) ? 1.0 : 0.0);
133 }
134
135 //----------------------------------------------------------------------------//
136
137 G4double G4StatAnalysis::GetStdDev() const
138 {
139 return ::sqrt(std::fabs(GetVariance()));
140 }
141
142 //----------------------------------------------------------------------------//
143
144 G4double G4StatAnalysis::GetVariance() const
145 {
146 // lambda for equation clarity (will be inlined)
147 auto compute_variance = [&]() {
148 return ((fSum2 - (std::pow(fSum1, 2.0) / fHits)) /
149 (((G4double) fHits) - 1.0));
150 };
151 return (fHits > 1) ? compute_variance() : 0.0;
152 }
153
154 //----------------------------------------------------------------------------//
155
156 G4double G4StatAnalysis::GetCoeffVariation() const
157 {
158 // lambda for equation clarity (will be inlined)
159 auto coefficient_of_variation = [&]() {
160 G4double hits = fHits;
161 return ::sqrt((hits / (hits - 1.0)) *
162 ((fSum2 / (fSum1 * fSum1)) - (1.0 / hits)));
163 };
164 return (fHits > 1) ? coefficient_of_variation() : 0.0;
165 // return (fHits > 0 && fabs(fSum1) > 0.0)
166 // ? (100.0*GetStdDev()/GetMean()) : 0.0;
167 }
168
169 //----------------------------------------------------------------------------//
170
171 G4double G4StatAnalysis::GetEfficiency() const
172 {
173 G4double hits = fHits;
174 G4double nzero = fHits - fZero;
175 return (fHits > 0) ? (nzero / hits) : 0.0;
176 }
177
178 //----------------------------------------------------------------------------//
179
180 G4double G4StatAnalysis::GetR2Int() const
181 {
182 G4double hits = fHits;
183 return (fHits > 0)
184 ? (fSum2 / (fSum1 * fSum1)) - 1.0 / (GetEfficiency() * hits)
185 : 0.0;
186 }
187
188 //----------------------------------------------------------------------------//
189
190 G4double G4StatAnalysis::GetR2Eff() const
191 {
192 G4double hits = fHits;
193 return (fHits > 0) ? (1.0 - GetEfficiency()) / (GetEfficiency() * hits) : 0.0;
194 }
195
196 //----------------------------------------------------------------------------//
197
198 G4StatAnalysis::operator G4double() const { return this->GetSum(); }
199
200 //----------------------------------------------------------------------------//
201
202 void G4StatAnalysis::Reset()
203 {
204 fHits = 0;
205 fZero = 0;
206 fSum1 = 0.0;
207 fSum2 = 0.0;
208 }
209
210 //----------------------------------------------------------------------------//
211
212 void G4StatAnalysis::Add(const G4double& val, const G4double& weight)
213 {
214 fHits += 1;
215 fSum1 += val * weight;
216 fSum2 += val * val * weight;
217 if(std::fabs(val * weight) <
218 std::fabs(GetMean() * std::numeric_limits<double>::epsilon()))
219 fZero += 1;
220 }
221
222 //----------------------------------------------------------------------------//
223
224 void G4StatAnalysis::Rescale(const G4double& factor)
225 {
226 fSum1 *= factor;
227 fSum2 *= factor * factor;
228 }
229
230 //----------------------------------------------------------------------------//
231
232 void G4StatAnalysis::PrintInfo(std::ostream& os, const std::string& tab) const
233 {
234 G4int _hits = this->GetHits();
235 G4double _sum = this->GetSum();
236 G4double _sigma = this->GetStdDev();
237 G4double _coeff = this->GetCoeffVariation();
238 G4double _error = this->GetRelativeError();
239 G4double _eff = this->GetEfficiency();
240 G4double _fom = this->GetFOM();
241 G4double _r2int = this->GetR2Int();
242 G4double _r2eff = this->GetR2Eff();
243
244 using std::fixed;
245 using std::ios;
246 using std::left;
247 using std::right;
248 using std::scientific;
249 using std::setprecision;
250 using std::setw;
251
252 std::stringstream ss;
253 ss << tab //<< scientific
254 << setprecision((G4int)os.precision()) << right << _sum << left
255 << " [sigma: " << right << _sigma << left << " | error: " << right
256 << _error << left << " | coeff: " << right << _coeff << left
257 << " | eff: " << right << _eff << left << " | fom: " << right << _fom
258 << left << " | r2int: " << right << _r2int << left << " | r2eff: " << right
259 << _r2eff << left << " | hits: " << right << _hits << left << " ]";
260
261 os << ss.str();
262 }
263
264 //----------------------------------------------------------------------------//
265
266 G4double G4StatAnalysis::GetCpuTime() const
267 {
268 tms* startTime = GetCpuClock();
269 tms endTime;
270 times(&endTime);
271 return ((endTime.tms_stime - startTime->tms_stime) +
272 (endTime.tms_utime - startTime->tms_utime)) /
273 sysconf(_SC_CLK_TCK);
274 }
275
276 //----------------------------------------------------------------------------//
277
278 G4StatAnalysis& G4StatAnalysis::operator+=(const G4double& _val)
279 {
280 this->Add(_val);
281 return *this;
282 }
283
284 //----------------------------------------------------------------------------//
285
286 G4StatAnalysis& G4StatAnalysis::operator/=(const G4double& _val)
287 {
288 fSum1 /= _val;
289 fSum2 /= (_val * _val);
290 return *this;
291 }
292
293 //----------------------------------------------------------------------------//
294
295 G4StatAnalysis& G4StatAnalysis::operator+=(const G4StatAnalysis& rhs)
296 {
297 fHits += rhs.fHits;
298 fSum1 += rhs.fSum1;
299 fSum2 += rhs.fSum2;
300 fZero += rhs.fZero;
301 return *this;
302 }
303
304 //----------------------------------------------------------------------------//
305
306 G4StatAnalysis& G4StatAnalysis::operator-=(const G4StatAnalysis& rhs)
307 {
308 fHits -= rhs.fHits;
309 fSum1 -= rhs.fSum1;
310 fSum2 -= rhs.fSum2;
311 fZero -= rhs.fZero;
312 return *this;
313 }
314
315 //----------------------------------------------------------------------------//
316
317 inline G4Allocator<G4StatAnalysis>*& _aStatAnalysisAllocator_G4MT_TLS_()
318 {
319 G4ThreadLocalStatic G4Allocator<G4StatAnalysis>* _instance =
320 new G4Allocator<G4StatAnalysis>();
321 return _instance;
322 }
323
324 //----------------------------------------------------------------------------//
325
326 void* G4StatAnalysis::operator new(std::size_t)
327 {
328 G4Allocator<G4StatAnalysis>& _allocator =
329 *_aStatAnalysisAllocator_G4MT_TLS_();
330 return (void*) _allocator.MallocSingle();
331 }
332
333 //----------------------------------------------------------------------------//
334
335 void G4StatAnalysis::operator delete(void* _ptr)
336 {
337 G4Allocator<G4StatAnalysis>& _allocator =
338 *_aStatAnalysisAllocator_G4MT_TLS_();
339 _allocator.FreeSingle((G4StatAnalysis*) _ptr);
340 }
341
342 //----------------------------------------------------------------------------//
343