Notebook

`abcd_pyhf`: Likelihood-based ABCD method for background estimation and hypothesis testing with `pyhf`¶

An introduction to the likelihood-based ABCD method and how to use the Python package abcd_pyhf to implement it.

Setup¶

There are some warnings generated by some of the packages used in this notebook that can safely be ignored in order to keep the output clean. These should disappear in future releases.

In [1]:

import warnings

warnings.filterwarnings(
    'ignore',
    r'Assigned errors must be positive\. Non-positive values are replaced by a heuristic\.',
    UserWarning,
    'iminuit'
)
warnings.filterwarnings(
    'ignore',
    'invalid value encountered in double_scalars',
    RuntimeWarning,
    'pyhf'
)

Introduction¶

In analyses that are searching for new physics processes, the general idea is usually to look for a potential excess of signal-like events beyond the expected number of events from background processes. The expectations of what signal events might look like and how the background events are distributed are often provided by Monte Carlo (MC) simulation. However, it's often impractical or impossible to generate MC simulations that accurately model the background in the most signal-like regions and/or provide enough statistics to be useful for very low-background analyses. This is especially common for QCD multijet backgrounds, for example. This necessitates the use of data-driven background models, where the background in the signal region is estimated directly from data outside the signal region. Without making any further assumptions about the background, this requires at least two independent discriminant variables.

We'll look at an artificial example of signal MC and data to understand this. (This example was randomly generated in the generate_example_data.ipynb notebook.) First, let's load the "signal MC":

In [2]:

import numpy as np

In [3]:

signal_mc = np.loadtxt('signal_mc.csv', delimiter=',')

This array consists of pairs of independent observables that I'll call $x$ and $y$. Each row represents one event. We can visualize this distribution in a 2D histogram:

In [4]:

import matplotlib.pyplot as plt

In [5]:

plt.hist2d(
    *signal_mc.T,
    bins=30,
    range=((0, 3), (0, 3))
)
plt.title('Signal Monte Carlo simulation')
plt.xlabel('$x$')
plt.xlim(0, 3)
plt.ylabel('$y$')
plt.ylim(0, 3)
plt.colorbar(label='Events')
plt.show()

As in a real analysis, we have a relatively high statistics sample of what this potential signal should look like in data if it exists. Of course we don't expect to see nearly this many signal events in real data, and we will have to deal with a distribution of background events in this plane.

Next, let's look at what our distribution of events in actual "data" looks like:

In [6]:

data = np.loadtxt('data.csv', delimiter=',')

In [7]:

plt.scatter(*data.T, s=10)
plt.title('Data')
plt.xlabel('$x$')
plt.xlim(0, 3)
plt.ylabel('$y$')
plt.ylim(0, 3)
plt.show()

Now the question is: how do we know if there are real signal events here or not? In order to answer, we need a method to estimate how many of these events are from background processes.

ABCD method¶

One of the simplest possible ways of estimating background from data is called the "ABCD method". In this method, we introduce threshold values in each of the two discriminants, separating this plane into four regions labeled $A$, $B$, $C$, and $D$:

In [8]:

x_threshold = 1
y_threshold = 1

In [9]:

plt.scatter(*data.T, s=10)
plt.plot((0, 3), (y_threshold, y_threshold), color='black')
plt.plot((x_threshold, x_threshold), (0, 3), color='black')
plt.title('Data')
plt.xlabel('$x$')
plt.xlim(0, 3)
plt.ylabel('$y$')
plt.ylim(0, 3)
region_label_color = 'red'
region_label_fontsize = 20
region_label_fontweight = 'bold'
plt.text(2.95, 2.95, 'A', color=region_label_color, fontsize=region_label_fontsize, fontweight=region_label_fontweight, horizontalalignment='right', verticalalignment='top')
plt.text(0.05, 2.95, 'B', color=region_label_color, fontsize=region_label_fontsize, fontweight=region_label_fontweight, verticalalignment='top')
plt.text(2.95, 0.05, 'C', color=region_label_color, fontsize=region_label_fontsize, fontweight=region_label_fontweight, horizontalalignment='right')
plt.text(0.05, 0.05, 'D', color=region_label_color, fontsize=region_label_fontsize, fontweight=region_label_fontweight)
plt.show()

