In [1]:
import cabinetry
We customize the output from cabinetry via a helper function. This is optional, and the logging module can be used directly as well to further customize the behavior.
In [2]:
cabinetry.set_logging()
Download a workspace from HEPData, extract it, pick a signal with pyhf. We use a workspace from an ATLAS search for bottom-squark pair production: JHEP 12 (2019) 060. The corresponding HEPData entry is 10.17182/hepdata.89408.v3.
In [3]:
import json
import pyhf
from pyhf.contrib.utils import download
download("https://www.hepdata.net/record/resource/1935437?view=true", "bottom-squarks")
ws = pyhf.Workspace(json.load(open("bottom-squarks/RegionC/BkgOnly.json")))
patchset = pyhf.PatchSet(json.load(open("bottom-squarks/RegionC/patchset.json")))
ws = patchset.apply(ws, "sbottom_600_280_150")
cabinetry.workspace.save(ws, "bottom-squarks.json")
INFO - pyhf.workspace - Validating spec against schema: workspace.json INFO - pyhf.patchset - Validating spec against schema: patchset.json INFO - pyhf.workspace - Validating spec against schema: workspace.json DEBUG - cabinetry.workspace - saving workspace to bottom-squarks.json
The bottom-squarks.json workspace is now ready to be used. We will run a maximum likelihood fit with cabinetry and visualize the results. First, we have a brief look at the content of the workspace:
In [4]:
!pyhf inspect bottom-squarks.json | head -n 6
Summary
------------------
channels 3
samples 9
parameters 63
modifiers 63
The fit model specified in the workspace is created next.
In [5]:
ws = cabinetry.workspace.load("bottom-squarks.json")
model, data = cabinetry.model_utils.model_and_data(ws)
INFO - pyhf.workspace - Validating spec against schema: workspace.json INFO - pyhf.pdf - Validating spec against schema: model.json INFO - pyhf.pdf - adding modifier EG_RESOLUTION_ALL (1 new nuisance parameters) INFO - pyhf.pdf - adding modifier EG_SCALE_ALL (1 new nuisance parameters) INFO - pyhf.pdf - adding modifier EL_EFF_ChargeIDSel_TOTAL_1NPCOR_PLUS_UNCOR (1 new nuisance parameters) INFO - pyhf.pdf - adding modifier EL_EFF_ID_TOTAL_1NPCOR_PLUS_UNCOR (1 new nuisance parameters) INFO - pyhf.pdf - adding modifier EL_EFF_Iso_TOTAL_1NPCOR_PLUS_UNCOR (1 new nuisance parameters) INFO - pyhf.pdf - adding modifier EL_EFF_Reco_TOTAL_1NPCOR_PLUS_UNCOR (1 new nuisance parameters) INFO - pyhf.pdf - adding modifier EL_EFF_TriggerEff_TOTAL_1NPCOR_PLUS_UNCOR (1 new nuisance parameters) INFO - pyhf.pdf - adding modifier EL_EFF_Trigger_TOTAL_1NPCOR_PLUS_UNCOR (1 new nuisance parameters) INFO - pyhf.pdf - adding modifier FT_EFF_B_systematics (1 new nuisance parameters) INFO - pyhf.pdf - adding modifier FT_EFF_C_systematics (1 new nuisance parameters) INFO - pyhf.pdf - adding modifier FT_EFF_Light_systematics (1 new nuisance parameters) INFO - pyhf.pdf - adding modifier FT_EFF_extrapolation (1 new nuisance parameters) INFO - pyhf.pdf - adding modifier FT_EFF_extrapolation_from_charm (1 new nuisance parameters) INFO - pyhf.pdf - adding modifier JET_EtaIntercalibration_NonClosure_highE (1 new nuisance parameters) INFO - pyhf.pdf - adding modifier JET_EtaIntercalibration_NonClosure_negEta (1 new nuisance parameters) INFO - pyhf.pdf - adding modifier JET_EtaIntercalibration_NonClosure_posEta (1 new nuisance parameters) INFO - pyhf.pdf - adding modifier JET_Flavor_Response (1 new nuisance parameters) INFO - pyhf.pdf - adding modifier JET_GroupedNP_1 (1 new nuisance parameters) INFO - pyhf.pdf - adding modifier JET_GroupedNP_2 (1 new nuisance parameters) INFO - pyhf.pdf - adding modifier JET_GroupedNP_3 (1 new nuisance parameters) INFO - pyhf.pdf - adding modifier JET_JER_DataVsMC (1 new nuisance parameters) INFO - pyhf.pdf - adding modifier JET_JER_EffectiveNP_1 (1 new nuisance parameters) INFO - pyhf.pdf - adding modifier JET_JER_EffectiveNP_2 (1 new nuisance parameters) INFO - pyhf.pdf - adding modifier JET_JER_EffectiveNP_3 (1 new nuisance parameters) INFO - pyhf.pdf - adding modifier JET_JER_EffectiveNP_4 (1 new nuisance parameters) INFO - pyhf.pdf - adding modifier JET_JER_EffectiveNP_5 (1 new nuisance parameters) INFO - pyhf.pdf - adding modifier JET_JER_EffectiveNP_6 (1 new nuisance parameters) INFO - pyhf.pdf - adding modifier JET_JER_EffectiveNP_7restTerm (1 new nuisance parameters) INFO - pyhf.pdf - adding modifier JET_JvtEfficiency (1 new nuisance parameters) INFO - pyhf.pdf - adding modifier MET_SoftTrk_ResoPara (1 new nuisance parameters) INFO - pyhf.pdf - adding modifier MET_SoftTrk_ResoPerp (1 new nuisance parameters) INFO - pyhf.pdf - adding modifier MET_SoftTrk_Scale (1 new nuisance parameters) INFO - pyhf.pdf - adding modifier MUON_EFF_BADMUON_STAT (1 new nuisance parameters) INFO - pyhf.pdf - adding modifier MUON_EFF_BADMUON_SYS (1 new nuisance parameters) INFO - pyhf.pdf - adding modifier MUON_EFF_ISO_STAT (1 new nuisance parameters) INFO - pyhf.pdf - adding modifier MUON_EFF_ISO_SYS (1 new nuisance parameters) INFO - pyhf.pdf - adding modifier MUON_EFF_RECO_STAT (1 new nuisance parameters) INFO - pyhf.pdf - adding modifier MUON_EFF_RECO_SYS (1 new nuisance parameters) INFO - pyhf.pdf - adding modifier MUON_EFF_TTVA_STAT (1 new nuisance parameters) INFO - pyhf.pdf - adding modifier MUON_EFF_TTVA_SYS (1 new nuisance parameters) INFO - pyhf.pdf - adding modifier MUON_EFF_TrigStatUncertainty (1 new nuisance parameters) INFO - pyhf.pdf - adding modifier MUON_EFF_TrigSystUncertainty (1 new nuisance parameters) INFO - pyhf.pdf - adding modifier MUON_ID (1 new nuisance parameters) INFO - pyhf.pdf - adding modifier MUON_MS (1 new nuisance parameters) INFO - pyhf.pdf - adding modifier MUON_SAGITTA_RESBIAS (1 new nuisance parameters) INFO - pyhf.pdf - adding modifier MUON_SAGITTA_RHO (1 new nuisance parameters) INFO - pyhf.pdf - adding modifier MUON_SCALE (1 new nuisance parameters) INFO - pyhf.pdf - adding modifier ttbar_FSR (1 new nuisance parameters) INFO - pyhf.pdf - adding modifier ttbar_Gen (1 new nuisance parameters) INFO - pyhf.pdf - adding modifier ttbar_ISR_Down (1 new nuisance parameters) INFO - pyhf.pdf - adding modifier ttbar_ISR_Up (1 new nuisance parameters) INFO - pyhf.pdf - adding modifier ttbar_PS (1 new nuisance parameters) INFO - pyhf.pdf - adding modifier Z_theory_SR (1 new nuisance parameters) INFO - pyhf.pdf - adding modifier lumi (1 new nuisance parameters) INFO - pyhf.pdf - adding modifier mu_z (1 new nuisance parameters) INFO - pyhf.pdf - adding modifier mu_SIG (1 new nuisance parameters) INFO - pyhf.pdf - adding modifier mu_ttbar (1 new nuisance parameters) INFO - pyhf.pdf - adding modifier SigRad (1 new nuisance parameters) INFO - pyhf.pdf - adding modifier ttH_theory (1 new nuisance parameters) INFO - pyhf.pdf - adding modifier ttZ_theory (1 new nuisance parameters) INFO - pyhf.pdf - adding modifier staterror_CRtt_cuts (1 new nuisance parameters) INFO - pyhf.pdf - adding modifier staterror_CRz_cuts (1 new nuisance parameters) INFO - pyhf.pdf - adding modifier staterror_SR_metsigST (4 new nuisance parameters)
We can take a look at a yield table for this model. We first generate the pre-fit model prediction, and then pass it to a function to produce a yield table from it.
In [6]:
model_prefit = cabinetry.model_utils.prediction(model)
_ = cabinetry.tabulate.yields(model_prefit, data)
DEBUG - cabinetry.model_utils - total stdev is [[7.07], [5.02], [4.04, 2.05, 0.964, 0.581]] DEBUG - cabinetry.model_utils - total stdev per channel is [7.07, 5.02, 5.42] INFO - cabinetry.tabulate - yields per bin for pre-fit model prediction: ╒═════════════════════╤═══════════════╤══════════════╤═══════════════╤══════════════╤═════════════╤═════════════╕ │ sample │ CRtt_cuts │ CRz_cuts │ SR_metsigST │ │ │ │ │ │ bin 1 │ bin 1 │ bin 1 │ bin 2 │ bin 3 │ bin 4 │ ╞═════════════════════╪═══════════════╪══════════════╪═══════════════╪══════════════╪═════════════╪═════════════╡ │ W │ 17.26 ± 2.90 │ 0.00 ± 0.00 │ 2.16 ± 0.47 │ 0.63 ± 0.16 │ 0.24 ± 0.08 │ 0.42 ± 0.15 │ ├─────────────────────┼───────────────┼──────────────┼───────────────┼──────────────┼─────────────┼─────────────┤ │ Z │ 0.32 ± 0.09 │ 53.18 ± 1.85 │ 8.77 ± 1.69 │ 5.91 ± 1.36 │ 1.98 ± 0.60 │ 1.12 ± 0.29 │ ├─────────────────────┼───────────────┼──────────────┼───────────────┼──────────────┼─────────────┼─────────────┤ │ diboson │ 2.77 ± 0.36 │ 5.92 ± 1.76 │ 0.94 ± 0.16 │ 0.59 ± 0.11 │ 0.17 ± 0.03 │ 0.13 ± 0.11 │ ├─────────────────────┼───────────────┼──────────────┼───────────────┼──────────────┼─────────────┼─────────────┤ │ sbottom_600_280_150 │ 14.76 ± 1.40 │ 0.34 ± 0.03 │ 0.98 ± 0.08 │ 0.00 ± 0.00 │ 0.13 ± 0.02 │ 0.00 ± 0.00 │ ├─────────────────────┼───────────────┼──────────────┼───────────────┼──────────────┼─────────────┼─────────────┤ │ st │ 30.75 ± 0.78 │ 0.00 ± 0.00 │ 3.46 ± 1.19 │ 1.56 ± 0.39 │ 0.87 ± 0.27 │ 0.56 ± 0.17 │ ├─────────────────────┼───────────────┼──────────────┼───────────────┼──────────────┼─────────────┼─────────────┤ │ ttH │ 2.88 ± 0.90 │ 0.10 ± 0.04 │ 0.16 ± 0.06 │ 0.04 ± 0.04 │ 0.08 ± 0.03 │ 0.00 ± 0.00 │ ├─────────────────────┼───────────────┼──────────────┼───────────────┼──────────────┼─────────────┼─────────────┤ │ ttW │ 3.58 ± 0.47 │ 0.00 ± 0.00 │ 0.29 ± 0.05 │ 0.16 ± 0.04 │ 0.14 ± 0.03 │ 0.04 ± 0.01 │ ├─────────────────────┼───────────────┼──────────────┼───────────────┼──────────────┼─────────────┼─────────────┤ │ ttZ │ 7.61 ± 1.16 │ 13.08 ± 2.26 │ 2.24 ± 0.36 │ 0.86 ± 0.17 │ 0.38 ± 0.10 │ 0.21 ± 0.05 │ ├─────────────────────┼───────────────┼──────────────┼───────────────┼──────────────┼─────────────┼─────────────┤ │ ttbar │ 108.26 ± 2.75 │ 3.09 ± 2.77 │ 4.96 ± 2.98 │ 1.38 ± 0.76 │ 0.43 ± 0.33 │ 0.16 ± 0.15 │ ├─────────────────────┼───────────────┼──────────────┼───────────────┼──────────────┼─────────────┼─────────────┤ │ total │ 188.19 ± 7.07 │ 75.71 ± 5.02 │ 23.95 ± 4.04 │ 11.11 ± 2.05 │ 4.42 ± 0.96 │ 2.64 ± 0.58 │ ├─────────────────────┼───────────────┼──────────────┼───────────────┼──────────────┼─────────────┼─────────────┤ │ data │ 143.00 │ 73.00 │ 28.00 │ 12.00 │ 4.00 │ 3.00 │ ╘═════════════════════╧═══════════════╧══════════════╧═══════════════╧══════════════╧═════════════╧═════════════╛ INFO - cabinetry.tabulate - saving table as tables/yields_per_bin_pre-fit.txt
We can also visualize the pre-fit model prediction and compare it to data. the visualize.data_mc function returns a list of dictionaries containing the matplotlib figures, which we could use to customize them as needed. We do not need to customize anything here, so we assign the return value to _.
In [7]:
_ = cabinetry.visualize.data_mc(model_prefit, data)
INFO - cabinetry.visualize.utils - saving figure as figures/CRtt_cuts_prefit.pdf INFO - cabinetry.visualize.utils - saving figure as figures/CRz_cuts_prefit.pdf INFO - cabinetry.visualize.utils - saving figure as figures/SR_metsigST_prefit.pdf
Next up is a maximum likelihood fit:
In [8]:
fit_results = cabinetry.fit.fit(model, data)
INFO - cabinetry.fit - performing maximum likelihood fit INFO - cabinetry.fit - Migrad status: ┌─────────────────────────────────────────────────────────────────────────┐ │ Migrad │ ├──────────────────────────────────┬──────────────────────────────────────┤ │ FCN = 106.2 │ Nfcn = 7421 │ │ EDM = 9.07e-05 (Goal: 0.0002) │ │ ├──────────────────────────────────┼──────────────────────────────────────┤ │ Valid Minimum │ SOME Parameters at limit │ ├──────────────────────────────────┼──────────────────────────────────────┤ │ Below EDM threshold (goal x 10) │ Below call limit │ ├───────────────┬──────────────────┼───────────┬─────────────┬────────────┤ │ Covariance │ Hesse ok │ Accurate │ Pos. def. │ Not forced │ └───────────────┴──────────────────┴───────────┴─────────────┴────────────┘ DEBUG - cabinetry.fit - -2 log(L) = 106.162804 at best-fit point INFO - cabinetry.fit - fit results (with symmetric uncertainties): INFO - cabinetry.fit - EG_RESOLUTION_ALL = -0.0009 +/- 0.9970 INFO - cabinetry.fit - EG_SCALE_ALL = 0.0025 +/- 0.9954 INFO - cabinetry.fit - EL_EFF_ChargeIDSel_TOTAL_1NPCOR_PLUS_UNCOR = -0.0013 +/- 0.9934 INFO - cabinetry.fit - EL_EFF_ID_TOTAL_1NPCOR_PLUS_UNCOR = -0.0553 +/- 0.9911 INFO - cabinetry.fit - EL_EFF_Iso_TOTAL_1NPCOR_PLUS_UNCOR = -0.0744 +/- 0.9900 INFO - cabinetry.fit - EL_EFF_Reco_TOTAL_1NPCOR_PLUS_UNCOR = -0.0063 +/- 0.9933 INFO - cabinetry.fit - EL_EFF_TriggerEff_TOTAL_1NPCOR_PLUS_UNCOR = -0.0013 +/- 0.9934 INFO - cabinetry.fit - EL_EFF_Trigger_TOTAL_1NPCOR_PLUS_UNCOR = -0.0050 +/- 0.9933 INFO - cabinetry.fit - FT_EFF_B_systematics = -0.0044 +/- 0.9932 INFO - cabinetry.fit - FT_EFF_C_systematics = 0.0092 +/- 0.9883 INFO - cabinetry.fit - FT_EFF_Light_systematics = 0.0003 +/- 0.9860 INFO - cabinetry.fit - FT_EFF_extrapolation = 0.0142 +/- 0.9921 INFO - cabinetry.fit - FT_EFF_extrapolation_from_charm = 0.0208 +/- 1.0045 INFO - cabinetry.fit - JET_EtaIntercalibration_NonClosure_highE = -0.0020 +/- 0.9934 INFO - cabinetry.fit - JET_EtaIntercalibration_NonClosure_negEta = -0.0048 +/- 0.9914 INFO - cabinetry.fit - JET_EtaIntercalibration_NonClosure_posEta = -0.0012 +/- 0.9939 INFO - cabinetry.fit - JET_Flavor_Response = 0.1369 +/- 1.0662 INFO - cabinetry.fit - JET_GroupedNP_1 = -0.2276 +/- 1.0015 INFO - cabinetry.fit - JET_GroupedNP_2 = -0.1853 +/- 1.0099 INFO - cabinetry.fit - JET_GroupedNP_3 = -0.0194 +/- 0.9981 INFO - cabinetry.fit - JET_JER_DataVsMC = -0.0552 +/- 0.9699 INFO - cabinetry.fit - JET_JER_EffectiveNP_1 = -0.1285 +/- 1.5951 INFO - cabinetry.fit - JET_JER_EffectiveNP_2 = 0.0292 +/- 2.1415 INFO - cabinetry.fit - JET_JER_EffectiveNP_3 = 0.1663 +/- 1.6090 INFO - cabinetry.fit - JET_JER_EffectiveNP_4 = 0.0144 +/- 1.4657 INFO - cabinetry.fit - JET_JER_EffectiveNP_5 = 0.0352 +/- 1.1874 INFO - cabinetry.fit - JET_JER_EffectiveNP_6 = -0.0142 +/- 1.0992 INFO - cabinetry.fit - JET_JER_EffectiveNP_7restTerm = 0.0461 +/- 1.5882 INFO - cabinetry.fit - JET_JvtEfficiency = -0.0052 +/- 0.8251 INFO - cabinetry.fit - MET_SoftTrk_ResoPara = 0.0075 +/- 1.0131 INFO - cabinetry.fit - MET_SoftTrk_ResoPerp = -0.0173 +/- 0.9183 INFO - cabinetry.fit - MET_SoftTrk_Scale = -0.0089 +/- 1.0006 INFO - cabinetry.fit - MUON_EFF_BADMUON_STAT = -0.0013 +/- 0.9934 INFO - cabinetry.fit - MUON_EFF_BADMUON_SYS = -0.0013 +/- 0.9934 INFO - cabinetry.fit - MUON_EFF_ISO_STAT = -0.0026 +/- 0.9934 INFO - cabinetry.fit - MUON_EFF_ISO_SYS = -0.0077 +/- 0.9933 INFO - cabinetry.fit - MUON_EFF_RECO_STAT = -0.0034 +/- 0.9933 INFO - cabinetry.fit - MUON_EFF_RECO_SYS = -0.0124 +/- 0.9933 INFO - cabinetry.fit - MUON_EFF_TTVA_STAT = -0.0022 +/- 0.9933 INFO - cabinetry.fit - MUON_EFF_TTVA_SYS = -0.0017 +/- 0.9934 INFO - cabinetry.fit - MUON_EFF_TrigStatUncertainty = -0.0032 +/- 0.9934 INFO - cabinetry.fit - MUON_EFF_TrigSystUncertainty = -0.0036 +/- 0.9933 INFO - cabinetry.fit - MUON_ID = 0.0078 +/- 0.9826 INFO - cabinetry.fit - MUON_MS = 0.0046 +/- 1.0095 INFO - cabinetry.fit - MUON_SAGITTA_RESBIAS = 0.0011 +/- 0.9376 INFO - cabinetry.fit - MUON_SAGITTA_RHO = -0.0013 +/- 0.9934 INFO - cabinetry.fit - MUON_SCALE = -0.0026 +/- 0.9916 INFO - cabinetry.fit - ttbar_FSR = 0.0836 +/- 1.0806 INFO - cabinetry.fit - ttbar_Gen = 0.4473 +/- 0.9446 INFO - cabinetry.fit - ttbar_ISR_Down = -0.0090 +/- 0.9766 INFO - cabinetry.fit - ttbar_ISR_Up = 0.3189 +/- 1.1357 INFO - cabinetry.fit - ttbar_PS = 0.2176 +/- 1.1645 INFO - cabinetry.fit - Z_theory_SR = 0.1810 +/- 0.9604 INFO - cabinetry.fit - lumi = 1.0000 +/- 0.0169 INFO - cabinetry.fit - mu_z = 1.0437 +/- 0.1717 INFO - cabinetry.fit - mu_SIG = 0.0001 +/- 8.0882 INFO - cabinetry.fit - mu_ttbar = 0.7995 +/- 0.0960 INFO - cabinetry.fit - SigRad = 0.0001 +/- 0.9933 INFO - cabinetry.fit - ttH_theory = -0.0041 +/- 0.9938 INFO - cabinetry.fit - ttZ_theory = -0.0148 +/- 0.9938 INFO - cabinetry.fit - staterror_CRtt_cuts[0] = 0.9990 +/- 0.0254 INFO - cabinetry.fit - staterror_CRz_cuts[0] = 0.9975 +/- 0.0342 INFO - cabinetry.fit - staterror_SR_metsigST[0] = 1.0125 +/- 0.0659 INFO - cabinetry.fit - staterror_SR_metsigST[1] = 1.0048 +/- 0.1097 INFO - cabinetry.fit - staterror_SR_metsigST[2] = 1.0023 +/- 0.1421 INFO - cabinetry.fit - staterror_SR_metsigST[3] = 1.0099 +/- 0.1611
We can now visualize the post-fit distributions. To do so, we need a post-fit model prediction. It is obtained like the pre-fit model prediction, but this time with an additional argument to pass in the fit results.
In [9]:
model_postfit = cabinetry.model_utils.prediction(model, fit_results=fit_results)
_ = cabinetry.visualize.data_mc(model_postfit, data)
DEBUG - cabinetry.model_utils - total stdev is [[120], [8.67], [8.46, 1.84, 1.36, 0.597]] DEBUG - cabinetry.model_utils - total stdev per channel is [120, 8.67, 10.1] INFO - cabinetry.visualize.utils - saving figure as figures/CRtt_cuts_postfit.pdf INFO - cabinetry.visualize.utils - saving figure as figures/CRz_cuts_postfit.pdf INFO - cabinetry.visualize.utils - saving figure as figures/SR_metsigST_postfit.pdf
The nuisance parameter pulls and correlations are visualized below.
In [10]:
cabinetry.visualize.pulls(fit_results, exclude="mu_SIG")
INFO - cabinetry.visualize.utils - saving figure as figures/pulls.pdf
Out[10]:
In [11]:
cabinetry.visualize.correlation_matrix(fit_results, pruning_threshold=0.2)
INFO - cabinetry.visualize.utils - saving figure as figures/correlation_matrix.pdf
Out[11]:
