On new physics searches with multidimensional differential shapes

TL;DR

This paper explores how multidimensional differential rates can enhance model discrimination in LHC new physics searches, despite limited impact on discovery potential, and discusses related systematic uncertainties.

Contribution

It provides a detailed analysis of the role of multidimensional shape information in LHC analyses and highlights subtle effects affecting their interpretability.

Findings

Multidimensional rates improve model discrimination.

Limited impact on discovery sensitivity.

Systematic uncertainties and interference effects are critical considerations.

Abstract

In the context of upcoming new physics searches at the LHC, we investigate the impact of multidimensional differential rates in typical LHC analyses. We discuss the properties of shape information, and argue that multidimensional rates bring limited information in the scope of a discovery, but can have a large impact on model discrimination. We also point out subtleties about systematic uncertainties cancellations and the Cauchy-Schwarz bound on interference terms.

Figures

Equations

$$L(\theta )=L_{\mathrm{tot}}(\theta )L_{\mathrm{shape}}(\theta ),$$

$$I_{\mathrm{shape}}^{2D}\le I_{\mathrm{shape}}^{1D,a}+I_{\mathrm{shape}}^{1D,b},$$

$$L(\theta )\approx L_{\max }\exp \left(-I\frac{(\theta -\bar{\theta }{)}^2}{2}\right),$$

$$\log B_0=\alpha n_{\mathrm{obs}}\frac{({\theta }_0-\bar{\theta }{)}^2}{2}-\frac{1}{2}\log \left(\frac{V\alpha n_{\mathrm{obs}}}{2\pi }\right),$$

$$\log B_0^{2D}\le \log B_0^{1D,a}+\log B_0^{1D,b}.$$

$$\log B_{12}=\log \left(\frac{L_{\max ,1}}{L_{\max ,2}}\right)-\log \left(\frac{{\alpha }_1}{{\alpha }_2}\right).$$

$$L_{\max ,1}/L_{\max ,2}\sim \exp (\beta n_{\mathrm{obs}}),$$

$$\mathcal{L}={\mathcal{L}}_{\mathrm{SM}}+c_{HW}{O}_{HW}+{\tilde{c}}_{HW}{\tilde{O}}_{HW},$$

$$O_{HW}(\tilde{O}_{HW})=\frac{2ig}{m_{W}^{2}}\big{[}D^{\mu}\Phi^{\dagger}T_{2k}D^{\nu}\Phi\big{]}W_{\mu\nu}^{k}(\tilde{W}_{\mu\nu}^{k}).$$

$$n_r=n{\int }_{{\mathcal{D}}_r}{f}_X(x)dx.$$

$$L(\theta )\equiv \mathrm{Pr}(\mathrm{data}|H,\theta ).$$

$$L=\prod _r{L}_r.$$

$$L_r(\theta )=n_r(\theta {)}^{{\hat{n}}_r}{e}^{-n_r(\theta )}.$$

$$L_{\mathrm{tot}}(\theta )=n(\theta {)}^{\hat{n}}{e}^{-n(\theta )}$$

$$L_{\mathrm{shape}}(\theta )=\prod _r{\left(\frac{n_r(\theta )}{n(\theta )}\right)}^{\hat{n}}.$$

$$\frac{{\int }_{{\Omega }^{1-\alpha }}d\theta p(\theta )}{{\int }_{\Omega }d\theta p(\theta )}=1-\alpha ,$$

$$B_{01}=\frac{{\int }_{{\Omega }_1}L({\theta }_1){\pi }_1({\theta }_1)}{{\int }_{{\Omega }_0}L({\theta }_0){\pi }_0({\theta }_0)},$$

$$B_0=\frac{{\int }_{\Omega }L(\theta )\pi (\theta )}{L({\theta }_0)},$$

$$\hat{\sigma }(c)={\hat{\sigma }}_{\mathrm{SM}}(1+{\delta }_{\mathrm{SM}})+c{\hat{\sigma }}_{\mathrm{int}}(1+{\delta }_{\mathrm{int}})+c^2{\hat{\sigma }}_{\mathrm{BSM}}(1+{\delta }_{\mathrm{BSM}})$$

$$\hat{\sigma }(c)={\hat{\sigma }}_0(c)(1+\delta \Delta (c)),$$