As long as $x$ and $y$ are statistically independent, it's easy to show that the number of events in region $A$, $N_A$ is approximately equal to $N_B N_C / N_D$, where the approximation is only due to finite statistics. Thus if $N_B$, $N_C$, and $N_D$ accurately reflect the number of background events in their respective regions, we have an unbiased estimator of the number of background events in region $A$ that does not rely on any information from region $A$.

Now let's look at the distribution of the simulated signal events:

In [10]:

plt.plot((0, 3), (y_threshold, y_threshold), color='black')
plt.plot((x_threshold, x_threshold), (0, 3), color='black')
plt.hist2d(
    *signal_mc.T,
    bins=30,
    range=((0, 3), (0, 3))
)

plt.title('Signal Monte Carlo simulaton')
plt.xlabel('$x$')
plt.xlim(0, 3)
plt.ylabel('$y$')
plt.ylim(0, 3)
region_label_color = 'red'
region_label_fontsize = 20
region_label_fontweight = 'bold'
plt.text(2.95, 2.95, 'A', color=region_label_color, fontsize=region_label_fontsize, fontweight=region_label_fontweight, horizontalalignment='right', verticalalignment='top')
plt.text(0.05, 2.95, 'B', color=region_label_color, fontsize=region_label_fontsize, fontweight=region_label_fontweight, verticalalignment='top')
plt.text(2.95, 0.05, 'C', color=region_label_color, fontsize=region_label_fontsize, fontweight=region_label_fontweight, horizontalalignment='right')
plt.text(0.05, 0.05, 'D', color=region_label_color, fontsize=region_label_fontsize, fontweight=region_label_fontweight)
plt.colorbar(label='Events')
plt.show()

There is a significant amount leakage of signal events into regions $B$, $C$, and $D$. So if there really are any signal events observed in data in region $A$, there very likely are signal events in data in the other regions as well. This unfortunately breaks the naive ABCD method, and if we tried to use it, we could end up totally missing our observation of signal in data.

Likelihood-based¶

The correct and much more robust way to handle this scenario is to use the "likelihood-based" ABCD method. In this method, we set up a system of equations.

Expected number of signal events in each region:

$n_X^\textrm{signal} = (\epsilon_X / \epsilon_A) \mu$

where $\epsilon_X$ is the signal efficiency of region $X$ ($X \in \{A, B, C, D\}$) estimated from MC and $\mu$ is the signal strength, defined as the number of signal events in data in region $A$.

Expected number of background events in each region:

$n_A^\textrm{bkg} = \mu_b$

$n_B^\textrm{bkg} = \tau_B \mu_b$

$n_C^\textrm{bkg} = \tau_C \mu_b$

$n_D^\textrm{bkg} = \tau_B \tau_C \mu_b$

where $\mu_b$ is defined to be the number of background events in region $A$ and $\tau_B$ and $\tau_C$ are nuisance parameters that enforce the standard ABCD relation $n_A^\textrm{bkg} = n_B^\textrm{bkg} n_C^\textrm{bkg} / n_D^\textrm{bkg}$.

Then the total expected number of events in region $X$ is $n_X = n_X^\textrm{signal} + n_X^\textrm{bkg}$. For particular values of all these parameters, we can calculate a Poisson likelihood based on the expected number of events and the observed number of events in each region. At this point, we can perform maximum likelihood fits with these parameters in order to do all statistical tests needed to get our analysis results.

`abcd_pyhf`¶

abcd_pyhf is a Python package that facilitates doing exactly these statistical tests. It's a convenient wrapper around pyhf, which is used to build the proper PDFs and run hypothesis tests for the likelihood-based ABCD method.

First we need to extract the yields (the number of events in each region for signal MC and data). I'll use a simple function for this counting:

In [11]:

def get_region_counts(x, y, x_cut, y_cut):
    return {
        'A': sum((x > x_cut) & (y > y_cut)),
        'B': sum((x > 0) & (x < x_cut) & (y > y_cut)),
        'C': sum((x > x_cut) & (y > 0) & (y < y_cut)),
        'D': sum((x > 0) & (x < x_cut) & (y > 0) & (y < y_cut))
    }

and then evaluate these yields in data and signal MC:

In [12]:

data_yields = get_region_counts(*data.T, x_threshold, y_threshold)
signal_mc_yields = get_region_counts(*signal_mc.T, x_threshold, y_threshold)

We also need the total systematic uncertainty for this analysis. This is a complicated number that would normally be the product of many dedicated studies. In this example, I'll just arbitrarily set it to 10%:

In [13]:

signal_uncertainty = 0.1

`ABCD`¶

The main class used in abcd_pyhf is called ABCD, which is an object that carries the information about a particular ABCD plane, including the yields from signal MC and data, as well as the systematic uncertainty. I'll pass in this information to make a new abcd object:

In [14]:

from abcd_pyhf import ABCD

In [15]:

abcd = ABCD(data_yields, signal_mc_yields, signal_uncertainty)

Parameter fitting¶

One of the lower-level tasks that abcd_pyhf can do is extracting the maximum likelihood estimators (MLEs) of the fit parameters and their uncertainty based on the observed data.

We can inspect the parameters used in fitting:

In [16]:

abcd.model.config.par_names()

Out[16]:

['mu', 'systematic_uncertainty', 'mu_b', 'tau_B', 'tau_C']

and we can run the fit with the signal strength fixed to zero (a background-only fit):

In [17]:

abcd.bkg_only_fit()

Out[17]:

array([[ 0.        ,  0.1       ],
       [ 0.        ,  0.1       ],
       [27.71897769,  3.56424424],
       [ 1.74226533,  0.1602753 ],
       [ 5.71076759,  0.71035935]])

or fit while allowing all parameters to float, so that the signal strength is at its MLE:

In [18]:

abcd.fit()

Out[18]:

array([[1.27676441e+01, 9.22126277e+00],
       [5.84459165e-05, 9.93343453e-01],
       [2.02291860e+01, 5.70029631e+00],
       [1.85885127e+00, 1.94267090e-01],
       [7.44703046e+00, 1.88646208e+00]])

The first column represents the parameter values, and the second column is the uncertainty. The first row is the signal strength, so $\mu = 13 \pm 9$ based on the fit.

Likelihood values¶

We can also extract the likelihood at various parameter values.

twice_nll() gets the likelihood at several values of $\mu$:

In [19]:

abcd.twice_nll()

Out[19]:

(array([ 0.,  1.,  2.,  3.,  4.,  5.,  6.,  7.,  8.,  9., 10., 11., 12.,
        13., 14., 15., 16., 17., 18., 19., 20., 21., 22., 23., 24., 25.,
        26., 27., 28., 29., 30., 31., 32., 33., 34., 35., 36., 37., 38.,
        39., 40., 41., 42., 43., 44., 45., 46., 47., 48., 49., 50., 51.,
        52., 53., 54., 55., 56., 57., 58., 59., 60., 61., 62.]),
 array([[1.82731843e+00],
        [1.55778510e+00],
        [1.30816135e+00],
        [1.07911913e+00],
        [8.71050942e-01],
        [6.84569046e-01],
        [5.20044362e-01],
        [3.77801183e-01],
        [2.58096057e-01],
        [1.61125297e-01],
        [8.68509051e-02],
        [3.54122377e-02],
        [6.71461821e-03],
        [6.59226440e-04],
        [1.68302247e-02],
        [5.51749899e-02],
        [1.15282598e-01],
        [1.96788824e-01],
        [2.99101368e-01],
        [4.21955519e-01],
        [5.64640277e-01],
        [7.26733245e-01],
        [9.07619960e-01],
        [1.10673018e+00],
        [1.32334152e+00],
        [1.55698571e+00],
        [1.80700479e+00],
        [2.07277877e+00],
        [2.35378082e+00],
        [2.64929129e+00],
        [2.95888915e+00],
        [3.28199469e+00],
        [3.61810438e+00],
        [3.96669413e+00],
        [4.32732254e+00],
        [4.69954877e+00],
        [5.08293958e+00],
        [5.47708464e+00],
        [5.88158300e+00],
        [6.29604812e+00],
        [6.72007981e+00],
        [7.15328930e+00],
        [7.59543012e+00],
        [8.04561382e+00],
        [8.50386620e+00],
        [8.96977453e+00],
        [9.44285711e+00],
        [9.92299959e+00],
        [1.04098231e+01],
        [1.09031155e+01],
        [1.14023413e+01],
        [1.19075053e+01],
        [1.24182998e+01],
        [1.29344730e+01],
        [1.34556696e+01],
        [1.39818532e+01],
        [1.45126826e+01],
        [1.50479855e+01],
        [1.55876448e+01],
        [1.61313433e+01],
        [1.66788970e+01],
        [1.72302603e+01],
        [1.77851699e+01]]))

