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 BayesianCalculator is based on Bayes's theorem and performs the integration using ROOT's numeric integration utilities
struct BayesianNumericalOptions {
double confLevel = 0.95 ;
TString integrationType = "";
int nToys = 10000;
bool scanPosterior = false;
int nScanPoints = 20;
int intervalType = 1;
double maxPOI = -999;
double nSigmaNuisance = -1;
};
BayesianNumericalOptions optBayes;
void StandardBayesianNumericalDemo(const char* infile = "",
const char* workspaceName = "combined",
const char* modelConfigName = "ModelConfig",
const char* dataName = "obsData") {
double confLevel = optBayes.confLevel;
TString integrationType = optBayes.integrationType;
int nToys = optBayes.nToys;
bool scanPosterior = optBayes.scanPosterior;
int nScanPoints = optBayes.nScanPoints;
int intervalType = optBayes.intervalType;
int maxPOI = optBayes.maxPOI;
double nSigmaNuisance = optBayes.nSigmaNuisance;
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;
}
if (nSigmaNuisance > 0) {
assert(pdf);
for (
int i = 0; i < nuisPar.
getSize(); ++i) {
assert( v);
std::cout <<
"setting interval for nuisance " << v->
GetName() <<
" : [ " << v->
getMin() <<
" , " << v->
getMax() <<
" ]" << std::endl;
}
}
bayesianCalc.SetConfidenceLevel(confLevel);
if (intervalType == 0) bayesianCalc.SetShortestInterval();
if (intervalType == 1) bayesianCalc.SetLeftSideTailFraction(0.5);
if (intervalType == 2) bayesianCalc.SetLeftSideTailFraction(0.);
if (!integrationType.IsNull() ) {
bayesianCalc.SetIntegrationType(integrationType);
bayesianCalc.SetNumIters(nToys);
}
if (integrationType.Contains("TOYMC") ) {
cout << "using TOYMC integration: make nuisance pdf from the model " << std::endl;
bayesianCalc.ForceNuisancePdf(*nuisPdf);
scanPosterior = true;
}
if (scanPosterior)
bayesianCalc.SetScanOfPosterior(nScanPoints);
if (maxPOI != -999 && maxPOI > poi->
getMin())
cout <<
"\n>>>> RESULT : " << confLevel*100 <<
"% interval on " << poi->
GetName()<<
" is : ["<<
cout << "\nDrawing plot of posterior function....." << endl;
bayesianCalc.SetScanOfPosterior(nScanPoints);
RooPlot * plot = bayesianCalc.GetPosteriorPlot();
}