Statistical issues in searches for new phenomena in High Energy Physics

There is a mathematical approach, developed largely by Russian mathematicians, which provides a series of axioms which probability must obey. They include concepts such as probability being represented by a number in the range zero to unity inclusive; zero probability means the relevant event never happens; the probability of something happening is 1 minus the probability of it not happening; etc. These axioms provide the necessary theoretical underpinning to the Bayesian and frequentist approaches to probability, but gives no real understanding of what it is.

The frequentist approach is to define the probability of a process as the limit as N tends to infinity of the fraction of times it occurs in N essentially identical independent trials. Since the definition involves a large number of repetitive trials, it does not apply to a unique event or to the value of a physical quantity. The following are examples of cases where a frequentist would not assign a probability:

• Rain in Manchester:Frequentists would not ascribe a probability about whether it was raining in Manchester on the day of the Brexit vote. They would say either it was or it was not, and although we personally do not know, we could find out by looking at meteorological records. There is no way we can have a large number of repetitions of that day, to produce an ensemble for us to see the fraction of days there was rain. Furthermore frequentists would similarly refuse to provide a probability estimate, even if 'Manchester' was replaced by 'The Sahara Desert'.
• Astronauts to Mars:Again, frequentists would not assign a probability to whether the first astronauts to Mars will return to Earth alive, because there is only one first set of astronauts going there. This does not mean it is an uninteresting question, especially if you are one of those first astronauts, but in that case do not ask a frequentist for an assessment of the probability.
• Physical constants:Similarly, frequentists will not speak about probabilities for numerical values of physical constants, or whether a particular theory is true or false. Examples would be whether the amount of dark matter in the Universe is above 25% of the critical density; or whether SUSY particles of mass below 1 TeV exist. At present, we do not know, but these are not issues where we can have a series of repeated Universes with differing answers, for us to see in what fraction the statement turns out to be true.

In contrast Bayesians regard probability as a personal issue, expressing the degree to which one believes a statement. Thus people with different information available or with different beliefs can assign different numerical values to a probability. For example, if I toss a coin, you would be likely to decide that the probability of heads is 0.5. However, I might have cheated and seen that it is tails, and so for me the probability is zero. Alternatively, perhaps I had a very quick glance so that you would not see that I was cheating, and I think that it is tails but am not completely sure, so maybe I assign a probability of 0.2 to heads.

Given that Bayesian probability is a measure of personal belief, it does not require an extended series of trials, and so can be applied to a single event, or to the numerical value of a physical constant. The decision of what numerical value to assign a probability is determined by the concept of a fair bet: if you think that there is a 20% probability of something specific happening, you should be prepared to offer 4 to 1 odds and allow a person to bet against you in whichever direction they want. If you assign the probability incorrectly, you are liable to lose money.

Bayes' theorem plays a central role in Bayesian analysis. It starts with the probability $P(A\ {\rm{and}}\ B)$ of A and B both happening. As can be readily verified, this satisfies

$$\begin{eqnarray}&&P(A\ {\rm{and}}\ {B})=P(A| {B})P({B})=P({B}| A)P(A),\end{eqnarray} \tag{ 2.1 }$$

where $P(A| B)$ is the conditional probability of A happening, given the fact that B has occurred; and the last part of the equation follows because $P(A\ {\rm{and}}\ B)$ must be the same as $P(B\ {\rm{and}}\ A)$. If we divide the second part of equation (2.1) by P(B) we obtain

$$\begin{eqnarray}&&P(A| B)=P(B| A)P(A)/P(B).\end{eqnarray} \tag{ 2.2 }$$

This is Bayes' theorem, and relates the conditional probabilities $P(A| B)$ and $P(B| A)$.

Not everyone appreciates that these probabilities are in general not equal⁵ . This sometimes manifests itself in statements about the probability of data, given the 'No New Physics' hypothesis, being incorrectly interpreted as the probability of 'No New Physics', given the data.

Bayes' theorem is acceptable not only to Bayesians, but also to frequentists, provided that all the probabilities are in accord with frequentist views of probability. The controversy begins when Bayesians replace A by 'hypothesis' or 'parameter value' and B by 'data'. Then, for example, Bayes theorem becomes

$$\begin{eqnarray}&&P({\rm{param}}| {\rm{data}})\propto P({\rm{data}}| {\rm{param}})P({\rm{param}}).\end{eqnarray} \tag{ 2.3 }$$

Here P(param) is the Bayesian prior which is the probability density that we assign to different parameter values before we perform our experiment; $P({\rm{data}}| {\rm{param}})$ is the likelihood function, encapsulating what we can say about the parameter as a result of our data; and $P({\rm{param}}| {\rm{data}})$ is the Bayesian posterior, combining the information from previous knowledge with that from our data. The problem is that equation (2.3) involves probability densities for parameter values which, as explained in section 2.1.1, are not acceptable for frequentists. Bayesians would retort that frequentists use too restrictive a definition of probability.

In addition, the choice of Bayesian priors can be a problem, especially in situations where we are sensitive to parameter values in a range where no or little previous information exists. The idea of choosing a constant prior, in order not to favour any particular value of the parameter μ, is fallacious; as well as being ignorant about μ, we could equally well say that we know nothing about ${\mu }^{2},1/\mu ,{ln}\mu$, etc, and priors which are constant in these different variables are not the same.

