#include "RooStats/SimpleInterval.h"
#include "RooStats/BayesianCalculator.h"
#include "RooStats/MCMCCalculator.h"
#include "RooStats/MCMCIntervalPlot.h"
#include "RooStats/SequentialProposal.h"
#include "RooStats/ProposalHelper.h"
#include "RooStats/ProfileLikelihoodCalculator.h"
#include "RooStats/LikelihoodInterval.h"
#include "RooStats/LikelihoodIntervalPlot.h"
#include "RooStats/HypoTestResult.h"

using namespace std;
using namespace RooFit;
using namespace RooStats;



void cbLimit_Minimal(){

	RooFit::MsgLevel msglevel = RooMsgService::instance().globalKillBelow();
	RooMsgService::instance().setGlobalKillBelow(RooFit::FATAL);	
	
	TFile *f = new TFile("dimuon_13TeV_2016.root");
	
	RooWorkspace * _myWS = (RooWorkspace*) f->Get("dimuon_13TeV_2016");

	_myWS->Print();

	_myWS->var("peak")->setVal(1000);
	
	ModelConfig* mc = (ModelConfig*) _myWS->obj("mc");
	
	RooDataSet* prunedData = (RooDataSet*)_myWS->data("data");

	_myWS->var("ratio")->setRange(0, 1.08e-5);
	_myWS->var("ratio")->setVal(1e-5);
	
	//getting estimate of limit from PLC
	ProfileLikelihoodCalculator plc(*prunedData, *mc);
	plc.SetConfidenceLevel(0.95);
	
	mpPlrInt = plc.GetInterval();	
	RooRealVar * _poi = (RooRealVar *)mc->GetParametersOfInterest()->first();
	double upper_limit = mpPlrInt->UpperLimit( *_poi );

	
	cout << "Upper limit from PCL: " << upper_limit << endl;

	cout << "Now testing MCMC" << endl;
	
	
	MCMCInterval * mcInt;

	// FIXME: testing: this proposal function seems fairly robust
	SequentialProposal sp(0.5);

	MCMCCalculator mcmc( *prunedData, *mc );
	mcmc.SetConfidenceLevel(0.95);
	mcmc.SetNumIters(10000);          // Metropolis-Hastings algorithm iterations

	mcmc.SetProposalFunction(sp);

	mcmc.SetNumBurnInSteps(1000); // first N steps to be ignored as burn-in
	mcmc.SetLeftSideTailFraction(0.0);
	mcmc.SetNumBins(100);

    mcInt = mcmc.GetInterval();	
	
    RooRealVar * p_first_poi = (RooRealVar*) mc->GetParametersOfInterest()->first();
    double poi_limit = mcInt->UpperLimit(*p_first_poi);
   
    cout << "limti from MCMC: " <<  poi_limit << endl;	
	
	}