Logo ROOT   6.18/05
Reference Guide
RooBernstein.cxx
Go to the documentation of this file.
1/*****************************************************************************
2 * Project: RooFit *
3 * Package: RooFitModels *
4 * @(#)root/roofit:$Id$
5 * Authors: *
6 * Kyle Cranmer *
7 * *
8 *****************************************************************************/
9
10/** \class RooBernstein
11 \ingroup Roofit
12
13Bernstein basis polynomials are positive-definite in the range [0,1].
14In this implementation, we extend [0,1] to be the range of the parameter.
15There are n+1 Bernstein basis polynomials of degree n:
16\f[
17 B_{i,n}(x) = \begin{pmatrix}n \\\ i \end{pmatrix} x^i \cdot (1-x)^{n-i}
18\f]
19Thus, by providing n coefficients that are positive-definite, there
20is a natural way to have well-behaved polynomial PDFs. For any n, the n+1 polynomials
21'form a partition of unity', i.e., they sum to one for all values of x.
22They can be used as a basis to span the space of polynomials with degree n or less:
23\f[
24 PDF(x, c_0, ..., c_n) = \mathcal{N} \cdot \sum_{i=0}^{n} c_i \cdot B_{i,n}(x).
25\f]
26By giving n+1 coefficients in the constructor, this class constructs the n+1
27polynomials of degree n, and sums them to form an element of the space of polynomials
28of degree n. \f$ \mathcal{N} \f$ is a normalisation constant that takes care of the
29cases where the \f$ c_i \f$ are not all equal to one.
30
31See also
32http://www.idav.ucdavis.edu/education/CAGDNotes/Bernstein-Polynomials.pdf
33**/
34
35#include "RooFit.h"
36
37#include "Riostream.h"
38#include "Riostream.h"
39#include <math.h>
40#include "TMath.h"
41#include "RooBernstein.h"
42#include "RooAbsReal.h"
43#include "RooRealVar.h"
44#include "RooArgList.h"
45
46#include "TError.h"
47
48using namespace std;
49
51
52////////////////////////////////////////////////////////////////////////////////
53
55{
56}
57
58////////////////////////////////////////////////////////////////////////////////
59/// Constructor
60
61RooBernstein::RooBernstein(const char* name, const char* title,
62 RooAbsReal& x, const RooArgList& coefList):
63 RooAbsPdf(name, title),
64 _x("x", "Dependent", this, x),
65 _coefList("coefficients","List of coefficients",this)
66{
67 TIterator* coefIter = coefList.createIterator() ;
68 RooAbsArg* coef ;
69 while((coef = (RooAbsArg*)coefIter->Next())) {
70 if (!dynamic_cast<RooAbsReal*>(coef)) {
71 cout << "RooBernstein::ctor(" << GetName() << ") ERROR: coefficient " << coef->GetName()
72 << " is not of type RooAbsReal" << endl ;
73 R__ASSERT(0) ;
74 }
75 _coefList.add(*coef) ;
76 }
77 delete coefIter ;
78}
79
80////////////////////////////////////////////////////////////////////////////////
81
82RooBernstein::RooBernstein(const RooBernstein& other, const char* name) :
83 RooAbsPdf(other, name),
84 _x("x", this, other._x),
85 _coefList("coefList",this,other._coefList)
86{
87}
88
89////////////////////////////////////////////////////////////////////////////////
90
92{
93 Double_t xmin = _x.min();
94 Double_t x = (_x - xmin) / (_x.max() - xmin); // rescale to [0,1]
95 Int_t degree = _coefList.getSize() - 1; // n+1 polys of degree n
97
98 if(degree == 0) {
99
100 return ((RooAbsReal *)iter.next())->getVal();
101
102 } else if(degree == 1) {
103
104 Double_t a0 = ((RooAbsReal *)iter.next())->getVal(); // c0
105 Double_t a1 = ((RooAbsReal *)iter.next())->getVal() - a0; // c1 - c0
106 return a1 * x + a0;
107
108 } else if(degree == 2) {
109
110 Double_t a0 = ((RooAbsReal *)iter.next())->getVal(); // c0
111 Double_t a1 = 2 * (((RooAbsReal *)iter.next())->getVal() - a0); // 2 * (c1 - c0)
112 Double_t a2 = ((RooAbsReal *)iter.next())->getVal() - a1 - a0; // c0 - 2 * c1 + c2
113 return (a2 * x + a1) * x + a0;
114
115 } else if(degree > 2) {
116
117 Double_t t = x;
118 Double_t s = 1 - x;
119
120 Double_t result = ((RooAbsReal *)iter.next())->getVal() * s;
121 for(Int_t i = 1; i < degree; i++) {
122 result = (result + t * TMath::Binomial(degree, i) * ((RooAbsReal *)iter.next())->getVal()) * s;
123 t *= x;
124 }
125 result += t * ((RooAbsReal *)iter.next())->getVal();
126
127 return result;
128 }
129
130 // in case list of arguments passed is empty
131 return TMath::SignalingNaN();
132}
133
134////////////////////////////////////////////////////////////////////////////////
135/// No analytical calculation available (yet) of integrals over subranges
136
137Int_t RooBernstein::getAnalyticalIntegral(RooArgSet& allVars, RooArgSet& analVars, const char* rangeName) const
138{
139 if (rangeName && strlen(rangeName)) {
140 return 0 ;
141 }
142
143 if (matchArgs(allVars, analVars, _x)) return 1;
144 return 0;
145}
146
147////////////////////////////////////////////////////////////////////////////////
148
149Double_t RooBernstein::analyticalIntegral(Int_t code, const char* rangeName) const
150{
151 R__ASSERT(code==1) ;
152 Double_t xmin = _x.min(rangeName); Double_t xmax = _x.max(rangeName);
153 Int_t degree= _coefList.getSize()-1; // n+1 polys of degree n
154 Double_t norm(0) ;
155
157 Double_t temp=0;
158 for (int i=0; i<=degree; ++i){
159 // for each of the i Bernstein basis polynomials
160 // represent it in the 'power basis' (the naive polynomial basis)
161 // where the integral is straight forward.
162 temp = 0;
163 for (int j=i; j<=degree; ++j){ // power basisŧ
164 temp += pow(-1.,j-i) * TMath::Binomial(degree, j) * TMath::Binomial(j,i) / (j+1);
165 }
166 temp *= ((RooAbsReal*)iter.next())->getVal(); // include coeff
167 norm += temp; // add this basis's contribution to total
168 }
169
170 norm *= xmax-xmin;
171 return norm;
172}
int Int_t
Definition: RtypesCore.h:41
double Double_t
Definition: RtypesCore.h:55
#define ClassImp(name)
Definition: Rtypes.h:365
#define R__ASSERT(e)
Definition: TError.h:96
char name[80]
Definition: TGX11.cxx:109
float xmin
Definition: THbookFile.cxx:93
float xmax
Definition: THbookFile.cxx:93
double pow(double, double)
RooAbsArg is the common abstract base class for objects that represent a value (of arbitrary type) an...
Definition: RooAbsArg.h:70
RooFIter fwdIterator() const R__SUGGEST_ALTERNATIVE("begin()
One-time forward iterator.
Int_t getSize() const
TIterator * createIterator(Bool_t dir=kIterForward) const R__SUGGEST_ALTERNATIVE("begin()
TIterator-style iteration over contained elements.
RooAbsReal is the common abstract base class for objects that represent a real value and implements f...
Definition: RooAbsReal.h:53
Bool_t matchArgs(const RooArgSet &allDeps, RooArgSet &numDeps, const RooArgProxy &a) const
Utility function for use in getAnalyticalIntegral().
Double_t getVal(const RooArgSet *normalisationSet=nullptr) const
Evaluate object.
Definition: RooAbsReal.h:81
RooArgList is a container object that can hold multiple RooAbsArg objects.
Definition: RooArgList.h:21
RooArgSet is a container object that can hold multiple RooAbsArg objects.
Definition: RooArgSet.h:28
Bernstein basis polynomials are positive-definite in the range [0,1].
Definition: RooBernstein.h:23
Double_t analyticalIntegral(Int_t code, const char *rangeName=0) const
Implements the actual analytical integral(s) advertised by getAnalyticalIntegral.
Int_t getAnalyticalIntegral(RooArgSet &allVars, RooArgSet &analVars, const char *rangeName=0) const
No analytical calculation available (yet) of integrals over subranges.
RooRealProxy _x
Definition: RooBernstein.h:39
Double_t evaluate() const
Evaluate this PDF / function / constant. Needs to be overridden by all derived classes.
RooListProxy _coefList
Definition: RooBernstein.h:40
A one-time forward iterator working on RooLinkedList or RooAbsCollection.
RooAbsArg * next()
Return next element or nullptr if at end.
virtual Bool_t add(const RooAbsArg &var, Bool_t silent=kFALSE)
Reimplementation of standard RooArgList::add()
Double_t min(const char *rname=0) const
Definition: RooRealProxy.h:56
Double_t max(const char *rname=0) const
Definition: RooRealProxy.h:57
Iterator abstract base class.
Definition: TIterator.h:30
virtual TObject * Next()=0
virtual const char * GetName() const
Returns name of object.
Definition: TNamed.h:47
Double_t x[n]
Definition: legend1.C:17
static constexpr double s
static constexpr double degree
Double_t Binomial(Int_t n, Int_t k)
Calculate the binomial coefficient n over k.
Definition: TMath.cxx:2090
Double_t SignalingNaN()
Returns a signaling NaN as defined by IEEE 754](http://en.wikipedia.org/wiki/NaN#Signaling_NaN)
Definition: TMath.h:896