The frequentist approach to parameter determination is very different, and uses only $P({\rm{data}}| {\rm{parameter}})$, the probability to observe the data given a specific parameter value, and avoids $P({\rm{parameter}}| {\rm{data}})$. Often, the data are summarised by a single quantity known as the 'test statistic' (see section 2.5.1) and the probability to obtain the observed test statistic value is used for parameter determination. For example, consider the case of estimating the temperature T of the fusion reactor at the centre of the Sun by using an estimate of the solar neutrino flux ϕ from one month's data from a large solar neutrino detector. The Neyman construction uses $P(\phi | T)$ at a given value of T to choose a region⁶ in ϕ for which the fraction of outcomes in that region is a pre-determined level $1-\alpha$, for example 90% (see figure 1). This is repeated for all values of T to produce the shaded confidence band, showing the likely values⁷ of the data for any value of the parameter of interest. Then we collect one month's data, from which we deduce a flux ${\phi }_{d}$. The intersection of a vertical line at ${\phi }_{d}$ with the confidence band then gives the frequentist range ${T}_{l}\leqslant T\leqslant {T}_{u}$ for the parameter T at the chosen confidence level.

Zoom In Zoom Out Reset image size

The Feldman–Cousins approach [1] is a frequentist method which uses an ordering rule for deciding which of the infinitely many possible data regions of probability $1-\alpha$ is chosen for the confidence band. It has many useful properties and is commonly used in HEP.

For the determination of a parameter μ, Bayesians and frequentists can end up with a statement of the form

$$\begin{eqnarray}&&{\rm{Prob}}\ ({\mu }_{l}\leqslant \mu \leqslant {\mu }_{u})=1-\alpha \end{eqnarray} \tag{ 2.4 }$$

(sometimes even with the same values of ${\mu }_{l},{\mu }_{u}$ and α), but their interpretations are very different. For Bayesians, the measurement has been performed once, and from the data, ${\mu }_{l}$ and ${\mu }_{u}$ have been determined and are constants. Although the physics parameter μ presumably has a unique value, this is unknown and they ascribe a probability (or credible) distribution to it. Then a fraction $1-\alpha$ of this distribution is in the range ${\mu }_{l}$ to ${\mu }_{u}$. In contrast, frequentists regard ${\mu }_{l}$ and ${\mu }_{u}$ as part of an ensemble of ranges for μ that would be obtained if the measurement were to be repeated many times. Then equation (2.4) is a statement that a fraction $1-\alpha$ of these ranges should contain the true value of the parameter.

This situation is where we are dealing with the result of the analysis e.g. Do we have convincing evidence for the existence of a SUSY particle? Is the parity of the Higgs boson positive (as expected) or negative? Is the hierarchy of neutrino masses normal or inverted? etc.

The different possible outcomes are summarised in table 1. A notable difference from event selector situations (see section 2.3.2) is that here we can make no decision regarding the two hypotheses. Also an incorrect answer is now more serious. In an Error of First Kind we reject H0 when it is in fact true. This results in a false discovery claim, which is extremely serious, and warrants a very stringent confidence level (e.g. the '$5\sigma$ criterion' for discovery claims). An Error of Second Kind means that we have missed a possible discovery, which is unfortunate, but is regarded as not so terrible as false discovery claims.

Hypothesis testing is used at an early stage of many analyses, where an attempt is made to select wanted events for further analysis, while having a strong rejection of the usually more copious background. Here the relevant hypotheses are 'signal' as H0 and 'background' as H1. An example could be that we want to select events containing b-quarks in a particular analysis. In these cases, there is no option of not making a decision; the events are either accepted or rejected. A multivariate method (for example: decision trees (DTs), neural networks or support vector machines) is usually adopted for this (see sections 2.3.3 and 4.5). Unlike the situation where our hypothesis testing is for the result of the experiment, incorrect classification is not a disaster, provided we can make a reliable estimate of the levels of loss of efficiency and background contamination, and correctly allow for them in the analysis. This will inevitably increase the statistical uncertainty on the result of the experiment as compared with an ideal situation of 100% efficiency and no background. An unfortunate possibility is that the result will be dominated by systematic effects associated with uncertainties in these estimates; or, even worse, that the systematics will be underestimated.

Multivariate methods usually allow the signal acceptance efficiency to be varied; thus with neural networks, the acceptance cut on the network's output can be altered. As the cut is loosened and the efficiency is increased (reducing Errors of First Kind), the background generally rises (more Errors of Second Kind) and the purity of the selected sample decreases. Some compromise is therefore needed in deciding where to make the cut. This is usually done by optimising the expected accuracy of the result (although some approaches build this optimisation into the multivariate method itself). Estimating the background is generally much more difficult than obtaining the signal efficiency, because the background could come from many sources, and it is often the hard-to-simulate tails of distributions that creep into the acceptance region⁹ . It may thus be prudent to reduce the systematic uncertainty at the expense of a somewhat larger statistical one.

The simplest, and one of the most common, multivariate techniques is a DT. In HEP, an event is typically summarised using a pertinent subset of the observables that are measured in the detector, or can be calculated from the measurements. A DT is essentially a sequence of thresholds which can be used to define a region of the observable space for which background events are rejected more strongly than those of the potential signal. A basic DT, with only two observables x and y, would look something like that shown in figure 2. The logic of this DT is shown in figure 2 (left). At the first node, the input value of x is compared to the threshold ${{x}}_{1}$. If the input value is smaller than ${{x}}_{1}$, the input value pair (x, y) is classified as background usually meaning an event with these input values would be rejected. If not, then the logic continues by following the direction indicated by the arrow to the second node. At the second node, the input value of y is compared to the threshold ${{y}}_{1}$ and a similar decision is made, this time based on the input value of y. Once the path ends, the input values (x, y) can be categorised into those which are selected as signal and those which are rejected as background. For each pair of input values, there is only a single, unique decision (or output), for example assigning −1 or +1 for background-like or signal-like, respectively. The values of x and y which end up as either signal or background are shown as the shaded blue or solid green regions of the observable space in figure 2 (right). The dashed lines which bound these regions are the values of the thresholds at each node in the DT.

Zoom In Zoom Out Reset image size

