4 #pragma warning(disable:1572)
17 #include "AIDA/IHistogram1D.h"
18 #include "AIDA/IAxis.h"
22 #include "GaudiUtils/HistoStats.h"
34 const double s_bad = -1 * std::numeric_limits<double>::max() ;
36 const long s_long_bad = std::numeric_limits<int>::min() ;
48 (
const AIDA::IHistogram1D* histo ,
49 const unsigned int order ,
52 if ( 0 == histo ) {
return s_bad ; }
53 if ( 0 == order ) {
return 1.0 ; }
54 if ( 1 == order ) {
return mean( histo ) -
value ; }
57 const double _r = rms ( histo ) ;
58 const double _d = value - mean ( histo ) ;
59 return _r *_r + _d * _d ;
61 const double n = nEff ( histo ) ;
62 if ( 0 >= n ) {
return 0.0 ; }
64 const int _order = (int) order ;
67 const AIDA::IAxis& axis = histo -> axis () ;
69 const int nBins = axis.bins() ;
73 for (
int i = 0 ;
i < nBins ; ++
i )
75 const double lE = axis . binLowerEdge (
i ) ;
76 const double uE = axis . binUpperEdge (
i ) ;
78 const double yBin = histo -> binHeight (
i ) ;
79 const double xBin = 0.5 * ( lE + uE ) ;
82 result += yBin * std::pow ( xBin - value , _order ) ;
84 if ( 0 != weight ) { result /= weight ; }
96 (
const AIDA::IHistogram1D* histo ,
97 const unsigned int order )
99 if ( 0 == histo ) {
return s_bad ; }
100 const double n = nEff ( histo ) ;
101 if ( 0 >= n ) {
return 0.0 ; }
102 const double a2o = moment ( histo , 2 * order ) ;
103 const double ao = moment ( histo , order ) ;
104 double result = a2o - ao*ao ;
106 result = std::max ( 0.0 , result ) ;
107 return std:: sqrt ( result ) ;
118 (
const AIDA::IHistogram1D* histo ,
119 const unsigned int order )
121 if ( 0 == histo ) {
return s_bad ; }
122 if ( 0 == order ) {
return 1.0 ; }
123 if ( 1 == order ) {
return 0.0 ; }
126 const double sigma = histo->rms() ;
127 return sigma * sigma ;
130 return moment ( histo , order , mean ( histo ) ) ;
143 (
const AIDA::IHistogram1D* histo ,
144 const unsigned int order )
146 if ( 0 == histo ) {
return s_bad ; }
147 const double n = nEff ( histo ) ;
148 if ( 0 >= n ) {
return 0.0 ; }
149 const double mu2 = centralMoment ( histo , 2 ) ;
150 const double muo = centralMoment ( histo , order ) ;
151 const double mu2o = centralMoment ( histo , 2 * order ) ;
152 const double muom = centralMoment ( histo , order - 1 ) ;
153 const double muop = centralMoment ( histo , order + 1 ) ;
154 double result = mu2o ;
155 result -= 2.0 * order * muom * muop ;
156 result -= muo * muo ;
157 result += order * order * mu2 * muom * muom ;
159 result = std::max ( 0.0 , result ) ;
160 return std:: sqrt ( result ) ;
166 (
const AIDA::IHistogram1D* histo )
168 if ( 0 == histo ) {
return s_bad ; }
169 const double mu3 = centralMoment ( histo , 3 ) ;
170 const double s3 = std::pow ( rms ( histo ) , 3 ) ;
171 return ( std::fabs(s3)>0 ? mu3/s3 : 0.0 );
177 (
const AIDA::IHistogram1D* histo )
179 if ( 0 == histo ) {
return s_bad ; }
180 const double n = nEff ( histo ) ;
181 if ( 2 > n ) {
return 0.0 ; }
183 result *= ( n - 2 ) ;
184 result /= ( n + 1 ) * ( n + 3 ) ;
185 return std::sqrt ( result ) ;
191 (
const AIDA::IHistogram1D* histo )
193 if ( 0 == histo ) {
return s_bad ; }
194 const double mu4 = centralMoment ( histo , 4 ) ;
195 const double s4 = std::pow ( rms ( histo ) , 4 ) ;
196 return ( std::fabs(s4)>0 ? mu4/s4 - 3.0 : 0.0 );
202 (
const AIDA::IHistogram1D* histo )
204 if ( 0 == histo ) {
return s_bad ; }
205 const double n = nEff ( histo ) ;
206 if ( 3 > n ) {
return 0.0 ; }
207 double result = 24 *
n ;
208 result *= ( n - 2 ) * ( n - 3 ) ;
209 result /= ( n + 1 ) * ( n + 1 ) ;
210 result /= ( n + 3 ) * ( n + 5 ) ;
211 return std::sqrt ( result ) ;
217 (
const AIDA::IHistogram1D* histo )
219 if ( 0 == histo ) {
return s_bad ; }
220 return histo -> equivalentBinEntries () ;
226 (
const AIDA::IHistogram1D* histo )
228 if ( 0 == histo ) {
return s_bad ; }
229 return histo -> mean() ;
235 (
const AIDA::IHistogram1D* histo )
237 if ( 0 == histo ) {
return s_bad ; }
238 const double n = nEff ( histo ) ;
239 return 0 >= n ? 0.0 : rms ( histo ) / std::sqrt ( n ) ;
245 (
const AIDA::IHistogram1D* histo )
247 if ( 0 == histo ) {
return s_bad ; }
248 return histo -> rms () ;
254 (
const AIDA::IHistogram1D* histo )
256 if ( 0 == histo ) {
return s_bad ; }
257 const double n = nEff ( histo ) ;
258 if ( 1 >= n ) {
return 0.0 ; }
259 double result = 2 + kurtosis ( histo ) ;
260 result += 2.0 /( n - 1 ) ;
262 result = std::max ( result , 0.0 ) ;
263 return histo -> rms() * std::sqrt ( result ) ;
269 (
const AIDA::IHistogram1D* histo )
271 if ( 0 == histo ) {
return s_bad ; }
275 const AIDA::IAxis& axis = histo -> axis () ;
277 const int nBins = axis.bins() ;
279 for (
int i = 0 ;
i < nBins ; ++
i )
282 const double berr = histo->binError (
i ) ;
284 error2 += berr * berr ;
286 return std::sqrt ( error2 ) ;
292 (
const AIDA::IHistogram1D* histo )
294 if ( 0 == histo ) {
return s_bad ; }
296 const double error = sumBinHeightErr( histo ) ;
297 if ( 0 > error ) {
return s_bad ; }
299 const double err1 = histo->binError ( AIDA::IAxis::UNDERFLOW_BIN ) ;
300 const double err2 = histo->binError ( AIDA::IAxis::OVERFLOW_BIN ) ;
302 return std::sqrt ( error * error + err1 * err1 + err2 * err2 ) ;
308 (
const AIDA::IHistogram1D* histo )
310 if ( 0 == histo ) {
return s_bad ; }
311 const long overflow = histo -> binEntries ( AIDA::IAxis::OVERFLOW_BIN ) ;
312 if ( 0 == overflow ) {
return 0 ; }
313 const long all = histo->allEntries() ;
314 if ( 0 == all ) {
return 0 ; }
315 if ( 0 > all ) {
return s_bad ; }
317 double _tmp = (double) overflow ;
326 (
const AIDA::IHistogram1D* histo )
328 if ( 0 == histo ) {
return s_bad ; }
329 const long underflow = histo -> binEntries ( AIDA::IAxis::UNDERFLOW_BIN ) ;
330 if ( 0 == underflow ) {
return 0 ; }
331 const long all = histo->allEntries() ;
332 if ( 0 == all ) {
return 0 ; }
333 if ( 0 > all ) {
return s_bad ; }
335 double _tmp = (double) underflow ;
344 (
const AIDA::IHistogram1D* histo )
346 if ( 0 == histo ) {
return s_bad ; }
347 const double overflow = histo -> binHeight ( AIDA::IAxis::OVERFLOW_BIN ) ;
348 if ( 0 == overflow ) {
return 0 ; }
349 const double all = histo -> sumAllBinHeights() ;
350 if ( 0 == all ) {
return 0 ; }
352 return overflow / all ;
358 (
const AIDA::IHistogram1D* histo )
360 if ( 0 == histo ) {
return s_bad ; }
361 const double underflow = histo -> binHeight ( AIDA::IAxis::UNDERFLOW_BIN ) ;
362 if ( 0 == underflow ) {
return 0 ; }
363 const double all = histo -> sumAllBinHeights() ;
364 if ( 0 == all ) {
return 0 ; }
366 return underflow / all ;
372 (
const AIDA::IHistogram1D* histo )
374 if ( 0 == histo ) {
return s_bad ; }
376 const double overflow = histo -> binEntries ( AIDA::IAxis::OVERFLOW_BIN );
377 const double all = histo -> allEntries () ;
379 if ( 0 > overflow || 0 >= all || overflow > all ) {
return s_bad ; }
381 const double n = std::max ( overflow , 1.0 ) ;
382 const double N = all ;
383 const double n1 = std::max ( N - n , 1.0 ) ;
385 return std::sqrt ( n * n1 / N ) /
N ;
391 (
const AIDA::IHistogram1D* histo )
393 if ( 0 == histo ) {
return s_bad ; }
395 const double underflow = histo -> binEntries ( AIDA::IAxis::UNDERFLOW_BIN );
396 const double all = histo -> allEntries () ;
398 if ( 0 > underflow || 0 >= all || underflow > all ) {
return s_bad ; }
400 const double n = std::max ( underflow , 1.0 ) ;
401 const double N = all ;
402 const double n1 = std::max ( N - n , 1.0 ) ;
404 return std::sqrt ( n * n1 / N ) /
N ;
410 (
const AIDA::IHistogram1D* histo )
412 if ( 0 == histo ) {
return s_bad ; }
414 const double all = histo -> sumAllBinHeights () ;
415 if ( 0 == all ) {
return s_bad ; }
416 const double overflow = histo -> binHeight ( AIDA::IAxis::OVERFLOW_BIN ) ;
417 const double oErr = histo -> binError ( AIDA::IAxis::OVERFLOW_BIN ) ;
418 if ( 0 > oErr ) {
return s_bad ; }
419 const double aErr = sumAllBinHeightErr ( histo ) ;
420 if ( 0 > aErr ) {
return s_bad ; }
422 double error2 = oErr * oErr ;
423 error2 += aErr * aErr * overflow * overflow / all / all ;
424 error2 /= all * all ;
426 return std::sqrt ( error2 ) ;
432 (
const AIDA::IHistogram1D* histo )
434 if ( 0 == histo ) {
return s_bad ; }
436 const double all = histo -> sumAllBinHeights () ;
437 if ( 0 == all ) {
return s_bad ; }
438 const double underflow = histo -> binHeight ( AIDA::IAxis::UNDERFLOW_BIN ) ;
439 const double oErr = histo -> binError ( AIDA::IAxis::UNDERFLOW_BIN ) ;
440 if ( 0 > oErr ) {
return s_bad ; }
441 const double aErr = sumAllBinHeightErr ( histo ) ;
442 if ( 0 > aErr ) {
return s_bad ; }
444 double error2 = oErr * oErr ;
445 error2 += aErr * aErr * underflow * underflow / all / all ;
446 error2 /= all * all ;
448 return std::sqrt ( error2 ) ;
460 (
const AIDA::IHistogram1D* histo ,
463 if ( 0 == histo ) {
return s_long_bad ; }
464 if ( 0 > imax ) {
return 0 ; }
465 long result = histo -> binEntries( AIDA::IAxis::UNDERFLOW_BIN ) ;
468 const AIDA::IAxis& axis = histo -> axis () ;
470 const int nBins = axis.bins() ;
472 for (
int i = 0 ;
i < nBins ; ++
i )
473 {
if (
i < imax ) { result += histo->binEntries (
i ) ; } }
476 { result += histo -> binEntries( AIDA::IAxis::OVERFLOW_BIN ) ; }
490 (
const AIDA::IHistogram1D* histo ,
494 if ( 0 == histo ) {
return s_long_bad ; }
495 if ( imin == imax ) {
return 0 ; }
496 if ( imax < imin ) {
return nEntries ( histo , imax ,imin ) ; }
500 { result += histo -> binEntries( AIDA::IAxis::UNDERFLOW_BIN ) ; }
502 const AIDA::IAxis& axis = histo -> axis () ;
504 const int nBins = axis.bins() ;
505 if ( nBins < imin ) {
return 0 ; }
507 for (
int i = 0 ;
i < nBins ; ++
i )
508 {
if ( imin <=
i &&
i < imax ) { result += histo->binEntries (
i ) ; } }
511 { result += histo -> binEntries( AIDA::IAxis::OVERFLOW_BIN ) ; }
525 (
const AIDA::IHistogram1D* histo ,
528 if ( 0 == histo ) {
return s_bad ; }
530 const double N = histo->allEntries () ;
531 if ( 0 >= N ) {
return s_bad ; }
532 const long n = nEntries ( histo , imax ) ;
533 if ( 0 > n ) {
return s_bad ; }
534 if ( n > N ) {
return s_bad ; }
548 (
const AIDA::IHistogram1D* histo ,
552 if ( 0 == histo ) {
return s_bad ; }
553 const double N = histo->allEntries () ;
554 if ( 0 >= N ) {
return s_bad ; }
555 const long n = nEntries ( histo , imin , imax ) ;
556 if ( 0 > n ) {
return s_bad ; }
557 if ( n > N ) {
return s_bad ; }
571 (
const AIDA::IHistogram1D* histo ,
574 if ( 0 == histo ) {
return s_bad ; }
576 const double N = histo->allEntries () ;
577 if ( 0 >= N ) {
return s_bad ; }
578 const long n = nEntries ( histo , imax ) ;
579 if ( 0 > n ) {
return s_bad ; }
580 if ( n > N ) {
return s_bad ; }
582 const double _n1 = std::max ( (
double) n , 1.0 ) ;
583 const double _n2 = std::max ( ( N - n ) , 1.0 ) ;
585 return std::sqrt ( _n1 * _n2 / N ) /
N ;
597 (
const AIDA::IHistogram1D* histo ,
601 if ( 0 == histo ) {
return s_bad ; }
603 const double N = histo->allEntries () ;
604 if ( 0 >= N ) {
return s_bad ; }
605 const long n = nEntries ( histo , imin , imax ) ;
606 if ( 0 > n ) {
return s_bad ; }
607 if ( n > N ) {
return s_bad ; }
609 const double _n1 = std::max ( (
double) n , 1.0 ) ;
610 const double _n2 = std::max ( ( N - n ) , 1.0 ) ;
612 return std::sqrt ( _n1 * _n2 / N ) /
N ;
static double underflowIntegralFracErr(const AIDA::IHistogram1D *histo)
the error on fraction of underflow integral
static double overflowIntegralFracErr(const AIDA::IHistogram1D *histo)
the error on fraction of overflow intergal
static double momentErr(const AIDA::IHistogram1D *histo, const unsigned int order)
evaluate the uncertanty for 'bin-by-bin'-moment
static double underflowEntriesFrac(const AIDA::IHistogram1D *histo)
the fraction of underflow entries (useful for shape comparison)
static long nEntries(const AIDA::IHistogram1D *histo, const int imax)
get number of entries in histogram up to the certain bin (not-included)
static double rmsErr(const AIDA::IHistogram1D *histo)
get an error in the rms value
static double underflowIntegralFrac(const AIDA::IHistogram1D *histo)
the fraction of underflow integral (useful for shape comparison)
static double moment(const AIDA::IHistogram1D *histo, const unsigned int order, const double value=0)
get the "bin-by-bin"-moment around the specified "value"
static double meanErr(const AIDA::IHistogram1D *histo)
get an error in the mean value
static double nEntriesFrac(const AIDA::IHistogram1D *histo, const int imax)
get the fraction of entries in histogram up to the certain bin (not-included)
static double overflowIntegralFrac(const AIDA::IHistogram1D *histo)
the fraction of overflow intergal (useful for shape comparison)
static double centralMoment(const AIDA::IHistogram1D *histo, const unsigned int order)
evaluate the 'bin-by-bin'-central moment (around the mean value)
static double centralMomentErr(const AIDA::IHistogram1D *histo, const unsigned int order)
evaluate the uncertanty for 'bin-by-bin'-central moment (around the mean value) ( the uncertanty is c...
static double rms(const AIDA::IHistogram1D *histo)
get the rms value for the histogram (just for completeness)
static double overflowEntriesFrac(const AIDA::IHistogram1D *histo)
the fraction of overflow entries (useful for shape comparison)
static double nEntriesFracErr(const AIDA::IHistogram1D *histo, const int imax)
get the (binominal) error for the fraction of entries in histogram up to the certain bin (not-include...
static double sumBinHeightErr(const AIDA::IHistogram1D *histo)
get an error in the sum bin height ("in-range integral")
static double overflowEntriesFracErr(const AIDA::IHistogram1D *histo)
error on fraction of overflow entries (useful for shape comparison)
static double sumAllBinHeightErr(const AIDA::IHistogram1D *histo)
get an error in the sum of all bin height ("integral")
static double kurtosisErr(const AIDA::IHistogram1D *histo)
get the error in kurtosis for the histogram
static double mean(const AIDA::IHistogram1D *histo)
get the mean value for the histogram (just for completeness)
static double skewnessErr(const AIDA::IHistogram1D *histo)
get the error in skewness for the histogram
static double underflowEntriesFracErr(const AIDA::IHistogram1D *histo)
the error on fraction of underflow entries (useful for shape comparison)
static double kurtosis(const AIDA::IHistogram1D *histo)
get the kurtosis for the histogram
static double nEff(const AIDA::IHistogram1D *histo)
get the effective entries (just for completeness)
static double skewness(const AIDA::IHistogram1D *histo)
get the skewness for the histogram