 ROOT   Reference Guide
Searching...
No Matches
RooPower.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
9 */
10
11/** \class RooPower
12 \ingroup Roofit
13
14RooPower implements a power law PDF of the form
15\f[ f(x) = \mathcal{N} \cdot \sum_{i} a_{i} * x^{b_i} \f]
16
17\image html RooPower.png
18**/
19
20#include <RooPower.h>
21
22#include <RooAbsReal.h>
23#include <RooArgList.h>
24#include <RooMsgService.h>
25
26#include <TError.h>
27
28#include <cmath>
29#include <cassert>
30#include <sstream>
31
33
34////////////////////////////////////////////////////////////////////////////////
35/// Create a power law in the variable x.
36/// \param[in] name Name of the PDF
37/// \param[in] title Title for plotting the PDF
38/// \param[in] x The variable of the polynomial
39/// \param[in] coefList The coefficients \f$a_i \f$
40/// \param[in] expList The exponentials \f$b_i \f$
41/// \f[
42/// \sum_{i=0}^{n} a_{i} * x^{b_{i}}
43/// \f]
44///
45/// This means that
46/// \code{.cpp}
47/// RooPower powl("pow", "pow", x, RooArgList(a1, a2), RooArgList(b1,b2))
48/// \endcode
49/// computes
50/// \f[
51/// \mathrm{pol}(x) = a1 * x^b1 + a2 * x^b2
52/// \f]
53
54RooPower::RooPower(const char *name, const char *title, RooAbsReal &x, const RooArgList &coefList,
55 const RooArgList &expList)
56 : RooAbsPdf(name, title),
57 _x("x", "Dependent", this, x),
58 _coefList("coefList", "List of coefficients", this),
59 _expList("expList", "List of exponents", this)
60{
61 if (coefList.size() != expList.size()) {
62 coutE(InputArguments) << "RooPower::ctor(" << GetName()
63 << ") ERROR: coefficient list and exponent list must be of same length" << std::endl;
64 return;
65 }
66 for (auto coef : coefList) {
67 if (!dynamic_cast<RooAbsReal *>(coef)) {
68 coutE(InputArguments) << "RooPower::ctor(" << GetName() << ") ERROR: coefficient " << coef->GetName()
69 << " is not of type RooAbsReal" << std::endl;
70 R__ASSERT(0);
71 }
73 }
74 for (auto exp : expList) {
75 if (!dynamic_cast<RooAbsReal *>(exp)) {
76 coutE(InputArguments) << "RooPower::ctor(" << GetName() << ") ERROR: coefficient " << exp->GetName()
77 << " is not of type RooAbsReal" << std::endl;
78 R__ASSERT(0);
79 }
81 }
82}
83
84////////////////////////////////////////////////////////////////////////////////
85/// Copy constructor
86
87RooPower::RooPower(const RooPower &other, const char *name)
88 : RooAbsPdf(other, name),
89 _x("x", this, other._x),
90 _coefList("coefList", this, other._coefList),
91 _expList("expList", this, other._expList)
92{
93}
94
95////////////////////////////////////////////////////////////////////////////////
96
97double RooPower::evaluate() const
98{
99 // Calculate and return value of polynomial
100
101 const unsigned sz = _coefList.size();
102 if (!sz) {
103 return 0.;
104 }
105
106 std::vector<double> coefs;
107 std::vector<double> exps;
108 coefs.reserve(sz);
109 exps.reserve(sz);
110 const RooArgSet *nset = _coefList.nset();
111 for (auto c : _coefList) {
112 coefs.push_back(static_cast<RooAbsReal *>(c)->getVal(nset));
113 }
114 for (auto c : _expList) {
115 exps.push_back(static_cast<RooAbsReal *>(c)->getVal(nset));
116 }
117 double x = this->_x;
118 double retval = 0;
119 for (unsigned int i = 0; i < sz; ++i) {
120 retval += coefs[i] * pow(x, exps[i]);
121 }
122 return retval;
123}
124
125////////////////////////////////////////////////////////////////////////////////
126/// Advertise to RooFit that this function can be analytically integrated.
127int RooPower::getAnalyticalIntegral(RooArgSet &allVars, RooArgSet &analVars, const char * /*rangeName*/) const
128{
129 if (matchArgs(allVars, analVars, _x))
130 return 1;
131 return 0;
132}
133
134////////////////////////////////////////////////////////////////////////////////
135/// Do the analytical integral according to the code that was returned by getAnalyticalIntegral().
136double RooPower::analyticalIntegral(int /*code*/, const char *rangeName) const
137{
138 const double xmin = _x.min(rangeName);
139 const double xmax = _x.max(rangeName);
140 const unsigned sz = _coefList.size();
141 if (!sz) {
142 return xmax - xmin;
143 }
144
145 std::vector<double> coefs;
146 std::vector<double> exps;
147 coefs.reserve(sz);
148 exps.reserve(sz);
149 const RooArgSet *nset = _coefList.nset();
150 for (auto c : _coefList) {
151 coefs.push_back(static_cast<RooAbsReal *>(c)->getVal(nset));
152 }
153 for (auto c : _expList) {
154 exps.push_back(static_cast<RooAbsReal *>(c)->getVal(nset));
155 }
156
157 double retval = 0;
158 for (unsigned int i = 0; i < sz; ++i) {
159 if (exps[i] == -1) {
160 retval += coefs[i] * (log(xmax) - log(xmin));
161 } else {
162 retval += coefs[i] / (exps[i] + 1) * (pow(xmax, (exps[i] + 1)) - pow(xmin, (exps[i] + 1)));
163 }
164 }
165 return retval;
166}
167
168std::string RooPower::getFormulaExpression(bool expand) const
169{
170 std::stringstream ss;
171 for (std::size_t i = 0; i < _coefList.size(); ++i) {
172 if (i != 0)
173 ss << "+";
174 if (expand)
175 ss << static_cast<RooAbsReal *>(_coefList.at(i))->getVal();
176 else
177 ss << _coefList.at(i)->GetName();
178 ss << "*pow(" << _x.GetName() << ",";
179 if (expand)
180 ss << static_cast<RooAbsReal *>(_expList.at(i))->getVal();
181 else
182 ss << _expList.at(i)->GetName();
183 ss << ")";
184 }
185 return ss.str().c_str();
186}
#define c(i)
Definition RSha256.hxx:101
#define coutE(a)
#define ClassImp(name)
Definition Rtypes.h:377
#define R__ASSERT(e)
Definition TError.h:118
char name
Definition TGX11.cxx:110
float xmin
float xmax
Storage_t::size_type size() const
const RooArgSet * nset() const
Definition RooAbsProxy.h:48
RooAbsReal is the common abstract base class for objects that represent a real value and implements f...
Definition RooAbsReal.h:61
double getVal(const RooArgSet *normalisationSet=nullptr) const
Evaluate object.
Definition RooAbsReal.h:105
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:55
bool add(const RooAbsArg &var, bool valueServer, bool shapeServer, bool silent)
Overloaded RooCollection_t::add() method insert object into set and registers object as server to own...
RooPower implements a power law PDF of the form.
Definition RooPower.h:20
RooListProxy _expList
Definition RooPower.h:45
std::string getFormulaExpression(bool expand) const
Definition RooPower.cxx:168
RooPower()
Definition RooPower.h:22
RooArgList const & coefList() const
Get the list of coefficients.
Definition RooPower.h:32
RooListProxy _coefList
Definition RooPower.h:44
RooRealProxy _x
Definition RooPower.h:43
int getAnalyticalIntegral(RooArgSet &allVars, RooArgSet &analVars, const char *rangeName=nullptr) const override
Advertise to RooFit that this function can be analytically integrated.
Definition RooPower.cxx:127
double analyticalIntegral(int code, const char *rangeName=nullptr) const override
Do the analytical integral according to the code that was returned by getAnalyticalIntegral().
Definition RooPower.cxx:136
double evaluate() const override
do not persist
Definition RooPower.cxx:97
RooArgList const & expList() const
Get the list of exponents.
Definition RooPower.h:35
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