staterror calculation

[zjkjsd]

Aloha @alexander-held,

Thank you for presenting at the Belle II pyhf workshop earlier; it was very informative! Would it be possible for you to answer the following questions?

1. Is it possible that cabinetry may support additional file formats for ntuple inputs in the future, such as parquet?
2. In the tutorials, I noticed that only signal samples have a normfactor modifier in the workspace.json file, while background samples do not. Is this because the normfactor for backgrounds is always expected to be 1? I am a bit confused by this, as I believe the amount of background in the model may depend on the amount of Monte Carlo used. Or is this just for demonstration purposes?
3. In the config_example.yml file, I noticed that some samples have disableStaterror turned on. Is this for samples that are not from MC? However, I thought even samples not from MC should still have statistical uncertainties. I may have misunderstood something here.
4. If disableStaterror == False, I noticed that a staterror modifier is automatically added to the sample. This is very helpful information as I was confused about what modifier to use for statistical uncertainty when studying pyhf. However, how are the numbers in staterror calculated? Is it $\sigma_{cb} = \sqrt{\sum_s\delta_{csb}}/\sum_s \nu^0_{csb}$? I tried to calculate the numbers by hand using the example from the tutorial, but couldn't get the same results.
   Attached, please find the workspace.json file from the tutorial.
   example_workspace.txt

