Logo ROOT  
Reference Guide
 
Loading...
Searching...
No Matches
RooPowerSum.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 RooPowerSum
12 \ingroup Roofit
13
14RooPowerSum 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 RooPowerSum.png
18**/
19
20#include <RooPowerSum.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 <array>
30#include <cassert>
31#include <cmath>
32#include <sstream>
33
35
36////////////////////////////////////////////////////////////////////////////////
37/// Create a power law in the variable `x`.
38/// \param[in] name Name of the PDF
39/// \param[in] title Title for plotting the PDF
40/// \param[in] x The variable of the polynomial
41/// \param[in] coefList The coefficients \f$ a_i \f$
42/// \param[in] expList The exponentials \f$ b_i \f$
43/// \f[
44/// \sum_{i=0}^{n} a_{i} * x^{b_{i}}
45/// \f]
46///
47/// This means that
48/// \code{.cpp}
49/// RooPowerSum powl("pow", "pow", x, RooArgList(a1, a2), RooArgList(b1,b2))
50/// \endcode
51/// computes
52/// \f[
53/// \mathrm{pol}(x) = a1 * x^b1 + a2 * x^b2
54/// \f]
55
56RooPowerSum::RooPowerSum(const char *name, const char *title, RooAbsReal &x, const RooArgList &coefList,
57 const RooArgList &expList)
58 : RooAbsPdf(name, title),
59 _x("x", "Dependent", this, x),
60 _coefList("coefList", "List of coefficients", this),
61 _expList("expList", "List of exponents", this)
62{
63 if (coefList.size() != expList.size()) {
64 coutE(InputArguments) << "RooPowerSum::ctor(" << GetName()
65 << ") ERROR: coefficient list and exponent list must be of same length" << std::endl;
66 return;
67 }
70}
71
72////////////////////////////////////////////////////////////////////////////////
73/// Copy constructor
74
75RooPowerSum::RooPowerSum(const RooPowerSum &other, const char *name)
76 : RooAbsPdf(other, name),
77 _x("x", this, other._x),
78 _coefList("coefList", this, other._coefList),
79 _expList("expList", this, other._expList)
80{
81}
82
83////////////////////////////////////////////////////////////////////////////////
84
85/// Compute multiple values of Power distribution.
87{
88 std::vector<std::span<const double>> vars;
89 vars.reserve(2 * _coefList.size() + 1);
90 vars.push_back(ctx.at(_x));
91
92 assert(_coefList.size() == _expList.size());
93
94 for (std::size_t i = 0; i < _coefList.size(); ++i) {
95 vars.push_back(ctx.at(&_coefList[i]));
96 vars.push_back(ctx.at(&_expList[i]));
97 }
98
99 std::array<double, 1> args{static_cast<double>(_coefList.size())};
100
101 RooBatchCompute::compute(ctx.config(this), RooBatchCompute::Power, ctx.output(), vars, args);
102}
103
104////////////////////////////////////////////////////////////////////////////////
105
107{
108 // Calculate and return value of polynomial
109
110 const unsigned sz = _coefList.size();
111 if (!sz) {
112 return 0.;
113 }
114
115 std::vector<double> coefs;
116 std::vector<double> exps;
117 coefs.reserve(sz);
118 exps.reserve(sz);
119 const RooArgSet *nset = _coefList.nset();
120 for (auto c : _coefList) {
121 coefs.push_back(static_cast<RooAbsReal *>(c)->getVal(nset));
122 }
123 for (auto c : _expList) {
124 exps.push_back(static_cast<RooAbsReal *>(c)->getVal(nset));
125 }
126 double x = this->_x;
127 double retval = 0;
128 for (unsigned int i = 0; i < sz; ++i) {
129 retval += coefs[i] * std::pow(x, exps[i]);
130 }
131 return retval;
132}
133
134////////////////////////////////////////////////////////////////////////////////
135/// Advertise to RooFit that this function can be analytically integrated.
136int RooPowerSum::getAnalyticalIntegral(RooArgSet &allVars, RooArgSet &analVars, const char * /*rangeName*/) const
137{
138 if (matchArgs(allVars, analVars, _x))
139 return 1;
140 return 0;
141}
142
143////////////////////////////////////////////////////////////////////////////////
144/// Do the analytical integral according to the code that was returned by getAnalyticalIntegral().
145double RooPowerSum::analyticalIntegral(int /*code*/, const char *rangeName) const
146{
147 const double xmin = _x.min(rangeName);
148 const double xmax = _x.max(rangeName);
149 const unsigned sz = _coefList.size();
150 if (!sz) {
151 return xmax - xmin;
152 }
153
154 std::vector<double> coefs;
155 std::vector<double> exps;
156 coefs.reserve(sz);
157 exps.reserve(sz);
158 const RooArgSet *nset = _coefList.nset();
159 for (auto c : _coefList) {
160 coefs.push_back(static_cast<RooAbsReal *>(c)->getVal(nset));
161 }
162 for (auto c : _expList) {
163 exps.push_back(static_cast<RooAbsReal *>(c)->getVal(nset));
164 }
165
166 double retval = 0;
167 for (unsigned int i = 0; i < sz; ++i) {
168 if (exps[i] == -1) {
169 retval += coefs[i] * (log(xmax) - log(xmin));
170 } else {
171 retval += coefs[i] / (exps[i] + 1) * (pow(xmax, (exps[i] + 1)) - pow(xmin, (exps[i] + 1)));
172 }
173 }
174 return retval;
175}
176
177std::string RooPowerSum::getFormulaExpression(bool expand) const
178{
179 std::stringstream ss;
180 for (std::size_t i = 0; i < _coefList.size(); ++i) {
181 if (i != 0)
182 ss << "+";
183 if (expand) {
184 ss << static_cast<RooAbsReal *>(_coefList.at(i))->getVal();
185 } else {
186 ss << _coefList.at(i)->GetName();
187 }
188 ss << "*pow(" << _x.GetName() << ",";
189 if (expand) {
190 ss << static_cast<RooAbsReal *>(_expList.at(i))->getVal();
191 } else {
192 ss << _expList.at(i)->GetName();
193 }
194 ss << ")";
195 }
196 return ss.str();
197}
#define c(i)
Definition RSha256.hxx:101
#define coutE(a)
#define ClassImp(name)
Definition Rtypes.h:382
char name[80]
Definition TGX11.cxx:110
float xmin
float xmax
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
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
RooPowerSum implements a power law PDF of the form.
Definition RooPowerSum.h:20
RooArgList const & coefList() const
Get the list of coefficients.
Definition RooPowerSum.h:32
double evaluate() const override
Evaluation.
int getAnalyticalIntegral(RooArgSet &allVars, RooArgSet &analVars, const char *rangeName=nullptr) const override
Advertise to RooFit that this function can be analytically integrated.
double analyticalIntegral(int code, const char *rangeName=nullptr) const override
Do the analytical integral according to the code that was returned by getAnalyticalIntegral().
RooRealProxy _x
Definition RooPowerSum.h:43
void doEval(RooFit::EvalContext &) const override
do not persist
RooListProxy _expList
Definition RooPowerSum.h:45
RooArgList const & expList() const
Get the list of exponents.
Definition RooPowerSum.h:35
std::string getFormulaExpression(bool expand) const
RooListProxy _coefList
Definition RooPowerSum.h:44
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={})