#include "TFile.h"
#include "TH1F.h"
#include "TCanvas.h"
#include "TH1.h"

#include "RooFit.h"
#include "RooRealVar.h"
#include "RooGaussian.h"
#include "RooDataHist.h"
#include "RooAddPdf.h"
#include "RooChebychev.h"
#include "RooDataSet.h"
#include "RooPlot.h"

using namespace RooFit;

void testIntervalShape (Float_t massPoint=7.5, Float_t sigmaRes=0.15) {

  float lowEdge = 5.5;
  float highEdge = 8.2;
  Int_t nevents = 300000;
  Int_t nbins = 27;
  

  RooRealVar m_roo("m_roo", "m_roo", lowEdge, highEdge);  
  RooRealVar a0("a0","a0", -1.42803e-01, -1.0, 1.0) ;
  RooRealVar a1("a1","a1", 2.87161e-02, -1.0, 1.0) ;
  RooRealVar a2("a2","a2", -1.45196e-03, -1.0, 1.0) ;
  RooRealVar a3("a3","a3", -6.00993e-04, -1.0, 1.0) ;

//  RooExponential func_bk("func_bk", "Background", m_roo, a0);
  RooChebychev func_bk("func_bk", "Background", m_roo, RooArgSet(a0,a1,a2,a3));

//  TH1F* inHist = func_bk.generate(m_roo, nevents);
  RooDataSet *tmpDH = func_bk.generate(m_roo, nevents);
  TH1F* inHist = tmpDH -> createHistogram ("inHist", m_roo, Binning(nbins));
  RooDataHist dh ("dh","dh",m_roo, inHist);



  if ( (massPoint-2*sigmaRes) > lowEdge) {    // has low mass sideband
     m_roo.setRange("lowBkg",lowEdge, massPoint-2*sigmaRes);
     if (massPoint+2*sigmaRes < highEdge) {   // has high mass sideband
       m_roo.setRange("highBkg", massPoint+2*sigmaRes, highEdge);
       RooFitResult *fitRes = func_bk.fitTo(dh, Range("lowBkg,highBkg"), Save());
     }
     else RooFitResult *fitRes = func_bk.fitTo(dh, Range("lowBkg"), Save());
  }
  else  {
    m_roo.setRange("highBkg", massPoint+2*sigmaRes, highEdge);
    RooFitResult *fitRes = func_bk.fitTo(dh, Range("highBkg"), Save());
  }

  // Plot data and PDF
  RooPlot* xframe = m_roo.frame(Title("Background")) ;
  dh.plotOn(xframe, Name("data"));
  func_bk.plotOn(xframe, Name("model"));
  xframe->Draw();

  cout << "Bkg Chi2 = " << xframe-> chiSquare("model", "data") << endl;
  cout << "The End!" << endl;
  
}