$$\bar{L}(c)=\int d\delta \pi (\delta )\frac{e^{-n(c)}n(c{)}^{n(c^{\prime })}}{\Gamma (n(c^{\prime })+1)},$$

Full text

I Introduction

In modern High Energy Physics, the use of large datasets has become commonplace. In two areas in particular, Particle Physics and Cosmology, the forefront of discoveries and characterization of new phenomena relies on extraction of information from complex datasets produced by experiments like Planck Ade et al. (2016) and the LHC de Florian et al. (2016). In both fields, a precise theoretical paradigm is used to interpret the data ($\Lambda$CDM and SM, respectively) and the search for new phenomena depends then on identifying subtle deviations within the data, often relying on machine learning techniques. For example, the discovery rare SM processes, like mono-top Aaltonen et al. (2009) and Higgs decays to tau-leptons Aad et al. (2016), has been achieved using this methodology.

On the theoretical side, these multivariate techniques obscure the physical understanding of which variables drive the analysis, making the re-interpretation of results very difficult and in general hindering the public use of the data. Yet more detailed information, in particular differential rates, is required to advance the programme of searching for a new paradigm beyond the standard one. For example, the use of differential information on Higgs production Ellis et al. (2014) has proven key to pushing the limits of understanding the impact of possible new phenomena in the Higgs boson properties.

In this paper we investigate the advantages and limitations of multidimensional shape information in searching for new physics and present two case studies, the new physics search in the context of the SM Effective Field Theory (SMEFT) and the characterization of the quantum numbers of a new resonance. These case studies, together with the material collected in the Appendix, can readily be used as guidelines on how experiments could provide data and how theorists would use it.

II Information content of multidimensional differential rates

In this section we study the information content of differential rates and their use in discovery and model comparison. The statistical formalism and details are provided in App. A.

The information content of a likelihood function with respect to a parameter $\theta$ is measured by the observed Fisher information $I_{\theta}[L]\equiv-\partial^{2}_{\theta}\log L$ 111The expected Fisher information has been recently advocated in Ref. Brehmer et al. (2016) as a way of visualizing the information content of differential distributions and the expected constraints on EFT operators. The results presented in this reference are orthogonal to those discussed in the present Letter, where we focus on the relative power of multidimensional analyses and use the observed Fisher information only for qualitative arguments. . The likelihood functions we focus on arise from event counting and can always be factored as

[TABLE]

where $L_{\rm tot}(\theta)$ contains the information on the total rates and $L_{\rm shape}(\theta)$ contains information on the shape of differential distributions 222This is true for Poisson statistics, but also in presence of systematic uncertainties, as these can always be split into as subset affecting only $L_{\rm tot}(\theta)$ and a subset affecting only $L_{\rm shape}(\theta)$. . The information content from total event number and from the shape ($I_{\rm tot}\equiv I[L_{\rm tot}]$, $I_{\rm shape}\equiv I[L_{\rm shape}]$) are thus independent from each other. $I_{\rm shape}$, therefore, could be arbitrarily large with respect to $I_{\rm tot}$, i.e. the amount of information contained in the shape could dominate over the amount of information contained in the total rate. It is thus fully justified to systematically take into account the shape information on top of the total rate information.

The information content of $L_{\rm shape}$ with respect to the dimensionality of the differential rate distribution is slightly more subtle. For concreteness let us consider the case of one kinematic variable (“1D”) versus two kinematic variables (“2D”). The variables are labelled $a$ and $b$. If the $a$ and $b$ variables are totally correlated, one has $I_{\rm shape}^{{\rm 2D}}=I_{\rm shape}^{{\rm 1D},a}=I_{\rm shape}^{{\rm 1D},b}$, and there is no gain in going from 1D distributions to 2D information. On the other extreme, if the two variables $a,b$ provide uncorrelated information, the likelihood factorises and the total information is given by $I_{\rm shape}^{{\rm 1D},a}+I_{\rm shape}^{{\rm 1D},b}$. This is the maximum information possible, thus one obtains that as

[TABLE]

the potential gain from 1D to 2D cannot be arbitrarily large.

The gain from 1D to 2D is maximal when the two 1D distributions are of the same order of magnitude, $I_{\rm shape}^{{\rm 1D},a}\sim I_{\rm shape}^{{\rm 1D},b}\equiv I_{\rm shape}^{{\rm 1D}}$, with the maximal value for $I_{\rm shape}^{{\rm 2D}}$ given by twice $I_{\rm shape}^{{\rm 1D}}$. In the rest of this paper, we will often refer to the information gain obtained from using a 1D differential distribution to using a 2D differential distribution as the “1D/2D” gain.

