| Geant4 Cross Reference |
1 // 1
2 // *******************************************
3 // * License and Disclaimer
4 // *
5 // * The Geant4 software is copyright of th
6 // * the Geant4 Collaboration. It is provided
7 // * conditions of the Geant4 Software License
8 // * LICENSE and available at http://cern.ch/
9 // * include a list of copyright holders.
10 // *
11 // * Neither the authors of this software syst
12 // * institutes,nor the agencies providing fin
13 // * work make any representation or warran
14 // * regarding this software system or assum
15 // * use. Please see the license in the file
16 // * for the full disclaimer and the limitatio
17 // *
18 // * This code implementation is the result
19 // * technical work of the GEANT4 collaboratio
20 // * By using, copying, modifying or distri
21 // * any work based on the software) you ag
22 // * use in resulting scientific publicati
23 // * acceptance of all terms of the Geant4 Sof
24 // *******************************************
25 //
26 // G4StatDouble class implementation
27 //
28 // Original Author: Giovanni Santin (ESA) - Oc
29 // Adapted by: John Apostolakis - November 201
30 // -------------------------------------------
31 #include "G4StatDouble.hh"
32
33 G4StatDouble::G4StatDouble() { reset(); }
34
35 G4StatDouble::G4StatDouble(G4double x) { fill(
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, G4doub
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: w
61 }
62 }
63
64 void G4StatDouble::scale(G4double value) { m_s
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
77 {
78 G4double factor = 0.;
79 // factor to rescale the Mean for the reques
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, G
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 par
100 // when w_i = w
101 // ((m_sum_w * m_sum_w) / (m_sum_w * m_s
102
103 // from NIST "DATAPLOT Reference manual"
104 // http://www.itl.nist.gov/div898/softwa
105 // rewritten based on: SUM[w(x-m)^2]/SUM
106 // and dividing it by sqrt[n] to go from
107 // rms of the mean value
108
109 (xn / (xn - 1)) * ((ssum_wx2 / ssum_w) -
110
111 tmp = std::max(tmp, 0.0); // this avoids
112 vrms = std::sqrt(tmp);
113 // G4cout << "[G4StatDoubleElement::rms]
114 // << " m_sum_wx2: " << m_sum_wx2
115 // << " m_n: " << m_n << " tmp:
116 // << G4endl;
117 // G4cout << "[G4StatDoubleElement::rms]
118 // 1))
119 // << " (m_sum_wx2 / m_sum_w): " << (
120 // << " (mean * mean): " << (mean * m
121 // << " ((m_sum_wx2 / m_sum_w) - (mea
122 // << ((m_sum_wx2 / m_sum_w) - (me
123 // << G4endl;
124 }
125 return vrms * m_scale;
126 }
127
128 G4double G4StatDouble::rms()
129 {
130 // this method computes the RMS with "all in
131 // all the sums are the internal ones: m_sum
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,
137 {
138 // this method computes the RMS with sum_w a
139 // ext_sum_w and ext_n:
140 // this means that the result is normalised
141 // it is useful when, given a number ext_n o
142 // ext_sum_w, only m_n (with sum of weights
143 // in the internal summation (e.g. for a dos
144 // only a few particles reach that volume)
145
146 return rms(m_sum_wx, m_sum_wx2, ext_sum_w, e
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