#ifndef __CINT__
#include "RooGlobalFunc.h"
#endif
#include "RooWorkspace.h"
#include "RooRealVar.h"
#include "RooAbsPdf.h"
#include "RooGenericPdf.h"
#include "RooArgList.h"
#include "RooDataSet.h"
#include "RooHist.h"
#include "RooDataHist.h"
#include "RooHistPdf.h"
#include "RooPlot.h"
#include "RooGlobalFunc.h"
#include "RooGaussian.h"
#include "RooConstVar.h"
#include "RooTruthModel.h"
#include "RooFormulaVar.h"
#include "RooRealSumPdf.h"
#include "RooPolyVar.h"
#include "RooProduct.h"
#include "RooFitResult.h"
#include "TH1.h"
#include "TCanvas.h"
#include "TAxis.h"
#include "RooBernstein.h"

#include "RooStats/ModelConfig.h"


using namespace RooFit;
using namespace RooStats;

////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
//load files
TFile *data_histo = new TFile("data.root");
TFile *Z_histo = new TFile("Z_mc.root");

//define histos for this macro
TH1F * control_bq_data =(TH1F*)data_histo->Get("graph_bckg_mass_CUT_control_bq"); //DATA (control region bq)
TH1F * signal_data = (TH1F*)data_histo->Get("graph_bckg_mass_CUT_signal"); //DATA (signal region)

TH1F * Z_signal = (TH1F*)Z_histo->Get("graph_Z_signal");
////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
	RooRealVar x("x","x",40,200); //invariant mass
	
	RooDataHist SignalData("SignalData","SignalData",x,Import(*signal_data));
	RooDataHist ControlData("ControlData","ControlData",x,Import(*control_bq_data));
	RooDataHist ZSignal("ZSignal","ZSignal",x,Import(*Z_signal));
////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
	
void NEW_ask() {	

	//get control region
	RooHistPdf ctrlpdf("ctrlpdf","ctrlpdf",x,ControlData,0);

	//get signal
	RooHistPdf Zresonance("Zresonance","Zresonance",x,ZSignal,0);

	//get transfer function
	RooRealVar a0("a0","a0", -0.167266);
	RooRealVar a1("a1","a1", 0.010425);
	RooRealVar a2("a2","a2", -0.000188002);
	RooRealVar a3("a3","a3", 1.56046e-06);
	RooRealVar a4("a4","a4", -6.21518e-09);
	RooRealVar a5("a5","a5", 9.60187e-12);
     RooRealVar p1("p1","p1",2.1514e+00,2.1,2.2);
     RooRealVar p2("p2","p2",2.6933e-02,2.6e-02,2.7e-02);
	
	RooGenericPdf tf_data("tf_data","tf_data","(a0+a1*x+a2*pow(x,2)+a3*pow(x,3)+a4*pow(x,4)+a5*pow(x,5))*(1+p1*exp(-p2*x))",RooArgSet(x,a0,a1,a2,a3,a4,a5,p1,p2));

	//get qcd-bckg signal region
	RooProdPdf prod("prod","prod",tf_data,ctrlpdf);

	
	//get full model
	RooRealVar n_qcd("n_{QCD}", "n_{QCD}",6202910,6150000, 6250000);		
	RooRealVar n_Z("n_{Z}", "n_{Z}", 25855,1000, 50000);		
	
	RooArgList shapes;
	RooArgList yields;
	shapes.add(Zresonance);
	shapes.add(prod);
	yields.add(n_Z);
	yields.add(n_qcd);
    
	RooAddPdf final_pdf ("final_pdf", "final_pdf", shapes, yields);

	RooFitResult * r3 = final_pdf.fitTo(SignalData,Optimize(1),Save(kTRUE),Extended(kTRUE),Strategy(1),Hesse(kTRUE));
	r3->Print();


	//plot graphs
	RooPlot* final = x.frame(Title("signal region"));
	SignalData.plotOn(final,DataError(RooAbsData::SumW2),Name("data"));
	SignalData.statOn(final);
	final_pdf.plotOn(final,Name("bckg"));
	final_pdf.paramOn(final);
	RooHist* hpull2 = final->pullHist();
	RooPlot* frame_pul2 = x.frame(Title(" ")) ;
	frame_pul2->addPlotable(hpull2,"P");
	TCanvas* c2=new TCanvas("c2","c2",900,600);
	c2->Divide(1,2);
	c2->cd(1);
	final->Draw();
	auto leg = new TLegend(0.6,0.6,0.8,0.8);
	leg->AddEntry(final->findObject("data"),"data","lp");
	leg->AddEntry(final->findObject("bckg"),"background","lp");
	leg->Draw("same");
	c2->cd(2);
	frame_pul2->Draw();

}