ROOT   6.10/09 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
13 Bernstein basis polynomials are positive-definite in the range [0,1].
14 In this implementation, we extend [0,1] to be the range of the parameter.
15 There are n+1 Bernstein basis polynomials of degree n.
16 Thus, by providing N coefficients that are positive-definite, there
17 is a natural way to have well behaved polynomial PDFs.
18 For any n, the n+1 basis polynomials 'form a partition of unity', eg.
19  they sum to one for all values of x. See
20 http://www.idav.ucdavis.edu/education/CAGDNotes/Bernstein-Polynomials.pdf
21 **/
22
23 #include "RooFit.h"
24
25 #include "Riostream.h"
26 #include "Riostream.h"
27 #include <math.h>
28 #include "TMath.h"
29 #include "RooBernstein.h"
30 #include "RooAbsReal.h"
31 #include "RooRealVar.h"
32 #include "RooArgList.h"
33
34 #include "TError.h"
35
36 using namespace std;
37
39
40 ////////////////////////////////////////////////////////////////////////////////
41
43 {
44 }
45
46 ////////////////////////////////////////////////////////////////////////////////
47 /// Constructor
48
49 RooBernstein::RooBernstein(const char* name, const char* title,
50  RooAbsReal& x, const RooArgList& coefList):
51  RooAbsPdf(name, title),
52  _x("x", "Dependent", this, x),
53  _coefList("coefficients","List of coefficients",this)
54 {
55  TIterator* coefIter = coefList.createIterator() ;
56  RooAbsArg* coef ;
57  while((coef = (RooAbsArg*)coefIter->Next())) {
58  if (!dynamic_cast<RooAbsReal*>(coef)) {
59  cout << "RooBernstein::ctor(" << GetName() << ") ERROR: coefficient " << coef->GetName()
60  << " is not of type RooAbsReal" << endl ;
61  R__ASSERT(0) ;
62  }
64  }
65  delete coefIter ;
66 }
67
68 ////////////////////////////////////////////////////////////////////////////////
69
70 RooBernstein::RooBernstein(const RooBernstein& other, const char* name) :
71  RooAbsPdf(other, name),
72  _x("x", this, other._x),
73  _coefList("coefList",this,other._coefList)
74 {
75 }
76
77 ////////////////////////////////////////////////////////////////////////////////
78
80 {
81  Double_t xmin = _x.min();
82  Double_t x = (_x - xmin) / (_x.max() - xmin); // rescale to [0,1]
83  Int_t degree = _coefList.getSize() - 1; // n+1 polys of degree n
85
86  if(degree == 0) {
87
88  return ((RooAbsReal *)iter.next())->getVal();
89
90  } else if(degree == 1) {
91
92  Double_t a0 = ((RooAbsReal *)iter.next())->getVal(); // c0
93  Double_t a1 = ((RooAbsReal *)iter.next())->getVal() - a0; // c1 - c0
94  return a1 * x + a0;
95
96  } else if(degree == 2) {
97
98  Double_t a0 = ((RooAbsReal *)iter.next())->getVal(); // c0
99  Double_t a1 = 2 * (((RooAbsReal *)iter.next())->getVal() - a0); // 2 * (c1 - c0)
100  Double_t a2 = ((RooAbsReal *)iter.next())->getVal() - a1 - a0; // c0 - 2 * c1 + c2
101  return (a2 * x + a1) * x + a0;
102
103  } else if(degree > 2) {
104
105  Double_t t = x;
106  Double_t s = 1 - x;
107
108  Double_t result = ((RooAbsReal *)iter.next())->getVal() * s;
109  for(Int_t i = 1; i < degree; i++) {
110  result = (result + t * TMath::Binomial(degree, i) * ((RooAbsReal *)iter.next())->getVal()) * s;
111  t *= x;
112  }
113  result += t * ((RooAbsReal *)iter.next())->getVal();
114
115  return result;
116  }
117
118  // in case list of arguments passed is empty
119  return TMath::SignalingNaN();
120 }
121
122 ////////////////////////////////////////////////////////////////////////////////
123 /// No analytical calculation available (yet) of integrals over subranges
124
125 Int_t RooBernstein::getAnalyticalIntegral(RooArgSet& allVars, RooArgSet& analVars, const char* rangeName) const
126 {
127  if (rangeName && strlen(rangeName)) {
128  return 0 ;
129  }
130
131  if (matchArgs(allVars, analVars, _x)) return 1;
132  return 0;
133 }
134
135 ////////////////////////////////////////////////////////////////////////////////
136
137 Double_t RooBernstein::analyticalIntegral(Int_t code, const char* rangeName) const
138 {
139  R__ASSERT(code==1) ;
140  Double_t xmin = _x.min(rangeName); Double_t xmax = _x.max(rangeName);
141  Int_t degree= _coefList.getSize()-1; // n+1 polys of degree n
142  Double_t norm(0) ;
143
144  RooFIter iter = _coefList.fwdIterator() ;
145  Double_t temp=0;
146  for (int i=0; i<=degree; ++i){
147  // for each of the i Bernstein basis polynomials
148  // represent it in the 'power basis' (the naive polynomial basis)
149  // where the integral is straight forward.
150  temp = 0;
151  for (int j=i; j<=degree; ++j){ // power basis≈ß
152  temp += pow(-1.,j-i) * TMath::Binomial(degree, j) * TMath::Binomial(j,i) / (j+1);
153  }
154  temp *= ((RooAbsReal*)iter.next())->getVal(); // include coeff
155  norm += temp; // add this basis's contribution to total
156  }
157
158  norm *= xmax-xmin;
159  return norm;
160 }
virtual const char * GetName() const
Returns name of object.
Definition: TNamed.h:47
RooRealProxy _x
Definition: RooBernstein.h:39
TIterator * createIterator(Bool_t dir=kIterForward) const
float xmin
Definition: THbookFile.cxx:93
Bool_t matchArgs(const RooArgSet &allDeps, RooArgSet &numDeps, const RooArgProxy &a) const
Utility function for use in getAnalyticalIntegral().
Double_t getVal(const RooArgSet *set=0) const
Definition: RooAbsReal.h:64
Bernstein basis polynomials are positive-definite in the range [0,1].
Definition: RooBernstein.h:23
RooListProxy _coefList
Definition: RooBernstein.h:40
#define R__ASSERT(e)
Definition: TError.h:96
int Int_t
Definition: RtypesCore.h:41
Int_t getAnalyticalIntegral(RooArgSet &allVars, RooArgSet &analVars, const char *rangeName=0) const
No analytical calculation available (yet) of integrals over subranges.
STL namespace.
Iterator abstract base class.
Definition: TIterator.h:30
Double_t analyticalIntegral(Int_t code, const char *rangeName=0) const
Implements the actual analytical integral(s) advertised by getAnalyticalIntegral. ...
Double_t x[n]
Definition: legend1.C:17
double pow(double, double)
virtual Bool_t add(const RooAbsArg &var, Bool_t silent=kFALSE)
Int_t getSize() const
Double_t evaluate() const
Double_t Binomial(Int_t n, Int_t k)
Calculate the binomial coefficient n over k.
Definition: TMath.cxx:2078
float xmax
Definition: THbookFile.cxx:93
Double_t SignalingNaN()
Definition: TMath.h:791
RooAbsArg * next()
#define ClassImp(name)
Definition: Rtypes.h:336
Double_t min(const char *rname=0) const
Definition: RooRealProxy.h:56
double Double_t
Definition: RtypesCore.h:55
RooAbsReal is the common abstract base class for objects that represent a real value and implements f...
Definition: RooAbsReal.h:53
RooFIter fwdIterator() const
Double_t max(const char *rname=0) const
Definition: RooRealProxy.h:57
RooAbsPdf is the abstract interface for all probability density functions The class provides hybrid a...
Definition: RooAbsPdf.h:41
virtual TObject * Next()=0
double result[121]
RooAbsArg is the common abstract base class for objects that represent a value (of arbitrary type) an...
Definition: RooAbsArg.h:66
double norm(double *x, double *p)
Definition: unuranDistr.cxx:40