diff --git a/.gitignore b/.gitignore
index bf8e780..ff6f4b0 100644
--- a/.gitignore
+++ b/.gitignore
@@ -133,3 +133,6 @@ _build
*.tar.gz
book/data/
book/1Lbb-likelihoods/
+
+# Mac artifacts
+.DS_Store
diff --git a/binder/requirements.txt b/binder/requirements.txt
index bd260bf..9274e23 100644
--- a/binder/requirements.txt
+++ b/binder/requirements.txt
@@ -4,6 +4,7 @@ ipywidgets~=8.0.7
ipympl~=0.9.3
pandas~=2.0.3
altair~=5.0.1
+rich==13.5.2
# jupyter notebooks
notebook>=7.0.0
jupyterlab>=4.0.3
diff --git a/book/UpperLimitsTable.ipynb b/book/UpperLimitsTable.ipynb
new file mode 100644
index 0000000..55aa5c9
--- /dev/null
+++ b/book/UpperLimitsTable.ipynb
@@ -0,0 +1,589 @@
+{
+ "cells": [
+ {
+ "cell_type": "code",
+ "execution_count": null,
+ "id": "960720c4",
+ "metadata": {},
+ "outputs": [],
+ "source": [
+ "import json\n",
+ "import tarfile\n",
+ "\n",
+ "import matplotlib.pyplot as plt\n",
+ "import numpy as np\n",
+ "import pyhf\n",
+ "import pyhf.contrib.utils\n",
+ "import scipy\n",
+ "from pyhf.contrib.viz import brazil\n",
+ "\n",
+ "pyhf.set_backend(\"numpy\", \"minuit\")\n",
+ "\n",
+ "from rich import box\n",
+ "from rich.console import Console\n",
+ "from rich.table import Table"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "id": "e2335cc8",
+ "metadata": {},
+ "source": [
+ "# Model-Independent Upper Limits\n",
+ "\n",
+ "This section is aiming to cover some common concepts about (model-independent) upper limits. While it won't necessarily be exhaustive, the aim is to summarize the current working knowledge of how High Energy Physics experiments such as ATLAS produce tables such as Table 20 from the Supersymmetry Electroweak 2-lepton, 2-jet analysis ([SUSY-2018-05](https://atlas.web.cern.ch/Atlas/GROUPS/PHYSICS/PAPERS/SUSY-2018-05/)):\n",
+ "\n",
+ "
\n",
+ " \n",
+ " \n",
+ " Table 20: Model-independent upper limits on the observed visible cross-section in the five electroweak search discovery regions, derived using pseudo-experiments. Left to right: background-only model post-fit total expected background, with the combined statistical and systematic uncertainties; observed data; 95 CL upper limits on the visible cross-section ($\\langle A\\varepsilon\\sigma\\rangle_\\mathrm{obs}^{95}$) and on the number of signal events ($\\mathrm{S}_\\mathrm{obs}^{95}$). The sixth column ($\\mathrm{S}_\\mathrm{exp}^{95}$) shows the expected 95 CL upper limit on the number of signal events, given the expected number (and $\\pm1\\sigma$ excursions of the expectation) of background events. The last two columns indicate the confidence level of the background-only hypothesis ($\\mathrm{CL}_\\mathrm{b}$) and discovery $p$-value with the corresponding Gaussian significance ($Z(s=0)$). $\\mathrm{CL}_\\mathrm{b}$ provides a measure of compatibility of the observed data with the signal strength hypothesis at the 95 CL limit relative to fluctuations of the background, and $p(s=0)$ measures compatibility of the observed data with the background-only hypothesis relative to fluctuations of the background. The $p$-value is capped at 0.5.\n",
+ " \n",
+ " \n",
+ "
\n",
+ "\n",
+ "Here, the ATLAS collaboration provided some text describing each of the individual columns. How does that translate to statistical fits and `pyhf`? Let's explore below. We will use the public probability models from the analysis published on [HEPData entry](https://www.hepdata.net/record/116034)."
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "id": "bcfb863f",
+ "metadata": {},
+ "source": [
+ "## Downloading the public probability models\n",
+ "\n",
+ "First, we'll use `pyhf.contrib` to download the models from HEPData:"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": null,
+ "id": "c411ed7c",
+ "metadata": {},
+ "outputs": [],
+ "source": [
+ "pyhf.contrib.utils.download(\n",
+ " \"https://doi.org/10.17182/hepdata.116034.v1/r34\", \"workspace_2L2J\"\n",
+ ")"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "id": "10f9e136",
+ "metadata": {},
+ "source": [
+ "Then use Python's `tarfile` (part of standard library) to extract out the JSONs:"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": null,
+ "id": "acae3295",
+ "metadata": {},
+ "outputs": [],
+ "source": [
+ "with tarfile.open(\"workspace_2L2J/ewk_discovery_bkgonly.json.tgz\") as fp:\n",
+ " bkgonly = pyhf.Workspace(json.load(fp.extractfile(\"ewk_discovery_bkgonly.json\")))\n",
+ "\n",
+ "with tarfile.open(\"workspace_2L2J/ewk_discovery_patchset.json.tgz\") as fp:\n",
+ " patchset = pyhf.PatchSet(json.load(fp.extractfile(\"ewk_discovery_patchset.json\")))"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "id": "05c9738f",
+ "metadata": {},
+ "source": [
+ "and then build our workspace for a particular signal region **DR-Int-EWK**:"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": null,
+ "id": "bfffecd3",
+ "metadata": {},
+ "outputs": [],
+ "source": [
+ "ws = patchset.apply(bkgonly, \"DR_Int_EWK\")\n",
+ "\n",
+ "model = ws.model()"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "id": "92851f11",
+ "metadata": {},
+ "source": [
+ "## Understanding expected and observed events\n",
+ "\n",
+ "The background estimate(s) in the signal region(s) are obtained from the background-only fit to control region(s) extrapolated to the signal region(s). \n",
+ "\n",
+ "This number, combined with the number of observed events in the signal region, allows us to redo the limit-setting or $p$-value determination without the public probability models.\n",
+ "\n",
+ "For the uncertainties on the background estimate(s) in the signal region(s), they're symmetric and calculated through standard error propagation described below."
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "id": "c5e284e2",
+ "metadata": {},
+ "source": [
+ "### Observed Events\n",
+ "\n",
+ "The observed events are just `ws.data(model)` which we will call `data`. These are ordered in the same way as the ordering of the channels in the model (`model.config.channels`):"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": null,
+ "id": "d0ee1d9d",
+ "metadata": {},
+ "outputs": [],
+ "source": [
+ "data = ws.data(model)\n",
+ "print(list(zip(model.config.channels, data)))"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "id": "e7068fd4",
+ "metadata": {},
+ "source": [
+ "### Expected Events\n",
+ "\n",
+ "The expected events are determined from a background-only fit. This is done by fixing the parameter of interest (i.e. for a search the non-SM signal strength) to zero. This fit will give us the fitted parameters $\\hat{\\theta}$. From this, we can print out the expected background events in each channel after this fit:"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": null,
+ "id": "f79763d5",
+ "metadata": {},
+ "outputs": [],
+ "source": [
+ "pars_uncrt_bonly, corr_bonly = pyhf.infer.mle.fixed_poi_fit(\n",
+ " 0.0, data, model, return_uncertainties=True, return_correlations=True\n",
+ ")\n",
+ "pars_bonly, uncrt_bonly = pars_uncrt_bonly.T\n",
+ "\n",
+ "table = Table(title=\"Per-region yields from background-only fit\", box=box.HORIZONTALS)\n",
+ "table.add_column(\"Region\", justify=\"left\", style=\"cyan\", no_wrap=True)\n",
+ "table.add_column(\"Total Bkg.\", justify=\"center\")\n",
+ "\n",
+ "for region, count in zip(model.config.channels, model.expected_actualdata(pars_bonly)):\n",
+ " table.add_row(region, f\"{count:0.2f}\")\n",
+ "\n",
+ "console = Console()\n",
+ "console.print(table)"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "id": "3940313d",
+ "metadata": {},
+ "source": [
+ "The uncertainty on the total background (expected event rate) is done by *linearly* propagating errors on parameters using the covariance matrix from a fit result. The error is calculated in the same way that ROOT's [`RooAbsReal::getPropagatedError()`](https://root.cern.ch/doc/master/classRooAbsReal.html#a4d7678837410aabcd93e6f0937de525b) calculates it. It is summarized as follows\n",
+ "\n",
+ "$$\n",
+ "\\mathrm{error}^2(x) = F_\\mathbf{a}(x) \\cdot \\mathrm{Corr}(\\mathbf{a},\\mathbf{a}') \\cdot F_{\\mathbf{a}'}^{\\mathrm{T}}(x)\n",
+ "$$\n",
+ "\n",
+ "where \n",
+ "\n",
+ "$$\n",
+ "F_\\mathbf{a}(x) = \\frac{ f(x, \\mathbf{a} + \\mathrm{d}\\mathbf{a}) - f(x, \\mathbf{a} - \\mathrm{d}\\mathbf{a}) }{2},\n",
+ "$$\n",
+ "\n",
+ "with $f(x)$ the probability model and $\\mathrm{d}\\mathbf{a}$ the vector of one-sigma uncertainties of all fit parameters taken from the fit result and $\\mathrm{Corr}(\\mathbf{a},\\mathbf{a}')$ the correlation matrix from the fit result.\n",
+ "\n",
+ "References:\n",
+ "- [Statistical Error Propagation (J. Phys. Chem. A 2001, 105, 15, 3917–3921)](https://web.stanford.edu/~kimth/www-mit/8.13/error_propagation.pdf)\n",
+ "- [root-project/root#13404](https://github.com/root-project/root/issues/13404)\n",
+ "\n",
+ "**Note**: be aware that some definitions of error propagation use covariance and some use correlation. You can translate between the two, turning covariance into correlation by multiplying by $\\partial a$ on both sides. If one uses covariance, you'll need information about the derivatives (e.g. through auto-differentiation).\n",
+ "\n",
+ "This could be done with a for-loop like so:\n",
+ "\n",
+ "```python\n",
+ "up_yields = []\n",
+ "dn_yields = []\n",
+ "for i, variation in enumerate(uncrt_bonly):\n",
+ " varied_up_pars = pars_bonly[:] # copy\n",
+ " varied_up_pars[i] += variation\n",
+ " \n",
+ " varied_dn_pars = pars_bonly[:] # copy\n",
+ " varied_dn_pars[i] -= variation\n",
+ "\n",
+ " up_yields.append(model.expected_actualdata(varied_up_pars))\n",
+ " dn_yields.append(model.expected_actualdata(varied_dn_pars))\n",
+ "```\n",
+ "\n",
+ "which gives us the up/down variations on the yields for each systematic. However, we'll batch this calculation to calculate all yield variations simultaneously. To do this, we set up our model to allow for batching."
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": null,
+ "id": "34c480b3",
+ "metadata": {},
+ "outputs": [],
+ "source": [
+ "npars = len(model.config.parameters)\n",
+ "model_batch = ws.model(batch_size=npars)\n",
+ "pars_batch = np.tile(pars_bonly, (npars, 1))"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "id": "ebf6a42f",
+ "metadata": {},
+ "source": [
+ "Then we calculate the `up_yields` and `dn_yields` which are used to calculate the $F_\\mathbf{a}(x)$ term above, denoted as `variations`:"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": null,
+ "id": "0a9f6848",
+ "metadata": {},
+ "outputs": [],
+ "source": [
+ "up_yields = model_batch.expected_actualdata(pars_batch + np.diag(uncrt_bonly))\n",
+ "dn_yields = model_batch.expected_actualdata(pars_batch - np.diag(uncrt_bonly))\n",
+ "variations = (up_yields - dn_yields) / 2"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "id": "4fd1b3b9",
+ "metadata": {},
+ "source": [
+ "Now we can calculate $\\mathrm{error}^2(x)$:"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": null,
+ "id": "bcc4aa16",
+ "metadata": {},
+ "outputs": [],
+ "source": [
+ "error_sq = np.einsum(\"il,ik,kl->l\", variations, corr_bonly, variations)\n",
+ "error_sq"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "id": "409d3365",
+ "metadata": {},
+ "source": [
+ "Taking the square root will give us the propagated uncertainty on the yields per-region for the background-only fit. We can update our table:"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": null,
+ "id": "0c3023bf",
+ "metadata": {},
+ "outputs": [],
+ "source": [
+ "table = Table(title=\"Per-region yields from background-only fit\", box=box.HORIZONTALS)\n",
+ "table.add_column(\"Region\", justify=\"left\", style=\"cyan\", no_wrap=True)\n",
+ "table.add_column(\"Total Bkg.\", justify=\"center\")\n",
+ "\n",
+ "for region, count, uncrt in zip(\n",
+ " model.config.channels, model.expected_actualdata(pars_bonly), np.sqrt(error_sq)\n",
+ "):\n",
+ " table.add_row(region, f\"{count:0.2f} ± {uncrt:0.2f}\")\n",
+ "\n",
+ "console = Console()\n",
+ "console.print(table)"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "id": "8123d0f3",
+ "metadata": {},
+ "source": [
+ "## Understanding visible cross-section\n",
+ "\n",
+ "_See [ATL-PHYS-PUB-2022-017](https://atlas.web.cern.ch/Atlas/GROUPS/PHYSICS/PUBNOTES/ATL-PHYS-PUB-2022-017/) for descriptions of detector efficiency and selection acceptance._\n",
+ "\n",
+ "**Production cross section upper limits** is another way of saying model-dependent limits. Without the detector efficiency (but keeping selection acceptance), the term **visible cross section upper limits** is used. \n",
+ "\n",
+ "For upper limits on numbers of events (and the corresponding cross section), after detector efficiency and selection acceptance effects, the term **model independent upper limits** is used. These upper limits are typically free from assumptions about the acceptance and efficiency of any model (and are model independent in that sense), but they do make the assumption that signal contamination in control regions is negligible.\n",
+ "\n",
+ "### Model-independent Upper Limits\n",
+ "\n",
+ "Model-independent upper limits (UL) at the 95% confidence level (CL) on the number of beyond the SM (BSM) events for each signal region are derived using the $\\mathrm{CL}_\\mathrm{s}$ prescription (DOI: [10.1088/0954-3899/28/10/313](https://doi.org/10.1088/0954-3899/28/10/313)) and neglecting any possible contamination in the control regions.\n",
+ "\n",
+ "To obtain the 95% CL model-independent UL / discovery $p_0$-value, the fit in the signal region proceeds in the same way as the background-only fit, except that an additional parameter for the non-SM signal strength, constrained to be non-negative, is fit."
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": null,
+ "id": "57fda65c",
+ "metadata": {},
+ "outputs": [],
+ "source": [
+ "pyhf.set_backend(\"numpy\", \"scipy\")\n",
+ "\n",
+ "mu_tests = np.linspace(15, 25, 10)\n",
+ "(\n",
+ " obs_limit,\n",
+ " exp_limit_and_bands,\n",
+ " (poi_tests, tests),\n",
+ ") = pyhf.infer.intervals.upper_limits.upper_limit(\n",
+ " data, model, scan=mu_tests, level=0.05, return_results=True\n",
+ ")"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": null,
+ "id": "f4190909",
+ "metadata": {},
+ "outputs": [],
+ "source": [
+ "print(f\"Observed UL: {obs_limit:0.2f}\")\n",
+ "print(\n",
+ " f\"Expected UL: {exp_limit_and_bands[2]:0.2f} [{exp_limit_and_bands[1]-exp_limit_and_bands[2]:0.2f}, +{exp_limit_and_bands[3]-exp_limit_and_bands[2]:0.2f}]\"\n",
+ ")"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "id": "27c5e954",
+ "metadata": {},
+ "source": [
+ "When calculating $S^{95}_\\mathrm{obs}$ (also for expected) it is typical to inject a _test_ signal (all expected event rates equal to unity) which means for this signal, we're assuming $A\\times\\varepsilon\\times\\sigma\\times\\mathcal{L}=1$. As $S^{95}_\\mathrm{obs}$ is an upper limit on $\\mu$ for this _test_ signal, normalizing this by the integrated luminosity of the data sample allows us to interpret as upper limits on the **visible (BSM) cross-section**, $\\langle A\\varepsilon\\sigma\\rangle_\\mathrm{obs}^{95}$. This is defined as the product of acceptance, reconstruction efficiency and production cross-section.\n",
+ "\n",
+ "Remember that we often call this the **visible cross-section**, where cross-section is $\\sigma$ and the visible cross-section is $A\\times\\varepsilon\\times\\sigma$. This is called a _model independent_ limit as we do not have a specific model in mind (hence the _test_ signal). We don't have specific values for selection acceptance $A$, detector efficiency $\\varepsilon$, or signal model cross-section $\\sigma$, so we can only say something about the product. If we have a specific signal model with $A\\varepsilon\\sigma$ that gives us a cross-section higher than the limit on the visible cross-section, we should be able to see that signal in the analysis."
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": null,
+ "id": "1f5b5e3a",
+ "metadata": {},
+ "outputs": [],
+ "source": [
+ "print(f\"Limit on visible cross-section: {obs_limit / 140:0.2f} ifb\") # 140 ifb"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "id": "13e055f6",
+ "metadata": {},
+ "source": [
+ "## Discovery $p$-values\n",
+ "\n",
+ "The significance of an excess can be quantified by the probability ($p_0$) that a background-only experiment is more signal-like than observed. The last two columns indicate the $\\mathrm{CL}_\\mathrm{b}$ value and the discovery $p$-value ($p(s)=0$) with the corresponding gaussian significance ($Z$$\\,$).\n",
+ "\n",
+ "These two numbers provide information about the compatibility of the observed data with a background-only hypothesis, albeit with different test statistics $q$. The $q_\\mu$ (or $\\tilde{q}_\\mu$) test statistic is optimized for hypothesis-testing the signal+background hypothesis. The $q_0$ test statistic is optimized for hypothesis-testing the background-only hypothesis.\n",
+ "\n",
+ "From Equation 52 in [arXiv:1503.07622](https://arxiv.org/abs/1503.07622):\n",
+ "- for $p_\\mu$: the null hypothesis is signal at a rate $\\mu' < \\mu$ and the alternative hypothesis is signal at a rate of $\\mu$\n",
+ "- for $p_0$: the null hypothesis is background-only (b) and alternative hypothesis is signal+background (s+b)\n",
+ "\n",
+ "For these hypotheses, we find that the $p$-values are\n",
+ "\n",
+ "$$\n",
+ "1 - p_\\mathrm{b} = P(q \\geq q_\\mathrm{obs} | \\mathrm{b}) = \\int_{q_\\mathrm{obs}}^\\infty f(q | \\mathrm{b} )\\ \\mathrm{d}q\n",
+ "$$\n",
+ "\n",
+ "and\n",
+ "\n",
+ "$$\n",
+ "p_\\mathrm{s+b} = P(q \\geq q_\\mathrm{obs} | \\mathrm{s+b}) = \\int_{q_\\mathrm{obs}}^\\infty f(q | \\mathrm{s+b} )\\ \\mathrm{d}q\n",
+ "$$\n",
+ "\n",
+ "Note that the same $q_\\mathrm{obs}$ value is used to determine $p_\\mu$ and $p_b$."
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "id": "b092d468",
+ "metadata": {},
+ "source": [
+ "### $p_\\mu$ (and $\\mathrm{CL}_\\mathrm{b}$)\n",
+ "\n",
+ "$\\mathrm{CL}_\\mathrm{b}$ provides a measure of compatibility of the observed data with the 95% CL signal strength hypothesis relative to fluctuations of the background. $\\mathrm{CL}_\\mathrm{b}$ is evaluating $p_\\mu$ in the case where $\\mu'=0$ and $\\mu=\\mu_{95}$ (the fitted $\\mu$ at the upper limit $S^{95}_\\mathrm{obs}$). This is saying that the tested $\\mu$ is the 95% CL signal strength ($\\mu_{95}$) but the assumed true $\\mu$ is the background-only hypothesis ($\\mu'=0$):\n",
+ "\n",
+ "$$\n",
+ "\\mathrm{CL}_\\mathrm{b} \\equiv 1 - p_b = \\int_{q_{\\mu_{95}}^{\\ \\mathrm{obs}}}^\\infty f(q_\\mu | \\mu' = 0)\\ \\mathrm{d} q_\\mu\n",
+ "$$\n",
+ "\n",
+ "Note that this calculation depends on the value of $\\mu$ at the fitted upper limit $S^{95}_\\mathrm{obs}$. You get a small $\\mathrm{CL}_\\mathrm{b}$ if you have a downward fluctuation in observed data relative to your signal+background hypothesis, preventing you from excluding your signal erroneously if you \"shouldn't be able to\" when using $\\mathrm{CL}_\\mathrm{s}$ (since you would be dividing by a small number). See the paragraph above plot on page two in [this document](https://twiki.cern.ch/twiki/pub/Sandbox/NotesOnStatistics/CLsInfo.pdf).\n",
+ "\n",
+ "### $p_0$\n",
+ "\n",
+ "$p(s)=0$ measures compatibility of the observed data with the background-only (zero signal strength) hypothesis relative to fluctuations of the background. Larger values indicate greater relative compatibility. $p(s)=0$ is not calculated in signal regions with a deficit with respect to the nominal background prediction.\n",
+ "\n",
+ "This is saying that the tested $\\mu$ is the background-only hypothesis ($\\mu=\\mu'=0$) and the assumed true $\\mu$ is the background-only hypothesis:\n",
+ "\n",
+ "$$\n",
+ "p_0 = P(q_0 \\geq q_0^\\mathrm{obs}) = \\int_{q_{\\mu_{95}}^{\\ \\mathrm{obs}}}^\\infty f(q_{\\mu=0} | \\mu'=0)\\ \\mathrm{d} q_{\\mu=0}\n",
+ "$$"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "id": "afa48732",
+ "metadata": {},
+ "source": [
+ "### Relationship to $\\mathrm{CL}_\\mathrm{s}$\n",
+ "\n",
+ "This will be covered through a learn notebook, which will be added as part of [scikit-hep/pyhf#2274](https://github.com/scikit-hep/pyhf/issues/2274)."
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "id": "005b6d1c",
+ "metadata": {},
+ "source": [
+ "### Calculating the $p$-values\n",
+ "\n",
+ "Armed with all of this knowledge, we can first calculate the discovery $p$-value evaluated at $\\mu=\\mu'=0$, as well as convert this to a Z-score assuming a \"standard\" Normal distribution:"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": null,
+ "id": "11011192",
+ "metadata": {},
+ "outputs": [],
+ "source": [
+ "p0 = pyhf.infer.hypotest(0.0, data, model, test_stat=\"q0\")\n",
+ "p0_Z = scipy.stats.norm.isf(p0)\n",
+ "print(f\"p(s)=0 (Z): {p0:0.2f} ({p0_Z:0.2f})\")"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "id": "1f10a3ca",
+ "metadata": {},
+ "source": [
+ "And similarly, we can extract $\\mathrm{CL}_\\mathrm{b}$ by evaluating at $\\mu=\\mu_{95}^\\mathrm{obs}$ :"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": null,
+ "id": "74072887",
+ "metadata": {},
+ "outputs": [],
+ "source": [
+ "_, (_, CLb) = pyhf.infer.hypotest(obs_limit, data, model, return_tail_probs=True)\n",
+ "print(f\"CLb: {CLb:0.2f}\")"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "id": "2586d091",
+ "metadata": {},
+ "source": [
+ "## Making the table\n",
+ "\n",
+ "Now that we have all of our numbers, we can make the table for the discovery region `DRInt_cuts`. We'll get the corresponding index of the channel so we can extract out the total background and uncertainty estimates from above."
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": null,
+ "id": "84be7d7b",
+ "metadata": {},
+ "outputs": [],
+ "source": [
+ "region_mask = model.config.channels.index(\"DRInt_cuts\")\n",
+ "bkg_count = model.expected_actualdata(pars_bonly)[region_mask]\n",
+ "bkg_uncrt = np.sqrt(error_sq)[region_mask]"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": null,
+ "id": "0d7ee562",
+ "metadata": {},
+ "outputs": [],
+ "source": [
+ "table = Table(\n",
+ " title=\"Model-independent upper limits on the observed visible cross-section in the five electroweak search discovery regions; derived using asymptotics, pyhf, and published likelihoods.\",\n",
+ " box=box.HORIZONTALS,\n",
+ ")\n",
+ "table.add_column(\"Signal Region\", justify=\"left\", style=\"cyan\", no_wrap=True)\n",
+ "table.add_column(\"Total Bkg.\", justify=\"center\")\n",
+ "table.add_column(\"Data\")\n",
+ "table.add_column(r\"XSec \\[fb]\")\n",
+ "table.add_column(\"S 95 Obs\")\n",
+ "table.add_column(\"S 95 Exp\")\n",
+ "table.add_column(\"CLb(S 95 Obs)\")\n",
+ "table.add_column(\"p(s=0) (Z)\")\n",
+ "\n",
+ "table.add_row(\n",
+ " model.config.channels[region_mask],\n",
+ " f\"{bkg_count:0.0f} ± {bkg_uncrt:0.0f}\",\n",
+ " f\"{data[region_mask]:0.0f}\",\n",
+ " f\"{obs_limit / 140:0.2f}\",\n",
+ " f\"{obs_limit:0.1f}\",\n",
+ " f\"{exp_limit_and_bands[2]:0.1f} + {exp_limit_and_bands[3] - exp_limit_and_bands[2]:0.1f} - {exp_limit_and_bands[2] - exp_limit_and_bands[1]:0.1f}\",\n",
+ " f\"{CLb:0.2f}\",\n",
+ " f\"{p0:0.2f} ({p0_Z:0.1f})\",\n",
+ ")\n",
+ "\n",
+ "console = Console()\n",
+ "console.print(table)"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "id": "334983cb",
+ "metadata": {},
+ "source": [
+ "
\n",
+ " \n",
+ " \n",
+ " Table 20: Model-independent upper limits on the observed visible cross-section in the five electroweak search discovery regions, derived using pseudo-experiments. Left to right: background-only model post-fit total expected background, with the combined statistical and systematic uncertainties; observed data; 95 CL upper limits on the visible cross-section ($\\langle A\\varepsilon\\sigma\\rangle_\\mathrm{obs}^{95}$) and on the number of signal events ($\\mathrm{S}_\\mathrm{obs}^{95}$). The sixth column ($\\mathrm{S}_\\mathrm{exp}^{95}$) shows the expected 95 CL upper limit on the number of signal events, given the expected number (and $\\pm1\\sigma$ excursions of the expectation) of background events. The last two columns indicate the confidence level of the background-only hypothesis ($\\mathrm{CL}_\\mathrm{b}$) and discovery $p$-value with the corresponding Gaussian significance ($Z(s=0)$). $\\mathrm{CL}_\\mathrm{b}$ provides a measure of compatibility of the observed data with the signal strength hypothesis at the 95 CL limit relative to fluctuations of the background, and $p(s=0)$ measures compatibility of the observed data with the background-only hypothesis relative to fluctuations of the background. The $p$-value is capped at 0.5.\n",
+ " \n",
+ " \n",
+ "
![](\"https://atlas.web.cern.ch/Atlas/GROUPS/PHYSICS/PAPERS/SUSY-2018-05/tab_20.png\")
\n",
+ "