14#ifndef RooFit_Detail_MathFuncs_h
15#define RooFit_Detail_MathFuncs_h
33 if (
n < 0 || k < 0 ||
n < k)
38 int k1 = std::min(k,
n - k);
41 for (
double i = k1; i > 1.; --i) {
51 int degree = nCoefs - 1;
56 }
else if (degree == 0) {
58 }
else if (degree == 1) {
61 double a1 = coefs[1] - a0;
62 return a1 * xScaled + a0;
64 }
else if (degree == 2) {
67 double a1 = 2 * (coefs[1] - a0);
68 double a2 = coefs[2] - a1 - a0;
69 return (a2 * xScaled + a1) * xScaled + a0;
73 double s = 1. - xScaled;
75 double result = coefs[0] * s;
76 for (
int i = 1; i < degree; i++) {
80 result += t * coefs[degree];
88 const double arg =
x - mean;
89 const double sig =
sigma;
90 return std::exp(-0.5 * arg * arg / (sig * sig));
93inline double product(
double const *factors, std::size_t nFactors)
96 for (std::size_t i = 0; i < nFactors; ++i) {
103inline double ratio(
double numerator,
double denominator)
105 return numerator / denominator;
108inline double bifurGauss(
double x,
double mean,
double sigmaL,
double sigmaR)
117inline double efficiency(
double effFuncVal,
int catIndex,
int sigCatIndex)
120 effFuncVal = std::clamp(effFuncVal, 0.0, 1.0);
122 if (catIndex == sigCatIndex)
125 return 1 - effFuncVal;
129template <
bool pdfMode = false>
130inline double polynomial(
double const *coeffs,
int nCoeffs,
int lowestOrder,
double x)
132 double retVal = coeffs[nCoeffs - 1];
133 for (
int i = nCoeffs - 2; i >= 0; i--)
134 retVal = coeffs[i] +
x * retVal;
135 retVal = retVal * std::pow(
x, lowestOrder);
136 return retVal + (pdfMode && lowestOrder > 0 ? 1.0 : 0.0);
139inline double chebychev(
double *coeffs,
unsigned int nCoeffs,
double x_in,
double xMin,
double xMax)
142 const double xPrime = (x_in - 0.5 * (xMax + xMin)) / (0.5 * (xMax - xMin));
147 double curr = xPrime;
148 double twox = 2 * xPrime;
150 double newval = twox * curr - last;
153 for (
unsigned int i = 0; nCoeffs != i; ++i) {
154 sum += last * coeffs[i];
155 newval = twox * curr - last;
166 for (
unsigned int i = 0; i < compSize; i++) {
167 sum -= std::log(comp[i]);
174 double binWidth = (high - low) / numBins;
175 return val >= high ? numBins - 1 : std::abs((val - low) / binWidth);
185 }
else if (
x == 0.0) {
186 return std::exp(-par);
189 return std::exp(out);
194 double paramVal,
double res)
199 return paramVal * (high - nominal);
201 return paramVal * (nominal - low);
203 }
else if (code == 1) {
206 return res * (std::pow(high / nominal, +paramVal) - 1);
208 return res * (std::pow(low / nominal, -paramVal) - 1);
210 }
else if (code == 2) {
212 double a = 0.5 * (high + low) - nominal;
213 double b = 0.5 * (high - low);
216 return (2 *
a +
b) * (paramVal - 1) + high - nominal;
217 }
else if (paramVal < -1) {
218 return -1 * (2 *
a -
b) * (paramVal + 1) + low - nominal;
220 return a * std::pow(paramVal, 2) +
b * paramVal +
c;
222 }
else if (code == 3) {
224 double a = 0.5 * (high + low) - nominal;
225 double b = 0.5 * (high - low);
228 return (2 *
a +
b) * (paramVal - 1) + high - nominal;
229 }
else if (paramVal < -1) {
230 return -1 * (2 *
a -
b) * (paramVal + 1) + low - nominal;
232 return a * std::pow(paramVal, 2) +
b * paramVal +
c;
234 }
else if (code == 4) {
237 return x * (high - nominal);
238 }
else if (
x <= -boundary) {
239 return x * (nominal - low);
243 double t =
x / boundary;
244 double eps_plus = high - nominal;
245 double eps_minus = nominal - low;
246 double S = 0.5 * (eps_plus + eps_minus);
247 double A = 0.0625 * (eps_plus - eps_minus);
249 return x * (S + t * A * (15 + t * t * (-10 + t * t * 3)));
250 }
else if (code == 5) {
254 mod = std::pow(high / nominal, +paramVal);
255 }
else if (
x <= -boundary) {
256 mod = std::pow(low / nominal, -paramVal);
259 double x0 = boundary;
262 double powUp = std::pow(high / nominal, x0);
263 double powDown = std::pow(low / nominal, x0);
264 double logHi = std::log(high);
265 double logLo = std::log(low);
266 double powUpLog = high <= 0.0 ? 0.0 : powUp * logHi;
267 double powDownLog = low <= 0.0 ? 0.0 : -powDown * logLo;
268 double powUpLog2 = high <= 0.0 ? 0.0 : powUpLog * logHi;
269 double powDownLog2 = low <= 0.0 ? 0.0 : -powDownLog * logLo;
271 double S0 = 0.5 * (powUp + powDown);
272 double A0 = 0.5 * (powUp - powDown);
273 double S1 = 0.5 * (powUpLog + powDownLog);
274 double A1 = 0.5 * (powUpLog - powDownLog);
275 double S2 = 0.5 * (powUpLog2 + powDownLog2);
276 double A2 = 0.5 * (powUpLog2 - powDownLog2);
280 double a = 1. / (8 * x0) * (15 * A0 - 7 * x0 *
S1 + x0 * x0 * A2);
281 double b = 1. / (8 * x0 * x0) * (-24 + 24 *
S0 - 9 * x0 * A1 + x0 * x0 * S2);
282 double c = 1. / (4 * std::pow(x0, 3)) * (-5 * A0 + 5 * x0 *
S1 - x0 * x0 * A2);
283 double d = 1. / (4 * std::pow(x0, 4)) * (12 - 12 *
S0 + 7 * x0 * A1 - x0 * x0 * S2);
284 double e = 1. / (8 * std::pow(x0, 5)) * (+3 * A0 - 3 * x0 *
S1 + x0 * x0 * A2);
285 double f = 1. / (8 * std::pow(x0, 6)) * (-8 + 8 *
S0 - 5 * x0 * A1 + x0 * x0 * S2);
291 return res * (mod - 1.0);
297inline double flexibleInterp(
unsigned int code,
double *params,
unsigned int n,
double *low,
double *high,
298 double boundary,
double nominal,
int doCutoff)
300 double total = nominal;
301 for (std::size_t i = 0; i <
n; ++i) {
305 return doCutoff && total <= 0 ? TMath::Limits<double>::Min() :
total;
330inline double nll(
double pdf,
double weight,
int binnedL,
int doBinOffset)
336 if (std::abs(pdf) < 1
e-10 && std::abs(weight) < 1
e-10) {
340 return pdf - weight - weight * (std::log(pdf) - std::log(weight));
344 return -weight * std::log(pdf);
352 for (
unsigned int i = 1; i <
n; ++i) {
361 double t = (
m - m0) /
sigma;
365 double absAlpha = std::abs((
double)alpha);
367 if (t >= -absAlpha) {
368 return std::exp(-0.5 * t * t);
370 double a = std::pow(
n / absAlpha,
n) * std::exp(-0.5 * absAlpha * absAlpha);
371 double b =
n / absAlpha - absAlpha;
373 return a / std::pow(
b - t,
n);
404 double scaledMin = 0.;
405 double scaledMax = 0.;
406 scaledMin = (xMin - mean) / xscale;
407 scaledMax = (xMax - mean) / xscale;
419 double prd = scaledMin * scaledMax;
421 cond = 2.0 - (ecmin + ecmax);
422 }
else if (scaledMax <= 0.0) {
423 cond = ecmax - ecmin;
425 cond = ecmin - ecmax;
427 return resultScale * cond;
435 const double resultScale = 0.5 * std::sqrt(
TMath::TwoPi());
438 return resultScale * (sigmaL * (
TMath::Erf((xMax - mean) / xscaleL) -
TMath::Erf((xMin - mean) / xscaleL)));
439 }
else if (xMin > mean) {
440 return resultScale * (sigmaR * (
TMath::Erf((xMax - mean) / xscaleR) -
TMath::Erf((xMin - mean) / xscaleR)));
449 if (constant == 0.0) {
453 return (std::exp(constant * xMax) - std::exp(constant * xMin)) / constant;
457template <
bool pdfMode = false>
458inline double polynomialIntegral(
double const *coeffs,
int nCoeffs,
int lowestOrder,
double xMin,
double xMax)
460 int denom = lowestOrder + nCoeffs;
461 double min = coeffs[nCoeffs - 1] /
double(denom);
462 double max = coeffs[nCoeffs - 1] /
double(denom);
464 for (
int i = nCoeffs - 2; i >= 0; i--) {
466 min = (coeffs[i] /
double(denom)) + xMin * min;
467 max = (coeffs[i] /
double(denom)) + xMax * max;
470 max = max * std::pow(xMax, 1 + lowestOrder);
471 min = min * std::pow(xMin, 1 + lowestOrder);
473 return max - min + (pdfMode && lowestOrder > 0.0 ? xMax - xMin : 0.0);
479#if defined(FP_FAST_FMA)
480 return std::fma(
x,
y, z);
483#if defined(__clang__)
484#pragma STDC FP_CONTRACT ON
490inline double chebychevIntegral(
double const *coeffs,
unsigned int nCoeffs,
double xMin,
double xMax,
double xMinFull,
493 const double halfrange = .5 * (xMax - xMin);
494 const double mid = .5 * (xMax + xMin);
497 const double b = (xMaxFull - mid) / halfrange;
498 const double a = (xMinFull - mid) / halfrange;
505 const unsigned int iend = nCoeffs;
509 const double c = coeffs[0];
514 double btwox = 2 *
b;
518 double atwox = 2 *
a;
521 double newval = atwox * acurr - alast;
525 newval = btwox * bcurr - blast;
529 for (
unsigned int i = 1; iend != i; ++i) {
531 const double c = coeffs[i];
532 const double term2 = (blast - alast) / nminus1;
534 newval = atwox * acurr - alast;
538 newval = btwox * bcurr - blast;
543 const double term1 = (bcurr - acurr) / (nminus1 + 1.);
544 const double intTn = 0.5 * (term1 - term2);
552 return halfrange *
sum;
557poissonIntegral(
int code,
double mu,
double x,
double integrandMin,
double integrandMax,
unsigned int protectNegative)
559 if (protectNegative && mu < 0.0) {
560 return std::exp(-2.0 * mu);
566 integrandMin = std::max(0., integrandMin);
568 if (integrandMax < 0. || integrandMax < integrandMin) {
571 const double delta = 100.0 * std::sqrt(mu);
575 if (integrandMin < std::max(mu - delta, 0.0) && integrandMax > mu + delta) {
580 const unsigned int ixMin = integrandMin;
581 const unsigned int ixMax = std::min(integrandMax + 1, (
double)std::numeric_limits<unsigned int>::max());
599 const double ix = 1 +
x;
606 const double root2 = std::sqrt(2.);
608 double ln_k = std::abs(std::log(k));
610 0.5 * (
TMath::Erf(std::log(xMax / m0) / (root2 * ln_k)) -
TMath::Erf(std::log(xMin / m0) / (root2 * ln_k)));
617 const double root2 = std::sqrt(2.);
619 double ln_k = std::abs(
sigma);
621 0.5 * (
TMath::Erf((std::log(xMax) - mu) / (root2 * ln_k)) -
TMath::Erf((std::log(xMin) - mu) / (root2 * ln_k)));
628 const double sqrtPiOver2 = 1.2533141373;
629 const double sqrt2 = 1.4142135624;
634 if (std::abs(
n - 1.0) < 1.0e-05)
637 double sig = std::abs(
sigma);
639 double tmin = (mMin - m0) / sig;
640 double tmax = (mMax - m0) / sig;
648 double absAlpha = std::abs(alpha);
650 if (tmin >= -absAlpha) {
652 }
else if (tmax <= -absAlpha) {
653 double a = std::pow(
n / absAlpha,
n) * std::exp(-0.5 * absAlpha * absAlpha);
654 double b =
n / absAlpha - absAlpha;
657 result +=
a * sig * (std::log(
b - tmin) - std::log(
b - tmax));
659 result +=
a * sig / (1.0 -
n) * (1.0 / (std::pow(
b - tmin,
n - 1.0)) - 1.0 / (std::pow(
b - tmax,
n - 1.0)));
662 double a = std::pow(
n / absAlpha,
n) * std::exp(-0.5 * absAlpha * absAlpha);
663 double b =
n / absAlpha - absAlpha;
667 term1 =
a * sig * (std::log(
b - tmin) - std::log(
n / absAlpha));
669 term1 =
a * sig / (1.0 -
n) * (1.0 / (std::pow(
b - tmin,
n - 1.0)) - 1.0 / (std::pow(
n / absAlpha,
n - 1.0)));
672 double term2 = sig * sqrtPiOver2 * (
approxErf(tmax / sqrt2) -
approxErf(-absAlpha / sqrt2));
687 int degree = nCoefs - 1;
690 for (
int i = 0; i <= degree; ++i) {
695 for (
int j = i; j <= degree; ++j) {
697 double oneOverJPlusOne = 1. / (j + 1.);
698 double powDiff = std::pow(xhiScaled, j + 1.) - std::pow(xloScaled, j + 1.);
699 temp += std::pow(-1., j - i) * binCoefs * powDiff * oneOverJPlusOne;
static unsigned int total
Option_t Option_t TPoint TPoint const char GetTextMagnitude GetFillStyle GetLineColor GetLineWidth GetMarkerStyle GetTextAlign GetTextColor GetTextSize void char Point_t Rectangle_t WindowAttributes_t Float_t Float_t Float_t Int_t Int_t UInt_t UInt_t Rectangle_t result
Option_t Option_t TPoint TPoint const char GetTextMagnitude GetFillStyle GetLineColor GetLineWidth GetMarkerStyle GetTextAlign GetTextColor GetTextSize void value
double lognormal_pdf(double x, double m, double s, double x0=0)
Probability density function of the lognormal distribution.
double landau_pdf(double x, double xi=1, double x0=0)
Probability density function of the Landau distribution:
double gamma_cdf_c(double x, double alpha, double theta, double x0=0)
Complement of the cumulative distribution function of the gamma distribution (upper tail).
double gamma_cdf(double x, double alpha, double theta, double x0=0)
Cumulative distribution function of the gamma distribution (lower tail).
double logNormalIntegral(double xMin, double xMax, double m0, double k)
double gaussianIntegral(double xMin, double xMax, double mean, double sigma)
Function to calculate the integral of an un-normalized RooGaussian over x.
double chebychevIntegral(double const *coeffs, unsigned int nCoeffs, double xMin, double xMax, double xMinFull, double xMaxFull)
double bifurGaussIntegral(double xMin, double xMax, double mean, double sigmaL, double sigmaR)
double cbShape(double m, double m0, double sigma, double alpha, double n)
double polynomial(double const *coeffs, int nCoeffs, int lowestOrder, double x)
In pdfMode, a coefficient for the constant term of 1.0 is implied if lowestOrder > 0.
double recursiveFraction(double *a, unsigned int n)
double constraintSum(double const *comp, unsigned int compSize)
double cbShapeIntegral(double mMin, double mMax, double m0, double sigma, double alpha, double n)
double flexibleInterp(unsigned int code, double *params, unsigned int n, double *low, double *high, double boundary, double nominal, int doCutoff)
double fast_fma(double x, double y, double z) noexcept
use fast FMA if available, fall back to normal arithmetic if not
double logNormalIntegralStandard(double xMin, double xMax, double mu, double sigma)
double landau(double x, double mu, double sigma)
double gaussian(double x, double mean, double sigma)
Function to evaluate an un-normalized RooGaussian.
double product(double const *factors, std::size_t nFactors)
double chebychev(double *coeffs, unsigned int nCoeffs, double x_in, double xMin, double xMax)
double poisson(double x, double par)
double binomial(int n, int k)
Calculates the binomial coefficient n over k.
unsigned int getUniformBinning(double low, double high, double val, unsigned int numBins)
double flexibleInterpSingle(unsigned int code, double low, double high, double boundary, double nominal, double paramVal, double res)
double logNormalStandard(double x, double sigma, double mu)
double bifurGauss(double x, double mean, double sigmaL, double sigmaR)
double ratio(double numerator, double denominator)
double polynomialIntegral(double const *coeffs, int nCoeffs, int lowestOrder, double xMin, double xMax)
In pdfMode, a coefficient for the constant term of 1.0 is implied if lowestOrder > 0.
double bernsteinIntegral(double xlo, double xhi, double xmin, double xmax, double *coefs, int nCoefs)
double approxErf(double arg)
double effProd(double eff, double pdf)
double poissonIntegral(int code, double mu, double x, double integrandMin, double integrandMax, unsigned int protectNegative)
double logNormal(double x, double k, double m0)
double nll(double pdf, double weight, int binnedL, int doBinOffset)
double bernstein(double x, double xmin, double xmax, double *coefs, int nCoefs)
The caller needs to make sure that there is at least one coefficient.
double efficiency(double effFuncVal, int catIndex, int sigCatIndex)
double exponentialIntegral(double xMin, double xMax, double constant)
The namespace RooFit contains mostly switches that change the behaviour of functions of PDFs (or othe...
Double_t Erf(Double_t x)
Computation of the error function erf(x).
Double_t QuietNaN()
Returns a quiet NaN as defined by IEEE 754.
constexpr Double_t Sqrt2()
Double_t Erfc(Double_t x)
Computes the complementary error function erfc(x).
Double_t LnGamma(Double_t z)
Computation of ln[gamma(z)] for all z.
Double_t SignalingNaN()
Returns a signaling NaN as defined by IEEE 754](http://en.wikipedia.org/wiki/NaN#Signaling_NaN).
constexpr Double_t TwoPi()
static uint64_t sum(uint64_t i)