By following this chain of logic for every event, those which are more likely to be a signal can be separated from those which are most likely to be a background. The threshold values chosen for ${{x}}_{1},\,{{x}}_{2},\,{{y}}_{1}$ and y₂ depend on how the DT is optimised. Typically in particle physics, one wants to reject as much of the background as possible while ensuring that the efficiency for selecting signal is large. While of course, in any real data event we do not know whether the event is a signal or a background, we can use simulated events, or data from samples not used in the analysis, to calculate the background and signal efficiencies of the selection. The threshold values can be optimised to yield both a high signal effiency and background rejection. This is known as 'training' the DT. Typically the samples used to assess the performance of the DT are different from those used to train it, to ensure that it has not learned merely to distinguish signal and background fluctuations in the particular training samples used. In this example, there are only a few nodes; however more nodes can be added to the tree which can pick out finer structure of the signal and background in the multivariate observable space. The number of nodes in a DT is usually referred to as its depth.

In particle physics, it is often the case that rather than using a single DT, multiple DTs (sometimes referred to as forests) are trained and their outputs are combined, for example by averaging the individual outputs. For a large enough forest, this can yield a nearly continuous output so that rather than having the decision of either signal or background as the output, the output is a real number between −1 and 1 indicating that the input (x, y) is more 'background like' or more 'signal like' respectively. A final cut(s) can be placed on this value to define one or more categories of events which will be used in the statistical analysis. It should be obvious that the DT shown in figure 2 (left) is not unique in defining this particular division of the observable space. A number of similar DTs could lead to the same division of the signal and background like regions. Part of the optimisation in terms of the efficiency of the DT would be to construct the most simple one which gives the same outcomes.

In a similar way, the procedure by which a neural network can select a particular region of multi-dimensional observable space can be readily understood in terms of the responses of the nodes of the network to their inputs. As for DTs, there is no need to regard a neural network as a mysterious black box.

The procedure for using either a forest of DTs or a neural network is as follows:

• Define samples for the signal and background.
• Use these samples for training.
• Use independent samples to measure the signal and background efficiencies, and potentially optimise one or more thresholds to select events.
• Run real data through the forest or neural network.

For a thorough review of the algorithms used in neural networks and DTs, see [4, 5], respectively. Fortunately, there are several software packages available to train and use DTs and neural networks. Those commonly used in HEP are detailed in [6, 7].

A test statistic t is a way of summarising the data, and is used in deciding between H0 and H1. In a counting experiment (e.g. searching for dark matter), it could be simply the number of observed events, but in analyses using more information for each event, a likelihood ratio, perhaps profiled over systematic nuisance parameters, might well be more appropriate (as described in the penultimate paragraph of section 2.5.6).

Some searches for new physics involve precision measurements of a parameter (e.g. the anomalous magnetic moment of the muon), which is then compared with the SM prediction. Since we are looking for evidence of corrections from the effect of virtual new particles, we would regard this as an indirect approach, rather than searching for their direct production. We do not consider such analyses further in this article.

As in the previous section, here we are also involved in parameter determination, but in this case the parameter σ is directly related to the production of the new particle or effect.

Searches for new particles may define special regions of observable space, where the new effects give rise to an excess of events e.g.

$$\begin{eqnarray}&&\kappa =A\epsilon L\sigma +b,\end{eqnarray} \tag{ 2.5 }$$

where κ is the expected number of events; A and are the geometrical and kinematic acceptance, and the efficiency for observing signal events, respectively; L is the integrated luminosity describing the intensities of the incident beams and the running time; σ is the cross-section of the process we are interested in, and usually constrained to be non-negative; and b is the expected background. In reality there will be statistical and/or systematic uncertainties associated with $A,\epsilon ,L,$ and b but for simplicity we will regard them as precisely known. In the absence of a signal, the observed number of events n_obs will be Poisson distributed with mean b; because of statistical fluctuations, n_obs can be smaller than b. For a given n_obs, we might hope to find that σ was conclusively larger than zero, and hence claim a discovery. A more modest aim is to set an upper limit on possible values of σ i.e. larger values of σ are excluded because they would likely have resulted in significantly more events than we saw.

The above description was in terms of a single region of interest for the data. This could be generalised to several bins of data, with the likelihood function using the number of events n_i in each bin, and the expected signal and background rates s_i and b_i. A special case is the use of just two bins, where one (the 'control region') is designed to contain no signal and hence can be used to estimate the background; and the other contains background and possibly signal too. This is sometimes referred to as the 'on-off' problem. A parameter determination method can be used to find a range for σ.

To check for discovery ($\sigma \ne 0$) when there is one parameter of interest σ, we need to use some sort of 2-sided interval. Upper limits are unsatisfactory in that they always include zero, and lower limits because they almost always exclude it. Maximum probability density intervals are not ideal because they are metric-dependent. In problems with one or more dimensions, the Feldman–Cousins approach [1] in principle provides a sensible way of defining the acceptance region; however, computing time may be an issue in several dimensions. By changing the confidence level for the interval until $\sigma =0$ is just excluded, we can determine the significance of the observed effect (see section 2.7).

Analyses often do involve more than one physics parameter of interest. For example, for a new particle its mass M and cross-section σ could be the parameters. Similarly for a 2-neutrino oscillation analysis, the parameters are ${\sin }^{2}2\theta$ and ${\rm{\Delta }}{m}^{2}$.¹¹ The procedure in such cases could use a single 2-dimensional approach; or alternatively a 'Raster Scan', where statements are made about one parameter for a variety of fixed values of the other (e.g. conclusions about σ for fixed values of M). Which is better depends on the circumstances [12]. For parameter determination, the 2-dimensional approach is preferable; the Raster Scan would provide a range for σ at each M, regardless of whether the fit to the data for those values of M and σ was acceptable or not. In contrast, the Raster Scan is recommended for hypothesis testing. It can find more than one mass peak; and it will not result in exclusion at other masses M₂ just because there is a peak at M₁. This would be especially undesirable if the analysis had little or no sensitivity to the presence of a peak at M₂. Finally Wilks' theorem (see section 2.10.1) does not apply when ($M,\sigma$) or (${{\rm{\sin }}}^{2}2\theta ,{\rm{\Delta }}{m}^{2}$) are the variables, whereas it may if a series of fixed values of M or ${\rm{\Delta }}{m}^{2}$ respectively are tested as in a Raster Scan. Thus in the Higgs search culminating in its discovery at 125 GeV (see section 4.6.1), and in analyses resulting in the exclusion of other Higgs-like bosons at higher masses, Raster Scans were performed.

