Contents

Likelihood Serialization and Patches

Using HEPData

Preserved on HEPData

As of this tutorial, ATLAS has published 4 full likelihoods to HEPData

Let’s explore the 1Lbb workspace a little bit shall we?

Getting the Data

We’ll use python requests and tarfile to download the tarball and extract the files we need.

import requests
import hashlib
import tarfile
import json

electroweakino_HEPData_URL = "https://www.hepdata.net/record/resource/1408476?view=true"
targz_filename = "data/1Lbb_workspaces.tar.gz"
response = requests.get(electroweakino_HEPData_URL, stream=True)
assert response.status_code == 200
with open(targz_filename, "wb") as file:
    file.write(response.content)
hashlib.sha256(open(targz_filename, "rb").read()).hexdigest()
# Open as a tarfile
tar_archive = tarfile.open(targz_filename, "r:gz")

Instantiate our objects

To make this easier on us, we’ll first define a nice helper function to extract out a file from the tarball.

def get_json_from_tarfile(tarfile, json_name):
    json_file = tarfile.extractfile(tarfile.getmember(json_name)).read().decode("utf8")
    return json.loads(json_file)

Let’s go ahead and import pyhf and get started.

import pyhf

spec = get_json_from_tarfile(tar_archive, "BkgOnly.json")
patchset = pyhf.PatchSet(get_json_from_tarfile(tar_archive, "patchset.json"))

So what did the analyzers give us for signal patches?

Patching in Signals

Let’s look at this pyhf.PatchSet object which provides a user-friendly way to interact with many signal patches at once.

PatchSet

patchset

<pyhf.patchset.PatchSet object with 125 patches at 0x7f4a28f833d0>

Oh wow, we’ve got 125 patches. What information does it have?

print(f'description: {patchset.description}')
print(f'    digests: {patchset.digests}')
print(f'     labels: {patchset.labels}')
print(f' references: {patchset.references}')
print(f'    version: {patchset.version}')

description: signal patchset for the SUSY EWK 1Lbb analysis
    digests: {'sha256': '2563672e1a165384faf49f1431e921d88c9c07ec10f150d5045576564f98f820'}
     labels: ['m1', 'm2']
 references: {'hepdata': 'ins1755298'}
    version: 1.0.0

So we’ve got a useful description of the signal patches… there’s a digest. Does that match the background-only workspace we have?

pyhf.utils.digest(spec)

'2563672e1a165384faf49f1431e921d88c9c07ec10f150d5045576564f98f820'

It does! In fact, this sort of verification check will be done automatically when applying patches using pyhf.PatchSet as we will see shortly. To manually verify, simply run pyhf.PatchSet.verify on the workspace. No error means everything is fine. It will loudly complain otherwise.

patchset.verify(spec)

No error, whew. Let’s move on.

The labels m1 and m2 tells us that we have the signal patches parametrized in 2-dimensional space, likely as \(m_1 = \tilde{\chi}_1^\pm\) and \(m_2 = \tilde{\chi}_1^0\)… but I guess we’ll see?

The references list the references for this dataset, which is pointing at the hepdata record for now.

Next, the version is the version of the schema set we’re using with pyhf (1.0.0).

And last, but certainly not least… its patches:

patchset.patches

