Geant4  10.03.p03
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
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  reset();
47  fill(x);
48 }
49 
51 {
52  m_sum_wx = 0.;
53  m_sum_wx2 = 0.;
54  m_n = 0;
55  m_sum_w = 0.;
56  m_sum_w2 = 0.;
57  m_scale = 1.;
58 }
59 
61 {}
62 
64 {
65  m_sum_wx += value * weight;
66  m_sum_wx2 += value * value * weight;
67  if(m_n < INT_MAX) { ++m_n; }
68  m_sum_w += weight;
69  m_sum_w2 += weight * weight;
70 
71  if (weight <= 0.)
72  {
73  G4cout << "[G4StatDouble::fill] WARNING: weight<=0. "
74  << weight << G4endl;
75  }
76 }
77 
79 {
80  m_scale = m_scale * value;
81 }
82 
84 {
85  G4double mean_val = 0.;
86  if (m_sum_w > 0.)
87  {
88  mean_val = m_sum_wx / m_sum_w;
89  }
90  return m_scale * mean_val;
91 }
92 
94 {
95  G4double factor = 0.;
96  // factor to rescale the Mean for the requested number
97  // of events (or sum of weights) ext_sum_w
98 
99  if (ext_sum_w > 0)
100  {
101  factor = m_sum_w;
102  factor /= ext_sum_w;
103  }
104  return mean() * factor;
105 
106 }
107 
109  G4double ssum_w, G4int nn)
110 {
111  G4double vrms = 0.0;
112  if (nn > 1)
113  {
114  G4double vmean = ssum_wx / ssum_w;
115  G4double xn = nn;
116  G4double tmp =
117  // from GNU Scientific Library. This part is equivalent to N/(N-1)
118  // when w_i = w
119  // ((m_sum_w * m_sum_w) / (m_sum_w * m_sum_w - m_sum_w2))
120 
121  // from NIST "DATAPLOT Reference manual", Page 2-66
122  // http://www.itl.nist.gov/div898/software/dataplot/refman2/ch2/weightsd.pdf
123  // rewritten based on: SUM[w(x-m)^2]/SUM[w] = SUM[wx^2]/SUM[w] - m^2
124  // and dividing it by sqrt[n] to go from rms of distribution to the
125  // rms of the mean value
126 
127  (xn / (xn - 1))
128  * ((ssum_wx2 / ssum_w) - (vmean * vmean));
129 
130  tmp = std::max(tmp, 0.0); // this avoids observed computation problem
131  vrms = std::sqrt( tmp );
132 // G4cout << "[G4StatDoubleElement::rms] m_sum_wx: " << m_sum_wx
133 // << " m_sum_wx2: " << m_sum_wx2 << " m_sum_w: " << m_sum_w
134 // << " m_n: " << m_n << " tmp: " << tmp<< " rms: " << rms
135 // << G4endl;
136 // G4cout << "[G4StatDoubleElement::rms] (m_n / (m_n - 1)): " << (xn/(xn - 1))
137 // << " (m_sum_wx2 / m_sum_w): " << (m_sum_wx2 / m_sum_w)
138 // << " (mean * mean): " << (mean * mean)
139 // << " ((m_sum_wx2 / m_sum_w) - (mean * mean)): "
140 // << ((m_sum_wx2 / m_sum_w) - (mean * mean))
141 // << G4endl;
142  }
143  return vrms * m_scale;
144 }
145 
147 {
148  // this method computes the RMS with "all internal" parameters:
149  // all the sums are the internal ones: m_sum_wx, m_sum_wx2, m_sum_w, m_n
150 
151  return rms(m_sum_wx, m_sum_wx2, m_sum_w, m_n);
152 }
153 
155 {
156  // this method computes the RMS with sum_w and n coming from outside:
157  // ext_sum_w and ext_n:
158  // this means that the result is normalised to the external events
159  // it is useful when, given a number ext_n of events with sum of the weights
160  // ext_sum_w, only m_n (with sum of weights m_sum_w) are actually accumulated
161  // in the internal summation (e.g. for a dose variable in a volume, because
162  // only a few particles reach that volume)
163 
164  return rms(m_sum_wx, m_sum_wx2, ext_sum_w, ext_n);
165 }
166 
168 {
169  m_n += ptr->n();
170  m_sum_w += ptr->sum_w();
171  m_sum_w2 += ptr->sum_w2();
172  m_sum_wx += ptr->sum_wx();
173  m_sum_wx2 += ptr->sum_wx2();
174 }
G4double sum_wx2() const
Definition: G4StatDouble.hh:99
G4double sum_w2() const
Definition: G4StatDouble.hh:97
void fill(G4double x, G4double weight=1.)
Definition: G4StatDouble.cc:63
void add(const G4StatDouble *)
tuple x
Definition: test.py:50
int G4int
Definition: G4Types.hh:78
G4double m_sum_wx
G4double sum_w() const
Definition: G4StatDouble.hh:96
G4double m_sum_wx2
virtual ~G4StatDouble()
Definition: G4StatDouble.cc:60
G4double m_sum_w2
G4GLOB_DLL std::ostream G4cout
const XML_Char int const XML_Char * value
Definition: expat.h:331
G4double sum_wx() const
Definition: G4StatDouble.hh:98
#define INT_MAX
Definition: templates.hh:111
G4double m_sum_w
T max(const T t1, const T t2)
brief Return the largest of the two arguments
void scale(G4double)
Definition: G4StatDouble.cc:78
G4double m_scale
G4double mean() const
Definition: G4StatDouble.cc:83
#define G4endl
Definition: G4ios.hh:61
G4double rms()
double G4double
Definition: G4Types.hh:76
G4int n() const
Definition: G4StatDouble.hh:95