Pursuing in our general considerations, let us evaluate the impact of the various pieces of information discussed above in the subsequent statistical analyses. In order to proceed, we need to consider statistical tests. We adopt the framework of Bayesian statistics, which allows a unified treatment of discovery and model comparison For hypothesis testing, the relevant quantity to use is the Bayes factor (see app. A for definitions).

We assume that the likelihood for each hypothesis can be approximated by a Gaussian with respect to the parameter of interest $\theta$. This limit tends to occur once at least $O(10)$ events are collected 333 This behaviour has first been described in Ref. Wald (1943).. The likelihood function then takes the form

[TABLE]

where $\bar{\theta}$ is the value of $\theta$ preferred by the data, $I$ is the Fisher information for $\theta$, and the constant $L_{\rm max}$ encodes the information about goodness-of-fit between the hypothesis and the data. In this Gaussian limit, the Bayes factors exhibit simple expressions with respect to the Fisher information(s). Moreover, the Fisher information depends *linearly * on the observed total event number $n_{\rm obs}$ to a good approximation. Hence we obtain $I_{\rm shape}^{{\rm 1D}}=\alpha^{{\rm 1D}}\,n_{\rm obs}$, $I_{\rm shape}^{{\rm 2D}}=\alpha^{{\rm 2D}}\,n_{\rm obs}$, with $\alpha^{{\rm 2D}}\leq\alpha^{{\rm 1D},a}+\alpha^{{\rm 1D},b}$. Note that the $\alpha^{{\rm 2D}}$ information coefficient can be at best $\alpha^{{\rm 1D},a}\sim\alpha^{{\rm 1D},b}$. This direct link of Fisher information to the event number is crucial to concretely quantify the impact of the various likelihoods.

To characterize discovery, we introduce the discovery Bayes factor, which compares a model hypothesis with a free parameter $\theta$ with the same hypothesis restricted to $\theta=\theta_{0}$ (see App. A). The discovery Bayes factor is given by

[TABLE]

with $\alpha=\alpha_{\rm tot}+\alpha_{\rm shape}$. Note that the constant $L_{\rm max}$ does not appear in this expression. The first term in Eq. (4) encodes the comparison of central values while the second term encodes prior information. This second term becomes quickly negligible once $n_{\rm obs}$ increases. Comparing the discovery Bayes factor for 1D and 2D distributions, and assuming $\bar{\theta}_{\rm 1D}\sim\bar{\theta}_{\rm 2D}$, we get that

[TABLE]

The bound is saturated when the 2D information is maximal, $\alpha_{\rm shape}^{\rm 2D}=\alpha_{\rm shape}^{{\rm 1D},a}+\alpha_{\rm shape}^{{\rm 1D},b}$, and for $\alpha_{\rm tot}\ll\alpha_{\rm shape}$. Finally, the information gain from 1D to 2D would be maximal when $\alpha_{\rm shape}^{{\rm 1D},a}\sim\alpha_{\rm shape}^{{\rm 1D},b}\equiv\alpha_{\rm shape}^{{\rm 1D}}$, in which case we obtain that $\log B_{0}^{\rm 2D}$ could be at most twice $\log B_{0}^{{\rm 1D}}$.

This bound on the 2D Bayes factor can be easily translated in terms of sample size and evidence strength. In terms of sample size, the bound can be translated using the fact that the 1D/2D gain amounts to at most doubling the $n_{\rm obs}$ from the 1D case. In terms of strength of evidence (see Jeffreys’ scale), we observe that moving from 1D to 2D can lead to a shift of at most one step in evidence strength. For instance, if the 1D Bayes factor would give moderate evidence ($\log B_{0}=2.5$), the 2D Bayes factor could at most reach strong evidence ($\log B_{0}=5$).

So far we have discussed how the 1D/2D information gain is bounded in the scope of a discovery. But what about model discrimination? Approximating the likelihoods as Gaussians in both hypotheses $H_{1}$, $H_{2}$, the Bayes factor comparing $H_{1}$ to $H_{2}$ reads

[TABLE]