and there's a convenience function (twice_nll_plot()) to plot this result:

In [20]:

abcd.twice_nll_plot()
plt.ylim(0, 4)
plt.show()

While the favored value of $\mu$ is around 13, we certainly can't exclude a signal-free hypothesis ($\mu = 0$)

$\textrm{CL}_s$ tests¶

We can easily extract $\textrm{CL}_b$, $\textrm{CL}_{s + b}$, and $\textrm{CL}_s$, which are used by LHC experiments for discovery and setting upper limits on signal strength [arXiv:1007.1727]:

In [21]:

abcd.clb()

Out[21]:

(array([ 0.,  1.,  2.,  3.,  4.,  5.,  6.,  7.,  8.,  9., 10., 11., 12.,
        13., 14., 15., 16., 17., 18., 19., 20., 21., 22., 23., 24., 25.,
        26., 27., 28., 29., 30., 31., 32., 33., 34., 35., 36., 37., 38.,
        39., 40., 41., 42., 43., 44., 45., 46., 47., 48., 49., 50., 51.,
        52., 53., 54., 55., 56., 57., 58., 59., 60., 61., 62.]),
 array([       nan, 0.54049462, 0.58114728, 0.62097115, 0.65943358,
        0.69637113, 0.7313514 , 0.76412622, 0.79449952, 0.82234034,
        0.84758169, 0.87022533, 0.89032565, 0.90368121, 0.90291818,
        0.90183495, 0.90070848, 0.89953496, 0.89834632, 0.89710649,
        0.89585016, 0.8945602 , 0.8932525 , 0.89193117, 0.89060104,
        0.88926322, 0.88792381, 0.8865893 , 0.88526165, 0.88393826,
        0.88262855, 0.88133203, 0.88004809, 0.87877957, 0.87752291,
        0.87627627, 0.8750297 , 0.87379191, 0.87255812, 0.87132787,
        0.87010166, 0.86888064, 0.86766039, 0.86646292, 0.86527101,
        0.86408463, 0.86291845, 0.8617614 , 0.86061656, 0.85948122,
        0.85836446, 0.85725737, 0.85616173, 0.85507638, 0.85400926,
        0.85294908, 0.8519015 , 0.85086892, 0.84984534, 0.84882924,
        0.84783104, 0.8468414 , 0.84586422]))

In [22]:

abcd.clsb()

Out[22]:

(array([ 0.,  1.,  2.,  3.,  4.,  5.,  6.,  7.,  8.,  9., 10., 11., 12.,
        13., 14., 15., 16., 17., 18., 19., 20., 21., 22., 23., 24., 25.,
        26., 27., 28., 29., 30., 31., 32., 33., 34., 35., 36., 37., 38.,
        39., 40., 41., 42., 43., 44., 45., 46., 47., 48., 49., 50., 51.,
        52., 53., 54., 55., 56., 57., 58., 59., 60., 61., 62.]),
 array([           nan, 5.00000000e-01, 5.00000000e-01, 5.00000000e-01,
        5.00000000e-01, 5.00000000e-01, 5.00000000e-01, 5.00000000e-01,
        5.00000000e-01, 5.00000000e-01, 5.00000000e-01, 5.00000000e-01,
        5.00000000e-01, 4.89758121e-01, 4.48389487e-01, 4.07145670e-01,
        3.67104154e-01, 3.28662698e-01, 2.92223166e-01, 2.57981444e-01,
        2.26198312e-01, 1.96972062e-01, 1.70373351e-01, 1.46396920e-01,
        1.24996615e-01, 1.06053558e-01, 8.94339759e-02, 7.49740490e-02,
        6.24897807e-02, 5.17980363e-02, 4.27032993e-02, 3.50218759e-02,
        2.85769803e-02, 2.32044003e-02, 1.87525990e-02, 1.50852687e-02,
        1.20811854e-02, 9.63368370e-03, 7.65002996e-03, 6.05037220e-03,
        4.76666803e-03, 3.74133149e-03, 2.92582167e-03, 2.28069801e-03,
        1.77196333e-03, 1.37241234e-03, 1.05986510e-03, 8.16130761e-04,
        6.26733742e-04, 4.80013294e-04, 3.66758007e-04, 2.79541115e-04,
        2.12573566e-04, 1.61293018e-04, 1.22133450e-04, 9.22917917e-05,
        6.96095556e-05, 5.24060320e-05, 3.93838006e-05, 2.95489210e-05,
        2.21353700e-05, 1.65559993e-05, 1.23651676e-05]))

