{ "cells": [ { "cell_type": "markdown", "id": "94c6edd3", "metadata": {}, "source": [ "# rs101_limitexample\n", "Limits: number counting experiment with uncertainty on both the background rate and signal efficiency.\n", "\n", "The usage of a Confidence Interval Calculator to set a limit on the signal is illustrated\n", "\n", "\n", "\n", "\n", "**Author:** Kyle Cranmer \n", "This notebook tutorial was automatically generated with ROOTBOOK-izer from the macro found in the ROOT repository on Tuesday, May 19, 2026 at 08:35 PM." ] }, { "cell_type": "code", "execution_count": 1, "id": "e376bc1b", "metadata": { "collapsed": false, "execution": { "iopub.execute_input": "2026-05-19T20:35:50.472019Z", "iopub.status.busy": "2026-05-19T20:35:50.471877Z", "iopub.status.idle": "2026-05-19T20:35:50.498765Z", "shell.execute_reply": "2026-05-19T20:35:50.498190Z" } }, "outputs": [], "source": [ "%%cpp -d\n", "#include \"RooProfileLL.h\"\n", "#include \"RooAbsPdf.h\"\n", "#include \"RooStats/HypoTestResult.h\"\n", "#include \"RooRealVar.h\"\n", "#include \"RooPlot.h\"\n", "#include \"RooDataSet.h\"\n", "#include \"RooTreeDataStore.h\"\n", "#include \"TTree.h\"\n", "#include \"TCanvas.h\"\n", "#include \"TLine.h\"\n", "#include \"TStopwatch.h\"\n", "\n", "#include \"RooStats/ProfileLikelihoodCalculator.h\"\n", "#include \"RooStats/MCMCCalculator.h\"\n", "#include \"RooStats/UniformProposal.h\"\n", "#include \"RooStats/FeldmanCousins.h\"\n", "#include \"RooStats/NumberCountingPdfFactory.h\"\n", "#include \"RooStats/ConfInterval.h\"\n", "#include \"RooStats/PointSetInterval.h\"\n", "#include \"RooStats/LikelihoodInterval.h\"\n", "#include \"RooStats/LikelihoodIntervalPlot.h\"\n", "#include \"RooStats/RooStatsUtils.h\"\n", "#include \"RooStats/ModelConfig.h\"\n", "#include \"RooStats/MCMCInterval.h\"\n", "#include \"RooStats/MCMCIntervalPlot.h\"\n", "#include \"RooStats/ProposalFunction.h\"\n", "#include \"RooStats/ProposalHelper.h\"\n", "#include \"RooFitResult.h\"\n", "#include \"TGraph2D.h\"\n", "\n", "#include \n", "\n", "using namespace RooFit;\n", "using namespace RooStats;" ] }, { "cell_type": "markdown", "id": "2c81b9a5", "metadata": {}, "source": [ "--------------------------------------\n", "An example of setting a limit in a number counting experiment with uncertainty on background and signal" ] }, { "cell_type": "markdown", "id": "1aa7a953", "metadata": {}, "source": [ "to time the macro" ] }, { "cell_type": "code", "execution_count": 2, "id": "98369a3e", "metadata": { "collapsed": false, "execution": { "iopub.execute_input": "2026-05-19T20:35:50.509190Z", "iopub.status.busy": "2026-05-19T20:35:50.509062Z", "iopub.status.idle": "2026-05-19T20:35:50.858981Z", "shell.execute_reply": "2026-05-19T20:35:50.858509Z" } }, "outputs": [], "source": [ "TStopwatch t;\n", "t.Start();" ] }, { "cell_type": "markdown", "id": "e1f46fc4", "metadata": {}, "source": [ "--------------------------------------\n", "The Model building stage\n", "--------------------------------------" ] }, { "cell_type": "code", "execution_count": 3, "id": "0840cf6b", "metadata": { "collapsed": false, "execution": { "iopub.execute_input": "2026-05-19T20:35:50.870593Z", "iopub.status.busy": "2026-05-19T20:35:50.870458Z", "iopub.status.idle": "2026-05-19T20:35:51.103136Z", "shell.execute_reply": "2026-05-19T20:35:51.097317Z" } }, "outputs": [], "source": [ "RooWorkspace wspace{};\n", "wspace.factory(\"Poisson::countingModel(obs[150,0,300], \"\n", " \"sum(s[50,0,120]*ratioSigEff[1.,0,3.],b[100]*ratioBkgEff[1.,0.,3.]))\"); // counting model" ] }, { "cell_type": "markdown", "id": "0aa99442", "metadata": {}, "source": [ "wspace.factory(\"Gaussian::sigConstraint(ratioSigEff,1,0.05)\"); // 5% signal efficiency uncertainty\n", "wspace.factory(\"Gaussian::bkgConstraint(ratioBkgEff,1,0.1)\"); // 10% background efficiency uncertainty" ] }, { "cell_type": "code", "execution_count": 4, "id": "b0c8e7c3", "metadata": { "collapsed": false, "execution": { "iopub.execute_input": "2026-05-19T20:35:51.115198Z", "iopub.status.busy": "2026-05-19T20:35:51.115063Z", "iopub.status.idle": "2026-05-19T20:35:51.328051Z", "shell.execute_reply": "2026-05-19T20:35:51.327436Z" } }, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "\n", "RooWorkspace() contents\n", "\n", "variables\n", "---------\n", "(b,gSigBkg,gSigEff,obs,ratioBkgEff,ratioSigEff,s)\n", "\n", "p.d.f.s\n", "-------\n", "RooGaussian::bkgConstraint[ x=gSigBkg mean=ratioBkgEff sigma=0.2 ] = 1\n", "RooPoisson::countingModel[ x=obs mean=countingModel_2 ] = 0.0325554\n", "RooProdPdf::modelWithConstraints[ countingModel * sigConstraint * bkgConstraint ] = 0.0325554\n", "RooGaussian::sigConstraint[ x=gSigEff mean=ratioSigEff sigma=0.05 ] = 1\n", "\n", "functions\n", "--------\n", "RooAddition::countingModel_2[ countingModel_2_[s_x_ratioSigEff] + countingModel_2_[b_x_ratioBkgEff] ] = 150\n", "RooProduct::countingModel_2_[b_x_ratioBkgEff][ b * ratioBkgEff ] = 100\n", "RooProduct::countingModel_2_[s_x_ratioSigEff][ s * ratioSigEff ] = 50\n", "\n" ] } ], "source": [ "wspace.factory(\"Gaussian::sigConstraint(gSigEff[1,0,3],ratioSigEff,0.05)\"); // 5% signal efficiency uncertainty\n", "wspace.factory(\"Gaussian::bkgConstraint(gSigBkg[1,0,3],ratioBkgEff,0.2)\"); // 10% background efficiency uncertainty\n", "wspace.factory(\"PROD::modelWithConstraints(countingModel,sigConstraint,bkgConstraint)\"); // product of terms\n", "wspace.Print();\n", "\n", "RooAbsPdf *modelWithConstraints = wspace.pdf(\"modelWithConstraints\"); // get the model\n", "RooRealVar *obs = wspace.var(\"obs\"); // get the observable\n", "RooRealVar *s = wspace.var(\"s\"); // get the signal we care about\n", "RooRealVar *b =\n", " wspace.var(\"b\"); // get the background and set it to a constant. Uncertainty included in ratioBkgEff\n", "b->setConstant();\n", "\n", "RooRealVar *ratioSigEff = wspace.var(\"ratioSigEff\"); // get uncertain parameter to constrain\n", "RooRealVar *ratioBkgEff = wspace.var(\"ratioBkgEff\"); // get uncertain parameter to constrain\n", "RooArgSet constrainedParams(*ratioSigEff,\n", " *ratioBkgEff); // need to constrain these in the fit (should change default behavior)\n", "\n", "RooRealVar *gSigEff = wspace.var(\"gSigEff\"); // global observables for signal efficiency\n", "RooRealVar *gSigBkg = wspace.var(\"gSigBkg\"); // global obs for background efficiency\n", "gSigEff->setConstant();\n", "gSigBkg->setConstant();" ] }, { "cell_type": "markdown", "id": "f83d85a6", "metadata": {}, "source": [ "Create an example dataset with 160 observed events" ] }, { "cell_type": "code", "execution_count": 5, "id": "10ed5f55", "metadata": { "collapsed": false, "execution": { "iopub.execute_input": "2026-05-19T20:35:51.329567Z", "iopub.status.busy": "2026-05-19T20:35:51.329451Z", "iopub.status.idle": "2026-05-19T20:35:51.538487Z", "shell.execute_reply": "2026-05-19T20:35:51.537857Z" } }, "outputs": [], "source": [ "obs->setVal(160.);\n", "RooDataSet dataOrig{\"exampleData\", \"exampleData\", *obs};\n", "dataOrig.add(*obs);\n", "\n", "RooArgSet all(*s, *ratioBkgEff, *ratioSigEff);" ] }, { "cell_type": "markdown", "id": "d06aa44d", "metadata": {}, "source": [ "not necessary" ] }, { "cell_type": "code", "execution_count": 6, "id": "e8603255", "metadata": { "collapsed": false, "execution": { "iopub.execute_input": "2026-05-19T20:35:51.540500Z", "iopub.status.busy": "2026-05-19T20:35:51.540379Z", "iopub.status.idle": "2026-05-19T20:35:51.750723Z", "shell.execute_reply": "2026-05-19T20:35:51.750268Z" } }, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "[#1] INFO:Minimization -- Including the following constraint terms in minimization: (sigConstraint,bkgConstraint)\n", "[#1] INFO:Minimization -- The global observables are not defined , normalize constraints with respect to the parameters (ratioSigEff,ratioBkgEff)\n", "[#1] INFO:Fitting -- RooAbsPdf::fitTo(modelWithConstraints) fixing normalization set for coefficient determination to observables in data\n", "[#1] INFO:Fitting -- using generic CPU library compiled with no vectorizations\n", "[#1] INFO:Fitting -- Creation of NLL object took 1.02969 ms\n", "[#1] INFO:Fitting -- RooAddition::defaultErrorLevel(nll_modelWithConstraints_exampleData) Summation contains a RooNLLVar, using its error level\n", "[#1] INFO:Minimization -- [fitFCN] No discrete parameters, performing continuous minimization only\n" ] } ], "source": [ "modelWithConstraints->fitTo(dataOrig, Constrain({*ratioSigEff, *ratioBkgEff}), PrintLevel(-1));" ] }, { "cell_type": "markdown", "id": "7ad3eaa2", "metadata": {}, "source": [ "Now let's make some confidence intervals for s, our parameter of interest" ] }, { "cell_type": "code", "execution_count": 7, "id": "74f4cb9e", "metadata": { "collapsed": false, "execution": { "iopub.execute_input": "2026-05-19T20:35:51.757588Z", "iopub.status.busy": "2026-05-19T20:35:51.757443Z", "iopub.status.idle": "2026-05-19T20:35:51.963415Z", "shell.execute_reply": "2026-05-19T20:35:51.963007Z" } }, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "[#1] INFO:ObjectHandling -- RooWorkspace::import() importing dataset exampleData\n" ] } ], "source": [ "RooArgSet paramOfInterest(*s);\n", "\n", "ModelConfig modelConfig(&wspace);\n", "modelConfig.SetPdf(*modelWithConstraints);\n", "modelConfig.SetParametersOfInterest(paramOfInterest);\n", "modelConfig.SetNuisanceParameters(constrainedParams);\n", "modelConfig.SetObservables(*obs);\n", "modelConfig.SetGlobalObservables(RooArgSet(*gSigEff, *gSigBkg));\n", "modelConfig.SetName(\"ModelConfig\");\n", "wspace.import(modelConfig);\n", "wspace.import(dataOrig);\n", "wspace.SetName(\"w\");" ] }, { "cell_type": "markdown", "id": "f5e87725", "metadata": {}, "source": [ "wspace.writeToFile(\"rs101_ws.root\");" ] }, { "cell_type": "markdown", "id": "f98acc55", "metadata": {}, "source": [ "Make sure we reference the data in the workspace from now on" ] }, { "cell_type": "code", "execution_count": 8, "id": "b44a4d6b", "metadata": { "collapsed": false, "execution": { "iopub.execute_input": "2026-05-19T20:35:51.964795Z", "iopub.status.busy": "2026-05-19T20:35:51.964674Z", "iopub.status.idle": "2026-05-19T20:35:52.173930Z", "shell.execute_reply": "2026-05-19T20:35:52.173471Z" } }, "outputs": [ { "name": "stderr", "output_type": "stream", "text": [ "input_line_77:2:2: warning: 'data' shadows a declaration with the same name in the 'std' namespace; use '::data' to reference this declaration\n", " RooDataSet &data = static_cast(*wspace.data(dataOrig.GetName()));\n", " ^\n" ] } ], "source": [ "RooDataSet &data = static_cast(*wspace.data(dataOrig.GetName()));" ] }, { "cell_type": "markdown", "id": "5594ed75", "metadata": {}, "source": [ "First, let's use a Calculator based on the Profile Likelihood Ratio\n", "ProfileLikelihoodCalculator plc(data, *modelWithConstraints, paramOfInterest);" ] }, { "cell_type": "code", "execution_count": 9, "id": "8f245600", "metadata": { "collapsed": false, "execution": { "iopub.execute_input": "2026-05-19T20:35:52.175184Z", "iopub.status.busy": "2026-05-19T20:35:52.175071Z", "iopub.status.idle": "2026-05-19T20:35:52.385360Z", "shell.execute_reply": "2026-05-19T20:35:52.384893Z" } }, "outputs": [ { "name": "stderr", "output_type": "stream", "text": [ "input_line_78:2:34: error: reference to 'data' is ambiguous\n", " ProfileLikelihoodCalculator plc(data, modelConfig);\n", " ^\n", "input_line_77:2:14: note: candidate found by name lookup is 'data'\n", " RooDataSet &data = static_cast(*wspace.data(dataOrig.GetName()));\n", " ^\n", "/usr/lib/gcc/x86_64-redhat-linux/14/../../../../include/c++/14/bits/range_access.h:344:5: note: candidate found by name lookup is 'std::data'\n", " data(initializer_list<_Tp> __il) noexcept\n", " ^\n", "/usr/lib/gcc/x86_64-redhat-linux/14/../../../../include/c++/14/bits/range_access.h:312:5: note: candidate found by name lookup is 'std::data'\n", " data(_Container& __cont) noexcept(noexcept(__cont.data()))\n", " ^\n", "/usr/lib/gcc/x86_64-redhat-linux/14/../../../../include/c++/14/bits/range_access.h:323:5: note: candidate found by name lookup is 'std::data'\n", " data(const _Container& __cont) noexcept(noexcept(__cont.data()))\n", " ^\n", "/usr/lib/gcc/x86_64-redhat-linux/14/../../../../include/c++/14/bits/range_access.h:334:5: note: candidate found by name lookup is 'std::data'\n", " data(_Tp (&__array)[_Nm]) noexcept\n", " ^\n", "input_line_78:2:34: error: unknown type name 'data'\n", " ProfileLikelihoodCalculator plc(data, modelConfig);\n", " ^\n", "input_line_78:2:40: error: unknown type name 'modelConfig'\n", " ProfileLikelihoodCalculator plc(data, modelConfig);\n", " ^\n" ] } ], "source": [ "ProfileLikelihoodCalculator plc(data, modelConfig);\n", "plc.SetTestSize(.05);\n", "std::unique_ptr lrinterval{static_cast(plc.GetInterval())};" ] }, { "cell_type": "markdown", "id": "60569314", "metadata": {}, "source": [ "Let's make a plot" ] }, { "cell_type": "code", "execution_count": 10, "id": "505544f9", "metadata": { "collapsed": false, "execution": { "iopub.execute_input": "2026-05-19T20:35:52.386717Z", "iopub.status.busy": "2026-05-19T20:35:52.386574Z", "iopub.status.idle": "2026-05-19T20:35:52.801865Z", "shell.execute_reply": "2026-05-19T20:35:52.792664Z" } }, "outputs": [ { "name": "stderr", "output_type": "stream", "text": [ "[runStaticInitializersOnce]: Failed to materialize symbols: { (main, { _ZN12__cling_N52710dataCanvasE, cling_module_320_.11, cling_module_320_.13, _ZN12__cling_N52716__cling_Un1Qu338EPv, $.cling-module-320.__inits.0, __vd_init_order__cling_Un1Qu39, _ZNK5cling7runtime8internal15LifetimeHandler9getMemoryEv, _GLOBAL__sub_I_cling_module_320, __orc_init_func.cling-module-320, _ZN5cling7runtime8internal15DynamicExprInfoC2EPKcPPvb, cling_module_320_.12, _ZN12__cling_N52724__dynamic__cling_Un1Qu30E, _ZN5cling7runtime8internal15DynamicExprInfoC1EPKcPPvb, _Z30__fd_init_order__cling_Un1Qu38v, cling_module_320_, _ZN12__cling_N5277plotIntE }) }\n", "IncrementalExecutor::executeFunction: symbol '_ZN5cling7runtime8internal15LifetimeHandlerC1EPNS1_15DynamicExprInfoEPN5clang11DeclContextEPKcPNS_11InterpreterE' unresolved while linking [cling interface function]!\n", "You are probably missing the definition of cling::runtime::internal::LifetimeHandler::LifetimeHandler(cling::runtime::internal::DynamicExprInfo*, clang::DeclContext*, char const*, cling::Interpreter*)\n", "Maybe you need to load the corresponding shared library?\n", "IncrementalExecutor::executeFunction: symbol '_ZN5cling7runtime8internal15LifetimeHandlerD1Ev' unresolved while linking [cling interface function]!\n", "You are probably missing the definition of cling::runtime::internal::LifetimeHandler::~LifetimeHandler()\n", "Maybe you need to load the corresponding shared library?\n" ] } ], "source": [ "auto dataCanvas = new TCanvas(\"dataCanvas\");\n", "dataCanvas->Divide(2, 1);\n", "\n", "dataCanvas->cd(1);\n", "LikelihoodIntervalPlot plotInt(lrinterval.get());\n", "plotInt.SetTitle(\"Profile Likelihood Ratio and Posterior for S\");\n", "plotInt.Draw();" ] }, { "cell_type": "markdown", "id": "5e01fcd2", "metadata": {}, "source": [ "Second, use a Calculator based on the Feldman Cousins technique" ] }, { "cell_type": "code", "execution_count": 11, "id": "c9d3da0d", "metadata": { "collapsed": false, "execution": { "iopub.execute_input": "2026-05-19T20:35:52.805038Z", "iopub.status.busy": "2026-05-19T20:35:52.804907Z", "iopub.status.idle": "2026-05-19T20:35:53.014975Z", "shell.execute_reply": "2026-05-19T20:35:53.014509Z" } }, "outputs": [ { "name": "stderr", "output_type": "stream", "text": [ "input_line_80:2:20: error: reference to 'data' is ambiguous\n", " FeldmanCousins fc(data, modelConfig);\n", " ^\n", "input_line_77:2:14: note: candidate found by name lookup is 'data'\n", " RooDataSet &data = static_cast(*wspace.data(dataOrig.GetName()));\n", " ^\n", "/usr/lib/gcc/x86_64-redhat-linux/14/../../../../include/c++/14/bits/range_access.h:344:5: note: candidate found by name lookup is 'std::data'\n", " data(initializer_list<_Tp> __il) noexcept\n", " ^\n", "/usr/lib/gcc/x86_64-redhat-linux/14/../../../../include/c++/14/bits/range_access.h:312:5: note: candidate found by name lookup is 'std::data'\n", " data(_Container& __cont) noexcept(noexcept(__cont.data()))\n", " ^\n", "/usr/lib/gcc/x86_64-redhat-linux/14/../../../../include/c++/14/bits/range_access.h:323:5: note: candidate found by name lookup is 'std::data'\n", " data(const _Container& __cont) noexcept(noexcept(__cont.data()))\n", " ^\n", "/usr/lib/gcc/x86_64-redhat-linux/14/../../../../include/c++/14/bits/range_access.h:334:5: note: candidate found by name lookup is 'std::data'\n", " data(_Tp (&__array)[_Nm]) noexcept\n", " ^\n", "input_line_80:2:20: error: unknown type name 'data'\n", " FeldmanCousins fc(data, modelConfig);\n", " ^\n", "input_line_80:2:26: error: unknown type name 'modelConfig'\n", " FeldmanCousins fc(data, modelConfig);\n", " ^\n" ] } ], "source": [ "FeldmanCousins fc(data, modelConfig);\n", "fc.UseAdaptiveSampling(true);\n", "fc.FluctuateNumDataEntries(false); // number counting analysis: dataset always has 1 entry with N events observed\n", "fc.SetNBins(100); // number of points to test per parameter\n", "fc.SetTestSize(.05);" ] }, { "cell_type": "markdown", "id": "a767c15d", "metadata": {}, "source": [ "fc.SaveBeltToFile(true); // optional" ] }, { "cell_type": "code", "execution_count": 12, "id": "45180aaa", "metadata": { "collapsed": false, "execution": { "iopub.execute_input": "2026-05-19T20:35:53.016922Z", "iopub.status.busy": "2026-05-19T20:35:53.016806Z", "iopub.status.idle": "2026-05-19T20:35:53.226850Z", "shell.execute_reply": "2026-05-19T20:35:53.226278Z" } }, "outputs": [ { "name": "stderr", "output_type": "stream", "text": [ "input_line_81:4:63: error: reference to 'data' is ambiguous\n", "std::unique_ptr fit{modelWithConstraints->fitTo(data, Save(true), PrintLevel(-1))};\n", " ^\n", "input_line_77:2:14: note: candidate found by name lookup is 'data'\n", " RooDataSet &data = static_cast(*wspace.data(dataOrig.GetName()));\n", " ^\n", "/usr/lib/gcc/x86_64-redhat-linux/14/../../../../include/c++/14/bits/range_access.h:344:5: note: candidate found by name lookup is 'std::data'\n", " data(initializer_list<_Tp> __il) noexcept\n", " ^\n", "/usr/lib/gcc/x86_64-redhat-linux/14/../../../../include/c++/14/bits/range_access.h:312:5: note: candidate found by name lookup is 'std::data'\n", " data(_Container& __cont) noexcept(noexcept(__cont.data()))\n", " ^\n", "/usr/lib/gcc/x86_64-redhat-linux/14/../../../../include/c++/14/bits/range_access.h:323:5: note: candidate found by name lookup is 'std::data'\n", " data(const _Container& __cont) noexcept(noexcept(__cont.data()))\n", " ^\n", "/usr/lib/gcc/x86_64-redhat-linux/14/../../../../include/c++/14/bits/range_access.h:334:5: note: candidate found by name lookup is 'std::data'\n", " data(_Tp (&__array)[_Nm]) noexcept\n", " ^\n" ] } ], "source": [ "std::unique_ptr fcint{static_cast(fc.GetInterval())};\n", "\n", "std::unique_ptr fit{modelWithConstraints->fitTo(data, Save(true), PrintLevel(-1))};" ] }, { "cell_type": "markdown", "id": "e38308d9", "metadata": {}, "source": [ "Third, use a Calculator based on Markov Chain monte carlo\n", "Before configuring the calculator, let's make a ProposalFunction\n", "that will achieve a high acceptance rate" ] }, { "cell_type": "code", "execution_count": 13, "id": "daad807c", "metadata": { "collapsed": false, "execution": { "iopub.execute_input": "2026-05-19T20:35:53.228425Z", "iopub.status.busy": "2026-05-19T20:35:53.228308Z", "iopub.status.idle": "2026-05-19T20:35:53.438935Z", "shell.execute_reply": "2026-05-19T20:35:53.438354Z" } }, "outputs": [ { "name": "stderr", "output_type": "stream", "text": [ "input_line_82:9:19: error: reference to 'data' is ambiguous\n", "MCMCCalculator mc(data, modelConfig);\n", " ^\n", "input_line_77:2:14: note: candidate found by name lookup is 'data'\n", " RooDataSet &data = static_cast(*wspace.data(dataOrig.GetName()));\n", " ^\n", "/usr/lib/gcc/x86_64-redhat-linux/14/../../../../include/c++/14/bits/range_access.h:344:5: note: candidate found by name lookup is 'std::data'\n", " data(initializer_list<_Tp> __il) noexcept\n", " ^\n", "/usr/lib/gcc/x86_64-redhat-linux/14/../../../../include/c++/14/bits/range_access.h:312:5: note: candidate found by name lookup is 'std::data'\n", " data(_Container& __cont) noexcept(noexcept(__cont.data()))\n", " ^\n", "/usr/lib/gcc/x86_64-redhat-linux/14/../../../../include/c++/14/bits/range_access.h:323:5: note: candidate found by name lookup is 'std::data'\n", " data(const _Container& __cont) noexcept(noexcept(__cont.data()))\n", " ^\n", "/usr/lib/gcc/x86_64-redhat-linux/14/../../../../include/c++/14/bits/range_access.h:334:5: note: candidate found by name lookup is 'std::data'\n", " data(_Tp (&__array)[_Nm]) noexcept\n", " ^\n", "input_line_82:9:19: error: unknown type name 'data'\n", "MCMCCalculator mc(data, modelConfig);\n", " ^\n", "input_line_82:9:25: error: unknown type name 'modelConfig'\n", "MCMCCalculator mc(data, modelConfig);\n", " ^\n" ] } ], "source": [ "ProposalHelper ph;\n", "ph.SetVariables((RooArgSet &)fit->floatParsFinal());\n", "ph.SetCovMatrix(fit->covarianceMatrix());\n", "ph.SetUpdateProposalParameters(true);\n", "ph.SetCacheSize(100);\n", "ProposalFunction *pdfProp = ph.GetProposalFunction();\n", "\n", "MCMCCalculator mc(data, modelConfig);\n", "mc.SetNumIters(20000); // steps to propose in the chain\n", "mc.SetTestSize(.05); // 95% CL\n", "mc.SetNumBurnInSteps(40); // ignore first N steps in chain as \"burn in\"\n", "mc.SetProposalFunction(*pdfProp);\n", "mc.SetLeftSideTailFraction(0.5); // find a \"central\" interval\n", "std::unique_ptr mcInt{static_cast(mc.GetInterval())};" ] }, { "cell_type": "markdown", "id": "70b861b6", "metadata": {}, "source": [ "Get Lower and Upper limits from Profile Calculator" ] }, { "cell_type": "code", "execution_count": 14, "id": "a1048e83", "metadata": { "collapsed": false, "execution": { "iopub.execute_input": "2026-05-19T20:35:53.450052Z", "iopub.status.busy": "2026-05-19T20:35:53.449913Z", "iopub.status.idle": "2026-05-19T20:35:53.664739Z", "shell.execute_reply": "2026-05-19T20:35:53.664386Z" } }, "outputs": [ { "name": "stderr", "output_type": "stream", "text": [ "[runStaticInitializersOnce]: Failed to materialize symbols: { (main, { __orc_init_func.cling-module-320 }) }\n", "cling JIT session error: Failed to materialize symbols: { (main, { _ZN5cling7runtime8internal15DynamicExprInfoC1EPKcPPvb }) }\n" ] } ], "source": [ "std::cout << \"Profile lower limit on s = \" << lrinterval->LowerLimit(*s) << std::endl;\n", "std::cout << \"Profile upper limit on s = \" << lrinterval->UpperLimit(*s) << std::endl;" ] }, { "cell_type": "markdown", "id": "972c1127", "metadata": {}, "source": [ "Get Lower and Upper limits from FeldmanCousins with profile construction" ] }, { "cell_type": "code", "execution_count": 15, "id": "9c396e88", "metadata": { "collapsed": false, "execution": { "iopub.execute_input": "2026-05-19T20:35:53.672444Z", "iopub.status.busy": "2026-05-19T20:35:53.672320Z", "iopub.status.idle": "2026-05-19T20:35:53.874888Z", "shell.execute_reply": "2026-05-19T20:35:53.874329Z" } }, "outputs": [ { "name": "stderr", "output_type": "stream", "text": [ "cling JIT session error: Failed to materialize symbols: { (main, { _ZN12__cling_N52710dataCanvasE }) }\n" ] } ], "source": [ "if (fcint) {\n", " double fcul = fcint->UpperLimit(*s);\n", " double fcll = fcint->LowerLimit(*s);\n", " std::cout << \"FC lower limit on s = \" << fcll << std::endl;\n", " std::cout << \"FC upper limit on s = \" << fcul << std::endl;\n", " auto fcllLine = new TLine(fcll, 0, fcll, 1);\n", " auto fculLine = new TLine(fcul, 0, fcul, 1);\n", " fcllLine->SetLineColor(kRed);\n", " fculLine->SetLineColor(kRed);\n", " fcllLine->Draw(\"same\");\n", " fculLine->Draw(\"same\");\n", " dataCanvas->Update();\n", "}" ] }, { "cell_type": "markdown", "id": "8bca5f12", "metadata": {}, "source": [ "Plot MCMC interval and print some statistics" ] }, { "cell_type": "code", "execution_count": 16, "id": "bf944ad2", "metadata": { "collapsed": false, "execution": { "iopub.execute_input": "2026-05-19T20:35:53.876190Z", "iopub.status.busy": "2026-05-19T20:35:53.876075Z", "iopub.status.idle": "2026-05-19T20:35:54.086040Z", "shell.execute_reply": "2026-05-19T20:35:54.085479Z" } }, "outputs": [ { "name": "stderr", "output_type": "stream", "text": [ "cling JIT session error: Failed to materialize symbols: { (main, { _ZNK5cling7runtime8internal15LifetimeHandler9getMemoryEv }) }\n", "[runStaticInitializersOnce]: Failed to materialize symbols: { (main, { cling_module_326_, _ZN12__cling_N53316__cling_Un1Qu344EPv, _ZN12__cling_N5334mculE, cling_module_326_.17, cling_module_326_.15, cling_module_326_.16, _ZN12__cling_N5334mcllE, cling_module_326_.18, _ZNK12TColorNumber6numberEv, _ZN12TColorNumberC1Ei, __vd_init_order__cling_Un1Qu313, __orc_init_func.cling-module-326, $.cling-module-326.__inits.0, _ZN8RooStats16MCMCIntervalPlot12SetLineColorE12TColorNumber, _ZN8RooStats16MCMCIntervalPlot12SetLineWidthEi, _ZN12TColorNumberC2Ei, _ZN12__cling_N53324__dynamic__cling_Un1Qu32E, _GLOBAL__sub_I_cling_module_326, _Z31__fd_init_order__cling_Un1Qu312v, _ZN12__cling_N5336mcPlotE }) }\n" ] } ], "source": [ "MCMCIntervalPlot mcPlot(*mcInt);\n", "mcPlot.SetLineColor(kMagenta);\n", "mcPlot.SetLineWidth(2);\n", "mcPlot.Draw(\"same\");\n", "\n", "double mcul = mcInt->UpperLimit(*s);\n", "double mcll = mcInt->LowerLimit(*s);\n", "std::cout << \"MCMC lower limit on s = \" << mcll << std::endl;\n", "std::cout << \"MCMC upper limit on s = \" << mcul << std::endl;\n", "std::cout << \"MCMC Actual confidence level: \" << mcInt->GetActualConfidenceLevel() << std::endl;" ] }, { "cell_type": "markdown", "id": "de26da70", "metadata": {}, "source": [ "3-d plot of the parameter points" ] }, { "cell_type": "code", "execution_count": 17, "id": "08305ac4", "metadata": { "collapsed": false, "execution": { "iopub.execute_input": "2026-05-19T20:35:54.087330Z", "iopub.status.busy": "2026-05-19T20:35:54.087216Z", "iopub.status.idle": "2026-05-19T20:35:54.297222Z", "shell.execute_reply": "2026-05-19T20:35:54.296687Z" } }, "outputs": [ { "name": "stderr", "output_type": "stream", "text": [ "[runStaticInitializersOnce]: Failed to materialize symbols: { (main, { __orc_init_func.cling-module-326 }) }\n", "cling JIT session error: Failed to materialize symbols: { (main, { _ZN12__cling_N52710dataCanvasE }) }\n" ] } ], "source": [ "dataCanvas->cd(2);" ] }, { "cell_type": "markdown", "id": "364a3b0d", "metadata": {}, "source": [ "also plot the points in the markov chain" ] }, { "cell_type": "code", "execution_count": 18, "id": "684437f7", "metadata": { "collapsed": false, "execution": { "iopub.execute_input": "2026-05-19T20:35:54.298488Z", "iopub.status.busy": "2026-05-19T20:35:54.298365Z", "iopub.status.idle": "2026-05-19T20:35:54.510975Z", "shell.execute_reply": "2026-05-19T20:35:54.510548Z" } }, "outputs": [ { "name": "stderr", "output_type": "stream", "text": [ "cling JIT session error: Failed to materialize symbols: { (main, { _ZNK5cling7runtime8internal15LifetimeHandler9getMemoryEv }) }\n", "[runStaticInitializersOnce]: Failed to materialize symbols: { (main, { _ZN12__cling_N5359chainDataE, $.cling-module-328.__inits.0, _ZSt3getILm0EJP10RooDataSetSt14default_deleteIS0_EEERKNSt13tuple_elementIXT_ESt5tupleIJDpT0_EEE4typeERKS8_, _ZN12__cling_N53516__cling_Un1Qu346EPv, __vd_init_order__cling_Un1Qu315, _ZN12__cling_N53524__dynamic__cling_Un1Qu33E, cling_module_328_.15, _ZNKSt10unique_ptrI10RooDataSetSt14default_deleteIS0_EEdeEv, _ZNKSt15__uniq_ptr_implI10RooDataSetSt14default_deleteIS0_EE6_M_ptrEv, _ZNSt11_Tuple_implILm0EJP10RooDataSetSt14default_deleteIS0_EEE7_M_headERKS4_, cling_module_328_.16, cling_module_328_.17, _GLOBAL__sub_I_cling_module_328, cling_module_328_, _ZNKSt10unique_ptrI10RooDataSetSt14default_deleteIS0_EEptEv, _ZNSt10_Head_baseILm0EP10RooDataSetLb0EE7_M_headERKS2_, _Z31__fd_init_order__cling_Un1Qu314v, _ZSt12__get_helperILm0EP10RooDataSetJSt14default_deleteIS0_EEERKT0_RKSt11_Tuple_implIXT_EJS4_DpT1_EE, _ZNKSt10unique_ptrI10RooDataSetSt14default_deleteIS0_EE3getEv, _ZN12__cling_N5355chainE, __orc_init_func.cling-module-328 }) }\n" ] } ], "source": [ "std::unique_ptr chainData{mcInt->GetChainAsDataSet()};\n", "\n", "assert(chainData);\n", "std::cout << \"plotting the chain data - nentries = \" << chainData->numEntries() << std::endl;\n", "TTree *chain = RooStats::GetAsTTree(\"chainTreeData\", \"chainTreeData\", *chainData);\n", "assert(chain);\n", "chain->SetMarkerStyle(6);\n", "chain->SetMarkerColor(kRed);\n", "\n", "chain->Draw(\"s:ratioSigEff:ratioBkgEff\", \"nll_MarkovChain_local_\", \"box\"); // 3-d box proportional to posterior" ] }, { "cell_type": "markdown", "id": "97a7ef5d", "metadata": {}, "source": [ "the points used in the profile construction" ] }, { "cell_type": "code", "execution_count": 19, "id": "04a61b87", "metadata": { "collapsed": false, "execution": { "iopub.execute_input": "2026-05-19T20:35:54.529230Z", "iopub.status.busy": "2026-05-19T20:35:54.529099Z", "iopub.status.idle": "2026-05-19T20:35:54.739832Z", "shell.execute_reply": "2026-05-19T20:35:54.739133Z" } }, "outputs": [ { "name": "stderr", "output_type": "stream", "text": [ "cling JIT session error: Failed to materialize symbols: { (main, { _ZNK5cling5Value21isPointerOrObjectTypeEv }) }\n", "[runStaticInitializersOnce]: Failed to materialize symbols: { (main, { $.cling-module-329.__inits.0, _ZN5cling5Value7CastFwdIP10RooDataSetE4castERKS0_, _ZNK5cling5Value6castAsIP10RooDataSetEET_v, _ZN12__cling_N53616__cling_Un1Qu347EPv, cling_module_329_, __orc_init_func.cling-module-329, _GLOBAL__sub_I_cling_module_329, _ZN5cling7runtime8internal9EvaluateTIP10RooDataSetEET_PNS1_15DynamicExprInfoEPN5clang11DeclContextE, _ZN12__cling_N53611parScanDataE }) }\n" ] } ], "source": [ "RooDataSet *parScanData = (RooDataSet *)fc.GetPointsToScan();\n", "assert(parScanData);\n", "std::cout << \"plotting the scanned points used in the frequentist construction - npoints = \"\n", " << parScanData->numEntries() << std::endl;" ] }, { "cell_type": "markdown", "id": "bb7bc203", "metadata": {}, "source": [ "getting the tree and drawing it -crashes (very strange....);\n", "TTree* parameterScan = RooStats::GetAsTTree(\"parScanTreeData\",\"parScanTreeData\",*parScanData);\n", "assert(parameterScan);\n", "parameterScan->Draw(\"s:ratioSigEff:ratioBkgEff\",\"\",\"goff\");" ] }, { "cell_type": "code", "execution_count": 20, "id": "19e0c57a", "metadata": { "collapsed": false, "execution": { "iopub.execute_input": "2026-05-19T20:35:54.741173Z", "iopub.status.busy": "2026-05-19T20:35:54.741056Z", "iopub.status.idle": "2026-05-19T20:35:54.955744Z", "shell.execute_reply": "2026-05-19T20:35:54.955338Z" } }, "outputs": [ { "name": "stderr", "output_type": "stream", "text": [ "cling JIT session error: Failed to materialize symbols: { (main, { _ZN12__cling_N53611parScanDataE }) }\n", "[runStaticInitializersOnce]: Failed to materialize symbols: { (main, { _ZN12__cling_N53716__cling_Un1Qu348EPv, $.cling-module-330.__inits.0, _GLOBAL__sub_I_cling_module_330, cling_module_330_, __orc_init_func.cling-module-330, _ZN12__cling_N5372grE }) }\n" ] } ], "source": [ "auto gr = new TGraph2D(parScanData->numEntries());\n", "for (int ievt = 0; ievt < parScanData->numEntries(); ++ievt) {\n", " const RooArgSet *evt = parScanData->get(ievt);\n", " double x = evt->getRealValue(\"ratioBkgEff\");\n", " double y = evt->getRealValue(\"ratioSigEff\");\n", " double z = evt->getRealValue(\"s\");\n", " gr->SetPoint(ievt, x, y, z);\n", " // std::cout << ievt << \" \" << x << \" \" << y << \" \" << z << std::endl;\n", "}\n", "gr->SetMarkerStyle(24);\n", "gr->Draw(\"P SAME\");" ] }, { "cell_type": "markdown", "id": "dccf24f5", "metadata": {}, "source": [ "print timing info" ] }, { "cell_type": "code", "execution_count": 21, "id": "b7820d36", "metadata": { "collapsed": false, "execution": { "iopub.execute_input": "2026-05-19T20:35:54.971276Z", "iopub.status.busy": "2026-05-19T20:35:54.971135Z", "iopub.status.idle": "2026-05-19T20:35:55.189007Z", "shell.execute_reply": "2026-05-19T20:35:55.188406Z" } }, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "Real time 0:00:04, CP time 0.990\n" ] }, { "name": "stderr", "output_type": "stream", "text": [ "[runStaticInitializersOnce]: Failed to materialize symbols: { (main, { __orc_init_func.cling-module-328 }) }\n" ] } ], "source": [ "t.Stop();\n", "t.Print();" ] }, { "cell_type": "markdown", "id": "f04aecba", "metadata": {}, "source": [ "Draw all canvases " ] }, { "cell_type": "code", "execution_count": 22, "id": "a0ccad3f", "metadata": { "collapsed": false, "execution": { "iopub.execute_input": "2026-05-19T20:35:55.190564Z", "iopub.status.busy": "2026-05-19T20:35:55.190437Z", "iopub.status.idle": "2026-05-19T20:35:55.399676Z", "shell.execute_reply": "2026-05-19T20:35:55.398990Z" } }, "outputs": [], "source": [ "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 }