#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 "RooDataSet.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 "RooStats/HistFactory/PiecewiseInterpolation.h"
#include "RooMCStudy.h"
#include "RooStats/ModelConfig.h"
#include "RooStats/BernsteinCorrection.h"
#include "TH1F.h"
#include "TPad.h"
#include "TLatex.h"
#include "TCanvas.h"
#include "TString.h"
using namespace RooFit;
using namespace RooStats;


void limit_calc_pdf()
{


  // open input file 
  TFile *file = TFile::Open("model_config_pdf.root");

  // get the workspace out of the file
  RooWorkspace* w = (RooWorkspace*) file->Get("w");
  const char * modelSBName="ModelConfig";

   RooAbsData * data = w->data("genData"); 
//   ModelConfig* bModel = (ModelConfig*) w->obj(modelBName);
   ModelConfig* sbModel = (ModelConfig*) w->obj(modelSBName);


ModelConfig* bModel = (ModelConfig*) sbModel->Clone();
RooRealVar* poi = (RooRealVar*) sbModel->GetParametersOfInterest()->first();
poi->setVal(0) ;
bModel->SetSnapshot( *poi  );

//from https://roostatsworkbook.readthedocs.io/en/latest/inferencetools.html#upper-limits


RooStats::AsymptoticCalculator  asympCalc(*data, *bModel, *sbModel);
asympCalc.SetOneSided(true);
RooStats::HypoTestInverter inverter(asympCalc);
inverter.SetConfidenceLevel(0.95);
inverter.UseCLs(true);
inverter.SetVerbose(true);
inverter.SetFixedScan(1000,0.0,10000); // set number of points for pdf case
RooStats::HypoTestInverterResult* result =  inverter.GetInterval();

std::cout << "Expected upper limits, using the B (alternate) model : " << std::endl;
std::cout << " expected limit (median) " << result->GetExpectedUpperLimit(0) << std::endl;
std::cout << " expected limit (-1 sig) " << result->GetExpectedUpperLimit(-1) << std::endl;
std::cout << " expected limit (+1 sig) " << result->GetExpectedUpperLimit(1) << std::endl;
std::cout << " expected limit (-2 sig) " << result->GetExpectedUpperLimit(-2) << std::endl;
std::cout << " expected limit (+2 sig) " << result->GetExpectedUpperLimit(2) << std::endl;

RooStats::HypoTestInverterPlot* plot = new RooStats::HypoTestInverterPlot("HTI_Result_Plot","HypoTest Scan Result",result);
plot->Draw("EXP");  // plot also CLb and CLs+b



}