talk.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "id": "d5089437",
   "metadata": {},
   "source": [
    "![cabinetry logo](_static/cabinetry_logo_small.png)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "1f0c4ba5-6bc4-42b9-ae82-f3e77a93f9a9",
   "metadata": {},
   "source": [
    "To run this notebook outside of the provided Binder environment, you need `cabinetry`. You can install it with `pip` by uncommenting the cell below. The `[contrib]` setup extra installs libraries that are not needed for the core `cabinetry` logic, but for `ROOT` file handling and plotting. We will also install the `tensorflow` backend for `pyhf` here for a demonstration below, as well as `hist` (again, not needed if running via Binder)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "9129f10d-9c5b-4b98-b463-24f2a489a0a7",
   "metadata": {},
   "outputs": [],
   "source": [
    "#! pip install cabinetry[contrib]\n",
    "#! pip install pyhf[tensorflow,contrib]\n",
    "#! pip install hist[plot]"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "ece924f2-6481-45e4-9132-f42cdbbe2002",
   "metadata": {},
   "source": [
    "<br>\n",
    "\n",
    "`cabinetry` outputs text with the Python `logging` module. We use a helper function to set some sensible defaults, but you are welcome to customize this to your own needs."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "743afc95",
   "metadata": {},
   "outputs": [],
   "source": [
    "import cabinetry\n",
    "cabinetry.set_logging()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "14cf960b-5605-48d2-910c-3359efaa0268",
   "metadata": {},
   "source": [
    "<br>\n",
    "We will also set up additional imports that are used later on in this notebook."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "a18df09f-0293-45ae-bd7f-e9498626ed8c",
   "metadata": {},
   "outputs": [],
   "source": [
    "import copy\n",
    "import json\n",
    "import pathlib\n",
    "\n",
    "import boost_histogram as bh\n",
    "import hist\n",
    "import numpy as np\n",
    "import pyhf\n",
    "from pyhf.contrib.utils import download"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "ce4b21f9-dcab-4efc-8d19-4b879d1544af",
   "metadata": {},
   "source": [
    "<br>\n",
    "\n",
    "We will have a look at three different ways of using `cabinetry` in this notebook:\n",
    "- creating a statistical model,\n",
    "- performing statistical inference with our model,\n",
    "- exploring the statistical model of an ATLAS analysis.\n",
    "\n",
    "There is also bonus material at the end for you to look through."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "96a5baec-6001-4f18-a33d-1929219807c5",
   "metadata": {},
   "source": [
    "<br>\n",
    "\n",
    "# Creating a statistical model\n",
    "Statistical models are constructed following instructions in a configuration. You can specify your model in a configuration in `YAML` or `JSON` format, or alternatively as a Python dictionary. Here we will go the dictionary route. Have a look at the `config_example.yml` file for a `YAML` example."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "375cc5f6-5a23-4c9f-a078-8ff0a8b4cbaa",
   "metadata": {},
   "source": [
    "<br>\n",
    "\n",
    "### Defining the model\n",
    "\n",
    "There are a few things we need. First up, some general settings:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "855f18bb",
   "metadata": {},
   "outputs": [],
   "source": [
    "config = {\n",
    "   \"General\":{\n",
    "      \"Measurement\": \"minimal_example\",\n",
    "      \"POI\": \"Signal_norm\",              # parameter of interest, which we want to measure \n",
    "      \"InputPath\": \"input/{SamplePath}\", # where to find input data\n",
    "      \"HistogramFolder\": \"histograms/\"\n",
    "\n",
    "   }\n",
    "}"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "86bbf3dc-ce51-4c94-ab76-c8291f9f3780",
   "metadata": {},
   "source": [
    "Note the `input/{SamplePaths}` structure, which we will get back to shortly.\n",
    "\n",
    "Now it is time to think more about physics, let's define a phase space region that contains data we want to fit to."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "d68a639c-32c7-44d6-854b-f362b3480d38",
   "metadata": {},
   "outputs": [],
   "source": [
    "config.update({\n",
    "   \"Regions\":[\n",
    "      {\n",
    "         \"Name\": \"Signal_region\",\n",
    "         \"Filter\": \"lep_charge > 0\",           # event selection \n",
    "         \"Variable\": \"jet_pt\",                 # which variable we bin histograms in\n",
    "         \"Binning\": [200, 300, 400, 500, 600]\n",
    "      }\n",
    "   ]\n",
    "})"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "e5a94400-d7f3-4ebe-9d9f-5ceb24347b85",
   "metadata": {},
   "source": [
    "`\"Regions\"` is a list, because we can use events from more than one phase space region.\n",
    "\n",
    "We also need to specify where our observed data is stored. Considerations about physics become important now: which types of processes do we expect to show up in this phase space region? We list them below, those samples are typically simulated with Monte Carlo methods."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "f2ccc354-efba-4616-8541-af22fd8244ce",
   "metadata": {},
   "outputs": [],
   "source": [
    "config.update({\n",
    "   \"Samples\":[\n",
    "      {\n",
    "         \"Name\": \"Data\",\n",
    "         \"Tree\": \"pseudodata\",\n",
    "         \"SamplePath\": \"data.root\",\n",
    "         \"Data\": True                       # observed data is handled differently, need to distinguish\n",
    "      },\n",
    "      {\n",
    "         \"Name\": \"Signal\",\n",
    "         \"Tree\": \"signal\",\n",
    "         \"SamplePath\": \"prediction.root\",\n",
    "         \"Weight\": \"weight\"                 # weights: Monte Carlo integration, simulation correction etc.\n",
    "      },\n",
    "      {\n",
    "         \"Name\": \"Background\",\n",
    "         \"Tree\": \"background\",\n",
    "         \"SamplePath\": \"prediction.root\",\n",
    "         \"Weight\": \"weight\"\n",
    "      }\n",
    "   ]\n",
    "})"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "3ad101b5-01ff-480e-8cb0-5d545cde1d95",
   "metadata": {},
   "source": [
    "`\"Samples\"` is again a list, allowing us to include arbitrarly many different types of processes we may need to consider. The `\"SamplePath\"` option replaces the placeholder in `input/{SamplePath}`, and specifies where each sample can be found (see the [documentation](https://cabinetry.readthedocs.io/en/latest/core.html#input-file-path-specification-for-ntuples)).\n",
    "\n",
    "We are almost done defining our statistical model. What about systematic uncertainties? For now, we won't define any."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "456acd07-d928-44a2-9065-e37037c568dd",
   "metadata": {},
   "outputs": [],
   "source": [
    "config.update({\"Systematics\": []})"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "d0948c92-7257-4104-ad19-863c4092c975",
   "metadata": {},
   "source": [
    "<br>\n",
    "\n",
    "There is one more thing we need though: we said our POI (parameter of interest) was `\"Signal_norm\"`. That is a normalization factor, describing the normalization of somthing. We still need to define what it does!"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "1f99436b-472f-49da-a8d3-1cc90ad5ec0b",
   "metadata": {},
   "outputs": [],
   "source": [
    "config.update({\n",
    "   \"NormFactors\":[\n",
    "      {\n",
    "         \"Name\": \"Signal_norm\",\n",
    "         \"Samples\": \"Signal\",    # we want this parameter to scale the signal\n",
    "         \"Nominal\": 1,\n",
    "         \"Bounds\": [-5, 10]\n",
    "      }\n",
    "   ]\n",
    "})"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "a39b019b-fd8f-48bf-98f5-7e167473d414",
   "metadata": {},
   "source": [
    "<br>\n",
    "\n",
    "It is a good idea to validate that our configuration satisfies the `cabinetry` configuration schema."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "b3ff3ef6-c30b-4eaa-b94a-2921b24cbb1f",
   "metadata": {},
   "outputs": [],
   "source": [
    "cabinetry.configuration.validate(config)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "0c30fc25-4660-4135-94f9-0bec1e2bacd5",
   "metadata": {},
   "source": [
    "<br>\n",
    "\n",
    "That looks good, but is not sufficient to be sure that everything is fine (some things will only show up at runtime). If we use some settings `cabinetry` does not know, or forget some it expects, we would see an error. We can also get some summary information about the configuration we defined."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "33827777-6a3b-4284-aa4f-14ff0b3cc20e",
   "metadata": {},
   "outputs": [],
   "source": [
    "cabinetry.configuration.print_overview(config)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "01b4ee6d-742f-4404-915a-6e228fa2e31b",
   "metadata": {},
   "source": [
    "<br>\n",
    "\n",
    "### Creating histograms\n",
    "\n",
    "We now need to create the required histograms for our statistical model. `cabinetry` implements the logic to figure out which histograms are needed, to create all instructions, and to send those off for execution. You can view our `config` as a generic specification and `cabinetry` as the reference implementation handling it, but you could just as well use your own code to interact with the specification. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "e7015f9b-e9d7-41bd-b9fe-5e8817481fc2",
   "metadata": {},
   "outputs": [],
   "source": [
    "cabinetry.templates.build(config, method=\"uproot\")"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "4e975c57-0deb-4fbf-985e-822fee298daf",
   "metadata": {},
   "source": [
    "`cabinetry` used `uproot` and `boost-histogram` to create three histograms for us: the distribution of our three samples in the one phase space region we defined. We also see a warning: there are no expected signal events in the first bin of the histogram.\n",
    "\n",
    "The histograms are saved to the folder specified under HistogramFolder in the General settings in the configuration file. In this case, this folder is histograms/:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "862f5cfb",
   "metadata": {},
   "outputs": [],
   "source": [
    "!ls histograms"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "9a60e180",
   "metadata": {},
   "source": [
    "A post-processing step can be run to apply optional operations like histogram smoothing. See the bonus material."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "b82677df-86a5-4f51-b515-ae2bee262984",
   "metadata": {},
   "outputs": [],
   "source": [
    "cabinetry.templates.postprocess(config)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "a7d25704-7724-47c9-b841-7cfcf70c978e",
   "metadata": {},
   "source": [
    "You can also provide existing histograms you built yourself for `cabinetry` to use, see the [cabinetry-tutorials](https://github.com/cabinetry/cabinetry-tutorials) repository for an example.\n",
    "\n",
    "<br>\n",
    "Let's visualize the histograms we have produced."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "e3260338-ccd7-4a4b-acbe-b29275677a5d",
   "metadata": {},
   "outputs": [],
   "source": [
    "_ = cabinetry.visualize.data_mc_from_histograms(config)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "1f197fa8-d706-4321-a94d-8c7075778224",
   "metadata": {},
   "source": [
    "We see our expected distributions of signal and background, as well as our observed (pseudo-) data."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "f9bec4cc-73e6-4b6a-8b34-4a8c5e768356",
   "metadata": {},
   "source": [
    "<br>\n",
    "\n",
    "### A more complex model: adding systematic uncertainties\n",
    "\n",
    "Let's make our model a bit more realistic and add a few systematic uncertainties:\n",
    "- a 5% luminosity uncertainty,\n",
    "- a `Modeling` uncertainty derived from comparing our nominal background prediction to that of a different simulator,\n",
    "- a `WeightBasedModeling` modeling uncertainty derived from varying the weights we apply to background events."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "a4ef9b7a-504b-4a59-a49a-fad9f9891c80",
   "metadata": {},
   "outputs": [],
   "source": [
    "config.update({\n",
    "   \"Systematics\":[\n",
    "      {\n",
    "         \"Name\": \"Luminosity\",\n",
    "         \"Up\": {\"Normalization\": 0.05},\n",
    "         \"Down\": {\"Normalization\": -0.05},\n",
    "         \"Type\": \"Normalization\"\n",
    "      },\n",
    "      {\n",
    "         \"Name\":\"Modeling\",\n",
    "         \"Up\": {\"Tree\": \"background_varied\"},\n",
    "         \"Down\": {\"Symmetrize\": True},\n",
    "         \"Samples\": \"Background\",\n",
    "         \"Type\": \"NormPlusShape\"\n",
    "      },\n",
    "      {\n",
    "         \"Name\": \"WeightBasedModeling\",\n",
    "         \"Up\": {\"Weight\": \"weight_up\"},\n",
    "         \"Down\": {\"Weight\": \"0.7*weight\"},\n",
    "         \"Samples\": \"Background\",\n",
    "         \"Type\": \"NormPlusShape\"\n",
    "      }\n",
    "   ],\n",
    "})"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "3b4eb839-836c-4c09-a35d-2b885a1450f6",
   "metadata": {},
   "source": [
    "<br>\n",
    "These new systematic uncertainties require new histograms, so let's build those."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "5c9d8c9e-0f07-4ba4-8fbc-99186b10800b",
   "metadata": {},
   "outputs": [],
   "source": [
    "cabinetry.templates.build(config, method=\"uproot\")\n",
    "cabinetry.templates.postprocess(config)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "229772c5-b76c-4dc0-bf6c-e237e22ab0e9",
   "metadata": {},
   "source": [
    "<br>\n",
    "We can visualize the template histograms corresponding to systematic variations of our model."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "915e27e2-9a11-4add-8a48-40cbcb4e94ef",
   "metadata": {},
   "outputs": [],
   "source": [
    "_ = cabinetry.visualize.templates(config)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "94aa9f17-6648-4feb-b078-31b6cc54ffe7",
   "metadata": {},
   "source": [
    "The top figure is our `Modeling` uncertainty, where we compare our nominal background prediction to that of another simulator, here called the \"up\" variation. At the bottom are the weight-based variations: you can see the \"down\" variation, defined by multiplying the nominal weight by `0.7`, is a factor `0.7` smaller than nominal."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "e155a6fe-74ed-4b8b-b270-d3e915281751",
   "metadata": {},
   "source": [
    "<br>\n",
    "\n",
    "### Building a workspace\n",
    "\n",
    "With all relevant histograms available, we can now construct a `pyhf` workspace. This is our serialized fit model. It contains everything needed to construct a likelihood function via `pyhf`, which is then used for inference."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "c2d2f57b-8cbe-49a3-bd74-d01e7e88942a",
   "metadata": {},
   "outputs": [],
   "source": [
    "workspace_path = \"example_workspace.json\"\n",
    "spec = cabinetry.workspace.build(config)\n",
    "cabinetry.workspace.save(spec, workspace_path)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "d36f9d4e",
   "metadata": {},
   "source": [
    "The `spec` object is a dictionary holding our workspace specification."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "da3ee968-125f-4e59-8540-50de9eeda76f",
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "print(json.dumps(spec, sort_keys=True, indent=4))"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "c1d02578-65b5-4da2-989d-112d653354e9",
   "metadata": {},
   "source": [
    "Saving is optional, we could directly continue with the `spec` object. `pyhf` has validated our workspace, and we are now ready for statistical inference."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "210e3dc9",
   "metadata": {},
   "source": [
    "<br>\n",
    "\n",
    "### Model structure\n",
    "\n",
    "It can be helpful to visualize the modifier structure of the statistical model we have built to catch potential issues. The `visualize.modifier_grid` function creates a figure showcasing the information about which modifiers (indicated by color) act on which region and sample when a given parameter (on the horizontal axis) is varied. To split this visualization from one table per region to one table per sample, use `split_by_sample=True`.\n",
    "\n",
    "We need the fit model (containing the probability density function) for the visualization, which we create from the workspace specification. We directly use the `pyhf` API for this here, in the next section you will see the `cabinetry` API used for the same task."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "7d1e93b3",
   "metadata": {},
   "outputs": [],
   "source": [
    "cabinetry.visualize.modifier_grid(pyhf.Workspace(spec).model())"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "6f477b12",
   "metadata": {},
   "source": [
    "Our model here is not very complex, for a more complex example and the accompanying visualization see [these figures](https://github.com/scikit-hep/cabinetry/issues/332#issuecomment-1214245346)."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "5535de2e-8ccc-4ba5-a298-b6e5c5c91129",
   "metadata": {},
   "source": [
    "<br>\n",
    "\n",
    "# Performing statistical inference with our model\n",
    "\n",
    "To perform inference, we need two things: a probability density function (pdf), or `model`, and data to fit it to. Both are derived from the workspace specification."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "b4e51e3c-0a61-4bce-961b-19742a28db81",
   "metadata": {},
   "outputs": [],
   "source": [
    "model, data = cabinetry.model_utils.model_and_data(spec)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "0ec123f2-f6c2-481d-a889-43b5d910f813",
   "metadata": {},
   "source": [
    "You can see our `Signal_norm` normalization showing up, in addition to parameters for the systematic uncertainties we defined.\n",
    "\n",
    "There is also `staterror_Signal_region`: these are parameters automatically created by `cabinetry` to encode systematic uncertainty due to the finite sample size of our predicted distributions for signal and background.\n",
    "\n",
    "`data` is a list, starting with the observed counts per bin."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "de9e18fc-88e6-41dc-9a34-7162d5a6ecd8",
   "metadata": {},
   "outputs": [],
   "source": [
    "data"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "23f65799-2f30-41a7-b390-16a24d5b4ad8",
   "metadata": {},
   "source": [
    "What comes after is so-called auxiliary data. Check out [this `pyhf` tutorial](https://pyhf.github.io/pyhf-tutorial/HelloWorld.html#auxiliary-data) to learn more about this important concept. You can view it as data observed in previous measurements, which inform our current model. In this case the auxiliary data is associated to the `staterror_Signal_region` parameter. "
   ]
  },
  {
   "cell_type": "markdown",
   "id": "6ec8e723-3c9e-4d31-8439-59919ad11d03",
   "metadata": {},
   "source": [
    "<br>\n",
    "\n",
    "### Maximum likelihood estimate (MLE)\n",
    "\n",
    "Let's fit our model to data to obtain the maximum likelihood estimate (MLE)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "4f023a2e-4599-4ee6-9dfd-6fa1c2ecef87",
   "metadata": {},
   "outputs": [],
   "source": [
    "fit_results = cabinetry.fit.fit(model, data)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "b6ad5704-2e04-4232-86b9-120bf5a2e4e2",
   "metadata": {},
   "source": [
    "<br>\n",
    "The fit converged, and we see the best-fit parameter results reported. The results are stored in a named tuple. This allows for easy access of the results."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "3d1a90b9-3005-4e53-aa99-ed754fa6af64",
   "metadata": {},
   "outputs": [],
   "source": [
    "for label, result, unc in zip(fit_results.labels, fit_results.bestfit, fit_results.uncertainty):\n",
    "    print(f\"{label}: {result:.3f} +/- {unc:.3f}\")"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "6a27c735-3084-4032-8fa3-5b38e58fc61b",
   "metadata": {},
   "source": [
    "<br>\n",
    "\n",
    "It is helpful to visualize the fit results. Let's start with the *pull plot* showing us best-fit parameter results."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "8fe42034-b8f4-482c-9fe3-06f64b98c5a1",
   "metadata": {},
   "outputs": [],
   "source": [
    "cabinetry.visualize.pulls(fit_results, exclude=\"Signal_norm\")"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "9ad13e89-9e6c-4c07-834f-401d2222f0e5",
   "metadata": {},
   "source": [
    "<br>\n",
    "\n",
    "The parameter correlation matrix has a handy `pruning_threshold` setting to filter out parameters that are not highly correlated with others."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "af402771-58ff-48b9-a185-1eb4c6c39f37",
   "metadata": {},
   "outputs": [],
   "source": [
    "cabinetry.visualize.correlation_matrix(fit_results, pruning_threshold=0.1)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "04850886-24fa-456f-9971-28fcb785d0a8",
   "metadata": {},
   "source": [
    "<br>\n",
    "How does the model look like after fit to data? Let's first look again at the model before the fit to data. We use information from the workspace for plotting, and from the config for cosmetics (axis label / binning). The results include the effect of systematic uncertainties."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "5c1efef8-1a0a-4c3a-b63f-a13eb50c4661",
   "metadata": {},
   "outputs": [],
   "source": [
    "model_pred = cabinetry.model_utils.prediction(model)\n",
    "figures = cabinetry.visualize.data_mc(model_pred, data, config=config)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "0053892b",
   "metadata": {},
   "source": [
    "It is possible to edit the figures created by `cabinetry` using the `matplotlib` API. The example below modifies an axis label to use $\\LaTeX$. `cabinetry` version 0.5.2 also added the possibility to more conveniently change the colors of histograms in the stack as shown below."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "050370bc",
   "metadata": {},
   "outputs": [],
   "source": [
    "figures = cabinetry.visualize.data_mc(model_pred, data, config=config,\n",
    "                                      colors={\"Signal\": \"salmon\", \"Background\": \"yellowgreen\"}, close_figure=True)\n",
    "ratio_panel = figures[0][\"figure\"].get_axes()[1]\n",
    "ratio_panel.set_xlabel(\"jet $p_T$\")\n",
    "figures[0][\"figure\"]  # show figure"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "e709e1e2",
   "metadata": {},
   "source": [
    "Yield tables can also be created from a model prediction, and compared to data. Optional keyword arguments control whether yields per bin are shown (`per_bin=True`, default) and whether bins summed per region are shown (`per_channel=True`, disabled by default). The yield table is also saved to disk by default, in a format customizable via the `table_format` argument."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "e91fd2ac",
   "metadata": {},
   "outputs": [],
   "source": [
    "_ = cabinetry.tabulate.yields(model_pred, data)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "5dd35e80",
   "metadata": {},
   "source": [
    "This table was also saved to disk."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "40795267",
   "metadata": {},
   "outputs": [],
   "source": [
    "!cat tables/yields_per_bin_pre-fit.txt"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "314d310a-758c-4e5a-bd73-5224181fb308",
   "metadata": {},
   "source": [
    "<br>\n",
    "\n",
    "We create the post-fit version of this plot by passing in `fit_results`."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "9233142b-eba5-4dc2-a5ba-60b3c0ae32e3",
   "metadata": {},
   "outputs": [],
   "source": [
    "model_pred_postfit = cabinetry.model_utils.prediction(model, fit_results=fit_results)\n",
    "_ = cabinetry.visualize.data_mc(model_pred_postfit, data, config=config)\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "82165476-4834-4f60-9511-d0f333f9c790",
   "metadata": {},
   "source": [
    "The contribution of signal in green has increased, consistent with the normalization `Signal_norm` having been fit to a value greater than `1`."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "436761f9-3390-41b6-b379-294cb677be0b",
   "metadata": {},
   "source": [
    "<br>\n",
    "\n",
    "### Expected sensitivity with the Asimov dataset\n",
    "\n",
    "We can evaluate the expected performance of our model with the so-called Asimov dataset (see [arXiv:1007.1727](https://arxiv.org/abs/1007.1727))."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "4f3aa7af-6896-4702-8ed9-d271efdd2635",
   "metadata": {},
   "outputs": [],
   "source": [
    "asimov_data = cabinetry.model_utils.asimov_data(model)\n",
    "_ = cabinetry.fit.fit(model, asimov_data)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "8e02d65b-a5a5-4bd4-9cc6-7c03d61f3918",
   "metadata": {},
   "source": [
    "By definition, none of the parameters are pulled away from their initial values in this fit."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "b6bd9650-3eba-4a29-bd94-8f98c57ab674",
   "metadata": {},
   "source": [
    "<br>\n",
    "\n",
    "### Beyond MLEs: discovery significance and parameter limits\n",
    "\n",
    "Now that we ran a basic fit, let's do something slightly more involved: calculate discovery significance."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "de33aecf-f941-4966-815a-f354454fefab",
   "metadata": {},
   "outputs": [],
   "source": [
    "significance_results = cabinetry.fit.significance(model, data)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "a33abd73-6e58-43a4-8a94-b3a054675170",
   "metadata": {},
   "source": [
    "The results are again packaged up in a named tuple. The observed significance is higher than the expected significance, consistent with fitting a `Signal_norm` value above one.\n",
    "\n",
    "We can also calculate expected and observed 95% confidence level upper parameter limits with the [CLs method](https://doi.org/10.1088%2F0954-3899%2F28%2F10%2F313). The implementation uses Brent bracketing to find `CLs=0.05` crossings. Let's use slightly different data for this: it is simple to switch out!"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "02b1f8d6-ec90-4273-9a42-d2edf9ad1c2f",
   "metadata": {},
   "outputs": [],
   "source": [
    "data_limit = [112, 129, 92, 63] + model.config.auxdata  # need auxiliary data as well\n",
    "limit_results = cabinetry.fit.limit(model, data_limit)\n",
    "cabinetry.visualize.limit(limit_results)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "1888691c-0932-4097-9b1d-51f0513b6c23",
   "metadata": {},
   "source": [
    "The observed limit follows the expectation quite closely with this new data we have used, which is very similar to the background-only prediction."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "0b1cacf9-93bb-4025-907d-10227344d132",
   "metadata": {},
   "source": [
    "<br>\n",
    "\n",
    "# Exploring the statistical model of an ATLAS analysis\n",
    "\n",
    "Let's explore the statistical model used in an ATLAS search for electroweakinos: [Eur.Phys.J.C 80 (2020) 8, 691](https://inspirehep.net/literature/1755298). The corresponding likelihoods are available on HEPData: [doi:10.17182/hepdata.90607.v2](https://doi.org/10.17182/hepdata.90607.v2).\n",
    "\n",
    "We can download it, pick one of the available signals with `pyhf`, and are ready to go."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "f8ce2a88-c33b-4bd0-9145-1c0939b8f864",
   "metadata": {},
   "outputs": [],
   "source": [
    "download(\"https://www.hepdata.net/record/resource/1267798?view=true\", \"electroweakinos\")\n",
    "ATLAS_ws = pyhf.Workspace(json.load(open(\"electroweakinos/1Lbb-likelihoods-hepdata/BkgOnly.json\")))\n",
    "patchset = pyhf.PatchSet(json.load(open(\"electroweakinos/1Lbb-likelihoods-hepdata/patchset.json\")))\n",
    "ATLAS_ws = patchset.apply(ATLAS_ws, \"C1N2_Wh_hbb_700_400\")\n",
    "cabinetry.workspace.save(ATLAS_ws, \"electroweakinos.json\")"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "74a8425c-d680-40c4-acd0-bf09b0d68ac8",
   "metadata": {},
   "source": [
    "<br>`pyhf` has a helpful command line interface utility to summarize the workspace content."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "4f0668ab-7874-478c-8820-4b32eb86b058",
   "metadata": {},
   "outputs": [],
   "source": [
    "!pyhf inspect electroweakinos.json | head -n 6"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "b5a5f722-8927-444f-92a7-6659c725d1e5",
   "metadata": {},
   "source": [
    "<br>\n",
    "\n",
    "Let's extract model and data with `cabinetry`."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "1eed7e5c-9494-49a6-9ec0-6b74ea36bfd0",
   "metadata": {
    "scrolled": true,
    "tags": []
   },
   "outputs": [],
   "source": [
    "ATLAS_model, ATLAS_data = cabinetry.model_utils.model_and_data(ATLAS_ws)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "9ed5710d-5294-42fa-92af-1cf3757d3f7f",
   "metadata": {},
   "source": [
    "<br>\n",
    "We are ready to take a closer look at the content. Which phase space regions are included?"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "26267735-2a7d-4c5c-81e6-49a0d1da9572",
   "metadata": {},
   "outputs": [],
   "source": [
    "ATLAS_model_pred = cabinetry.model_utils.prediction(ATLAS_model)\n",
    "_ = cabinetry.visualize.data_mc(ATLAS_model_pred, ATLAS_data)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "ed5cbcd9-6171-4d3e-8889-cd9fba546641",
   "metadata": {},
   "source": [
    "<br>\n",
    "Let's fit this model to data."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "3548fe2a-320f-467f-a0e4-eb861f823f2f",
   "metadata": {
    "scrolled": true,
    "tags": []
   },
   "outputs": [],
   "source": [
    "ATLAS_fit_results = cabinetry.fit.fit(ATLAS_model, ATLAS_data)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "969e0c54-f7ed-48a2-9839-edfce1a1a5ba",
   "metadata": {},
   "source": [
    "<br>\n",
    "To parse the results more easily, let's visualize them."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "add7f5dc-8291-4382-ad65-71350d870e9b",
   "metadata": {},
   "outputs": [],
   "source": [
    "cabinetry.visualize.pulls(ATLAS_fit_results)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "8b639bc8-9a9d-4be7-9f07-192da13fc206",
   "metadata": {},
   "source": [
    "<br>\n",
    "\n",
    "# Bonus material\n",
    "\n",
    "Additional material that probably does not fit into the main talk. Check it out for more information!"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "3860400b-60ad-4977-85a3-7c9c48be2425",
   "metadata": {},
   "source": [
    "<br>\n",
    "\n",
    "### Command line interface\n",
    "\n",
    "`cabinetry` also provides a [command line interface](https://cabinetry.readthedocs.io/en/latest/cli.html). Feel free to explore it!"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "99f0af6f-d7f8-439e-9611-3bbea0fcdc2c",
   "metadata": {},
   "outputs": [],
   "source": [
    "! cabinetry --help"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "e50ce638-12a7-467d-a571-c7750cf0602b",
   "metadata": {},
   "source": [
    "<br>\n",
    "\n",
    "### Input histograms via a custom function\n",
    "\n",
    "It is possible to inject custom code into `cabinetry`, which is used for template histogram creation. You can find more details in the [documentation](https://cabinetry.readthedocs.io/en/latest/advanced.html#overrides-for-template-building)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "46e10381-eb2c-4833-8d20-f60a4064601a",
   "metadata": {},
   "outputs": [],
   "source": [
    "my_router = cabinetry.route.Router()\n",
    "\n",
    "# define a custom template builder function that is executed for data samples\n",
    "@my_router.register_template_builder(sample_name=\"Data\")\n",
    "def build_data_hist(region: dict, sample: dict, systematic: dict, template: str):\n",
    "    hist = bh.Histogram(\n",
    "        bh.axis.Variable(region[\"Binning\"], underflow=False, overflow=False),\n",
    "        storage=bh.storage.Weight(),\n",
    "    )\n",
    "    yields = np.asarray([17, 12, 25, 20])\n",
    "    variance = np.asarray([1.5, 1.2, 1.8, 1.6])\n",
    "    hist[...] = np.stack([yields, variance], axis=-1)\n",
    "    return hist  # return a boost-histogram histogram\n",
    "\n",
    "custom_config = copy.deepcopy(config)\n",
    "custom_config[\"General\"][\"HistogramFolder\"] = \"histograms_custom/\"\n",
    "cabinetry.templates.build(custom_config, router=my_router)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "9505c159-3019-4b26-bdbf-2059ce756b8f",
   "metadata": {},
   "source": [
    "The histogram creation called our function to create the data histogram.\n",
    "\n",
    "We can load the histogram and check the yields to verify they were correctly picked up."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "682a5c6b-3934-4908-8994-ebb3c675f88b",
   "metadata": {},
   "outputs": [],
   "source": [
    "h = cabinetry.histo.Histogram.from_path(pathlib.Path(\"histograms_custom/Signal_region_Data.npz\"))\n",
    "h.values()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "e1bf2189-4c7a-4351-9f41-8f4aeb4eeaaa",
   "metadata": {},
   "source": [
    "<br>\n",
    "\n",
    "### Histogram smoothing\n",
    "\n",
    "This shows how to apply smoothing to a histogram after producing it. We copy our config and save the resulting histograms in a new location to not interfere with the model we have used so far."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "4c0e8a13-4055-4ecb-a1b9-45b75608e270",
   "metadata": {},
   "outputs": [],
   "source": [
    "smoothing_config = copy.deepcopy(config)\n",
    "smoothing_config[\"General\"][\"HistogramFolder\"] = \"histograms_smoothing/\"\n",
    "smoothing_config.update({\n",
    "   \"Systematics\":[\n",
    "      {\n",
    "         \"Name\":\"Modeling\",\n",
    "         \"Up\": {\"SamplePath\": \"prediction.root\", \"Tree\": \"background_varied\"},\n",
    "         \"Down\": {\"Symmetrize\": True},\n",
    "         \"Samples\": \"Background\",\n",
    "         \"Smoothing\": {\"Algorithm\": \"353QH, twice\"},  # smoothing applied\n",
    "         \"Type\": \"NormPlusShape\"\n",
    "      }\n",
    "   ]\n",
    "})\n",
    "cabinetry.templates.build(smoothing_config, method=\"uproot\")\n",
    "cabinetry.templates.postprocess(smoothing_config)\n",
    "_ = cabinetry.visualize.templates(smoothing_config)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "9a95e8e7-c601-4fcd-938f-34097b37e135",
   "metadata": {},
   "source": [
    "The original histogram is shown with the dotted line, and the dashed line shows the histogram after applying the [353QH, twice](https://cds.cern.ch/record/186223/) algorithm (same as `TH1::SmoothArray` in `ROOT`)."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "f2317609-be7d-4053-938c-531f18dea3f7",
   "metadata": {},
   "source": [
    "<br>\n",
    "\n",
    "### Exploring histograms with `hist`\n",
    "\n",
    "`cabinetry` uses `boost-histogram` internally to handle histograms. That means we can use `hist` for visualization. There are multiple ways for loading histograms, the method below uses information from the config."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "ae1b3b94-ed37-4047-a73e-ac93226768ff",
   "metadata": {},
   "outputs": [],
   "source": [
    "cabinetry_histogram = cabinetry.histo.Histogram.from_config(\n",
    "    config[\"General\"][\"HistogramFolder\"],\n",
    "    config[\"Regions\"][0],  # we only have one region: Signal_region\n",
    "    config[\"Samples\"][0],  # data\n",
    "    {\"Name\": \"Nominal\"},   # no systematics\n",
    ")\n",
    "_ = hist.Hist(cabinetry_histogram).plot()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "231f902c-e1e4-46c5-a220-e7239060bb06",
   "metadata": {},
   "source": [
    "<br>\n",
    "\n",
    "### Using a different `pyhf` backend"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "c3d1b722-cfc8-4245-9727-60d090a6c9c6",
   "metadata": {},
   "source": [
    "If your fits take a long time, or you have a GPU available and would like to make use of it, you can switch the computational backend used by `pyhf`. Set the backend via `pyhf` as you would in standalone use, and `cabinetry` will respect that choice."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "3bd2316a-ac94-4a3c-8ca3-93935b100867",
   "metadata": {},
   "outputs": [],
   "source": [
    "pyhf.set_backend(\"tensorflow\")\n",
    "_ = cabinetry.fit.fit(model, data)\n",
    "pyhf.set_backend(\"numpy\")  # switch back to numpy"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "26477db9-1348-4b6e-91d3-3ac9c6db52b1",
   "metadata": {},
   "source": [
    "<br>\n",
    "\n",
    "### Going back and forth between `cabinetry` and `pyhf`\n",
    "\n",
    "You can go back and forth between `cabinetry` and `pyhf`. Here is a simple example that builds the Asimov dataset with `cabinetry` and performs a fit with `pyhf`."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "a77dc37b-092c-41c0-bc0c-e01229758f26",
   "metadata": {},
   "outputs": [],
   "source": [
    "asimov_data = cabinetry.model_utils.asimov_data(model)\n",
    "pyhf.infer.mle.fit(asimov_data, model)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "df8836da-b0bc-4ca9-992c-1c10494474df",
   "metadata": {},
   "source": [
    "<br>\n",
    "\n",
    "You can also take advantage of the `pyhf` workspace format and edit parameters directly and then build a new model. Let's fix the `Modeling` nuisance parameter to `-0.5` and repeat the fit, this time with `cabinetry`. The model is obtained directly via `pyhf` here instead of using `cabinetry.model_utils.model_and_data`."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "8d85cc1f-c911-4212-bfa8-522130250ac8",
   "metadata": {},
   "outputs": [],
   "source": [
    "spec_edited = copy.deepcopy(spec)\n",
    "spec_edited[\"measurements\"][0][\"config\"][\"parameters\"].append({\"name\": \"Modeling\", \"inits\": [-0.5], \"fixed\": True})\n",
    "model_edited = pyhf.Workspace(spec_edited).model()\n",
    "_ = cabinetry.fit.fit(model_edited, data)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "b019b430-af3f-4e3c-a31b-7b8b36fc608e",
   "metadata": {},
   "source": [
    "The `Modeling` parameter shows up at its fixed value `-0.5` in the fit results, with an uncertainty of `0` since it is set to constant."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "415db153-d873-4502-86d0-bb160cc1b20a",
   "metadata": {},
   "source": [
    "<br>\n",
    "\n",
    "### Further MLE options: MINOS and goodness-of-fit\n",
    "\n",
    "The `cabinetry.fit.fit` API allows use of the [MINOS algorithm](https://iminuit.readthedocs.io/en/stable/reference.html#iminuit.Minuit.minos) to compute confidence intervals via parameter scans. All parameters for which the algorithm should be run have to be listed."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "e0e43f11-064e-48a2-806e-825fcf540aee",
   "metadata": {},
   "outputs": [],
   "source": [
    "_ = cabinetry.fit.fit(model, data, minos=[\"Signal_norm\", \"Modeling\"])"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "f19b130a-4fba-4b1d-9cfc-af6d73cd9422",
   "metadata": {},
   "source": [
    "MINOS uncertainties are not restricted to be symmetric by construction.\n",
    "\n",
    "`cabinetry` also implementes a goodness-of-fit calculation using the [saturated model](http://www.physics.ucla.edu/~cousins/stats/cousins_saturated.pdf)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "d6ffa762-cd64-473d-a3c4-8ba4cea21a93",
   "metadata": {},
   "outputs": [],
   "source": [
    "_ = cabinetry.fit.fit(model, data, goodness_of_fit=True)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "8c15819e-891d-44da-992f-5760c7c8bc73",
   "metadata": {},
   "source": [
    "<br>\n",
    "\n",
    "### Likelihood scan\n",
    "\n",
    "We can perform a likelihood scan, where we hold a parameter fixed at different values and perform a MLE for the other parameters, recording the likelihood at each step in the scan."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "31410b9c-40e3-4bae-9bbd-c3572263877c",
   "metadata": {
    "scrolled": true,
    "tags": []
   },
   "outputs": [],
   "source": [
    "scan_results = cabinetry.fit.scan(model, data, \"Signal_norm\", n_steps=5)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "f1df702b-eb1f-4a71-8cde-0c5f1dd41d07",
   "metadata": {},
   "source": [
    "<br>\n",
    "The scan in this example shows relatively good agreement with a naive Gaussian approximation."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "2113f255-7e77-4635-b625-51ef440dd208",
   "metadata": {},
   "outputs": [],
   "source": [
    "cabinetry.visualize.scan(scan_results)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "fcd6ceb7-63f1-4bad-a62a-fb8c46abbb47",
   "metadata": {},
   "source": [
    "<br>\n",
    "\n",
    "### Nuisance parameter ranking\n",
    "\n",
    "We can rank nuisance parameters (NPs) by their impact on the POI: how much does the POI change if the NP varies within its uncertainty? This requires a lot of MLE fits."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "c154b331-66dc-472a-8307-ec9b6c527895",
   "metadata": {
    "scrolled": true,
    "tags": []
   },
   "outputs": [],
   "source": [
    "ranking_results = cabinetry.fit.ranking(model, data, fit_results=fit_results)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "e914d135-8e56-4992-b275-60aae493c579",
   "metadata": {},
   "source": [
    "<br>\n",
    "The figure visualizes the impact of nuisance parameters on the POI in order of decreasing impact."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "ec1737c3-84a4-4469-a8ab-ad3a454d8513",
   "metadata": {},
   "outputs": [],
   "source": [
    "cabinetry.visualize.ranking(ranking_results)"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "stats",
   "language": "python",
   "name": "stats"
  },
  "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.16"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}