ROOT   Reference Guide
Searching...
No Matches
AnalyticalIntegrals.h
Go to the documentation of this file.
1/*
2 * Project: RooFit
3 * Authors:
4 * Jonas Rembser, CERN 2023
5 * Garima Singh, CERN 2023
6 *
7 * Copyright (c) 2023, CERN
8 *
9 * Redistribution and use in source and binary forms,
10 * with or without modification, are permitted according to the terms
12 */
13
14#ifndef RooFit_Detail_AnalyticalIntegrals_h
15#define RooFit_Detail_AnalyticalIntegrals_h
16
18
19#include <TMath.h>
21
22#include <cmath>
23
24namespace RooFit {
25
26namespace Detail {
27
28// For RooCBShape
29inline double approxErf(double arg)
30{
31 if (arg > 5.0)
32 return 1.0;
33 if (arg < -5.0)
34 return -1.0;
35
36 return TMath::Erf(arg);
37}
38
39namespace AnalyticalIntegrals {
40
41/// @brief Function to calculate the integral of an un-normalized RooGaussian over x. To calculate the integral over
42/// mean, just interchange the respective values of x and mean.
43/// @param xMin Minimum value of variable to integrate wrt.
44/// @param xMax Maximum value of of variable to integrate wrt.
45/// @param mean Mean.
46/// @param sigma Sigma.
47/// @return The integral of an un-normalized RooGaussian over the value in x.
48inline double gaussianIntegral(double xMin, double xMax, double mean, double sigma)
49{
50 // The normalisation constant 1./sqrt(2*pi*sigma^2) is left out in evaluate().
51 // Therefore, the integral is scaled up by that amount to make RooFit normalise
52 // correctly.
53 double resultScale = 0.5 * std::sqrt(TMath::TwoPi()) * sigma;
54
55 // Here everything is scaled and shifted into a standard normal distribution:
56 double xscale = TMath::Sqrt2() * sigma;
57 double scaledMin = 0.;
58 double scaledMax = 0.;
59 scaledMin = (xMin - mean) / xscale;
60 scaledMax = (xMax - mean) / xscale;
61
62 // Here we go for maximum precision: We compute all integrals in the UPPER
63 // tail of the Gaussian, because erfc has the highest precision there.
64 // Therefore, the different cases for range limits in the negative hemisphere are mapped onto
65 // the equivalent points in the upper hemisphere using erfc(-x) = 2. - erfc(x)
66 double ecmin = TMath::Erfc(std::abs(scaledMin));
67 double ecmax = TMath::Erfc(std::abs(scaledMax));
68
69 double cond = 0.0;
70 // Don't put this "prd" inside the "if" because clad will not be able to differentiate the code correctly (as of
71 // v1.1)!
72 double prd = scaledMin * scaledMax;
73 if (prd < 0.0) {
74 cond = 2.0 - (ecmin + ecmax);
75 } else if (scaledMax <= 0.0) {
76 cond = ecmax - ecmin;
77 } else {
78 cond = ecmin - ecmax;
79 }
80 return resultScale * cond;
81}
82
83inline double bifurGaussIntegral(double xMin, double xMax, double mean, double sigmaL, double sigmaR)
84{
85 const double xscaleL = TMath::Sqrt2() * sigmaL;
86 const double xscaleR = TMath::Sqrt2() * sigmaR;
87
88 const double resultScale = 0.5 * std::sqrt(TMath::TwoPi());
89
90 if (xMax < mean) {
91 return resultScale * (sigmaL * (TMath::Erf((xMax - mean) / xscaleL) - TMath::Erf((xMin - mean) / xscaleL)));
92 } else if (xMin > mean) {
93 return resultScale * (sigmaR * (TMath::Erf((xMax - mean) / xscaleR) - TMath::Erf((xMin - mean) / xscaleR)));
94 } else {
95 return resultScale *
96 (sigmaR * TMath::Erf((xMax - mean) / xscaleR) - sigmaL * TMath::Erf((xMin - mean) / xscaleL));
97 }
98}
99
100inline double exponentialIntegral(double xMin, double xMax, double constant)
101{
102 if (constant == 0.0) {
103 return xMax - xMin;
104 }
105
106 return (std::exp(constant * xMax) - std::exp(constant * xMin)) / constant;
107}
108
109/// In pdfMode, a coefficient for the constant term of 1.0 is implied if lowestOrder > 0.
110template <bool pdfMode = false>
111inline double polynomialIntegral(double const *coeffs, int nCoeffs, int lowestOrder, double xMin, double xMax)
112{
113 int denom = lowestOrder + nCoeffs;
114 double min = coeffs[nCoeffs - 1] / double(denom);
115 double max = coeffs[nCoeffs - 1] / double(denom);
116
117 for (int i = nCoeffs - 2; i >= 0; i--) {
118 denom--;
119 min = (coeffs[i] / double(denom)) + xMin * min;
120 max = (coeffs[i] / double(denom)) + xMax * max;
121 }
122
123 max = max * std::pow(xMax, 1 + lowestOrder);
124 min = min * std::pow(xMin, 1 + lowestOrder);
125
126 return max - min + (pdfMode && lowestOrder > 0.0 ? xMax - xMin : 0.0);
127}
128
129/// use fast FMA if available, fall back to normal arithmetic if not
130inline double fast_fma(double x, double y, double z) noexcept
131{
132#if defined(FP_FAST_FMA) // check if std::fma has fast hardware implementation
133 return std::fma(x, y, z);
134#else // defined(FP_FAST_FMA)
135 // std::fma might be slow, so use a more pedestrian implementation
136#if defined(__clang__)
137#pragma STDC FP_CONTRACT ON // hint clang that using an FMA is okay here
138#endif // defined(__clang__)
139 return (x * y) + z;
140#endif // defined(FP_FAST_FMA)
141}
142
143inline double chebychevIntegral(double const *coeffs, unsigned int nCoeffs, double xMin, double xMax, double xMinFull,
144 double xMaxFull)
145{
146 const double halfrange = .5 * (xMax - xMin);
147 const double mid = .5 * (xMax + xMin);
148
149 // the full range of the function is mapped to the normalised [-1, 1] range
150 const double b = (xMaxFull - mid) / halfrange;
151 const double a = (xMinFull - mid) / halfrange;
152
153 // coefficient for integral(T_0(x)) is 1 (implicit), integrate by hand
154 // T_0(x) and T_1(x), and use for n > 1: integral(T_n(x) dx) =
155 // (T_n+1(x) / (n + 1) - T_n-1(x) / (n - 1)) / 2
156 double sum = b - a; // integrate T_0(x) by hand
157
158 const unsigned int iend = nCoeffs;
159 if (iend > 0) {
160 {
161 // integrate T_1(x) by hand...
162 const double c = coeffs[0];
163 sum = fast_fma(0.5 * (b + a) * (b - a), c, sum);
164 }
165 if (1 < iend) {
166 double bcurr = b;
167 double btwox = 2 * b;
168 double blast = 1;
169
170 double acurr = a;
171 double atwox = 2 * a;
172 double alast = 1;
173
174 double newval = atwox * acurr - alast;
175 alast = acurr;
176 acurr = newval;
177
178 newval = btwox * bcurr - blast;
179 blast = bcurr;
180 bcurr = newval;
181 double nminus1 = 1.;
182 for (unsigned int i = 1; iend != i; ++i) {
183 // integrate using recursion relation
184 const double c = coeffs[i];
185 const double term2 = (blast - alast) / nminus1;
186
187 newval = atwox * acurr - alast;
188 alast = acurr;
189 acurr = newval;
190
191 newval = btwox * bcurr - blast;
192 blast = bcurr;
193 bcurr = newval;
194
195 ++nminus1;
196 const double term1 = (bcurr - acurr) / (nminus1 + 1.);
197 const double intTn = 0.5 * (term1 - term2);
198 sum = fast_fma(intTn, c, sum);
199 }
200 }
201 }
202
203 // take care to multiply with the right factor to account for the mapping to
204 // normalised range [-1, 1]
205 return halfrange * sum;
206}
207
208// Clad does not like std::max and std::min so redefined here for simplicity.
209inline double max(double x, double y)
210{
211 return x >= y ? x : y;
212}
213
214inline double min(double x, double y)
215{
216 return x <= y ? x : y;
217}
218
219// The last param should be of type bool but it is not as that causes some issues with Cling for some reason...
220inline double
221poissonIntegral(int code, double mu, double x, double integrandMin, double integrandMax, unsigned int protectNegative)
222{
223 if (protectNegative && mu < 0.0) {
224 return std::exp(-2.0 * mu); // make it fall quickly
225 }
226
227 if (code == 1) {
228 // Implement integral over x as summation. Add special handling in case
229 // range boundaries are not on integer values of x
230 integrandMin = max(0, integrandMin);
231
232 if (integrandMax < 0. || integrandMax < integrandMin) {
233 return 0;
234 }
235 const double delta = 100.0 * std::sqrt(mu);
236 // If the limits are more than many standard deviations away from the mean,
237 // we might as well return the integral of the full Poisson distribution to
238 // save computing time.
239 if (integrandMin < max(mu - delta, 0.0) && integrandMax > mu + delta) {
240 return 1.;
241 }
242
243 // The range as integers. ixMin is included, ixMax outside.
244 const unsigned int ixMin = integrandMin;
245 const unsigned int ixMax = min(integrandMax + 1, (double)std::numeric_limits<unsigned int>::max());
246
247 // Sum from 0 to just before the bin outside of the range.
248 if (ixMin == 0) {
249 return ROOT::Math::gamma_cdf_c(mu, ixMax, 1);
250 } else {
251 // If necessary, subtract from 0 to the beginning of the range
252 if (ixMin <= mu) {
253 return ROOT::Math::gamma_cdf_c(mu, ixMax, 1) - ROOT::Math::gamma_cdf_c(mu, ixMin, 1);
254 } else {
255 // Avoid catastrophic cancellation in the high tails:
256 return ROOT::Math::gamma_cdf(mu, ixMin, 1) - ROOT::Math::gamma_cdf(mu, ixMax, 1);
257 }
258 }
259 }
260
261 // the integral with respect to the mean is the integral of a gamma distribution
262 // negative ix does not need protection (gamma returns 0.0)
263 const double ix = 1 + x;
264
265 return ROOT::Math::gamma_cdf(integrandMax, ix, 1.0) - ROOT::Math::gamma_cdf(integrandMin, ix, 1.0);
266}
267
268inline double logNormalIntegral(double xMin, double xMax, double m0, double k)
269{
270 const double root2 = std::sqrt(2.);
271
272 double ln_k = std::abs(std::log(k));
273 double ret =
274 0.5 * (TMath::Erf(std::log(xMax / m0) / (root2 * ln_k)) - TMath::Erf(std::log(xMin / m0) / (root2 * ln_k)));
275
276 return ret;
277}
278
279inline double logNormalIntegralStandard(double xMin, double xMax, double mu, double sigma)
280{
281 const double root2 = std::sqrt(2.);
282
283 double ln_k = std::abs(sigma);
284 double ret =
285 0.5 * (TMath::Erf((std::log(xMax) - mu) / (root2 * ln_k)) - TMath::Erf((std::log(xMin) - mu) / (root2 * ln_k)));
286
287 return ret;
288}
289
290inline double cbShapeIntegral(double mMin, double mMax, double m0, double sigma, double alpha, double n)
291{
292 const double sqrtPiOver2 = 1.2533141373;
293 const double sqrt2 = 1.4142135624;
294
295 double result = 0.0;
296 bool useLog = false;
297
298 if (std::abs(n - 1.0) < 1.0e-05)
299 useLog = true;
300
301 double sig = std::abs((double)sigma);
302
303 double tmin = (mMin - m0) / sig;
304 double tmax = (mMax - m0) / sig;
305
306 if (alpha < 0) {
307 double tmp = tmin;
308 tmin = -tmax;
309 tmax = -tmp;
310 }
311
312 double absAlpha = std::abs((double)alpha);
313
314 if (tmin >= -absAlpha) {
315 result += sig * sqrtPiOver2 * (RooFit::Detail::approxErf(tmax / sqrt2) - RooFit::Detail::approxErf(tmin / sqrt2));
316 } else if (tmax <= -absAlpha) {
317 double a = std::pow(n / absAlpha, n) * std::exp(-0.5 * absAlpha * absAlpha);
318 double b = n / absAlpha - absAlpha;
319
320 if (useLog) {
321 result += a * sig * (std::log(b - tmin) - std::log(b - tmax));
322 } else {
323 result += a * sig / (1.0 - n) * (1.0 / (std::pow(b - tmin, n - 1.0)) - 1.0 / (std::pow(b - tmax, n - 1.0)));
324 }
325 } else {
326 double a = std::pow(n / absAlpha, n) * std::exp(-0.5 * absAlpha * absAlpha);
327 double b = n / absAlpha - absAlpha;
328
329 double term1 = 0.0;
330 if (useLog) {
331 term1 = a * sig * (std::log(b - tmin) - std::log(n / absAlpha));
332 } else {
333 term1 = a * sig / (1.0 - n) * (1.0 / (std::pow(b - tmin, n - 1.0)) - 1.0 / (std::pow(n / absAlpha, n - 1.0)));
334 }
335
336 double term2 =
337 sig * sqrtPiOver2 * (RooFit::Detail::approxErf(tmax / sqrt2) - RooFit::Detail::approxErf(-absAlpha / sqrt2));
338
339 result += term1 + term2;
340 }
341
342 if (result == 0)
343 return 1.E-300;
344 return result;
345}
346
347inline double bernsteinIntegral(double xlo, double xhi, double xmin, double xmax, double *coefs, int nCoefs)
348{
349 double xloScaled = (xlo - xmin) / (xmax - xmin);
350 double xhiScaled = (xhi - xmin) / (xmax - xmin);
351
352 int degree = nCoefs - 1; // n+1 polys of degree n
353 double norm = 0.;
354
355 for (int i = 0; i <= degree; ++i) {
356 // for each of the i Bernstein basis polynomials
357 // represent it in the 'power basis' (the naive polynomial basis)
358 // where the integral is straight forward.
359 double temp = 0.;
360 for (int j = i; j <= degree; ++j) { // power basis≈ß
361 double binCoefs = binomial(degree, j) * binomial(j, i);
362 double oneOverJPlusOne = 1. / (j + 1.);
363 double powDiff = std::pow(xhiScaled, j + 1.) - std::pow(xloScaled, j + 1.);
364 temp += std::pow(-1., j - i) * binCoefs * powDiff * oneOverJPlusOne;
365 }
366 temp *= coefs[i]; // include coeff
367 norm += temp; // add this basis's contribution to total
368 }
369
370 return norm * (xmax - xmin);
371}
372
373} // namespace AnalyticalIntegrals
374
375} // namespace Detail
376
377} // namespace RooFit
378
379#endif
#define b(i)
Definition RSha256.hxx:100
#define c(i)
Definition RSha256.hxx:101
#define a(i)
Definition RSha256.hxx:99
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
float xmin
float xmax
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).
const Double_t sigma
Double_t y[n]
Definition legend1.C:17
Double_t x[n]
Definition legend1.C:17
const Int_t n
Definition legend1.C:16
double cbShapeIntegral(double mMin, double mMax, double m0, double sigma, double alpha, double n)
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 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 bernsteinIntegral(double xlo, double xhi, double xmin, double xmax, double *coefs, int nCoefs)
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.
double binomial(int n, int k)
Calculates the binomial coefficient n over k.
double approxErf(double arg)
The namespace RooFit contains mostly switches that change the behaviour of functions of PDFs (or othe...
Definition JSONIO.h:26
Double_t Erf(Double_t x)
Computation of the error function erf(x).
Definition TMath.cxx:190
constexpr Double_t Sqrt2()
Definition TMath.h:86
Double_t Erfc(Double_t x)
Computes the complementary error function erfc(x).
Definition TMath.cxx:199
constexpr Double_t TwoPi()
Definition TMath.h:44
static uint64_t sum(uint64_t i)
Definition Factory.cxx:2345