Logo ROOT  
Reference Guide
 
Loading...
Searching...
No Matches
RooPolyFunc.cxx
Go to the documentation of this file.
1/*****************************************************************************
2 * Project: RooFit *
3 * Package: RooFitModels *
4 * @(#)root/roofit:$Id$
5 * Authors: *
6 * WV, Wouter Verkerke, UC Santa Barbara, verkerke@slac.stanford.edu *
7 * DK, David Kirkby, UC Irvine, dkirkby@uci.edu *
8 * *
9 * Copyright (c) 2021: *
10 * CERN, Switzerland *
11 * *
12 * Redistribution and use in source and binary forms, *
13 * with or without modification, are permitted according to the terms *
14 * listed in LICENSE (http://roofit.sourceforge.net/license.txt) *
15 *****************************************************************************/
16
17/** \class RooPolyFunc
18 \ingroup Roofit
19
20RooPolyFunc implements a polynomial function in multi-variables.
21The polynomial coefficients are implemented as doubles and are not
22part of the RooFit computation graph.
23**/
24
25#include "RooPolyFunc.h"
26
27#include "RooAbsReal.h"
28#include "RooArgList.h"
29#include "RooArgSet.h"
30#include "RooDerivative.h"
31#include "RooMsgService.h"
32#include "RooRealVar.h"
33
34using namespace std;
35using namespace RooFit;
36
38
39////////////////////////////////////////////////////////////////////////////////
40/// coverity[UNINIT_CTOR]
41
42void RooPolyFunc::addTerm(double coefficient)
43{
44 int n_terms = _terms.size();
45 std::string coeff_name = Form("%s_c%d", GetName(), n_terms);
46 std::string term_name = Form("%s_t%d", GetName(), n_terms);
47 auto termList = std::make_unique<RooListProxy>(term_name.c_str(), term_name.c_str(), this);
48 auto coeff = new RooRealVar(coeff_name.c_str(), coeff_name.c_str(), coefficient);
49 auto exponents = new RooArgList();
50
51 for (const auto &var : _vars) {
52 std::string exponent_name = Form("%s_%s^%d", GetName(), var->GetName(), 0);
53 auto exponent = new RooRealVar(exponent_name.c_str(), exponent_name.c_str(), 0);
54 exponents->add(*exponent);
55 }
56
57 termList->addOwned(*exponents);
58 termList->addOwned(*coeff);
59 this->_terms.push_back(move(termList));
60}
61
62void RooPolyFunc::addTerm(double coefficient, const RooAbsReal &var1, int exp1)
63{
64 int n_terms = _terms.size();
65 std::string coeff_name = Form("%s_c%d", GetName(), n_terms);
66 std::string term_name = Form("%s_t%d", GetName(), n_terms);
67 auto termList = std::make_unique<RooListProxy>(term_name.c_str(), term_name.c_str(), this);
68 auto coeff = new RooRealVar(coeff_name.c_str(), coeff_name.c_str(), coefficient);
69 auto exponents = new RooArgList();
70
71 // linear iterate over all the variables, create var1^exp1 ..vark^0
72 for (const auto &var : _vars) {
73 int exp = 0;
74 if (strcmp(var1.GetName(), var->GetName()) == 0)
75 exp += exp1;
76 std::string exponent_name = Form("%s_%s^%d", GetName(), var->GetName(), exp);
77 auto exponent = new RooRealVar(exponent_name.c_str(), exponent_name.c_str(), exp);
78 exponents->add(*exponent);
79 }
80
81 termList->addOwned(*exponents);
82 termList->addOwned(*coeff);
83 this->_terms.push_back(move(termList));
84}
85
86void RooPolyFunc::addTerm(double coefficient, const RooAbsReal &var1, int exp1, const RooAbsReal &var2, int exp2)
87{
88
89 int n_terms = _terms.size();
90 std::string coeff_name = Form("%s_c%d", GetName(), n_terms);
91 std::string term_name = Form("%s_t%d", GetName(), n_terms);
92 auto termList = std::make_unique<RooListProxy>(term_name.c_str(), term_name.c_str(), this);
93 auto coeff = new RooRealVar(coeff_name.c_str(), coeff_name.c_str(), coefficient);
94 auto exponents = new RooArgList();
95
96 for (const auto &var : _vars) {
97 int exp = 0;
98 if (strcmp(var1.GetName(), var->GetName()) == 0)
99 exp += exp1;
100 if (strcmp(var2.GetName(), var->GetName()) == 0)
101 exp += exp2;
102 std::string exponent_name = Form("%s_%s^%d", GetName(), var->GetName(), exp);
103 auto exponent = new RooRealVar(exponent_name.c_str(), exponent_name.c_str(), exp);
104 exponents->add(*exponent);
105 }
106 termList->addOwned(*exponents);
107 termList->addOwned(*coeff);
108 this->_terms.push_back(move(termList));
109}
110
111void RooPolyFunc::addTerm(double coefficient, const RooAbsCollection &exponents)
112{
113 if (exponents.size() != _vars.size()) {
114 coutE(InputArguments) << "RooPolyFunc::addTerm(" << GetName() << ") WARNING: number of exponents ("
115 << exponents.size() << ") provided do not match the number of variables (" << _vars.size()
116 << ")" << endl;
117 }
118 int n_terms = _terms.size();
119 std::string coeff_name = Form("%s_c%d", GetName(), n_terms);
120 std::string term_name = Form("%s_t%d", GetName(), n_terms);
121 auto termList = std::make_unique<RooListProxy>(term_name.c_str(), term_name.c_str(), this);
122 auto coeff = new RooRealVar(coeff_name.c_str(), coeff_name.c_str(), coefficient);
123 termList->addOwned(exponents);
124 termList->addOwned(*coeff);
125 this->_terms.push_back(move(termList));
126}
127
128////////////////////////////////////////////////////////////////////////////////
129/// Default constructor
130
132
133////////////////////////////////////////////////////////////////////////////////
134/// Parameterised constructor
135
136RooPolyFunc::RooPolyFunc(const char *name, const char *title, const RooAbsCollection &vars)
137 : RooAbsReal(name, title), _vars("x", "list of dependent variables", this)
138{
139 for (const auto &var : vars) {
140 if (!dynamic_cast<RooAbsReal *>(var)) {
141 std::stringstream ss;
142 ss << "RooPolyFunc::ctor(" << GetName() << ") ERROR: coefficient " << var->GetName()
143 << " is not of type RooAbsReal";
144 const std::string errorMsg = ss.str();
145 coutE(InputArguments) << errorMsg << std::endl;
146 throw std::runtime_error(errorMsg);
147 }
148 _vars.add(*var);
149 }
150}
151
152////////////////////////////////////////////////////////////////////////////////
153/// Copy constructor
154
155RooPolyFunc::RooPolyFunc(const RooPolyFunc &other, const char *name)
156 : RooAbsReal(other, name), _vars("vars", this, other._vars)
157{
158 for (auto const &term : other._terms) {
159 this->_terms.emplace_back(std::make_unique<RooListProxy>(term->GetName(), this, *term));
160 }
161}
162
163////////////////////////////////////////////////////////////////////////////////
164/// Assignment operator
165
167{
169 _vars = other._vars;
170
171 for (auto const &term : other._terms) {
172 this->_terms.emplace_back(std::make_unique<RooListProxy>(term->GetName(), this, *term));
173 }
174 return *this;
175}
176
177////////////////////////////////////////////////////////////////////////////////
178/// Evaluate value of Polynomial.
180{
181 // Calculate and return value of polynomial
182 double poly_sum(0.0);
183 for (const auto &term : _terms) {
184 double poly_term(1.0);
185 size_t n_vars = term->size() - 1;
186 for (size_t i_var = 0; i_var < n_vars; ++i_var) {
187 auto var = dynamic_cast<RooRealVar *>(_vars.at(i_var));
188 auto exp = dynamic_cast<RooRealVar *>(term->at(i_var));
189 poly_term *= pow(var->getVal(), exp->getVal());
190 }
191 auto coef = dynamic_cast<RooRealVar *>(term->at(n_vars));
192 poly_sum += coef->getVal() * poly_term;
193 }
194 return poly_sum;
195}
196
197void setCoordinates(const RooAbsCollection &observables, std::vector<double> const &observableValues)
198// set all observables to expansion co-ordinate
199{
200 std::size_t iObs = 0;
201 for (auto *var : static_range_cast<RooRealVar *>(observables)) {
202 var->setVal(observableValues[iObs++]);
203 }
204}
205
206void fixObservables(const RooAbsCollection &observables)
207{
208 for (auto *var : static_range_cast<RooRealVar *>(observables)) {
209 var->setConstant(true);
210 }
211}
212////////////////////////////////////////////////////////////////////////////////
213///// Taylor expanding given function in terms of observables around
214///// observableValues. Supports expansions upto order 2.
215///// \param[in] function of variables that is taylor expanded.
216///// \param[in] observables set of variables to perform the expansion.
217///// \param[in] observableValues co-ordinates around which expansion is performed.
218///// \param[in] order order of the expansion (0,1,2 supported).
219///// \param[in] eps1 precision for first derivative and second derivative.
220///// \param[in] eps2 precision for second partial derivative of cross-derivative.
221std::unique_ptr<RooAbsReal>
222RooPolyFunc::taylorExpand(const char *name, const char *title, RooAbsReal &func, const RooAbsCollection &observables,
223 std::vector<double> const &observableValues, int order, double eps1, double eps2)
224{
225 // create the taylor expansion polynomial
226 auto taylor_poly = std::make_unique<RooPolyFunc>(name, title, observables);
227
228 // taylor expansion can be performed only for order 0, 1, 2 currently
229 if (order >= 3 || order <= 0) {
230 std::stringstream errorMsgStream;
231 errorMsgStream << "RooPolyFunc::taylorExpand(" << name << ") ERROR: order must be 0, 1, or 2";
232 const auto errorMsg = errorMsgStream.str();
233 oocoutE(taylor_poly.get(), InputArguments) << errorMsg << std::endl;
234 throw std::invalid_argument(errorMsg);
235 }
236
237 setCoordinates(observables, observableValues);
238
239 // estimate taylor expansion polynomial for different orders
240 // f(x) = f(x=x0)
241 // + \sum_i + \frac{df}{dx_i}|_{x_i=x_i^0}(x - x_i^0)
242 // + \sum_{i,j} 0.5 * \frac{d^f}{dx_i x_j}
243 // ( x_i - x_i^0 )( x_j - x_j^0 )
244 // note in the polynomial the () brackets are expanded out
245 for (int i_order = 0; i_order <= order; ++i_order) {
246 switch (i_order) {
247 case 0: {
248 taylor_poly->addTerm(func.getVal());
249 break;
250 }
251 case 1: {
252 for (auto *var : static_range_cast<RooRealVar *>(observables)) {
253 double var1_val = var->getVal();
254 auto deriv = func.derivative(*var, 1, eps2);
255 double deriv_val = deriv->getVal();
256 setCoordinates(observables, observableValues);
257 taylor_poly->addTerm(deriv_val, *var, 1);
258 if (var1_val != 0.0) {
259 taylor_poly->addTerm(deriv_val * var1_val * -1.0);
260 }
261 }
262 break;
263 }
264 case 2: {
265 for (auto *var1 : static_range_cast<RooRealVar *>(observables)) {
266 double var1_val = var1->getVal();
267 auto deriv1 = func.derivative(*var1, 1, eps1);
268 for (auto *var2 : static_range_cast<RooRealVar *>(observables)) {
269 double var2_val = var2->getVal();
270 double deriv_val = 0.0;
271 if (strcmp(var1->GetName(), var2->GetName()) == 0) {
272 auto deriv2 = func.derivative(*var2, 2, eps2);
273 deriv_val = 0.5 * deriv2->getVal();
274 setCoordinates(observables, observableValues);
275 } else {
276 auto deriv2 = deriv1->derivative(*var2, 1, eps2);
277 deriv_val = 0.5 * deriv2->getVal();
278 setCoordinates(observables, observableValues);
279 }
280 taylor_poly->addTerm(deriv_val, *var1, 1, *var2, 1);
281 if (var1_val != 0.0 || var2_val != 0.0) {
282 taylor_poly->addTerm(deriv_val * var1_val * var2_val);
283 taylor_poly->addTerm(deriv_val * var2_val * -1.0, *var1, 1);
284 taylor_poly->addTerm(deriv_val * var1_val * -1.0, *var2, 1);
285 }
286 }
287 }
288 break;
289 }
290 }
291 }
292 return taylor_poly;
293}
294
295////////////////////////////////////////////////////////////////////////////////
296/// Taylor expanding given function in terms of observables around
297/// defaultValue for all observables.
298std::unique_ptr<RooAbsReal> RooPolyFunc::taylorExpand(const char *name, const char *title, RooAbsReal &func,
299 const RooAbsCollection &observables, double observablesValue,
300 int order, double eps1, double eps2)
301{
302 return RooPolyFunc::taylorExpand(name, title, func, observables,
303 std::vector<double>(observables.size(), observablesValue), order, eps1, eps2);
304}
#define oocoutE(o, a)
#define coutE(a)
void fixObservables(const RooAbsCollection &observables)
void setCoordinates(const RooAbsCollection &observables, std::vector< double > const &observableValues)
#define ClassImp(name)
Definition Rtypes.h:364
char name[80]
Definition TGX11.cxx:110
char * Form(const char *fmt,...)
RooAbsCollection is an abstract container object that can hold multiple RooAbsArg objects.
Storage_t::size_type size() const
RooAbsReal is the common abstract base class for objects that represent a real value and implements f...
Definition RooAbsReal.h:64
RooDerivative * derivative(RooRealVar &obs, Int_t order=1, Double_t eps=0.001)
Return function representing first, second or third order derivative of this function.
RooAbsReal & operator=(const RooAbsReal &other)
Assign values, name and configs from another RooAbsReal.
Double_t getVal(const RooArgSet *normalisationSet=nullptr) const
Evaluate object.
Definition RooAbsReal.h:94
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
virtual Bool_t add(const RooAbsArg &var, Bool_t silent=kFALSE) override
Reimplementation of standard RooArgList::add()
RooPolyFunc implements a polynomial function in multi-variables.
Definition RooPolyFunc.h:27
static std::unique_ptr< RooAbsReal > taylorExpand(const char *name, const char *title, RooAbsReal &func, const RooAbsCollection &observables, std::vector< double > const &observableValues, int order=1, double eps1=1e-6, double eps2=1e-3)
RooListProxy _vars
Definition RooPolyFunc.h:51
RooPolyFunc()
Default constructor.
void addTerm(double coefficient)
coverity[UNINIT_CTOR]
RooPolyFunc & operator=(const RooPolyFunc &other)
Assignment operator.
std::vector< std::unique_ptr< RooListProxy > > _terms
Definition RooPolyFunc.h:52
double evaluate() const
Evaluation.
RooRealVar represents a variable that can be changed from the outside.
Definition RooRealVar.h:39
virtual const char * GetName() const
Returns name of object.
Definition TNamed.h:47
The namespace RooFit contains mostly switches that change the behaviour of functions of PDFs (or othe...
Definition Common.h:18
@ InputArguments