Likelihood Serialization and Preservation¶

Preserved on HEPData¶

As of PyHEP 2020, ATLAS has published 4 full likelihoods to HEPData

import requests
import hashlib
import tarfile
import json

electroweakino_HEPData_URL = "https://www.hepdata.net/record/resource/1267798?view=true"
targz_filename = "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)
assert (
    hashlib.sha256(open(targz_filename, "rb").read()).hexdigest()
    == "64bbbef9f1aaf9e30d75c8975de4789484329b2b825d89331a6f2081661aa728"
)
hashlib.sha256(open(targz_filename, "rb").read()).hexdigest()
# Open as a tarfile
tar_archive = tarfile.open(targz_filename, "r:gz")

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


def get_bkg_and_signal(tarfile, directory_name, model_point):
    background_only = get_json_from_tarfile(tarfile, directory_name + "/BkgOnly.json",)
    patchset = pyhf.PatchSet(
        get_json_from_tarfile(tarfile, directory_name + "/patchset.json")
    )
    signal_patch = patchset[model_point]
    return background_only, signal_patch

import pyhf

def calculate_CLs(bkgonly_json, signal_patch_json):
    """
    Calculate the observed CLs and the expected CLs band from a background only
    and signal patch.
    Args:
        bkgonly_json: The JSON for the background only model
        signal_patch_json: The JSON Patch for the signal model
    Returns:
        CLs_obs: The observed CLs value
        CLs_exp: List of the expected CLs value band
    """
    workspace = pyhf.workspace.Workspace(bkgonly_json)
    model = workspace.model(
        measurement_name=None,
        patches=[signal_patch_json],
        modifier_settings={
            "normsys": {"interpcode": "code4"},
            "histosys": {"interpcode": "code4p"},
        },
    )
    result = pyhf.infer.hypotest(
        1.0, workspace.data(model), model, qtilde=True, return_expected_set=True
    )
    return result[0].tolist()[0], result[-1].ravel().tolist()

oneLbb_background, oneLbb_Wh_hbb_750_100_signal_patch = get_bkg_and_signal(
    tar_archive, "1Lbb-likelihoods-hepdata",(750, 100),  # "C1N2_Wh_hbb_750_100(750, 100)"
    )

CLs_obs, CLs_exp_band = calculate_CLs(oneLbb_background, oneLbb_Wh_hbb_750_100_signal_patch)

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

print('Observed: {}, Expected band: {}'.format(CLs_obs, CLs_exp_band))

Observed: 0.0662858973723076, Expected band: [9.521240025279684e-05, 0.0010232044661932368, 0.009465318735204287, 0.0658248752163788, 0.2824224872612863]

Model Inference Visualization