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


void RooFit_test(){

  RooRealVar mass("mass","mass",4700,5900,"MeV/c^{2}");
  //
  // ** GAUSSIAN **
  RooRealVar  mean     ("mean"      ,"mean"  ,4600 ,"MeV/c^{2}"); 
  RooRealVar  sigma    ("sigma"     ,"sigma" ,270  ,"MeV/c^{2}"); 
  RooGaussian gauss    ("gauss"     ,"gauss" ,mass,mean,sigma);
  RooRealVar  nEv_gauss("nEv_gauss" ,"nEv_gauss",1500,0,10000);
  
  // ** CHEBYCHEV **
  RooRealVar   slope   ("slope"    ,"slope" ,-0.0001,-10.0,0.0,"MeV/c^{2}"); // -0.05
  RooChebychev cheb    ("cheb"     ,"cheb"  ,mass,slope);
  //RooExponential cheb    ("cheb"     ,"cheb"  ,mass,slope); // If you try using the RooExponential instead of RooChebychev, 
                                                              // the fit works in the two subranges. 
  RooRealVar   nEv_cheb("nEv_cheb" ,"nEv_cheb",1500,0,10000);

  // ** Set "sidebands" ranges **
  mass.setRange("SB1",4700,5100);
  mass.setRange("SB2",5500,5900); 
  

  // -------------------------------------------------------
  // ** FIRST TEST: generate and fit the CHEBYCHEV only **
  // -------------------------------------------------------
  RooExtendPdf cheb_pdf ("cheb_pdf", "cheb_pdf", cheb ,nEv_cheb);
  //
  // ** Generate **
  RooDataSet *data1 = cheb_pdf.generate(mass,3000,Extended(kTRUE));

  // ** Make plot **
  TCanvas *c1 = new TCanvas("c1","c1",800,800);
  c1->Divide(1,2);
  //
  c1->cd(1);
  // ** Fit and plot in the FULL range **
  cheb_pdf.fitTo(*data1,Extended(kTRUE));
  RooPlot *frame1a = mass.frame(Title("Invariant mass"),Bins(60));
  data1   ->plotOn(frame1a); 
  frame1a ->Draw();
  cheb_pdf.plotOn(frame1a);
  frame1a ->Draw();

  c1->cd(2);
  // ** Fit and plot in the SB1 and SB2 RANGES **
  cheb_pdf.fitTo(*data1,Extended(kTRUE),Range("SB1,SB2"));
  RooPlot *frame1b = mass.frame(Title("Invariant mass"),Bins(60));
  data1   ->plotOn(frame1b); 
  frame1b ->Draw();
  cheb_pdf.plotOn(frame1b);
  frame1b ->Draw();


  // --------------------------------------------------------------
  // ** SECOND TEST: generate and fit the CHEBYCHEV + GAUSSIAN **
  // --------------------------------------------------------------
  RooAddPdf model("model","model",
		  RooArgSet(gauss,cheb),
		  RooArgSet(nEv_gauss,nEv_cheb));
  
  // ** Generate **
  RooDataSet *data2 = model.generate(mass,3000);
  
  // ** Make plot **
  TCanvas *c2 = new TCanvas("c2","c2",800,800);
  c2->Divide(1,2);
  //
  c2->cd(1);
  // ** Fit and plot in the FULL range **
  model.fitTo(*data2,Extended(kTRUE));
  RooPlot *frame2a = mass.frame(Title("Invariant mass"),Bins(60));
  data2  ->plotOn(frame2a); 
  frame2a->Draw();
  model  .plotOn(frame2a);
  model  .plotOn(frame2a,Components(gauss),LineColor(kViolet),LineStyle(kDashed));
  model  .plotOn(frame2a,Components(cheb),LineColor(kGreen+1),LineStyle(kDashed));
  frame2a->Draw();

  c2->cd(2);
  // ** Fit and plot in the SB1 and SB2 RANGES **
  model.fitTo(*data2,Extended(kTRUE),Range("SB1,SB2"));
  RooPlot *frame2b = mass.frame(Title("Invariant mass"),Bins(60));
  data2  ->plotOn(frame2b); 
  frame2b->Draw();
  model  .plotOn(frame2b);
  model  .plotOn(frame2b,Components(gauss),LineColor(kViolet),LineStyle(kDashed));
  model  .plotOn(frame2b,Components(cheb),LineColor(kGreen+1),LineStyle(kDashed));
  frame2b->Draw();

}