Gaudi Framework, version v20r3
Generated: 24 Nov 2008
|
|
|
Generated: 24 Nov 2008 |
00001 // $Id: HistoStats.cpp,v 1.5 2008/06/04 08:18:22 marcocle Exp $ 00002 // ============================================================================ 00003 // Include files 00004 // ============================================================================ 00005 // STD & STL 00006 // ============================================================================ 00007 #include <cmath> 00008 #include <limits> 00009 #include <climits> 00010 // ============================================================================ 00011 // AIDA 00012 // ============================================================================ 00013 #include "AIDA/IHistogram1D.h" 00014 #include "AIDA/IAxis.h" 00015 // ============================================================================ 00016 // GaudiUtils 00017 // ============================================================================ 00018 #include "GaudiUtils/HistoStats.h" 00019 // ============================================================================ 00025 // ============================================================================ 00026 namespace 00027 { 00028 // ========================================================================== 00030 const double s_bad = -1 * std::numeric_limits<double>::max() ; 00031 // ========================================================================== 00032 } 00033 // ============================================================================ 00034 /* get the moment of the certain around the specified "value" 00035 * @param histo histogram 00036 * @param order the momentm order 00037 * @param value central value 00038 * @return the evaluated moment 00039 */ 00040 // ============================================================================ 00041 double Gaudi::Utils::HistoStats::moment 00042 ( const AIDA::IHistogram1D* histo , 00043 const unsigned int order , 00044 const double value ) 00045 { 00046 if ( 0 == histo ) { return s_bad ; } // RETURN 00047 if ( 0 == order ) { return 1.0 ; } // RETURN 00048 if ( 1 == order ) { return mean( histo ) - value ; } // RETURN 00049 if ( 2 == order ) 00050 { 00051 const double _r = rms ( histo ) ; 00052 const double _d = value - mean ( histo ) ; 00053 return _r *_r + _d * _d ; // RETURN 00054 } 00055 const double n = nEff ( histo ) ; 00056 if ( 0 >= n ) { return 0.0 ; } // RETURN 00057 00058 const int _order = (int) order ; 00059 00060 // get the exis 00061 const AIDA::IAxis& axis = histo -> axis () ; 00062 // number of bins 00063 const int nBins = axis.bins() ; 00064 double result = 0 ; 00065 double weight = 0 ; 00066 // loop over all bins 00067 for ( int i = 0 ; i < nBins ; ++i ) 00068 { 00069 const double lE = axis . binLowerEdge ( i ) ; 00070 const double uE = axis . binUpperEdge ( i ) ; 00071 // 00072 const double yBin = histo -> binHeight ( i ) ; // bin weight 00073 const double xBin = 0.5 * ( lE + uE ) ; // bin center 00074 // 00075 weight += yBin ; 00076 result += yBin * std::pow ( xBin - value , _order ) ; 00077 } 00078 if ( 0 != weight ) { result /= weight ; } 00079 return result ; 00080 } 00081 // ====================================================================== 00088 // ====================================================================== 00089 double Gaudi::Utils::HistoStats::momentErr 00090 ( const AIDA::IHistogram1D* histo , 00091 const unsigned int order ) 00092 { 00093 if ( 0 == histo ) { return s_bad ; } // RETURN 00094 const double n = nEff ( histo ) ; 00095 if ( 0 >= n ) { return 0.0 ; } // RETURN 00096 const double a2o = moment ( histo , 2 * order ) ; // a(2o) 00097 const double ao = moment ( histo , 2 * order ) ; // a(o) 00098 double result = a2o - ao*ao ; 00099 result /= n ; 00100 result = std::max ( 0.0 , result ) ; 00101 return std:: sqrt ( result ) ; // RETURN 00102 } 00103 // ====================================================================== 00104 /* evaluate the central momentum (around the mean value) 00105 * @param histo histogram 00106 * @param order the momentm order 00107 * @param value central value 00108 * @return the evaluated central moment 00109 */ 00110 // ====================================================================== 00111 double Gaudi::Utils::HistoStats::centralMoment 00112 ( const AIDA::IHistogram1D* histo , 00113 const unsigned int order ) 00114 { 00115 if ( 0 == histo ) { return s_bad ; } // RETURN 00116 if ( 0 == order ) { return 1.0 ; } // RETURN 00117 if ( 1 == order ) { return 0.0 ; } // RETURN 00118 if ( 2 == order ) 00119 { 00120 const double sigma = histo->rms() ; 00121 return sigma * sigma ; // RETURN 00122 } 00123 // delegate the actual evaluation to another method: 00124 return moment ( histo , order , mean ( histo ) ) ; 00125 } 00126 // ====================================================================== 00127 /* evaluate the uncertanty for 'bin-by-bin'-central momentum 00128 * (around the mean value) 00129 * ( the uncertanty is calculated with O(1/n2) precision) 00130 * @param histo histogram 00131 * @param order the moment parameter 00132 * @param value central value 00133 * @return the evaluated uncertanty in the central moment 00134 */ 00135 // ====================================================================== 00136 double Gaudi::Utils::HistoStats::centralMomentErr 00137 ( const AIDA::IHistogram1D* histo , 00138 const unsigned int order ) 00139 { 00140 if ( 0 == histo ) { return s_bad ; } // RETURN 00141 const double n = nEff ( histo ) ; 00142 if ( 0 >= n ) { return 0.0 ; } // RETURN 00143 const double mu2 = centralMoment ( histo , 2 ) ; // mu(2) 00144 const double muo = centralMoment ( histo , order ) ; // mu(o) 00145 const double mu2o = centralMoment ( histo , 2 * order ) ; // mu(2o) 00146 const double muom = centralMoment ( histo , order - 1 ) ; // mu(o-1) 00147 const double muop = centralMoment ( histo , order + 1 ) ; // mu(o+1) 00148 double result = mu2o ; 00149 result -= 2.0 * order * muom * muop ; 00150 result -= muo * muo ; 00151 result += order * order * mu2 * muom * muom ; 00152 result /= n ; 00153 result = std::max ( 0.0 , result ) ; 00154 return std:: sqrt ( result ) ; // RETURN 00155 } 00156 // ====================================================================== 00157 // get the skewness for the histogram 00158 // ====================================================================== 00159 double Gaudi::Utils::HistoStats::skewness 00160 ( const AIDA::IHistogram1D* histo ) 00161 { 00162 if ( 0 == histo ) { return s_bad ; } // RETURN 00163 const double mu3 = centralMoment ( histo , 3 ) ; 00164 const double s3 = std::pow ( rms ( histo ) , 3 ) ; 00165 return ( std::fabs(s3)>0 ? mu3/s3 : 0.0 ); 00166 } 00167 // ====================================================================== 00168 // get the error in skewness 00169 // ====================================================================== 00170 double Gaudi::Utils::HistoStats::skewnessErr 00171 ( const AIDA::IHistogram1D* histo ) 00172 { 00173 if ( 0 == histo ) { return s_bad ; } // RETURN 00174 const double n = nEff ( histo ) ; 00175 if ( 2 > n ) { return 0.0 ; } // RETURN 00176 double result = 6 ; 00177 result *= ( n - 2 ) ; 00178 result /= ( n + 1 ) * ( n + 3 ) ; 00179 return std::sqrt ( result ) ; 00180 } 00181 // ====================================================================== 00182 // get the kurtosis for the histogram 00183 // ====================================================================== 00184 double Gaudi::Utils::HistoStats::kurtosis 00185 ( const AIDA::IHistogram1D* histo ) 00186 { 00187 if ( 0 == histo ) { return s_bad ; } // RETURN 00188 const double mu4 = centralMoment ( histo , 4 ) ; 00189 const double s4 = std::pow ( rms ( histo ) , 4 ) ; 00190 return ( std::fabs(s4)>0 ? mu4/s4 - 3.0 : 0.0 ); 00191 } 00192 // ====================================================================== 00193 // get the error in kurtosis 00194 // ====================================================================== 00195 double Gaudi::Utils::HistoStats::kurtosisErr 00196 ( const AIDA::IHistogram1D* histo ) 00197 { 00198 if ( 0 == histo ) { return s_bad ; } // RETURN 00199 const double n = nEff ( histo ) ; 00200 if ( 3 > n ) { return 0.0 ; } // RETURN 00201 double result = 24 ; 00202 result *= ( n - 2 ) * ( n - 3 ) ; 00203 result /= ( n + 1 ) * ( n + 1 ) ; 00204 result /= ( n + 3 ) * ( n + 5 ) ; 00205 return std::sqrt ( result ) ; 00206 } 00207 // ====================================================================== 00208 // get the effective entries 00209 // ====================================================================== 00210 double Gaudi::Utils::HistoStats::nEff 00211 ( const AIDA::IHistogram1D* histo ) 00212 { 00213 if ( 0 == histo ) { return s_bad ; } 00214 return histo -> equivalentBinEntries () ; 00215 } 00216 // ====================================================================== 00217 // get the mean value for the histogram 00218 // ====================================================================== 00219 double Gaudi::Utils::HistoStats::mean 00220 ( const AIDA::IHistogram1D* histo ) 00221 { 00222 if ( 0 == histo ) { return s_bad ; } 00223 return histo -> mean() ; 00224 } 00225 // ====================================================================== 00226 // get an error in the mean value 00227 // ====================================================================== 00228 double Gaudi::Utils::HistoStats::meanErr 00229 ( const AIDA::IHistogram1D* histo ) 00230 { 00231 if ( 0 == histo ) { return s_bad ; } 00232 const double n = nEff ( histo ) ; 00233 return 0 >= n ? 0.0 : rms ( histo ) / std::sqrt ( n ) ; 00234 } 00235 // ====================================================================== 00236 // get the rms value for the histogram 00237 // ====================================================================== 00238 double Gaudi::Utils::HistoStats::rms 00239 ( const AIDA::IHistogram1D* histo ) 00240 { 00241 if ( 0 == histo ) { return s_bad ; } 00242 return histo -> rms () ; 00243 } 00244 // ====================================================================== 00245 // get the error in rms 00246 // ====================================================================== 00247 double Gaudi::Utils::HistoStats::rmsErr 00248 ( const AIDA::IHistogram1D* histo ) 00249 { 00250 if ( 0 == histo ) { return s_bad ; } 00251 const double n = nEff ( histo ) ; 00252 if ( 1 >= n ) { return 0.0 ; } 00253 double result = 2 + kurtosis ( histo ) ; 00254 result += 2.0 /( n - 1 ) ; 00255 result /= 4.0 * n ; 00256 result = std::max ( result , 0.0 ) ; 00257 return histo -> rms() * std::sqrt ( result ) ; 00258 } 00259 // ============================================================================ 00260 // The END 00261 // ============================================================================