ROOT   Reference Guide
Searching...
No Matches
EvaluateFuncs.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_EvaluateFuncs_h
15#define RooFit_Detail_EvaluateFuncs_h
16
17#include <TMath.h>
19
20#include <cmath>
21
22namespace RooFit {
23
24namespace Detail {
25
26/// Calculates the binomial coefficient n over k.
27/// Equivalent to TMath::Binomial, but inlined.
28inline double binomial(int n, int k)
29{
30 if (n < 0 || k < 0 || n < k)
31 return TMath::SignalingNaN();
32 if (k == 0 || n == k)
33 return 1;
34
35 int k1 = std::min(k, n - k);
36 int k2 = n - k1;
37 double fact = k2 + 1;
38 for (double i = k1; i > 1.; --i) {
39 fact *= (k2 + i) / i;
40 }
41 return fact;
42}
43
44namespace EvaluateFuncs {
45
46/// The caller needs to make sure that there is at least one coefficient.
47inline double bernsteinEvaluate(double x, double xmin, double xmax, double *coefs, int nCoefs)
48{
49 double xScaled = (x - xmin) / (xmax - xmin); // rescale to [0,1]
50 int degree = nCoefs - 1; // n+1 polys of degree n
51
52 // in case list of arguments passed is empty
53 if(degree < 0) {
54 return TMath::SignalingNaN();
55 } else if (degree == 0) {
56 return coefs[0];
57 } else if (degree == 1) {
58
59 double a0 = coefs[0]; // c0
60 double a1 = coefs[1] - a0; // c1 - c0
61 return a1 * xScaled + a0;
62
63 } else if (degree == 2) {
64
65 double a0 = coefs[0]; // c0
66 double a1 = 2 * (coefs[1] - a0); // 2 * (c1 - c0)
67 double a2 = coefs[2] - a1 - a0; // c0 - 2 * c1 + c2
68 return (a2 * xScaled + a1) * xScaled + a0;
69 }
70
71 double t = xScaled;
72 double s = 1. - xScaled;
73
74 double result = coefs[0] * s;
75 for (int i = 1; i < degree; i++) {
76 result = (result + t * binomial(degree, i) * coefs[i]) * s;
77 t *= xScaled;
78 }
79 result += t * coefs[degree];
80
81 return result;
82}
83
84/// @brief Function to evaluate an un-normalized RooGaussian.
85inline double gaussianEvaluate(double x, double mean, double sigma)
86{
87 const double arg = x - mean;
88 const double sig = sigma;
89 return std::exp(-0.5 * arg * arg / (sig * sig));
90}
91
92// RooRatio evaluate function.
93inline double ratioEvaluate(double numerator, double denominator)
94{
95 return numerator / denominator;
96}
97
98inline double bifurGaussEvaluate(double x, double mean, double sigmaL, double sigmaR)
99{
100 // Note: this simplification does not work with Clad as of v1.1!
101 // return gaussianEvaluate(x, mean, x < mean ? sigmaL : sigmaR);
102 if (x < mean)
103 return gaussianEvaluate(x, mean, sigmaL);
104 return gaussianEvaluate(x, mean, sigmaR);
105}
106
107inline double efficiencyEvaluate(double effFuncVal, int catIndex, int sigCatIndex)
108{
109 // Truncate efficiency function in range 0.0-1.0
110 effFuncVal = std::clamp(effFuncVal, 0.0, 1.0);
111
112 if (catIndex == sigCatIndex)
113 return effFuncVal; // Accept case
114 else
115 return 1 - effFuncVal; // Reject case
116}
117
118/// In pdfMode, a coefficient for the constant term of 1.0 is implied if lowestOrder > 0.
119template <bool pdfMode = false>
120inline double polynomialEvaluate(double const *coeffs, int nCoeffs, int lowestOrder, double x)
121{
122 double retVal = coeffs[nCoeffs - 1];
123 for (int i = nCoeffs - 2; i >= 0; i--)
124 retVal = coeffs[i] + x * retVal;
125 retVal = retVal * std::pow(x, lowestOrder);
126 return retVal + (pdfMode && lowestOrder > 0 ? 1.0 : 0.0);
127}
128
129inline double chebychevEvaluate(double *coeffs, unsigned int nCoeffs, double x_in, double xMin, double xMax)
130{
131 // transform to range [-1, +1]
132 const double xPrime = (x_in - 0.5 * (xMax + xMin)) / (0.5 * (xMax - xMin));
133
134 // extract current values of coefficients
135 double sum = 1.;
136 if (nCoeffs > 0) {
137 double curr = xPrime;
138 double twox = 2 * xPrime;
139 double last = 1;
140 double newval = twox * curr - last;
141 last = curr;
142 curr = newval;
143 for (unsigned int i = 0; nCoeffs != i; ++i) {
144 sum += last * coeffs[i];
145 newval = twox * curr - last;
146 last = curr;
147 curr = newval;
148 }
149 }
150 return sum;
151}
152
153inline double constraintSumEvaluate(double const *comp, unsigned int compSize)
154{
155 double sum = 0;
156 for (unsigned int i = 0; i < compSize; i++) {
157 sum -= std::log(comp[i]);
158 }
159 return sum;
160}
161
162inline unsigned int getUniformBinning(double low, double high, double val, unsigned int numBins)
163{
164 double binWidth = (high - low) / numBins;
165 return val >= high ? numBins - 1 : std::abs((val - low) / binWidth);
166}
167
168inline double poissonEvaluate(double x, double par)
169{
170 if (par < 0)
171 return TMath::QuietNaN();
172
173 if (x < 0) {
174 return 0;
175 } else if (x == 0.0) {
176 return std::exp(-par);
177 } else {
178 double out = x * std::log(par) - TMath::LnGamma(x + 1.) - par;
179 return std::exp(out);
180 }
181}
182
183inline double interpolate6thDegree(double x, double low, double high, double nominal, double boundary)
184{
185 double t = x / boundary;
186 double eps_plus = high - nominal;
187 double eps_minus = nominal - low;
188 double S = 0.5 * (eps_plus + eps_minus);
189 double A = 0.0625 * (eps_plus - eps_minus);
190
191 return x * (S + t * A * (15 + t * t * (-10 + t * t * 3)));
192}
193
194inline double interpolate6thDegreeExp(double x, double low, double high, double nominal, double boundary)
195{
196 double x0 = boundary;
197
198 // GHL: Swagato's suggestions
199 double powUp = std::pow(high / nominal, x0);
200 double powDown = std::pow(low / nominal, x0);
201 double logHi = std::log(high);
202 double logLo = std::log(low);
203 double powUpLog = high <= 0.0 ? 0.0 : powUp * logHi;
204 double powDownLog = low <= 0.0 ? 0.0 : -powDown * logLo;
205 double powUpLog2 = high <= 0.0 ? 0.0 : powUpLog * logHi;
206 double powDownLog2 = low <= 0.0 ? 0.0 : -powDownLog * logLo;
207
208 double S0 = 0.5 * (powUp + powDown);
209 double A0 = 0.5 * (powUp - powDown);
210 double S1 = 0.5 * (powUpLog + powDownLog);
211 double A1 = 0.5 * (powUpLog - powDownLog);
212 double S2 = 0.5 * (powUpLog2 + powDownLog2);
213 double A2 = 0.5 * (powUpLog2 - powDownLog2);
214
215 // fcns+der+2nd_der are eq at bd
216
217 double a = 1. / (8 * x0) * (15 * A0 - 7 * x0 * S1 + x0 * x0 * A2);
218 double b = 1. / (8 * x0 * x0) * (-24 + 24 * S0 - 9 * x0 * A1 + x0 * x0 * S2);
219 double c = 1. / (4 * std::pow(x0, 3)) * (-5 * A0 + 5 * x0 * S1 - x0 * x0 * A2);
220 double d = 1. / (4 * std::pow(x0, 4)) * (12 - 12 * S0 + 7 * x0 * A1 - x0 * x0 * S2);
221 double e = 1. / (8 * std::pow(x0, 5)) * (+3 * A0 - 3 * x0 * S1 + x0 * x0 * A2);
222 double f = 1. / (8 * std::pow(x0, 6)) * (-8 + 8 * S0 - 5 * x0 * A1 + x0 * x0 * S2);
223
224 // evaluate the 6-th degree polynomial using Horner's method
225 double value = 1. + x * (a + x * (b + x * (c + x * (d + x * (e + x * f)))));
226 return value;
227}
228
229inline double
230flexibleInterp(unsigned int code, double low, double high, double boundary, double nominal, double paramVal, double res)
231{
232 if (code == 0) {
233 // piece-wise linear
234 if (paramVal > 0) {
235 return paramVal * (high - nominal);
236 } else {
237 return paramVal * (nominal - low);
238 }
239 } else if (code == 1) {
240 // piece-wise log
241 if (paramVal >= 0) {
242 return res * (std::pow(high / nominal, +paramVal) - 1);
243 } else {
244 return res * (std::pow(low / nominal, -paramVal) - 1);
245 }
246 } else if (code == 2) {
247 // parabolic with linear
248 double a = 0.5 * (high + low) - nominal;
249 double b = 0.5 * (high - low);
250 double c = 0;
251 if (paramVal > 1) {
252 return (2 * a + b) * (paramVal - 1) + high - nominal;
253 } else if (paramVal < -1) {
254 return -1 * (2 * a - b) * (paramVal + 1) + low - nominal;
255 } else {
256 return a * std::pow(paramVal, 2) + b * paramVal + c;
257 }
258 } else if (code == 3) {
259 // parabolic version of log-normal
260 double a = 0.5 * (high + low) - nominal;
261 double b = 0.5 * (high - low);
262 double c = 0;
263 if (paramVal > 1) {
264 return (2 * a + b) * (paramVal - 1) + high - nominal;
265 } else if (paramVal < -1) {
266 return -1 * (2 * a - b) * (paramVal + 1) + low - nominal;
267 } else {
268 return a * std::pow(paramVal, 2) + b * paramVal + c;
269 }
270 } else if (code == 4) {
271 double x = paramVal;
272 if (x >= boundary) {
273 return x * (high - nominal);
274 } else if (x <= -boundary) {
275 return x * (nominal - low);
276 }
277
278 return interpolate6thDegree(x, low, high, nominal, boundary);
279 } else if (code == 5) {
280 double x = paramVal;
281 double mod = 1.0;
282 if (x >= boundary) {
283 mod = std::pow(high / nominal, +paramVal);
284 } else if (x <= -boundary) {
285 mod = std::pow(low / nominal, -paramVal);
286 } else {
287 mod = interpolate6thDegreeExp(x, low, high, nominal, boundary);
288 }
289 return res * (mod - 1.0);
290 }
291
292 return 0.0;
293}
294
295inline double flexibleInterpEvaluate(unsigned int code, double *params, unsigned int n, double *low, double *high,
296 double boundary, double nominal)
297{
298 double total = nominal;
299 for (std::size_t i = 0; i < n; ++i) {
300 total += flexibleInterp(code, low[i], high[i], boundary, nominal, params[i], total);
301 }
302
304}
305
306inline double piecewiseInterpolationEvaluate(unsigned int code, double *low, double *high, double nominal,
307 double *params, unsigned int n)
308{
309 double total = nominal;
310 for (std::size_t i = 0; i < n; ++i) {
311 total += flexibleInterp(code, low[i], high[i], 1.0, nominal, params[i], total);
312 }
313
315}
316
317inline double landauEvaluate(double x, double mu, double sigma)
318{
319 if (sigma <= 0.)
320 return 0.;
321 return ROOT::Math::landau_pdf((x - mu) / sigma);
322}
323
324inline double logNormalEvaluate(double x, double k, double m0)
325{
326 return ROOT::Math::lognormal_pdf(x, std::log(m0), std::abs(std::log(k)));
327}
328
329inline double logNormalEvaluateStandard(double x, double sigma, double mu)
330{
331 return ROOT::Math::lognormal_pdf(x, mu, std::abs(sigma));
332}
333
334inline double effProdEvaluate(double eff, double pdf)
335{
336 return eff * pdf;
337}
338
339inline double nllEvaluate(double pdf, double weight, int binnedL, int doBinOffset)
340{
341 if (binnedL) {
342 // Special handling of this case since std::log(Poisson(0,0)=0 but can't be
343 // calculated with usual log-formula since std::log(mu)=0. No update of result
344 // is required since term=0.
345 if (std::abs(pdf) < 1e-10 && std::abs(weight) < 1e-10) {
346 return 0.0;
347 }
348 if (doBinOffset) {
349 return pdf - weight - weight * (std::log(pdf) - std::log(weight));
350 }
351 return pdf - weight * std::log(pdf) + TMath::LnGamma(weight + 1);
352 } else {
353 return -weight * std::log(pdf);
354 }
355}
356
357inline double recursiveFractionEvaluate(double *a, unsigned int n)
358{
359 double prod = a[0];
360
361 for (unsigned int i = 1; i < n; ++i) {
362 prod *= 1.0 - a[i];
363 }
364
365 return prod;
366}
367
368inline double cbShapeEvaluate(double m, double m0, double sigma, double alpha, double n)
369{
370 double t = (m - m0) / sigma;
371 if (alpha < 0)
372 t = -t;
373
374 double absAlpha = std::abs((double)alpha);
375
376 if (t >= -absAlpha) {
377 return std::exp(-0.5 * t * t);
378 } else {
379 double a = std::pow(n / absAlpha, n) * std::exp(-0.5 * absAlpha * absAlpha);
380 double b = n / absAlpha - absAlpha;
381
382 return a / std::pow(b - t, n);
383 }
384}
385
386} // namespace EvaluateFuncs
387
388} // namespace Detail
389
390} // namespace RooFit
391
392#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:
const Double_t sigma
Double_t x[n]
Definition legend1.C:17
const Int_t n
Definition legend1.C:16
double poissonEvaluate(double x, double par)
double bifurGaussEvaluate(double x, double mean, double sigmaL, double sigmaR)
double constraintSumEvaluate(double const *comp, unsigned int compSize)
double gaussianEvaluate(double x, double mean, double sigma)
Function to evaluate an un-normalized RooGaussian.
double piecewiseInterpolationEvaluate(unsigned int code, double *low, double *high, double nominal, double *params, unsigned int n)
double flexibleInterp(unsigned int code, double low, double high, double boundary, double nominal, double paramVal, double res)
double landauEvaluate(double x, double mu, double sigma)
double effProdEvaluate(double eff, double pdf)
double polynomialEvaluate(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 flexibleInterpEvaluate(unsigned int code, double *params, unsigned int n, double *low, double *high, double boundary, double nominal)
double recursiveFractionEvaluate(double *a, unsigned int n)
double logNormalEvaluate(double x, double k, double m0)
double logNormalEvaluateStandard(double x, double sigma, double mu)
double efficiencyEvaluate(double effFuncVal, int catIndex, int sigCatIndex)
double interpolate6thDegreeExp(double x, double low, double high, double nominal, double boundary)
double bernsteinEvaluate(double x, double xmin, double xmax, double *coefs, int nCoefs)
The caller needs to make sure that there is at least one coefficient.
unsigned int getUniformBinning(double low, double high, double val, unsigned int numBins)
double cbShapeEvaluate(double m, double m0, double sigma, double alpha, double n)
double ratioEvaluate(double numerator, double denominator)
double interpolate6thDegree(double x, double low, double high, double nominal, double boundary)
double nllEvaluate(double pdf, double weight, int binnedL, int doBinOffset)
double chebychevEvaluate(double *coeffs, unsigned int nCoeffs, double x_in, double xMin, double xMax)
double binomial(int n, int k)
Calculates the binomial coefficient n over k.
The namespace RooFit contains mostly switches that change the behaviour of functions of PDFs (or othe...
Definition JSONIO.h:26
Double_t QuietNaN()
Returns a quiet NaN as defined by IEEE 754.
Definition TMath.h:902
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
TMarker m
Definition textangle.C:8
static uint64_t sum(uint64_t i)
Definition Factory.cxx:2345