Note that the structure of this Bayes factor is different from the discovery Bayes factor. The first term encodes the relative goodness-of-fit of the models with respect to data, whereas the second term in Eq. (6) is a ratio of Fisher information, and should be understood as a measure of the relative fine-tuning of the two models, see Fichet (2012). This second, “naturalness” term is independent of $n_{\rm obs}$. In contrast, the ratio of maximum likelihoods depends in general on $n_{\rm obs}$, as goodness-of-fit is in general different in both hypotheses. In fact, in the large sample limit, one expects

[TABLE]

where $\beta$ a positive or negative constant. The case $\beta>0$ corresponds to the $H_{1}$ model being a better fit than $H_{2}$, and conversely. The absolute value of $\log\left(L_{{\rm max},1}/L_{{\rm max},2}\right)$ is thus expected to grow with $n_{\rm obs}$, reducing the relative impact of the naturalness term.

We can now compare a Bayes factor based on a 1D distribution, $\log B_{12}^{\rm 1D}$, with a Bayes factor based on a 2D distribution, $\log B_{12}^{\rm 2D}$. Neglecting the naturalness terms (which are however different for 1D and 2D), we are left with comparing the goodness-of-fit terms of the 1D and 2D cases, roughly given by $\beta_{\rm 1D}n_{\rm obs}$ and $\beta_{\rm 2D}n_{\rm obs}$. We have found no bound on the $\beta_{\rm 2D}/\beta_{\rm 1D}$ ratio based on general information considerations. This suggests that the 1D/2D information gain can be arbitrarily large in case of model comparison. This is related to the fact that for model comparison the goodness-of-fit matters, while in case of discovery it does not, i.e. $B_{0}$ does not depend on $L_{\rm max}$.

III Case studies for discovery and characterization

In this section we evaluate the impact of multidimensional differential information using realistic LHC simulations and apply the procedure to various hypotheses of physics beyond the SM (BSM).

III.1 Simulation and analysis setup

To simulate the conditions in the LHC for different model hypothese, we use FeynRules to implement the BSM models, and the UFO Degrande et al. (2012) output to interface with MadGraph5 aMC@NLO platform Alwall et al. (2014). The parton events are then passed through Pythia Sjostrand et al. (2006) for parton-showering and hadronization. Finally, the hadrons are reconstructed via the anti-$k_{T}$ algorithm Cacciari et al. (2008) with an R-parameter set to 0.4 using the FastJet Cacciari et al. (2012) interface of MadAnalysis 5 Conte et al. (2013). A jet is tagged as arising from a $b$-quark when a B-hadron is present within a cone of radius $R=0.4$ centred on the jet momentum direction. A private pyROOT script has been developed in order to automatize and monitor the whole analysis in the framework of MadAnalysis.

As our focus is on evaluation and comparison of analyses based on future data prospects, we have introduced “projected” data in our likelihoods, see App. A for details. No statistical fluctuations in the projected data are assumed, and these are thus directly given by the expected rates.

III.2 Case I: CP-violating and -conserving SMEFT

In the scenario where new particles are too heavy to be produced on-shell at the LHC, their observables effects are better described by a low-energy effective theory, in the so-called SMEFT framework. For our case study we assume the presence of two characteristic dimension-six operators,

[TABLE]

where the operators $O_{HW}$ and $\tilde{O}_{HW}$ are defined as

[TABLE]

Here $\Phi$ is the Higgs doublet and $W_{\mu\nu}$ the $SU(2)_{L}$ field strength. The Wilson coefficients $c_{i}$ are normalized following the SILH basis conventionsContino et al. (2013) and their current bounds can be found in Refs. Ferreira et al. (2016); Ellis et al. (2015). Note we use the implementation of these operators provided in Ref. Alloul et al. (2014).

As an example of the use of differential information, we consider Higgs production in the Vector Boson Fusion (VBF) mechanism, and generate samples of $450$K events. Basic selection cuts require the presence of two jets with a transverse momentum $p^{j}_{T}>20$ GeV, pseudo-rapidity $|\eta_{j}|<4.5$, as well as typical VBF cuts: the dijet invariant mass is required to be larger than $400$ GeV and the jet separation in pseudo-rapidity to be above $2.8$. The analysis selects two high-momentum jets $j_{1},j_{2}$ and two photons $\gamma_{1},\gamma_{2}$ from the Higgs decay. The indexes $1$ and $2$ denote the leading and sub-leading particles.

