ROOT   Reference Guide
Searching...
No Matches
MathFuncs.h
Go to the documentation of this file.
1/*
2 * Project: RooFit
3 * Authors:
4 * Jonas Rembser, CERN 2024
5 * Garima Singh, CERN 2023
6 *
7 * Copyright (c) 2024, 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_MathFuncs_h
15#define RooFit_Detail_MathFuncs_h
16
17#include <TMath.h>
20
21#include <cmath>
22
23namespace RooFit {
24
25namespace Detail {
26
27namespace MathFuncs {
28
29/// Calculates the binomial coefficient n over k.
30/// Equivalent to TMath::Binomial, but inlined.
31inline double binomial(int n, int k)
32{
33 if (n < 0 || k < 0 || n < k)
34 return TMath::SignalingNaN();
35 if (k == 0 || n == k)
36 return 1;
37
38 int k1 = std::min(k, n - k);
39 int k2 = n - k1;
40 double fact = k2 + 1;
41 for (double i = k1; i > 1.; --i) {
42 fact *= (k2 + i) / i;
43 }
44 return fact;
45}
46
47/// The caller needs to make sure that there is at least one coefficient.
48inline double bernstein(double x, double xmin, double xmax, double *coefs, int nCoefs)
49{
50 double xScaled = (x - xmin) / (xmax - xmin); // rescale to [0,1]
51 int degree = nCoefs - 1; // n+1 polys of degree n
52
53 // in case list of arguments passed is empty
54 if (degree < 0) {
55 return TMath::SignalingNaN();
56 } else if (degree == 0) {
57 return coefs[0];
58 } else if (degree == 1) {
59
60 double a0 = coefs[0]; // c0
61 double a1 = coefs[1] - a0; // c1 - c0
62 return a1 * xScaled + a0;
63
64 } else if (degree == 2) {
65
66 double a0 = coefs[0]; // c0
67 double a1 = 2 * (coefs[1] - a0); // 2 * (c1 - c0)
68 double a2 = coefs[2] - a1 - a0; // c0 - 2 * c1 + c2
69 return (a2 * xScaled + a1) * xScaled + a0;
70 }
71
72 double t = xScaled;
73 double s = 1. - xScaled;
74
75 double result = coefs[0] * s;
76 for (int i = 1; i < degree; i++) {
77 result = (result + t * binomial(degree, i) * coefs[i]) * s;
78 t *= xScaled;
79 }
80 result += t * coefs[degree];
81
82 return result;
83}
84
85/// @brief Function to evaluate an un-normalized RooGaussian.
86inline double gaussian(double x, double mean, double sigma)
87{
88 const double arg = x - mean;
89 const double sig = sigma;
90 return std::exp(-0.5 * arg * arg / (sig * sig));
91}
92
93// RooRatio evaluate function.
94inline double ratio(double numerator, double denominator)
95{
96 return numerator / denominator;
97}
98
99inline double bifurGauss(double x, double mean, double sigmaL, double sigmaR)
100{
101 // Note: this simplification does not work with Clad as of v1.1!
102 // return gaussian(x, mean, x < mean ? sigmaL : sigmaR);
103 if (x < mean)
104 return gaussian(x, mean, sigmaL);
105 return gaussian(x, mean, sigmaR);
106}
107
108inline double efficiency(double effFuncVal, int catIndex, int sigCatIndex)
109{
110 // Truncate efficiency function in range 0.0-1.0
111 effFuncVal = std::clamp(effFuncVal, 0.0, 1.0);
112
113 if (catIndex == sigCatIndex)
114 return effFuncVal; // Accept case
115 else
116 return 1 - effFuncVal; // Reject case
117}
118
119/// In pdfMode, a coefficient for the constant term of 1.0 is implied if lowestOrder > 0.
120template <bool pdfMode = false>
121inline double polynomial(double const *coeffs, int nCoeffs, int lowestOrder, double x)
122{
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);
128}
129
130inline double chebychev(double *coeffs, unsigned int nCoeffs, double x_in, double xMin, double xMax)
131{
132 // transform to range [-1, +1]
133 const double xPrime = (x_in - 0.5 * (xMax + xMin)) / (0.5 * (xMax - xMin));
134
135 // extract current values of coefficients
136 double sum = 1.;
137 if (nCoeffs > 0) {
138 double curr = xPrime;
139 double twox = 2 * xPrime;
140 double last = 1;
141 double newval = twox * curr - last;
142 last = curr;
143 curr = newval;
144 for (unsigned int i = 0; nCoeffs != i; ++i) {
145 sum += last * coeffs[i];
146 newval = twox * curr - last;
147 last = curr;
148 curr = newval;
149 }
150 }
151 return sum;
152}
153
154inline double constraintSum(double const *comp, unsigned int compSize)
155{
156 double sum = 0;
157 for (unsigned int i = 0; i < compSize; i++) {
158 sum -= std::log(comp[i]);
159 }
160 return sum;
161}
162
163inline unsigned int getUniformBinning(double low, double high, double val, unsigned int numBins)
164{
165 double binWidth = (high - low) / numBins;
166 return val >= high ? numBins - 1 : std::abs((val - low) / binWidth);
167}
168
169inline double poisson(double x, double par)
170{
171 if (par < 0)
172 return TMath::QuietNaN();
173
174 if (x < 0) {
175 return 0;
176 } else if (x == 0.0) {
177 return std::exp(-par);
178 } else {
179 double out = x * std::log(par) - TMath::LnGamma(x + 1.) - par;
180 return std::exp(out);
181 }
182}
183
184inline double flexibleInterpSingle(unsigned int code, double low, double high, double boundary, double nominal,
185 double paramVal, double res)
186{
187 if (code == 0) {
188 // piece-wise linear
189 if (paramVal > 0) {
190 return paramVal * (high - nominal);
191 } else {
192 return paramVal * (nominal - low);
193 }
194 } else if (code == 1) {
195 // piece-wise log
196 if (paramVal >= 0) {
197 return res * (std::pow(high / nominal, +paramVal) - 1);
198 } else {
199 return res * (std::pow(low / nominal, -paramVal) - 1);
200 }
201 } else if (code == 2) {
202 // parabolic with linear
203 double a = 0.5 * (high + low) - nominal;
204 double b = 0.5 * (high - low);
205 double c = 0;
206 if (paramVal > 1) {
207 return (2 * a + b) * (paramVal - 1) + high - nominal;
208 } else if (paramVal < -1) {
209 return -1 * (2 * a - b) * (paramVal + 1) + low - nominal;
210 } else {
211 return a * std::pow(paramVal, 2) + b * paramVal + c;
212 }
213 } else if (code == 3) {
214 // parabolic version of log-normal
215 double a = 0.5 * (high + low) - nominal;
216 double b = 0.5 * (high - low);
217 double c = 0;
218 if (paramVal > 1) {
219 return (2 * a + b) * (paramVal - 1) + high - nominal;
220 } else if (paramVal < -1) {
221 return -1 * (2 * a - b) * (paramVal + 1) + low - nominal;
222 } else {
223 return a * std::pow(paramVal, 2) + b * paramVal + c;
224 }
225 } else if (code == 4) {
226 double x = paramVal;
227 if (x >= boundary) {
228 return x * (high - nominal);
229 } else if (x <= -boundary) {
230 return x * (nominal - low);
231 }
232
233 // interpolate 6th degree
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);
239
240 return x * (S + t * A * (15 + t * t * (-10 + t * t * 3)));
241 } else if (code == 5) {
242 double x = paramVal;
243 double mod = 1.0;
244 if (x >= boundary) {
245 mod = std::pow(high / nominal, +paramVal);
246 } else if (x <= -boundary) {
247 mod = std::pow(low / nominal, -paramVal);
248 } else {
249 // interpolate 6th degree exp
250 double x0 = boundary;
251
252 // GHL: Swagato's suggestions
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;
261
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);
268
269 // fcns+der+2nd_der are eq at bd
270
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);
277
278 // evaluate the 6-th degree polynomial using Horner's method
279 double value = 1. + x * (a + x * (b + x * (c + x * (d + x * (e + x * f)))));
280 mod = value;
281 }
282 return res * (mod - 1.0);
283 }
284
285 return 0.0;
286}
287
288inline double flexibleInterp(unsigned int code, double *params, unsigned int n, double *low, double *high,
289 double boundary, double nominal, int doCutoff)
290{
291 double total = nominal;
292 for (std::size_t i = 0; i < n; ++i) {
293 total += flexibleInterpSingle(code, low[i], high[i], boundary, nominal, params[i], total);
294 }
295
296 return doCutoff && total <= 0 ? TMath::Limits<double>::Min() : total;
297}
298
299inline double landau(double x, double mu, double sigma)
300{
301 if (sigma <= 0.)
302 return 0.;
303 return ROOT::Math::landau_pdf((x - mu) / sigma);
304}
305
306inline double logNormal(double x, double k, double m0)
307{
308 return ROOT::Math::lognormal_pdf(x, std::log(m0), std::abs(std::log(k)));
309}
310
311inline double logNormalStandard(double x, double sigma, double mu)
312{
313 return ROOT::Math::lognormal_pdf(x, mu, std::abs(sigma));
314}
315
316inline double effProd(double eff, double pdf)
317{
318 return eff * pdf;
319}
320
321inline double nll(double pdf, double weight, int binnedL, int doBinOffset)
322{
323 if (binnedL) {
324 // Special handling of this case since std::log(Poisson(0,0)=0 but can't be
325 // calculated with usual log-formula since std::log(mu)=0. No update of result
326 // is required since term=0.
327 if (std::abs(pdf) < 1e-10 && std::abs(weight) < 1e-10) {
328 return 0.0;
329 }
330 if (doBinOffset) {
331 return pdf - weight - weight * (std::log(pdf) - std::log(weight));
332 }
333 return pdf - weight * std::log(pdf) + TMath::LnGamma(weight + 1);
334 } else {
335 return -weight * std::log(pdf);
336 }
337}
338
339inline double recursiveFraction(double *a, unsigned int n)
340{
341 double prod = a[0];
342
343 for (unsigned int i = 1; i < n; ++i) {
344 prod *= 1.0 - a[i];
345 }
346
347 return prod;
348}
349
350inline double cbShape(double m, double m0, double sigma, double alpha, double n)
351{
352 double t = (m - m0) / sigma;
353 if (alpha < 0)
354 t = -t;
355
356 double absAlpha = std::abs((double)alpha);
357
358 if (t >= -absAlpha) {
359 return std::exp(-0.5 * t * t);
360 } else {
361 double a = std::pow(n / absAlpha, n) * std::exp(-0.5 * absAlpha * absAlpha);
362 double b = n / absAlpha - absAlpha;
363
364 return a / std::pow(b - t, n);
365 }
366}
367
368// For RooCBShape
369inline double approxErf(double arg)
370{
371 if (arg > 5.0)
372 return 1.0;
373 if (arg < -5.0)
374 return -1.0;
375
376 return TMath::Erf(arg);
377}
378
379/// @brief Function to calculate the integral of an un-normalized RooGaussian over x. To calculate the integral over
380/// mean, just interchange the respective values of x and mean.
381/// @param xMin Minimum value of variable to integrate wrt.
382/// @param xMax Maximum value of of variable to integrate wrt.
383/// @param mean Mean.
384/// @param sigma Sigma.
385/// @return The integral of an un-normalized RooGaussian over the value in x.
386inline double gaussianIntegral(double xMin, double xMax, double mean, double sigma)
387{
388 // The normalisation constant 1./sqrt(2*pi*sigma^2) is left out in evaluate().
389 // Therefore, the integral is scaled up by that amount to make RooFit normalise
390 // correctly.
391 double resultScale = 0.5 * std::sqrt(TMath::TwoPi()) * sigma;
392
393 // Here everything is scaled and shifted into a standard normal distribution:
394 double xscale = TMath::Sqrt2() * sigma;
395 double scaledMin = 0.;
396 double scaledMax = 0.;
397 scaledMin = (xMin - mean) / xscale;
398 scaledMax = (xMax - mean) / xscale;
399
400 // Here we go for maximum precision: We compute all integrals in the UPPER
401 // tail of the Gaussian, because erfc has the highest precision there.
402 // Therefore, the different cases for range limits in the negative hemisphere are mapped onto
403 // the equivalent points in the upper hemisphere using erfc(-x) = 2. - erfc(x)
404 double ecmin = TMath::Erfc(std::abs(scaledMin));
405 double ecmax = TMath::Erfc(std::abs(scaledMax));
406
407 double cond = 0.0;
408 // Don't put this "prd" inside the "if" because clad will not be able to differentiate the code correctly (as of
409 // v1.1)!
410 double prd = scaledMin * scaledMax;
411 if (prd < 0.0) {
412 cond = 2.0 - (ecmin + ecmax);
413 } else if (scaledMax <= 0.0) {
414 cond = ecmax - ecmin;
415 } else {
416 cond = ecmin - ecmax;
417 }
418 return resultScale * cond;
419}
420
421inline double bifurGaussIntegral(double xMin, double xMax, double mean, double sigmaL, double sigmaR)
422{
423 const double xscaleL = TMath::Sqrt2() * sigmaL;
424 const double xscaleR = TMath::Sqrt2() * sigmaR;
425
426 const double resultScale = 0.5 * std::sqrt(TMath::TwoPi());
427
428 if (xMax < mean) {
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)));
432 } else {
433 return resultScale *
434 (sigmaR * TMath::Erf((xMax - mean) / xscaleR) - sigmaL * TMath::Erf((xMin - mean) / xscaleL));
435 }
436}
437
438inline double exponentialIntegral(double xMin, double xMax, double constant)
439{
440 if (constant == 0.0) {
441 return xMax - xMin;
442 }
443
444 return (std::exp(constant * xMax) - std::exp(constant * xMin)) / constant;
445}
446
447/// In pdfMode, a coefficient for the constant term of 1.0 is implied if lowestOrder > 0.
448template <bool pdfMode = false>
449inline double polynomialIntegral(double const *coeffs, int nCoeffs, int lowestOrder, double xMin, double xMax)
450{
451 int denom = lowestOrder + nCoeffs;
452 double min = coeffs[nCoeffs - 1] / double(denom);
453 double max = coeffs[nCoeffs - 1] / double(denom);
454
455 for (int i = nCoeffs - 2; i >= 0; i--) {
456 denom--;
457 min = (coeffs[i] / double(denom)) + xMin * min;
458 max = (coeffs[i] / double(denom)) + xMax * max;
459 }
460
461 max = max * std::pow(xMax, 1 + lowestOrder);
462 min = min * std::pow(xMin, 1 + lowestOrder);
463
464 return max - min + (pdfMode && lowestOrder > 0.0 ? xMax - xMin : 0.0);
465}
466
467/// use fast FMA if available, fall back to normal arithmetic if not
468inline double fast_fma(double x, double y, double z) noexcept
469{
470#if defined(FP_FAST_FMA) // check if std::fma has fast hardware implementation
471 return std::fma(x, y, z);
472#else // defined(FP_FAST_FMA)
473 // std::fma might be slow, so use a more pedestrian implementation
474#if defined(__clang__)
475#pragma STDC FP_CONTRACT ON // hint clang that using an FMA is okay here
476#endif // defined(__clang__)
477 return (x * y) + z;
478#endif // defined(FP_FAST_FMA)
479}
480
481inline double chebychevIntegral(double const *coeffs, unsigned int nCoeffs, double xMin, double xMax, double xMinFull,
482 double xMaxFull)
483{
484 const double halfrange = .5 * (xMax - xMin);
485 const double mid = .5 * (xMax + xMin);
486
487 // the full range of the function is mapped to the normalised [-1, 1] range
488 const double b = (xMaxFull - mid) / halfrange;
489 const double a = (xMinFull - mid) / halfrange;
490
491 // coefficient for integral(T_0(x)) is 1 (implicit), integrate by hand
492 // T_0(x) and T_1(x), and use for n > 1: integral(T_n(x) dx) =
493 // (T_n+1(x) / (n + 1) - T_n-1(x) / (n - 1)) / 2
494 double sum = b - a; // integrate T_0(x) by hand
495
496 const unsigned int iend = nCoeffs;
497 if (iend > 0) {
498 {
499 // integrate T_1(x) by hand...
500 const double c = coeffs[0];
501 sum = fast_fma(0.5 * (b + a) * (b - a), c, sum);
502 }
503 if (1 < iend) {
504 double bcurr = b;
505 double btwox = 2 * b;
506 double blast = 1;
507
508 double acurr = a;
509 double atwox = 2 * a;
510 double alast = 1;
511
512 double newval = atwox * acurr - alast;
513 alast = acurr;
514 acurr = newval;
515
516 newval = btwox * bcurr - blast;
517 blast = bcurr;
518 bcurr = newval;
519 double nminus1 = 1.;
520 for (unsigned int i = 1; iend != i; ++i) {
521 // integrate using recursion relation
522 const double c = coeffs[i];
523 const double term2 = (blast - alast) / nminus1;
524
525 newval = atwox * acurr - alast;
526 alast = acurr;
527 acurr = newval;
528
529 newval = btwox * bcurr - blast;
530 blast = bcurr;
531 bcurr = newval;
532
533 ++nminus1;
534 const double term1 = (bcurr - acurr) / (nminus1 + 1.);
535 const double intTn = 0.5 * (term1 - term2);
536 sum = fast_fma(intTn, c, sum);
537 }
538 }
539 }
540
541 // take care to multiply with the right factor to account for the mapping to
542 // normalised range [-1, 1]
543 return halfrange * sum;
544}
545
546// The last param should be of type bool but it is not as that causes some issues with Cling for some reason...
547inline double
548poissonIntegral(int code, double mu, double x, double integrandMin, double integrandMax, unsigned int protectNegative)
549{
550 if (protectNegative && mu < 0.0) {
551 return std::exp(-2.0 * mu); // make it fall quickly
552 }
553
554 if (code == 1) {
555 // Implement integral over x as summation. Add special handling in case
556 // range boundaries are not on integer values of x
557 integrandMin = std::max(0., integrandMin);
558
559 if (integrandMax < 0. || integrandMax < integrandMin) {
560 return 0;
561 }
562 const double delta = 100.0 * std::sqrt(mu);
563 // If the limits are more than many standard deviations away from the mean,
564 // we might as well return the integral of the full Poisson distribution to
565 // save computing time.
566 if (integrandMin < std::max(mu - delta, 0.0) && integrandMax > mu + delta) {
567 return 1.;
568 }
569
570 // The range as integers. ixMin is included, ixMax outside.
571 const unsigned int ixMin = integrandMin;
572 const unsigned int ixMax = std::min(integrandMax + 1, (double)std::numeric_limits<unsigned int>::max());
573
574 // Sum from 0 to just before the bin outside of the range.
575 if (ixMin == 0) {
576 return ROOT::Math::gamma_cdf_c(mu, ixMax, 1);
577 } else {
578 // If necessary, subtract from 0 to the beginning of the range
579 if (ixMin <= mu) {
580 return ROOT::Math::gamma_cdf_c(mu, ixMax, 1) - ROOT::Math::gamma_cdf_c(mu, ixMin, 1);
581 } else {
582 // Avoid catastrophic cancellation in the high tails:
583 return ROOT::Math::gamma_cdf(mu, ixMin, 1) - ROOT::Math::gamma_cdf(mu, ixMax, 1);
584 }
585 }
586 }
587
588 // the integral with respect to the mean is the integral of a gamma distribution
589 // negative ix does not need protection (gamma returns 0.0)
590 const double ix = 1 + x;
591
592 return ROOT::Math::gamma_cdf(integrandMax, ix, 1.0) - ROOT::Math::gamma_cdf(integrandMin, ix, 1.0);
593}
594
595inline double logNormalIntegral(double xMin, double xMax, double m0, double k)
596{
597 const double root2 = std::sqrt(2.);
598
599 double ln_k = std::abs(std::log(k));
600 double ret =
601 0.5 * (TMath::Erf(std::log(xMax / m0) / (root2 * ln_k)) - TMath::Erf(std::log(xMin / m0) / (root2 * ln_k)));
602
603 return ret;
604}
605
606inline double logNormalIntegralStandard(double xMin, double xMax, double mu, double sigma)
607{
608 const double root2 = std::sqrt(2.);
609
610 double ln_k = std::abs(sigma);
611 double ret =
612 0.5 * (TMath::Erf((std::log(xMax) - mu) / (root2 * ln_k)) - TMath::Erf((std::log(xMin) - mu) / (root2 * ln_k)));
613
614 return ret;
615}
616
617inline double cbShapeIntegral(double mMin, double mMax, double m0, double sigma, double alpha, double n)
618{
619 const double sqrtPiOver2 = 1.2533141373;
620 const double sqrt2 = 1.4142135624;
621
622 double result = 0.0;
623 bool useLog = false;
624
625 if (std::abs(n - 1.0) < 1.0e-05)
626 useLog = true;
627
628 double sig = std::abs(sigma);
629
630 double tmin = (mMin - m0) / sig;
631 double tmax = (mMax - m0) / sig;
632
633 if (alpha < 0) {
634 double tmp = tmin;
635 tmin = -tmax;
636 tmax = -tmp;
637 }
638
639 double absAlpha = std::abs(alpha);
640
641 if (tmin >= -absAlpha) {
642 result += sig * sqrtPiOver2 * (approxErf(tmax / sqrt2) - approxErf(tmin / sqrt2));
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;
646
647 if (useLog) {
648 result += a * sig * (std::log(b - tmin) - std::log(b - tmax));
649 } else {
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)));
651 }
652 } else {
653 double a = std::pow(n / absAlpha, n) * std::exp(-0.5 * absAlpha * absAlpha);
654 double b = n / absAlpha - absAlpha;
655
656 double term1 = 0.0;
657 if (useLog) {
658 term1 = a * sig * (std::log(b - tmin) - std::log(n / absAlpha));
659 } else {
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)));
661 }
662
663 double term2 = sig * sqrtPiOver2 * (approxErf(tmax / sqrt2) - approxErf(-absAlpha / sqrt2));
664
665 result += term1 + term2;
666 }
667
668 if (result == 0)
669 return 1.E-300;
670 return result;
671}
672
673inline double bernsteinIntegral(double xlo, double xhi, double xmin, double xmax, double *coefs, int nCoefs)
674{
675 double xloScaled = (xlo - xmin) / (xmax - xmin);
676 double xhiScaled = (xhi - xmin) / (xmax - xmin);
677
678 int degree = nCoefs - 1; // n+1 polys of degree n
679 double norm = 0.;
680
681 for (int i = 0; i <= degree; ++i) {
682 // for each of the i Bernstein basis polynomials
683 // represent it in the 'power basis' (the naive polynomial basis)
684 // where the integral is straight forward.
685 double temp = 0.;
686 for (int j = i; j <= degree; ++j) { // power basis≈ß
687 double binCoefs = binomial(degree, j) * binomial(j, i);
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;
691 }
692 temp *= coefs[i]; // include coeff
693 norm += temp; // add this basis's contribution to total
694 }
695
696 return norm * (xmax - xmin);
697}
698
699} // namespace MathFuncs
700
701} // namespace Detail
702
703} // namespace RooFit
704
705#endif
#define d(i)
Definition RSha256.hxx:102
#define b(i)
Definition RSha256.hxx:100
#define f(i)
Definition RSha256.hxx:104
#define S0(x)
Definition RSha256.hxx:88
#define S1(x)
Definition RSha256.hxx:89
#define c(i)
Definition RSha256.hxx:101
#define a(i)
Definition RSha256.hxx:99
#define e(i)
Definition RSha256.hxx:103
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
float xmin
float xmax
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).
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 logNormalIntegral(double xMin, double xMax, double m0, double k)
Definition MathFuncs.h:595
double gaussianIntegral(double xMin, double xMax, double mean, double sigma)
Function to calculate the integral of an un-normalized RooGaussian over x.
Definition MathFuncs.h:386
double chebychevIntegral(double const *coeffs, unsigned int nCoeffs, double xMin, double xMax, double xMinFull, double xMaxFull)
Definition MathFuncs.h:481
double bifurGaussIntegral(double xMin, double xMax, double mean, double sigmaL, double sigmaR)
Definition MathFuncs.h:421
double cbShape(double m, double m0, double sigma, double alpha, double n)
Definition MathFuncs.h:350
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.
Definition MathFuncs.h:121
double recursiveFraction(double *a, unsigned int n)
Definition MathFuncs.h:339
double constraintSum(double const *comp, unsigned int compSize)
Definition MathFuncs.h:154
double cbShapeIntegral(double mMin, double mMax, double m0, double sigma, double alpha, double n)
Definition MathFuncs.h:617
double flexibleInterp(unsigned int code, double *params, unsigned int n, double *low, double *high, double boundary, double nominal, int doCutoff)
Definition MathFuncs.h:288
double fast_fma(double x, double y, double z) noexcept
use fast FMA if available, fall back to normal arithmetic if not
Definition MathFuncs.h:468
double logNormalIntegralStandard(double xMin, double xMax, double mu, double sigma)
Definition MathFuncs.h:606
double landau(double x, double mu, double sigma)
Definition MathFuncs.h:299
double gaussian(double x, double mean, double sigma)
Function to evaluate an un-normalized RooGaussian.
Definition MathFuncs.h:86
double chebychev(double *coeffs, unsigned int nCoeffs, double x_in, double xMin, double xMax)
Definition MathFuncs.h:130
double poisson(double x, double par)
Definition MathFuncs.h:169
double binomial(int n, int k)
Calculates the binomial coefficient n over k.
Definition MathFuncs.h:31
unsigned int getUniformBinning(double low, double high, double val, unsigned int numBins)
Definition MathFuncs.h:163
double flexibleInterpSingle(unsigned int code, double low, double high, double boundary, double nominal, double paramVal, double res)
Definition MathFuncs.h:184
double logNormalStandard(double x, double sigma, double mu)
Definition MathFuncs.h:311
double bifurGauss(double x, double mean, double sigmaL, double sigmaR)
Definition MathFuncs.h:99
double ratio(double numerator, double denominator)
Definition MathFuncs.h:94
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.
Definition MathFuncs.h:449
double bernsteinIntegral(double xlo, double xhi, double xmin, double xmax, double *coefs, int nCoefs)
Definition MathFuncs.h:673
double approxErf(double arg)
Definition MathFuncs.h:369
double effProd(double eff, double pdf)
Definition MathFuncs.h:316
double poissonIntegral(int code, double mu, double x, double integrandMin, double integrandMax, unsigned int protectNegative)
Definition MathFuncs.h:548
double logNormal(double x, double k, double m0)
Definition MathFuncs.h:306
double nll(double pdf, double weight, int binnedL, int doBinOffset)
Definition MathFuncs.h:321
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.
Definition MathFuncs.h:48
double efficiency(double effFuncVal, int catIndex, int sigCatIndex)
Definition MathFuncs.h:108
double exponentialIntegral(double xMin, double xMax, double constant)
Definition MathFuncs.h:438
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
Double_t QuietNaN()
Returns a quiet NaN as defined by IEEE 754.
Definition TMath.h:902
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
Double_t LnGamma(Double_t z)
Computation of ln[gamma(z)] for all z.
Definition TMath.cxx:509
Double_t SignalingNaN()
Returns a signaling NaN as defined by IEEE 754](http://en.wikipedia.org/wiki/NaN#Signaling_NaN).
Definition TMath.h:910
constexpr Double_t TwoPi()
Definition TMath.h:44
TMarker m
Definition textangle.C:8
static uint64_t sum(uint64_t i)
Definition Factory.cxx:2345