diff --git a/book/UpperLimitsTable.ipynb b/book/UpperLimitsTable.ipynb index 1df6cf1..bf9b40d 100644 --- a/book/UpperLimitsTable.ipynb +++ b/book/UpperLimitsTable.ipynb @@ -7,21 +7,22 @@ "metadata": {}, "outputs": [], "source": [ - "import pyhf\n", - "import pyhf.contrib.utils\n", - "from pyhf.contrib.viz import brazil\n", - "import tarfile\n", "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", "# for making tables\n", "from rich.console import Console\n", - "from rich.table import Table\n", - "from rich import box" + "from rich.table import Table" ] }, { @@ -61,7 +62,9 @@ "metadata": {}, "outputs": [], "source": [ - "pyhf.contrib.utils.download(\"https://doi.org/10.17182/hepdata.116034.v1/r34\", \"workspace_2L2J\")" + "pyhf.contrib.utils.download(\n", + " \"https://doi.org/10.17182/hepdata.116034.v1/r34\", \"workspace_2L2J\"\n", + ")" ] }, { @@ -81,7 +84,7 @@ "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", + "\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\")))" ] @@ -200,16 +203,18 @@ } ], "source": [ - "pars_uncrt_bonly, corr_bonly = pyhf.infer.mle.fixed_poi_fit(0.0, data, model, return_uncertainties=True, return_correlations=True)\n", + "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", + "for region, count in zip(model.config.channels, model.expected_actualdata(pars_bonly)):\n", " table.add_row(region, f\"{count:0.2f}\")\n", - " \n", + "\n", "console = Console()\n", "console.print(table)" ] @@ -267,7 +272,7 @@ "source": [ "npars = len(model.config.parameters)\n", "model_batch = ws.model(batch_size=npars)\n", - "pars_batch = np.tile(pars_bonly, (npars,1))" + "pars_batch = np.tile(pars_bonly, (npars, 1))" ] }, { @@ -287,7 +292,7 @@ "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" + "variations = (up_yields - dn_yields) / 2" ] }, { @@ -317,7 +322,7 @@ } ], "source": [ - "error_sq = np.diag(np.einsum('ij,ik,kl->jl', variations, corr_bonly, variations))\n", + "error_sq = np.diag(np.einsum(\"ij,ik,kl->jl\", variations, corr_bonly, variations))\n", "error_sq" ] }, @@ -374,9 +379,11 @@ "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(model.config.channels, model.expected_actualdata(pars_bonly), np.sqrt(error_sq)):\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", + "\n", "console = Console()\n", "console.print(table)" ] @@ -411,7 +418,11 @@ "pyhf.set_backend(\"numpy\", \"scipy\")\n", "\n", "mu_tests = np.linspace(15, 25, 10)\n", - "obs_limit, exp_limit_and_bands, (poi_tests, tests) = pyhf.infer.intervals.upper_limits.upper_limit(\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, return_results=True\n", ")" ] @@ -432,8 +443,10 @@ } ], "source": [ - "print(f'Observed UL: {obs_limit:0.2f}')\n", - "print(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}]')" + "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", + ")" ] }, { @@ -461,7 +474,7 @@ } ], "source": [ - "print(f'Limit on visible cross-section: {obs_limit / 140:0.2f}') # 140 ifb" + "print(f\"Limit on visible cross-section: {obs_limit / 140:0.2f}\") # 140 ifb" ] }, { @@ -566,9 +579,11 @@ } ], "source": [ - "p0 = pyhf.infer.hypotest(0.0, data, model, test_stat='q0')#, calctype='toybased', ntoys=10000)\n", + "p0 = pyhf.infer.hypotest(\n", + " 0.0, data, model, test_stat=\"q0\"\n", + ") # , calctype='toybased', ntoys=10000)\n", "p0_Z = scipy.stats.norm.isf(p0)\n", - "print(f'p(s)=0 (Z): {p0:0.2f} ({p0_Z:0.2f})')" + "print(f\"p(s)=0 (Z): {p0:0.2f} ({p0_Z:0.2f})\")" ] }, { @@ -595,7 +610,7 @@ ], "source": [ "_, (_, CLb) = pyhf.infer.hypotest(obs_limit, data, model, return_tail_probs=True)\n", - "print(f'CLb: {CLb:0.2f}')" + "print(f\"CLb: {CLb:0.2f}\")" ] }, { @@ -615,7 +630,7 @@ "metadata": {}, "outputs": [], "source": [ - "region_mask = model.config.channels.index('DRInt_cuts')\n", + "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]" ] @@ -653,7 +668,10 @@ } ], "source": [ - "table = Table(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.\", box=box.HORIZONTALS)\n", + "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", @@ -664,14 +682,14 @@ "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", + " 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", + " f\"{p0:0.2f} ({p0_Z:0.1f})\",\n", ")\n", "\n", "console = Console()\n",