ROOT logo
/*****************************************************************************
 * Project: RooFit                                                           *
 * Package: RooFitModels                                                     *
 * @(#)root/roofit:$Id: RooBernstein.cxx 28259 2009-04-16 16:21:16Z wouter $
 * Authors:                                                                  *
 *   Kyle Cranmer
 *                                                                           *
 *****************************************************************************/

//////////////////////////////////////////////////////////////////////////////
//
// BEGIN_HTML
// 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)
{
  // Constructor
  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; // n+1 polys of degree n

  Double_t temp=0, tempx=0;
  for (int i=0; i<=degree; ++i){
    tempx = (_x-xmin)/(xmax-xmin); // rescale to [0,1]
    temp += ((RooAbsReal&) _coefList[i]).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 
{
  // No analytical calculation available (yet) of integrals over subranges
  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; // n+1 polys of degree n
  Double_t norm(0) ;

  Double_t temp=0;
  for (int i=0; i<=degree; ++i){
    // for each of the i Bernstein basis polynomials
    // represent it in the 'power basis' (the naive polynomial basis)
    // where the integral is straight forward.
    temp = 0;
    for (int j=i; j<=degree; ++j){ // power basisŧ
      temp += pow(-1.,j-i) * TMath::Binomial(degree, j) * TMath::Binomial(j,i) / (j+1);
    }
    temp *= ((RooAbsReal&) _coefList[i]).getVal(); // include coeff
    norm += temp; // add this basis's contribution to total
  }

  norm *= xmax-xmin;
  return norm;
}
 RooBernstein.cxx:1
 RooBernstein.cxx:2
 RooBernstein.cxx:3
 RooBernstein.cxx:4
 RooBernstein.cxx:5
 RooBernstein.cxx:6
 RooBernstein.cxx:7
 RooBernstein.cxx:8
 RooBernstein.cxx:9
 RooBernstein.cxx:10
 RooBernstein.cxx:11
 RooBernstein.cxx:12
 RooBernstein.cxx:13
 RooBernstein.cxx:14
 RooBernstein.cxx:15
 RooBernstein.cxx:16
 RooBernstein.cxx:17
 RooBernstein.cxx:18
 RooBernstein.cxx:19
 RooBernstein.cxx:20
 RooBernstein.cxx:21
 RooBernstein.cxx:22
 RooBernstein.cxx:23
 RooBernstein.cxx:24
 RooBernstein.cxx:25
 RooBernstein.cxx:26
 RooBernstein.cxx:27
 RooBernstein.cxx:28
 RooBernstein.cxx:29
 RooBernstein.cxx:30
 RooBernstein.cxx:31
 RooBernstein.cxx:32
 RooBernstein.cxx:33
 RooBernstein.cxx:34
 RooBernstein.cxx:35
 RooBernstein.cxx:36
 RooBernstein.cxx:37
 RooBernstein.cxx:38
 RooBernstein.cxx:39
 RooBernstein.cxx:40
 RooBernstein.cxx:41
 RooBernstein.cxx:42
 RooBernstein.cxx:43
 RooBernstein.cxx:44
 RooBernstein.cxx:45
 RooBernstein.cxx:46
 RooBernstein.cxx:47
 RooBernstein.cxx:48
 RooBernstein.cxx:49
 RooBernstein.cxx:50
 RooBernstein.cxx:51
 RooBernstein.cxx:52
 RooBernstein.cxx:53
 RooBernstein.cxx:54
 RooBernstein.cxx:55
 RooBernstein.cxx:56
 RooBernstein.cxx:57
 RooBernstein.cxx:58
 RooBernstein.cxx:59
 RooBernstein.cxx:60
 RooBernstein.cxx:61
 RooBernstein.cxx:62
 RooBernstein.cxx:63
 RooBernstein.cxx:64
 RooBernstein.cxx:65
 RooBernstein.cxx:66
 RooBernstein.cxx:67
 RooBernstein.cxx:68
 RooBernstein.cxx:69
 RooBernstein.cxx:70
 RooBernstein.cxx:71
 RooBernstein.cxx:72
 RooBernstein.cxx:73
 RooBernstein.cxx:74
 RooBernstein.cxx:75
 RooBernstein.cxx:76
 RooBernstein.cxx:77
 RooBernstein.cxx:78
 RooBernstein.cxx:79
 RooBernstein.cxx:80
 RooBernstein.cxx:81
 RooBernstein.cxx:82
 RooBernstein.cxx:83
 RooBernstein.cxx:84
 RooBernstein.cxx:85
 RooBernstein.cxx:86
 RooBernstein.cxx:87
 RooBernstein.cxx:88
 RooBernstein.cxx:89
 RooBernstein.cxx:90
 RooBernstein.cxx:91
 RooBernstein.cxx:92
 RooBernstein.cxx:93
 RooBernstein.cxx:94
 RooBernstein.cxx:95
 RooBernstein.cxx:96
 RooBernstein.cxx:97
 RooBernstein.cxx:98
 RooBernstein.cxx:99
 RooBernstein.cxx:100
 RooBernstein.cxx:101
 RooBernstein.cxx:102
 RooBernstein.cxx:103
 RooBernstein.cxx:104
 RooBernstein.cxx:105
 RooBernstein.cxx:106
 RooBernstein.cxx:107
 RooBernstein.cxx:108
 RooBernstein.cxx:109
 RooBernstein.cxx:110
 RooBernstein.cxx:111
 RooBernstein.cxx:112
 RooBernstein.cxx:113
 RooBernstein.cxx:114
 RooBernstein.cxx:115
 RooBernstein.cxx:116
 RooBernstein.cxx:117
 RooBernstein.cxx:118
 RooBernstein.cxx:119
 RooBernstein.cxx:120
 RooBernstein.cxx:121
 RooBernstein.cxx:122
 RooBernstein.cxx:123
 RooBernstein.cxx:124
 RooBernstein.cxx:125
 RooBernstein.cxx:126
 RooBernstein.cxx:127
 RooBernstein.cxx:128
 RooBernstein.cxx:129
 RooBernstein.cxx:130
 RooBernstein.cxx:131