#include <TBranch.h>
#include <TCanvas.h>
#include <TCut.h>
#include <TFile.h>
#include <TH1.h>
#include <TH2.h>
#include <TMath.h>
#include <TTree.h>

#ifndef __CINT__
#include <RooGlobalFunc.h>
#endif//__CINT__
#include <RooAbsPdf.h>
#include <RooDataHist.h>
#include <RooHistPdf.h>
#include <RooPlot.h>
#include <RooProdPdf.h>
#include <RooRealVar.h>

void HistPdfTest(void)
{
	using namespace RooFit;

	TFile *fin = TFile::Open("./pp510Run12/S2BGpDSTs/mc_high/merged_data.root");
	TTree *tin = (TTree*)fin->Get("stobg");

	TCut Cut_SC   = "eta < 0 && charge < 0";
	TCut Cut_DCA  = "Rpc1DcaSt1 < 100 && Rpc3DcaSt3 < 100";
	TCut Cut_FVTX = "fvtx_dr > -100 && fvtx_dphi > -100 && fvtx_dtheta > -100";
	TCut Cut_Comb = (Cut_SC + Cut_DCA + Cut_FVTX);

	TH2F *TestH = new TH2F("DG0_DDG0", "", 100, 0, 20, 100, 0, 9);

	tin->Project(TestH->GetName(), "DDG0:DG0", Cut_Comb);

	RooRealVar DDG0("DDG0", "DDG0", 0, 9.0);
	RooRealVar DG0("DG0", "DG0", 0,	20.0);

	//Compose PDFs
	//Projected into each direction
	RooDataHist *TestD1x = new RooDataHist("TestD1x", "", RooArgSet(DG0),  TestH->ProjectionX());
	RooDataHist *TestD1y = new RooDataHist("TestD1y", "", RooArgSet(DDG0), TestH->ProjectionY());
	RooHistPdf *TestP1x = new RooHistPdf("TestP1x", "", RooArgSet(DG0),  *TestD1x);
	RooHistPdf *TestP1y = new RooHistPdf("TestP1y", "", RooArgSet(DDG0), *TestD1y);
	RooAbsPdf *AP1 = new RooProdPdf("AP1", "", RooArgSet(*TestP1x, *TestP1y));
	//Whole 2D histogram
	RooDataHist *TestD2 = new RooDataHist("TestD2", "", RooArgSet(DG0, DDG0), TestH);
	RooHistPdf *TestP2 = new RooHistPdf("TestP2", "", RooArgSet(DG0, DDG0), *TestD2);
	RooAbsPdf *AP2 = new RooProdPdf("AP2", "", RooArgSet(*TestP2));

	//Draw
	RooPlot *f1_x =  DG0.frame();
	RooPlot *f1_y = DDG0.frame();
	AP1->plotOn(f1_x);
	AP1->plotOn(f1_y);
	RooPlot *f2_x =  DG0.frame();
	RooPlot *f2_y = DDG0.frame();
	AP2->plotOn(f2_x);
	AP2->plotOn(f2_y);

	TCanvas *c1 = new TCanvas("c1", "", 1280, 800);
	c1->Divide(2, 2);
	c1->cd(1);
	f1_x->Draw();
	c1->cd(2);
	f1_y->Draw();
	c1->cd(3);
	f2_x->Draw();
	c1->cd(4);
	f2_y->Draw();

	return;
}