In comparing two hypotheses, the weighted sums of squared deviations¹² S₀ and S₁ may be calculated by comparing the data with predictions of each hypothesis in turn (see equation (2.7)). In deciding whether or not it is possible to choose between the two hypotheses, it is generally better to use ${\rm{\Delta }}S={S}_{0}-{S}_{1}$, rather than the individual S₀ and S₁, especially if the numbers of degrees of freedom involved in calculating S₀ and S₁ are large.

Many methods of distinguishing between hypotheses involve their p-values, which are the probabilities of obtaining a test statistic as deviant from the predicted value or more so, i.e. they are the tail probabilities of the relevant pdf of the test statistic (see figure 3). These are discussed in more detail in section 2.6.

Zoom In Zoom Out Reset image size

It is important to realise that there are several different types of likelihood ratio that are used in particle physics analyses:

• Comparison with Saturated Model [13]: The maximum value of an unbinned likelihood is not a measure of goodness of fit. However for a histogram, −2 times the Poisson $\mathrm{ln}({\rm{LR}})$
  $$\begin{eqnarray}-2\ast \mathrm{ln}{{\rm{LR}}}_{{\rm{Poisson}}} & = & -2\ast {\rm{\Sigma }}\mathrm{ln}[P({n}_{i}| {\kappa }_{i})/P({n}_{i}| {n}_{i})]\\ & = & -2\ast {\rm{\Sigma }}[{n}_{i}\mathrm{ln}({\kappa }_{i}/{n}_{i})+{n}_{i}-{\kappa }_{i}]\end{eqnarray} \tag{ 2.6 }$$
  asymptotically approximates to the weighted sum of squares (see equation (2.7) below). Here, $P({n}_{i}| {\kappa }_{i})$ is the Poisson probability for observing n_i events in the ith bin, when the predicted value is ${\kappa }_{i};$ $P({n}_{i}| {n}_{i})$ is the Poisson probability, evaluated for the best value of ${\kappa }_{i}$ i.e. n_i; and the summation is over the bins of the histogram. This provides a likelihood method of estimating goodness of fit for histograms. For small values of the ${\kappa }_{i}$, it is better than the standard weighted sum of squares
  $$\begin{eqnarray}&&S={\rm{\Sigma }}{({\kappa }_{i}-{n}_{i})}^{2}/{n}_{i}\ \ \ {\rm{or}}\ \ \ \ {\rm{\Sigma }}{({\kappa }_{i}-{n}_{i})}^{2}/{\kappa }_{i}.\end{eqnarray} \tag{ 2.7 }$$
  These are respectively the Neyman or Pearson versions. They rely on the Gaussian approximation for the Poisson.
• The Feldman–Cousins method [1] also makes use of the likelihood ratio ${ \mathcal L }(\kappa | x)/{ \mathcal L }({\kappa }_{{\rm{best}}}| x)$ as its ordering rule for the Neyman construction of the confidence belt for the likely values of the observable x for each value of the parameter κ. Here ${\kappa }_{{\rm{best}}}$ is the best physical value of the parameter κ for the particular x.
• ${{ \mathcal L }}_{1}/{{ \mathcal L }}_{0}$: For choosing between two hypotheses H0 and H1, the ratio of their corresponding likelihoods $\lambda ={{ \mathcal L }}_{1}/{{ \mathcal L }}_{0}$ can be used as a test statistic. For example, in the search for the SM Higgs production at the Fermilab Tevatron, this was used by setting $\mu =0$ for H0 and $\mu =1$ for H1. Here μ is the ratio between the observed cross-section and that predicted by the SM. Analyses trying to distinguish between two different orderings of neutrino masses (normal and inverted hierarchies) also use this form of LR [14, 15].In the absence of nuisance parameters, this LR is that used in implementing the Neyman–Pearson lemma for the optimal way of choosing between two simple hypotheses.When the expected distributions of the test statistic for the H0 and H1 are Gaussians of equal width, $-2\ast \mathrm{ln}({{ \mathcal L }}_{1}/{{ \mathcal L }}_{0})={S}_{1}-{S}_{0}$. There is an offset between them, however, when the Gaussian widths are different.
• In contrast, for the LHC Higgs search, two LR test statistics were used: ${ \mathcal L }(\mu =0)/{ \mathcal L }({\mu }_{{\rm{best}}})$ for testing discovery, and ${ \mathcal L }(\mu =1)/{ \mathcal L }({\mu }_{{\rm{best}}})$ for exclusion of H1. i.e. each hypothesis was tested separately, using its likelihood compared with the best possible one (as μ was varied). There are also variants of these test statistics where the LR is set equal to unity when the best value of μ is in certain ranges. These are to avoid excluding H₀ when there is a downward fluctuation causing ${\mu }_{{\rm{best}}}$ to be negative; or an upward fluctuation that might result in the exclusion of H1 [16]. An extension of these test statistics defined with a continuous value of μ (i.e ${ \mathcal L }(\mu )/{ \mathcal L }({\mu }_{{\rm{best}}})$) is used to determine an interval in the parameter μ. With ${\mu }_{{\rm{best}}}$ restricted to be greater than or equal to zero, the Feldman–Cousins method is reproduced.An advantage of these test statistics (as compared with a generic ratio of likelihoods) is that their asymptotic distributions are easier to calculate, thus reducing the need for Monte Carlo simulation. The expected distributions of various log likelihood ratios are discussed in section 2.10.