In the case of $\textrm{CL}_s$, we get the observed values as well as the expected band, which is the median value (assuming the true value of $\mu$ is zero), $\pm 1 \sigma$, and $\pm 2 \sigma$.

In [23]:

abcd.cls()

Out[23]:

(array([ 0.,  1.,  2.,  3.,  4.,  5.,  6.,  7.,  8.,  9., 10., 11., 12.,
        13., 14., 15., 16., 17., 18., 19., 20., 21., 22., 23., 24., 25.,
        26., 27., 28., 29., 30., 31., 32., 33., 34., 35., 36., 37., 38.,
        39., 40., 41., 42., 43., 44., 45., 46., 47., 48., 49., 50., 51.,
        52., 53., 54., 55., 56., 57., 58., 59., 60., 61., 62.]),
 array([           nan, 9.25078580e-01, 8.60367105e-01, 8.05190390e-01,
        7.58226473e-01, 7.18007935e-01, 6.83665883e-01, 6.54342162e-01,
        6.29326997e-01, 6.08020759e-01, 5.89913641e-01, 5.74563829e-01,
        5.61592268e-01, 5.41958951e-01, 4.96600353e-01, 4.51463616e-01,
        4.07572664e-01, 3.65369565e-01, 3.25290103e-01, 2.87570591e-01,
        2.52495697e-01, 2.20188716e-01, 1.90733698e-01, 1.64134772e-01,
        1.40350853e-01, 1.19260030e-01, 1.00722579e-01, 8.45645769e-02,
        7.05890519e-02, 5.85991561e-02, 4.83819603e-02, 3.97374369e-02,
        3.24720669e-02, 2.64052568e-02, 2.13699252e-02, 1.72151971e-02,
        1.38066004e-02, 1.10251464e-02, 8.76735863e-03, 6.94385251e-03,
        5.47828865e-03, 4.30592113e-03, 3.37208163e-03, 2.63219345e-03,
        2.04787091e-03, 1.58828464e-03, 1.22823321e-03, 9.47049566e-04,
        7.28238065e-04, 5.58491892e-04, 4.27275387e-04, 3.26087737e-04,
        2.48286695e-04, 1.88629954e-04, 1.43011856e-04, 1.08203167e-04,
        8.17108028e-05, 6.15911932e-05, 4.63423152e-05, 3.48113845e-05,
        2.61082327e-05, 1.95502951e-05, 1.46183836e-05]),
 array([[1.00000000e+00, 7.82002601e-01, 6.03646409e-01, 4.61477246e-01,
         3.49721834e-01, 2.62357190e-01, 1.95011965e-01, 1.43632814e-01,
         1.04846253e-01, 7.58696200e-02, 5.44429256e-02, 3.87526636e-02,
         2.73723539e-02, 1.91905906e-02, 1.33590602e-02, 9.23822541e-03,
         6.34752329e-03, 4.33509799e-03, 2.94386618e-03, 1.98840777e-03,
         1.33627902e-03, 8.93853145e-04, 5.95223759e-04, 3.94677700e-04,
         2.60716492e-04, 1.71579728e-04, 1.12527021e-04, 7.35552460e-05,
         4.79277723e-05, 3.11469603e-05, 2.01835240e-05, 1.30452897e-05,
         8.41133259e-06, 5.41111817e-06, 3.47382101e-06, 2.22587740e-06,
         1.42410376e-06, 9.09656865e-07, 5.80251639e-07, 3.69668640e-07,
         2.35245159e-07, 1.49549415e-07, 9.49843504e-08, 6.02773000e-08,
         3.82266466e-08, 2.42309406e-08, 1.53495114e-08, 9.71907498e-09,
         6.15168582e-09, 3.89257817e-09, 2.46252656e-09, 1.55758285e-09,
         9.85090247e-10, 6.23016578e-10, 3.93996304e-10, 2.49198496e-10,
         1.57636934e-10, 9.97248240e-11, 6.31025012e-11, 3.99462084e-11,
         2.52910859e-11, 1.60181109e-11, 1.01485850e-11],
        [1.00000000e+00, 8.52795126e-01, 7.19387406e-01, 6.01500194e-01,
         4.98782292e-01, 4.09781332e-01, 3.33738484e-01, 2.69446367e-01,
         2.15671205e-01, 1.71171086e-01, 1.34735016e-01, 1.05203370e-01,
         8.15071079e-02, 6.26705503e-02, 4.78349028e-02, 3.62581041e-02,
         2.72967329e-02, 2.04175624e-02, 1.51777229e-02, 1.12160801e-02,
         8.24186149e-03, 6.02431200e-03, 4.38081428e-03, 3.17002559e-03,
         2.28357862e-03, 1.63770650e-03, 1.16959310e-03, 8.31926448e-04,
         5.89442716e-04, 4.16208665e-04, 2.92836758e-04, 2.05351498e-04,
         1.43550829e-04, 1.00048053e-04, 6.95329913e-05, 4.81977774e-05,
         3.33328490e-05, 2.29983991e-05, 1.58344621e-05, 1.08804477e-05,
         7.46250675e-06, 5.10934438e-06, 3.49251698e-06, 2.38364350e-06,
         1.62460502e-06, 1.10595191e-06, 7.51892918e-07, 5.10609374e-07,
         3.46394769e-07, 2.34769752e-07, 1.58976196e-07, 1.07565182e-07,
         7.27265071e-08, 4.91406033e-08, 3.31815653e-08, 2.23948252e-08,
         1.51077116e-08, 1.01867051e-08, 6.86618391e-09, 4.62728037e-09,
         3.11720265e-09, 2.09949968e-09, 1.41378230e-09],
        [1.00000000e+00, 9.19010750e-01, 8.37705447e-01, 7.58057706e-01,
         6.81132833e-01, 6.07257732e-01, 5.37297202e-01, 4.71747568e-01,
         4.11000950e-01, 3.55319314e-01, 3.04836622e-01, 2.59549333e-01,
         2.19348704e-01, 1.84015815e-01, 1.53265551e-01, 1.26768049e-01,
         1.04134453e-01, 8.49762484e-02, 6.88982057e-02, 5.55159399e-02,
         4.44650510e-02, 3.54101293e-02, 2.80416446e-02, 2.20865441e-02,
         1.73080986e-02, 1.34957692e-02, 1.04729611e-02, 8.08971692e-03,
         6.22077340e-03, 4.76402579e-03, 3.63315949e-03, 2.75978806e-03,
         2.08843798e-03, 1.57464143e-03, 1.18313700e-03, 8.86039890e-04,
         6.61564843e-04, 4.92472643e-04, 3.65573813e-04, 2.70650888e-04,
         1.99867837e-04, 1.47240018e-04, 1.08220749e-04, 7.93665815e-05,
         5.80869446e-05, 4.24333860e-05, 3.09381344e-05, 2.25172133e-05,
         1.63609726e-05, 1.18691248e-05, 8.59761529e-06, 6.21900842e-06,
         4.49244750e-06, 3.24122847e-06, 2.33557595e-06, 1.68118481e-06,
         1.20888931e-06, 8.68362315e-07, 6.23186969e-07, 4.46904551e-07,
         3.20193704e-07, 2.29239039e-07, 1.64003133e-07],
        [1.00000000e+00, 9.69272753e-01, 9.35101669e-01, 8.97992383e-01,
         8.58265232e-01, 8.15978624e-01, 7.71592423e-01, 7.25505087e-01,
         6.78183483e-01, 6.30137453e-01, 5.81904346e-01, 5.34013405e-01,
         4.86983915e-01, 4.41281729e-01, 3.97330926e-01, 3.55510179e-01,
         3.16091877e-01, 2.79301584e-01, 2.45284087e-01, 2.14114427e-01,
         1.85802878e-01, 1.60309233e-01, 1.37531155e-01, 1.17336631e-01,
         9.95768353e-02, 8.40620496e-02, 7.06042029e-02, 5.90071251e-02,
         4.90758788e-02, 4.06304176e-02, 3.34839920e-02, 2.74730810e-02,
         2.24452947e-02, 1.82618685e-02, 1.47991219e-02, 1.19470767e-02,
         9.61025980e-03, 7.70301550e-03, 6.15346691e-03, 4.89968509e-03,
         3.88920814e-03, 3.07786085e-03, 2.42875105e-03, 1.91119853e-03,
         1.49997329e-03, 1.17431649e-03, 9.17058742e-04, 7.14479582e-04,
         5.55398357e-04, 4.30808011e-04, 3.33475837e-04, 2.57621916e-04,
         1.98644235e-04, 1.52893681e-04, 1.17469848e-04, 9.01070503e-05,
         6.90088728e-05, 5.27676354e-05, 4.02905135e-05, 3.07241521e-05,
         2.33960785e-05, 1.77935022e-05, 1.35160562e-05],
        [1.00000000e+00, 9.93781700e-01, 9.86121127e-01, 9.76898497e-01,
         9.65954682e-01, 9.53043388e-01, 9.38022216e-01, 9.20737885e-01,
         9.01074662e-01, 8.78961103e-01, 8.54380038e-01, 8.27366603e-01,
         7.98020316e-01, 7.66489437e-01, 7.32983521e-01, 6.97778049e-01,
         6.61162239e-01, 6.23480358e-01, 5.85093136e-01, 5.46372385e-01,
         5.07688209e-01, 4.69407827e-01, 4.31854622e-01, 3.95331896e-01,
         3.60129845e-01, 3.26457461e-01, 2.94504778e-01, 2.64410114e-01,
         2.36267361e-01, 2.10157072e-01, 1.86074084e-01, 1.64013344e-01,
         1.43934545e-01, 1.25770984e-01, 1.09439428e-01, 9.48401228e-02,
         8.18685134e-02, 7.03976309e-02, 6.03084546e-02, 5.14782444e-02,
         4.37866363e-02, 3.71172277e-02, 3.13596217e-02, 2.64098997e-02,
         2.21727713e-02, 1.85605085e-02, 1.54908737e-02, 1.28925435e-02,
         1.07008343e-02, 8.85836291e-03, 7.31446671e-03, 6.02478250e-03,
         4.95069396e-03, 4.05882168e-03, 3.32011902e-03, 2.71013203e-03,
         2.20766218e-03, 1.79469293e-03, 1.45617616e-03, 1.17941203e-03,
         9.53474473e-04, 7.69494523e-04, 6.19974382e-04]]))

Brazil plot¶

There's another convenience function for plotting the $\textrm{CL}_s$ band, sometimes called a "Brazil plot":

In [24]:

abcd.brazil_plot()
plt.xlim(0)
plt.show()

Upper limits¶

And finally we can extract the observed and expected upper limit band for the signal strength (where the lines in the plot above cross 0.05, for a 95% confidence level):

In [25]:

abcd.upper_limit()

Out[25]:

(29.84163563742085,
 [10.283164525558307,
  13.854061161638194,
  19.499139929855087,
  27.90694811529501,
  39.19218925297826])

Summary and outlook¶

The ABCD method
- One of the simplest and most common data-driven background estimation methods
- Really a family of many methods using very similar techniques, but almost all that are actually used in modern analyses are some variant of the likelihood-based method described here
abcd_pyhf
- A small, pure-Python package
- Depends only on NumPy, Matplotlib, and pyhf
  - No dependence on experiment-specific formats or code
- Works out of the box for most simple cases
  - Intended to be extensible to accommodate more complicated cases
- Using this for my own ATLAS analysis, and a few other ATLAS analyses now also using or trying out this package

abcd_pyhf: Likelihood-based ABCD method for background estimation and hypothesis testing with pyhf¶