#include "RooRealVar.h"
#include "RooGaussian.h"
#include "RooUniform.h"
#include "RooDataSet.h"
#include "RooDataHist.h"
#include "RooHistFunc.h"
#include "RooRealSumPdf.h"
#include "RooParamHistFunc.h"
#include "RooHistConstraint.h"
#include "RooProdPdf.h"
#include "RooPlot.h"
#include "TCanvas.h"
#include "TPaveText.h"
#include <iostream>
#include <memory>
using namespace RooFit;
void rf709_BarlowBeeston()
{
  // First, construct a likelihood model with a Gaussian signal on top of a uniform background
	RooRealVar x("x", "x", -20, 20);
	x.setBins(25);
	RooRealVar meanG("meanG", "meanG", 1, -10, 10);
	RooRealVar sigG("sigG", "sigG", 1.5, -10, 10);
	RooGaussian g("g", "Gauss", x, meanG, sigG);
	RooUniform u("u", "Uniform", x);
  // Generate the data to be fitted
	std::unique_ptr<RooDataSet> sigData(g.generate(x, 50));
	std::unique_ptr<RooDataSet> bkgData(u.generate(x, 1000));
	RooDataSet sumData("sumData", "Gauss + Uniform", x);
	sumData.append(*sigData);
	sumData.append(*bkgData);
  // Make histogram templates for signal and background.
  // Let's take a signal distribution with low statistics and a more accurate
  // background distribution.
  // Normally, these come from Monte Carlo simulations, but we will just generate them.
	std::unique_ptr<RooDataHist> dh_sig( g.generateBinned(x, 50) );
	std::unique_ptr<RooDataHist> dh_bkg( u.generateBinned(x, 10000) );
  // ***** Case 0 - 'Rigid templates' *****
  // Construct histogram shapes for signal and background
	RooHistFunc p_h_sig("p_h_sig","p_h_sig",x,*dh_sig);
	RooHistFunc p_h_bkg("p_h_bkg","p_h_bkg",x,*dh_bkg);

  // ***** Case 1 - 'Barlow Beeston' *****
  // Construct parameterized histogram shapes for signal and background
	RooParamHistFunc p_ph_sig1("p_ph_sig","p_ph_sig",*dh_sig);
	RooParamHistFunc p_ph_bkg1("p_ph_bkg","p_ph_bkg",*dh_bkg);
	RooRealVar Asig1("Asig","Asig",1,0.01,5000);
	RooRealVar Abkg1("Abkg","Abkg",1,0.01,5000);
  // Construct the sum of these
	RooRealSumPdf model_tmp("sp_ph", "sp_ph",
		RooArgList(p_ph_sig1,p_ph_bkg1),
		RooArgList(Asig1,Abkg1),
		true);
  // Construct the subsidiary poisson measurements constraining the histogram parameters
  // These ensure that the bin contents of the histograms are only allowed to vary within
  // the statistical uncertainty of the Monte Carlo.
	RooHistConstraint hc_sig("hc_sig","hc_sig",p_ph_sig1);
	RooHistConstraint hc_bkg("hc_bkg","hc_bkg",p_ph_bkg1);
  // Construct the joint model with template PDFs and constraints
	RooProdPdf model1("model1","model1",RooArgSet(hc_sig,hc_bkg),Conditional(model_tmp,x));


	auto result1 = model1.fitTo(sumData, PrintLevel(0), Save());

	TCanvas* can = new TCanvas("can", "", 1500, 600);
	can->cd();

	auto frame = x.frame(Title(""));
	sumData.plotOn(frame);
	model1.plotOn(frame, LineColor(kBlue), VisualizeError(*result1));
  // Plot data again to show it on top of error bands:
	sumData.plotOn(frame);
	model1.plotOn(frame, LineColor(kBlue));
	model1.plotOn(frame, Components(p_ph_sig1), LineColor(kAzure));
	model1.plotOn(frame, Components(p_ph_bkg1), LineColor(kRed));
	model1.paramOn(frame, Parameters(RooArgSet(Asig1, Abkg1)));
	sigData->plotOn(frame, MarkerColor(kBlue));
	frame->Draw();

	can->SaveAs("test.png");


	x.setRange("window",-20,20);
	RooAbsReal * integral_sig1 = p_ph_sig1.createIntegral(RooArgSet(x), RooArgSet(x), "window");
	Double_t Asig1val = Asig1.getVal();
	Double_t Asig1val_err = Asig1.getError();
	Double_t integval_sig = integral_sig1->getVal();
	Double_t dintegval_sig = integral_sig1->getPropagatedError(*result1, RooArgSet(x));
	cout << "Asig1 : " << Asig1val << " +/- " << Asig1val_err << endl;
	cout << "Signal Integral in window: " << integval_sig  << " +/- " << dintegval_sig  << endl;


	x.setRange("window1",-5,0);
	RooAbsReal * integral_sig2 = p_ph_sig1.createIntegral(RooArgSet(x), RooArgSet(x), "window1");
	integval_sig = integral_sig2->getVal();
	dintegval_sig = integral_sig2->getPropagatedError(*result1, RooArgSet(x));
	cout << "Signal Integral in window1: " << integval_sig  << " +/- " << dintegval_sig  << endl;
}