As well as depending on the parameter(s) of interest ϕ, in the presence of systematic effects, each of these likelihoods also depends on the corresponding nuisance parameters θ. Within the likelihood paradigm, the profile likelihood is often used. This is obtained for each value of the physics parameter(s) by re-maximising the likelihood by suitable choice of θ, i.e. ${{ \mathcal L }}_{{\rm{profile}}}(\phi )={ \mathcal L }(\phi ,{\theta }_{{\rm{best}}}(\phi ))$. For the LHC version of the test statistic, an advantage is that asymptotically it is independent of any nuisance parameters θ, and so it is not necessary to consider the effect of varying the θ. In contrast the Tevatron approach deals with the systematics by integrating over them.

For a comparison of p-values and LR for hypothesis testing, see section 2.6.1. The method used in Particle Physics of quoting p-values with the test statistic being a LR can be regarded as either:

• A p-value method where the test statistic just happens to be a LR; or
• A LR method where p-values are used merely to calibrate the significance of the values of the LR.

There are various Bayesian methods available (e.g. posterior probability ratio, Bayes factor, AIC [17], BIC [18], etc). Largely because of their dependence on the choice of priors occurring in one hypothesis and not in the other, these methods tend not to be used in Particle Physics searches for new phenomena.

As already mentioned, this is not formally used much in searches for New Physics. It is, however, one of the reasons behind the different cut-offs on the p-values used for claiming a discovery ($5\sigma$) and for excluding some form of New Physics (5%). Here the relevant 'costs' in the discovery case include the embarrassment of making a discovery claim which turns out to be incorrect, as compared with the reward for making the discovery when it is true.

There is sometimes discussion of why an approach using the ${\rm{LR}}={{ \mathcal L }}_{1}/{{ \mathcal L }}_{0}$ can give a very different numerical answer from a p-value calculation. A reason some agreement might be expected is that they are both addressing the question of whether there is evidence in the data for new physics.

In fact they measure very different things. Thus p₀ (see figure 3 (b)) simply measures the consistency with the null hypothesis, while the LR takes both hypotheses into account. Furthermore a p-value is the tail area of a pdf beyond the observed value of the test statistic, while a likelihood uses the value of the pdf at the observed test statistic. There is thus no reason to expect them to bear any particular relationship to each other. Figure 4(b) shows that even at a given value of p₀, the LR can have any value within a range.

Zoom In Zoom Out Reset image size

Sometimes the effect of New Physics could appear in several different reactions in a given experiment; or in different experiments. It would then be sensible to combine the information, in order to improve the sensitivity of the search. The best way to do this is to perform a joint analysis, but this is not always possible, so an alternative is to combine the individual p-values. Problems with this are:

• Different effects: Two analyses might each have small p₀-values because they both disagree with the SM, but have different discrepancies e.g. peaks at quite different mass values.
• Non-uniqueness: As reviewed by Cousins [21], in the combination of $n\ p$-values, we are trying to find a transformation from the n (presumed uniform) 1-dimensional distributions to just one uniform distribution; this can clearly be achieved in many ways, thus yielding a variety of different possible combined p-values. Which is best requires extra information about the possible alternative hypotheses we think might be relevant; and more details of the individual analyses beyond just their p-values. Again the desirability of a combined analysis is demonstrated.
• Selection bias: Care must be taken to combine all relevant analyses, and not just the ones which give small individual p-values.

One possibility is to calculate the product ${\pi }_{{\rm{obs}}}$ of the n individual p-values. Then the probability P of the product of n independent uniformly distributed p-values being smaller than ${\pi }_{{\rm{obs}}}$ is

$$\begin{eqnarray}&&P={\pi }_{{\rm{obs}}}{\rm{\Sigma }}{(-{\rm{ln}}{\pi }_{{\rm{obs}}})}^{j}/j!\end{eqnarray} \tag{ 2.8 }$$

where the summation extends from j  =  0 to $n-1$. Thus P is larger than ${\pi }_{{\rm{obs}}}$. For two measurements $P={p}_{1}{p}_{2}(1-{\rm{ln}}({p}_{1}{p}_{2}))$.

If the p_i-values were obtained from weighted sums of squares S_i which are expected to have ${\chi }^{2}$ distributions with numbers of degrees of freedom ${\nu }_{i}$, an alternative is to use ${S}_{{\rm{comb}}}={\rm{\Sigma }}{S}_{i}$ and ${\nu }_{{\rm{comb}}}={\rm{\Sigma }}{\nu }_{i}$ to obtain the overall P. In the special case where the individual ${\nu }_{i}$ are all 2, this becomes equivalent to using equation (2.8).

Another, known as the Stouffer method [22], uses

$$\begin{eqnarray}&&{z}_{{\rm{comb}}}={\rm{\Sigma }}{z}_{i}/\sqrt{n},\end{eqnarray} \tag{ 2.9 }$$

where z_i are the significances¹⁴ or signed z-scores (i.e the number of standard deviations corresponding to the one-sided Gaussian p-value, with p = 0.5 corresponding to z = 0) and z_comb is the combined value.

When comparing data with two hypotheses, a common test statistic t is −2 times the logarithm of the likelihood ratio, or almost equivalently the difference in S_min between the two hypotheses, where S_min is the usual weighted sum of squared discrepancies between theory and data, minimised with respect to any free parameters in each case. Wilks' theorem [28] states that under special circumstances, the distribution of t is the mathematical ${\chi }^{2}$, with the number of degrees of freedom ν equal to the difference in the numbers of free parameters in the two hypotheses.

The conditions are:

• The null hypothesis must be true.
• The hypotheses are nested.
• When we reduce H1 to H0, the values of all the extra parameters must be well defined, and not at an extreme end of their possible ranges.
• There are enough data for asymptotic approximations to apply.

