4 # pragma warning( disable : 1572 ) 17 #include "AIDA/IAxis.h" 18 #include "AIDA/IHistogram1D.h" 19 #include "AIDA/IProfile1D.h" 42 double _nEff(
const AIDA::IHistogram1D* histo ) {
return ( histo ? histo->equivalentBinEntries() : s_bad ); }
43 double _nEff(
const AIDA::IProfile1D* histo ) {
45 return ( histo ? histo->sumBinHeights() : s_bad );
50 template <
typename HISTO>
51 double _rms(
const HISTO* histo ) {
52 return ( histo ? histo->rms() : s_bad );
57 template <
typename HISTO>
58 double _mean(
const HISTO* histo ) {
59 return ( histo ? histo->mean() : s_bad );
64 template <
typename HISTO>
65 double _moment(
const HISTO* histo,
const unsigned int order,
const double value = 0 ) {
66 if (
UNLIKELY( !histo ) ) {
return s_bad; }
67 if ( 0 == order ) {
return 1.0; }
68 if ( 1 == order ) {
return _mean( histo ) - value; }
70 const auto _r = _rms( histo );
71 const auto _d = value - _mean( histo );
74 const auto n = _nEff( histo );
75 if ( 0 >=
n ) {
return 0.0; }
77 const auto& axis = histo->axis();
79 const auto nBins = axis.bins();
80 double result{0}, weight{0};
82 for (
int i = 0; i < nBins; ++i ) {
83 const auto lE = axis.binLowerEdge( i );
84 const auto uE = axis.binUpperEdge( i );
86 const auto yBin = histo->binHeight( i );
87 const auto xBin = 0.5 * double( lE + uE );
90 result += yBin *
std::pow( xBin - value, order );
92 if ( 0 != weight ) { result /= weight; }
98 template <
typename HISTO>
99 double _momentErr(
const HISTO* histo,
const unsigned int order ) {
100 if (
UNLIKELY( !histo ) ) {
return s_bad; }
101 const auto n = _nEff( histo );
102 if ( 0 >=
n ) {
return 0.0; }
103 const auto a2o = _moment( histo, 2 * order );
104 const auto ao = _moment( histo, order );
111 template <
typename HISTO>
112 double _centralMoment(
const HISTO* histo,
const unsigned int order ) {
113 if (
UNLIKELY( !histo ) ) {
return s_bad; }
114 if ( 0 == order ) {
return 1.0; }
115 if ( 1 == order ) {
return 0.0; }
117 return std::pow( _rms( histo ), 2 );
120 return _moment( histo, order, _mean( histo ) );
125 template <
typename HISTO>
126 double _centralMomentErr(
const HISTO* histo,
const unsigned int order ) {
127 if (
UNLIKELY( !histo ) ) {
return s_bad; }
128 const auto n = _nEff( histo );
129 if ( 0 >=
n ) {
return 0.0; }
130 const auto mu2 = _centralMoment( histo, 2 );
131 const auto muo = _centralMoment( histo, order );
132 const auto mu2o = _centralMoment( histo, 2 * order );
133 const auto muom = _centralMoment( histo, order - 1 );
134 const auto muop = _centralMoment( histo, order + 1 );
136 ( mu2o - ( 2.0 * order * muom * muop ) -
std::pow( muo, 2 ) + ( order * order * mu2 * muom * muom ) ) /
n;
142 template <
typename HISTO>
143 double _skewness(
const HISTO* histo ) {
144 if (
UNLIKELY( !histo ) ) {
return s_bad; }
145 const auto mu3 = _centralMoment( histo, 3 );
146 const auto s3 =
std::pow( _rms( histo ), 3 );
147 return (
std::fabs( s3 ) > 0 ? mu3 / s3 : 0.0 );
152 template <
typename HISTO>
153 double _skewnessErr(
const HISTO* histo ) {
154 if (
UNLIKELY( !histo ) ) {
return s_bad; }
155 const auto n = _nEff( histo );
156 if ( 2 >
n ) {
return 0.0; }
157 const auto result = 6.0 * (
n - 2 ) / ( (
n + 1 ) * (
n + 3 ) );
163 template <
typename HISTO>
164 double _kurtosis(
const HISTO* histo ) {
165 if (
UNLIKELY( !histo ) ) {
return s_bad; }
166 const auto mu4 = _centralMoment( histo, 4 );
167 const auto s4 =
std::pow( _rms( histo ), 4 );
168 return (
std::fabs( s4 ) > 0 ? mu4 / s4 - 3.0 : 0.0 );
173 template <
typename HISTO>
174 double _kurtosisErr(
const HISTO* histo ) {
175 if (
UNLIKELY( !histo ) ) {
return s_bad; }
176 const auto n = _nEff( histo );
177 if ( 3 >
n ) {
return 0.0; }
178 auto result = 24.0 *
n;
179 result *= (
n - 2 ) * (
n - 3 );
180 result /= (
n + 1 ) * (
n + 1 );
181 result /= (
n + 3 ) * (
n + 5 );
187 template <
typename HISTO>
188 double _meanErr(
const HISTO* histo ) {
189 if (
UNLIKELY( !histo ) ) {
return s_bad; }
190 const auto n = _nEff( histo );
191 return ( 0 >=
n ? 0.0 : _rms( histo ) /
std::sqrt(
n ) );
196 template <
typename HISTO>
197 double _rmsErr(
const HISTO* histo ) {
198 if (
UNLIKELY( !histo ) ) {
return s_bad; }
199 const auto n = _nEff( histo );
200 if ( 1 >=
n ) {
return 0.0; }
201 auto result = 2.0 + _kurtosis( histo );
202 result += 2.0 / (
n - 1 );
209 template <
typename HISTO>
210 double _sumBinHeightErr(
const HISTO* histo ) {
211 if (
UNLIKELY( !histo ) ) {
return s_bad; }
215 const auto& axis = histo->axis();
217 const auto nBins = axis.bins();
219 for (
int i = 0; i < nBins; ++i ) {
221 error2 +=
std::pow( histo->binError( i ), 2 );
228 template <
typename HISTO>
229 double _sumAllBinHeightErr(
const HISTO* histo ) {
230 if (
UNLIKELY( !histo ) ) {
return s_bad; }
232 const auto error = _sumBinHeightErr( histo );
233 if ( 0 > error ) {
return s_bad; }
235 const auto err1 = histo->binError( AIDA::IAxis::UNDERFLOW_BIN );
236 const auto err2 = histo->binError( AIDA::IAxis::OVERFLOW_BIN );
243 template <
typename HISTO>
244 double _overflowEntriesFrac(
const HISTO* histo ) {
245 if (
UNLIKELY( !histo ) ) {
return s_bad; }
246 const auto overflow = histo->binEntries( AIDA::IAxis::OVERFLOW_BIN );
247 if ( 0 == overflow ) {
return 0; }
248 const auto all = histo->allEntries();
249 if ( 0 == all ) {
return 0; }
250 if ( 0 > all ) {
return s_bad; }
252 return (
double)overflow / (double)all;
257 template <
typename HISTO>
258 double _underflowEntriesFrac(
const HISTO* histo ) {
259 if (
UNLIKELY( !histo ) ) {
return s_bad; }
260 const auto underflow = histo->binEntries( AIDA::IAxis::UNDERFLOW_BIN );
261 if ( 0 == underflow ) {
return 0; }
262 const auto all = histo->allEntries();
263 if ( 0 == all ) {
return 0; }
264 if ( 0 > all ) {
return s_bad; }
266 return (
double)underflow / (double)all;
271 template <
typename HISTO>
272 double _overflowIntegralFrac(
const HISTO* histo ) {
273 if (
UNLIKELY( !histo ) ) {
return s_bad; }
274 const auto overflow = histo->binHeight( AIDA::IAxis::OVERFLOW_BIN );
275 if ( 0 == overflow ) {
return 0; }
276 const auto all = histo->sumAllBinHeights();
277 if ( 0 == all ) {
return 0; }
279 return (
double)overflow / (double)all;
284 template <
typename HISTO>
285 double _underflowIntegralFrac(
const HISTO* histo ) {
286 if (
UNLIKELY( !histo ) ) {
return s_bad; }
287 const auto underflow = histo->binHeight( AIDA::IAxis::UNDERFLOW_BIN );
288 if ( 0 == underflow ) {
return 0; }
289 const auto all = histo->sumAllBinHeights();
290 if ( 0 == all ) {
return 0; }
292 return (
double)underflow / (double)all;
297 template <
typename HISTO>
298 double _overflowEntriesFracErr(
const HISTO* histo ) {
299 if (
UNLIKELY( !histo ) ) {
return s_bad; }
301 const auto overflow = histo->binEntries( AIDA::IAxis::OVERFLOW_BIN );
302 const auto all = histo->allEntries();
304 if ( 0 > overflow || 0 >= all || overflow > all ) {
return s_bad; }
306 const double n =
std::max( (
double)overflow, 1.0 );
307 const double N = all;
308 const double n1 =
std::max( N - n, 1.0 );
315 template <
typename HISTO>
316 double _underflowEntriesFracErr(
const HISTO* histo ) {
317 if (
UNLIKELY( !histo ) ) {
return s_bad; }
319 const auto underflow = histo->binEntries( AIDA::IAxis::UNDERFLOW_BIN );
320 const auto all = histo->allEntries();
322 if ( 0 > underflow || 0 >= all || underflow > all ) {
return s_bad; }
324 const double n =
std::max( (
double)underflow, 1.0 );
325 const double N = all;
326 const double n1 =
std::max( N - n, 1.0 );
333 template <
typename HISTO>
334 double _overflowIntegralFracErr(
const HISTO* histo ) {
335 if (
UNLIKELY( !histo ) ) {
return s_bad; }
337 const auto all = histo->sumAllBinHeights();
338 if ( 0 == all ) {
return s_bad; }
339 const auto overflow = histo->binHeight( AIDA::IAxis::OVERFLOW_BIN );
340 const auto oErr = histo->binError( AIDA::IAxis::OVERFLOW_BIN );
341 if ( 0 > oErr ) {
return s_bad; }
342 const auto aErr = _sumAllBinHeightErr( histo );
343 if ( 0 > aErr ) {
return s_bad; }
345 auto error2 =
std::pow( (
double)oErr, 2 );
347 error2 /=
std::pow( (
double)all, 2 );
354 template <
typename HISTO>
355 double _underflowIntegralFracErr(
const HISTO* histo ) {
356 if (
UNLIKELY( !histo ) ) {
return s_bad; }
358 const auto all = histo->sumAllBinHeights();
359 if ( 0 == all ) {
return s_bad; }
360 const auto underflow = histo->binHeight( AIDA::IAxis::UNDERFLOW_BIN );
361 const auto oErr = histo->binError( AIDA::IAxis::UNDERFLOW_BIN );
362 if ( 0 > oErr ) {
return s_bad; }
363 const auto aErr = _sumAllBinHeightErr( histo );
364 if ( 0 > aErr ) {
return s_bad; }
366 auto error2 =
std::pow( (
double)oErr, 2 );
368 error2 /=
std::pow( (
double)all, 2 );
375 template <
typename HISTO>
376 long _nEntries(
const HISTO* histo,
const int imax ) {
377 if (
UNLIKELY( !histo ) ) {
return s_long_bad; }
378 if ( 0 > imax ) {
return 0; }
379 long result = histo->binEntries( AIDA::IAxis::UNDERFLOW_BIN );
382 const auto& axis = histo->axis();
384 const auto nBins = axis.bins();
386 for (
int i = 0; i < nBins && i < imax; ++i ) { result += histo->binEntries( i ); }
388 if ( nBins < imax ) { result += histo->binEntries( AIDA::IAxis::OVERFLOW_BIN ); }
395 template <
typename HISTO>
396 long _nEntries(
const HISTO* histo,
400 if (
UNLIKELY( !histo ) ) {
return s_long_bad; }
401 if ( imin == imax ) {
return 0; }
402 if ( imax < imin ) {
return _nEntries( histo, imax, imin ); }
405 if ( 0 > imin ) { result += histo->binEntries( AIDA::IAxis::UNDERFLOW_BIN ); }
407 const auto& axis = histo->axis();
409 const auto nBins = axis.bins();
410 if ( nBins < imin ) {
return 0; }
412 for (
int i = imin; i < nBins && imin <= i && i < imax; ++i ) { result += histo->binEntries( i ); }
414 if ( nBins < imax ) { result += histo->binEntries( AIDA::IAxis::OVERFLOW_BIN ); }
421 template <
typename HISTO>
422 double _nEntriesFrac(
const HISTO* histo,
const int imax ) {
423 if (
UNLIKELY( !histo ) ) {
return s_bad; }
425 const auto N = histo->allEntries();
426 if ( 0 >= N ) {
return s_bad; }
427 const auto n = _nEntries( histo, imax );
428 if ( 0 > n ) {
return s_bad; }
429 if ( n > N ) {
return s_bad; }
431 return (
double)n / (double)N;
437 template <
typename HISTO>
438 double _nEntriesFrac(
const HISTO* histo,
442 if (
UNLIKELY( !histo ) ) {
return s_bad; }
443 const auto N = histo->allEntries();
444 if ( 0 >= N ) {
return s_bad; }
445 const auto n = _nEntries( histo, imin, imax );
446 if ( 0 > n ) {
return s_bad; }
447 if ( n > N ) {
return s_bad; }
449 return (
double)n / (double)N;
454 template <
typename HISTO>
455 double _nEntriesFracErr(
const HISTO* histo,
const int imax ) {
456 if (
UNLIKELY( !histo ) ) {
return s_bad; }
458 const auto N = histo->allEntries();
459 if ( 0 >= N ) {
return s_bad; }
460 const auto n = _nEntries( histo, imax );
461 if ( 0 > n ) {
return s_bad; }
462 if ( n > N ) {
return s_bad; }
464 const auto _n1 =
std::max( (
double)n, 1.0 );
465 const auto _n2 =
std::max( (
double)( N - n ), 1.0 );
472 template <
typename HISTO>
473 double _nEntriesFracErr(
const HISTO* histo,
477 if (
UNLIKELY( !histo ) ) {
return s_bad; }
479 const auto N = histo->allEntries();
480 if ( 0 >= N ) {
return s_bad; }
481 const auto n = _nEntries( histo, imin, imax );
482 if ( 0 > n ) {
return s_bad; }
483 if ( n > N ) {
return s_bad; }
485 const auto _n1 =
std::max( (
double)n, 1.0 );
486 const auto _n2 =
std::max( (
double)( N - n ), 1.0 );
500 const double value ) {
501 return _moment( histo, order, value );
512 return _moment( histo, order, value );
523 return _momentErr( histo, order );
534 return _momentErr( histo, order );
545 return _centralMoment( histo, order );
556 return _centralMoment( histo, order );
569 return _centralMomentErr( histo, order );
582 return _centralMomentErr( histo, order );
660 return _sumBinHeightErr( histo );
670 return _sumAllBinHeightErr( histo );
676 return _sumAllBinHeightErr( histo );
682 return _overflowEntriesFrac( histo );
688 return _overflowEntriesFrac( histo );
694 return _underflowEntriesFrac( histo );
700 return _underflowEntriesFrac( histo );
706 return _overflowIntegralFrac( histo );
712 return _overflowIntegralFrac( histo );
718 return _underflowIntegralFrac( histo );
724 return _underflowIntegralFrac( histo );
730 return _overflowEntriesFracErr( histo );
736 return _overflowEntriesFracErr( histo );
742 return _underflowEntriesFracErr( histo );
748 return _underflowEntriesFracErr( histo );
754 return _overflowIntegralFracErr( histo );
760 return _overflowIntegralFracErr( histo );
766 return _underflowIntegralFracErr( histo );
772 return _underflowIntegralFracErr( histo );
784 return _nEntries( histo, imax );
796 return _nEntries( histo, imax );
811 return _nEntries( histo, imin, imax );
826 return _nEntries( histo, imin, imax );
838 return _nEntriesFrac( histo, imax );
850 return _nEntriesFrac( histo, imax );
865 return _nEntriesFrac( histo, imin, imax );
880 return _nEntriesFrac( histo, imin, imax );
892 return _nEntriesFracErr( histo, imax );
904 return _nEntriesFracErr( histo, imax );
919 return _nEntriesFracErr( histo, imin, imax );
934 return _nEntriesFracErr( histo, imin, imax );
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