#include "RooRealVar.h"
#include "RooDataSet.h"
#include "RooGaussian.h"
//#ifndef __CINT__
#include "RooGlobalFunc.h"
//#endif
#include "RooRealVar.h"
#include "RooArgList.h"
#include "RooFormulaVar.h"
#include "RooDataSet.h"
#include "RooGaussian.h"
#include "RooPolyVar.h"
#include "RooExponential.h"
#include "RooAbsData.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 multibody_2D_model_data()
{
  /*---------------------GET THE DATA & DEFINE X,Y DATASETS ----------------------------------*/
  
  //The root file 
  TFile *f = new TFile("kpipi0_semileptonic.root");
  
  //X dataset
  TTree *tree = (TTree*)f->Get("d0pi_distributions");
  RooRealVar x("d0pi_kpipi0_semileptonic","", 2.005, 2.03);//d0pi
  //RooDataSet datax("datax","data", tree, RooArgSet(x));

  //Y dataset
  TTree *tree1 = (TTree*)f->Get("d0_distributions");
  RooRealVar y("d0_kpipi0_semileptonic","", 1.75, 2.05);//d0
  //RooDataSet datay("datay","data", tree1, RooArgSet(y));

  //XY dataset
  // RooDataSet dataxy("dataxy","dataxy",RooArgSet(x,y));

  //TTree *tree = (TTree*)f1 ->Get("myTree");
  Double_t xval,yval;
  tree->SetBranchAddress("d0pi_kpipi0_semileptonic",&xval);
  tree1->SetBranchAddress("d0_kpipi0_semileptonic",&yval);
  
  RooDataSet* datax=new RooDataSet("datax","",RooArgSet(x));
  RooDataSet* datay=new RooDataSet("datay","",RooArgSet(y));
  
  Int_t nevt=tree->GetEntries();
  
  for(int i=0;i<nevt;i++)
    {
      tree->GetEntry(i);
      x.setVal(xval);
      datax->add(RooArgSet(x));
    }
  Int_t nevt1=tree1->GetEntries();
  for(int j=0;j<nevt1;j++)
    {
      tree1->GetEntry(j);
      y.setVal(yval);
      datay->add(RooArgSet(y));
    }
  datay->merge(datax);

  /*----------------------------------- MODEL THE DATA-----------------------------*/

  //Create function f(y) = a0 + a1*y (this is to consider the variation of the delta as a function of m(KK))
  RooRealVar a0("a0", "a0", -6.55247, -10, 10);
  RooRealVar a1("a1", "a1", 3.04389, 1, 10);
  RooPolyVar fy("fy", "fy", y, RooArgSet(a0, a1));// y = m(KK)
  
  //Create the johnson for d0pi i.e. x 
  //RooRealVar mu_J("#mu_{J}","mean_johnson",2.01,2.005,2.03);
  RooRealVar mu_J("#mu_{J}","mean_johnson",2.008427);
  RooRealVar sigma_J ("#sigma_{J}", "sigma_johnson", 0.0024,0,1);
  RooRealVar gamma_J("#gamma_{J}", "gamma", 1);
  RooRealVar delta_J("#delta_{J}", "delta", -10, -20., 20.);
  myRooJohnsonSU johnsonx("johnson", "Johnson PDF", x, mu_J, sigma_J, fy, gamma_J);//fy in place of delta

  RooRealVar mpi("mpi", "nominal pion mass", 0.13957039);
  RooRealVar alpha3("#alpha", "alpha", 0.);
  RooRealVar Beta3("#beta", "beta", 0.);
  myRooPolBG bkg1("bkg","background component", x, mpi, alpha3, Beta3);
  //////////////////////////////////////////////////////////////////////////////////
  
  //Create the double exponential for d0 i.e. y
  //RooRealVar f1("f_{e}", "",0.2, 0, 1);
  RooRealVar f1("f_{e}", "",0);
  //1st exponential 
  RooRealVar alpha1("a_{1}", "A2", -9.275,-10.,0.);
  RooExponential exp1("exp1", "1st exponential", y, alpha1);
  //2st exponential 
  RooRealVar alpha2("a_{2}", "A2",-15,-20.,0. );
  RooExponential exp2("exp2", "2nd exponential", y, alpha2);
  //Double exponential
  RooAddPdf expy("expy","exp1 + exp2",exp1,exp2,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", x, mpi, alpha, Beta);
  /////////////////////////////////////////////////////////////////////////
  
  RooProdPdf kpipi0_sig("kpipi0_sig_model", "kpipi0_model", RooArgSet(johnsonx,expy));
  RooProdPdf kpipi0_bkg("kpipi0_bkg_model", "kpipi0_model", RooArgSet(bkg1,bkg));

  RooRealVar n_sig("N_{sig}", "n_{s}", 1000, 0, 20000);
  RooRealVar n_bkg("N_{bkg}", "n_{b}", 0);
  RooAddPdf total("total","",RooArgList(kpipi0_sig,kpipi0_bkg),RooArgList(n_sig,n_bkg));//final signal + bkg
    
  //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));
  
  //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 THE DATA-------------------------------------*/
  //RooFitResult *fitresultx = johnsonx.fitTo(datax, Save(true), Strategy(2), Extended(true));
  //RooFitResult *fitresulty = expy.fitTo(datay, Save(true), Strategy(2), Extended(true));
  //RooFitResult *fitresultxy = total.fitTo(datay, Save(true), Strategy(2), Extended(true));    

  total.fitTo(*datax);
  
  //The Product
  RooProdPdf model("model", "johnsonx(x|y)*expy(y)", expy, Conditional(johnsonx, x));

  /*-----------------------------------PLOT THE RESULTS-------------------------------------*/   
  //Plot x distribution of data and projection of model on x
  RooPlot *xframe = x.frame();
  //datax.plotOn(xframe);
  model.plotOn(xframe);
  
  //Plot y distribution of data and projection of model on y 
  RooPlot *yframe = y.frame();
  //datay.plotOn(yframe);
  model.plotOn(yframe);

  //Plot xy distribution of data and projection of model on y 
  RooPlot *xyframe = x.frame();
  datay->plotOn(xyframe);
  total.plotOn(xyframe);
  total.paramOn(xyframe);
  //total.plotOn(xyframe,Components(johnsonx),LineStyle(kDotted),LineColor(kRed));
  //total.plotOn(xyframe,Components(expy),LineStyle(kDashed),LineColor(kCyan));
  
    
  
  // Make two-dimensional plot in x vs y
  TH1 *hh_model = model.createHistogram("hh_model", x, Binning(100), YVar(y, Binning(100)));
  hh_model->SetLineColor(kBlue);
  
  // Make canvas and draw RooPlots
  TCanvas *c = new TCanvas("rf305_condcorrprod", "rf05_condcorrprod", 1600, 400);
  c->Divide(4);
  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();
  c->cd(3);
  gPad->SetLeftMargin(0.20);
  hh_model->GetZaxis()->SetTitleOffset(2.5);
  hh_model->Draw("surf");

  c->cd(4);
  xyframe->GetYaxis()->SetTitleOffset(1.6);
  xyframe->Draw();
 
  
}