When H0 is true, we still expect its S_min to be larger than that for H1 (because H1 has extra parameters), but the difference should be small. However, when H1 is true, the difference can be large. Wilks' theorem provides an easy way of assessing how small is small: the value of t can immediately be turned into a p-value for the null hypothesis, via the known distribution of t. When the theorem does not apply, t may well still be a useful test statistic, but it will usually be necessary to calculate its expected distribution under H0 by simulation.

The examples below illustrate when the conditions do or do not apply:

• Polynomials: e.g. H0 is a straight line, and H1 is a polynomial of degree 3. These are nested since by setting the squared and cubic terms in H1 to zero, we obtain H0. Furthermore, these special values are not at the end of mathematically or physically allowed ranges. Wilks' theorem applies, with ν equal to the difference in the number of free parameters in the two hypotheses; in this example $\nu =2$.
• Presence of signal. Here H0 is just a background distribution, while H1 also has a signal bump of mass M and cross-section σ. H1 becomes H0 if σ is set to zero, but in that case M is undefined. Also unless σ is allowed to be negative, its special value is at the end of its physical range. Thus the theorem does not apply. However, if M is kept fixed (and negative values for σ are allowed), either because its value is known or because we are performing a Raster Scan, Wilks' theorem may apply.
• Spin of resonance: We try to determine whether the spin of a resonance is 1 or 2. These hypotheses are not nested, so the theorem does not apply.

We can think of each of these categories as having a likelihood function ${{ \mathcal L }}_{c}(\mu ,\theta );$ with μ defined as the ratio of the assumed value to the theoretical cross-section as in section 2.5.6. The statistical combination of the data in all of the independent categories, is performed by taking the product,

$$\begin{eqnarray}&&{ \mathcal L }(\mu ,\theta )=\displaystyle \prod _{c}{{ \mathcal L }}_{c}(\mu ,\theta ),\end{eqnarray} \tag{ 4.1 }$$

where θ are the nuisance parameters which may be correlated between different categories. An important point to notice is that the parameter μ is continuous. One could imagine only being interested in the value $\mu =1$ since that would yield the signal rates, in every category, as predicted by the SM. However, if nature contained some new physics which altered that rate in some way, or perhaps the SM calculations were incorrect in some fashion, we would not want to miss discovering something, just because it was not exactly the Higgs boson predicted by the SM. In fact, one can think of the μ parameter as representing a whole family of alternate hypotheses $H(\mu )$ which predict the existence of a 'Higgs-like' boson, with similar properties, decay rates etc, but being produced with a different rate. Instead then of comparing $\mu =0$ with $\mu =1$, we compare it with $\mu =\hat{\mu }$ which denotes the value of μ which maximises the likelihood ${ \mathcal L }(\mu ,\theta )$ for the data we observed, thereby giving the 'most-likely' value for μ. The value of ${m}_{{{\rm{H}}}^{0}}$ is implicitly included in the definition of $H(\mu )$ so typically, a comparison between $H(\mu )$ and H0 is made for different values of ${m}_{{{\rm{H}}}^{0}}$, within some reasonable range in which the search is expected to be sensitive.

The quantity ${q}_{0}=-2\mathrm{ln}({ \mathcal L }(\mu =0,{\theta }_{{\rm{p}}{\rm{r}}{\rm{o}}{\rm{f}}}^{0})/{ \mathcal L }({\mu }_{{\rm{best}}},{\theta }_{{\rm{best}}}))$, was used as the test statistic for the discovery of the Higgs boson. Here, ${\theta }_{{\rm{p}}{\rm{r}}{\rm{o}}{\rm{f}}}^{0}$ denotes the values of the nuisance parameters resulting from profiling the likelihood function when $\mu =0$; similarly ${\theta }_{{\rm{best}}}$ denotes the corresponding values for $\mu ={\mu }_{{\rm{best}}}$. Clearly, large values of q₀ would indicate that the data would favour strongly the presence of a Higgs boson signal, which we would call a 'significant' result. By comparing the observed q₀ with its expected distribution (usually determined by simulation, but sometimes available analytically), we can determine the p-value p₀; this is the probability of a statistical fluctuation, consistent with a particular value of ${m}_{{{\rm{H}}}^{0}}$, resulting in q₀ being at least as extreme as its observed value, assuming that H0 is true. An example distribution of q₀ under H0 is given in figure 9. The observed value of q₀ is indicated by the arrow.

Zoom In Zoom Out Reset image size

The probability to observe a value of q₀, under H0, which is larger than the one we observed in the real data is called the local p-value p₀, for the assumed value of ${m}_{{{\rm{H}}}^{0}}$. It can be obtained from the fraction of these pseudo-data sets for which ${q}_{0}\gt {q}_{0}^{{\rm{obs}}}$,

$$\begin{eqnarray}&&{p}_{0}={\rm{Prob}}({q}_{0}\gt {q}_{0}^{{\rm{obs}}}| H0)={\int }_{{q}_{0}^{{\rm{obs}}}}^{+\infty }f({q}_{0}| H0){{\rm{d}}{q}}_{0},\end{eqnarray} \tag{ 4.2 }$$

where $f({q}_{0}| H0)$ is the pdf for q₀ under H0, for the particular mass being investigated. As the significance of an observation increases, it becomes more and more computationally challenging to produce these distributions with pseudo-data. This issue is overcome through the use of asymptotic approximations which are valid when there are a large number of events. In figure 9, the expected distribution of q₀ using these approximations, shown as the green curve, agrees very well with that from the toys and hence the local p-value can be determined analytically. ATLAS and CMS use these asymptotic approximations when calculating p-values for different values of ${m}_{{{\rm{H}}}^{0}}$. The local p-value for the data collected by CMS is shown in figure 10 using data from each decay mode separately and the combination of all modes. The smallest such value occurs at the large dip with a p-value of ≈3 × 10⁻⁷, corresponding to a z-score of about 5 standard deviations, the predefined requirement to announce a discovery. The dashed line in the plot shows what the median expected value of p₀ would be at each value of ${m}_{{{\rm{H}}}^{0}}$ if the Higgs boson was ${m}_{{{\rm{H}}}^{0}}$, according to the SM predictions.

