{
 "cells": [
  {
   "cell_type": "markdown",
   "id": "217a2358",
   "metadata": {},
   "source": [
    "# rs401c_FeldmanCousins\n",
    "Produces an interval on the mean signal in a number counting experiment with known background using the\n",
    "Feldman-Cousins technique.\n",
    "\n",
    "Using the RooStats FeldmanCousins tool with 200 bins\n",
    "it takes 1 min and the interval is [0.2625, 10.6125]\n",
    "with a step size of 0.075.\n",
    "The interval in Feldman & Cousins's original paper is [.29, 10.81] Phys.Rev.D57:3873-3889,1998.\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:35 PM.</small></i>"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "975c7554",
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "%%cpp -d\n",
    "#include \"RooGlobalFunc.h\"\n",
    "#include \"RooStats/ConfInterval.h\"\n",
    "#include \"RooStats/PointSetInterval.h\"\n",
    "#include \"RooStats/ConfidenceBelt.h\"\n",
    "#include \"RooStats/FeldmanCousins.h\"\n",
    "#include \"RooStats/ModelConfig.h\"\n",
    "\n",
    "#include \"RooWorkspace.h\"\n",
    "#include \"RooDataSet.h\"\n",
    "#include \"RooRealVar.h\"\n",
    "#include \"RooConstVar.h\"\n",
    "#include \"RooAddition.h\"\n",
    "\n",
    "#include \"RooDataHist.h\"\n",
    "\n",
    "#include \"RooPoisson.h\"\n",
    "#include \"RooPlot.h\"\n",
    "\n",
    "#include \"TCanvas.h\"\n",
    "#include \"TTree.h\"\n",
    "#include \"TH1F.h\"\n",
    "#include \"TMarker.h\"\n",
    "#include \"TStopwatch.h\"\n",
    "\n",
    "#include <iostream>\n",
    "\n",
    "using namespace RooFit;\n",
    "using namespace RooStats;"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "85c4d6c3",
   "metadata": {},
   "source": [
    "to time the macro... about 30 s"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "6361259a",
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "TStopwatch t;\n",
    "t.Start();"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "48b15dd9",
   "metadata": {},
   "source": [
    "make a simple model"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "4ffc53a0",
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "RooRealVar x(\"x\", \"\", 1, 0, 50);\n",
    "RooRealVar mu(\"mu\", \"\", 2.5, 0, 15); // with a limit on mu>=0\n",
    "RooConstVar b(\"b\", \"\", 3.);\n",
    "RooAddition mean(\"mean\", \"\", RooArgList(mu, b));\n",
    "RooPoisson pois(\"pois\", \"\", x, mean);\n",
    "RooArgSet parameters(mu);"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "38c17ff6",
   "metadata": {},
   "source": [
    "create a toy dataset"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "0c6d9d0f",
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "std::unique_ptr<RooDataSet> data{pois.generate({x}, 1)};\n",
    "data->Print(\"v\");\n",
    "\n",
    "TCanvas *dataCanvas = new TCanvas(\"dataCanvas\");\n",
    "RooPlot *frame = x.frame();\n",
    "data->plotOn(frame);\n",
    "frame->Draw();\n",
    "dataCanvas->Update();\n",
    "\n",
    "RooWorkspace *w = new RooWorkspace();\n",
    "ModelConfig modelConfig(\"poissonProblem\", w);\n",
    "modelConfig.SetPdf(pois);\n",
    "modelConfig.SetParametersOfInterest(parameters);\n",
    "modelConfig.SetObservables(RooArgSet(x));\n",
    "w->Print();"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "983a70b9",
   "metadata": {},
   "source": [
    "show use of Feldman-Cousins"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "e4022b43",
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "RooStats::FeldmanCousins fc(*data, modelConfig);\n",
    "fc.SetTestSize(.05); // set size of test\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"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "eda85fed",
   "metadata": {},
   "source": [
    "use the Feldman-Cousins tool"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "41d24d2d",
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "PointSetInterval *interval = (PointSetInterval *)fc.GetInterval();"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "4eac7c89",
   "metadata": {},
   "source": [
    "make a canvas for plots"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "e1acb104",
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "TCanvas *intervalCanvas = new TCanvas(\"intervalCanvas\");\n",
    "\n",
    "std::cout << \"is this point in the interval? \" << interval->IsInInterval(parameters) << std::endl;\n",
    "\n",
    "std::cout << \"interval is [\" << interval->LowerLimit(mu) << \", \" << interval->UpperLimit(mu) << \"]\" << endl;"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "55ed23f1",
   "metadata": {},
   "source": [
    "using 200 bins it takes 1 min and the answer is\n",
    "interval is [0.2625, 10.6125] with a step size of .075\n",
    "The interval in Feldman & Cousins's original paper is [.29, 10.81]\n",
    "Phys.Rev.D57:3873-3889,1998."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "e5b42460",
   "metadata": {},
   "source": [
    "No dedicated plotting class yet, so do it by hand:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "3a8a3882",
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "RooDataHist *parameterScan = (RooDataHist *)fc.GetPointsToScan();\n",
    "TH1F *hist = (TH1F *)parameterScan->createHistogram(\"mu\", Binning(30));\n",
    "hist->Draw();\n",
    "\n",
    "RooArgSet *tmpPoint;"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "0b1df0d9",
   "metadata": {},
   "source": [
    "loop over points to test"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "98679c06",
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "for (Int_t i = 0; i < parameterScan->numEntries(); ++i) {\n",
    "   //    cout << \"on parameter point \" << i << \" out of \" << parameterScan->numEntries() << endl;\n",
    "   // get a parameter point from the list of points to test.\n",
    "   tmpPoint = (RooArgSet *)parameterScan->get(i)->clone(\"temp\");\n",
    "\n",
    "   TMarker *mark = new TMarker(tmpPoint->getRealValue(\"mu\"), 1, 25);\n",
    "   if (interval->IsInInterval(*tmpPoint))\n",
    "      mark->SetMarkerColor(kBlue);\n",
    "   else\n",
    "      mark->SetMarkerColor(kRed);\n",
    "\n",
    "   mark->Draw(\"s\");\n",
    "   // delete tmpPoint;\n",
    "   //    delete mark;\n",
    "}\n",
    "t.Stop();\n",
    "t.Print();"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "ae530a0e",
   "metadata": {},
   "source": [
    "Draw all canvases "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "350bee16",
   "metadata": {
    "collapsed": false
   },
   "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
}