In order to determine the differential rates for arbitrary values of the effective operator coefficients, we use the reconstruction method described in Brenner et al. (2016); ATL (2015); Fichet et al. (2016) - which has been dubbed “morphing” in experimental references. The optimal version of the reconstruction method has been described in Fichet et al. (2016) and is used in our analyses. The reconstruction provides estimates of the various components of the rate on a given bin, $\hat{\sigma}_{r}(c)=\hat{\sigma}^{\rm SM}_{r}+c\hat{\sigma}^{\rm int}_{r}+c^{2}\hat{\sigma}^{\rm BSM}_{r}$ with $c=c_{HW},$ $\tilde{c}_{HW}$. An important subtlety related to the estimation of the interference component in regions with low event rates is described in App. B

The projected data are directly given by the $\hat{\sigma}(c^{\prime})$ event rates, where $c^{\prime}$ is the value operator coefficient assumed to be present in these projected data. The fact that we use the same rates $\hat{\sigma}(c)$ for both projected and expected rates leads to an interesting simplification. It turns out that the main Monte Carlo uncertainties cancel out from the likelihood, leaving the maximum likelihood rigorously unchanged (see app. C). Rather, the uncertainties are only changing the Fisher information part, and more generally the likelihood line-shape. This simplification implies that in practice, the number of Monte Carlo events to perform the simulation needs to be only mildly larger than the nominal number of events. Having for example $n^{\rm MC}_{r}>3n^{\rm obs}_{r}$ gives a systematic uncertainty of $\sim 33\%$ on the Fisher information, and thus of $\sim 16.5\%$ on the projected statistical uncertainty.

Having described all the aspects of our simulations and likelihoods, let us proceed with the data analysis focusing on 1D and 2D differential distributions. We have tested the constraining power of a set of basic kinematic variables including transverse momenta, azimuthal angles and longitudinal rapidity differences between final state objects. Throughout our study we found that the pair of variables with the best 1D/2D gain for discovery are $p_{T}^{\gamma_{1}}$, $\Delta\phi_{\gamma\gamma}$, hence we present the analysis with respect to these two variables, see Fig. 1. Note, though, that the analysis does not include detector effects which could change the set of optimal variables.

We first compute the posterior distributions for $L_{\rm tot}$, $L_{\rm 1D}$, $L_{\rm 2D}$, assuming $c^{\prime}_{HW}=\tilde{c}^{\prime}_{HW}=0$ in the projected data. The preferred regions are shown in Fig. 2. We can see that taking into acccount the 1D distributions is crucial in order to lift the degeneracy in the $c_{HW}-\tilde{c}_{HW}$ plane. In contrast, the gain from 1D to 2D differential information turns out to be mild.

Assuming that the $O_{HW}$ operator with $c^{\prime}_{HW}=-0.01$ is present in the data, we compute the discovery Bayes factor for $\ O_{HW}$ as a function of the sample size, as shown in Fig. 3. A mild gain between $B_{0}^{\rm tot}$ and $B_{0}^{\rm 1D}$ and between $B_{0}^{\rm 1D}$ and $B_{0}^{\rm 2D}$ is observed. That the two 1D Bayes factors have almost the same value is apparently a mere coincidence. The 1D/2D gain for this pair of kinematic variables is the best we found among the kinematic variables considered. A positive result is also obtained when assuming the existence of the operator $\tilde{O}_{HW}$ instead of $O_{HW}$.

Still assuming the presence of the $O_{HW}$ in the data, we then use a Bayes factor comparing the $c_{HW}\neq 0$ hypothesis with the $\tilde{c}_{HW}\neq 0$ hypothesis, see Eq. (6). The result is shown in Fig. 4. We observe that the 1D/2D gain in this case is much larger than for discovery. For example we can see that the 1D/2D gain in sample size is about $90\%$, which corresponds to almost doubling the sample size 444We define the “1D/2D gain in sample size” as $n^{\rm 2D}|_{B=150}/n^{\rm 1D}|_{B=150}-1$.. For comparison, for the discovery, the gain in sample size is of $20\%$.

We should stress that certain kinematic variables such as $m_{jj}$ have a better discriminating power than the variables we consider but are not as good for discovery, hence we present results based on the $p_{T}^{\gamma_{1}}$, $\Delta\phi_{\gamma\gamma}$ variables in order to have a direct comparison with the discovery Bayes factor. Nevertheless, the large 1D/2D gain persists for these other combinations of variables, the pair $m_{jj}-p_{T}^{\gamma_{1}}$, for instance, has 1D/2D gain in sample size is about $\sim 100\%$ .