[<pyhf.patchset.Patch object 'C1N2_Wh_hbb_1000_0(1000, 0)' at 0x7f4a547da0a0>,
 <pyhf.patchset.Patch object 'C1N2_Wh_hbb_1000_100(1000, 100)' at 0x7f4a57867b80>,
 <pyhf.patchset.Patch object 'C1N2_Wh_hbb_1000_150(1000, 150)' at 0x7f4a547cea60>,
 <pyhf.patchset.Patch object 'C1N2_Wh_hbb_1000_200(1000, 200)' at 0x7f4a547ce940>,
 <pyhf.patchset.Patch object 'C1N2_Wh_hbb_1000_250(1000, 250)' at 0x7f4a547ce2b0>,
 <pyhf.patchset.Patch object 'C1N2_Wh_hbb_1000_300(1000, 300)' at 0x7f4a29a17af0>,
 <pyhf.patchset.Patch object 'C1N2_Wh_hbb_1000_350(1000, 350)' at 0x7f4a547e2d60>,
 <pyhf.patchset.Patch object 'C1N2_Wh_hbb_1000_400(1000, 400)' at 0x7f4a547e2220>,
 <pyhf.patchset.Patch object 'C1N2_Wh_hbb_1000_50(1000, 50)' at 0x7f4a28f04280>,
 <pyhf.patchset.Patch object 'C1N2_Wh_hbb_150_0(150, 0)' at 0x7f4a28f042b0>,
 <pyhf.patchset.Patch object 'C1N2_Wh_hbb_165_35(165, 35)' at 0x7f4a28f042e0>,
 <pyhf.patchset.Patch object 'C1N2_Wh_hbb_175_0(175, 0)' at 0x7f4a28f04400>,
 <pyhf.patchset.Patch object 'C1N2_Wh_hbb_175_25(175, 25)' at 0x7f4a28f04430>,
 <pyhf.patchset.Patch object 'C1N2_Wh_hbb_190_60(190, 60)' at 0x7f4a28f04490>,
 <pyhf.patchset.Patch object 'C1N2_Wh_hbb_200_0(200, 0)' at 0x7f4a28f044f0>,
 <pyhf.patchset.Patch object 'C1N2_Wh_hbb_200_25(200, 25)' at 0x7f4a28f04550>,
 <pyhf.patchset.Patch object 'C1N2_Wh_hbb_200_50(200, 50)' at 0x7f4a28f045b0>,
 <pyhf.patchset.Patch object 'C1N2_Wh_hbb_225_0(225, 0)' at 0x7f4a28f04610>,
 <pyhf.patchset.Patch object 'C1N2_Wh_hbb_225_25(225, 25)' at 0x7f4a28f04670>,
 <pyhf.patchset.Patch object 'C1N2_Wh_hbb_225_50(225, 50)' at 0x7f4a28f046d0>,
 <pyhf.patchset.Patch object 'C1N2_Wh_hbb_250_0(250, 0)' at 0x7f4a28f04730>,
 <pyhf.patchset.Patch object 'C1N2_Wh_hbb_250_100(250, 100)' at 0x7f4a28f04790>,
 <pyhf.patchset.Patch object 'C1N2_Wh_hbb_250_25(250, 25)' at 0x7f4a28f047f0>,
 <pyhf.patchset.Patch object 'C1N2_Wh_hbb_250_50(250, 50)' at 0x7f4a28f04850>,
 <pyhf.patchset.Patch object 'C1N2_Wh_hbb_250_75(250, 75)' at 0x7f4a28f048b0>,
 <pyhf.patchset.Patch object 'C1N2_Wh_hbb_275_0(275, 0)' at 0x7f4a28f04910>,
 <pyhf.patchset.Patch object 'C1N2_Wh_hbb_275_25(275, 25)' at 0x7f4a28f04970>,
 <pyhf.patchset.Patch object 'C1N2_Wh_hbb_275_50(275, 50)' at 0x7f4a28f049d0>,
 <pyhf.patchset.Patch object 'C1N2_Wh_hbb_275_75(275, 75)' at 0x7f4a28f04a30>,
 <pyhf.patchset.Patch object 'C1N2_Wh_hbb_300_0(300, 0)' at 0x7f4a28f04a90>,
 <pyhf.patchset.Patch object 'C1N2_Wh_hbb_300_150(300, 150)' at 0x7f4a28f04af0>,
 <pyhf.patchset.Patch object 'C1N2_Wh_hbb_300_25(300, 25)' at 0x7f4a28f04b50>,
 <pyhf.patchset.Patch object 'C1N2_Wh_hbb_300_50(300, 50)' at 0x7f4a28f04bb0>,
 <pyhf.patchset.Patch object 'C1N2_Wh_hbb_300_75(300, 75)' at 0x7f4a28f04c10>,
 <pyhf.patchset.Patch object 'C1N2_Wh_hbb_325_0(325, 0)' at 0x7f4a28f04c70>,
 <pyhf.patchset.Patch object 'C1N2_Wh_hbb_325_50(325, 50)' at 0x7f4a28f04cd0>,
 <pyhf.patchset.Patch object 'C1N2_Wh_hbb_350_0(350, 0)' at 0x7f4a28f04d30>,
 <pyhf.patchset.Patch object 'C1N2_Wh_hbb_350_100(350, 100)' at 0x7f4a547d3af0>,
 <pyhf.patchset.Patch object 'C1N2_Wh_hbb_350_150(350, 150)' at 0x7f4a547daee0>,
 <pyhf.patchset.Patch object 'C1N2_Wh_hbb_350_200(350, 200)' at 0x7f4a28f04dc0>,
 <pyhf.patchset.Patch object 'C1N2_Wh_hbb_350_25(350, 25)' at 0x7f4a28f04e20>,
 <pyhf.patchset.Patch object 'C1N2_Wh_hbb_350_50(350, 50)' at 0x7f4a28f04e80>,
 <pyhf.patchset.Patch object 'C1N2_Wh_hbb_350_75(350, 75)' at 0x7f4a28f04ee0>,
 <pyhf.patchset.Patch object 'C1N2_Wh_hbb_375_0(375, 0)' at 0x7f4a28f04f40>,
 <pyhf.patchset.Patch object 'C1N2_Wh_hbb_375_50(375, 50)' at 0x7f4a28f04fa0>,
 <pyhf.patchset.Patch object 'C1N2_Wh_hbb_400_0(400, 0)' at 0x7f4a24c4e040>,
 <pyhf.patchset.Patch object 'C1N2_Wh_hbb_400_100(400, 100)' at 0x7f4a24c4e0a0>,
 <pyhf.patchset.Patch object 'C1N2_Wh_hbb_400_150(400, 150)' at 0x7f4a24c4e100>,
 <pyhf.patchset.Patch object 'C1N2_Wh_hbb_400_200(400, 200)' at 0x7f4a24c4e160>,
 <pyhf.patchset.Patch object 'C1N2_Wh_hbb_400_25(400, 25)' at 0x7f4a24c4e1c0>,
 <pyhf.patchset.Patch object 'C1N2_Wh_hbb_400_250(400, 250)' at 0x7f4a24c4e220>,
 <pyhf.patchset.Patch object 'C1N2_Wh_hbb_400_50(400, 50)' at 0x7f4a24c4e280>,
 <pyhf.patchset.Patch object 'C1N2_Wh_hbb_425_0(425, 0)' at 0x7f4a24c4e2e0>,
 <pyhf.patchset.Patch object 'C1N2_Wh_hbb_450_0(450, 0)' at 0x7f4a24c4e340>,
 <pyhf.patchset.Patch object 'C1N2_Wh_hbb_450_100(450, 100)' at 0x7f4a24c4e3a0>,
 <pyhf.patchset.Patch object 'C1N2_Wh_hbb_450_150(450, 150)' at 0x7f4a24c4e400>,
 <pyhf.patchset.Patch object 'C1N2_Wh_hbb_450_200(450, 200)' at 0x7f4a24c4e460>,
 <pyhf.patchset.Patch object 'C1N2_Wh_hbb_450_250(450, 250)' at 0x7f4a24c4e4c0>,
 <pyhf.patchset.Patch object 'C1N2_Wh_hbb_450_300(450, 300)' at 0x7f4a24c4e520>,
 <pyhf.patchset.Patch object 'C1N2_Wh_hbb_450_50(450, 50)' at 0x7f4a24c4e580>,
 <pyhf.patchset.Patch object 'C1N2_Wh_hbb_500_0(500, 0)' at 0x7f4a24c4e5e0>,
 <pyhf.patchset.Patch object 'C1N2_Wh_hbb_500_100(500, 100)' at 0x7f4a24c4e640>,
 <pyhf.patchset.Patch object 'C1N2_Wh_hbb_500_150(500, 150)' at 0x7f4a24c4e6a0>,
 <pyhf.patchset.Patch object 'C1N2_Wh_hbb_500_200(500, 200)' at 0x7f4a24c4e700>,
 <pyhf.patchset.Patch object 'C1N2_Wh_hbb_500_250(500, 250)' at 0x7f4a24c4e760>,
 <pyhf.patchset.Patch object 'C1N2_Wh_hbb_500_300(500, 300)' at 0x7f4a24c4e7c0>,
 <pyhf.patchset.Patch object 'C1N2_Wh_hbb_500_350(500, 350)' at 0x7f4a24c4e820>,
 <pyhf.patchset.Patch object 'C1N2_Wh_hbb_500_50(500, 50)' at 0x7f4a24c4e880>,
 <pyhf.patchset.Patch object 'C1N2_Wh_hbb_535_400(535, 400)' at 0x7f4a24c4e8e0>,
 <pyhf.patchset.Patch object 'C1N2_Wh_hbb_550_0(550, 0)' at 0x7f4a24c4e940>,
 <pyhf.patchset.Patch object 'C1N2_Wh_hbb_550_100(550, 100)' at 0x7f4a24c4e9a0>,
 <pyhf.patchset.Patch object 'C1N2_Wh_hbb_550_150(550, 150)' at 0x7f4a24c4ea00>,
 <pyhf.patchset.Patch object 'C1N2_Wh_hbb_550_200(550, 200)' at 0x7f4a24c4ea60>,
 <pyhf.patchset.Patch object 'C1N2_Wh_hbb_550_250(550, 250)' at 0x7f4a24c4eac0>,
 <pyhf.patchset.Patch object 'C1N2_Wh_hbb_550_300(550, 300)' at 0x7f4a24c4eb20>,
 <pyhf.patchset.Patch object 'C1N2_Wh_hbb_550_50(550, 50)' at 0x7f4a24c4eb80>,
 <pyhf.patchset.Patch object 'C1N2_Wh_hbb_600_0(600, 0)' at 0x7f4a24c4ebe0>,
 <pyhf.patchset.Patch object 'C1N2_Wh_hbb_600_100(600, 100)' at 0x7f4a24c4ec40>,
 <pyhf.patchset.Patch object 'C1N2_Wh_hbb_600_150(600, 150)' at 0x7f4a24c4eca0>,
 <pyhf.patchset.Patch object 'C1N2_Wh_hbb_600_200(600, 200)' at 0x7f4a24c4ed00>,
 <pyhf.patchset.Patch object 'C1N2_Wh_hbb_600_250(600, 250)' at 0x7f4a24c4ed60>,
 <pyhf.patchset.Patch object 'C1N2_Wh_hbb_600_300(600, 300)' at 0x7f4a24c4edc0>,
 <pyhf.patchset.Patch object 'C1N2_Wh_hbb_600_350(600, 350)' at 0x7f4a24c4ee20>,
 <pyhf.patchset.Patch object 'C1N2_Wh_hbb_600_400(600, 400)' at 0x7f4a24c4ee80>,
 <pyhf.patchset.Patch object 'C1N2_Wh_hbb_600_50(600, 50)' at 0x7f4a24c4eee0>,
 <pyhf.patchset.Patch object 'C1N2_Wh_hbb_650_0(650, 0)' at 0x7f4a24c4ef40>,
 <pyhf.patchset.Patch object 'C1N2_Wh_hbb_650_100(650, 100)' at 0x7f4a24c4efa0>,
 <pyhf.patchset.Patch object 'C1N2_Wh_hbb_650_150(650, 150)' at 0x7f4a24c54040>,
 <pyhf.patchset.Patch object 'C1N2_Wh_hbb_650_200(650, 200)' at 0x7f4a24c540a0>,
 <pyhf.patchset.Patch object 'C1N2_Wh_hbb_650_250(650, 250)' at 0x7f4a24c54100>,
 <pyhf.patchset.Patch object 'C1N2_Wh_hbb_650_300(650, 300)' at 0x7f4a24c54160>,
 <pyhf.patchset.Patch object 'C1N2_Wh_hbb_650_50(650, 50)' at 0x7f4a24c541c0>,
 <pyhf.patchset.Patch object 'C1N2_Wh_hbb_700_0(700, 0)' at 0x7f4a24c54220>,
 <pyhf.patchset.Patch object 'C1N2_Wh_hbb_700_100(700, 100)' at 0x7f4a24c54280>,
 <pyhf.patchset.Patch object 'C1N2_Wh_hbb_700_150(700, 150)' at 0x7f4a24c542e0>,
 <pyhf.patchset.Patch object 'C1N2_Wh_hbb_700_200(700, 200)' at 0x7f4a24c54340>,
 <pyhf.patchset.Patch object 'C1N2_Wh_hbb_700_250(700, 250)' at 0x7f4a24c543a0>,
 <pyhf.patchset.Patch object 'C1N2_Wh_hbb_700_300(700, 300)' at 0x7f4a24c54400>,
 <pyhf.patchset.Patch object 'C1N2_Wh_hbb_700_350(700, 350)' at 0x7f4a24c54460>,
 <pyhf.patchset.Patch object 'C1N2_Wh_hbb_700_400(700, 400)' at 0x7f4a24c544c0>,
 <pyhf.patchset.Patch object 'C1N2_Wh_hbb_700_50(700, 50)' at 0x7f4a24c54520>,
 <pyhf.patchset.Patch object 'C1N2_Wh_hbb_750_100(750, 100)' at 0x7f4a24c54580>,
 <pyhf.patchset.Patch object 'C1N2_Wh_hbb_750_150(750, 150)' at 0x7f4a24c545e0>,
 <pyhf.patchset.Patch object 'C1N2_Wh_hbb_750_200(750, 200)' at 0x7f4a24c54640>,
 <pyhf.patchset.Patch object 'C1N2_Wh_hbb_750_250(750, 250)' at 0x7f4a24c546a0>,
 <pyhf.patchset.Patch object 'C1N2_Wh_hbb_750_300(750, 300)' at 0x7f4a24c54700>,
 <pyhf.patchset.Patch object 'C1N2_Wh_hbb_750_50(750, 50)' at 0x7f4a24c54760>,
 <pyhf.patchset.Patch object 'C1N2_Wh_hbb_800_0(800, 0)' at 0x7f4a24c547c0>,
 <pyhf.patchset.Patch object 'C1N2_Wh_hbb_800_100(800, 100)' at 0x7f4a24c54820>,
 <pyhf.patchset.Patch object 'C1N2_Wh_hbb_800_150(800, 150)' at 0x7f4a24c54880>,
 <pyhf.patchset.Patch object 'C1N2_Wh_hbb_800_200(800, 200)' at 0x7f4a24c548e0>,
 <pyhf.patchset.Patch object 'C1N2_Wh_hbb_800_250(800, 250)' at 0x7f4a24c54940>,
 <pyhf.patchset.Patch object 'C1N2_Wh_hbb_800_300(800, 300)' at 0x7f4a24c549a0>,
 <pyhf.patchset.Patch object 'C1N2_Wh_hbb_800_350(800, 350)' at 0x7f4a24c54a00>,
 <pyhf.patchset.Patch object 'C1N2_Wh_hbb_800_400(800, 400)' at 0x7f4a24c54a60>,
 <pyhf.patchset.Patch object 'C1N2_Wh_hbb_800_50(800, 50)' at 0x7f4a24c54ac0>,
 <pyhf.patchset.Patch object 'C1N2_Wh_hbb_900_0(900, 0)' at 0x7f4a24c54b20>,
 <pyhf.patchset.Patch object 'C1N2_Wh_hbb_900_100(900, 100)' at 0x7f4a24c54b80>,
 <pyhf.patchset.Patch object 'C1N2_Wh_hbb_900_150(900, 150)' at 0x7f4a24c54be0>,
 <pyhf.patchset.Patch object 'C1N2_Wh_hbb_900_200(900, 200)' at 0x7f4a24c54c40>,
 <pyhf.patchset.Patch object 'C1N2_Wh_hbb_900_250(900, 250)' at 0x7f4a24c54ca0>,
 <pyhf.patchset.Patch object 'C1N2_Wh_hbb_900_300(900, 300)' at 0x7f4a24c54d00>,
 <pyhf.patchset.Patch object 'C1N2_Wh_hbb_900_350(900, 350)' at 0x7f4a24c54d60>,
 <pyhf.patchset.Patch object 'C1N2_Wh_hbb_900_400(900, 400)' at 0x7f4a24c54dc0>,
 <pyhf.patchset.Patch object 'C1N2_Wh_hbb_900_50(900, 50)' at 0x7f4a24c54e20>]

