{
"cells": [
{
"cell_type": "markdown",
"id": "df40a5eb",
"metadata": {},
"source": [
"# rs401d_FeldmanCousins\n",
"Neutrino Oscillation Example from Feldman & Cousins\n",
"\n",
"This tutorial shows a more complex example using the FeldmanCousins utility\n",
"to create a confidence interval for a toy neutrino oscillation experiment.\n",
"The example attempts to faithfully reproduce the toy example described in Feldman & Cousins'\n",
"original paper, Phys.Rev.D57:3873-3889,1998.\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": "a2d35b24",
"metadata": {
"collapsed": false,
"execution": {
"iopub.execute_input": "2026-05-19T20:36:02.553991Z",
"iopub.status.busy": "2026-05-19T20:36:02.553876Z",
"iopub.status.idle": "2026-05-19T20:36:02.570956Z",
"shell.execute_reply": "2026-05-19T20:36:02.570370Z"
}
},
"outputs": [],
"source": [
"%%cpp -d\n",
"#include \"RooGlobalFunc.h\"\n",
"#include \"RooStats/ConfInterval.h\"\n",
"#include \"RooStats/FeldmanCousins.h\"\n",
"#include \"RooStats/ProfileLikelihoodCalculator.h\"\n",
"#include \"RooStats/MCMCCalculator.h\"\n",
"#include \"RooStats/UniformProposal.h\"\n",
"#include \"RooStats/LikelihoodIntervalPlot.h\"\n",
"#include \"RooStats/MCMCIntervalPlot.h\"\n",
"#include \"RooStats/MCMCInterval.h\"\n",
"\n",
"#include \"RooDataSet.h\"\n",
"#include \"RooDataHist.h\"\n",
"#include \"RooRealVar.h\"\n",
"#include \"RooConstVar.h\"\n",
"#include \"RooAddition.h\"\n",
"#include \"RooProduct.h\"\n",
"#include \"RooProdPdf.h\"\n",
"#include \"RooAddPdf.h\"\n",
"\n",
"#include \"TROOT.h\"\n",
"#include \"RooPolynomial.h\"\n",
"#include \"RooRandom.h\"\n",
"\n",
"#include \"RooProfileLL.h\"\n",
"\n",
"#include \"RooPlot.h\"\n",
"\n",
"#include \"TCanvas.h\"\n",
"#include \"TH1F.h\"\n",
"#include \"TH2F.h\"\n",
"#include \"TTree.h\"\n",
"#include \"TMarker.h\"\n",
"#include \"TStopwatch.h\"\n",
"\n",
"#include \n",
"\n",
"using namespace RooFit;\n",
"using namespace RooStats;"
]
},
{
"cell_type": "markdown",
"id": "fae37fb4",
"metadata": {},
"source": [
" Arguments are defined. "
]
},
{
"cell_type": "code",
"execution_count": 2,
"id": "ee5f3921",
"metadata": {
"collapsed": false,
"execution": {
"iopub.execute_input": "2026-05-19T20:36:02.572696Z",
"iopub.status.busy": "2026-05-19T20:36:02.572563Z",
"iopub.status.idle": "2026-05-19T20:36:02.920078Z",
"shell.execute_reply": "2026-05-19T20:36:02.918396Z"
}
},
"outputs": [],
"source": [
"bool doFeldmanCousins = false;\n",
"bool doMCMC = true;"
]
},
{
"cell_type": "code",
"execution_count": 3,
"id": "0c585097",
"metadata": {
"collapsed": false,
"execution": {
"iopub.execute_input": "2026-05-19T20:36:02.932423Z",
"iopub.status.busy": "2026-05-19T20:36:02.932286Z",
"iopub.status.idle": "2026-05-19T20:36:03.159023Z",
"shell.execute_reply": "2026-05-19T20:36:03.152361Z"
}
},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"[#1] INFO:InputArguments -- Silent clipping of values to range in `RooRealVar::setVal()` enabled.\n"
]
}
],
"source": [
"RooRealVar::enableSilentClipping();"
]
},
{
"cell_type": "markdown",
"id": "dde21f4a",
"metadata": {},
"source": [
"to time the macro"
]
},
{
"cell_type": "code",
"execution_count": 4,
"id": "6e44154b",
"metadata": {
"collapsed": false,
"execution": {
"iopub.execute_input": "2026-05-19T20:36:03.160799Z",
"iopub.status.busy": "2026-05-19T20:36:03.160673Z",
"iopub.status.idle": "2026-05-19T20:36:03.369450Z",
"shell.execute_reply": "2026-05-19T20:36:03.369012Z"
}
},
"outputs": [],
"source": [
"TStopwatch t;\n",
"t.Start();"
]
},
{
"cell_type": "markdown",
"id": "f2a3a613",
"metadata": {},
"source": [
"Taken from Feldman & Cousins paper, Phys.Rev.D57:3873-3889,1998.\n",
"e-Print: physics/9711021 (see page 13.)\n",
"\n",
"Quantum mechanics dictates that the probability of such a transformation is given by the formula\n",
"$P (\\nu\\mu \\rightarrow \\nu e ) = sin^2 (2\\theta) sin^2 (1.27 \\Delta m^2 L /E )$\n",
"where P is the probability for a $\\nu\\mu$ to transform into a $\\nu e$ , L is the distance in km between\n",
"the creation of the neutrino from meson decay and its interaction in the detector, E is the\n",
"neutrino energy in GeV, and $\\Delta m^2 = |m^2 - m^2 |$ in $(eV/c^2 )^2$ .\n",
"\n",
"To demonstrate how this works in practice, and how it compares to alternative approaches\n",
"that have been used, we consider a toy model of a typical neutrino oscillation experiment.\n",
"The toy model is defined by the following parameters: Mesons are assumed to decay to\n",
"neutrinos uniformly in a region 600 m to 1000 m from the detector. The expected background\n",
"from conventional $\\nu e$ interactions and misidentified $\\nu\\mu$ interactions is assumed to be 100\n",
"events in each of 5 energy bins which span the region from 10 to 60 GeV. We assume that\n",
"the $\\nu\\mu$ flux is such that if $P (\\nu\\mu \\rightarrow \\nu e ) = 0.01$ averaged over any bin, then that bin\n",
"would\n",
"have an expected additional contribution of 100 events due to $\\nu\\mu \\rightarrow \\nu e$ oscillations."
]
},
{
"cell_type": "markdown",
"id": "fcd58e8b",
"metadata": {},
"source": [
"Make signal model model"
]
},
{
"cell_type": "code",
"execution_count": 5,
"id": "566ebeb9",
"metadata": {
"collapsed": false,
"execution": {
"iopub.execute_input": "2026-05-19T20:36:03.375976Z",
"iopub.status.busy": "2026-05-19T20:36:03.375854Z",
"iopub.status.idle": "2026-05-19T20:36:03.594102Z",
"shell.execute_reply": "2026-05-19T20:36:03.592700Z"
}
},
"outputs": [],
"source": [
"RooRealVar E(\"E\", \"\", 15, 10, 60, \"GeV\");\n",
"RooRealVar L(\"L\", \"\", .800, .600, 1.0, \"km\"); // need these units in formula\n",
"RooRealVar deltaMSq(\"deltaMSq\", \"#Delta m^{2}\", 40, 1, 300, \"eV/c^{2}\");\n",
"RooRealVar sinSq2theta(\"sinSq2theta\", \"sin^{2}(2#theta)\", .006, .0, .02);"
]
},
{
"cell_type": "markdown",
"id": "b4cf5bd1",
"metadata": {},
"source": [
"RooRealVar deltaMSq(\"deltaMSq\",\"#Delta m^{2}\",40,20,70,\"eV/c^{2}\");\n",
"RooRealVar sinSq2theta(\"sinSq2theta\",\"sin^{2}(2#theta)\", .006,.001,.01);\n",
"PDF for oscillation only describes deltaMSq dependence, sinSq2theta goes into sigNorm"
]
},