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).
Processing /mnt/vdb/lsf/workspace/root-makedoc-v608/rootspi/rdoc/src/v6-08-00-patches/tutorials/roostats/StandardProfileLikelihoodDemo.C...
[1mRooFit v3.60 -- Developed by Wouter Verkerke and David Kirkby[0m
Copyright (C) 2000-2013 NIKHEF, University of California & Stanford University
All rights reserved, please read http://roofit.sourceforge.net/license.txt
[#1] INFO:Minization -- p.d.f. provides expected number of events, including extended term in likelihood.
[#1] INFO:Minization -- Including the following contraint terms in minimization: (alpha_syst2Constraint,alpha_syst3Constraint,gamma_stat_channel1_bin_0_constraint,gamma_stat_channel1_bin_1_constraint)
[#0] PROGRESS:Minization -- ProfileLikelihoodCalcultor::DoGLobalFit - find MLE
[#0] PROGRESS:Minization -- ProfileLikelihoodCalcultor::DoMinimizeNLL - using Minuit / Migrad with strategy 1
[#1] INFO:Fitting -- RooAddition::defaultErrorLevel(nll_simPdf_obsData_with_constr) Summation contains a RooNLLVar, using its error level
[#1] INFO:Minization -- RooMinimizer::optimizeConst: activating const optimization
RooAbsTestStatistic::initSimMode: creating slave calculator #0 for state channel1 (2 dataset entries)
[#1] INFO:NumericIntegration -- RooRealIntegral::init(channel1_model_Int[obs_x_channel1]) using numeric integrator RooBinIntegrator to calculate Int(obs_x_channel1)
[#1] INFO:Fitting -- RooAbsTestStatistic::initSimMode: created 1 slave calculators.
[#1] INFO:NumericIntegration -- RooRealIntegral::init(channel1_model_Int[obs_x_channel1]) using numeric integrator RooBinIntegrator to calculate Int(obs_x_channel1)
[#1] INFO:Minization -- The following expressions have been identified as constant and will be precalculated and cached: (signal_channel1_nominal,background1_channel1_nominal,background2_channel1_nominal)
[#1] INFO:Minization -- The following expressions will be evaluated in cache-and-track mode: (mc_stat_channel1)
[#1] INFO:NumericIntegration -- RooRealIntegral::init(gamma_stat_channel1_bin_0_constraint_Int[gamma_stat_channel1_bin_0]) using numeric integrator RooIntegrator1D to calculate Int(gamma_stat_channel1_bin_0)
[#1] INFO:NumericIntegration -- RooRealIntegral::init(gamma_stat_channel1_bin_1_constraint_Int[gamma_stat_channel1_bin_1]) using numeric integrator RooIntegrator1D to calculate Int(gamma_stat_channel1_bin_1)
[#1] INFO:Minization --
RooFitResult: minimized FCN value: -1044.38, estimated distance to minimum: 2.40386e-08
covariance matrix quality: Full, accurate covariance matrix
Status : MINIMIZE=0
Floating Parameter FinalValue +/- Error
-------------------- --------------------------
SigXsecOverSM 1.1153e+00 +/- 5.86e-01
alpha_syst2 -8.9017e-03 +/- 9.82e-01
alpha_syst3 1.7918e-02 +/- 9.48e-01
gamma_stat_channel1_bin_0 9.9955e-01 +/- 4.93e-02
gamma_stat_channel1_bin_1 1.0036e+00 +/- 8.01e-02
>>>> RESULT : 95% interval on SigXsecOverSM is : [Warning: lower value for SigXsecOverSM is at limit 0
0, 2.32744]
making a plot of the profile likelihood function ....(if it is taking a lot of time use less points or the TF1 drawing option)
[#1] INFO:Minization -- RooProfileLL::evaluate(nll_simPdf_obsData_with_constr_Profile[SigXsecOverSM]) Creating instance of MINUIT
[#1] INFO:Fitting -- RooAddition::defaultErrorLevel(nll_simPdf_obsData_with_constr) Summation contains a RooNLLVar, using its error level
[#1] INFO:Minization -- RooProfileLL::evaluate(nll_simPdf_obsData_with_constr_Profile[SigXsecOverSM]) determining minimum likelihood for current configurations w.r.t all observable
[#1] INFO:Minization -- RooProfileLL::evaluate(nll_simPdf_obsData_with_constr_Profile[SigXsecOverSM]) minimum found at (SigXsecOverSM=1.11511)
.
[#1] INFO:Minization -- RooProfileLL::evaluate(nll_simPdf_obsData_with_constr_Profile[SigXsecOverSM]) Creating instance of MINUIT
[#1] INFO:Fitting -- RooAddition::defaultErrorLevel(nll_simPdf_obsData_with_constr) Summation contains a RooNLLVar, using its error level
[#1] INFO:Minization -- RooProfileLL::evaluate(nll_simPdf_obsData_with_constr_Profile[SigXsecOverSM]) determining minimum likelihood for current configurations w.r.t all observable
RooAbsTestStatistic::initSimMode: creating slave calculator #0 for state channel1 (2 dataset entries)
[#1] INFO:NumericIntegration -- RooRealIntegral::init(channel1_model_Int[obs_x_channel1]) using numeric integrator RooBinIntegrator to calculate Int(obs_x_channel1)
[#1] INFO:Fitting -- RooAbsTestStatistic::initSimMode: created 1 slave calculators.
[#1] INFO:NumericIntegration -- RooRealIntegral::init(channel1_model_Int[obs_x_channel1]) using numeric integrator RooBinIntegrator to calculate Int(obs_x_channel1)
[#1] INFO:NumericIntegration -- RooRealIntegral::init(gamma_stat_channel1_bin_0_constraint_Int[gamma_stat_channel1_bin_0]) using numeric integrator RooIntegrator1D to calculate Int(gamma_stat_channel1_bin_0)
[#1] INFO:NumericIntegration -- RooRealIntegral::init(gamma_stat_channel1_bin_1_constraint_Int[gamma_stat_channel1_bin_1]) using numeric integrator RooIntegrator1D to calculate Int(gamma_stat_channel1_bin_1)
[#0] ERROR:InputArguments -- RooArgSet::checkForDup: ERROR argument with name SigXsecOverSM is already in this set
[#1] INFO:Minization -- RooProfileLL::evaluate(nll_simPdf_obsData_with_constr_Profile[SigXsecOverSM]) minimum found at (SigXsecOverSM=1.11509)
......................................................................................................
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);
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";
}
}