III.3 Case II: Testing the spin of a resonance

We consider now the discovery of a new resonance with either spin zero or two, and how our analysis would help on the characterization of the resonance using a final state of a pair of $Z$ bosons further decaying leptonically, $pp\rightarrow\phi+X\rightarrow 2l^{+}2l^{-}+X$. The spin-0 and spin-2 resonance behaviour is simulated using existing FeynRules models. Samples of $450$K events were generated, following the hadronization procedure previously described. The two pairs of opposite-sign leptons are required to have an invariant mass close to the $Z$ boson mass, $75\mathrm{GeV}<m_{ll}<105\mathrm{GeV}$, and are sorted by their transverse momentum, with $l_{1}$ being the hardest lepton. We chose $p_{T}^{l_{1}}$, $\Delta\phi_{ll}$ as kinetic variables for analysis. Widths and production rates of the two resonances are assumed to be the same, so that only differential distributions may be used to distinguish the spin of the resonance.

Following the same approach as for Case I, we assume that the projected data arises from a spin-2 resonance, and compute the Bayes factor $B_{20}$ comparing the spin-2 hypothesis to spin-0 hypothesis. The 1D/2D gain in sample size is found to be $\sim 50\%$. Similar results are obtained when assuming a spin 0 resonance and computing $B_{02}$. We observe thus again a substantial gain of information when using the 2D distribution instead of individual 1D distributions.

IV Conclusions

In the view of future new physics searches and characterization at the LHC, we investigate the impact of multidimensional differential rates in typical LHC analyses. Through general observations based on Bayes factors and Fisher information, we find that in the occurrence of a discovery, the gain from using 2D differential distributions instead of 1D is fundamentally bounded. In contrast, for model discrimination, no such bound is found, thus the gain from 1D to 2D can be much higher. To illustrate these features and show realistic values of the 1D/2D information gain, we study two new physics scenarios: operators from the SMEFT, and bosonic resonances.

We carried out discovery and discrimination tests in the VBF channel in presence of CP-even and CP-odd operators. The best 1D/2D gain for discovery is found for the combination of variables ($p_{T}^{\gamma_{1}},\Delta\Phi_{\gamma\gamma}$). As expected, the 1D/2D gain for CP discrimination is found to be much larger than for discovery. This observation also holds for various other choices of variables. In the presence of a heavy bosonic resonance, we evaluate the discrimination power of spin-0 versus spin-2 using the ($p_{T}^{l_{1}},\Delta\Phi_{ll}$) variables, and observe a 1D/2D gain of about $50\%$. Note, though, that the procedure of adding more differential information saturates as the kinematic information in a given final state is limited. Hence the gain from 2D to higher-dimensional distributions will be restricted due to correlations among of the variables involved.

All details needed to reproduce our analysis are provided, and important subtleties generally present in these analyses are pointed out. First, in the reconstruction method of the differential rates in the SMEFT, we point out that in the low-event regions of phase space the Cauchy-Schwarz bound on the interference can be violated by the large systematic uncertainties, resulting in unphysical results. Second, we find that when using projected data, a cancellation between the leading uncertainties of expected and projected rates naturally occurs, implying that the maximum likelihood would remain unchanged and hence the MonteCarlo sample size would only have to be mildly larger than the nominal sample size to provide meaningful results.

We hope this study serves as a guide to experiments to provide differential information to theoretical collaborations, and as how to use this information for model discrimination.

Acknowledgements

The S.F. work was supported by the São Paulo Research Foundation (FAPESP) under grants #2011/11973 and #2014/21477-2. V. S. is supported by the Science Technology and Facilities Council (STFC) under grant number ST/L000504/1. F.F. is supported by Coordenação de Aperfeiçoamento de Pessoal de Nível Superior (CAPES).

Appendix A Statistical basics

In this Appendix we set the notation of the statistical analysis. We denote phase space by $\cal D$, and consider a binning of $\cal D$ in $d$ dimensions. The bins are set along a dimension $i\in(1\ldots d)$ and labelled by $r_{i}$, with the coordinates $(r_{1},\ldots,r_{d})$ of a bin denoted $r$, and the associated piece of phase space ${\cal D}_{r}$.

