//#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 "TPaveLabel.h"
#include "RooBinning.h"
#include "RooDataHist.h"
#include "RooJohnson.h"
#include "RooHypatia2.h"
#include "/home/chanchal/Documents/chanchal phd/analysis_new/bkg_study_mc15rib/WS_sample/myRooJohnsonSU.cpp"
#include "/home/chanchal/Documents/chanchal phd/analysis_new/bkg_study_mc15rib/WS_sample/myRooPolBG.cpp"

using namespace RooFit ;
using namespace std;

void fit_d0_RS_binned_expo()
{

  

  TFile *f = new TFile("RS_fit_signal_genMC.root");


  TH1 *h1 = (TH1*)f-> Get("Dz_M");

  RooRealVar d0("d0","", 1.8, 1.95,"");
  RooDataHist data("data","data", RooArgList(d0),h1);
  

   //signal best check ---------------
  RooRealVar f_J("f_{J}", "", 0.695409, 0, 1);                                                                                               
  RooRealVar mu_J("#mu_{J}","mean_johnson", 1.86580, 1.8, 1.95);                                                      
  RooRealVar sigma_J ("#sigma_{J}", "sigma_johnson", 0.0166931, 0, 0.1);                                                                      
  // RooRealVar gamma_J("#gamma_{J}", "gamma", 0.69361);                                                                               
  // RooRealVar delta_J("#delta_{J}", "delta", 1.24021); 

  RooRealVar gamma_J("#gamma_{J}", "gamma", 0.340421);                                                                               
  RooRealVar delta_J("#delta_{J}", "delta", 1.18599); 
  myRooJohnsonSU johnson("johnson", "Johnson PDF", d0, mu_J, sigma_J, delta_J, gamma_J );
  
  RooRealVar f_G1("f_{G}", "", 0.466759, 0, 1);                                                                                              
  RooRealVar mu_G1("#mu_{G1}", "mean_gauss1", 1.86457, 1.8, 1.95);                                                     
  // RooRealVar mu_G2("#mu_{G2}", "mean_gauss2", 1.86457, 1.8, 1.95);                                                         
  RooRealVar sigma_G1("#sigma_{G1}", "sigma_gauss1", 0.0100087, 0, 0.1);                                                                      
  RooRealVar sigma_G2("#sigma_{G2}", "sigma_gauss2", 0.00513264, 0, 0.01); 

  RooGaussian gauss1("gauss1", "1st gaussian PDF", d0, mu_G1, sigma_G1);                                                           
  RooGaussian gauss2("gauss2", "2nd gaussian PDF", d0, mu_G1, sigma_G2);                                                           
  RooAddPdf doublegaussian("doublegaussian", "double_gaussian", gauss1, gauss2, f_G1);                                                 
  RooAddPdf sig("sig", "johnson+double_gaussian", johnson, doublegaussian, f_J); 


  // // // background: exponential
  RooRealVar alpha1("a_{1}", "A2", -3.2452, -10.,10.);
  RooExponential bkg("exp1", "1st exponential", d0, alpha1);


  RooRealVar n_sig("N_{sig}", "n_{s}", 4212333, 0, 5000000);
  RooRealVar n_bkg("N_{bkg}", "n_{b}", 745536, 0, 900000);
  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* d0plot = d0.frame();
  data.plotOn(d0plot);
  //sig.plotOn(d0plot);
  //sig.paramOn(d0plot);

  pdf.plotOn(d0plot);
  pdf.paramOn(d0plot);
  d0plot->SetYTitle("Candidates per 1.5 MeV/#it{c}^{2}");
  d0plot->GetYaxis()->CenterTitle();
  d0plot->GetYaxis()->SetTitleSize(0.05);
  
  // chi^2 value
  cout << "chi^2 = " << d0plot->chiSquare() << endl;
  double nbins = 9*d0plot->chiSquare(9) / (d0plot->chiSquare(9) - d0plot->chiSquare());// 10 is the no. of fit parameters
  cout << "chi^2/ndf = " << d0plot->chiSquare()*nbins << "/" << nbins-9 << endl;
  cout << "chi^2/ndf = " << d0plot->chiSquare(9)<< endl;
  cout << "100bin_histo: chi^2/ndf = " << d0plot->chiSquare()*100 << "/" << 100-9 << endl;

  double d0plot_chi2 = d0plot->chiSquare();
  TPaveLabel *t1 = new TPaveLabel(0.7,0.6,0.9,0.66, Form("#chi^{2}/ndf = %f", d0plot_chi2),"brNDC");
  t1->SetBorderSize(1);
  t1->SetFillStyle(0);
  //t1->SetTextFont(40);
  t1->SetTextSize(0.70);
  t1->SetLineColor(0);

  // pull distribution
  RooHist* hpull = d0plot->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 = d0.frame(Title(" "));
  pullplot->addPlotable(hpull,"B");
  // pullplot->SetYTitle("Pull");
  // pullplot->SetXTitle("#it{m(K^{+}#pi^{-}#pi^{0})} [GeV/#it{c}^{2}]");
  pullplot->SetXTitle("#it{m(K^{-}#pi^{+}#pi^{0})} [GeV/#it{c}^{2}]");

  pullplot->SetMinimum(-4.);
  pullplot->SetMaximum(4.);
  pullplot->GetXaxis()->SetLabelSize(0.1);
  pullplot->GetXaxis()->SetTitleSize(0.1);
  pullplot->GetYaxis()->SetLabelSize(0.07);
  pullplot->GetXaxis()->CenterTitle();
  // pullplot->GetXaxis()->SetTitleSize(0.05);
  
  // chi^2 value
  
  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(d0plot,Components(RooArgSet(sig)),LineColor(kRed),LineStyle(kDashed));
  // pdf.plotOn(d0plot,Components(RooArgSet(doublegaussian)),LineColor(kBlack),LineStyle(kDashed));
  // pdf.plotOn(d0plot,Components(RooArgSet(johnson)),LineColor(kGreen),LineStyle(kDashed));
  pdf.plotOn(d0plot,Components(RooArgSet(bkg)),LineColor(kBlue),LineStyle(kDashed));



  d0plot->Draw();
  t1->Draw();

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

   
  }