Logo ROOT  
Reference Guide
 
All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Properties Friends Macros Modules Pages
Loading...
Searching...
No Matches
RooLegacyExpPoly.cxx
Go to the documentation of this file.
1/*
2 * Project: RooFit
3 *
4 * Copyright (c) 2022, CERN
5 *
6 * Redistribution and use in source and binary forms,
7 * with or without modification, are permitted according to the terms
8 * listed in LICENSE (http://roofit.sourceforge.net/license.txt)
9 */
10
11/** \class RooLegacyExpPoly
12 \ingroup Roofit
13
14RooLegacyExpPoly implements a polynomial PDF of the form \f[ f(x) =
15\mathcal{N} \cdot \exp( \sum_{i} a_{i} * x^{i} ) \f] \f$ \mathcal{N}
16\f$ is a normalisation constant that is automatically calculated when
17the function is used in computations.
18
19The sum can be truncated at the low end. See the main constructor
20RooLegacyExpPoly::RooLegacyExpPoly(const char*, const char*, RooAbsReal&, const RooArgList&, int)
21
22\image html RooLegacyExpPoly.png
23
24**/
25
26#include <RooLegacyExpPoly.h>
27
28#include <RooAbsReal.h>
29#include <RooArgList.h>
30#include <RooMath.h>
31#include <RooMsgService.h>
32#include <RooRealVar.h>
33#include "RooBatchCompute.h"
34
35#include <TMath.h>
36#include <TError.h>
37
38#include <array>
39#include <cassert>
40#include <cmath>
41#include <complex>
42#include <sstream>
43
44
45////////////////////////////////////////////////////////////////////////////////
46/// Create a polynomial in the variable `x`.
47/// \param[in] name Name of the PDF
48/// \param[in] title Title for plotting the PDF
49/// \param[in] x The variable of the polynomial
50/// \param[in] coefList The coefficients \f$ a_i \f$
51/// \param[in] lowestOrder [optional] Truncate the sum such that it skips the lower orders:
52/// \f[
53/// 1. + \sum_{i=0}^{\mathrm{coefList.size()}} a_{i} * x^{(i + \mathrm{lowestOrder})}
54/// \f]
55///
56/// This means that
57/// \code{.cpp}
58/// RooLegacyExpPoly pol("pol", "pol", x, RooArgList(a, b), lowestOrder = 2)
59/// \endcode
60/// computes
61/// \f[
62/// \mathrm{pol}(x) = 1 * x^0 + (0 * x^{\ldots}) + a * x^2 + b * x^3.
63/// \f]
64
65RooLegacyExpPoly::RooLegacyExpPoly(const char *name, const char *title, RooAbsReal &x, const RooArgList &coefList, int lowestOrder)
66 : RooAbsPdf(name, title),
67 _x("x", "Dependent", this, x),
68 _coefList("coefList", "List of coefficients", this),
69 _lowestOrder(lowestOrder)
70{
71 // Check lowest order
72 if (_lowestOrder < 0) {
73 coutE(InputArguments) << "RooLegacyExpPoly::ctor(" << GetName()
74 << ") WARNING: lowestOrder must be >=0, setting value to 0" << std::endl;
75 _lowestOrder = 0;
76 }
77
79}
80
81////////////////////////////////////////////////////////////////////////////////
82/// Copy constructor
83
86 _x("x", this, other._x),
87 _coefList("coefList", this, other._coefList),
88 _lowestOrder(other._lowestOrder)
89{
90}
91
92////////////////////////////////////////////////////////////////////////////////
93
95{
96 // Calculate and return value of polynomial
97
98 const unsigned sz = _coefList.size();
99 const int lowestOrder = _lowestOrder;
100 if (!sz)
101 return lowestOrder ? 1. : 0.;
102 std::vector<double> coefs;
103 coefs.reserve(sz);
104
105 const RooArgSet *nset = _coefList.nset();
106 for (auto coef : _coefList) {
107 coefs.push_back(static_cast<RooAbsReal *>(coef)->getVal(nset));
108 };
109 const double x = _x;
110 double xpow = std::pow(x, lowestOrder);
111 double retval = 0;
112 for (size_t i = 0; i < sz; ++i) {
113 retval += coefs[i] * xpow;
114 xpow *= x;
115 }
116
117 if (std::numeric_limits<double>::max_exponent < retval) {
118 coutE(InputArguments) << "RooLegacyExpPoly::evaluateLog(" << GetName() << ") ERROR: exponent at " << x
119 << " larger than allowed maximum, result will be infinite! " << retval << " > "
120 << std::numeric_limits<double>::max_exponent << " in " << this->getFormulaExpression(true)
121 << std::endl;
122 }
123 return retval;
124}
125
126////////////////////////////////////////////////////////////////////////////////
127
128/// Compute multiple values of ExpPoly distribution.
130{
131 std::vector<std::span<const double>> vars;
132 vars.reserve(_coefList.size() + 1);
133 vars.push_back(ctx.at(_x));
134
135 std::vector<double> coefVals;
136 for (RooAbsArg *coef : _coefList) {
137 vars.push_back(ctx.at(coef));
138 }
139
140 std::array<double, 2> args{static_cast<double>(_lowestOrder), static_cast<double>(_coefList.size())};
141
142 RooBatchCompute::compute(ctx.config(this), RooBatchCompute::ExpPoly, ctx.output(), vars, args);
143}
144
145////////////////////////////////////////////////////////////////////////////////
146
148{
149 // Adjust the limits of all the coefficients to reflect the numeric boundaries
150
151 const unsigned sz = _coefList.size();
152 double max = std::numeric_limits<double>::max_exponent / sz;
153 const int lowestOrder = _lowestOrder;
154 std::vector<double> coefs;
155 coefs.reserve(sz);
156
157 RooRealVar *x = dynamic_cast<RooRealVar *>(&(*_x));
158 if (x) {
159 const double xmax = x->getMax();
160 double xmaxpow = std::pow(xmax, lowestOrder);
161 for (size_t i = 0; i < sz; ++i) {
162 double thismax = max / xmaxpow;
163 RooRealVar *coef = dynamic_cast<RooRealVar *>(this->_coefList.at(i));
164 if (coef) {
165 coef->setVal(thismax);
166 coef->setMax(thismax);
167 }
168 xmaxpow *= xmax;
169 }
170 }
171}
172
173////////////////////////////////////////////////////////////////////////////////
174
176{
177 // Calculate and return value of function
178
179 const double logval = this->evaluateLog();
180 const double val = std::exp(logval);
181 if (std::isinf(val)) {
182 coutE(InputArguments) << "RooLegacyExpPoly::evaluate(" << GetName()
183 << ") ERROR: result of exponentiation is infinite! exponent was " << logval << std::endl;
184 }
185 return val;
186}
187
188////////////////////////////////////////////////////////////////////////////////
189
190double RooLegacyExpPoly::getLogVal(const RooArgSet *nset) const
191{
192 return RooAbsPdf::getLogVal(nset);
193 // return this->evaluateLog();
194}
195
196////////////////////////////////////////////////////////////////////////////////
197
198int RooLegacyExpPoly::getAnalyticalIntegral(RooArgSet &allVars, RooArgSet &analVars, const char * /*rangeName*/) const
199{
200
201 if ((_coefList.size() + _lowestOrder < 4) &&
202 ((_coefList.size() + _lowestOrder < 3) ||
203 (static_cast<RooAbsReal *>(_coefList.at(2 - _lowestOrder))->getVal() <= 0)) &&
204 matchArgs(allVars, analVars, _x)) {
205 return 0;
206 }
207 return 0;
208}
209
210////////////////////////////////////////////////////////////////////////////////
211
212#define PI TMath::Pi()
213
214namespace {
215double deltaerf(long double x1, long double x2)
216{
217 // several things happening here
218 // 1. use "erf(x) = -erf(-x)" to evaluate the function only ever for the positive side (higher precision)
219 // 2. use "erf(x) = 1.-erfc(x)" and, instead of "erf(x1) - erf(x2)", do "(1.-erfc(x1)) - (1.-erfc(x2)) = erfc(x2) -
220 // erfc(x1)"
221 double y2 = x1 > 0 ? erfc(x1) : -erfc(-x1);
222 double y1 = x2 > 0 ? erfc(x2) : -erfc(-x2);
223 if (y1 == y2) {
224 std::cout << "WARNING in calculation of analytical integral limited by numerical precision" << std::endl;
225 std::cout << "x: " << x1 << " , " << x2 << std::endl;
226 std::cout << "y: " << y1 << " , " << y2 << std::endl;
227 }
228 return y1 - y2;
229}
230
231double deltaerfi(double x1, double x2)
232{
233 std::complex<double> u1 = {x1, 0.};
234 std::complex<double> u2 = {x2, 0.};
235
236 std::complex<double> y2 = x2 > 0 ? RooMath::faddeeva(u2) : RooMath::faddeeva(-u2);
237 std::complex<double> y1 = x1 > 0 ? RooMath::faddeeva(u1) : RooMath::faddeeva(-u1);
238 if (y1 == y2) {
239 std::cout << "WARNING in calculation of analytical integral limited by numerical precision" << std::endl;
240 std::cout << "x: " << x1 << " , " << x2 << std::endl;
241 std::cout << "y: " << y1 << " , " << y2 << std::endl;
242 }
243 return y1.imag() - y2.imag();
244}
245} // namespace
246
247double RooLegacyExpPoly::analyticalIntegral(int /*code*/, const char *rangeName) const
248{
249 const double xmin = _x.min(rangeName);
250 const double xmax = _x.max(rangeName);
251 const unsigned sz = _coefList.size();
252 if (!sz)
253 return xmax - xmin;
254
255 std::vector<double> coefs;
256 coefs.reserve(sz);
257 const RooArgSet *nset = _coefList.nset();
258 for (auto c : _coefList) {
259 coefs.push_back(static_cast<RooAbsReal *>(c)->getVal(nset));
260 }
261
262 switch (_coefList.size() + _lowestOrder) {
263 case 1: return xmax - xmin;
264 case 2: {
265 const double a = coefs[1 - _lowestOrder];
266 if (a != 0) {
267 return 1. / a * (exp(a * xmax) - exp(a * xmin)) * (_lowestOrder == 0 ? exp(coefs[0]) : 1);
268 } else {
269 return xmax - xmin;
270 }
271 }
272 case 3: {
273 const double a = coefs[2 - _lowestOrder];
274 const double b = _lowestOrder == 2 ? 0. : coefs[1 - _lowestOrder];
275 const double c = _lowestOrder == 0 ? coefs[0] : 0.;
276 const double absa = std::abs(a);
277 const double sqrta = std::sqrt(absa);
278 if (a < 0) {
279 double d = ::deltaerf((-b + 2 * absa * xmax) / (2 * sqrta), (-b + 2 * absa * xmin) / (2 * sqrta));
280 double retval = exp(b * b / (4 * absa) + c) * std::sqrt(PI) * d / (2 * sqrta);
281 return retval;
282 } else if (a > 0) {
283 double d = ::deltaerfi((b + 2 * absa * xmax) / (2 * sqrta), (b + 2 * absa * xmin) / (2 * sqrta));
284 double retval = exp(-b * b / (4 * absa) + c) * std::sqrt(PI) * d / (2 * sqrta);
285 return retval;
286 } else if (b != 0) {
287 return 1. / b * (std::exp(b * xmax) - exp(b * xmin)) * exp(c);
288 } else {
289 return xmax - xmin;
290 }
291 }
292 }
293 return 0.;
294}
295
296////////////////////////////////////////////////////////////////////////////////
297
299{
300 std::stringstream ss;
301 ss << "exp(";
302 int order = _lowestOrder;
303 for (auto coef : _coefList) {
304 if (order != _lowestOrder)
305 ss << "+";
306 if (expand) {
308 } else {
309 ss << coef->GetName();
310 }
311 ss << "*pow(" << _x.GetName() << "," << order << ")";
312 ++order;
313 }
314 ss << ")";
315 return ss.str();
316}
#define d(i)
Definition RSha256.hxx:102
#define b(i)
Definition RSha256.hxx:100
#define c(i)
Definition RSha256.hxx:101
#define a(i)
Definition RSha256.hxx:99
#define coutE(a)
#define PI
ROOT::Detail::TRangeCast< T, true > TRangeDynCast
TRangeDynCast is an adapter class that allows the typed iteration through a TCollection.
Option_t Option_t TPoint TPoint const char x2
Option_t Option_t TPoint TPoint const char x1
Option_t Option_t TPoint TPoint const char y2
Option_t Option_t TPoint TPoint const char y1
char name[80]
Definition TGX11.cxx:110
float xmin
float xmax
Common abstract base class for objects that represent a value and a "shape" in RooFit.
Definition RooAbsArg.h:77
Storage_t::size_type size() const
bool addTyped(const RooAbsCollection &list, bool silent=false)
Adds elements of a given RooAbsCollection to the container if they match the specified type.
Abstract interface for all probability density functions.
Definition RooAbsPdf.h:40
virtual double getLogVal(const RooArgSet *set=nullptr) const
Return the log of the current value with given normalization An error message is printed if the argum...
const RooArgSet * nset() const
Definition RooAbsProxy.h:52
Abstract base class for objects that represent a real value and implements functionality common to al...
Definition RooAbsReal.h:59
double getVal(const RooArgSet *normalisationSet=nullptr) const
Evaluate object.
Definition RooAbsReal.h:103
bool matchArgs(const RooArgSet &allDeps, RooArgSet &numDeps, const RooArgProxy &a) const
Utility function for use in getAnalyticalIntegral().
RooArgList is a container object that can hold multiple RooAbsArg objects.
Definition RooArgList.h:22
RooAbsArg * at(Int_t idx) const
Return object at given index, or nullptr if index is out of range.
Definition RooArgList.h:110
RooArgSet is a container object that can hold multiple RooAbsArg objects.
Definition RooArgSet.h:24
std::span< const double > at(RooAbsArg const *arg, RooAbsArg const *caller=nullptr)
std::span< double > output()
RooBatchCompute::Config config(RooAbsArg const *arg) const
RooLegacyExpPoly implements a polynomial PDF of the form.
int getAnalyticalIntegral(RooArgSet &allVars, RooArgSet &analVars, const char *rangeName=nullptr) const override
Interface function getAnalyticalIntergral advertises the analytical integrals that are supported.
double evaluateLog() const
RooArgList const & coefList() const
Get the coefficient list.
RooListProxy _coefList
RooAbsReal const & x() const
Get the x variable.
double getLogVal(const RooArgSet *nset) const override
Return the log of the current value with given normalization An error message is printed if the argum...
std::string getFormulaExpression(bool expand) const
double evaluate() const override
Evaluation.
void doEval(RooFit::EvalContext &) const override
Compute multiple values of ExpPoly distribution.
double analyticalIntegral(int code, const char *rangeName=nullptr) const override
Implements the actual analytical integral(s) advertised by getAnalyticalIntegral.
int lowestOrder() const
Return the order for the first coefficient in the list.
static std::complex< double > faddeeva(std::complex< double > z)
evaluate Faddeeva function for complex argument
Definition RooMath.cxx:30
Variable that can be changed from the outside.
Definition RooRealVar.h:37
void setVal(double value) override
Set value of variable to 'value'.
void setMax(const char *name, double value)
Set maximum of name range to given value.
double max(const char *rname=nullptr) const
Query upper limit of range. This requires the payload to be RooAbsRealLValue or derived.
double min(const char *rname=nullptr) const
Query lower limit of range. This requires the payload to be RooAbsRealLValue or derived.
const char * GetName() const override
Returns name of object.
Definition TNamed.h:47
Double_t x[n]
Definition legend1.C:17
void compute(Config cfg, Computer comp, std::span< double > output, VarSpan vars, ArgSpan extraArgs={})