The observed event number in the bin $r$ is denoted $\hat{n}_{r}$, and the expected event number for a given value of the underlying parameter $\theta$ is denoted $n_{r}(\theta)$. Total number of observed events is $\hat{n}=\sum_{r}\hat{n}_{r}$ and total number of observed events is $n=\sum_{r}n_{r}$.

For further convenience one also introduces the $d$-dimensional density of expected events $f_{X}(x)$, where $X=(X_{i})$ denotes the set of binned variables. $f_{X}(x)$ is simply the differential event rate normalized by the total event rate. The expected event number in a bin $r$ is then given by

[TABLE]

A.1 Likelihood

The likelihood function $L$ is defined as the conditional probability of obtaining the observed data given a hypothesis, taken as a function of this hypothesis. For a hypothesis $H$ with a set of parameters $\theta$,

[TABLE]

The likelihood function can be defined up to an overall constant factor.

The events counted in each of the bins are statistically independent, hence the likelihood factorises as

[TABLE]

The event number in every bin follows a Poisson statistics, so that the likelihood function in the bin $r$ is given by

[TABLE]

For a given integrated luminosity $\cal L$, $n_{r}(\theta)$ is given by the event rate on the bin, $n_{r}(\theta)={\cal L}\sigma_{r}(\theta)$.

The likelihood can be formally factored in a Poisson term $L_{\rm tot}$ containing the information about the total rate and a term $L_{\rm shape}$ containing the information about the shape of the differential distribution, so that $L(\theta)=L_{\rm tot}(\theta)L_{\rm shape}(\theta)$ with

[TABLE]

[TABLE]

A.2 Credible regions and hypothesis testing

We adopt the framework of Bayesian statistics 555The Bayesian framework is consistent with the “likelihood principle”, which states that all experimental information is encoded in the likelihood function. This is not the case of, for example, frequentist $p$-values.. The model parameters are given an a-priori probability density $\pi(\theta)$, called “prior”, that can encode both subjective and objective information. The “posterior” density is defined as $p(\theta)\propto L(\theta)\pi(\theta)$, it provides the preferred regions of $\theta$ ones data are taken into account. The shape of the prior becomes irrelevant once enough data are accumulated, i.e. when the posterior is data-dominated.

A so-called $1-\alpha$ credible region of highest density is defined by the domain $\Omega^{1-\alpha}=\{\theta\,|\,p(\theta)>p_{1-\alpha}\}$, where $p_{1-\alpha}$ is determined by the fraction of integrated posterior

[TABLE]

$\Omega$ being the whole parameter space. We will use the credible regions associated with $1-\alpha=\{68.27\%\,,\,95.45\%\,,\,99.73\%\}$ 666Note that confidence regions are not uniquely defined, but the method of highest density is the most commonly used and arguably the most natural..

Comparison between two hypotheses $H_{0}$ and $H_{1}$ is done by means of the Bayes factor

[TABLE]

where the $\pi_{0,1}$ are the priors for hypotheses $H_{0,1}$ respectively. The Bayes factor is interpreted using the Jeffreys’ scale Trotta (2008), which associates weak, moderate and strong evidence in favour of $H_{0}$ to the threshold values $\log B_{01}\sim 1,2.5,5$ (i.e. $B_{01}\sim 3,12,150$ ).

The Bayes factor framework can be used in the context of new physics searches. In order to assess that the data favour a hypothesis where a parameter $\theta$ is different from a given value $\theta_{0}$ one has to compare the $H_{1}$ hypothesis to $H_{0}\equiv H_{1}|\theta=\theta_{0}$. In the context of effective operators, $H_{1}$ can for instance be the SM deformed by higher dimensional operators (the SMEFT), while $H_{0}$ is the SM. Defining $B_{0}\equiv 1/B_{01}$, we have

[TABLE]

that we refer to as the discovery Bayes factor. The test assesses that $\theta\neq\theta_{0}$ for $B_{0}>1$, using the thresholds given above.

A.3 Projected data

In order to evaluate the sensitivity of a future analysis, measurement, or experiment, one can rely on imaginary, speculative data. That is, instead of introducing actual observed data in the likelihood Eq. (11), one can instead introduce speculative data coming for instance from a simulation of the experiment. We refer to these as projected data.

