//////////////////////////////////////////////////////////////////////////////
//
// BEGIN_HTML
// p.d.f shape that can be used to model the background in
// the D*-D0 mass-difference distributions
// END_HTML
//

#include "myRooPolBG.h"

#include "RooFit.h"

using namespace std;

//ClassImp(myRooPolBG)

//_____________________________________________________________________________
myRooPolBG::myRooPolBG(const char *name, const char *title,
							  RooAbsReal& _dm, RooAbsReal& _dm0,
							  RooAbsReal& _b, RooAbsReal& _c) :
  RooAbsPdf(name,title),
  dm("dm","Dstar-D0 Mass Diff",this, _dm),
  dm0("dm0","Threshold",this, _dm0),
  B("B","Shape Parameter 2",this, _b),
  C("C","Shape Parameter 3",this, _c)
{;}


//_____________________________________________________________________________
myRooPolBG::myRooPolBG(const myRooPolBG& other, const char *name) :
	RooAbsPdf(other,name), dm("dm",this,other.dm), dm0("dm0",this,other.dm0),
	B("B",this,other.B), C("C",this,other.C)
{;}


//_____________________________________________________________________________
Double_t myRooPolBG::evaluate() const
{
  Double_t arg = (dm - dm0);
  if (arg <= 0.) return 0.;
  Double_t val = std::pow(arg,0.5) + B*std::pow(arg,1.5)+ C*std::pow(arg,2.5);
  return val;
}

//_____________________________________________________________________________
Double_t myRooPolBG::analyticalIntegral(Int_t code, const char* rangeName) const
{
	assert(code==1);
	
	Double_t x_min = dm.min(rangeName) - dm0;
	if (x_min<0.) x_min = 0.;
	Double_t x_max = dm.max(rangeName) - dm0;
	if (x_max<x_min) return 0.;
	
	Double_t val_min = std::pow(x_min,1.5)/1.5 + B*std::pow(x_min,2.5)/2.5 + C*std::pow(x_min,3.5)/3.5;
	Double_t val_max = std::pow(x_max,1.5)/1.5 + B*std::pow(x_max,2.5)/2.5 + C*std::pow(x_max,3.5)/3.5;
	
	return (val_max-val_min);
}