<jupyter_start><jupyter_text>![cabinetry logo](_static/cabinetry_logo_small.png) 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 `JAX` backend for `pyhf` here for a demonstration below, as well as `hist` (again, not needed if running via Binder).<jupyter_code>#! pip install cabinetry[contrib]
#! pip install pyhf[jax,contrib]
#! pip install hist[plot]<jupyter_output><empty_output><jupyter_text>`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.<jupyter_code>import cabinetry
cabinetry.set_logging()<jupyter_output><empty_output><jupyter_text>We will also set up additional imports that are used later on in this notebook.<jupyter_code>import copy
import json
import pathlib

import boost_histogram as bh
import hist
import numpy as np
import pyhf
from pyhf.contrib.utils import download<jupyter_output><empty_output><jupyter_text>We will have a look at three different ways of using `cabinetry` in this notebook:- creating a statistical model,- performing statistical inference with our model,- exploring the statistical model of an ATLAS analysis.There is also bonus material at the end for you to look through.  Creating a statistical modelStatistical 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.  Defining the modelThere are a few things we need. First up, some general settings:<jupyter_code>config = {
   "General":{
      "Measurement": "PyHEP_cabinetry",
      "POI": "Signal_norm",               # parameter of interest, which we want to measure 
      "InputPath": "input/{SamplePaths}", # where to find input data
      "HistogramFolder": "histograms/"

   }
}<jupyter_output><empty_output><jupyter_text>Note the `input/{SamplePaths}` structure, which we will get back to shortly.Now it is time to think more about physics, let's define a phase space region that contains data we want to fit to.<jupyter_code>config.update({
   "Regions":[
      {
         "Name": "Signal_region",
         "Filter": "lep_charge > 0",           # event selection 
         "Variable": "jet_pt",                 # which variable we bin histograms in
         "Binning": [200, 300, 400, 500, 600]
      }
   ]
})<jupyter_output><empty_output><jupyter_text>`"Regions"` is a list, because we can use events from more than one phase space region.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.<jupyter_code>config.update({
   "Samples":[
      {
         "Name": "Data",
         "Tree": "pseudodata",
         "SamplePaths": "data.root",
         "Data": True                       # observed data is handled differently, need to distinguish
      },
      {
         "Name": "Signal",
         "Tree": "signal",
         "SamplePaths": "prediction.root",
         "Weight": "weight"                 # weights: Monte Carlo integration, simulation correction etc.
      },
      {
         "Name": "Background",
         "Tree": "background",
         "SamplePaths": "prediction.root",
         "Weight": "weight"
      }
   ]
})<jupyter_output><empty_output><jupyter_text>`"Samples"` is again a list, allowing us to include arbitrarly many different types of processes we may need to consider. The `"SamplePaths"` option replaces the placeholder in `input/{SamplePaths}`, and specifies where each sample can be found (see the [documentation](https://cabinetry.readthedocs.io/en/latest/core.htmlinput-file-path-specification)).We are almost done defining our statistical model. What about systematic uncertainties? For now, we won't define any.<jupyter_code>config.update({"Systematics": []})<jupyter_output><empty_output><jupyter_text>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!<jupyter_code>config.update({
   "NormFactors":[
      {
         "Name": "Signal_norm",
         "Samples": "Signal",    # we want this parameter to scale the signal
         "Nominal": 1,
         "Bounds": [-5, 10]
      }
   ]
})<jupyter_output><empty_output><jupyter_text>It is a good idea to validate that our configuration satisfies the `cabinetry` configuration schema.<jupyter_code>cabinetry.configuration.validate(config)<jupyter_output><empty_output><jupyter_text>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.<jupyter_code>cabinetry.configuration.print_overview(config)<jupyter_output><empty_output><jupyter_text>Creating histogramsWe 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.<jupyter_code>cabinetry.template_builder.create_histograms(config, method="uproot")<jupyter_output><empty_output><jupyter_text>`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.A post-processing step can be run to apply optional operations like histogram smoothing. See the bonus material.<jupyter_code>cabinetry.template_postprocessor.run(config)<jupyter_output><empty_output><jupyter_text>Let's visualize the histograms we have produced.<jupyter_code>cabinetry.visualize.data_MC_from_histograms(config)<jupyter_output><empty_output><jupyter_text>We see our expected distributions of signal and background, as well as our observed (pseudo-) data.  A more complex model: adding systematic uncertaintiesLet's make our model a bit more realistic and add a few systematic uncertainties:- a 5% luminosity uncertainty,- a `Modeling` uncertainty derived from comparing our nominal background prediction to that of a different simulator,- a `WeightBasedModeling` modeling uncertainty derived from varying the weights we apply to background events.<jupyter_code>config.update({
   "Systematics":[
      {
         "Name": "Luminosity",
         "Up": {"Normalization": 0.05},
         "Down": {"Normalization": -0.05},
         "Type": "Normalization"
      },
      {
         "Name":"Modeling",
         "Up": {"Tree": "background_varied"},
         "Down": {"Symmetrize": True},
         "Samples": "Background",
         "Type": "NormPlusShape"
      },
      {
         "Name": "WeightBasedModeling",
         "Up": {"Weight": "weight_up"},
         "Down": {"Weight": "0.7*weight"},
         "Samples": "Background",
         "Type": "NormPlusShape"
      }
   ],
})<jupyter_output><empty_output><jupyter_text>These new systematic uncertainties require new histograms, so let's build those.<jupyter_code>cabinetry.template_builder.create_histograms(config, method="uproot")
cabinetry.template_postprocessor.run(config)<jupyter_output><empty_output><jupyter_text>We can visualize the template histograms corresponding to systematic variations of our model.<jupyter_code>cabinetry.visualize.templates(config)<jupyter_output><empty_output><jupyter_text>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.  Building a workspaceWith 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.<jupyter_code>workspace_path = "PyHEP_workspace.json"
ws = cabinetry.workspace.build(config)
cabinetry.workspace.save(ws, workspace_path)<jupyter_output><empty_output><jupyter_text>Saving is optional, we could directly continue with the `ws` object. `pyhf` has validated our workspace, and we are now ready for statistical inference.  Performing statistical inference with our modelTo perform inference, we need two things: a probability density function (pdf), or `model`, and data to fit it to. Both are contained in the workspace.<jupyter_code>model, data = cabinetry.model_utils.model_and_data(ws)<jupyter_output><empty_output><jupyter_text>You can see our `Signal_norm` normalization showing up, in addition to parameters for the systematic uncertainties we defined.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.`data` is a list, starting with the observed counts per bin.<jupyter_code>data<jupyter_output><empty_output><jupyter_text>What comes after is so-called auxiliary data. Check out [this `pyhf` tutorial](https://pyhf.github.io/pyhf-tutorial/HelloWorld.htmlauxiliary-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.   Maximum likelihood estimate (MLE)Let's fit our model to data to obtain the maximum likelihood estimate (MLE).<jupyter_code>fit_results = cabinetry.fit.fit(model, data)<jupyter_output><empty_output><jupyter_text>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.<jupyter_code>for label, result, unc in zip(fit_results.labels, fit_results.bestfit, fit_results.uncertainty):
    print(f"{label}: {result:.3f} +/- {unc:.3f}")<jupyter_output><empty_output><jupyter_text>It is usually helpful to visualize the fit results. Let's start with the *pull plot* showing us best-fit parameter results.<jupyter_code>cabinetry.visualize.pulls(fit_results, exclude="Signal_norm")<jupyter_output><empty_output><jupyter_text>The parameter correlation matrix has a handy `pruning_threshold` setting to filter out parameters that are not highly correlated with others.<jupyter_code>cabinetry.visualize.correlation_matrix(fit_results, pruning_threshold=0.1)<jupyter_output><empty_output><jupyter_text>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.<jupyter_code>cabinetry.visualize.data_MC(model, data, config=config)<jupyter_output><empty_output><jupyter_text>This gave us not only a plot, but also yield tables with model prediction and uncertainties. We create the post-fit version of this plot by passing in `fit_results`, and this time we skip the table to fit both figures onto the same screen.<jupyter_code>cabinetry.visualize.data_MC(model, data, config=config, fit_results=fit_results, include_table=False)<jupyter_output><empty_output><jupyter_text>The contribution of signal in green has increased.  Expected sensitivity with the Asimov datasetWe can evaluate the expected performance of our model with the so-called Asimov dataset (see [arXiv:1007.1727](https://arxiv.org/abs/1007.1727)).<jupyter_code>asimov_data = cabinetry.model_utils.build_Asimov_data(model)
_ = cabinetry.fit.fit(model, asimov_data)<jupyter_output><empty_output><jupyter_text>By definition, none of the parameters are pulled away from their initial values in this fit.  Beyond MLEs: discovery significance and parameter limitsNow that we ran a basic fit, let's do something slightly more involved: calculate discovery significance.<jupyter_code>significance_results = cabinetry.fit.significance(model, data)<jupyter_output><empty_output><jupyter_text>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.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!<jupyter_code>data_limit = [112, 129, 92, 63] + model.config.auxdata  # need auxiliary data as well
limit_results = cabinetry.fit.limit(model, data_limit)
cabinetry.visualize.limit(limit_results)<jupyter_output><empty_output><jupyter_text>The observed limit follows the expectation quite closely with this new data we have used, which is very similar to the background-only prediction.  Exploring the statistical model of an ATLAS analysisLet'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).We can download it, pick one of the available signals with `pyhf`, and are ready to go.<jupyter_code>download("https://www.hepdata.net/record/resource/1267798?view=true", "electroweakinos")
ATLAS_ws = pyhf.Workspace(json.load(open("electroweakinos/1Lbb-likelihoods-hepdata/BkgOnly.json")))
patchset = pyhf.PatchSet(json.load(open("electroweakinos/1Lbb-likelihoods-hepdata/patchset.json")))
ATLAS_ws = patchset.apply(ATLAS_ws, "C1N2_Wh_hbb_700_400")
cabinetry.workspace.save(ATLAS_ws, "electroweakinos.json")<jupyter_output><empty_output><jupyter_text>`pyhf` has a helpful command line interface utility to summarize the workspace content.<jupyter_code>!pyhf inspect electroweakinos.json | head -n 6<jupyter_output><empty_output><jupyter_text>Let's extract model and data with `cabinetry`.<jupyter_code>ATLAS_model, ATLAS_data = cabinetry.model_utils.model_and_data(ATLAS_ws)<jupyter_output><empty_output><jupyter_text>We are ready to take a closer look at the content. Which phase space regions are included?<jupyter_code>cabinetry.visualize.data_MC(ATLAS_model, ATLAS_data, include_table=False)<jupyter_output><empty_output><jupyter_text>Let's fit this model to data.<jupyter_code>ATLAS_fit_results = cabinetry.fit.fit(ATLAS_model, ATLAS_data)<jupyter_output><empty_output><jupyter_text>To parse the results more easily, let's visualize them.<jupyter_code>cabinetry.visualize.pulls(ATLAS_fit_results)<jupyter_output><empty_output><jupyter_text>You can now also visualize the post-fit model distribution if you wish. It takes a few seconds to calculate, and is therefore commented out below.<jupyter_code>#cabinetry.visualize.data_MC(ATLAS_model, ATLAS_data, fit_results=ATLAS_fit_results, include_table=False)<jupyter_output><empty_output><jupyter_text>Bonus materialAdditional material that does not fit into the main talk. Check it out for more information!  Command line interface`cabinetry` also provides a [command line interface](https://cabinetry.readthedocs.io/en/latest/cli.html). Feel free to explore it!<jupyter_code>! cabinetry --help<jupyter_output><empty_output><jupyter_text>Input histograms via a custom functionIt 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.htmloverrides-for-template-building).<jupyter_code>my_router = cabinetry.route.Router()

# define a custom template builder function that is executed for data samples
@my_router.register_template_builder(sample_name="Data")
def build_data_hist(region: dict, sample: dict, systematic: dict, template: str):
    hist = bh.Histogram(
        bh.axis.Variable(region["Binning"], underflow=False, overflow=False),
        storage=bh.storage.Weight(),
    )
    yields = np.asarray([17, 12, 25, 20])
    variance = np.asarray([1.5, 1.2, 1.8, 1.6])
    hist[...] = np.stack([yields, variance], axis=-1)
    return hist  # return a boost-histogram histogram

custom_config = copy.deepcopy(config)
custom_config["General"]["HistogramFolder"] = "histograms_custom/"
cabinetry.template_builder.create_histograms(custom_config, router=my_router)<jupyter_output><empty_output><jupyter_text>The histogram creation called our function to create the data histogram.We can load the histogram and check the yields to verify they were correctly picked up.<jupyter_code>h = cabinetry.histo.Histogram.from_path(pathlib.Path("histograms_custom/Signal_region_Data_Nominal.npz"))
h.values()<jupyter_output><empty_output><jupyter_text>Histogram smoothingThis 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.<jupyter_code>smoothing_config = copy.deepcopy(config)
smoothing_config["General"]["HistogramFolder"] = "histograms_smoothing/"
smoothing_config.update({
   "Systematics":[
      {
         "Name":"Modeling",
         "Up": {"SamplePaths": "prediction.root", "Tree": "background_varied"},
         "Down": {"Symmetrize": True},
         "Samples": "Background",
         "Smoothing": {"Algorithm": "353QH, twice"},  # smoothing applied
         "Type": "NormPlusShape"
      }
   ]
})
cabinetry.template_builder.create_histograms(smoothing_config, method="uproot")
cabinetry.template_postprocessor.run(smoothing_config)
cabinetry.visualize.templates(smoothing_config)<jupyter_output><empty_output><jupyter_text>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`).  Exploring histograms with `hist``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.<jupyter_code>cabinetry_histogram = cabinetry.histo.Histogram.from_config(
    config["General"]["HistogramFolder"],
    config["Regions"][0],  # we only have one region: Signal_region
    config["Samples"][0],  # data
    {"Name": "Nominal"},   # no systematics
)
hist.Hist(cabinetry_histogram).plot()<jupyter_output><empty_output><jupyter_text>Using a different `pyhf` backend 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.<jupyter_code>pyhf.set_backend("jax")
_ = cabinetry.fit.fit(model, data)<jupyter_output><empty_output><jupyter_text>Going back and forth between `cabinetry` and `pyhf`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`.<jupyter_code>asimov_data = cabinetry.model_utils.build_Asimov_data(model)
pyhf.infer.mle.fit(asimov_data, model)<jupyter_output><empty_output><jupyter_text>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`.<jupyter_code>ws_edited = copy.deepcopy(ws)
ws_edited["measurements"][0]["config"]["parameters"].append({"name": "Modeling", "inits": [-0.5], "fixed": True})
model_edited = pyhf.Workspace(ws_edited).model()
_ = cabinetry.fit.fit(model_edited, data)<jupyter_output><empty_output><jupyter_text>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.  Further MLE options: MINOS and goodness-of-fitThe `cabinetry.fit.fit` API allows use of the [MINOS algorithm](https://iminuit.readthedocs.io/en/stable/reference.htmliminuit.Minuit.minos) to compute confidence intervals via parameter scans. All parameters for which the algorithm should be run have to be listed.<jupyter_code>_ = cabinetry.fit.fit(model, data, minos=["Signal_norm", "Modeling"])<jupyter_output><empty_output><jupyter_text>MINOS uncertainties are not restricted to be symmetric by construction.`cabinetry` also implementes a goodness-of-fit calculation using the [saturated model](http://www.physics.ucla.edu/~cousins/stats/cousins_saturated.pdf).<jupyter_code>_ = cabinetry.fit.fit(model, data, goodness_of_fit=True)<jupyter_output><empty_output><jupyter_text>Likelihood scanWe 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.<jupyter_code>scan_results = cabinetry.fit.scan(model, data, "Signal_norm", n_steps=5)<jupyter_output><empty_output><jupyter_text>The scan in this example shows relatively good agreement with a naive Gaussian approximation.<jupyter_code>cabinetry.visualize.scan(scan_results)<jupyter_output><empty_output><jupyter_text>Nuisance parameter rankingWe 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.<jupyter_code>ranking_results = cabinetry.fit.ranking(model, data, fit_results=fit_results)<jupyter_output><empty_output><jupyter_text>The figure visualizes the impact of nuisance parameters on the POI in order of decreasing impact.<jupyter_code>cabinetry.visualize.ranking(ranking_results)<jupyter_output><empty_output>
{'fraction_non_alphanumeric': 0.09071739225936827, 'fraction_numerical': 0.008566533409480296, 'mean_word_length': 5.602088167053364, 'pattern_counts': {'":': 65, '<': 197, '<?xml version=': 0, '>': 198, 'https://': 11, 'lorem ipsum': 0, 'www.': 2, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 3, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'chain_length': '49', 'dir': 'jupyter-structured-clean-dedup', 'id': '119891', 'license': 'BSD-3-Clause', 'n_tokens_mistral': 6915, 'n_tokens_neox': 6440, 'n_words': 2228, 'path': 'talk.ipynb', 'repo_name': 'alexander-held/PyHEP-2021-cabinetry-talk', 'size': '22763'}