So we can see all the patches listed both by name such as C1N2_Wh_hbb_900_250 as well as a pair of points (900, 250). Why is this useful? The PatchSet object acts like a special dictionary look-up where it will grab the patch you need based on the unique key you provide it.

For example, we can look up by name

patchset['C1N2_Wh_hbb_900_250']

<pyhf.patchset.Patch object 'C1N2_Wh_hbb_900_250(900, 250)' at 0x7f4a24c54ca0>

or by the pair of points

patchset[(900, 250)]

<pyhf.patchset.Patch object 'C1N2_Wh_hbb_900_250(900, 250)' at 0x7f4a24c54ca0>

Patches

A pyhf.PatchSet is a collection of pyhf.Patch objects. What is a patch indeed? It contains enough information about how to apply the signal patch to the corresponding background-only workspace (matched by digest).

patch = patchset['C1N2_Wh_hbb_900_250']
print(f'  name: {patch.name}')
print(f'values: {patch.values}')

  name: C1N2_Wh_hbb_900_250
values: (900, 250)

Most importantly, it contains the patch information itself. Specifically, this inherits from the jsonpatch.JsonPatch object, which is a 3rd party module providing native support for json patching in python. That means we can simply apply the patch to our workspace directly!

print(f' samples pre-patch: {pyhf.Workspace(spec).samples}')
print(f'samples post-patch: {pyhf.Workspace(patch.apply(spec)).samples}')

 samples pre-patch: ['diboson', 'multiboson', 'singletop', 'ttbar', 'tth', 'ttv', 'vh', 'wjets', 'zjets']

samples post-patch: ['C1N2_Wh_hbb_900_250', 'diboson', 'multiboson', 'singletop', 'ttbar', 'tth', 'ttv', 'vh', 'wjets', 'zjets']

Or, more quickly, from the PatchSet object:

print(f' samples pre-patch: {pyhf.Workspace(spec).samples}')
print(f'samples post-patch: {pyhf.Workspace(patchset.apply(spec, (900, 250))).samples}')

 samples pre-patch: ['diboson', 'multiboson', 'singletop', 'ttbar', 'tth', 'ttv', 'vh', 'wjets', 'zjets']

samples post-patch: ['C1N2_Wh_hbb_900_250', 'diboson', 'multiboson', 'singletop', 'ttbar', 'tth', 'ttv', 'vh', 'wjets', 'zjets']

Patching via Model Creation

One last way to apply the patching is to, instead of patching workspaces, we patch the models as we build them from the background-only workspace. This maybe makes it easier to treat the background-only workspace as immutable, and patch in signal models when grabbing the model. Check it out.

workspace = pyhf.Workspace(spec)

First, load up our background-only spec into the workspace. Then let’s create a model.

model = workspace.model(patches=[patchset['C1N2_Wh_hbb_900_250']])
print(f'samples (workspace): {workspace.samples}')
print(f'samples (  model  ): {model.config.samples}')

