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));
94inline double ratio(
double numerator,
double denominator)
96 return numerator / denominator;
99inline double bifurGauss(
double x,
double mean,
double sigmaL,
double sigmaR)
108inline double efficiency(
double effFuncVal,
int catIndex,
int sigCatIndex)
111 effFuncVal = std::clamp(effFuncVal, 0.0, 1.0);
113 if (catIndex == sigCatIndex)
116 return 1 - effFuncVal;
120template <
bool pdfMode = false>
121inline double polynomial(
double const *coeffs,
int nCoeffs,
int lowestOrder,
double x)
123 double retVal = coeffs[nCoeffs - 1];
124 for (
int i = nCoeffs - 2; i >= 0; i--)
125 retVal = coeffs[i] +
x * retVal;
126 retVal = retVal * std::pow(
x, lowestOrder);
127 return retVal + (pdfMode && lowestOrder > 0 ? 1.0 : 0.0);
130inline double chebychev(
double *coeffs,
unsigned int nCoeffs,
double x_in,
double xMin,
double xMax)
133 const double xPrime = (x_in - 0.5 * (xMax + xMin)) / (0.5 * (xMax - xMin));
138 double curr = xPrime;
139 double twox = 2 * xPrime;
141 double newval = twox * curr - last;
144 for (
unsigned int i = 0; nCoeffs != i; ++i) {
145 sum += last * coeffs[i];
146 newval = twox * curr - last;
157 for (
unsigned int i = 0; i < compSize; i++) {
158 sum -= std::log(comp[i]);
165 double binWidth = (high - low) / numBins;
166 return val >= high ? numBins - 1 : std::abs((val - low) / binWidth);
176 }
else if (
x == 0.0) {
177 return std::exp(-par);
180 return std::exp(out);
185 double paramVal,
double res)
190 return paramVal * (high - nominal);
192 return paramVal * (nominal - low);
194 }
else if (code == 1) {
197 return res * (std::pow(high / nominal, +paramVal) - 1);
199 return res * (std::pow(low / nominal, -paramVal) - 1);
201 }
else if (code == 2) {
203 double a = 0.5 * (high + low) - nominal;
204 double b = 0.5 * (high - low);
207 return (2 *
a +
b) * (paramVal - 1) + high - nominal;
208 }
else if (paramVal < -1) {
209 return -1 * (2 *
a -
b) * (paramVal + 1) + low - nominal;
211 return a * std::pow(paramVal, 2) +
b * paramVal +
c;
213 }
else if (code == 3) {
215 double a = 0.5 * (high + low) - nominal;
216 double b = 0.5 * (high - low);
219 return (2 *
a +
b) * (paramVal - 1) + high - nominal;
220 }
else if (paramVal < -1) {
221 return -1 * (2 *
a -
b) * (paramVal + 1) + low - nominal;
223 return a * std::pow(paramVal, 2) +
b * paramVal +
c;
225 }
else if (code == 4) {
228 return x * (high - nominal);
229 }
else if (
x <= -boundary) {
230 return x * (nominal - low);
234 double t =
x / boundary;
235 double eps_plus = high - nominal;
236 double eps_minus = nominal - low;
237 double S = 0.5 * (eps_plus + eps_minus);
238 double A = 0.0625 * (eps_plus - eps_minus);
240 return x * (S + t * A * (15 + t * t * (-10 + t * t * 3)));
241 }
else if (code == 5) {
245 mod = std::pow(high / nominal, +paramVal);
246 }
else if (
x <= -boundary) {
247 mod = std::pow(low / nominal, -paramVal);
250 double x0 = boundary;
253 double powUp = std::pow(high / nominal, x0);
254 double powDown = std::pow(low / nominal, x0);
255 double logHi = std::log(high);
256 double logLo = std::log(low);
257 double powUpLog = high <= 0.0 ? 0.0 : powUp * logHi;
258 double powDownLog = low <= 0.0 ? 0.0 : -powDown * logLo;
259 double powUpLog2 = high <= 0.0 ? 0.0 : powUpLog * logHi;
260 double powDownLog2 = low <= 0.0 ? 0.0 : -powDownLog * logLo;
262 double S0 = 0.5 * (powUp + powDown);
263 double A0 = 0.5 * (powUp - powDown);
264 double S1 = 0.5 * (powUpLog + powDownLog);
265 double A1 = 0.5 * (powUpLog - powDownLog);
266 double S2 = 0.5 * (powUpLog2 + powDownLog2);
267 double A2 = 0.5 * (powUpLog2 - powDownLog2);
271 double a = 1. / (8 * x0) * (15 * A0 - 7 * x0 *
S1 + x0 * x0 * A2);
272 double b = 1. / (8 * x0 * x0) * (-24 + 24 *
S0 - 9 * x0 * A1 + x0 * x0 * S2);
273 double c = 1. / (4 * std::pow(x0, 3)) * (-5 * A0 + 5 * x0 *
S1 - x0 * x0 * A2);
274 double d = 1. / (4 * std::pow(x0, 4)) * (12 - 12 *
S0 + 7 * x0 * A1 - x0 * x0 * S2);
275 double e = 1. / (8 * std::pow(x0, 5)) * (+3 * A0 - 3 * x0 *
S1 + x0 * x0 * A2);
276 double f = 1. / (8 * std::pow(x0, 6)) * (-8 + 8 *
S0 - 5 * x0 * A1 + x0 * x0 * S2);
282 return res * (mod - 1.0);
288inline double flexibleInterp(
unsigned int code,
double *params,
unsigned int n,
double *low,
double *high,
289 double boundary,
double nominal,
int doCutoff)
291 double total = nominal;
292 for (std::size_t i = 0; i <
n; ++i) {
296 return doCutoff && total <= 0 ? TMath::Limits<double>::Min() :
total;
321inline double nll(
double pdf,
double weight,
int binnedL,
int doBinOffset)
327 if (std::abs(pdf) < 1
e-10 && std::abs(weight) < 1
e-10) {
331 return pdf - weight - weight * (std::log(pdf) - std::log(weight));
335 return -weight * std::log(pdf);
343 for (
unsigned int i = 1; i <
n; ++i) {
352 double t = (
m - m0) /
sigma;
356 double absAlpha = std::abs((
double)alpha);
358 if (t >= -absAlpha) {
359 return std::exp(-0.5 * t * t);
361 double a = std::pow(
n / absAlpha,
n) * std::exp(-0.5 * absAlpha * absAlpha);
362 double b =
n / absAlpha - absAlpha;
364 return a / std::pow(
b - t,
n);
395 double scaledMin = 0.;
396 double scaledMax = 0.;
397 scaledMin = (xMin - mean) / xscale;
398 scaledMax = (xMax - mean) / xscale;
410 double prd = scaledMin * scaledMax;
412 cond = 2.0 - (ecmin + ecmax);
413 }
else if (scaledMax <= 0.0) {
414 cond = ecmax - ecmin;
416 cond = ecmin - ecmax;
418 return resultScale * cond;
426 const double resultScale = 0.5 * std::sqrt(
TMath::TwoPi());
429 return resultScale * (sigmaL * (
TMath::Erf((xMax - mean) / xscaleL) -
TMath::Erf((xMin - mean) / xscaleL)));
430 }
else if (xMin > mean) {
431 return resultScale * (sigmaR * (
TMath::Erf((xMax - mean) / xscaleR) -
TMath::Erf((xMin - mean) / xscaleR)));
440 if (constant == 0.0) {
444 return (std::exp(constant * xMax) - std::exp(constant * xMin)) / constant;
448template <
bool pdfMode = false>
449inline double polynomialIntegral(
double const *coeffs,
int nCoeffs,
int lowestOrder,
double xMin,
double xMax)
451 int denom = lowestOrder + nCoeffs;
452 double min = coeffs[nCoeffs - 1] /
double(denom);
453 double max = coeffs[nCoeffs - 1] /
double(denom);
455 for (
int i = nCoeffs - 2; i >= 0; i--) {
457 min = (coeffs[i] /
double(denom)) + xMin * min;
458 max = (coeffs[i] /
double(denom)) + xMax * max;
461 max = max * std::pow(xMax, 1 + lowestOrder);
462 min = min * std::pow(xMin, 1 + lowestOrder);
464 return max - min + (pdfMode && lowestOrder > 0.0 ? xMax - xMin : 0.0);
470#if defined(FP_FAST_FMA)
471 return std::fma(
x,
y, z);
474#if defined(__clang__)
475#pragma STDC FP_CONTRACT ON
481inline double chebychevIntegral(
double const *coeffs,
unsigned int nCoeffs,
double xMin,
double xMax,
double xMinFull,
484 const double halfrange = .5 * (xMax - xMin);
485 const double mid = .5 * (xMax + xMin);
488 const double b = (xMaxFull - mid) / halfrange;
489 const double a = (xMinFull - mid) / halfrange;
496 const unsigned int iend = nCoeffs;
500 const double c = coeffs[0];
505 double btwox = 2 *
b;
509 double atwox = 2 *
a;
512 double newval = atwox * acurr - alast;
516 newval = btwox * bcurr - blast;
520 for (
unsigned int i = 1; iend != i; ++i) {
522 const double c = coeffs[i];
523 const double term2 = (blast - alast) / nminus1;
525 newval = atwox * acurr - alast;
529 newval = btwox * bcurr - blast;
534 const double term1 = (bcurr - acurr) / (nminus1 + 1.);
535 const double intTn = 0.5 * (term1 - term2);
543 return halfrange *
sum;
548poissonIntegral(
int code,
double mu,
double x,
double integrandMin,
double integrandMax,
unsigned int protectNegative)
550 if (protectNegative && mu < 0.0) {
551 return std::exp(-2.0 * mu);
557 integrandMin = std::max(0., integrandMin);
559 if (integrandMax < 0. || integrandMax < integrandMin) {
562 const double delta = 100.0 * std::sqrt(mu);
566 if (integrandMin < std::max(mu - delta, 0.0) && integrandMax > mu + delta) {
571 const unsigned int ixMin = integrandMin;
572 const unsigned int ixMax = std::min(integrandMax + 1, (
double)std::numeric_limits<unsigned int>::max());
590 const double ix = 1 +
x;
597 const double root2 = std::sqrt(2.);
599 double ln_k = std::abs(std::log(k));
601 0.5 * (
TMath::Erf(std::log(xMax / m0) / (root2 * ln_k)) -
TMath::Erf(std::log(xMin / m0) / (root2 * ln_k)));
608 const double root2 = std::sqrt(2.);
610 double ln_k = std::abs(
sigma);
612 0.5 * (
TMath::Erf((std::log(xMax) - mu) / (root2 * ln_k)) -
TMath::Erf((std::log(xMin) - mu) / (root2 * ln_k)));
619 const double sqrtPiOver2 = 1.2533141373;
620 const double sqrt2 = 1.4142135624;
625 if (std::abs(
n - 1.0) < 1.0e-05)
628 double sig = std::abs(
sigma);
630 double tmin = (mMin - m0) / sig;
631 double tmax = (mMax - m0) / sig;
639 double absAlpha = std::abs(alpha);
641 if (tmin >= -absAlpha) {
643 }
else if (tmax <= -absAlpha) {
644 double a = std::pow(
n / absAlpha,
n) * std::exp(-0.5 * absAlpha * absAlpha);
645 double b =
n / absAlpha - absAlpha;
648 result +=
a * sig * (std::log(
b - tmin) - std::log(
b - tmax));
650 result +=
a * sig / (1.0 -
n) * (1.0 / (std::pow(
b - tmin,
n - 1.0)) - 1.0 / (std::pow(
b - tmax,
n - 1.0)));
653 double a = std::pow(
n / absAlpha,
n) * std::exp(-0.5 * absAlpha * absAlpha);
654 double b =
n / absAlpha - absAlpha;
658 term1 =
a * sig * (std::log(
b - tmin) - std::log(
n / absAlpha));
660 term1 =
a * sig / (1.0 -
n) * (1.0 / (std::pow(
b - tmin,
n - 1.0)) - 1.0 / (std::pow(
n / absAlpha,
n - 1.0)));
663 double term2 = sig * sqrtPiOver2 * (
approxErf(tmax / sqrt2) -
approxErf(-absAlpha / sqrt2));
678 int degree = nCoefs - 1;
681 for (
int i = 0; i <= degree; ++i) {
686 for (
int j = i; j <= degree; ++j) {
688 double oneOverJPlusOne = 1. / (j + 1.);
689 double powDiff = std::pow(xhiScaled, j + 1.) - std::pow(xloScaled, j + 1.);
690 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 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)