// Bernstein basis polynomials are positive-definite in the range [0,1].
// In this implementation, we extend [0,1] to be the range of the parameter.
// There are n+1 Bernstein basis polynomials of degree n.
// Thus, by providing N coefficients that are positive-definite, there
// is a natural way to have well bahaved polynomail PDFs.
// For any n, the n+1 basis polynomials 'form a partition of unity', eg.
// they sum to one for all values of x. See
// http://www.idav.ucdavis.edu/education/CAGDNotes/Bernstein-Polynomials.pdf
// END_HTML
#include "RooFit.h"
#include "Riostream.h"
#include "Riostream.h"
#include <math.h>
#include "TMath.h"
#include "RooBernstein.h"
#include "RooAbsReal.h"
#include "RooRealVar.h"
#include "RooArgList.h"
ClassImp(RooBernstein)
;
RooBernstein::RooBernstein()
{
}
RooBernstein::RooBernstein(const char* name, const char* title,
RooAbsReal& x, const RooArgList& coefList):
RooAbsPdf(name, title),
_x("x", "Dependent", this, x),
_coefList("coefficients","List of coefficients",this)
{
TIterator* coefIter = coefList.createIterator() ;
RooAbsArg* coef ;
while((coef = (RooAbsArg*)coefIter->Next())) {
if (!dynamic_cast<RooAbsReal*>(coef)) {
cout << "RooBernstein::ctor(" << GetName() << ") ERROR: coefficient " << coef->GetName()
<< " is not of type RooAbsReal" << endl ;
assert(0) ;
}
_coefList.add(*coef) ;
}
delete coefIter ;
}
RooBernstein::RooBernstein(const RooBernstein& other, const char* name) :
RooAbsPdf(other, name),
_x("x", this, other._x),
_coefList("coefList",this,other._coefList)
{
}
Double_t RooBernstein::evaluate() const
{
Double_t xmin = _x.min(); Double_t xmax = _x.max();
Int_t degree= _coefList.getSize()-1;
Double_t temp=0, tempx=0;
RooFIter iter = _coefList.fwdIterator() ;
for (int i=0; i<=degree; ++i){
tempx = (_x-xmin)/(xmax-xmin);
temp += ((RooAbsReal*)iter.next())->getVal() *
TMath::Binomial(degree, i) * pow(tempx,i) * pow(1-tempx,degree-i);
}
return temp;
}
Int_t RooBernstein::getAnalyticalIntegral(RooArgSet& allVars, RooArgSet& analVars, const char* rangeName) const
{
if (rangeName && strlen(rangeName)) {
return 0 ;
}
if (matchArgs(allVars, analVars, _x)) return 1;
return 0;
}
Double_t RooBernstein::analyticalIntegral(Int_t code, const char* rangeName) const
{
assert(code==1) ;
Double_t xmin = _x.min(rangeName); Double_t xmax = _x.max(rangeName);
Int_t degree= _coefList.getSize()-1;
Double_t norm(0) ;
RooFIter iter = _coefList.fwdIterator() ;
Double_t temp=0;
for (int i=0; i<=degree; ++i){
temp = 0;
for (int j=i; j<=degree; ++j){
temp += pow(-1.,j-i) * TMath::Binomial(degree, j) * TMath::Binomial(j,i) / (j+1);
}
temp *= ((RooAbsReal*)iter.next())->getVal();
norm += temp;
}
norm *= xmax-xmin;
return norm;
}