20 #include <AIDA/IAxis.h>
21 #include <AIDA/IHistogram1D.h>
22 #include <AIDA/IProfile1D.h>
45 double _nEff(
const AIDA::IHistogram1D* histo ) {
return ( histo ? histo->equivalentBinEntries() : s_bad ); }
46 double _nEff(
const AIDA::IProfile1D* histo ) {
48 return ( histo ? histo->sumBinHeights() : s_bad );
53 template <
typename HISTO>
54 double _rms(
const HISTO* histo ) {
55 return ( histo ? histo->rms() : s_bad );
60 template <
typename HISTO>
61 double _mean(
const HISTO* histo ) {
62 return ( histo ? histo->mean() : s_bad );
67 template <
typename HISTO>
68 double _moment(
const HISTO* histo,
const unsigned int order,
const double value = 0 ) {
69 if ( !histo ) {
return s_bad; }
70 if ( 0 == order ) {
return 1.0; }
71 if ( 1 == order ) {
return _mean( histo ) - value; }
73 const auto _r = _rms( histo );
74 const auto _d = value - _mean( histo );
77 const auto n = _nEff( histo );
78 if ( 0 >=
n ) {
return 0.0; }
80 const auto& axis = histo->axis();
82 const auto nBins = axis.bins();
83 double result{ 0 }, weight{ 0 };
85 for (
int i = 0; i < nBins; ++i ) {
86 const auto lE = axis.binLowerEdge( i );
87 const auto uE = axis.binUpperEdge( i );
89 const auto yBin = histo->binHeight( i );
90 const auto xBin = 0.5 * double( lE + uE );
93 result += yBin *
std::pow( xBin - value, order );
101 template <
typename HISTO>
102 double _momentErr(
const HISTO* histo,
const unsigned int order ) {
103 if ( !histo ) {
return s_bad; }
104 const auto n = _nEff( histo );
105 if ( 0 >=
n ) {
return 0.0; }
106 const auto a2o = _moment( histo, 2 * order );
107 const auto ao = _moment( histo, order );
114 template <
typename HISTO>
115 double _centralMoment(
const HISTO* histo,
const unsigned int order ) {
116 if ( !histo ) {
return s_bad; }
117 if ( 0 == order ) {
return 1.0; }
118 if ( 1 == order ) {
return 0.0; }
120 return std::pow( _rms( histo ), 2 );
123 return _moment( histo, order, _mean( histo ) );
128 template <
typename HISTO>
129 double _centralMomentErr(
const HISTO* histo,
const unsigned int order ) {
130 if ( !histo ) {
return s_bad; }
131 const auto n = _nEff( histo );
132 if ( 0 >=
n ) {
return 0.0; }
133 const auto mu2 = _centralMoment( histo, 2 );
134 const auto muo = _centralMoment( histo, order );
135 const auto mu2o = _centralMoment( histo, 2 * order );
136 const auto muom = _centralMoment( histo, order - 1 );
137 const auto muop = _centralMoment( histo, order + 1 );
139 ( mu2o - ( 2.0 * order * muom * muop ) -
std::pow( muo, 2 ) + ( order * order * mu2 * muom * muom ) ) /
n;
145 template <
typename HISTO>
146 double _skewness(
const HISTO* histo ) {
147 if ( !histo ) {
return s_bad; }
148 const auto mu3 = _centralMoment( histo, 3 );
149 const auto s3 =
std::pow( _rms( histo ), 3 );
150 return (
std::fabs( s3 ) > 0 ? mu3 / s3 : 0.0 );
155 template <
typename HISTO>
156 double _skewnessErr(
const HISTO* histo ) {
157 if ( !histo ) {
return s_bad; }
158 const auto n = _nEff( histo );
159 if ( 2 >
n ) {
return 0.0; }
160 const auto result = 6.0 * (
n - 2 ) / ( (
n + 1 ) * (
n + 3 ) );
166 template <
typename HISTO>
167 double _kurtosis(
const HISTO* histo ) {
168 if ( !histo ) {
return s_bad; }
169 const auto mu4 = _centralMoment( histo, 4 );
170 const auto s4 =
std::pow( _rms( histo ), 4 );
171 return (
std::fabs( s4 ) > 0 ? mu4 / s4 - 3.0 : 0.0 );
176 template <
typename HISTO>
177 double _kurtosisErr(
const HISTO* histo ) {
178 if ( !histo ) {
return s_bad; }
179 const auto n = _nEff( histo );
180 if ( 3 >
n ) {
return 0.0; }
181 auto result = 24.0 *
n;
182 result *= (
n - 2 ) * (
n - 3 );
183 result /= (
n + 1 ) * (
n + 1 );
184 result /= (
n + 3 ) * (
n + 5 );
190 template <
typename HISTO>
191 double _meanErr(
const HISTO* histo ) {
192 if ( !histo ) {
return s_bad; }
193 const auto n = _nEff( histo );
194 return ( 0 >=
n ? 0.0 : _rms( histo ) /
std::sqrt(
n ) );
199 template <
typename HISTO>
200 double _rmsErr(
const HISTO* histo ) {
201 if ( !histo ) {
return s_bad; }
202 const auto n = _nEff( histo );
203 if ( 1 >=
n ) {
return 0.0; }
204 auto result = 2.0 + _kurtosis( histo );
205 result += 2.0 / (
n - 1 );
212 template <
typename HISTO>
213 double _sumBinHeightErr(
const HISTO* histo ) {
214 if ( !histo ) {
return s_bad; }
218 const auto& axis = histo->axis();
220 const auto nBins = axis.bins();
222 for (
int i = 0; i < nBins; ++i ) {
224 error2 +=
std::pow( histo->binError( i ), 2 );
231 template <
typename HISTO>
232 double _sumAllBinHeightErr(
const HISTO* histo ) {
233 if ( !histo ) {
return s_bad; }
235 const auto error = _sumBinHeightErr( histo );
236 if ( 0 > error ) {
return s_bad; }
238 const auto err1 = histo->binError( AIDA::IAxis::UNDERFLOW_BIN );
239 const auto err2 = histo->binError( AIDA::IAxis::OVERFLOW_BIN );
246 template <
typename HISTO>
247 double _overflowEntriesFrac(
const HISTO* histo ) {
248 if ( !histo ) {
return s_bad; }
249 const auto overflow = histo->binEntries( AIDA::IAxis::OVERFLOW_BIN );
250 if ( 0 == overflow ) {
return 0; }
251 const auto all = histo->allEntries();
252 if ( 0 ==
all ) {
return 0; }
253 if ( 0 >
all ) {
return s_bad; }
255 return (
double)overflow / (double)
all;
260 template <
typename HISTO>
261 double _underflowEntriesFrac(
const HISTO* histo ) {
262 if ( !histo ) {
return s_bad; }
263 const auto underflow = histo->binEntries( AIDA::IAxis::UNDERFLOW_BIN );
264 if ( 0 == underflow ) {
return 0; }
265 const auto all = histo->allEntries();
266 if ( 0 ==
all ) {
return 0; }
267 if ( 0 >
all ) {
return s_bad; }
269 return (
double)underflow / (double)
all;
274 template <
typename HISTO>
275 double _overflowIntegralFrac(
const HISTO* histo ) {
276 if ( !histo ) {
return s_bad; }
277 const auto overflow = histo->binHeight( AIDA::IAxis::OVERFLOW_BIN );
278 const auto all = histo->sumAllBinHeights();
281 return (
double)overflow / (double)
all;
286 template <
typename HISTO>
287 double _underflowIntegralFrac(
const HISTO* histo ) {
288 if ( !histo ) {
return s_bad; }
289 const auto underflow = histo->binHeight( AIDA::IAxis::UNDERFLOW_BIN );
290 const auto all = histo->sumAllBinHeights();
293 return (
double)underflow / (double)
all;
298 template <
typename HISTO>
299 double _overflowEntriesFracErr(
const HISTO* histo ) {
300 if ( !histo ) {
return s_bad; }
302 const auto overflow = histo->binEntries( AIDA::IAxis::OVERFLOW_BIN );
303 const auto all = histo->allEntries();
305 if ( 0 > overflow || 0 >=
all || overflow >
all ) {
return s_bad; }
307 const double n =
std::max( (
double)overflow, 1.0 );
308 const double N =
all;
316 template <
typename HISTO>
317 double _underflowEntriesFracErr(
const HISTO* histo ) {
318 if ( !histo ) {
return s_bad; }
320 const auto underflow = histo->binEntries( AIDA::IAxis::UNDERFLOW_BIN );
321 const auto all = histo->allEntries();
323 if ( 0 > underflow || 0 >=
all || underflow >
all ) {
return s_bad; }
325 const double n =
std::max( (
double)underflow, 1.0 );
326 const double N =
all;
334 template <
typename HISTO>
335 double _overflowIntegralFracErr(
const HISTO* histo ) {
336 if ( !histo ) {
return s_bad; }
338 const auto all = histo->sumAllBinHeights();
340 const auto overflow = histo->binHeight( AIDA::IAxis::OVERFLOW_BIN );
341 const auto oErr = histo->binError( AIDA::IAxis::OVERFLOW_BIN );
342 if ( 0 > oErr ) {
return s_bad; }
343 const auto aErr = _sumAllBinHeightErr( histo );
344 if ( 0 > aErr ) {
return s_bad; }
346 auto error2 =
std::pow( (
double)oErr, 2 );
355 template <
typename HISTO>
356 double _underflowIntegralFracErr(
const HISTO* histo ) {
357 if ( !histo ) {
return s_bad; }
359 const auto all = histo->sumAllBinHeights();
361 const auto underflow = histo->binHeight( AIDA::IAxis::UNDERFLOW_BIN );
362 const auto oErr = histo->binError( AIDA::IAxis::UNDERFLOW_BIN );
363 if ( 0 > oErr ) {
return s_bad; }
364 const auto aErr = _sumAllBinHeightErr( histo );
365 if ( 0 > aErr ) {
return s_bad; }
367 auto error2 =
std::pow( (
double)oErr, 2 );
376 template <
typename HISTO>
377 long _nEntries(
const HISTO* histo,
const int imax ) {
378 if ( !histo ) {
return s_long_bad; }
379 if ( 0 > imax ) {
return 0; }
380 long result = histo->binEntries( AIDA::IAxis::UNDERFLOW_BIN );
383 const auto& axis = histo->axis();
385 const auto nBins = axis.bins();
387 for (
int i = 0; i < nBins && i < imax; ++i ) { result += histo->binEntries( i ); }
389 if ( nBins < imax ) { result += histo->binEntries( AIDA::IAxis::OVERFLOW_BIN ); }
396 template <
typename HISTO>
397 long _nEntries(
const HISTO* histo,
401 if ( !histo ) {
return s_long_bad; }
402 if ( imin == imax ) {
return 0; }
403 if ( imax < imin ) {
return _nEntries( histo, imax, imin ); }
406 if ( 0 > imin ) { result += histo->binEntries( AIDA::IAxis::UNDERFLOW_BIN ); }
408 const auto& axis = histo->axis();
410 const auto nBins = axis.bins();
411 if ( nBins < imin ) {
return 0; }
413 for (
int i = imin; i < nBins && imin <= i && i < imax; ++i ) { result += histo->binEntries( i ); }
415 if ( nBins < imax ) { result += histo->binEntries( AIDA::IAxis::OVERFLOW_BIN ); }
422 template <
typename HISTO>
423 double _nEntriesFrac(
const HISTO* histo,
const int imax ) {
424 if ( !histo ) {
return s_bad; }
426 const auto N = histo->allEntries();
427 if ( 0 >=
N ) {
return s_bad; }
428 const auto n = _nEntries( histo, imax );
429 if ( 0 >
n ) {
return s_bad; }
430 if (
n >
N ) {
return s_bad; }
432 return (
double)
n / (double)
N;
438 template <
typename HISTO>
439 double _nEntriesFrac(
const HISTO* histo,
443 if ( !histo ) {
return s_bad; }
444 const auto N = histo->allEntries();
445 if ( 0 >=
N ) {
return s_bad; }
446 const auto n = _nEntries( histo, imin, imax );
447 if ( 0 >
n ) {
return s_bad; }
448 if (
n >
N ) {
return s_bad; }
450 return (
double)
n / (double)
N;
455 template <
typename HISTO>
456 double _nEntriesFracErr(
const HISTO* histo,
const int imax ) {
457 if ( !histo ) {
return s_bad; }
459 const auto N = histo->allEntries();
460 if ( 0 >=
N ) {
return s_bad; }
461 const auto n = _nEntries( histo, imax );
462 if ( 0 >
n ) {
return s_bad; }
463 if (
n >
N ) {
return s_bad; }
465 const auto _n1 =
std::max( (
double)
n, 1.0 );
466 const auto _n2 =
std::max( (
double)(
N -
n ), 1.0 );
473 template <
typename HISTO>
474 double _nEntriesFracErr(
const HISTO* histo,
478 if ( !histo ) {
return s_bad; }
480 const auto N = histo->allEntries();
481 if ( 0 >=
N ) {
return s_bad; }
482 const auto n = _nEntries( histo, imin, imax );
483 if ( 0 >
n ) {
return s_bad; }
484 if (
n >
N ) {
return s_bad; }
486 const auto _n1 =
std::max( (
double)
n, 1.0 );
487 const auto _n2 =
std::max( (
double)(
N -
n ), 1.0 );
501 const double value ) {
502 return _moment( histo, order, value );
513 return _moment( histo, order, value );
524 return _momentErr( histo, order );
535 return _momentErr( histo, order );
546 return _centralMoment( histo, order );
557 return _centralMoment( histo, order );
570 return _centralMomentErr( histo, order );
583 return _centralMomentErr( histo, order );
661 return _sumBinHeightErr( histo );
671 return _sumAllBinHeightErr( histo );
677 return _sumAllBinHeightErr( histo );
683 return _overflowEntriesFrac( histo );
689 return _overflowEntriesFrac( histo );
695 return _underflowEntriesFrac( histo );
701 return _underflowEntriesFrac( histo );
707 return _overflowIntegralFrac( histo );
713 return _overflowIntegralFrac( histo );
719 return _underflowIntegralFrac( histo );
725 return _underflowIntegralFrac( histo );
731 return _overflowEntriesFracErr( histo );
737 return _overflowEntriesFracErr( histo );
743 return _underflowEntriesFracErr( histo );
749 return _underflowEntriesFracErr( histo );
755 return _overflowIntegralFracErr( histo );
761 return _overflowIntegralFracErr( histo );
767 return _underflowIntegralFracErr( histo );
773 return _underflowIntegralFracErr( histo );
785 return _nEntries( histo, imax );
797 return _nEntries( histo, imax );
812 return _nEntries( histo, imin, imax );
827 return _nEntries( histo, imin, imax );
839 return _nEntriesFrac( histo, imax );
851 return _nEntriesFrac( histo, imax );
866 return _nEntriesFrac( histo, imin, imax );
881 return _nEntriesFrac( histo, imin, imax );
893 return _nEntriesFracErr( histo, imax );
905 return _nEntriesFracErr( histo, imax );
920 return _nEntriesFracErr( histo, imin, imax );
935 return _nEntriesFracErr( histo, imin, imax );