Zoom In Zoom Out Reset image size

As explained in section 2.8, we need to convert the local p-value into a global one; the former gives the chance of a statistical fluctuation giving an effect at 125 GeV at least as significant as the observed one, while the global one is for a range of possible Higgs masses. This is chosen to be a reasonable range over which all of the channels used in the combination are sensitive. Above this range, for example, the branching fractions for the Higgs decays to $\gamma \gamma ,\,{bb}$ and $\tau \tau$ channels become negligible compared to that for WW and ZZ. For some channels, like the high resolution channel $H\to \gamma \gamma$, there are several wiggles in the local p-value. This indicates the fact that the peak could occur in a number of different places (or specifically, different values of ${m}_{{{\rm{H}}}^{0}}$) so that a large look elsewhere effect is expected in this channel. Conversely, for channels like the $H\to {\mathrm{WW}}^{* }\to 2l2\nu$, with poorer resolution, the p-value is rather smooth with respect to ${m}_{{{\rm{H}}}^{0}}$ meaning if an excess is seen at one value of ${m}_{{{\rm{H}}}^{0}}$, the same excess will be compatible with other, nearby values of ${m}_{{{\rm{H}}}^{0}}$ too. Effectively, this suggests that there are fewer independent values of ${m}_{{{\rm{H}}}^{0}}$ which can be tested. For the combination, for which outcomes are more significant if the excess appears around the same mass in the two high-resolution channels, the look elsewhere effect is somewhere in between. The global p-value (using values of ${m}_{{{\rm{H}}}^{0}}$ in the range shown in figure 10) is $3.4\times {10}^{-6}$, corresponding to a z-score of 4.5.

Many analyses at ATLAS and CMS rely on theoretical calculations and Monte Carlo simulation in order to calculate the number of background events, b in equation (2.5), from a particular process. This could be determined using

$$\begin{eqnarray}&&b={A}_{b}{\epsilon }_{b}L{\sigma }_{b},\end{eqnarray} \tag{ 4.3 }$$

where ${\sigma }_{b}$ is the production cross-section for the background process, A_b and ${\epsilon }_{b}$ are the overall acceptance and efficiency for background events to pass all of the selections and L is the integrated luminosity. Typically, quantities such as the efficiency, acceptance and luminosity can be measured with an accuracy of a few percent at the LHC. The predicted cross-section can, in some cases, only be known to within 20% to 30%! In ATLAS and CMS, these kind of uncertainties are modelled assuming a log-normal distribution meaning the number of background events becomes

$$\begin{eqnarray}&&b\to b(\theta )=b\cdot {(1.2)}^{\theta },\end{eqnarray} \tag{ 4.4 }$$

where θ is the nuisance parameter for the cross-section uncertainty²⁵ (of size 20%). The nuisance parameter is assigned a Gaussian prior probability (or in the frequentist case the likelihood is multiplied by a Gaussian term) centred at 0 with a width of 1. This means that a value of $\theta =1$ corresponds to a positive shift of the nuisance parameter by one standard deviation and yields a background rate of $1.2b$. Instead, a value of $\theta =-1$ is a negative shift of one standard deviation and yields a value of $b/1.2\approx 0.83b$. This seems a little odd in that the value of $b(\theta )$ is not symmetric about b. However, this is an intentional feature of using log-normals since the number of expected background events can never go below 0.

Suppose we have an experiment where b  =  900, while the number of signal expected is 90. If we observe 990 events in our data, with no systematic uncertainty on the background, this would constitute a $3\sigma$ evidence for our signal. However, when considering a 20% uncertainty, this excess of events could be explained by assuming $\theta \approx 0.52$, only around half a standard deviation from its assumed value. The significance of the excess would drop to roughly $0.5\sigma$ which is not worth getting excited over.

If the rate of expected signal is very small compared to the backgrounds as in this case, uncertainties of this size will make it very difficult to determine whether or not any excess seen in the data is some hint of new physics or just some inaccuracy in the background rate, which can easily be explained by shifting the nuisance parameters a little. Worse yet, adding more data will not solve the problem since typically systematic uncertainties are not statistical in nature so remain a constant frustration as the luminosity increases. In the Higgs boson search at ATLAS and CMS, the majority of the work is understanding these background rates and reducing their systematic uncertainties. Moreover, in real analyses, there are typically several contributing backgrounds, which must be summed to obtain the total, potentially having different sources of systematic uncertainties.

One common way to reduce these kinds of systematic uncertainties is to directly use the data to estimate the background rates.

For searches which use kinematic variables to separate the signal from the background, both the rate and the distribution (shape) of the background must be accurately modelled. This often means accounting for additional experimental uncertainties related to the reconstruction of the decay products, such as the momenta of the reconstructed particles in the event. Typically, searches are performed in the tails of the background distributions which make their estimation from MC unreliable.

In the ATLAS and CMS ${{\rm{H}}}^{0}\to \gamma \gamma$ analyses, the backgrounds are determined by fitting functions to the invariant mass spectra directly in the data. This means that there are no uncertainties associated with the luminosity, cross-section or acceptance and efficiency of the backgrounds. However, there is no such thing as a (systematics) free lunch in particle physics. The problem still arises that one has to select the appropriate function to fit the background. Without relying heavily on theoretical calculations, a number of different smooth functions will provide good fits to the background in this analysis. This choice itself is a systematic uncertainty which can be handled in a number of different ways.

