%%cpp -d

#include "TFile.h"
#include "TROOT.h"
#include "TSystem.h"
#include "RooWorkspace.h"
#include "RooAbsData.h"
#include "RooRealVar.h"

#include "RooStats/ModelConfig.h"
#include "RooStats/ProfileLikelihoodCalculator.h"
#include "RooStats/LikelihoodInterval.h"
#include "RooStats/LikelihoodIntervalPlot.h"

using namespace RooFit;
using namespace RooStats;

struct ProfileLikelihoodOptions {

   double confLevel = 0.95;
   int nScanPoints = 50;
   bool plotAsTF1 = false;
   double poiMinPlot = 1;
   double poiMaxPlot = 0;
   bool doHypoTest = false;
   double nullValue = 0;
};

ProfileLikelihoodOptions optPL;

const char *infile = "";
const char *workspaceName = "combined";
const char *modelConfigName = "ModelConfig";
const char *dataName = "obsData";

double confLevel = optPL.confLevel;
double nScanPoints = optPL.nScanPoints;
bool plotAsTF1 = optPL.plotAsTF1;
double poiXMin = optPL.poiMinPlot;
double poiXMax = optPL.poiMaxPlot;
bool doHypoTest = optPL.doHypoTest;
double nullParamValue = optPL.nullValue;

const char *filename = "";
if (!strcmp(infile, "")) {
   filename = "results/example_combined_GaussExample_model.root";
   bool fileExist = !gSystem->AccessPathName(filename); // note opposite return code
   // if file does not exists generate with histfactory
   if (!fileExist) {
      // Normally this would be run on the command line
      cout << "will run standard hist2workspace example" << endl;
      gROOT->ProcessLine(".! prepareHistFactory .");
      gROOT->ProcessLine(".! hist2workspace config/example.xml");
      cout << "\n\n---------------------" << endl;
      cout << "Done creating example input" << endl;
      cout << "---------------------\n\n" << endl;
   }

} else
   filename = infile;

TFile *file = TFile::Open(filename);

RooWorkspace *w = (RooWorkspace *)file->Get(workspaceName);

ModelConfig *mc = (ModelConfig *)w->obj(modelConfigName);

RooAbsData *data = w->data(dataName);

ProfileLikelihoodCalculator pl(*data, *mc);
pl.SetConfidenceLevel(confLevel); // 95% interval
LikelihoodInterval *interval = pl.GetInterval();

RooRealVar *firstPOI = (RooRealVar *)mc->GetParametersOfInterest()->first();
cout << "\n>>>> RESULT : " << confLevel * 100 << "% interval on " << firstPOI->GetName() << " is : ["
     << interval->LowerLimit(*firstPOI) << ", " << interval->UpperLimit(*firstPOI) << "]\n " << endl;

cout << "making a plot of the profile likelihood function ....(if it is taking a lot of time use less points or the "
        "TF1 drawing option)\n";
LikelihoodIntervalPlot plot(interval);
plot.SetNPoints(nScanPoints); // do not use too many points, it could become very slow for some models
if (poiXMin < poiXMax)
   plot.SetRange(poiXMin, poiXMax);
TString opt;
if (plotAsTF1)
   opt += TString("tf1");
plot.Draw(opt); // use option TF1 if too slow (plot.Draw("tf1")

// if requested perform also an hypothesis test for the significance
if (doHypoTest) {
   RooArgSet nullparams("nullparams");
   nullparams.addClone(*firstPOI);
   nullparams.setRealValue(firstPOI->GetName(), nullParamValue);
   pl.SetNullParameters(nullparams);
   std::cout << "Perform Test of Hypothesis : null Hypothesis is " << firstPOI->GetName() << nullParamValue
             << std::endl;
   auto result = pl.GetHypoTest();
   std::cout << "\n>>>> Hypotheis Test Result \n";
   result->Print();
}

gROOT->GetListOfCanvases()->Draw()