With default parameters the macro will attempt to run the standard hist2workspace example and read the ROOT file that it produces.
The actual heart of the demo is only about 10 lines long.
The ProfileLikelihoodCalculator is based on Wilks's theorem and the asymptotic properties of the profile likelihood ratio (eg. that it is chi-square distributed for the true value).
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;
void StandardProfileLikelihoodDemo(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";
if (!fileExist) {
#ifdef _WIN32
cout << "HistFactory file cannot be generated on Windows - exit" << endl;
return;
#endif
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;
if(!file ){
cout <<"StandardRooStatsDemoMacro: Input file " << filename << " is not found" << endl;
return;
}
if(!w){
cout <<"workspace not found" << endl;
return;
}
if(!data || !mc){
cout << "data or ModelConfig was not found" <<endl;
return;
}
pl.SetConfidenceLevel(confLevel);
cout <<
"\n>>>> RESULT : " << confLevel*100 <<
"% interval on " <<firstPOI->
GetName()<<
" is : ["<<
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";
plot.SetNPoints(nScanPoints);
if (poiXMin < poiXMax) plot.SetRange(poiXMin, poiXMax);
TString opt;
if (plotAsTF1) opt += TString("tf1");
plot.Draw(opt);
if (doHypoTest) {
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";
}
}