RooFit/Detail/MathFuncs.h

0001 /*
0002  * Project: RooFit
0003  * Authors:
0004  *   Jonas Rembser, CERN 2024
0005  *   Garima Singh, CERN 2023
0006  *
0007  * Copyright (c) 2024, CERN
0008  *
0009  * Redistribution and use in source and binary forms,
0010  * with or without modification, are permitted according to the terms
0011  * listed in LICENSE (http://roofit.sourceforge.net/license.txt)
0012  */
0013
0014 #ifndef RooFit_Detail_MathFuncs_h
0015 #define RooFit_Detail_MathFuncs_h
0016
0017 #include <TMath.h>
0018 #include <Math/PdfFuncMathCore.h>
0019 #include <Math/ProbFuncMathCore.h>
0020
0021 #include <cmath>
0022
0023 namespace RooFit {
0024
0025 namespace Detail {
0026
0027 namespace MathFuncs {
0028
0029 /// Calculates the binomial coefficient n over k.
0030 /// Equivalent to TMath::Binomial, but inlined.
0031 inline double binomial(int n, int k)
0032 {
0033    if (n < 0 || k < 0 || n < k)
0034       return TMath::SignalingNaN();
0035    if (k == 0 || n == k)
0036       return 1;
0037
0038    int k1 = std::min(k, n - k);
0039    int k2 = n - k1;
0040    double fact = k2 + 1;
0041    for (double i = k1; i > 1.; --i) {
0042       fact *= (k2 + i) / i;
0043    }
0044    return fact;
0045 }
0046
0047 /// The caller needs to make sure that there is at least one coefficient.
0048 inline double bernstein(double x, double xmin, double xmax, double *coefs, int nCoefs)
0049 {
0050    double xScaled = (x - xmin) / (xmax - xmin); // rescale to [0,1]
0051    int degree = nCoefs - 1;                     // n+1 polys of degree n
0052
0053    // in case list of arguments passed is empty
0054    if (degree < 0) {
0055       return TMath::SignalingNaN();
0056    } else if (degree == 0) {
0057       return coefs[0];
0058    } else if (degree == 1) {
0059
0060       double a0 = coefs[0];      // c0
0061       double a1 = coefs[1] - a0; // c1 - c0
0062       return a1 * xScaled + a0;
0063
0064    } else if (degree == 2) {
0065
0066       double a0 = coefs[0];            // c0
0067       double a1 = 2 * (coefs[1] - a0); // 2 * (c1 - c0)
0068       double a2 = coefs[2] - a1 - a0;  // c0 - 2 * c1 + c2
0069       return (a2 * xScaled + a1) * xScaled + a0;
0070    }
0071
0072    double t = xScaled;
0073    double s = 1. - xScaled;
0074
0075    double result = coefs[0] * s;
0076    for (int i = 1; i < degree; i++) {
0077       result = (result + t * binomial(degree, i) * coefs[i]) * s;
0078       t *= xScaled;
0079    }
0080    result += t * coefs[degree];
0081
0082    return result;
0083 }
0084
0085 /// @brief Function to evaluate an un-normalized RooGaussian.
0086 inline double gaussian(double x, double mean, double sigma)
0087 {
0088    const double arg = x - mean;
0089    const double sig = sigma;
0090    return std::exp(-0.5 * arg * arg / (sig * sig));
0091 }
0092
0093 inline double product(double const *factors, std::size_t nFactors)
0094 {
0095    double out = 1.0;
0096    for (std::size_t i = 0; i < nFactors; ++i) {
0097       out *= factors[i];
0098    }
0099    return out;
0100 }
0101
0102 // RooRatio evaluate function.
0103 inline double ratio(double numerator, double denominator)
0104 {
0105    return numerator / denominator;
0106 }
0107
0108 inline double bifurGauss(double x, double mean, double sigmaL, double sigmaR)
0109 {
0110    // Note: this simplification does not work with Clad as of v1.1!
0111    // return gaussian(x, mean, x < mean ? sigmaL : sigmaR);
0112    if (x < mean)
0113       return gaussian(x, mean, sigmaL);
0114    return gaussian(x, mean, sigmaR);
0115 }
0116
0117 inline double efficiency(double effFuncVal, int catIndex, int sigCatIndex)
0118 {
0119    // Truncate efficiency function in range 0.0-1.0
0120    effFuncVal = std::clamp(effFuncVal, 0.0, 1.0);
0121
0122    if (catIndex == sigCatIndex)
0123       return effFuncVal; // Accept case
0124    else
0125       return 1 - effFuncVal; // Reject case
0126 }
0127
0128 /// In pdfMode, a coefficient for the constant term of 1.0 is implied if lowestOrder > 0.
0129 template <bool pdfMode = false>
0130 inline double polynomial(double const *coeffs, int nCoeffs, int lowestOrder, double x)
0131 {
0132    double retVal = coeffs[nCoeffs - 1];
0133    for (int i = nCoeffs - 2; i >= 0; i--)
0134       retVal = coeffs[i] + x * retVal;
0135    retVal = retVal * std::pow(x, lowestOrder);
0136    return retVal + (pdfMode && lowestOrder > 0 ? 1.0 : 0.0);
0137 }
0138
0139 inline double chebychev(double *coeffs, unsigned int nCoeffs, double x_in, double xMin, double xMax)
0140 {
0141    // transform to range [-1, +1]
0142    const double xPrime = (x_in - 0.5 * (xMax + xMin)) / (0.5 * (xMax - xMin));
0143
0144    // extract current values of coefficients
0145    double sum = 1.;
0146    if (nCoeffs > 0) {
0147       double curr = xPrime;
0148       double twox = 2 * xPrime;
0149       double last = 1;
0150       double newval = twox * curr - last;
0151       last = curr;
0152       curr = newval;
0153       for (unsigned int i = 0; nCoeffs != i; ++i) {
0154          sum += last * coeffs[i];
0155          newval = twox * curr - last;
0156          last = curr;
0157          curr = newval;
0158       }
0159    }
0160    return sum;
0161 }
0162
0163 inline double constraintSum(double const *comp, unsigned int compSize)
0164 {
0165    double sum = 0;
0166    for (unsigned int i = 0; i < compSize; i++) {
0167       sum -= std::log(comp[i]);
0168    }
0169    return sum;
0170 }
0171
0172 inline unsigned int getUniformBinning(double low, double high, double val, unsigned int numBins)
0173 {
0174    double binWidth = (high - low) / numBins;
0175    return val >= high ? numBins - 1 : std::abs((val - low) / binWidth);
0176 }
0177
0178 inline double interpolate1d(double low, double high, double val, unsigned int numBins, double const* vals)
0179 {
0180    double binWidth = (high - low) / numBins;
0181    int idx = val >= high ? numBins - 1 : std::abs((val - low) / binWidth);
0182
0183    // interpolation
0184    double central = low + (idx + 0.5) * binWidth;
0185    if (val > low + 0.5 * binWidth && val < high - 0.5 * binWidth) {
0186       double slope;
0187       if (val < central) {
0188           slope = vals[idx] - vals[idx - 1];
0189       } else {
0190           slope = vals[idx + 1] - vals[idx];
0191       }
0192       return vals[idx] + slope * (val - central) / binWidth;
0193    }
0194
0195    return vals[idx];
0196 }
0197
0198 inline double poisson(double x, double par)
0199 {
0200    if (par < 0)
0201       return TMath::QuietNaN();
0202
0203    if (x < 0) {
0204       return 0;
0205    } else if (x == 0.0) {
0206       return std::exp(-par);
0207    } else {
0208       double out = x * std::log(par) - TMath::LnGamma(x + 1.) - par;
0209       return std::exp(out);
0210    }
0211 }
0212
0213 inline double flexibleInterpSingle(unsigned int code, double low, double high, double boundary, double nominal,
0214                                    double paramVal, double res)
0215 {
0216    if (code == 0) {
0217       // piece-wise linear
0218       if (paramVal > 0) {
0219          return paramVal * (high - nominal);
0220       } else {
0221          return paramVal * (nominal - low);
0222       }
0223    } else if (code == 1) {
0224       // piece-wise log
0225       if (paramVal >= 0) {
0226          return res * (std::pow(high / nominal, +paramVal) - 1);
0227       } else {
0228          return res * (std::pow(low / nominal, -paramVal) - 1);
0229       }
0230    } else if (code == 2) {
0231       // parabolic with linear
0232       double a = 0.5 * (high + low) - nominal;
0233       double b = 0.5 * (high - low);
0234       double c = 0;
0235       if (paramVal > 1) {
0236          return (2 * a + b) * (paramVal - 1) + high - nominal;
0237       } else if (paramVal < -1) {
0238          return -1 * (2 * a - b) * (paramVal + 1) + low - nominal;
0239       } else {
0240          return a * std::pow(paramVal, 2) + b * paramVal + c;
0241       }
0242    } else if (code == 3) {
0243       // parabolic version of log-normal
0244       double a = 0.5 * (high + low) - nominal;
0245       double b = 0.5 * (high - low);
0246       double c = 0;
0247       if (paramVal > 1) {
0248          return (2 * a + b) * (paramVal - 1) + high - nominal;
0249       } else if (paramVal < -1) {
0250          return -1 * (2 * a - b) * (paramVal + 1) + low - nominal;
0251       } else {
0252          return a * std::pow(paramVal, 2) + b * paramVal + c;
0253       }
0254    } else if (code == 4) {
0255       double x = paramVal;
0256       if (x >= boundary) {
0257          return x * (high - nominal);
0258       } else if (x <= -boundary) {
0259          return x * (nominal - low);
0260       }
0261
0262       // interpolate 6th degree
0263       double t = x / boundary;
0264       double eps_plus = high - nominal;
0265       double eps_minus = nominal - low;
0266       double S = 0.5 * (eps_plus + eps_minus);
0267       double A = 0.0625 * (eps_plus - eps_minus);
0268
0269       return x * (S + t * A * (15 + t * t * (-10 + t * t * 3)));
0270    } else if (code == 5) {
0271       double x = paramVal;
0272       double mod = 1.0;
0273       if (x >= boundary) {
0274          mod = std::pow(high / nominal, +paramVal);
0275       } else if (x <= -boundary) {
0276          mod = std::pow(low / nominal, -paramVal);
0277       } else {
0278          // interpolate 6th degree exp
0279          double x0 = boundary;
0280
0281          high /= nominal;
0282          low /= nominal;
0283
0284          // GHL: Swagato's suggestions
0285          double powUp = std::pow(high, x0);
0286          double powDown = std::pow(low, x0);
0287          double logHi = std::log(high);
0288          double logLo = std::log(low);
0289          double powUpLog = high <= 0.0 ? 0.0 : powUp * logHi;
0290          double powDownLog = low <= 0.0 ? 0.0 : -powDown * logLo;
0291          double powUpLog2 = high <= 0.0 ? 0.0 : powUpLog * logHi;
0292          double powDownLog2 = low <= 0.0 ? 0.0 : -powDownLog * logLo;
0293
0294          double S0 = 0.5 * (powUp + powDown);
0295          double A0 = 0.5 * (powUp - powDown);
0296          double S1 = 0.5 * (powUpLog + powDownLog);
0297          double A1 = 0.5 * (powUpLog - powDownLog);
0298          double S2 = 0.5 * (powUpLog2 + powDownLog2);
0299          double A2 = 0.5 * (powUpLog2 - powDownLog2);
0300
0301          // fcns+der+2nd_der are eq at bd
0302
0303          double a = 1. / (8 * x0) * (15 * A0 - 7 * x0 * S1 + x0 * x0 * A2);
0304          double b = 1. / (8 * x0 * x0) * (-24 + 24 * S0 - 9 * x0 * A1 + x0 * x0 * S2);
0305          double c = 1. / (4 * std::pow(x0, 3)) * (-5 * A0 + 5 * x0 * S1 - x0 * x0 * A2);
0306          double d = 1. / (4 * std::pow(x0, 4)) * (12 - 12 * S0 + 7 * x0 * A1 - x0 * x0 * S2);
0307          double e = 1. / (8 * std::pow(x0, 5)) * (+3 * A0 - 3 * x0 * S1 + x0 * x0 * A2);
0308          double f = 1. / (8 * std::pow(x0, 6)) * (-8 + 8 * S0 - 5 * x0 * A1 + x0 * x0 * S2);
0309
0310          // evaluate the 6-th degree polynomial using Horner's method
0311          double value = 1. + x * (a + x * (b + x * (c + x * (d + x * (e + x * f)))));
0312          mod = value;
0313       }
0314       return res * (mod - 1.0);
0315    }
0316
0317    return 0.0;
0318 }
0319
0320 inline double flexibleInterp(unsigned int code, double const *params, unsigned int n, double const *low,
0321                              double const *high, double boundary, double nominal, int doCutoff)
0322 {
0323    double total = nominal;
0324    for (std::size_t i = 0; i < n; ++i) {
0325       total += flexibleInterpSingle(code, low[i], high[i], boundary, nominal, params[i], total);
0326    }
0327
0328    return doCutoff && total <= 0 ? TMath::Limits<double>::Min() : total;
0329 }
0330
0331 inline double landau(double x, double mu, double sigma)
0332 {
0333    if (sigma <= 0.)
0334       return 0.;
0335    return ROOT::Math::landau_pdf((x - mu) / sigma);
0336 }
0337
0338 inline double logNormal(double x, double k, double m0)
0339 {
0340    return ROOT::Math::lognormal_pdf(x, std::log(m0), std::abs(std::log(k)));
0341 }
0342
0343 inline double logNormalStandard(double x, double sigma, double mu)
0344 {
0345    return ROOT::Math::lognormal_pdf(x, mu, std::abs(sigma));
0346 }
0347
0348 inline double effProd(double eff, double pdf)
0349 {
0350    return eff * pdf;
0351 }
0352
0353 inline double nll(double pdf, double weight, int binnedL, int doBinOffset)
0354 {
0355    if (binnedL) {
0356       // Special handling of this case since std::log(Poisson(0,0)=0 but can't be
0357       // calculated with usual log-formula since std::log(mu)=0. No update of result
0358       // is required since term=0.
0359       if (std::abs(pdf) < 1e-10 && std::abs(weight) < 1e-10) {
0360          return 0.0;
0361       }
0362       if (doBinOffset) {
0363          return pdf - weight - weight * (std::log(pdf) - std::log(weight));
0364       }
0365       return pdf - weight * std::log(pdf) + TMath::LnGamma(weight + 1);
0366    } else {
0367       return -weight * std::log(pdf);
0368    }
0369 }
0370
0371 inline double recursiveFraction(double *a, unsigned int n)
0372 {
0373    double prod = a[0];
0374
0375    for (unsigned int i = 1; i < n; ++i) {
0376       prod *= 1.0 - a[i];
0377    }
0378
0379    return prod;
0380 }
0381
0382 inline double cbShape(double m, double m0, double sigma, double alpha, double n)
0383 {
0384    double t = (m - m0) / sigma;
0385    if (alpha < 0)
0386       t = -t;
0387
0388    double absAlpha = std::abs((double)alpha);
0389
0390    if (t >= -absAlpha) {
0391       return std::exp(-0.5 * t * t);
0392    } else {
0393       double a = std::pow(n / absAlpha, n) * std::exp(-0.5 * absAlpha * absAlpha);
0394       double b = n / absAlpha - absAlpha;
0395
0396       return a / std::pow(b - t, n);
0397    }
0398 }
0399
0400 // For RooCBShape
0401 inline double approxErf(double arg)
0402 {
0403    if (arg > 5.0)
0404       return 1.0;
0405    if (arg < -5.0)
0406       return -1.0;
0407
0408    return TMath::Erf(arg);
0409 }
0410
0411 /// @brief Function to calculate the integral of an un-normalized RooGaussian over x. To calculate the integral over
0412 /// mean, just interchange the respective values of x and mean.
0413 /// @param xMin Minimum value of variable to integrate wrt.
0414 /// @param xMax Maximum value of of variable to integrate wrt.
0415 /// @param mean Mean.
0416 /// @param sigma Sigma.
0417 /// @return The integral of an un-normalized RooGaussian over the value in x.
0418 inline double gaussianIntegral(double xMin, double xMax, double mean, double sigma)
0419 {
0420    // The normalisation constant 1./sqrt(2*pi*sigma^2) is left out in evaluate().
0421    // Therefore, the integral is scaled up by that amount to make RooFit normalise
0422    // correctly.
0423    double resultScale = 0.5 * std::sqrt(TMath::TwoPi()) * sigma;
0424
0425    // Here everything is scaled and shifted into a standard normal distribution:
0426    double xscale = TMath::Sqrt2() * sigma;
0427    double scaledMin = 0.;
0428    double scaledMax = 0.;
0429    scaledMin = (xMin - mean) / xscale;
0430    scaledMax = (xMax - mean) / xscale;
0431
0432    // Here we go for maximum precision: We compute all integrals in the UPPER
0433    // tail of the Gaussian, because erfc has the highest precision there.
0434    // Therefore, the different cases for range limits in the negative hemisphere are mapped onto
0435    // the equivalent points in the upper hemisphere using erfc(-x) = 2. - erfc(x)
0436    double ecmin = TMath::Erfc(std::abs(scaledMin));
0437    double ecmax = TMath::Erfc(std::abs(scaledMax));
0438
0439    double cond = 0.0;
0440    // Don't put this "prd" inside the "if" because clad will not be able to differentiate the code correctly (as of
0441    // v1.1)!
0442    double prd = scaledMin * scaledMax;
0443    if (prd < 0.0) {
0444       cond = 2.0 - (ecmin + ecmax);
0445    } else if (scaledMax <= 0.0) {
0446       cond = ecmax - ecmin;
0447    } else {
0448       cond = ecmin - ecmax;
0449    }
0450    return resultScale * cond;
0451 }
0452
0453 inline double bifurGaussIntegral(double xMin, double xMax, double mean, double sigmaL, double sigmaR)
0454 {
0455    const double xscaleL = TMath::Sqrt2() * sigmaL;
0456    const double xscaleR = TMath::Sqrt2() * sigmaR;
0457
0458    const double resultScale = 0.5 * std::sqrt(TMath::TwoPi());
0459
0460    if (xMax < mean) {
0461       return resultScale * (sigmaL * (TMath::Erf((xMax - mean) / xscaleL) - TMath::Erf((xMin - mean) / xscaleL)));
0462    } else if (xMin > mean) {
0463       return resultScale * (sigmaR * (TMath::Erf((xMax - mean) / xscaleR) - TMath::Erf((xMin - mean) / xscaleR)));
0464    } else {
0465       return resultScale *
0466              (sigmaR * TMath::Erf((xMax - mean) / xscaleR) - sigmaL * TMath::Erf((xMin - mean) / xscaleL));
0467    }
0468 }
0469
0470 inline double exponentialIntegral(double xMin, double xMax, double constant)
0471 {
0472    if (constant == 0.0) {
0473       return xMax - xMin;
0474    }
0475
0476    return (std::exp(constant * xMax) - std::exp(constant * xMin)) / constant;
0477 }
0478
0479 /// In pdfMode, a coefficient for the constant term of 1.0 is implied if lowestOrder > 0.
0480 template <bool pdfMode = false>
0481 inline double polynomialIntegral(double const *coeffs, int nCoeffs, int lowestOrder, double xMin, double xMax)
0482 {
0483    int denom = lowestOrder + nCoeffs;
0484    double min = coeffs[nCoeffs - 1] / double(denom);
0485    double max = coeffs[nCoeffs - 1] / double(denom);
0486
0487    for (int i = nCoeffs - 2; i >= 0; i--) {
0488       denom--;
0489       min = (coeffs[i] / double(denom)) + xMin * min;
0490       max = (coeffs[i] / double(denom)) + xMax * max;
0491    }
0492
0493    max = max * std::pow(xMax, 1 + lowestOrder);
0494    min = min * std::pow(xMin, 1 + lowestOrder);
0495
0496    return max - min + (pdfMode && lowestOrder > 0.0 ? xMax - xMin : 0.0);
0497 }
0498
0499 /// use fast FMA if available, fall back to normal arithmetic if not
0500 inline double fast_fma(double x, double y, double z) noexcept
0501 {
0502 #if defined(FP_FAST_FMA) // check if std::fma has fast hardware implementation
0503    return std::fma(x, y, z);
0504 #else // defined(FP_FAST_FMA)
0505    // std::fma might be slow, so use a more pedestrian implementation
0506 #if defined(__clang__)
0507 #pragma STDC FP_CONTRACT ON // hint clang that using an FMA is okay here
0508 #endif                      // defined(__clang__)
0509    return (x * y) + z;
0510 #endif                      // defined(FP_FAST_FMA)
0511 }
0512
0513 inline double chebychevIntegral(double const *coeffs, unsigned int nCoeffs, double xMin, double xMax, double xMinFull,
0514                                 double xMaxFull)
0515 {
0516    const double halfrange = .5 * (xMax - xMin);
0517    const double mid = .5 * (xMax + xMin);
0518
0519    // the full range of the function is mapped to the normalised [-1, 1] range
0520    const double b = (xMaxFull - mid) / halfrange;
0521    const double a = (xMinFull - mid) / halfrange;
0522
0523    // coefficient for integral(T_0(x)) is 1 (implicit), integrate by hand
0524    // T_0(x) and T_1(x), and use for n > 1: integral(T_n(x) dx) =
0525    // (T_n+1(x) / (n + 1) - T_n-1(x) / (n - 1)) / 2
0526    double sum = b - a; // integrate T_0(x) by hand
0527
0528    const unsigned int iend = nCoeffs;
0529    if (iend > 0) {
0530       {
0531          // integrate T_1(x) by hand...
0532          const double c = coeffs[0];
0533          sum = fast_fma(0.5 * (b + a) * (b - a), c, sum);
0534       }
0535       if (1 < iend) {
0536          double bcurr = b;
0537          double btwox = 2 * b;
0538          double blast = 1;
0539
0540          double acurr = a;
0541          double atwox = 2 * a;
0542          double alast = 1;
0543
0544          double newval = atwox * acurr - alast;
0545          alast = acurr;
0546          acurr = newval;
0547
0548          newval = btwox * bcurr - blast;
0549          blast = bcurr;
0550          bcurr = newval;
0551          double nminus1 = 1.;
0552          for (unsigned int i = 1; iend != i; ++i) {
0553             // integrate using recursion relation
0554             const double c = coeffs[i];
0555             const double term2 = (blast - alast) / nminus1;
0556
0557             newval = atwox * acurr - alast;
0558             alast = acurr;
0559             acurr = newval;
0560
0561             newval = btwox * bcurr - blast;
0562             blast = bcurr;
0563             bcurr = newval;
0564
0565             ++nminus1;
0566             const double term1 = (bcurr - acurr) / (nminus1 + 1.);
0567             const double intTn = 0.5 * (term1 - term2);
0568             sum = fast_fma(intTn, c, sum);
0569          }
0570       }
0571    }
0572
0573    // take care to multiply with the right factor to account for the mapping to
0574    // normalised range [-1, 1]
0575    return halfrange * sum;
0576 }
0577
0578 // The last param should be of type bool but it is not as that causes some issues with Cling for some reason...
0579 inline double
0580 poissonIntegral(int code, double mu, double x, double integrandMin, double integrandMax, unsigned int protectNegative)
0581 {
0582    if (protectNegative && mu < 0.0) {
0583       return std::exp(-2.0 * mu); // make it fall quickly
0584    }
0585
0586    if (code == 1) {
0587       // Implement integral over x as summation. Add special handling in case
0588       // range boundaries are not on integer values of x
0589       integrandMin = std::max(0., integrandMin);
0590
0591       if (integrandMax < 0. || integrandMax < integrandMin) {
0592          return 0;
0593       }
0594       const double delta = 100.0 * std::sqrt(mu);
0595       // If the limits are more than many standard deviations away from the mean,
0596       // we might as well return the integral of the full Poisson distribution to
0597       // save computing time.
0598       if (integrandMin < std::max(mu - delta, 0.0) && integrandMax > mu + delta) {
0599          return 1.;
0600       }
0601
0602       // The range as integers. ixMin is included, ixMax outside.
0603       const unsigned int ixMin = integrandMin;
0604       const unsigned int ixMax = std::min(integrandMax + 1, (double)std::numeric_limits<unsigned int>::max());
0605
0606       // Sum from 0 to just before the bin outside of the range.
0607       if (ixMin == 0) {
0608          return ROOT::Math::gamma_cdf_c(mu, ixMax, 1);
0609       } else {
0610          // If necessary, subtract from 0 to the beginning of the range
0611          if (ixMin <= mu) {
0612             return ROOT::Math::gamma_cdf_c(mu, ixMax, 1) - ROOT::Math::gamma_cdf_c(mu, ixMin, 1);
0613          } else {
0614             // Avoid catastrophic cancellation in the high tails:
0615             return ROOT::Math::gamma_cdf(mu, ixMin, 1) - ROOT::Math::gamma_cdf(mu, ixMax, 1);
0616          }
0617       }
0618    }
0619
0620    // the integral with respect to the mean is the integral of a gamma distribution
0621    // negative ix does not need protection (gamma returns 0.0)
0622    const double ix = 1 + x;
0623
0624    return ROOT::Math::gamma_cdf(integrandMax, ix, 1.0) - ROOT::Math::gamma_cdf(integrandMin, ix, 1.0);
0625 }
0626
0627 inline double logNormalIntegral(double xMin, double xMax, double m0, double k)
0628 {
0629    const double root2 = std::sqrt(2.);
0630
0631    double ln_k = std::abs(std::log(k));
0632    double ret =
0633       0.5 * (TMath::Erf(std::log(xMax / m0) / (root2 * ln_k)) - TMath::Erf(std::log(xMin / m0) / (root2 * ln_k)));
0634
0635    return ret;
0636 }
0637
0638 inline double logNormalIntegralStandard(double xMin, double xMax, double mu, double sigma)
0639 {
0640    const double root2 = std::sqrt(2.);
0641
0642    double ln_k = std::abs(sigma);
0643    double ret =
0644       0.5 * (TMath::Erf((std::log(xMax) - mu) / (root2 * ln_k)) - TMath::Erf((std::log(xMin) - mu) / (root2 * ln_k)));
0645
0646    return ret;
0647 }
0648
0649 inline double cbShapeIntegral(double mMin, double mMax, double m0, double sigma, double alpha, double n)
0650 {
0651    const double sqrtPiOver2 = 1.2533141373;
0652    const double sqrt2 = 1.4142135624;
0653
0654    double result = 0.0;
0655    bool useLog = false;
0656
0657    if (std::abs(n - 1.0) < 1.0e-05)
0658       useLog = true;
0659
0660    double sig = std::abs(sigma);
0661
0662    double tmin = (mMin - m0) / sig;
0663    double tmax = (mMax - m0) / sig;
0664
0665    if (alpha < 0) {
0666       double tmp = tmin;
0667       tmin = -tmax;
0668       tmax = -tmp;
0669    }
0670
0671    double absAlpha = std::abs(alpha);
0672
0673    if (tmin >= -absAlpha) {
0674       result += sig * sqrtPiOver2 * (approxErf(tmax / sqrt2) - approxErf(tmin / sqrt2));
0675    } else if (tmax <= -absAlpha) {
0676       double a = std::pow(n / absAlpha, n) * std::exp(-0.5 * absAlpha * absAlpha);
0677       double b = n / absAlpha - absAlpha;
0678
0679       if (useLog) {
0680          result += a * sig * (std::log(b - tmin) - std::log(b - tmax));
0681       } else {
0682          result += a * sig / (1.0 - n) * (1.0 / (std::pow(b - tmin, n - 1.0)) - 1.0 / (std::pow(b - tmax, n - 1.0)));
0683       }
0684    } else {
0685       double a = std::pow(n / absAlpha, n) * std::exp(-0.5 * absAlpha * absAlpha);
0686       double b = n / absAlpha - absAlpha;
0687
0688       double term1 = 0.0;
0689       if (useLog) {
0690          term1 = a * sig * (std::log(b - tmin) - std::log(n / absAlpha));
0691       } else {
0692          term1 = a * sig / (1.0 - n) * (1.0 / (std::pow(b - tmin, n - 1.0)) - 1.0 / (std::pow(n / absAlpha, n - 1.0)));
0693       }
0694
0695       double term2 = sig * sqrtPiOver2 * (approxErf(tmax / sqrt2) - approxErf(-absAlpha / sqrt2));
0696
0697       result += term1 + term2;
0698    }
0699
0700    if (result == 0)
0701       return 1.E-300;
0702    return result;
0703 }
0704
0705 inline double bernsteinIntegral(double xlo, double xhi, double xmin, double xmax, double *coefs, int nCoefs)
0706 {
0707    double xloScaled = (xlo - xmin) / (xmax - xmin);
0708    double xhiScaled = (xhi - xmin) / (xmax - xmin);
0709
0710    int degree = nCoefs - 1; // n+1 polys of degree n
0711    double norm = 0.;
0712
0713    for (int i = 0; i <= degree; ++i) {
0714       // for each of the i Bernstein basis polynomials
0715       // represent it in the 'power basis' (the naive polynomial basis)
0716       // where the integral is straight forward.
0717       double temp = 0.;
0718       for (int j = i; j <= degree; ++j) { // power basis≈ß
0719          double binCoefs = binomial(degree, j) * binomial(j, i);
0720          double oneOverJPlusOne = 1. / (j + 1.);
0721          double powDiff = std::pow(xhiScaled, j + 1.) - std::pow(xloScaled, j + 1.);
0722          temp += std::pow(-1., j - i) * binCoefs * powDiff * oneOverJPlusOne;
0723       }
0724       temp *= coefs[i]; // include coeff
0725       norm += temp;     // add this basis's contribution to total
0726    }
0727
0728    return norm * (xmax - xmin);
0729 }
0730
0731 } // namespace MathFuncs
0732
0733 } // namespace Detail
0734
0735 } // namespace RooFit
0736
0737 #endif