An important subtlety, well discussed in Cowan et al. (2011), is that an assumption has to be made on the statistical fluctuations present in the projected data. Along this paper, we will simply consider the case where no statistical fluctuations are present in the projected data. A dataset satisfying this condition is sometimes referred to as an “Asimov” dataset Cowan et al. (2011).

The event numbers in the projected dataset assuming no statistical fluctuations and the presence of an operator with coefficient $c^{\prime}$ are then simply given by $\mathscr{L}\sigma_{r}(c^{\prime})$. In practice, these rates have to be estimated by MonteCarlo simulations, just like the expected ones.

Appendix B Violation of the interference’s upper bound

The interference component in the true rate $\sigma_{r}(c)=\sigma^{\rm SM}_{r}+c\sigma^{\rm int}_{r}+c^{2}\sigma^{\rm BSM}_{r}$ satisfies a bound based on Cauchy-Schwarz inequality, $|\sigma_{\rm int}|<2\sqrt{\sigma_{\rm SM}\sigma_{\rm BSM}}$. From a concrete viewpoint, this bound prevents the cross-section for becoming negative for any value of $c$. Now, on bins featuring few events, the uncertainty is large enough so that there is a non-negligible probability that the Cauchy-Schwarz bound be violated, i.e. $|\hat{\sigma}_{\rm int}|<2\sqrt{\hat{\sigma}_{\rm SM}\hat{\sigma}_{\rm BSM}}$ can be false. This then results in a negative rate on the bin for some interval of $c$.

If the simulation has been done with enough events compared to the nominal event number, for any bin and any $c$, these rare events regions have a negligible impact on the subsequent analysis. Violation of the Cauchy-Schwarz constitutes then only a practical obstruction, and it is convenient to simply remove the bins on which the Cauchy-Schwarz bound is not respected 777Another way of reconstructing the differential rate is to use the recently implemented option in Madgraph 5 that removes the $\sigma_{\rm BSM}$ term from the amplitudes. We found a large amount of bins violating the Cauchy-Schwarz bound when using this prescription, and thus did not pursue with this method..

Appendix C Cancellation of systematic uncertainties

When confronting the expected event numbers to observed event numbers, one has to make sure that the systematic uncertainty on the expected event numbers has to be negligible with respect to statistical uncertainty, for every bin and for both SM and BSM hypotheses (for any relevant value of $c$, for example). In practice, this means that the number of MC events has to be larger than the number of observed events in all of these situations.

However, when the analyses involve projected data instead of actual data, a very nice and useful feature appears. It turns out that, provided one uses the same estimates for the projected rates and for the expected ones, the respective systematic uncertainties present in these rates will approximately cancel each other.

Let us detail how this occurs. The uncertainties on the reconstructed rates take the form

[TABLE]

where the $\delta$’s are the nuisance parameters -which are in general correlated (see Fichet et al. (2016) for their correlation matrix). Let us first adopt the Gaussian limit for simplicity. The likelihood using projected data takes the form $\exp(-\mathscr{L}(\hat{\sigma}(c^{\prime})-\hat{\sigma}(c))^{2}/2\hat{\sigma}(c))$, in which it is clear that the SM uncertainty in the numerator cancels out exactly, and the other ones are suppressed by $c-c^{\prime}$ and $c^{2}-(c^{\prime})^{2}$. As a result, the maximum of the likelihood remains unchanged and the main effect of the uncertainty is to distort the Gaussian and the Fisher information.

But these features turn out to be rigorously valid beyond the Gaussian limit, for Poisson likelihoods with any number of events. To see this, one first combines the three nuisance parameters into a single one. This operation is rigorously defined, and has already lead to useful developments in the context of LHC analyses Fichet and Moreau (2016); Fichet (2016); Arbey et al. (2016). After combination, the event rate takes the form

[TABLE]

where $\delta$ is the nuisance parameter and $\Delta(c)$ controls the relative magnitude of the uncertainty. The marginal Poisson likelihood with projected data is

[TABLE]

where $n(c)={\mathscr{L}\hat{\sigma}(c)}$ and $\pi(\delta)$ is the prior for $\delta$. Computing the derivate of $\bar{L}(c)$, it turns out that the maximum of $\bar{L}(c)$ still occurs for $c=c^{\prime}$, in spite of the deformations induced by the uncertainties, and for any $\pi(\delta)$.