//#ifndef __CINT__
#include "RooGlobalFunc.h"
//#endif
#include "RooRealVar.h"
#include "RooArgList.h"
#include "RooFormulaVar.h"
#include "RooDataSet.h"
#include "RooGaussian.h"
#include "RooChebychev.h"
#include "RooExponential.h"
#include "RooBifurGauss.h"
#include "RooAddModel.h"
#include "RooProdPdf.h"
#include "TCanvas.h"
#include "RooPlot.h"
#include "RooHist.h"
#include "RooCBShape.h"
#include "RooPolynomial.h"
#include "RooBinning.h"
#include "TH1.h"
#include "TH2.h"
#include "RooAddPdf.h"
#include "RooProdPdf.h"
#include "RooFitResult.h"
#include "RooGenericPdf.h"
#include "RooLandau.h"
#include "TChain.h"
#include "TCanvas.h"
#include "TAxis.h"
#include "RooPlot.h"
#include "RooCategory.h"
#include "RooSuperCategory.h"
#include "RooSimultaneous.h"
#include "RooNLLVar.h"
#include "TFile.h"
//#include "myRooJohnsonSU.cpp"
#include "myRooPolBG.cpp"

using namespace RooFit ;
using namespace std;

void kpipi0_semileptonic_bkg_fit_d0pi()
{
  //TFile *f = new TFile("kpipi0_semileptonic.root");
  TFile *f = new TFile("kpipi0_semileptonic_1.75to1.76.root");
  //TFile *f = new TFile("kpipi0_semileptonic_1.76to1.77.root");
  //TFile *f = new TFile("kpipi0_semileptonic_1.77to1.78.root");
  //TFile *f = new TFile("kpipi0_semileptonic_1.78to1.795.root");
  //TFile *f = new TFile("kpipi0_semileptonic_1.795to1.81.root");
  //TFile *f = new TFile("kpipi0_semileptonic_1.81to1.83.root");
  //TFile *f = new TFile("kpipi0_semileptonic_1.83to1.86.root");
  //TFile *f = new TFile("kpipi0_semileptonic_1.86to1.9.root");
  //TFile *f = new TFile("kpipi0_semileptonic_1.9to1.96.root");
  //TFile *f = new TFile("kpipi0_semileptonic_1.96to2.05.root");
    
  TTree *tree = (TTree*)f->Get("d0pi_distributions");
  cout << "Entries = " << tree->GetEntries() << endl;
  
  RooRealVar d0pi("d0pi_kpipi0_semileptonic","", 2.005, 2.03, "");
  RooDataSet data("data","data", tree, RooArgSet(d0pi));

  RooRealVar f1("f_{L}", "",0.89);
  RooRealVar delta("#delta", "shift in mean",0.01,0,5);
  RooRealVar s("#sigma", "width",0.0001,0,5);
  
  //Landau
  RooRealVar mean_L("#mu_{L}", "mean",2.01056);
  RooFormulaVar mean_L_shifted("#mu_L_shifted", "mean_{L} + delta",RooArgList(mean_L, delta));

  RooRealVar sigma_L("#sigma_{L}", "sigma", 0.00181);
  RooFormulaVar sigma_L_shifted("#sigma_{L}_shifted", "sigma_{L} + s*sigma_{L}",RooArgList(sigma_L, s));

  RooLandau landau("landau", "landau component", d0pi, mean_L_shifted, sigma_L_shifted);

  //Gaussian
  RooRealVar mean_G("#mu_{G}", "mean",2.0063);
  RooFormulaVar mean_G_shifted("#mu_G_shifted", "mean_{G} + delta",RooArgList(mean_G, delta));

  RooRealVar sigma_G("#sigma_{G}", "sigma", 0.0039);
  RooFormulaVar sigma_G_shifted("#sigma_{G}_shifted", "sigma_{G} + s*sigma_{G}",RooArgList(sigma_G, s));

  RooGaussian gauss("gauss", "gauss component", d0pi, mean_G_shifted, sigma_G_shifted);
  
  RooAddPdf sig("sig","landau + gaussian",landau,gauss,f1);
  
  //background
  RooRealVar mpi("mpi", "nominal pion mass", 0.13957039);
  RooRealVar alpha("#alpha", "alpha", 0.);
  RooRealVar beta("#beta", "beta", 0.);
  myRooPolBG bkg("bkg","background component", d0pi, mpi, alpha, beta);
  
  //total
  RooRealVar n_sig("N_{sig}", "n_{s}", 2000, 0, 20000);
  RooRealVar n_bkg("N_{bkg}", "n_{b}", 0);
  RooAddPdf pdf("pdf", "two component model",RooArgList(sig, bkg), RooArgList(n_sig, n_bkg));

  //fit
  RooFitResult *fitresult = pdf.fitTo(data, Save(true), Strategy(2), Extended(true));
      
  RooPlot* d0piplot = d0pi.frame();
  data.plotOn(d0piplot);
  //sig.plotOn(d0piplot);
  //sig.paramOn(d0piplot);
  pdf.plotOn(d0piplot);
  pdf.paramOn(d0piplot);

  cout << "chi^2 = " << d0piplot->chiSquare() << endl;

  RooHist* hpull = d0piplot->pullHist();
  hpull->SetFillStyle(1001);
  hpull->SetFillColor(1);
  for(int i=0;i<hpull->GetN();++i) hpull->SetPointError(i,0.0,0.0,0.0,0.0);
  RooPlot* pullplot = d0pi.frame(Title(""));
  pullplot->addPlotable(hpull,"B");
  // pullplot->SetYTitle("Pull");
  //pullplot->SetXTitle("#it{#Deltam} [GeV/c^{2}]");
  pullplot->SetMinimum(-4.);
  pullplot->SetMaximum(4.);
  pullplot->GetXaxis()->SetLabelSize(0.1);
  pullplot->GetXaxis()->SetTitleSize(0.1);
  pullplot->GetYaxis()->SetLabelSize(0.07);
  
    
  TCanvas *canvas = new TCanvas("canvas","canvas", 600, 600);

  Double_t xlow, ylow, xup, yup;
  canvas->GetPad(0)->GetPadPar(xlow,ylow,xup,yup);
  canvas->Divide(1,2);

  TVirtualPad *upPad = canvas->GetPad(1);
  upPad->SetPad(xlow,ylow+0.25*(yup-ylow),xup,yup);
  
  TVirtualPad *dwPad = canvas->GetPad(2);
  dwPad->SetPad(xlow,ylow,xup,ylow+0.25*(yup-ylow));

  canvas->Update();
  canvas->cd(1);
  
  pdf.plotOn(d0piplot,Components(RooArgSet(sig)),LineColor(kRed),LineStyle(kDashed));

  d0piplot->Draw();
  canvas->cd(2);
  pullplot->Draw();
  
    
  canvas->Update();
   
    
  }