pyhf for upper limit on BF

[GiuliaFrau]

Dear experts,
I hope you can help with this. I am trying to use pyhf to set an upper limit on the branching fraction of a very rare decay. The idea is to use the function pyhf.infer.intervals.upperlimit which takes as input some data and a model. Given that we are performing a blinded analysis, our observed data corresponds to the background, defined in the model together with the signal. Again, as the analysis is blinded I don’t have a value for the signal. What I have is alpha or the so-called single event sensitivity which is related to the number of signal events by the formula: BF = N_sig*alpha. I read the documentation about the modifiers, which seems to be the way to add a normalization factor to the signal events, but it’s not clear to me how to pass these alphas so that they can act as I want. Moreover, the limit is to be evaluated in 16x8 bins for each data-taking year. Do you happen to have any examples about my case? Thanks for your help.
Giulia

[alexander-held]

Hi Giulia, I am not sure I understand the setup correctly, so please correct me where I might be wrong. Let's start with a single bin: you have some expected contribution from background(s) and signal. You know N_sig in that bin, and now you want to calculate a limit for BF = N_sig*alpha. The alpha here would be a normfactor modifier in pyhf, a free-floating multiplicative parameter acting on the signal sample. You could then either:

• set N_sig = 1 in your workspace, effectively re-scaling your alpha by N_sig and the thing called "alpha" in your model then is equal to BF and you can calculate the limit for it,
• use alpha as your POI (N_sig is still the actual number of expected signal events for some benchmark setup), you calculate the limit for alpha and can re-scale that limit manually by multiplying it by N_sig.

I am assuming that there is some sensible value for N_sig that can be obtained, and you could then test your sensitivity with an Asimov dataset (where your expectation is background + that benchmark signal setup).

I am less sure about the specific implications from now including multiple bins. Do you know the signal acceptance across bins, and know N_sig for each bin (up to some constant overall normalization factor)? If yes, then nothing changes about the setup, there are just more bins and channels. Are there multiple potential signal models which may predict different distributions? If there are specific models, they could be tested independently. Are the potential signal contributions per bin not related to each other? Then multiple normalization factors could be used.

[alexander-held]

Sorry, I think I am still not following the setup here. Let's take one step back and talk about this in more general terms. You have some bins, in which you expect contributions from backgrounds and a signal process. You know the number of background events in each bin. Do you have any estimate for the distribution (up to the normalization) of signal events across your bins (e.g. from MC simulation)?

pyhf implements a template fit, so you need to build a template for signal in one way or another. Perhaps it would be helpful if you could describe how you would proceed after unblinding, because it is not clear to me at the moment how the fact that the analysis is blinded matters conceptually when designing the model. I understand that you may not want to fit observed data at the moment, but your model will still need to include signal if you want to evaluate its sensitivity (and you can use it to fit synthetic data instead).

Assuming that your parameter of interest is the branching fraction BF, and that you know alpha in each bin, then you can calculate how much signal you expect in each bin for a given value of BF: N_sig = BF / alpha (assuming I am correctly understanding what is meant by N_sig here, I am not sure about that). Another way to put this is that the signal in each bin scales linearly with BF (so you use that as a normfactor), and the number of signal events across different bins varies like 1/alpha.

Let's say we have two bins and alpha=1 in one bin, alpha=5 in the other bin. The signal template will then have [1.0, 0.2] events, and be scaled by some normfactor called BF. That would then give you N_sig as a function of BF (which is a parameter determined when fitting your model). I am not confident that this is the setup you have in mind, so please let me know how this may differ from what your case requires.

[pkoppenb]

I presume that what Guilia wants (correct me if I am wrong) is the limit given the case of N_sig=0. You see a number of events in your box (that could be 16x8 boxes, not sure how to deal with that) which is consistent with the background expectation, but Poisson-fluctuates. Based on that you set a limit on N_sig, which you can linearly scale to a limit on BF by the factor alpha.

@alexander-held : I presume the expression is more like alpha = 10^-9 and then (measured BF) = (observed N_sig) * alpha, which is the same as what you wrote, but in general we do not know BF a priori. It could even be 0 in the SM.

[alexander-held]

I understand that we may not see signal events in reality, and that the observed data could be completely compatible with background-only. What I am trying to understand is how to build a template for the signal distribution across all bins. That is needed to set a limit on the BF, regardless of whether the best-fit value indicates zero contribution from that template or not (because in the process of determining the limit, we test different values for the signal normalization to find the upper limit).

Am I understanding correctly that alpha can in general vary per bin, but is known for all bins? If it includes some acceptance effects that differ across bins, my proposal would be to include that effect already in the signal template, by re-scaling each signal bin content like I mention in the previous comment.

Even if the potential number of signal events is completely unknown, for the purpose of model-building, we need to create a template with some number of signal events per bin. We can set the amount of signal events in the template to 1 in a given bin for example, and then attach some signal normalization factor mu to it. We can calculate the limit on mu, which translates to a limit on BF when multiplied by alpha. The limit on mu times the yield in that bin should also give the limit on the amount of signal events in that bin. The choice of a yield of 1 is arbitrary here, setting it to 2 will just change the limit on mu by a factor 2 accordingly, resulting in the same number for the limit on signal events.

My proposal in the previous comment was to directly work with BF as the parameter of interest, and build a signal template from the knowledge of alpha in every bin. If I assume that BF = 1.0, can you tell me how many signal events would be expected in a given bin? If yes, then that is the signal template to use. If not, as long as we know the number in all bins up to some overall constant normalization factor, we can calculate everything up to that normalization factor, which would rescale our results. If the latter is true, I would need more details to understand how this kind of analysis can be done at all, unless you would normally quote the result as a function of some currently unknown rescaling factor.

If we do not know at all how potential signal contribution in different bins may relate to each other, we can also scale them by different normalization factors, but I am not sure that this is what you are looking for?

[GiuliaFrau]

Yes, @pkoppenb got the point. The decay we are studying is not allowed in the SM so we don't expect to see any signal. The aim of our study is to improve the current limit on the BF which is of the order of 10e-8. I will try to make my setup more clear.

The analysis is performed in bins of the output of a multivariate classifier (XGB) and a PID variable, 4x4 bins. In each of these 16 bins we have a mass template. For each XGB and PID bin we evaluate the expected number of bkg events n_bkg_i and the fraction of signal events that falls in that particular bin, epsilon_i. The sum of the fractions over the total number of bins is 1.

In order to evaluate the total BF we would do: BF = n_sig * alpha where n_sig corresponds to the number of signal events measured at LHCb while alpha, our overall normalization factor, is evaluated from the MC and brings its own uncertainty. In each XGB and PID we expect to measure a number of signal events n_sig_i which can be written as epsilon_i * n_sig.

So yes, in principle we could assume that alpha can vary per bin and that alpha_i would be given by alpha / epsilon_i.

From what I understood you suggest either to

• assume that our n_sig_i = 1, include alpha_i as normalization factor, obtain a limit on the BF and then scale it by alpha
  or
• assume that BF = 1 so that n_sig_i = 1 / alpha_i, obtain a limit on the BF and then again scale it by alpha

I think this might work. The problem would be that we split our data sample in the three data-taking years, which means that we have three overall normalization factors and we of course want a combined limit for the three years. So my question is, if I use one of the two methods you suggested, for which normalization factor should I rescale my result?

[matthewfeickert]

@GiuliaFrau 👋 I'm trying to go back through some old Discussions and check if things were resolved. It looks like you still had a question. Did that ever get resolved, or is this still something that is unclear?