//#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 "RooGamma.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 d0pi_fit()
{
  
  TFile *f = new TFile("kpipi0_semileptonic_1.75to1.78.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_{G}", "",0.3,0,1);
   
  //Johnson
  RooRealVar mu_J("#mu_{J}","mean_johnson",2.01,2.006, 2.03);                                                          
  RooRealVar sigma_J ("#sigma_{J}", "sigma_johnson",  0.001, 0, 0.1);                                                                        
  RooRealVar gamma_J("#gamma_{J}", "gamma", 1, -20., 20.);                                                                                 
  RooRealVar delta_J("#delta_{J}", "delta", 10, -20., 20.);                                                                                 
  myRooJohnsonSU johnson("johnson", "Johnson PDF", d0pi, mu_J, sigma_J, delta_J, gamma_J );

  //RooRealVar x("x","x", 1.0, 0.0, 10.0);
  //RooRealVar x0("x0","x0", 0);
  //RooRealVar a("a","a", 1.1);
  //RooRealVar b("b","b", 1);

  RooRealVar x0("x0","x0", 2.005);
  RooRealVar a("a","a", 2.006);
  RooRealVar b("b","b", 2.005);
  
  RooFormulaVar gamma("#gamma", "a + 1", RooArgList(a));
  RooFormulaVar beta("#beta", "1./b", RooArgList(b));
  RooGamma model("model", "Gamma pdf", d0pi, gamma, beta, x0);

  RooAddPdf sig("sig","gamma + johnson",model,johnson,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}", 10000, 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();
   
    
  }