//#ifndef __CINT__
#include "RooGlobalFunc.h"
//#endif
#include "RooRealVar.h"
#include "RooGenericPdf.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 "RooFit.h"

using namespace RooFit ;
using namespace std;

class RooSelfNormalizedGaussian : public RooAbsPdf {
public:
  RooSelfNormalizedGaussian(const char *name, const char *title,
                            RooAbsReal& _x, RooAbsReal& _mean, RooAbsReal& _sigma) :
  RooAbsPdf(name,title),
  x("x","Observable",this,_x),
  mean("mean","Mean",this,_mean),
  sigma("sigma","Width",this,_sigma)
{}

  RooSelfNormalizedGaussian(const RooSelfNormalizedGaussian& other, const char* name=0) :
    RooAbsPdf(other,name), x("x",this,other.x), mean("mean",this,other.mean),
    sigma("sigma",this,other.sigma) {}

  virtual TObject* clone(const char* newname) const override {
    return new RooSelfNormalizedGaussian(*this,newname);
  }

  Bool_t selfNormalized() const override { return true; }

  Int_t getAnalyticalIntegral(RooArgSet& allVars, RooArgSet& analVars, const char* rangeName=0) const override {
    if (matchArgs(allVars,analVars,x)) return 1 ;
    if (matchArgs(allVars,analVars,mean)) return 2 ;
    return 0 ;
  }

  Double_t analyticalIntegral(Int_t code, const char* rangeName=0) const override {
    assert(code==1 || code==2);
    return 1.0;
  }

protected:

  Double_t evaluate() const override {
    const double arg = x - mean;
    const double sig = sigma;
    return 1./(sig * std::sqrt(TMath::TwoPi())) * std::exp(-0.5*arg*arg/(sig*sig));
  }

private:
  RooRealProxy x ;
  RooRealProxy mean ;
  RooRealProxy sigma ;
};


void unbinned_2D_fit_Kpipi0_RS_pi0M_conv_new_histo_v2()  
{
	

  TFile *f = new TFile("DtoKpipi0_RS_mc14a_100fb_total_light_2012_minus_with_cut_2D_fit_isSignal_tree_new_criteria.root");
 
  TTree *tree = (TTree*)f->Get("distributions");
  cout << "Entries = " << tree->GetEntries() << endl;
	
  
  //create both fit observables and get the data in RooFit land
  RooRealVar pi0_M_prefit("pi0_m_prefit","#it{m}(#it{#pi}^{0})", 0.1, 0.18, "GeV/#it{c}^{2}");
  RooRealVar pi0_ErrM("pi0_errm","#it{m}(#it{#pi}^{0}) uncertainty", 0, 0.03, "GeV/#it{c}^{2}");
  RooDataSet* dataxy = new RooDataSet("data","data",tree,RooArgSet(pi0_M_prefit, pi0_ErrM));



  //*  D o u b l e    G a u s s i o n     R e s o l u t i o n    f u n c t i o n
  //  ------------------------------------------------------------------------

  // create the different fit parameters to be used in the fit pdf
  
  

  RooRealVar f_G1("f_{G1}", "frac_gauss1", 0.16);  
  RooRealVar mpi0("mpi0", "mpi0", 0.1349768); // PDG pi0 mass

  RooRealVar b1("b1", "b1", 0.001, -0.01, 0.01);  
  RooRealVar b2("b2", "b2", 0.0021, -0.01, 0.01); 
  RooRealVar s1("s1", "s1", 1.025, 0, 3);     
  RooRealVar s2("s2", "s2", 1.015, 0, 3);  

 
  //combine the parameters, as they enter in different combinations in the fit pdf
  RooFormulaVar mu_G1("mu_G1", "mpi0+b1", RooArgList(mpi0,b1));
  RooFormulaVar sigma_G1("sigma_G1", "s1*pi0_errm", RooArgList(s1,pi0_ErrM));

  RooFormulaVar mu_G2("mu_G2", "mpi0+b2", RooArgList(mpi0,b2));
  RooFormulaVar sigma_G2("sigma_G2", "s2*pi0_errm", RooArgList(s2,pi0_ErrM));



  // Resolution function (Double Gaussion)
  RooSelfNormalizedGaussian gauss1("gauss1", "1st gaussian PDF", pi0_M_prefit, mu_G1, sigma_G1);
  RooSelfNormalizedGaussian gauss2("gauss2", "2nd gaussian PDF", pi0_M_prefit, mu_G2, sigma_G2);
  RooAddPdf doublegaussian("doublegaussian", "doublegaussian_sig", gauss1, gauss2, f_G1);
  //RooAddPdf doublegaussian("doublegaussian", "g1+g2", RooArgList(gauss1 ,gauss2) , RooArgList(f_G1));


  // C o n s t r u c t  p d f   f o r   e r r o r   m a s s
  // -----------------------------------------------------------------
  

  // Construct a histogram pdf to describe the shape of the ErrM distribution
  RooDataSet datamerr("datamerr","datamerr", tree, RooArgSet(pi0_ErrM));

  RooDataHist *expHistDterr = datamerr.binnedClone();
  RooHistPdf pdfErr("pdfErr", "pdfErr", pi0_ErrM, *expHistDterr);
 

  // C o n s t r u c t   c o n d i t i o n a l   p r o d u c t   d e c a y 
   // -----------------------------------------------------------------------------
 

  // Construct production of conditional decay_dm(dt|dterr) with empirical pdfErr(dterr)
  RooProdPdf pdf("pdf", "pdf", pdfErr, Conditional(doublegaussian, pi0_M_prefit));
  
  //perform the fit
  pdf.fitTo(*dataxy);


  /*-----------------------------------PLOT THE RESULTS-------------------------------------*/

  //Plot pi0_M_prefit projections
  RooPlot *xframe = pi0_M_prefit.frame();
  dataxy->plotOn(xframe);
  pdf.plotOn(xframe, ProjWData(*dataxy));
  pdf.paramOn(xframe,dataxy);
  pdf.plotOn(xframe,Components(gauss1),LineStyle(kDashed),LineColor(kRed), ProjWData(*dataxy));
  pdf.plotOn(xframe,Components(gauss2),LineStyle(kDashed),LineColor(kGreen), ProjWData(*dataxy));
  pdf.plotOn(xframe,Components(doublegaussian),LineStyle(kDashed),LineColor(kMagenta), ProjWData(*dataxy));
  
 
  //Plot pi0_ErrM projections
  RooPlot *yframe = pi0_ErrM.frame();
  dataxy->plotOn(yframe);
  pdf.plotOn(yframe);
  //pdf.plotOn(yframe,Components(johnson),LineStyle(kDashed),LineColor(kRed));
 
 
  // Make canvas and draw RooPlots
  TCanvas *c = new TCanvas("c", "c", 1200, 600);
  c->Divide(2);
  c->cd(1);
  gPad->SetLeftMargin(0.15);
  xframe->GetYaxis()->SetTitleOffset(1.6);
  xframe->Draw();
  c->cd(2);
  gPad->SetLeftMargin(0.15);
  yframe->GetYaxis()->SetTitleOffset(1.6);
  yframe->Draw();

}