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