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 );