1/exercise_1b.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# File: exercise_1b.py\n",
    "## Analyze the following Run I CMS $H \\rightarrow 4 \\ell$ data\n",
    "\n",
    " - Created: 18-Dec-2015 CMSDAS 2016, LPC Fermilab HBP\n",
    " - Updated: 10-Jan-2019 CMSDAS 2019, LPC Fermilab JMD\n",
    " - Updated: 03-Jan-2022 CMSDAS 2019, LPC Fermilab DRY\n",
    " - Updated: 07-Jan-2024 CMSDAS 2024, LPC Fermilab HBP\n",
    "\n",
    "```\n",
    "   N        = 25\n",
    "   b_hat_zz =  6.8 +/- 0.3\n",
    "   b_hat_zx =  2.6 +/- 0.4\n",
    "   s_hat    = 17.3 +/- 1.3 (mH = 125 GeV)\n",
    "            = 19.6 +/- 1.3 (mH = 126 GeV)\n",
    "```\n",
    " which pertains to data in the range $121.5 \\le m_{4\\ell} \\le 130.5\\text{ GeV}$\n",
    " for 7 and 8 TeV data. The main backgrounds are $p p \\rightarrow ZZ \\rightarrow 4l$ and $p p \\rightarrow Z + X \\rightarrow 4l$, where $X$ is typically one or ore jets. \n",
    " \n",
    "## Introduction\n",
    "\n",
    "See **exercise_1a**.\n",
    "\n",
    "## Statistical Model\n",
    "\n",
    "As before, we'll model event counts $X$ with a $\\text{Poisson}(X, n)$ distribution, where the parameter $n$ is the mean count. In this version of **exercise_1**, we'll keep things simple and model the signal and background estimates with Gaussians though this is potentially problematic. Using Gaussians here is potentially problematic because signal and background estimates are constrained to be $\\ge 0$, while the domain of a $\\text{Gaussian}(x, \\mu, \\sigma)$ is $-\\infty < x < +\\infty$. \n",
    "What saves us is that the estimated signal and background uncertainties are small and, therefore, the supports of the Gaussians are far enough from zero to avoid problems with negative estimates. \n",
    "\n",
    "The statistical model in this version of the single-count analysis is given by\n",
    "\n",
    "\n",
    "\\begin{align}\n",
    "    p(X, Y| \\mu, \\nu) & = \\text{Poisson}(X, \\mu  s + b_{ZZ} + b_{ZX}) \\nonumber\\\\\n",
    "                    & \\times \\text{Gaussian}(\\hat{s}, s, \\delta s) \\nonumber\\\\\n",
    "                    & \\times \\text{Gaussian}(\\hat{b}_{ZZ}, b_{ZZ}, \\delta b_{ZZ}) \\nonumber\\\\\n",
    "                    & \\times \\text{Gaussian}(\\hat{b}_{ZX}, b_{ZX}, \\delta b_{ZX}) ,\n",
    "\\end{align}"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {
    "scrolled": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Welcome to JupyROOT 6.28/06\n"
     ]
    }
   ],
   "source": [
    "import os, sys\n",
    "import ROOT"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {
    "scrolled": false
   },
   "outputs": [],
   "source": [
    "def createWorkspace(wsname, wsfilename):\n",
    "    # The most convenient way to use RooFit/RooStats is to \n",
    "    # make a workspace so that we can use its factory method\n",
    "    wspace = ROOT.RooWorkspace(wsname)\n",
    "\n",
    "    #-----------------------------------------------------\n",
    "    # Create parameters\n",
    "    #\n",
    "    # Use the factory method of the RooWorkspace to create\n",
    "    # parameters\n",
    "    #\n",
    "    # syntax:\n",
    "    #        <name>[value, min-value, max-value]\n",
    "    #-----------------------------------------------------\n",
    "    # observations\n",
    "    params = [('N',       25,     0,  50),\n",
    "    # ZZ background estimate        \n",
    "              ('b_hat_zz', 6.8,   0,  15),\n",
    "              ('db_zz',    0.3,   0,   5),\n",
    "    # Z+X background estimate\n",
    "              ('b_hat_zx', 2.6,   0,  15),\n",
    "              ('db_zx',    0.3,   0,   5),\n",
    "    # signal estimate (mH=125 GeV)       \n",
    "              ('s_hat',   17.3,   0,  25),\n",
    "              ('ds',       1.3,   0,   5),\n",
    "    # nuisance parameters\n",
    "              ('b_zx',    2.6,    0,  10),\n",
    "              ('b_zz',    6.8,    0,  15),\n",
    "              ('s',      17.3,    0,  25),\n",
    "    # parameter of interest\n",
    "              ('mu',      1.0,    0,   4)]\n",
    "\n",
    "    for t in params:\n",
    "        cmd = '{}[{}, {}, {}]'.format(t[0], t[1], t[2], t[3]) # Could also do .format(*t)?\n",
    "        wspace.factory(cmd)        \n",
    "    wspace.var('mu').SetTitle('#mu')\n",
    "\n",
    "    # fix all background and signal parameters\n",
    "    for t in params[1:-4]:\n",
    "        name = t[0]\n",
    "        print('=> make {:8s} = {:5.1f} constant'.format(name, wspace.var(name).getVal()))\n",
    "        wspace.var(name).setConstant()\n",
    "\n",
    "    #-----------------------------------------------------\n",
    "    # Create expression for the mean count n of the\n",
    "    # Poisson distribution\n",
    "    #\n",
    "    # syntax:\n",
    "    #        expr::<name>(\"expression\", var1, var2, ...)\n",
    "    #-----------------------------------------------------\n",
    "    cmd  = 'expr::n(\"mu*s + b_zz + b_zx\", mu, s, b_zz, b_zx)'\n",
    "    wspace.factory(cmd)\n",
    "\n",
    "    #-----------------------------------------------------\n",
    "    # Create pdfs\n",
    "    #\n",
    "    # syntax:\n",
    "    #        pdf_name::<name>(var1, var2, ...)\n",
    "    #\n",
    "    # where the \"Roo\" prefix is dropped in pdf_name, e.g.\n",
    "    #-----------------------------------------------------\n",
    "    pdfs = [('Poisson', 'pN',    '(N, n)'),\n",
    "            # Truncated Gaussian constraints (in general, bad! See comments above.)\n",
    "            ('Gaussian','pB_zz', '(b_zz, b_hat_zz, db_zz)'),\n",
    "            ('Gaussian','pB_zx', '(b_zx, b_hat_zx, db_zx)'),\n",
    "            ('Gaussian','pS', '(s, s_hat, ds)'), \n",
    "           ]\n",
    "    \n",
    "    prodpdf = ''\n",
    "    for pdfargs in pdfs:\n",
    "        wspace.factory('{}::{}{}'.format(pdfargs[0], pdfargs[1], pdfargs[2]))\n",
    "        name = pdfargs[1]\n",
    "        prodpdf += \"{}, \".format(name)\n",
    "    prodpdf = prodpdf[:-2] # remove last \", \"\n",
    "    \n",
    "    # multiply the pdfs together. use upper case PROD to\n",
    "    # do this\n",
    "    wspace.factory('PROD::model({})'.format(prodpdf))\n",
    "\n",
    "    # create a prior, since one is needed in Bayesian\n",
    "    # calculations\n",
    "    wspace.factory('Uniform::prior({mu, s, b_zz, b_zx})')\n",
    "\n",
    "    #-----------------------------------------------------\n",
    "    # Define a few useful sets. We'll treat the signal\n",
    "    # and background estimates as constants rather than\n",
    "    # auxiliary\n",
    "    #-----------------------------------------------------\n",
    "    sets = [('obs',  'N'),           # observations\n",
    "            ('poi',  'mu'),          # parameter of interest\n",
    "            ('nuis', 's,b_zz,b_zx')] # nuisance parameters (leave no spaces)\n",
    "    for t in sets:\n",
    "        name, parlist = t\n",
    "        wspace.defineSet(name, parlist)\n",
    "    \n",
    "    #-----------------------------------------------------        \n",
    "    # create a dataset\n",
    "    #-----------------------------------------------------    \n",
    "    data = ROOT.RooDataSet('data', 'data', wspace.set('obs'))\n",
    "    data.add(wspace.set('obs'))\n",
    "    \n",
    "    # import dataset into workspace\n",
    "    wspace.Import(data)\n",
    "        \n",
    "    #-----------------------------------------------------\n",
    "    # Create model configuration. This is needed for the\n",
    "    # statistical analyses\n",
    "    #-----------------------------------------------------\n",
    "    cfg = ROOT.RooStats.ModelConfig('cfg')\n",
    "    cfg.SetWorkspace(wspace)\n",
    "    cfg.SetPdf(wspace.pdf('model'))\n",
    "    cfg.SetPriorPdf(wspace.pdf('prior'))\n",
    "    cfg.SetParametersOfInterest(wspace.set('poi'))\n",
    "    cfg.SetNuisanceParameters(wspace.set('nuis'))\n",
    "\n",
    "    # import model configuration into workspace\n",
    "    wspace.Import(cfg)\n",
    "\n",
    "    wspace.Print()\n",
    "    \n",
    "    # write out workspace\n",
    "    wspace.writeToFile(wsfilename)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "createWorkspace('CMSDAS', 'single_count_1.2.root')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {
    "scrolled": false
   },
   "outputs": [],
   "source": [
    "def analyzeWorkspace(wsname, wsfilename):\n",
    "\n",
    "    # Open workspace file\n",
    "    wsfile = ROOT.TFile.Open(wsfilename)\n",
    "\n",
    "    # Get workspace\n",
    "    wspace = wsfile.Get(wsname) \n",
    "\n",
    "    # Get data\n",
    "    data = wspace.data('data')\n",
    "\n",
    "    # Get model configuration    \n",
    "    cfg  = wspace.obj('cfg')\n",
    "\n",
    "    #-----------------------------------------------------    \n",
    "    # Fit model to data\n",
    "    #-----------------------------------------------------\n",
    "    results = wspace.pdf('model').fitTo(data, ROOT.RooFit.Save())\n",
    "    results.Print()\n",
    "    \n",
    "    #-----------------------------------------------------    \n",
    "    # Compute interval based on profile likelihood\n",
    "    #-----------------------------------------------------\n",
    "    # suppress some (apparently) innocuous warnings\n",
    "    msgservice = ROOT.RooMsgService.instance()\n",
    "    msgservice.setGlobalKillBelow(ROOT.RooFit.FATAL)\n",
    "        \n",
    "    print('compute interval using profile likelihood')\n",
    "    plc = ROOT.RooStats.ProfileLikelihoodCalculator(data, cfg)\n",
    "    CL  = 0.683\n",
    "    plc.SetConfidenceLevel(CL)\n",
    "    plcInterval= plc.GetInterval()\n",
    "    lowerLimit = plcInterval.LowerLimit(wspace.var('mu'))\n",
    "    upperLimit = plcInterval.UpperLimit(wspace.var('mu'))\n",
    "\n",
    "    print('\\tPL {:4.1f}% CL interval = [{:5.2f}, {:5.2f}]'.format(100*CL, lowerLimit, upperLimit))\n",
    "\n",
    "    plcplot = ROOT.RooStats.LikelihoodIntervalPlot(plcInterval)      \n",
    "    plccanvas = ROOT.TCanvas('fig_PL_1.2', 'PL', 800, 400)\n",
    "    plccanvas.Divide(2, 1)\n",
    "    plccanvas.cd(1)\n",
    "    plcplot.SetRange(0,4)\n",
    "    plcplot.SetMaximum(3)\n",
    "    plcplot.Draw()\n",
    "    \n",
    "    #----------------------------------------------------- \n",
    "    # compute an 95% upper limit on mu by\n",
    "    # computing a 90% central interval and\n",
    "    # ignoring the lower limit\n",
    "    #----------------------------------------------------- \n",
    "    CL = 0.90\n",
    "    plc.SetConfidenceLevel(CL)\n",
    "    plcInterval = plc.GetInterval()\n",
    "    upperLimit  = plcInterval.UpperLimit(wspace.var('mu'))\n",
    "\n",
    "    CL = 0.95\n",
    "    print('\\tPL {:4.1f}% upper limit = {:5.2f}\\n'.format(100*CL, upperLimit))\n",
    "      \n",
    "    plccanvas.cd(2)\n",
    "    plcplot2 = ROOT.RooStats.LikelihoodIntervalPlot(plcInterval)\n",
    "    plcplot2.SetRange(0,4)\n",
    "    plcplot2.SetMaximum(3)\n",
    "    plcplot2.Draw()\n",
    "    plccanvas.Update()\n",
    "    \n",
    "    #-----------------------------------------------------    \n",
    "    # Compute interval based on Bayesian calculator\n",
    "    #-----------------------------------------------------\n",
    "    print('compute interval using Bayesian calculator')\n",
    "    bc = ROOT.RooStats.BayesianCalculator(data, cfg)\n",
    "    CL  = 0.683\n",
    "    bc.SetConfidenceLevel(CL)\n",
    "    bcInterval = bc.GetInterval()\n",
    "    lowerLimit = bcInterval.LowerLimit()\n",
    "    upperLimit = bcInterval.UpperLimit()\n",
    "\n",
    "    print('\\tBayes {:4.1f}% CL interval = [{:5.2f}, {:5.2f}]'.format(100*CL, lowerLimit, upperLimit))\n",
    "\n",
    "    # calculate posterior density at 50 points\n",
    "    print(\"\\t\\t...be patient...!\")\n",
    "    bc.SetScanOfPosterior(50)\n",
    "    bcplot = bc.GetPosteriorPlot()\n",
    "    bccanvas = ROOT.TCanvas('fig_Bayes_1.2', 'Bayes', 800, 400)\n",
    "    bccanvas.Divide(2, 1)\n",
    "    bccanvas.cd(1)\n",
    "    bcplot.Draw()\n",
    "    bccanvas.Update()\n",
    "\n",
    "    # compute an 95% upper limit on mu\n",
    "    CL  = 0.950\n",
    "    bc.SetConfidenceLevel(CL)\n",
    "    # 0   = upper limit\n",
    "    # 0.5 = central limits (default)\n",
    "    # 1   = lower limit\n",
    "    bc.SetLeftSideTailFraction(0)\n",
    "    bcInterval = bc.GetInterval()\n",
    "    upperLimit = bcInterval.UpperLimit()\n",
    "\n",
    "    print('\\tBayes {:4.1f}% upper limit = {:5.2f}\\n'.format(100*CL, upperLimit))\n",
    "\n",
    "    # calculate posterior density at 50 points\n",
    "    bc.SetScanOfPosterior(50)\n",
    "    bcplot2 = bc.GetPosteriorPlot()\n",
    "    bccanvas.cd(2)\n",
    "    bcplot2.Draw()\n",
    "    bccanvas.Update()\n",
    "\n",
    "    # save canvases\n",
    "    plccanvas.Draw()\n",
    "    bccanvas.Draw()\n",
    "    plccanvas.SaveAs('.png')\n",
    "    bccanvas.SaveAs('.png')\n",
    "    return plccanvas, bccanvas"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "plccanvas, bccanvas = analyzeWorkspace('CMSDAS', 'single_count_1.2.root')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": false
   },
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "@webio": {
   "lastCommId": null,
   "lastKernelId": null
  },
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.9.12"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}