00001 // 00002 // ******************************************************************** 00003 // * License and Disclaimer * 00004 // * * 00005 // * The Geant4 software is copyright of the Copyright Holders of * 00006 // * the Geant4 Collaboration. It is provided under the terms and * 00007 // * conditions of the Geant4 Software License, included in the file * 00008 // * LICENSE and available at http://cern.ch/geant4/license . These * 00009 // * include a list of copyright holders. * 00010 // * * 00011 // * Neither the authors of this software system, nor their employing * 00012 // * institutes,nor the agencies providing financial support for this * 00013 // * work make any representation or warranty, express or implied, * 00014 // * regarding this software system or assume any liability for its * 00015 // * use. Please see the license in the file LICENSE and URL above * 00016 // * for the full disclaimer and the limitation of liability. * 00017 // * * 00018 // * This code implementation is the result of the scientific and * 00019 // * technical work of the GEANT4 collaboration. * 00020 // * By using, copying, modifying or distributing the software (or * 00021 // * any work based on the software) you agree to acknowledge its * 00022 // * use in resulting scientific publications, and indicate your * 00023 // * acceptance of all terms of the Geant4 Software license. * 00024 // ******************************************************************** 00025 // 00026 // 00027 // $Id:$ 00028 // GEANT4 tag $Name: not supported by cvs2svn $ 00029 // 00030 // 00031 // ---------------------------------------------------------------------- 00032 // class G4StatDouble 00033 // 00034 // Implementation. 00035 // Original Author: Giovanni Santin (ESA) - October 2005 in GRAS tool 00036 // Adapted by: John Apostolakis - November 2011 00037 00038 #include "G4StatDouble.hh" 00039 00040 G4StatDouble::G4StatDouble() 00041 { 00042 reset(); 00043 } 00044 00045 void G4StatDouble::reset() 00046 { 00047 m_sum_wx = 0.; 00048 m_sum_wx2 = 0.; 00049 m_n = 0; 00050 m_sum_w = 0.; 00051 m_sum_w2 = 0.; 00052 m_scale = 1.; 00053 } 00054 00055 G4StatDouble::~G4StatDouble() 00056 {} 00057 00058 void G4StatDouble::fill(G4double value, G4double weight) 00059 { 00060 m_sum_wx += value * weight; 00061 m_sum_wx2 += value * value * weight; 00062 m_n++; 00063 m_sum_w += weight; 00064 m_sum_w2 += weight * weight; 00065 00066 if (weight <= 0.) 00067 { 00068 G4cout << "[G4StatDouble::fill] WARNING: weight<=0. " 00069 << weight << G4endl; 00070 } 00071 } 00072 00073 void G4StatDouble::scale(G4double value) 00074 { 00075 m_scale = m_scale * value; 00076 } 00077 00078 G4double G4StatDouble::mean() const 00079 { 00080 G4double mean_val = 0.; 00081 if (m_sum_w > 0.) 00082 { 00083 mean_val = m_sum_wx / m_sum_w; 00084 } 00085 return m_scale * mean_val; 00086 } 00087 00088 G4double G4StatDouble::mean(G4double ext_sum_w) const 00089 { 00090 G4double factor = 0.; 00091 // factor to rescale the Mean for the requested number 00092 // of events (or sum of weights) ext_sum_w 00093 00094 if (ext_sum_w > 0) 00095 { 00096 factor = m_sum_w; 00097 factor /= ext_sum_w; 00098 } 00099 return mean() * factor; 00100 00101 } 00102 00103 G4double G4StatDouble::rms(G4double ssum_wx, G4double ssum_wx2, 00104 G4double ssum_w, G4int nn) 00105 { 00106 G4double vrms; 00107 if (nn > 1) 00108 { 00109 G4double vmean = ssum_wx / ssum_w; 00110 G4double xn = nn; 00111 G4double tmp = 00112 // from GNU Scientific Library. This part is equivalent to N/(N-1) 00113 // when w_i = w 00114 // ((m_sum_w * m_sum_w) / (m_sum_w * m_sum_w - m_sum_w2)) 00115 00116 // from NIST "DATAPLOT Reference manual", Page 2-66 00117 // http://www.itl.nist.gov/div898/software/dataplot/refman2/ch2/weightsd.pdf 00118 // rewritten based on: SUM[w(x-m)^2]/SUM[w] = SUM[wx^2]/SUM[w] - m^2 00119 // and dividing it by sqrt[n] to go from rms of distribution to the 00120 // rms of the mean value 00121 00122 (1. / (xn - 1)) 00123 * ((ssum_wx2 / ssum_w) - (vmean * vmean)); 00124 00125 if (tmp < 0.) tmp=0.; // this avoids observed computation problem 00126 vrms = std::sqrt( tmp ); 00127 // G4cout << "[G4StatDoubleElement::rms] m_sum_wx: " << m_sum_wx 00128 // << " m_sum_wx2: " << m_sum_wx2 << " m_sum_w: " << m_sum_w 00129 // << " m_n: " << m_n << " tmp: " << tmp<< " rms: " << rms 00130 // << G4endl; 00131 // G4cout << "[G4StatDoubleElement::rms] (m_n / (m_n - 1)): " << (xn/(xn - 1)) 00132 // << " (m_sum_wx2 / m_sum_w): " << (m_sum_wx2 / m_sum_w) 00133 // << " (mean * mean): " << (mean * mean) 00134 // << " ((m_sum_wx2 / m_sum_w) - (mean * mean)): " 00135 // << ((m_sum_wx2 / m_sum_w) - (mean * mean)) 00136 // << G4endl; 00137 } 00138 else 00139 { 00140 vrms = -1.; 00141 } 00142 return vrms * m_scale; 00143 } 00144 00145 G4double G4StatDouble::rms() 00146 { 00147 // this method computes the RMS with "all internal" parameters: 00148 // all the sums are the internal ones: m_sum_wx, m_sum_wx2, m_sum_w, m_n 00149 00150 return rms(m_sum_wx, m_sum_wx2, m_sum_w, m_n); 00151 } 00152 00153 G4double G4StatDouble::rms(G4double ext_sum_w, G4int ext_n) 00154 { 00155 // this method computes the RMS with sum_w and n coming from outside: 00156 // ext_sum_w and ext_n: 00157 // this means that the result is normalised to the external events 00158 // it is useful when, given a number ext_n of events with sum of the weights 00159 // ext_sum_w, only m_n (with sum of weights m_sum_w) are actually accumulated 00160 // in the internal summation (e.g. for a dose variable in a volume, because 00161 // only a few particles reach that volume) 00162 00163 return rms(m_sum_wx, m_sum_wx2, ext_sum_w, ext_n); 00164 }