Computing main model yields

[APMDSLHC]

Dear pyhf developers,

I want to check that all the yields on my main model are positive. I might be wrong but I understood that the model.expected_actualdata method returns the yields of the main model. So I checked with the following example. I compare the yields computed in each bin by simply doing signal+background with the ones obtained using model.expected_actualdata at the point of suggested initial parameters. I checked that the initial parameters give the nominal values of the data, i.e. I checked that all multiplicative modifiers ('normfactor', 'lumi' and 'staterror') parameters are set to 1 and that all additive modifiers ('normsys' and 'histosys') parameters are set to 0. So in principle the two computations should give the same results, but they don't (without taking the yields order into account).
Can someone please tell me what I did wrong?

Just in case, please find below the small code I use to make the computations.
Here I use the ATLAS-SUSY-2019-08 BkgOnly.json file and I patch the signals produced by SModelS
BkgOnly.zip

import pyhf
import json
import jsonpatch

# ATLAS-SUSY-2019-08 json file
with open("./BkgOnly.json", "r") as f:
    bkg = json.load(f)

# Signals produced by SModelS
nsig = [1.4349825452486582, 1.0533138982733217, 1.1158810930352483, 
        0.8276376394629201, 0.8551014982035671, 1.6913336713087228, 
        0.719039660507198, 0.9461983338046028, 1.5376094024329356]

# The patch produced by SModelS
patch = [
  {
    "op": "add",
    "path": "/channels/5/samples/0",
    "value": {
      "data": [
        nsig[0],
        nsig[1],
        nsig[2]
      ],
      "modifiers": [
        {
          "data": None,
          "type": "normfactor",
          "name": "mu_SIG"
        },
        {
          "data": None,
          "type": "lumi",
          "name": "lumi"
        }
      ],
      "name": "bsm"
    }
  },
  {
    "op": "add",
    "path": "/channels/6/samples/0",
    "value": {
      "data": [
        nsig[3],
        nsig[4],
        nsig[5]
      ],
      "modifiers": [
        {
          "data": None,
          "type": "normfactor",
          "name": "mu_SIG"
        },
        {
          "data": None,
          "type": "lumi",
          "name": "lumi"
        }
      ],
      "name": "bsm"
    }
  },
  {
    "op": "add",
    "path": "/channels/7/samples/0",
    "value": {
      "data": [
        nsig[6],
        nsig[7],
        nsig[8]
      ],
      "modifiers": [
        {
          "data": None,
          "type": "normfactor",
          "name": "mu_SIG"
        },
        {
          "data": None,
          "type": "lumi",
          "name": "lumi"
        }
      ],
      "name": "bsm"
    }
  },
  {
    "op": "remove",
    "path": "/channels/4"
  },
  {
    "op": "remove",
    "path": "/channels/3"
  },
  {
    "op": "remove",
    "path": "/channels/2"
  },
  {
    "op": "remove",
    "path": "/channels/1"
  },
  {
    "op": "remove",
    "path": "/channels/0"
  }
]

# Trying to reproduce the issue
llhdSpec = jsonpatch.apply_patch(bkg, patch)
msettings = {'normsys': {'interpcode': 'code4'}, 'histosys': {'interpcode': 'code4p'}}
workspace = pyhf.Workspace(llhdSpec)
model = workspace.model(modifier_settings=msettings)

# Background yields (corresponds to sum of each sample data)
nbg = [4.1, 2.9, 1.1, 4.6, 2.6, 1.4, 8.8, 11.3, 7.3]

print([sig+bg for sig,bg in zip(nsig,nbg)])

init_pars = model.config.suggested_init()
print(model.expected_actualdata(init_pars))

Output:

[5.534982545248658, 3.9533138982733216, 2.2158810930352484, 5.42763763946292, 3.4551014982035673, 3.091333671308723, 9.519039660507199, 12.246198333804603, 8.837609402432935]
[5.6248706751907775, 3.9536840686356154, 2.1237284515078936, 11.682015762293096, 11.26837128568917, 8.468149572953479, 6.204882733857231, 3.66526047735124, 3.078161378599372]

I searched for similar questions and didn't find anything, sorry if it has already been answered.

[alexander-held]

Hi @APMDSLHC, how are you getting this list

nbg = [4.1, 2.9, 1.1, 4.6, 2.6, 1.4, 8.8, 11.3, 7.3]

of background yields? I think we'll need to understand that in order to address why the results do not match.

The background-only yields in the workspace are:

[4.18988813, 2.90037017, 1.00784736, 10.9629761, 10.32217295, 6.93054017, 5.37724509, 2.81015898, 1.38682771]

[APMDSLHC]

Hi @alexander-held,

Thank you for your reply.
I got these from table 5 of the ATLAS-SUSY-2019-08 paper: https://atlas.web.cern.ch/Atlas/GROUPS/PHYSICS/PAPERS/SUSY-2019-08/.

The discrepancy indeed came from the yields in this table not matching what one obtains by summing the samples data of the json file.
If I replace nbg by the background-only yields in the workspace, the difference vanishes.

So I should use the json file data instead of the yields table and the probelm is solved.
Thanks a lot!

[alexander-held]

I cannot comment on what may cause the change between the table and the workspace without a closer look at the paper, but it seems like this is not just due to rounding. Your understanding of expected_actualdata is correct though.

[APMDSLHC]

Ok, thank you very much.

[kratsg]

It looks like this is resolved here. If there is a discrepancy between published HEPData and what the paper quotes, please bring it up with the corresponding experimental collaboration / authors to investigate.

For future reference, there are some details in the pyhf tutorial relating to the upper limit tables (as they also report nbkg there too).