{
 "cells": [
  {
   "cell_type": "markdown",
   "id": "84422f22",
   "metadata": {},
   "source": [
    "# StandardBayesianNumericalDemo\n",
    "Standard demo of the numerical Bayesian calculator\n",
    "\n",
    "This is a standard demo that can be used with any ROOT file\n",
    "prepared in the standard way.  You specify:\n",
    " - name for input ROOT file\n",
    " - name of workspace inside ROOT file that holds model and data\n",
    " - name of ModelConfig that specifies details for calculator tools\n",
    " - name of dataset\n",
    "\n",
    "With default parameters the macro will attempt to run the\n",
    "standard hist2workspace example and read the ROOT file\n",
    "that it produces.\n",
    "\n",
    "The actual heart of the demo is only about 10 lines long.\n",
    "\n",
    "The BayesianCalculator is based on Bayes's theorem\n",
    "and performs the integration using ROOT's numeric integration utilities\n",
    "\n",
    "\n",
    "\n",
    "\n",
    "**Author:** Kyle Cranmer  \n",
    "<i><small>This notebook tutorial was automatically generated with <a href= \"https://github.com/root-project/root/blob/master/documentation/doxygen/converttonotebook.py\">ROOTBOOK-izer</a> from the macro found in the ROOT repository  on Tuesday, May 19, 2026 at 08:36 PM.</small></i>"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "df1e41f1",
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "%%cpp -d\n",
    "\n",
    "#include \"TFile.h\"\n",
    "#include \"TROOT.h\"\n",
    "#include \"RooWorkspace.h\"\n",
    "#include \"RooAbsData.h\"\n",
    "#include \"RooRealVar.h\"\n",
    "\n",
    "#include \"RooUniform.h\"\n",
    "#include \"RooStats/ModelConfig.h\"\n",
    "#include \"RooStats/BayesianCalculator.h\"\n",
    "#include \"RooStats/SimpleInterval.h\"\n",
    "#include \"RooStats/RooStatsUtils.h\"\n",
    "#include \"RooPlot.h\"\n",
    "#include \"TSystem.h\"\n",
    "\n",
    "#include <cassert>\n",
    "\n",
    "using namespace RooFit;\n",
    "using namespace RooStats;\n",
    "\n",
    "struct BayesianNumericalOptions {\n",
    "\n",
    "   double confLevel = 0.95;      // interval CL\n",
    "   TString integrationType = \"\"; // integration Type (default is adaptive (numerical integration)\n",
    "   // possible values are \"TOYMC\" (toy MC integration, work when nuisances have a constraints pdf)\n",
    "   //  \"VEGAS\" , \"MISER\", or \"PLAIN\"  (these are all possible MC integration)\n",
    "   int nToys =\n",
    "      10000; // number of toys used for the MC integrations - for Vegas should be probably set to an higher value\n",
    "   bool scanPosterior =\n",
    "      false; // flag to compute interval by scanning posterior (it is more robust but maybe less precise)\n",
    "   bool plotPosterior = false; // plot posterior function after having computed the interval\n",
    "   int nScanPoints = 50; // number of points for scanning the posterior (if scanPosterior = false it is used only for\n",
    "                         // plotting). Use by default a low value to speed-up tutorial\n",
    "   int intervalType = 1; // type of interval (0 is shortest, 1 central, 2 upper limit)\n",
    "   double maxPOI = -999; // force a different value of POI for doing the scan (default is given value)\n",
    "   double nSigmaNuisance = -1; // force integration of nuisance parameters to be within nSigma of their error (do first\n",
    "                               // a model fit to find nuisance error)\n",
    "};\n",
    "\n",
    "BayesianNumericalOptions optBayes;"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "1bca93f0",
   "metadata": {},
   "source": [
    " Arguments are defined. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "018b8608",
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "const char *infile = \"\";\n",
    "const char *workspaceName = \"combined\";\n",
    "const char *modelConfigName = \"ModelConfig\";\n",
    "const char *dataName = \"obsData\";"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "337d6a96",
   "metadata": {},
   "source": [
    "option definitions"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "df337527",
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "double confLevel = optBayes.confLevel;\n",
    "TString integrationType = optBayes.integrationType;\n",
    "int nToys = optBayes.nToys;\n",
    "bool scanPosterior = optBayes.scanPosterior;\n",
    "bool plotPosterior = optBayes.plotPosterior;\n",
    "int nScanPoints = optBayes.nScanPoints;\n",
    "int intervalType = optBayes.intervalType;\n",
    "int maxPOI = optBayes.maxPOI;\n",
    "double nSigmaNuisance = optBayes.nSigmaNuisance;"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "1bb7925c",
   "metadata": {},
   "source": [
    "-------------------------------------------------------\n",
    "First part is just to access a user-defined file\n",
    "or create the standard example file if it doesn't exist"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "89e31fd7",
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "const char *filename = \"\";\n",
    "if (!strcmp(infile, \"\")) {\n",
    "   filename = \"results/example_combined_GaussExample_model.root\";\n",
    "   bool fileExist = !gSystem->AccessPathName(filename); // note opposite return code\n",
    "   // if file does not exists generate with histfactory\n",
    "   if (!fileExist) {\n",
    "      // Normally this would be run on the command line\n",
    "      cout << \"will run standard hist2workspace example\" << endl;\n",
    "      gROOT->ProcessLine(\".! prepareHistFactory .\");\n",
    "      gROOT->ProcessLine(\".! hist2workspace config/example.xml\");\n",
    "      cout << \"\\n\\n---------------------\" << endl;\n",
    "      cout << \"Done creating example input\" << endl;\n",
    "      cout << \"---------------------\\n\\n\" << endl;\n",
    "   }\n",
    "\n",
    "} else\n",
    "   filename = infile;"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "9467e2b6",
   "metadata": {},
   "source": [
    "Try to open the file"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "fd2f97e4",
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "TFile *file = TFile::Open(filename);"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "1b3448c3",
   "metadata": {},
   "source": [
    "if input file was specified but not found, quit"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "61086be2",
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "if (!file) {\n",
    "   cout << \"StandardRooStatsDemoMacro: Input file \" << filename << \" is not found\" << endl;\n",
    "   return;\n",
    "}"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "5b379f57",
   "metadata": {},
   "source": [
    "-------------------------------------------------------\n",
    "Tutorial starts here\n",
    "-------------------------------------------------------"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "7c48ff1e",
   "metadata": {},
   "source": [
    "get the workspace out of the file"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "66fbeda8",
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "RooWorkspace *w = (RooWorkspace *)file->Get(workspaceName);\n",
    "if (!w) {\n",
    "   cout << \"workspace not found\" << endl;\n",
    "   return;\n",
    "}"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "6bcc99c5",
   "metadata": {},
   "source": [
    "get the modelConfig out of the file"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "07fa9809",
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "ModelConfig *mc = (ModelConfig *)w->obj(modelConfigName);"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "71b09a44",
   "metadata": {},
   "source": [
    "get the modelConfig out of the file"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "2d24e664",
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "RooAbsData *data = w->data(dataName);"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "1b16af00",
   "metadata": {},
   "source": [
    "make sure ingredients are found"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "45a02691",
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "if (!data || !mc) {\n",
    "   w->Print();\n",
    "   cout << \"data or ModelConfig was not found\" << endl;\n",
    "   return;\n",
    "}"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "97364049",
   "metadata": {},
   "source": [
    "------------------------------------------\n",
    "create and use the BayesianCalculator\n",
    "to find and plot the 95% credible interval\n",
    "on the parameter of interest as specified\n",
    "in the model config"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "8fa73658",
   "metadata": {},
   "source": [
    "before we do that, we must specify our prior\n",
    "it belongs in the model config, but it may not have\n",
    "been specified"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "9ccafe99",
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "RooUniform prior(\"prior\", \"\", *mc->GetParametersOfInterest());\n",
    "w->import(prior);\n",
    "mc->SetPriorPdf(*w->pdf(\"prior\"));"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "c99986c1",
   "metadata": {},
   "source": [
    "do without systematics\n",
    "mc->SetNuisanceParameters(RooArgSet() );"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "6dc47623",
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "if (nSigmaNuisance > 0) {\n",
    "   RooAbsPdf *pdf = mc->GetPdf();\n",
    "   assert(pdf);\n",
    "   std::unique_ptr<RooFitResult> res{\n",
    "      pdf->fitTo(*data, Save(true), Minimizer(ROOT::Math::MinimizerOptions::DefaultMinimizerType().c_str()),\n",
    "                 Hesse(true), PrintLevel(ROOT::Math::MinimizerOptions::DefaultPrintLevel() - 1))};\n",
    "\n",
    "   res->Print();\n",
    "   RooArgList nuisPar(*mc->GetNuisanceParameters());\n",
    "   for (int i = 0; i < nuisPar.getSize(); ++i) {\n",
    "      RooRealVar *v = dynamic_cast<RooRealVar *>(&nuisPar[i]);\n",
    "      assert(v);\n",
    "      v->setMin(TMath::Max(v->getMin(), v->getVal() - nSigmaNuisance * v->getError()));\n",
    "      v->setMax(TMath::Min(v->getMax(), v->getVal() + nSigmaNuisance * v->getError()));\n",
    "      std::cout << \"setting interval for nuisance  \" << v->GetName() << \" : [ \" << v->getMin() << \" , \"\n",
    "                << v->getMax() << \" ]\" << std::endl;\n",
    "   }\n",
    "}\n",
    "\n",
    "BayesianCalculator bayesianCalc(*data, *mc);\n",
    "bayesianCalc.SetConfidenceLevel(confLevel); // 95% interval"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "00861da0",
   "metadata": {},
   "source": [
    "default of the calculator is central interval.  here use shortest , central or upper limit depending on input\n",
    "doing a shortest interval might require a longer time since it requires a scan of the posterior function"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "f9142dd8",
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "if (intervalType == 0)\n",
    "   bayesianCalc.SetShortestInterval(); // for shortest interval\n",
    "if (intervalType == 1)\n",
    "   bayesianCalc.SetLeftSideTailFraction(0.5); // for central interval\n",
    "if (intervalType == 2)\n",
    "   bayesianCalc.SetLeftSideTailFraction(0.); // for upper limit\n",
    "\n",
    "if (!integrationType.IsNull()) {\n",
    "   bayesianCalc.SetIntegrationType(integrationType); // set integrationType\n",
    "   bayesianCalc.SetNumIters(nToys); // set number of iterations (i.e. number of toys for MC integrations)\n",
    "}"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "f713c594",
   "metadata": {},
   "source": [
    "in case of toyMC make a nuisance pdf"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "03a52277",
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "if (integrationType.Contains(\"TOYMC\")) {\n",
    "   RooAbsPdf *nuisPdf = RooStats::MakeNuisancePdf(*mc, \"nuisance_pdf\");\n",
    "   cout << \"using TOYMC integration: make nuisance pdf from the model \" << std::endl;\n",
    "   nuisPdf->Print();\n",
    "   bayesianCalc.ForceNuisancePdf(*nuisPdf);\n",
    "   scanPosterior = true; // for ToyMC the posterior is scanned anyway so used given points\n",
    "}"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "54cde4a1",
   "metadata": {},
   "source": [
    "compute interval by scanning the posterior function"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "0534ddb5",
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "if (scanPosterior)\n",
    "   bayesianCalc.SetScanOfPosterior(nScanPoints);\n",
    "\n",
    "RooRealVar *poi = (RooRealVar *)mc->GetParametersOfInterest()->first();\n",
    "if (maxPOI != -999 && maxPOI > poi->getMin())\n",
    "   poi->setMax(maxPOI);\n",
    "\n",
    "SimpleInterval *interval = bayesianCalc.GetInterval();"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "983a7ac1",
   "metadata": {},
   "source": [
    "print out the interval on the first Parameter of Interest"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "716ae8df",
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "cout << \"\\n>>>> RESULT : \" << confLevel * 100 << \"% interval on \" << poi->GetName() << \" is : [\"\n",
    "     << interval->LowerLimit() << \", \" << interval->UpperLimit() << \"] \" << endl;"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "318122dc",
   "metadata": {},
   "source": [
    "end in case plotting is not requested"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "493fe830",
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "if (!plotPosterior)\n",
    "   return;"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "32ea1d31",
   "metadata": {},
   "source": [
    "make a plot\n",
    "since plotting may take a long time (it requires evaluating\n",
    "the posterior in many points) this command will speed up\n",
    "by reducing the number of points to plot - do 50"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "ea212f08",
   "metadata": {},
   "source": [
    "ignore errors of PDF if is zero"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "467bd53f",
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "RooAbsReal::setEvalErrorLoggingMode(RooAbsReal::Ignore);\n",
    "\n",
    "cout << \"\\nDrawing plot of posterior function.....\" << endl;"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "1a852ca0",
   "metadata": {},
   "source": [
    "always plot using number of scan points"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "4667c0bc",
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "bayesianCalc.SetScanOfPosterior(nScanPoints);\n",
    "\n",
    "RooPlot *plot = bayesianCalc.GetPosteriorPlot();\n",
    "plot->Draw();"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "18d81c40",
   "metadata": {},
   "source": [
    "Draw all canvases "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "96dc7e74",
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "%jsroot on\n",
    "gROOT->GetListOfCanvases()->Draw()"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "ROOT C++",
   "language": "c++",
   "name": "root"
  },
  "language_info": {
   "codemirror_mode": "text/x-c++src",
   "file_extension": ".C",
   "mimetype": " text/x-c++src",
   "name": "c++"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}