| 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 // G4StatDouble class implementation
27 //
28 // Original Author: Giovanni Santin (ESA) - October 2005 in GRAS tool
29 // Adapted by: John Apostolakis - November 2011
30 // --------------------------------------------------------------------
31 #include "G4StatDouble.hh"
32
33 G4StatDouble::G4StatDouble() { reset(); }
34
35 G4StatDouble::G4StatDouble(G4double x) { fill(x); }
36
37 void G4StatDouble::reset()
38 {
39 m_sum_wx = 0.;
40 m_sum_wx2 = 0.;
41 m_n = 0;
42 m_sum_w = 0.;
43 m_sum_w2 = 0.;
44 m_scale = 1.;
45 }
46
47 void G4StatDouble::fill(G4double value, G4double weight)
48 {
49 m_sum_wx += value * weight;
50 m_sum_wx2 += value * value * weight;
51 if(m_n < INT_MAX)
52 {
53 ++m_n;
54 }
55 m_sum_w += weight;
56 m_sum_w2 += weight * weight;
57
58 if(weight <= 0.)
59 {
60 G4cout << "[G4StatDouble::fill] WARNING: weight<=0. " << weight << G4endl;
61 }
62 }
63
64 void G4StatDouble::scale(G4double value) { m_scale = m_scale * value; }
65
66 G4double G4StatDouble::mean() const
67 {
68 G4double mean_val = 0.;
69 if(m_sum_w > 0.)
70 {
71 mean_val = m_sum_wx / m_sum_w;
72 }
73 return m_scale * mean_val;
74 }
75
76 G4double G4StatDouble::mean(G4double ext_sum_w) const
77 {
78 G4double factor = 0.;
79 // factor to rescale the Mean for the requested number
80 // of events (or sum of weights) ext_sum_w
81
82 if(ext_sum_w > 0)
83 {
84 factor = m_sum_w;
85 factor /= ext_sum_w;
86 }
87 return mean() * factor;
88 }
89
90 G4double G4StatDouble::rms(G4double ssum_wx, G4double ssum_wx2, G4double ssum_w,
91 G4int nn)
92 {
93 G4double vrms = 0.0;
94 if(nn > 1)
95 {
96 G4double vmean = ssum_wx / ssum_w;
97 G4double xn = nn;
98 G4double tmp =
99 // from GNU Scientific Library. This part is equivalent to N/(N-1)
100 // when w_i = w
101 // ((m_sum_w * m_sum_w) / (m_sum_w * m_sum_w - m_sum_w2))
102
103 // from NIST "DATAPLOT Reference manual", Page 2-66
104 // http://www.itl.nist.gov/div898/software/dataplot/refman2/ch2/weightsd.pdf
105 // rewritten based on: SUM[w(x-m)^2]/SUM[w] = SUM[wx^2]/SUM[w] - m^2
106 // and dividing it by sqrt[n] to go from rms of distribution to the
107 // rms of the mean value
108
109 (xn / (xn - 1)) * ((ssum_wx2 / ssum_w) - (vmean * vmean));
110
111 tmp = std::max(tmp, 0.0); // this avoids observed computation problem
112 vrms = std::sqrt(tmp);
113 // G4cout << "[G4StatDoubleElement::rms] m_sum_wx: " << m_sum_wx
114 // << " m_sum_wx2: " << m_sum_wx2 << " m_sum_w: " << m_sum_w
115 // << " m_n: " << m_n << " tmp: " << tmp<< " rms: " << rms
116 // << G4endl;
117 // G4cout << "[G4StatDoubleElement::rms] (m_n / (m_n - 1)): " << (xn/(xn -
118 // 1))
119 // << " (m_sum_wx2 / m_sum_w): " << (m_sum_wx2 / m_sum_w)
120 // << " (mean * mean): " << (mean * mean)
121 // << " ((m_sum_wx2 / m_sum_w) - (mean * mean)): "
122 // << ((m_sum_wx2 / m_sum_w) - (mean * mean))
123 // << G4endl;
124 }
125 return vrms * m_scale;
126 }
127
128 G4double G4StatDouble::rms()
129 {
130 // this method computes the RMS with "all internal" parameters:
131 // all the sums are the internal ones: m_sum_wx, m_sum_wx2, m_sum_w, m_n
132
133 return rms(m_sum_wx, m_sum_wx2, m_sum_w, m_n);
134 }
135
136 G4double G4StatDouble::rms(G4double ext_sum_w, G4int ext_n)
137 {
138 // this method computes the RMS with sum_w and n coming from outside:
139 // ext_sum_w and ext_n:
140 // this means that the result is normalised to the external events
141 // it is useful when, given a number ext_n of events with sum of the weights
142 // ext_sum_w, only m_n (with sum of weights m_sum_w) are actually accumulated
143 // in the internal summation (e.g. for a dose variable in a volume, because
144 // only a few particles reach that volume)
145
146 return rms(m_sum_wx, m_sum_wx2, ext_sum_w, ext_n);
147 }
148
149 void G4StatDouble::add(const G4StatDouble* ptr)
150 {
151 m_n += ptr->n();
152 m_sum_w += ptr->sum_w();
153 m_sum_w2 += ptr->sum_w2();
154 m_sum_wx += ptr->sum_wx();
155 m_sum_wx2 += ptr->sum_wx2();
156 }
157