14#ifndef RooFit_Detail_AnalyticalIntegrals_h
15#define RooFit_Detail_AnalyticalIntegrals_h
26namespace AnalyticalIntegrals {
44 double scaledMin = 0.;
45 double scaledMax = 0.;
46 scaledMin = (xMin - mean) / xscale;
47 scaledMax = (xMax - mean) / xscale;
59 double prd = scaledMin * scaledMax;
61 cond = 2.0 - (ecmin + ecmax);
62 }
else if (scaledMax <= 0.0) {
67 return resultScale * cond;
70inline double bifurGaussIntegral(
double xMin,
double xMax,
double mean,
double sigmaL,
double sigmaR)
75 const double resultScale = 0.5 * std::sqrt(
TMath::TwoPi());
78 return resultScale * (sigmaL * (
TMath::Erf((xMax - mean) / xscaleL) -
TMath::Erf((xMin - mean) / xscaleL)));
79 }
else if (xMin > mean) {
80 return resultScale * (sigmaR * (
TMath::Erf((xMax - mean) / xscaleR) -
TMath::Erf((xMin - mean) / xscaleR)));
82 return resultScale * (sigmaR *
TMath::Erf((xMax - mean) / xscaleR) - sigmaL *
TMath::Erf((xMin - mean) / xscaleL));
88 if (constant == 0.0) {
92 return (std::exp(constant * xMax) - std::exp(constant * xMin)) / constant;
96template <
bool pdfMode = false>
97inline double polynomialIntegral(
double const *coeffs,
int nCoeffs,
int lowestOrder,
double xMin,
double xMax)
99 int denom = lowestOrder + nCoeffs;
100 double min = coeffs[nCoeffs - 1] /
double(denom);
101 double max = coeffs[nCoeffs - 1] /
double(denom);
103 for (
int i = nCoeffs - 2; i >= 0; i--) {
106 max = (coeffs[i] /
double(denom)) + xMax * max;
109 max = max * std::pow(xMax, 1 + lowestOrder);
110 min =
min * std::pow(xMin, 1 + lowestOrder);
112 return max -
min + (pdfMode && lowestOrder > 0.0 ? xMax - xMin : 0.0);
118#if defined(FP_FAST_FMA)
119 return std::fma(
x,
y, z);
122#if defined(__clang__)
123#pragma STDC FP_CONTRACT ON
129inline double chebychevIntegral(
double const *coeffs,
unsigned int nCoeffs,
double xMin,
double xMax,
double xMinFull,
132 const double halfrange = .5 * (xMax - xMin);
133 const double mid = .5 * (xMax + xMin);
136 const double b = (xMaxFull - mid) / halfrange;
137 const double a = (xMinFull - mid) / halfrange;
144 const unsigned int iend = nCoeffs;
148 const double c = coeffs[0];
153 double btwox = 2 *
b;
157 double atwox = 2 *
a;
160 double newval = atwox * acurr - alast;
164 newval = btwox * bcurr - blast;
168 for (
unsigned int i = 1; iend != i; ++i) {
170 const double c = coeffs[i];
171 const double term2 = (blast - alast) / nminus1;
173 newval = atwox * acurr - alast;
177 newval = btwox * bcurr - blast;
182 const double term1 = (bcurr - acurr) / (nminus1 + 1.);
183 const double intTn = 0.5 * (term1 - term2);
191 return halfrange *
sum;
195inline double max(
double x,
double y)
197 return x >=
y ?
x :
y;
200inline double min(
double x,
double y)
202 return x <=
y ?
x :
y;
207poissonIntegral(
int code,
double mu,
double x,
double integrandMin,
double integrandMax,
unsigned int protectNegative)
209 if (protectNegative && mu < 0.0) {
210 return std::exp(-2.0 * mu);
216 integrandMin = max(0, integrandMin);
218 if (integrandMax < 0. || integrandMax < integrandMin) {
221 const double delta = 100.0 * std::sqrt(mu);
225 if (integrandMin < max(mu - delta, 0.0) && integrandMax > mu + delta) {
230 const unsigned int ixMin = integrandMin;
231 const unsigned int ixMax =
min(integrandMax + 1, (
double)std::numeric_limits<unsigned int>::max());
249 const double ix = 1 +
x;
256 const double root2 = std::sqrt(2.);
260 0.5 * (
TMath::Erf(std::log(xMax / m0) / (root2 * ln_k)) -
TMath::Erf(std::log(xMin / m0) / (root2 * ln_k)));
267 const double root2 = std::sqrt(2.);
269 double ln_k = std::abs(
sigma);
271 0.5 * (
TMath::Erf((std::log(xMax) - mu) / (root2 * ln_k)) -
TMath::Erf((std::log(xMin) - mu) / (root2 * ln_k)));
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 chebychevIntegral(double const *coeffs, unsigned int nCoeffs, double xMin, double xMax, double xMinFull, double xMaxFull)
double poissonIntegral(int code, double mu, double x, double integrandMin, double integrandMax, unsigned int protectNegative)
double fast_fma(double x, double y, double z) noexcept
use fast FMA if available, fall back to normal arithmetic if not
double min(double x, double y)
double gaussianIntegral(double xMin, double xMax, double mean, double sigma)
Function to calculate the integral of an un-normalized RooGaussian over x.
double bifurGaussIntegral(double xMin, double xMax, double mean, double sigmaL, double sigmaR)
double exponentialIntegral(double xMin, double xMax, double constant)
double logNormalIntegralStandard(double xMin, double xMax, double mu, double sigma)
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.
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).
constexpr Double_t Sqrt2()
Double_t Erfc(Double_t x)
Computes the complementary error function erfc(x).
Short_t Abs(Short_t d)
Returns the absolute value of parameter Short_t d.
constexpr Double_t TwoPi()
static uint64_t sum(uint64_t i)