samples (workspace): ['diboson', 'multiboson', 'singletop', 'ttbar', 'tth', 'ttv', 'vh', 'wjets', 'zjets']
samples (  model  ): ['C1N2_Wh_hbb_900_250', 'diboson', 'multiboson', 'singletop', 'ttbar', 'tth', 'ttv', 'vh', 'wjets', 'zjets']

Doing Physics

So we want to try and reproduce part of the contour. At least convince ourselves we’re doing physics and not fauxsics. … Anyway… Let’s remind ourselves of the 1Lbb contour as we don’t have the photographic memory of the ATLAS SUSY conveners

1Lbb exclusion contour

So let’s work around the 700-900 GeV \(\tilde{\chi}_1^\pm, \tilde{\chi}_2^0\) region. We’ll look at two points here:

  • C1N2_Wh_hbb_650_0(650, 0) which is below the contour and excluded

  • C1N2_Wh_hbb_1000_0(1000, 0) which is above the contour and not excluded

Let’s perform a “standard” hypothesis test (with \(\mu = 1\) null BSM hypothesis) on both of these and use the \(\text{CL}_\text{s}\) values to convince ourselves that we just did reproducible physics!?!

Interpolators?

Ok, one last thing I should mention for the HistFitter users. We use nonlinear interpolators by default for the normsys and histosys modifiers. These correspond to code4 and code4p in the HistFactory paper. We’ll additionally configure the modifier settings to switch to these interpolators by default.

