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