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