#ifndef __CINT__
#include "RooGlobalFunc.h"
#endif
#include "RooRealVar.h"
#include "RooDataSet.h"
#include "RooGaussian.h"
#include "RooChebychev.h"
#include "RooAddPdf.h"
#include "RooExtendPdf.h"
#include "RooBernstein.h"
#include "TCanvas.h"
#include "TAxis.h"
#include "TH1.h"
#include "RooPlot.h"
using namespace RooFit ;


void CheckJPsi()
{

// drawing options
gStyle->SetOptStat(0);
gROOT->SetStyle("Plain");
gStyle->SetCanvasBorderMode(0);
gStyle->SetCanvasColor(kWhite);
gStyle->SetFrameBorderMode(0);
gStyle->SetFrameFillColor(kWhite);
gStyle->SetPalette(1); 



TFile *myfile = new TFile("./PDC09.root", "READ");
TTree *ATuple = (TTree *)myfile->Get("fTwoMuonNtuple");

//___________________________________________________________________________________
  // S e t u p   c o m p o n e n t   p d f s 
  // ---------------------------------------

  // Declare observable x
  RooRealVar x("x","x",2,5) ; //2.6-4.4

  RooRealVar mean("mean","mean_jpsi",3.096,3,3.2);
  RooRealVar sigma1("sigma1","sigma_jpsi",0.070,0.050,0.080);
  RooRealVar alpha("alpha","alpha",1,0.5,2);
  RooRealVar n("n","n",4,3,5);
  RooCBShape sig1("jpsi","crystal ball PDF",x,mean,sigma1,alpha,n);

  // Build Bernstein polynomial p.d.f.
  RooRealVar a0("a0","a0",0,-1.5,1);
  RooRealVar a1("a1","a1",0,-2,1) ;
  RooRealVar a2("a2","a2",0,-2,1) ;
  RooRealVar a3("a3","a3",0,-2,1) ;
  RooChebychev bkg("bkg","Background",x,RooArgSet(a0,a1,a2));

  // C o n s t r u c t   e x t e n d e d   c o m p o s i t e   m o d e l 
  // -------------------------------------------------------------------

  // Sum the composite signal and background into an extended pdf nsig*sig+nbkg*bkg
  RooRealVar nsig("nsig","number of signal events",1000,0.,15000) ;
  RooRealVar nbkg("nbkg","number of background events",1000,0.,5000);
  RooAddPdf  *model = new RooAddPdf("model","(CB+a)",RooArgList(bkg,sig1),RooArgList(nbkg,nsig)) ;


  // S a m p l e ,   f i t   a n d   p l o t   e x t e n d e d   m o d e l 
  // ---------------------------------------------------------------------
  
  
  TH1F *t = new TH1F("t","t",50,2.0,5.0);
  ATuple->Project("t","m");
  
  RooDataHist *data = new RooDataHist("data","data x",x,t);
  
  x.setRange(2.5,4.0) ;
  RooFitResult *fr = model.fitTo(*data,Save()) ;

  x.setRange(2,5.0) ;
  RooPlot *xframe_tmp = x.frame(Title("extended ML fit example")) ;
  data->plotOn(xframe_tmp,Name("data"));
  model->plotOn(xframe_tmp,Name("model"),Range(2.5,4.0)) ;
  model->plotOn(xframe_tmp ,Components(bkg),LineStyle(kDashed),Range(2.5,4.0)) ;
  model.plotOn(xframe_tmp ,Components(RooArgSet(bkg,sigma1)),LineStyle(kDotted),Range(2.5,4.0)) ;
  
  nDOF = fr->floatParsFinal().getSize();
  Double_t chi2 = xframe_tmp->chiSquare("model", "data", nDOF);
  cout << "The chi2/nDOF is " << chi2 << endl;
  
  //2.7,3.5  
  x.setRange("window",2.8,3.4) ;
  RooAbsReal* fracSigRange = sig1.createIntegral(x,x,"window") ;
  Double_t nsigWindow = nsig.getVal();
  cout << "The signal event yield is " << nsigWindow << endl;
  
  RooAbsReal* fracBkgRange = bkg.createIntegral(x,x,"window") ;
  Double_t nbkgWindow = nbkg.getVal() * fracBkgRange->getVal() ;
  cout << "The Background event yield is " << nbkgWindow << endl;
  
  xframe_tmp->Draw();

}