examples/combine_limits.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Reading combine limits\n",
    "\n",
    "The CMS Higgs Physics Analysis Group supports a widely used tool called [combine](https://cms-hcomb.gitbooks.io/combine/content/). `hepdata_lib` has the capability to directly import the limits from the ntuple created, see e.g. [Asymptotic Frequentists Limits](https://cms-hcomb.gitbooks.io/combine/content/part3/commonstatsmethods.html#asymptotic-frequentist-limits).\n",
    "\n",
    "As explained in the [Getting started notebook](Getting_started.ipynb), a `Submission` needs to exist or be created. Here, we'll just create one without any additional information:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Welcome to JupyROOT 6.28/06\n"
     ]
    }
   ],
   "source": [
    "from hepdata_lib import Submission\n",
    "submission = Submission()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The limit plot will again be a `Table`, in the example Figure 6a from page 14 of the publication. Let's add all this, some more details as well as the actual plot (for thumbnail creation) to the `Table`:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [],
   "source": [
    "from hepdata_lib import Table\n",
    "table = Table(\"Figure 6a\")\n",
    "table.description = \"Exclusion limits on the product of the production cross section and the branching fraction for a new spin-2 resonance decaying to WW, as a function of the resonance mass hypothesis.\"\n",
    "\n",
    "table.location = \"Data from Figure 6a, located on page 14.\"\n",
    "table.keywords[\"observables\"] = [\"SIG\"]\n",
    "table.keywords[\"reactions\"] = [\"P P --> GRAVITON --> W+ W-\"]\n",
    "\n",
    "table.add_image(\"example_inputs/Figure_006-a.pdf\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now we can read in the `combine` ntuple, which in our case has the name `WWmerged.higgsCombineTest.Asymptotic.root`, and contains the upper limits on the production cross section of a resonance as a function of the mass of the resonance (\"bulk graviton\"). The file is read in with the `RootFileReader` using `read_limit_tree()`:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[[1.00000000e+03 1.82468891e-02 2.45227814e-02 3.55224609e-02\n",
      "  5.20889647e-02 7.50288591e-02 2.16936472e-02]\n",
      " [1.05000000e+03 1.50303841e-02 2.02972870e-02 2.89306641e-02\n",
      "  4.24229726e-02 6.07501119e-02 2.36925946e-02]\n",
      " [1.10000000e+03 1.26993656e-02 1.71843916e-02 2.45361328e-02\n",
      "  3.55877690e-02 5.09835742e-02 2.60327992e-02]\n",
      " [1.15000000e+03 1.08987093e-02 1.47477984e-02 2.10571289e-02\n",
      "  3.07096094e-02 4.38561961e-02 2.58909682e-02]\n",
      " [1.20000000e+03 9.50872898e-03 1.28669199e-02 1.83715820e-02\n",
      "  2.67930217e-02 3.82629447e-02 2.56530152e-02]\n",
      " [1.25000000e+03 8.37147236e-03 1.13280192e-02 1.61743164e-02\n",
      "  2.37174835e-02 3.37647162e-02 2.82829272e-02]\n",
      " [1.30000000e+03 7.48693943e-03 1.01310965e-02 1.44653320e-02\n",
      "  2.10961662e-02 3.01272999e-02 2.87586050e-02]\n",
      " [1.35000000e+03 6.72876835e-03 9.10516270e-03 1.30004883e-02\n",
      "  1.89598463e-02 2.70764343e-02 2.79213802e-02]\n",
      " [1.40000000e+03 6.04957342e-03 8.18609633e-03 1.16882324e-02\n",
      "  1.70460586e-02 2.43433677e-02 2.63606315e-02]\n",
      " [1.45000000e+03 5.48094511e-03 7.41664600e-03 1.05895996e-02\n",
      "  1.55282393e-02 2.21063327e-02 2.37855563e-02]\n",
      " [1.50000000e+03 5.00708818e-03 6.77543785e-03 9.67407227e-03\n",
      "  1.41857406e-02 2.01951228e-02 1.95946506e-02]\n",
      " [1.55000000e+03 4.57906723e-03 6.19214773e-03 8.88061523e-03\n",
      "  1.30222412e-02 1.85387395e-02 1.24887261e-02]\n",
      " [1.60000000e+03 4.20141220e-03 5.69687225e-03 8.14819336e-03\n",
      "  1.19482419e-02 1.71099994e-02 8.62416925e-03]\n",
      " [1.65000000e+03 3.86983156e-03 5.23653720e-03 7.47680664e-03\n",
      "  1.10233482e-02 1.57358106e-02 7.97497603e-03]\n",
      " [1.70000000e+03 3.55845690e-03 4.83500492e-03 6.92749023e-03\n",
      "  1.02134701e-02 1.45797106e-02 7.67257264e-03]\n",
      " [1.75000000e+03 3.27628851e-03 4.45161201e-03 6.37817383e-03\n",
      "  9.45443846e-03 1.35316607e-02 7.03882435e-03]\n",
      " [1.80000000e+03 3.04898620e-03 4.12021624e-03 5.93566895e-03\n",
      "  8.75118840e-03 1.25649497e-02 6.02582596e-03]\n",
      " [1.85000000e+03 2.82952189e-03 3.82364541e-03 5.50842285e-03\n",
      "  8.12128279e-03 1.16605321e-02 4.99097279e-03]\n",
      " [1.90000000e+03 2.62132287e-03 3.55680892e-03 5.14221191e-03\n",
      "  7.58136390e-03 1.08853169e-02 4.32900421e-03]\n",
      " [1.95000000e+03 2.43464112e-03 3.32179712e-03 4.77600098e-03\n",
      "  7.07951980e-03 1.01907123e-02 3.98417229e-03]\n",
      " [2.00000000e+03 2.27907300e-03 3.09241842e-03 4.47082520e-03\n",
      "  6.62715407e-03 9.53954831e-03 3.69187159e-03]\n",
      " [2.05000000e+03 2.13906169e-03 2.90244073e-03 4.19616699e-03\n",
      "  6.22002548e-03 8.95350147e-03 3.31831580e-03]\n",
      " [2.10000000e+03 1.99916959e-03 2.72386242e-03 3.95202637e-03\n",
      "  5.85813308e-03 8.43257084e-03 2.90806403e-03]\n",
      " [2.15000000e+03 1.88291073e-03 2.56727636e-03 3.70788574e-03\n",
      "  5.52579993e-03 7.97376037e-03 2.59142203e-03]\n",
      " [2.20000000e+03 1.76760554e-03 2.42184661e-03 3.49426270e-03\n",
      "  5.20744128e-03 7.55669596e-03 2.45318057e-03]\n",
      " [2.25000000e+03 1.67497993e-03 2.28215521e-03 3.31115723e-03\n",
      "  4.93456237e-03 7.16071203e-03 2.35423506e-03]\n",
      " [2.30000000e+03 1.57013535e-03 2.14826828e-03 3.12805176e-03\n",
      "  4.66168299e-03 6.76472764e-03 2.26915726e-03]\n",
      " [2.35000000e+03 1.47822499e-03 2.02251598e-03 2.94494629e-03\n",
      "  4.41228133e-03 6.41753245e-03 2.16532549e-03]\n",
      " [2.40000000e+03 1.39780343e-03 1.91248278e-03 2.78472900e-03\n",
      "  4.17223526e-03 6.06839266e-03 2.03717706e-03]\n",
      " [2.45000000e+03 1.32121146e-03 1.80768920e-03 2.63214111e-03\n",
      "  3.94361932e-03 5.76759921e-03 1.89680624e-03]\n",
      " [2.50000000e+03 1.24740601e-03 1.70543790e-03 2.49481201e-03\n",
      "  3.75775434e-03 5.47763519e-03 1.76238368e-03]\n",
      " [2.55000000e+03 1.17710233e-03 1.61612616e-03 2.37274170e-03\n",
      "  3.57388868e-03 5.20961639e-03 1.65494961e-03]\n",
      " [2.60000000e+03 1.11654401e-03 1.53298141e-03 2.25067139e-03\n",
      "  3.39002302e-03 4.96858172e-03 1.57756203e-03]\n",
      " [2.65000000e+03 1.06355548e-03 1.46022975e-03 2.14385986e-03\n",
      "  3.22914054e-03 4.73278435e-03 1.52508232e-03]\n",
      " [2.70000000e+03 1.01056695e-03 1.38747809e-03 2.03704834e-03\n",
      "  3.06825805e-03 4.52140998e-03 1.54582785e-03]\n",
      " [2.75000000e+03 9.57578421e-04 1.31472643e-03 1.93023682e-03\n",
      "  2.92276358e-03 4.31558909e-03 1.62568427e-03]\n",
      " [2.80000000e+03 8.97794962e-04 1.24327769e-03 1.83868408e-03\n",
      "  2.78413459e-03 4.11089743e-03 1.68444329e-03]\n",
      " [2.85000000e+03 8.53091478e-04 1.18137174e-03 1.74713135e-03\n",
      "  2.65943375e-03 3.93427955e-03 1.70904340e-03]\n",
      " [2.90000000e+03 8.08387995e-04 1.11946580e-03 1.65557861e-03\n",
      "  2.53327307e-03 3.75451311e-03 1.68491204e-03]\n",
      " [2.95000000e+03 7.71135092e-04 1.06787751e-03 1.57928467e-03\n",
      "  2.41653249e-03 3.58149409e-03 1.63289677e-03]\n",
      " [3.00000000e+03 7.33882189e-04 1.01628923e-03 1.50299072e-03\n",
      "  2.29979167e-03 3.42621398e-03 1.56086459e-03]\n",
      " [3.05000000e+03 6.96629286e-04 9.64700943e-04 1.42669678e-03\n",
      "  2.20579864e-03 3.28041799e-03 1.46800676e-03]\n",
      " [3.10000000e+03 6.56157732e-04 9.16678691e-04 1.36566162e-03\n",
      "  2.10054568e-03 3.15071130e-03 1.35949652e-03]\n",
      " [3.15000000e+03 6.26832247e-04 8.75709811e-04 1.30462646e-03\n",
      "  2.00666650e-03 3.00989719e-03 1.24343458e-03]\n",
      " [3.20000000e+03 5.92648983e-04 8.31666868e-04 1.24359131e-03\n",
      "  1.92765845e-03 2.89807026e-03 1.13871649e-03]\n",
      " [3.25000000e+03 5.63561916e-04 7.86012970e-04 1.18255615e-03\n",
      "  1.83304946e-03 2.78337952e-03 1.06529953e-03]\n",
      " [3.30000000e+03 5.41746616e-04 7.55586661e-04 1.13677979e-03\n",
      "  1.76209270e-03 2.64915614e-03 1.00478655e-03]\n",
      " [3.35000000e+03 5.12659550e-04 7.15018250e-04 1.07574463e-03\n",
      "  1.68463553e-03 2.56467308e-03 9.42639734e-04]\n",
      " [3.40000000e+03 4.90844250e-04 6.84591942e-04 1.02996826e-03\n",
      "  1.61294895e-03 2.45553814e-03 8.82151778e-04]\n",
      " [3.45000000e+03 4.69028950e-04 6.54165633e-04 9.84191895e-04\n",
      "  1.55695446e-03 2.37582577e-03 8.13336504e-04]\n",
      " [3.50000000e+03 4.39584255e-04 6.24335371e-04 9.53674316e-04\n",
      "  1.49347121e-03 2.27364642e-03 7.49498249e-04]\n",
      " [3.55000000e+03 4.18484211e-04 5.94367273e-04 9.07897949e-04\n",
      "  1.43626030e-03 2.19165324e-03 6.95054358e-04]\n",
      " [3.60000000e+03 4.04417515e-04 5.74388541e-04 8.77380371e-04\n",
      "  1.38798263e-03 2.11798423e-03 6.52207793e-04]\n",
      " [3.65000000e+03 3.83317471e-04 5.44420443e-04 8.31604004e-04\n",
      "  1.32882537e-03 2.05067685e-03 6.30999156e-04]\n",
      " [3.70000000e+03 3.69250774e-04 5.24441712e-04 8.01086426e-04\n",
      "  1.28006120e-03 1.97542273e-03 6.19517683e-04]\n",
      " [3.75000000e+03 3.43143940e-04 4.96749766e-04 7.70568848e-04\n",
      "  1.23129692e-03 1.91754173e-03 6.14266436e-04]\n",
      " [3.80000000e+03 3.29554081e-04 4.77076508e-04 7.40051270e-04\n",
      "  1.19433214e-03 1.84611650e-03 6.08301083e-04]\n",
      " [3.85000000e+03 3.15964222e-04 4.57403250e-04 7.09533691e-04\n",
      "  1.15639425e-03 1.78996241e-03 5.92740747e-04]\n",
      " [3.90000000e+03 3.09169292e-04 4.35532071e-04 6.94274902e-04\n",
      "  1.10938633e-03 1.72768615e-03 5.78964283e-04]\n",
      " [3.95000000e+03 2.95579433e-04 4.27893363e-04 6.63757324e-04\n",
      "  1.07120513e-03 1.68539467e-03 5.57112107e-04]\n",
      " [4.00000000e+03 2.86936760e-04 4.11389396e-04 6.33239746e-04\n",
      "  1.04214728e-03 1.62880821e-03 5.32287541e-04]\n",
      " [4.05000000e+03 2.60710716e-04 3.89104709e-04 6.17980957e-04\n",
      "  1.00225548e-03 1.58459472e-03 4.97839554e-04]\n",
      " [4.10000000e+03 2.54273415e-04 3.79497185e-04 6.02722168e-04\n",
      "  9.77508491e-04 1.54546893e-03 4.66706482e-04]\n",
      " [4.15000000e+03 2.47836113e-04 3.69889662e-04 5.87463379e-04\n",
      "  9.52761446e-04 1.50634313e-03 4.35238981e-04]\n",
      " [4.20000000e+03 2.34961510e-04 3.50674614e-04 5.56945801e-04\n",
      "  9.21027386e-04 1.47027173e-03 4.01716424e-04]\n",
      " [4.25000000e+03 2.34961510e-04 3.40612605e-04 5.56945801e-04\n",
      "  9.03267355e-04 1.44044030e-03 3.73395350e-04]\n",
      " [4.30000000e+03 2.28524208e-04 3.36173922e-04 5.41687012e-04\n",
      "  8.95793783e-04 1.40650908e-03 3.47121374e-04]\n",
      " [4.35000000e+03 2.22086906e-04 3.26704234e-04 5.26428223e-04\n",
      "  8.70560179e-04 1.38970884e-03 3.32671776e-04]\n",
      " [4.40000000e+03 2.15649605e-04 3.17234546e-04 5.11169434e-04\n",
      "  8.61626875e-04 1.36496371e-03 3.21948230e-04]]\n"
     ]
    }
   ],
   "source": [
    "from hepdata_lib import RootFileReader\n",
    "reader = RootFileReader(\"example_inputs/WWmerged.higgsCombineTest.Asymptotic.root\")\n",
    "data = reader.read_limit_tree()\n",
    "print(data)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The `RootFileReader` returned an array of tuples, where each tuple is for a mass point, here starting at 1000 (GeV). If you have one file per mass point, you can easily chain together the ntuples in a loop. If you have text files, use `np.loadtxt(\"limits.txt\", skiprows=0)` as in the [Getting started notebook](Getting_started.ipynb).\n",
    "\n",
    "Let's define the variables and uncertainties that are in the ntuple. The format, i.e. their position in the tuple (observed, expected, ±1/±2 sigma variations) is a standard in `combine`. Each `Uncertainty` is associated to a `Variable` by adding it using `add_uncertainty()`:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [],
   "source": [
    "from hepdata_lib import Variable, Uncertainty\n",
    "d = Variable(\"Bulk graviton mass\", is_independent=True, is_binned=False, units=\"GeV\")\n",
    "d.values = data[:,0]\n",
    "\n",
    "obs = Variable(\"Cross section upper limit at 95% CL\", is_independent=False, is_binned=False, units=\"pb\")\n",
    "obs.values = data[:,6]\n",
    "obs.add_qualifier(\"Limit\", \"Observed\")\n",
    "obs.add_qualifier(\"SQRT(S)\", 13, \"TeV\")\n",
    "obs.add_qualifier(\"LUMINOSITY\", 35.9, \"fb$^{-1}$\")\n",
    "\n",
    "exp = Variable(\"Cross section upper limit at 95% CL\", is_independent=False, is_binned=False, units=\"pb\")\n",
    "exp.values = data[:,3]\n",
    "exp.add_qualifier(\"Limit\", \"Expected\")\n",
    "exp.add_qualifier(\"SQRT(S)\", 13, \"TeV\")\n",
    "exp.add_qualifier(\"LUMINOSITY\", 35.9, \"fb$^{-1}$\")\n",
    "\n",
    "# +/- 1 sigma\n",
    "unc_1s = Uncertainty(\"1 s.d.\", is_symmetric=False)\n",
    "unc_1s.set_values_from_intervals(zip(data[:,2], data[:,4]), nominal=exp.values)\n",
    "exp.add_uncertainty(unc_1s)\n",
    "\n",
    "# +/- 2 sigma\n",
    "unc_2s = Uncertainty(\"2 s.d.\", is_symmetric=False)\n",
    "unc_2s.set_values_from_intervals(zip(data[:,1], data[:,5]), nominal=exp.values)\n",
    "exp.add_uncertainty(unc_2s)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "This is it, all that needs to be done is to add the variables to the `Table` and the `Table` to the `Submission`, and create the files. Again, please refer to the [Getting started notebook](Getting_started.ipynb) for a complete example."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [],
   "source": [
    "table.add_variable(d)\n",
    "table.add_variable(obs)\n",
    "table.add_variable(exp)\n",
    "table.add_additional_resource(\"ROOT file\", \"example_inputs/WWmerged.higgsCombineTest.Asymptotic.root\", copy_file=True)  # optional\n",
    "submission.add_table(table)\n",
    "submission.create_files(\"example_output\",remove_old=True)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Figure_006-a.png\n",
      "WWmerged.higgsCombineTest.Asymptotic.root\n",
      "figure_6a.yaml\n",
      "submission.yaml\n",
      "thumb_Figure_006-a.png\n"
     ]
    }
   ],
   "source": [
    "!ls example_output"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {
    "scrolled": true
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "dependent_variables:\n",
      "- header:\n",
      "    name: Cross section upper limit at 95% CL\n",
      "    units: pb\n",
      "  qualifiers:\n",
      "  - name: Limit\n",
      "    value: Observed\n",
      "  - name: SQRT(S)\n",
      "    units: TeV\n",
      "    value: 13\n"
     ]
    }
   ],
   "source": [
    "!head example_output/figure_6a.yaml"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3 (ipykernel)",
   "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.11.6"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}