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#include "RooBatchCompute.h"
26
27#include <TError.h>
28
29#include <cmath>
30#include <cassert>
31#include <sstream>
32
34
35////////////////////////////////////////////////////////////////////////////////
36/// Create a power law in the variable x.
37/// \param[in] name Name of the PDF
38/// \param[in] title Title for plotting the PDF
39/// \param[in] x The variable of the polynomial
40/// \param[in] coefList The coefficients \f$a_i \f$
41/// \param[in] expList The exponentials \f$b_i \f$
42/// \f[
43/// \sum_{i=0}^{n} a_{i} * x^{b_{i}}
44/// \f]
45///
46/// This means that
47/// \code{.cpp}
48/// RooPower powl("pow", "pow", x, RooArgList(a1, a2), RooArgList(b1,b2))
49/// \endcode
50/// computes
51/// \f[
52/// \mathrm{pol}(x) = a1 * x^b1 + a2 * x^b2
53/// \f]
54
55RooPower::RooPower(const char *name, const char *title, RooAbsReal &x, const RooArgList &coefList,
56 const RooArgList &expList)
57 : RooAbsPdf(name, title),
58 _x("x", "Dependent", this, x),
59 _coefList("coefList", "List of coefficients", this),
60 _expList("expList", "List of exponents", this)
61{
62 if (coefList.size() != expList.size()) {
63 coutE(InputArguments) << "RooPower::ctor(" << GetName()
64 << ") ERROR: coefficient list and exponent list must be of same length" << std::endl;
65 return;
66 }
67 for (auto coef : coefList) {
68 if (!dynamic_cast<RooAbsReal *>(coef)) {
69 coutE(InputArguments) << "RooPower::ctor(" << GetName() << ") ERROR: coefficient " << coef->GetName()
70 << " is not of type RooAbsReal" << std::endl;
71 R__ASSERT(0);
72 }
74 }
75 for (auto exp : expList) {
76 if (!dynamic_cast<RooAbsReal *>(exp)) {
77 coutE(InputArguments) << "RooPower::ctor(" << GetName() << ") ERROR: coefficient " << exp->GetName()
78 << " is not of type RooAbsReal" << std::endl;
79 R__ASSERT(0);
80 }
82 }
83}
84
85////////////////////////////////////////////////////////////////////////////////
86/// Copy constructor
87
88RooPower::RooPower(const RooPower &other, const char *name)
89 : RooAbsPdf(other, name),
90 _x("x", this, other._x),
91 _coefList("coefList", this, other._coefList),
92 _expList("expList", this, other._expList)
93{
94}
95
96////////////////////////////////////////////////////////////////////////////////
97
98/// Compute multiple values of Power distribution.
99void RooPower::computeBatch(double *output, size_t nEvents, RooFit::Detail::DataMap const &dataMap) const
100{
102 vars.reserve(2 * _coefList.size() + 1);
103 vars.push_back(dataMap.at(_x));
104
105 assert(_coefList.size() == _expList.size());
106
107 for (std::size_t i = 0; i < _coefList.size(); ++i) {
108 vars.push_back(dataMap.at(&_coefList[i]));
109 vars.push_back(dataMap.at(&_expList[i]));
110 }
111
113 args.push_back(_coefList.size());
114
115 RooBatchCompute::compute(dataMap.config(this), RooBatchCompute::Power, output, nEvents, vars, args);
116}
117
118////////////////////////////////////////////////////////////////////////////////
119
120double RooPower::evaluate() const
121{
122 // Calculate and return value of polynomial
123
124 const unsigned sz = _coefList.size();
125 if (!sz) {
126 return 0.;
127 }
128
129 std::vector<double> coefs;
130 std::vector<double> exps;
131 coefs.reserve(sz);
132 exps.reserve(sz);
133 const RooArgSet *nset = _coefList.nset();
134 for (auto c : _coefList) {
135 coefs.push_back(static_cast<RooAbsReal *>(c)->getVal(nset));
136 }
137 for (auto c : _expList) {
138 exps.push_back(static_cast<RooAbsReal *>(c)->getVal(nset));
139 }
140 double x = this->_x;
141 double retval = 0;
142 for (unsigned int i = 0; i < sz; ++i) {
143 retval += coefs[i] * std::pow(x, exps[i]);
144 }
145 return retval;
146}
147
148////////////////////////////////////////////////////////////////////////////////
149/// Advertise to RooFit that this function can be analytically integrated.
150int RooPower::getAnalyticalIntegral(RooArgSet &allVars, RooArgSet &analVars, const char * /*rangeName*/) const
151{
152 if (matchArgs(allVars, analVars, _x))
153 return 1;
154 return 0;
155}
156
157////////////////////////////////////////////////////////////////////////////////
158/// Do the analytical integral according to the code that was returned by getAnalyticalIntegral().
159double RooPower::analyticalIntegral(int /*code*/, const char *rangeName) const
160{
161 const double xmin = _x.min(rangeName);
162 const double xmax = _x.max(rangeName);
163 const unsigned sz = _coefList.size();
164 if (!sz) {
165 return xmax - xmin;
166 }
167
168 std::vector<double> coefs;
169 std::vector<double> exps;
170 coefs.reserve(sz);
171 exps.reserve(sz);
172 const RooArgSet *nset = _coefList.nset();
173 for (auto c : _coefList) {
174 coefs.push_back(static_cast<RooAbsReal *>(c)->getVal(nset));
175 }
176 for (auto c : _expList) {
177 exps.push_back(static_cast<RooAbsReal *>(c)->getVal(nset));
178 }
179
180 double retval = 0;
181 for (unsigned int i = 0; i < sz; ++i) {
182 if (exps[i] == -1) {
183 retval += coefs[i] * (log(xmax) - log(xmin));
184 } else {
185 retval += coefs[i] / (exps[i] + 1) * (pow(xmax, (exps[i] + 1)) - pow(xmin, (exps[i] + 1)));
186 }
187 }
188 return retval;
189}
190
191std::string RooPower::getFormulaExpression(bool expand) const
192{
193 std::stringstream ss;
194 for (std::size_t i = 0; i < _coefList.size(); ++i) {
195 if (i != 0)
196 ss << "+";
197 if (expand)
198 ss << static_cast<RooAbsReal *>(_coefList.at(i))->getVal();
199 else
200 ss << _coefList.at(i)->GetName();
201 ss << "*pow(" << _x.GetName() << ",";
202 if (expand)
203 ss << static_cast<RooAbsReal *>(_expList.at(i))->getVal();
204 else
205 ss << _expList.at(i)->GetName();
206 ss << ")";
207 }
208 return ss.str();
209}
#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[80]
Definition TGX11.cxx:110
float xmin
float xmax
Storage_t::size_type size() const
const RooArgSet * nset() const
Definition RooAbsProxy.h:52
RooAbsReal is the common abstract base class for objects that represent a real value and implements f...
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: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...
RooBatchCompute::Config config(RooAbsArg const *arg) const
Definition DataMap.cxx:40
std::span< const double > at(RooAbsArg const *arg, RooAbsArg const *caller=nullptr)
Definition DataMap.cxx:22
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:191
void computeBatch(double *output, size_t size, RooFit::Detail::DataMap const &) const override
do not persist
Definition RooPower.cxx:99
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:150
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:159
double evaluate() const override
Evaluation.
Definition RooPower.cxx:120
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
std::vector< std::span< const double > > VarVector
std::vector< double > ArgVector
void compute(Config cfg, Computer comp, RestrictArr output, size_t size, const VarVector &vars, ArgVector &extraArgs)
static void output()