
{
"cell_type": "code",
"execution_count": 6,
"id": "352814db",
"metadata": {
"collapsed": false,
"execution": {
"iopub.execute_input": "2026-05-19T20:36:03.613412Z",
"iopub.status.busy": "2026-05-19T20:36:03.613269Z",
"iopub.status.idle": "2026-05-19T20:37:08.074407Z",
"shell.execute_reply": "2026-05-19T20:37:08.074109Z"
}
},
"outputs": [
{
"data": {
"text/html": [
"\n",
"\n",
"\n",
"
\n",
"\n",
"\n",
"\n"
],
"text/plain": [
""
]
},
"metadata": {},
"output_type": "display_data"
},
{
"name": "stdout",
"output_type": "stream",
"text": [
"[#1] INFO:NumericIntegration -- RooRealIntegral::init(PnmuTone_Int[L]) using numeric integrator RooIntegrator1D to calculate Int(L)\n",
"generate toy data with nEvents = 692\n",
"Minuit2Minimizer: Minimize with max-calls 1000 convergence for edm < 1 strategy 1\n",
"Minuit2Minimizer : Valid minimum - status = 0\n",
"FVAL = -1143.51305519945186\n",
"Edm = 2.62634022232556406e-06\n",
"Nfcn = 85\n",
"deltaMSq\t = 1.00008\t +/- 20.6319\t(limited)\n",
"sinSq2theta\t = 0.02\t +/- 0.000213424\t(limited)\n",
"[#1] INFO:Minimization -- p.d.f. provides expected number of events, including extended term in likelihood.\n",
"[#1] INFO:Fitting -- RooAbsPdf::fitTo(model) 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 396.962 μs\n",
"Metropolis-Hastings progress: ....................................................................................................\n",
"[#1] INFO:Eval -- Proposal acceptance rate: 0.4%\n",
"[#1] INFO:Eval -- Number of steps in chain: 20\n",
"[#0] ERROR:InputArguments -- MCMCInterval::CreateKeysPdf: creation of Keys PDF failed: Number of burn-in steps (num steps to ignore) >= number of steps in Markov chain.\n",
"[#0] WARNING:NumericIntegration -- RooAdaptiveIntegratorND::dtor(PnmuTone) WARNING: Number of suppressed warningings about integral evaluations where target precision was not reached is 633\n",
"Real time 0:00:36, CP time 36.040\n",
"MCMC actual confidence level: 0\n",
"[#0] ERROR:InputArguments -- MCMCInterval::CreateKeysPdf: creation of Keys PDF failed: Number of burn-in steps (num steps to ignore) >= number of steps in Markov chain.\n",
"[#0] ERROR:InputArguments -- MCMCIntervalPlot::DrawKeysPdfInterval: Couldn't get posterior Keys PDF.\n",
"[#1] INFO:Minimization -- RooProfileLL::evaluate(RooEvaluatorWrapper_Profile[deltaMSq,sinSq2theta]) Creating instance of MINUIT\n",
"[#1] INFO:Fitting -- RooAddition::defaultErrorLevel(nll_model_modelData) Summation contains a RooNLLVar, using its error level\n",
"[#1] INFO:Minimization -- RooProfileLL::evaluate(RooEvaluatorWrapper_Profile[deltaMSq,sinSq2theta]) determining minimum likelihood for current configurations w.r.t all observable\n",
"[#1] INFO:Minimization -- [fitFCN] No discrete parameters, performing continuous minimization only\n",
"[#1] INFO:Minimization -- RooProfileLL::evaluate(RooEvaluatorWrapper_Profile[deltaMSq,sinSq2theta]) minimum found at (deltaMSq=1, sinSq2theta=0.02)\n",
".[#1] INFO:Minimization -- [fitFCN] No discrete parameters, performing continuous minimization only\n",
".[#1] INFO:Minimization -- [fitFCN] No discrete parameters, performing continuous minimization only\n",
"[#1] INFO:Minimization -- LikelihoodInterval - Finding the contour of deltaMSq ( 0 ) and sinSq2theta ( 1 ) \n",
"Real time 0:01:04, CP time 64.310\n"
]
},
{
"name": "stderr",
"output_type": "stream",
"text": [
"input_line_53:62:1: warning: 'data' shadows a declaration with the same name in the 'std' namespace; use '::data' to reference this declaration\n",
"RooDataSet *data = model.generate(RooArgSet(E), nEventsData);\n",
"^\n",
"Info in : MnSeedGenerator Computing seed using NumericalGradient calculator\n",
"Info in : MnSeedGenerator Evaluated function and gradient in 7.80153 ms\n",
"Info in : MnSeedGenerator Initial state: FCN = -1090.831711 Edm = 1095.465905 NCalls = 9\n",
"Info in : NegativeG2LineSearch Doing a NegativeG2LineSearch since one of the G2 component is negative\n",
"Info in : NegativeG2LineSearch Done after 18.4429 ms\n",
"Info in : MnSeedGenerator Negative G2 found - new state: \n",
" Minimum value : -1141.929575\n",
" Edm : 0.8557230398\n",
" Internal parameters:\t[ -1.218687397 1.558064949]\t\n",
" Internal gradient :\t[ 20.04216908 -0.2880868692]\t\n",
" Internal covariance matrix:\n",
"[[ 0.0085029968 0]\n",
" [ 0 0.088389344]]]\n",
"Info in : MnSeedGenerator Initial state \n",
" Minimum value : -1141.929575\n",
" Edm : 0.8557230398\n",
" Internal parameters:\t[ -1.218687397 1.558064949]\t\n",
" Internal gradient :\t[ 20.04216908 -0.2880868692]\t\n",
" Internal covariance matrix:\n",
"[[ 0.0085029968 0]\n",
" [ 0 0.088389344]]]\n",
"Info in : VariableMetricBuilder Start iterating until Edm is < 0.001 with call limit = 1000\n",
"Info in : VariableMetricBuilder 0 - FCN = -1141.929575 Edm = 0.8557230398 NCalls = 39\n",
"Info in : VariableMetricBuilder 1 - FCN = -1143.495594 Edm = 0.01234303762 NCalls = 53\n",
"Info in : VariableMetricBuilder 2 - FCN = -1143.508215 Edm = 0.0003439496695 NCalls = 58\n",
"Info in : MnHesse Done after 6.53555 ms\n",
"Info in : VariableMetricBuilder After Hessian\n",
"Info in : VariableMetricBuilder 3 - FCN = -1143.508215 Edm = 0.003981095725 NCalls = 68\n",
"Info in : VariableMetricBuilder Tolerance not sufficient, continue minimization; Edm 0.0039811 Required 0.001\n",
"Info in : VariableMetricBuilder 4 - FCN = -1143.513055 Edm = 1.940940537e-06 NCalls = 75\n",
"Info in : MnHesse Done after 6.54549 ms\n",
"Info in : VariableMetricBuilder After Hessian\n",
"Info in : VariableMetricBuilder 5 - FCN = -1143.513055 Edm = 2.626340222e-06 NCalls = 85\n",
"Info in : VariableMetricBuilder Stop iterating after 34.8939 ms\n",
"Info in : Minuit2Minimizer::Hesse Using max-calls 1000\n",
"Info in : MnHesse Done after 7.45666 ms\n",
"Info in : Minuit2Minimizer::Hesse Hesse is valid - matrix is accurate\n"
]
}
],
"source": [
"auto oscillationFormula = \"std::pow(std::sin(1.27 * x[2] * x[0] / x[1]), 2)\";\n",
"RooGenericPdf PnmuTone(\"PnmuTone\", \"P(#nu_{#mu} #rightarrow #nu_{e}\", oscillationFormula, {L, E, deltaMSq});\n",
"\n",
"// only E is observable, so create the signal model by integrating out L\n",
"RooAbsPdf *sigModel = PnmuTone.createProjection(L);\n",
"\n",
"// create $ \\int dE' dL' P(E',L' | \\Delta m^2)$.\n",
"// Given RooFit will renormalize the PDF in the range of the observables,\n",
"// the average probability to oscillate in the experiment's acceptance\n",
"// needs to be incorporated into the extended term in the likelihood.\n",
"// Do this by creating a RooAbsReal representing the integral and divide by\n",
"// the area in the E-L plane.\n",
"// The integral should be over \"primed\" observables, so we need\n",
"// an independent copy of PnmuTone not to interfere with the original.\n",
"\n",
"// Independent copy for Integral\n",
"RooRealVar EPrime(\"EPrime\", \"\", 15, 10, 60, \"GeV\");\n",
"RooRealVar LPrime(\"LPrime\", \"\", .800, .600, 1.0, \"km\"); // need these units in formula\n",
"RooGenericPdf PnmuTonePrime(\"PnmuTonePrime\", \"P(#nu_{#mu} #rightarrow #nu_{e}\", oscillationFormula, {LPrime, EPrime, deltaMSq});\n",
"RooAbsReal *intProbToOscInExp = PnmuTonePrime.createIntegral(RooArgSet(EPrime, LPrime));\n",
"\n",
"// Getting the flux is a bit tricky. It is more clear to include a cross section term that is not\n",
"// explicitly referred to in the text, eg.\n",
"// number events in bin = flux * cross-section for nu_e interaction in E bin * average prob nu_mu osc. to nu_e in bin\n",
"// let maxEventsInBin = flux * cross-section for nu_e interaction in E bin\n",
"// maxEventsInBin * 1% chance per bin = 100 events / bin\n",
"// therefore maxEventsInBin = 10,000.\n",
"// for 5 bins, this means maxEventsTot = 50,000\n",
"RooConstVar maxEventsTot(\"maxEventsTot\", \"maximum number of sinal events\", 50000);\n",
"RooConstVar inverseArea(\"inverseArea\", \"1/(#Delta E #Delta L)\",\n",
" 1. / (EPrime.getMax() - EPrime.getMin()) / (LPrime.getMax() - LPrime.getMin()));\n",
"\n",
"// $sigNorm = maxEventsTot \\cdot \\int dE dL \\frac{P_{oscillate\\ in\\ experiment}}{Area} \\cdot {sin}^2(2\\theta)$\n",
"RooProduct sigNorm(\"sigNorm\", \"\", RooArgSet(maxEventsTot, *intProbToOscInExp, inverseArea, sinSq2theta));\n",
"// bkg = 5 bins * 100 events / bin\n",
"RooConstVar bkgNorm(\"bkgNorm\", \"normalization for background\", 500);\n",
"\n",
"// flat background (0th order polynomial, so no arguments for coefficients)\n",
"RooPolynomial bkgEShape(\"bkgEShape\", \"flat bkg shape\", E);\n",
"\n",
"// total model\n",
"RooAddPdf model(\"model\", \"\", RooArgList(*sigModel, bkgEShape), RooArgList(sigNorm, bkgNorm));\n",
"\n",
"// for debugging, check model tree\n",
"// model.printCompactTree();\n",
"// model.graphVizTree(\"model.dot\");\n",
"\n",
"// turn off some messages\n",
"RooMsgService::instance().setStreamStatus(0, false);\n",
"RooMsgService::instance().setStreamStatus(1, false);\n",
"RooMsgService::instance().setStreamStatus(2, false);\n",
"\n",
"// --------------------------------------\n",
"// n events in data to data, simply sum of sig+bkg\n",
"Int_t nEventsData = bkgNorm.getVal() + sigNorm.getVal();\n",
"cout << \"generate toy data with nEvents = \" << nEventsData << endl;\n",
"// adjust random seed to get a toy dataset similar to one in paper.\n",
"// Found by trial and error (3 trials, so not very \"fine tuned\")\n",
"RooRandom::randomGenerator()->SetSeed(3);\n",
"// create a toy dataset\n",
"RooDataSet *data = model.generate(RooArgSet(E), nEventsData);\n",
"\n",
"// --------------------------------------\n",
"// make some plots\n",
"TCanvas *dataCanvas = new TCanvas(\"dataCanvas\");\n",
"dataCanvas->Divide(2, 2);\n",
"\n",
"// plot the PDF\n",
"dataCanvas->cd(1);\n",
"TH1 *hh = PnmuTone.createHistogram(\"hh\", E, Binning(40), YVar(L, Binning(40)), Scaling(false));\n",
"hh->SetLineColor(kBlue);\n",
"hh->SetTitle(\"True Signal Model\");\n",
"hh->Draw(\"surf\");\n",
"\n",
"// plot the data with the best fit\n",
"dataCanvas->cd(2);\n",
"RooPlot *Eframe = E.frame();\n",
"data->plotOn(Eframe);\n",
"model.fitTo(*data, Extended());\n",
"model.plotOn(Eframe);\n",
"model.plotOn(Eframe, Components(*sigModel), LineColor(kRed));\n",
"model.plotOn(Eframe, Components(bkgEShape), LineColor(kGreen));\n",
"model.plotOn(Eframe);\n",
"Eframe->SetTitle(\"toy data with best fit model (and sig+bkg components)\");\n",
"Eframe->Draw();\n",
"\n",
"// plot the likelihood function\n",
"dataCanvas->cd(3);\n",
"std::unique_ptr nll{model.createNLL(*data, Extended(true))};\n",
"RooProfileLL pll(\"pll\", \"\", *nll, RooArgSet(deltaMSq, sinSq2theta));\n",
"// TH1* hhh = nll.createHistogram(\"hhh\",sinSq2theta,Binning(40),YVar(deltaMSq,Binning(40))) ;\n",
"TH1 *hhh = pll.createHistogram(\"hhh\", sinSq2theta, Binning(40), YVar(deltaMSq, Binning(40)), Scaling(false));\n",
"hhh->SetLineColor(kBlue);\n",
"hhh->SetTitle(\"Likelihood Function\");\n",
"hhh->Draw(\"surf\");\n",
"\n",
"dataCanvas->Update();\n",
"\n",
"// --------------------------------------------------------------\n",
"// show use of Feldman-Cousins utility in RooStats\n",
"// set the distribution creator, which encodes the test statistic\n",
"RooArgSet parameters(deltaMSq, sinSq2theta);\n",
"RooWorkspace *w = new RooWorkspace();\n",
"\n",
"ModelConfig modelConfig;\n",
"modelConfig.SetWorkspace(*w);\n",
"modelConfig.SetPdf(model);\n",
"modelConfig.SetParametersOfInterest(parameters);\n",
"\n",
"RooStats::FeldmanCousins fc(*data, modelConfig);\n",
"fc.SetTestSize(.1); // set size of test\n",
"fc.UseAdaptiveSampling(true);\n",
"fc.SetNBins(10); // number of points to test per parameter\n",
"\n",
"// use the Feldman-Cousins tool\n",
"ConfInterval *interval = 0;\n",
"if (doFeldmanCousins)\n",
" interval = fc.GetInterval();\n",
"\n",
"// ---------------------------------------------------------\n",
"// show use of ProfileLikeihoodCalculator utility in RooStats\n",
"RooStats::ProfileLikelihoodCalculator plc(*data, modelConfig);\n",
"plc.SetTestSize(.1);\n",
"\n",
"ConfInterval *plcInterval = plc.GetInterval();\n",
"\n",
"// --------------------------------------------\n",
"// show use of MCMCCalculator utility in RooStats\n",
"MCMCInterval *mcInt = NULL;\n",
"\n",
"if (doMCMC) {\n",
" // turn some messages back on\n",
" RooMsgService::instance().setStreamStatus(0, true);\n",
" RooMsgService::instance().setStreamStatus(1, true);\n",
"\n",
" TStopwatch mcmcWatch;\n",
" mcmcWatch.Start();\n",
"\n",
" RooArgList axisList(deltaMSq, sinSq2theta);\n",
" MCMCCalculator mc(*data, modelConfig);\n",
" mc.SetNumIters(5000);\n",
" mc.SetNumBurnInSteps(100);\n",
" mc.SetUseKeys(true);\n",
" mc.SetTestSize(.1);\n",
" mc.SetAxes(axisList); // set which is x and y axis in posterior histogram\n",
" // mc.SetNumBins(50);\n",
" mcInt = (MCMCInterval *)mc.GetInterval();\n",
"\n",
" mcmcWatch.Stop();\n",
" mcmcWatch.Print();\n",
"}\n",
"// -------------------------------\n",
"// make plot of resulting interval\n",
"\n",
"dataCanvas->cd(4);\n",
"\n",
"// first plot a small dot for every point tested\n",
"if (doFeldmanCousins) {\n",
" RooDataHist *parameterScan = (RooDataHist *)fc.GetPointsToScan();\n",
" TH2F *hist = parameterScan->createHistogram(deltaMSq,sinSq2theta, 30, 30);\n",
" // hist->Draw();\n",
" TH2F *forContour = (TH2F *)hist->Clone();\n",
"\n",
" // now loop through the points and put a marker if it's in the interval\n",
" RooArgSet *tmpPoint;\n",
" // loop over points to test\n",
" for (Int_t i = 0; i < parameterScan->numEntries(); ++i) {\n",
" // get a parameter point from the list of points to test.\n",
" tmpPoint = (RooArgSet *)parameterScan->get(i)->clone(\"temp\");\n",
"\n",
" if (interval) {\n",
" if (interval->IsInInterval(*tmpPoint)) {\n",
" forContour->SetBinContent(\n",
" hist->FindBin(tmpPoint->getRealValue(\"sinSq2theta\"), tmpPoint->getRealValue(\"deltaMSq\")), 1);\n",
" } else {\n",
" forContour->SetBinContent(\n",
" hist->FindBin(tmpPoint->getRealValue(\"sinSq2theta\"), tmpPoint->getRealValue(\"deltaMSq\")), 0);\n",
" }\n",
" }\n",
"\n",
" delete tmpPoint;\n",
" }\n",
"\n",
" if (interval) {\n",
" Double_t level = 0.5;\n",
" forContour->SetContour(1, &level);\n",
" forContour->SetLineWidth(2);\n",
" forContour->SetLineColor(kRed);\n",
" forContour->Draw(\"cont2,same\");\n",
" }\n",
"}\n",
"\n",
"MCMCIntervalPlot *mcPlot = NULL;\n",
"if (mcInt) {\n",
" cout << \"MCMC actual confidence level: \" << mcInt->GetActualConfidenceLevel() << endl;\n",
" mcPlot = new MCMCIntervalPlot(*mcInt);\n",
" mcPlot->SetLineColor(kMagenta);\n",
" mcPlot->Draw();\n",
"}\n",
"dataCanvas->Update();\n",
"\n",
"LikelihoodIntervalPlot plotInt((LikelihoodInterval *)plcInterval);\n",
"plotInt.SetTitle(\"90% Confidence Intervals\");\n",
"if (mcInt)\n",
" plotInt.Draw(\"same\");\n",
"else\n",
" plotInt.Draw();\n",
"dataCanvas->Update();\n",
"\n",
"/// print timing info\n",
"t.Stop();\n",
"t.Print();"
]
},
{
"cell_type": "markdown",
"id": "83ad37df",
"metadata": {},
"source": [
"Draw all canvases "
]
},
{
"cell_type": "code",
"execution_count": 7,
"id": "ba57fad4",
"metadata": {
"collapsed": false,
"execution": {
"iopub.execute_input": "2026-05-19T20:37:08.084854Z",
"iopub.status.busy": "2026-05-19T20:37:08.084720Z",
"iopub.status.idle": "2026-05-19T20:37:08.294073Z",
"shell.execute_reply": "2026-05-19T20:37:08.293540Z"
}
},
"outputs": [
{
"data": {
"text/html": [
"\n",
"\n",
"\n",
"
\n",
"\n",
"\n",
"\n"
],
"text/plain": [
""
]
},
"metadata": {},
"output_type": "display_data"
}
],
"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
}