At the time of discovery, the CMS analysis used polynomial functions to model the background whose number of parameters were chosen such that any potential bias from choosing such a function would be much smaller than the statistical uncertainty, thereby rendering it negligible. Since then, the analysis has adopted a method in which this choice of function is a nuisance parameter which can take on discrete values, each of which represents a particular choice from a list of suitable functional forms [44]. This discrete nuisance parameter is profiled, along with all the other nuisance parameters, effectively leaving the choice of function itself free in the fit²⁶ .

Accurately determining both the rate and the shape of the background contributions was not only an issue during the search for the Higgs boson but continues to present a challenge for measuring its properties. As the precision of these measurements increases, due to inclusion of new data, the issue of reducing systematic uncertainties becomes ever more difficult and the techniques used to overcome them become increasingly sophisticated.

An important property of the Higgs boson in the SM is that it should have zero intrinsic spin. If the spin of the discovered particle were something else, it could not be identified with the SM Higgs boson. Since the spin of a particle is not a continuous quantity, a hypothesis test is performed to determine whether or not the hypothesis of a spin-0 particle can be rejected in favour of an alternative hypothesis such as spin-2. Experimentally, these two hypotheses would yield different angular distributions of the Higgs decay products. For example, in the ${{\rm{H}}}^{0}\to {\mathrm{ZZ}}^{* }\to 4l$ channel, the angles defined using the directions of the leptons will depend on whether or not the decaying boson has spin-0 or spin-2. Using the distribution of these angles in the data the two hypotheses, spin-0 or spin-2, can be distinguished. A test statistic is defined using the ratio of the likelihoods to observe the distribution in the data under each of these two hypotheses. As usual, the convention is to take the log of this ratio so that the numbers are of the order 1, and multiply by a factor of −2 so that the quantity asymptotically approximates to the difference of the weighted sum of squares for the two hypotheses (see section 2.5.6). Figure 11 shows the distribution of the test statistic under each of the two hypotheses, generated using pseudo-experiments. The shaded blue or solid yellow histograms show the distribution that we would expect if the decay products originated from a spin-2 boson or the SM Higgs boson (with spin-0), respectively. The red arrow indicates the value of the test statistic evaluated on the actual angular distribution which is observed in the data collected by CMS [45].

Zoom In Zoom Out Reset image size

From these distributions two p-values are calculated. The first (${p}_{{2}^{+}}$) is the probability to observe a value which is at least as large as the one observed, assuming the spin-2 hypothesis. The value of ${p}_{{2}^{+}}$ in this case is around 0.008, which corresponds to around 3 standard deviations from the mean of a normal distribution. The other p-value (${p}_{{0}^{+}}$) is defined as the probability to observe a value which is smaller than the one observed, assuming the spin-0 hypothesis. In this case, ${p}_{{0}^{+}}$ is around 0.62. From this we can conclude that the hypothesis of a spin-0 particle is certainly more favourable than a spin-2. Note however, that in this kind of test, had the observed data yielded a value of the test statistic around -20 or smaller, neither the spin-2, nor the spin-0 hypothesis would be a very good candidate and both could have been rejected. In such a case, one would have to revisit the list of possible models to test and it would have probably sparked a lot of interest in the particle physics community. When performing hypothesis tests of this kind, ATLAS and CMS use the CL_s criterion (defined as ${p}_{{2}^{+}}/(1-{p}_{{0}^{+}})$) to avoid excluding in cases where the two distributions lie almost on top of one another and the observed value of the test statistic lies somewhere in the tails of the distributions. In the example shown in figure 11, ${{CL}}_{s}=0.02$ and so with the usual 95% confidence level criterion, this particular spin-2 model is excluded.

A major task was to precisely measure the mass of the Higgs boson since, in the SM, once this is known, all of the other properties of the Higgs boson can be predicted. In order to get the best measurement of the mass, ATLAS and CMS used their datasets from the two high resolution decay channels, $H\to \gamma \gamma$ and $H\to {\mathrm{ZZ}}^{* }\to 4l$. The likelihoods from ATLAS and CMS were combined to give an overall likelihood for the combined dataset [46]. Figure 12 shows the result of that combination. The figure shows $-2{\rm{ln}}{\rm{\Lambda }}({m}_{{{\rm{H}}}^{0}})$ (where ${\rm{\Lambda }}({m}_{{{\rm{H}}}^{0}})$ is the ratio of the profile likelihood at a particular value of ${m}_{{{\rm{H}}}^{0}}$ to the value obtained at ${m}_{{{\rm{H}}}^{0}}={m}_{{{\rm{H}}}^{0},{\rm{best}}}$) for the sub-combinations of each decay mode separately and the overall combination. In each case the minimum of the curve indicates the value of ${m}_{{{\rm{H}}}^{0}}$ which best fits that particular dataset. The widths of the curves indicate the level of uncertainty in each of the measurements. By appealing to Wilks' theorem, the uncertainty of the measurement is estimated from the values at which the curve increases by 1 from the minimum i.e the regions in ${m}_{{{\rm{H}}}^{0}}$ for which $-2{\rm{ln}}{\rm{\Lambda }}({m}_{{{\rm{H}}}^{0}})\leqslant 1$. Naturally, the width of the combined measurement is smaller than in either of the two channels alone. This is expected since the data used in each of the two channels are statistically independent from one another. Since the statistical component of the uncertainties is much larger than the systematic component, as shown by the small differences between the dashed and solid lines, the systematics play a small role in the combination.

Zoom In Zoom Out Reset image size

Many other properties of the Higgs boson are currently under study at ATLAS and CMS. As more and more data are collected, the precision to which these properties can be determined and the range in the types of properties which can be measured will increase. Up to this point however, the measured properties of the discovered particle are consistent with those that are predicted for the Higgs boson in the SM. However, particle physicists are ever hopeful that at some point we will see something we did not expect and that is when the real challenge and excitement will begin.