# Signal      "data": [
                        0.0,
                        1.585369130539275,
                        23.61642679992917,
                        24.548922232659436
                    ],
                    "modifiers": [
                        {
                            "data": [
                                0.0,
                                0.19931165969755615,
                                0.7741071346184436,
                                0.7883050276148931
                            ],
                            "name": "staterror_Signal_region",
                            "type": "staterror"

# Bkg          "data": [
                        112.73896936192415,
                        128.62169538923106,
                        88.107008384283,
                        55.246070720721384
                    ],
                    "modifiers": [
                        {
                            "data": [
                                4.761366784144052,
                                5.106450362569687,
                                4.21104366965744,
                                3.34933334869164
                            ],

If I want to calculate the last bin in staterror_Signal_region, would I do the following:
$$\frac{\sqrt{\sqrt{24.548922232659436}+\sqrt{55.246070720721384}}}{24.548922232659436+55.246070720721384}\times24.548922232659436=1.082799$$
which is different from the cabinetry calculation result of 0.7883050276148931. I must have misunderstood something.

5. How does an incorrect staterror affect the fitting? Would the algorithm mistakenly search for the minimum in a larger space if the staterror is larger than the correct values?

Thank you again for your time and assistance!
Mahalo,
Boyang

[alexander-held]

Hi, sure, please find some answers below.

1. Certainly possible, no current concrete plans for it. It certainly would not be a very high priority item for me right now, but happy to talk about this if somebody is willing to provide an implementation. I have not looked into how difficult that would be to implement (which will mostly come down to how many things do not trivially map from ROOT files to Parquet). You can currently work around this either via https://cabinetry.readthedocs.io/en/latest/advanced.html#overrides-for-template-building or by building the histograms yourself externally and using cabinetry.templates.collect.
2. This is just an example for demonstration. The choice here only depends on the analysis you are interested in. It is also not uncommon to have multiple normalization factors that affect different samples.
3. This is also just for demonstration. One somewhat common use case is that your signal sample may contain a lot of simulated events (and therefore have a very small MC statistical uncertainty) while your backgrounds have larger uncertainties. In such a case you would probably not want to scale your signal if the (large, due to the backgrounds) staterror parameter gets pulled and therefore may choose to exclude the signal.
4. The statistical uncertainty calculation here is more complicated, as it depends on the distribution of weights in your samples. If all those weights are 1 then you arrive at the sqrt(N) calculation, but that is typically not the case. The uncertainty on the yield in a bin, filled with events of weight w_i, is sqrt(sum_i(w_i^2)) with the yield being sum_i w_i. In the case of cabinetry, that calculation is handled by boost-histogram automatically. The uncertainty sums in quadrature when you have contributions from multiple samples.
5. The fit result only makes physical sense if your model does. I cannot answer that in general, as it will be highly case-dependent how much of a difference this makes. Make sure your staterror terms make sense to avoid this problem.

[zjkjsd]

Thank you for your helpful answer! I have some follow-up questions, thank you in advance!

2. Yeah, I remember your answer here demonstrated multiple normfactor acting on samples. Is it correct to assume that there should be at least one normfactor for a sample in practice? If not, in what case, would a sample not need a normfactor?

3. The example is very helpful!

4. I am able to reproduce the correct staterror now! You answered this question in the workshop but I didn't understand until now, thank you very much! This question has puzzled me for months!

   1. Does $\delta_{csb}$ miss a ^2 here in the numerator?
   2. Could you give me some hints about how auxiliary data is calculated for staterror? I noticed that the auxi_data are 1 when parameters are 1, but how can we get that?
   3. I have never thought to assign weights to events, why are the weights typically not 1 and different between events? I guess this probably depends on the analysis, would it be possible for you to give an example or point me to relevant documents?

Mahalo,
Boyang

[alexander-held]

2. No, you do not always need a normfactor per sample. If you have a prediction for the normalization (and potentially an uncertainty for that prediction) that you can trust, you can just use that directly. This is a physics question that depends a lot on your use case and I cannot answer this in general.

4.i. Yes, that looks like a mistake, I'll check with the pyhf team.

4.ii. The auxiliary data is 1 by convention for staterror. What you're looking for is probably the width of the constraint term, that is what relates to the actual statistical uncertainties. You can get that via model.config.param_set("parameter name").width().

4.iii. Have a look at Monte Carlo event generators (MadGraph for example), it can be a practical consequence of how the event generation is done (including weights being negative). In addition to that, in the ATLAS/CMS applications I am most used to, it is also very common to pick up additional event-level weights that relate to calibration and systematic uncertainties.

[zjkjsd]

Thank you very much! The examples you gave are very inspiring and helpful! I have a few follow-up questions for Q.4:

1. I don't fully understand the terminology in the table below (quoted from the pyhf doc). For example, the staterror modifier for a given sample in a workspace will have data, name, and type 3 sections. What is the data of a staterror modifier? Is it the Modification $\gamma_b$? or the Input $\delta_b^2 = \sum_s\delta^2_{sb}$? or the squared root of Input, i.e. $\delta_b$?
2. Is it correct that the $\gamma_b$ and $\delta_b$ are the mean and std of the constraint Gaussian term? If the data of a staterror is $\delta_b$, how is $\gamma_b$ calculated in cabinetry or pyhf?
3. I noticed that for many modifiers, the mean of the constraint term is the Modification. Is this actually how we define the Modification?

Description	Modification	Constraint Term	Input
MC Stat. Uncertainty	$\kappa_{scb}(\gamma_b) = \gamma_b$	$\prod_b \mathrm{Gaus}\left(a_{\gamma_b} = 1\middle|,\gamma_b,\delta_b\right)$	$\delta_b^2 = \sum_s\delta^2_{sb}$

4. Thanks for directing me to MadGraph, @nlo generator. The interference could cause the weights to be negative if I understand correctly. I am wondering if it makes sense that all events have a generator weight of 1, for example, for B meson decays with EvtGen where events are produced according to correct SM branching fractions and decay models. I am still learning, but I feel the only weight people talk about in Belle II is the additional event-level weights that relate to calibration and systematic uncertainties, and all MC events have a default generator weight of 1 and this weight is currently not changeable.
5. I also want to ask how to draw toy MC in cabinetry or pyhf. I have seen code like the below with pyhf:

toys = model.make_pdf(pyhf.tensorlib.astensor(np.asarray(model.config.suggested_init()))).sample((n_toys,))

If there is only 1 channel, my understanding of the code is, to sum all samples first and draw samples from the Poisson distribution defined at each bin with the mean of the expected bin count. Is this correct?
If I get an independent MC sample and want to draw toys from it, can I do the following?

toys = []
def compute_toys(itry):
    rng = np.random.default_rng(itry)
    toy = rng.poisson(MC) # MC is the histogram of a MC sample
    toys.append(toy)

6. Other people have suggested I use inverse transform sampling to make toys, i.e. calculate a cdf from the histogram and use a uniform distribution together to sample events one by one (i.e. a pseudo experiment). Would this be different from the toy methods above? I would greatly appreciate it if you could share some insights! (some of my questions might be a little general, but I couldn't think of anywhere to ask that could give me a better answer than here. Thank you!)