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