modifier_settings={
    "normsys": {"interpcode": "code4"},
    "histosys": {"interpcode": "code4p"},
}

Doing Physics, for real now

model_below = workspace.model(
    patches=[patchset['C1N2_Wh_hbb_650_0']],
    modifier_settings=modifier_settings
)

model_above = workspace.model(
    patches=[patchset['C1N2_Wh_hbb_1000_0']],
    modifier_settings=modifier_settings
)

We’ve made our models. Let’s test hypotheses!

Note: this will not be as instantaneous as our simple models.. but it should still be pretty fast!

result_below = pyhf.infer.hypotest(
    1.0,
    workspace.data(model_below),
    model_below,
    qtilde=True,
    return_expected_set=True
)
print(f"Observed: {result_below[0]}, Expected band: {[exp.tolist() for exp in result_below[1]]}")

/opt/hostedtoolcache/Python/3.8.5/x64/lib/python3.8/site-packages/pyhf/tensor/numpy_backend.py:271: RuntimeWarning: invalid value encountered in log
  return n * np.log(lam) - lam - gammaln(n + 1.0)

Observed: 0.01713750229590585, Expected band: [4.841885092942874e-06, 9.133212662778044e-05, 0.001465995592955356, 0.01732873165693218, 0.12149957599602453]

result_above = pyhf.infer.hypotest(
    1.0,
    workspace.data(model_above),
    model_above,
    qtilde=True,
    return_expected_set=True
)
print(f"Observed: {result_above[0]}, Expected band: {[exp.tolist() for exp in result_above[1]]}")

Observed: 0.5856783708143126, Expected band: [0.04988036057970551, 0.12644646848485364, 0.2925782309330889, 0.5694123817162843, 0.8475953316084923]

And as you can see, we’re getting results that we generally expect. Excluded models are those for which \(\text{CL}_\text{s} < 0.05\). Additionally, you can see that the expected bands \(-2\sigma\) for the \((1000, 0)\) point is just slightly below the observed result for the \((650, 0)\) point which is what we observe in the figure above.

Modifiers Making a Pull Plot