examples/reading_histograms.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Reading histograms\n",
    "\n",
    "One of the most common tasks will be translating [ROOT](https://root.cern.ch) histograms into the HEPData format. `hepdata_lib` will help you with that, and this notebook will demonstrate how to do that.\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 plot will be a `Table`, in this example Figure 4a from page 12 (upper left) of the [publication](https://cms-results.web.cern.ch/cms-results/public-results/publications/B2G-17-009/index.html), which shows the distribution of the reconstructed B quark mass for the data as well as the individual simulated processes. 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 4a\")\n",
    "table.description = \"Distribution in the reconstructed B quark mass, after applying all selections to events with no forward jet, compared to the background distributions estimated before fitting. The plot refers to the low-mass mB analysis. The expectations for signal MC events are given by the blue histogram lines. Different contributions to background are indicated by the colour-filled histograms. The grey-hatched error band shows total uncertainties in the background expectation. The ratio of observations to background expectations is given in the lower panel, together with the total uncertainties prior to fitting, indicated by the grey-hatched band.\"\n",
    "table.location = \"Data from Figure 4 (upper left), located on page 12.\"\n",
    "table.keywords[\"observables\"] = [\"N\"]\n",
    "table.add_image(\"example_inputs/CMS-B2G-17-009_Figure_004-a.pdf\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The individual plot components are stored in different ROOT files, one for the individual background processes (one histogram per process plus the total), another one for the data, and the third for the signal process. All histograms here are of type [TH1](https://root.cern.ch/doc/master/classTH1.html), but you can also read in 2-dimensional [TH2](https://root.cern.ch/doc/master/classTH2.html) using `read_hist_2d(...)` instead of `read_hist_1d(...)`:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [],
   "source": [
    "from hepdata_lib import RootFileReader\n",
    "\n",
    "reader = RootFileReader(\"example_inputs/mlfit_lm_1000.root\")\n",
    "reader_data = RootFileReader(\"example_inputs/Data_cat0_singleH.root\")\n",
    "reader_signal = RootFileReader(\"example_inputs/BprimeBToHB1000_cat0_singleH.root\")\n",
    "\n",
    "TotalBackground = reader.read_hist_1d(\"shapes_prefit/cat0_singleH/total_background\")\n",
    "TT = reader.read_hist_1d(\"shapes_prefit/cat0_singleH/TT\")\n",
    "QCD = reader.read_hist_1d(\"shapes_prefit/cat0_singleH/QCDTT\")\n",
    "WJets = reader.read_hist_1d(\"shapes_prefit/cat0_singleH/WJets\")\n",
    "ZJets = reader.read_hist_1d(\"shapes_prefit/cat0_singleH/ZJets\")\n",
    "\n",
    "Data = reader_data.read_hist_1d(\"h_bprimemass_SRlm\")\n",
    "\n",
    "signal = reader_signal.read_hist_1d(\"h_bprimemass_SRlm\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The content of the histograms is stored as a dictionary, with keys `x` (bin center), `y` (bin value or for `TH2` the bin center of the 2nd dimension), `z` (`TH2` only: bin value), as well as the bin errors `dy` (`dz` for `TH2`). Furthermore, the lower and upper bin edges (`x_edges`, for `TH2` also `y_edges`) are stored for each bin:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "dict_keys(['x', 'y', 'x_edges', 'x_labels', 'dy'])"
      ]
     },
     "execution_count": 4,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "TotalBackground.keys()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The `RootFileReader` automatically recognises if the histogram has symmetric or assymmetric errors based on [TH1::GetBinErrorOption()](https://root.cern.ch/doc/master/classTH1.html#ac6e38c12259ab72c0d574614ee5a61c7). Symmetric errors are returned if this returns `TH1::kNormal`, in this case (as for the example here) the errors are a plain `float` per bin, otherwise a `tuple` of `float`. The bin edges are always stored as `tuple`:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "x <class 'list'> <class 'float'>\n",
      "y <class 'list'> <class 'float'>\n",
      "x_edges <class 'list'> <class 'tuple'>\n",
      "x_labels <class 'list'> <class 'str'>\n",
      "dy <class 'list'> <class 'float'>\n"
     ]
    }
   ],
   "source": [
    "from __future__ import print_function\n",
    "for key in TotalBackground.keys():\n",
    "    print(key, type(TotalBackground[key]), type(TotalBackground[key][0]))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now define the variables:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [],
   "source": [
    "from hepdata_lib import Variable, Uncertainty\n",
    "\n",
    "# x-axis: B quark mass\n",
    "mmed = Variable(\"$M_{bH}$\", is_independent=True, is_binned=False, units=\"GeV\")\n",
    "mmed.values = signal[\"x\"]\n",
    "\n",
    "# y-axis: N events\n",
    "sig = Variable(\"Number of signal events\", is_independent=False, is_binned=False, units=\"\")\n",
    "sig.values = signal[\"y\"]\n",
    "\n",
    "totalbackground = Variable(\"Number of background events\", is_independent=False, is_binned=False, units=\"\")\n",
    "totalbackground.values = TotalBackground[\"y\"]\n",
    "\n",
    "tt = Variable(\"Number of ttbar events\", is_independent=False, is_binned=False, units=\"\")\n",
    "tt.values = TT[\"y\"]\n",
    "\n",
    "qcd = Variable(\"Number of qcd events\", is_independent=False, is_binned=False, units=\"\")\n",
    "qcd.values = QCD[\"y\"]\n",
    "\n",
    "wjets = Variable(\"Number of wjets events\", is_independent=False, is_binned=False, units=\"\")\n",
    "wjets.values = WJets[\"y\"]\n",
    "\n",
    "zjets = Variable(\"Number of zjets events\", is_independent=False, is_binned=False, units=\"\")\n",
    "zjets.values = ZJets[\"y\"]\n",
    "\n",
    "data = Variable(\"Number of data events\", is_independent=False, is_binned=False, units=\"\")\n",
    "data.values = Data[\"y\"]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "For the data as well as the background total, we will also provide the associated uncertainties:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [],
   "source": [
    "from hepdata_lib import Uncertainty\n",
    "\n",
    "unc_totalbackground = Uncertainty(\"total uncertainty\", is_symmetric=True)\n",
    "unc_totalbackground.values = TotalBackground[\"dy\"]\n",
    "\n",
    "unc_data = Uncertainty(\"Poisson errors\", is_symmetric=True)\n",
    "unc_data.values = Data[\"dy\"]\n",
    "\n",
    "totalbackground.add_uncertainty(unc_totalbackground)\n",
    "data.add_uncertainty(unc_data)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now we can add the variables to the `Table` and the `Table` to the `Submission`, and create the files. Please refer to the [Getting started notebook](Getting_started.ipynb) for a complete example."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [],
   "source": [
    "table.add_variable(mmed)\n",
    "table.add_variable(sig)\n",
    "table.add_variable(totalbackground)\n",
    "table.add_variable(tt)\n",
    "table.add_variable(qcd)\n",
    "table.add_variable(zjets)\n",
    "table.add_variable(wjets)\n",
    "table.add_variable(data)\n",
    "\n",
    "table.add_additional_resource(\"Background ROOT file\", \"example_inputs/mlfit_lm_1000.root\", copy_file=True)  # optional\n",
    "table.add_additional_resource(\"Data ROOT file\", \"example_inputs/Data_cat0_singleH.root\", copy_file=True)  # optional\n",
    "table.add_additional_resource(\"Signal ROOT file\", \"example_inputs/BprimeBToHB1000_cat0_singleH.root\", copy_file=True)  # optional\n",
    "\n",
    "submission.add_table(table)\n",
    "\n",
    "submission.create_files(\"example_output\",remove_old=True)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "dependent_variables:\n",
      "- header:\n",
      "    name: Number of signal events\n",
      "  values:\n",
      "  - value: 0.95335\n",
      "  - value: 0.83057\n",
      "  - value: 3.5395\n",
      "  - value: 1.8212\n",
      "  - value: 0.04857\n",
      "  - value: 0.047358\n"
     ]
    }
   ],
   "source": [
    "!head example_output/figure_4a.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
}