# 1 Introduction

## 1 Introduction

A consistently recurring topic at LEP2 has been the interpretation and combination of results from searches for new particles. The fundamental task is to interpret the collected dataset in the context of two complementary hypotheses. The first hypothesis – the null hypothesis – is that the dataset is compatible with non-signal Standard Model background production alone, and the second is that the dataset is compatible with the sum of signal and Standard Model background production. In most cases, the search for new particles proceeds via several parallel searches for final states. The results from all of these subchannels are then combined to produce a final result.

All existing confidence level calculations follow the same general strategy . A test statistic or estimator is constructed to quantify the “signal-ness” of a real or simulated experiment. The “signal-ness” of a single observed experiment leads to the confidence level on, for example, the null hypothesis that the observed experiment is incompatible with signal and background both being produced. Most calculation methods use an ensemble of toy Monte Carlo experiments to generate the estimator distribution against which the observed experiment is compared. This generation can be rather time-consuming when the number of toy Monte Carlo experiments is great (as it must be for high precision calculations) or if the number of signal and background expected for each experiment is great (as it is for the case of searches optimized to use background subtraction).

In this note, we present an improved method for calculating confidence levels in the context of searches for new particles. Specifically, when the likelihood ratio is used as an estimator, the experiment estimator distribution may be calculated analytically with the Fourier transform. With this approach, the disadvantage of toy Monte Carlo experiments is avoided. The analytic method offers several advantages over existing methods, the most dramatic of which is the increase in calculation speed and precision.

## 2 Likelihood ratio estimator for searches

The likelihood ratio estimator is the ratio of the probabilities of observing an event under two search hypotheses. The estimator for a single experiment is

$$E=C\frac{_{s+b}}{_b}.$$
(1)   
Here $`_{s+b}`$ is the probability density function for signal+background experiments and $`_b`$ is the probability density function for background-only experiments. Because the constant factor $`C`$ appears in each event’s estimator, it does not affect the ordering of the estimators – an event cannot become more signal-like by choosing a different $`C`$. For clarity in this note, the constant is chosen to be $`e^s`$, where $`s`$ is the expected number of signal events. <sup>1</sup><sup>1</sup>1When considering the two production hypotheses and calculating an exclusion, the expected signal $`s`$ is uniquely determined by the cross section. If the cross section is not fixed, then $`e^s`$ is not constant, and $`C`$ may be set to unity. For the simplest case of event counting with no discriminant variables (or, equivalently, with perfectly non-discriminating variables), the estimator can be calculated with Poisson probabilities alone. In practice, not every event is equally signal-like. Each search may have one or more event variables that discriminate between signal-like and background-like events. For the general case, the probabilities $`_{s+b}`$ and $`_b`$ are functions of the observed events’ measured variables.

As an example, consider a search using one discriminant variable $`m`$, the reconstructed Higgs mass. The signal and background have different probability density functions of $`m`$, defined as $`f_s(m)`$ and $`f_b(m)`$, respectively. (For searches with more than one discriminant variable, $`m`$ is replaced by a vector of discriminant variables $`\stackrel{}{x}`$.) It is then straightforward to calculate $`_{s+b}`$ and $`_b`$ for a single event, taking into account the event weighting coming from the discriminant variables:

$$E=e^s\frac{P_{s+b}}{P_b}=e^s\frac{e^{(s+b)}\left[sf_s(m)+bf_b(m)\right]}{e^b\left[bf_b(m)\right]}.$$
(2)   
The likelihood ratio estimator can be shown to maximize the discovery potential and exclusion potential of a search for new particles . Such an estimator, both with and without discriminant variables, has been used successfully by the LEP2 collaborations to calculate confidence levels for searches .

## 3 Ensemble estimator distributions via Fast Fourier Transform (FFT)

One way to form an estimator for an ensemble of events is to generate a large number of toy Monte Carlo experiments, each experiment having a number of events generated from a Poisson distribution. Another way is to analytically compute the probability density function of the ensemble estimator given the probability density function of the event estimator. The discussion of this section pursues the latter approach.

The likelihood ratio estimator is a multiplicative estimator. This means the estimator for an ensemble of events is formed by multiplying the individual event estimators. Alternatively, the logarithms of the estimators may be summed. In the following derivation, $`F=\mathrm{ln}E`$, where $`E`$ is the likelihood ratio estimator.

For an experiment with 0 events observed, the estimator is trivial:

$`E`$ $`=`$ $`e^s{\displaystyle \frac{e^{(s+b)}}{e^b}}=1`$ (3)   
$`F`$ $`=`$ $`0`$ (4)   
$`\rho _0(F)`$ $`=`$ $`\delta (F),`$ (5)   
where $`\rho _0(F)`$ is the probability density function of $`F`$ for experiments with 0 observed events.

For an experiment with exactly one event, the estimator is, again using the reconstructed Higgs mass $`m`$,

$`E`$ $`=`$ $`e^s{\displaystyle \frac{e^{(s+b)}\left[sf_s(m)+bf_b(m)\right]}{e^b\left[bf_b(m)\right]}},`$ (6)   
$`F`$ $`=`$ $`\mathrm{ln}{\displaystyle \frac{sf_s(m)+bf_b(m)}{bf_b(m)}},`$ (7)   
and the probability density function of $`F`$ is defined as $`\rho _1(F)`$.

For an experiment with exactly two events, the estimators of the two events are multiplied to form an event estimator. If the reconstructed Higgs masses of the two events are $`m_1`$ and $`m_2`$, then

$`E`$ $`=`$ $`{\displaystyle \frac{\left[sf_s(m_1)+bf_b(m_1)\right]\left[sf_s(m_2)+bf_b(m_2)\right]}{\left[bf_b(m_1)\right]\left[bf_b(m_2)\right]}}`$ (8)   
$`F`$ $`=`$ $`\mathrm{ln}{\displaystyle \frac{sf_s(m_1)+bf_b(m_1)}{bf_b(m_1)}}+\mathrm{ln}{\displaystyle \frac{sf_s(m_2)+bf_b(m_2)}{bf_b(m_2)}}.`$ (9)   
The probability density function for exactly two particles $`\rho _2(F)`$ is simply the convolution of $`\rho _1(F)`$ with itself:

$`\rho _2(F)`$ $`=`$ $`{\displaystyle \rho _1(F_1)\rho _1(F_2)\delta (FF_1F_2)𝑑F_1𝑑F_2}`$ (10)   
$`=`$ $`\rho _1(F)\rho _1(F).`$ (11)   
The generalization to the case of $`n`$ events is straightforward and encouraging:

$`E`$ $`=`$ $`{\displaystyle \underset{i=1}{\overset{n}{}}}{\displaystyle \frac{sf_s(m_i)+bf_b(m_i)}{bf_b(m_i)}}`$ (12)   
$`F`$ $`=`$ $`{\displaystyle \underset{i=1}{\overset{n}{}}}\mathrm{ln}{\displaystyle \frac{sf_s(m_i)+bf_b(m_i)}{bf_b(m_i)}}`$ (13)   
$`\rho _n(F)`$ $`=`$ $`{\displaystyle \int \cdots \int \underset{i=1}{\overset{n}{}}\left[\rho _1(F_i)dF_i\right]\delta \left(F\underset{i=1}{\overset{n}{}}F_i\right)}`$ (14)   
$`=`$ $`\underset{n\text{times}}{\underset{}{\rho _1(F)\mathrm{}\rho _1(F)}}.`$ (15)   
Next, the convolution of $`\rho _1(F)`$ is rendered manageable by an application of the relationship between the convolution and the Fourier transform.

If $`A(F)=B(F)C(F)`$, then the Fourier transforms of $`A`$, $`B`$, and $`C`$ satisfy

$$\overline{A(G)}=\overline{B(G)}\overline{C(G)}.$$
(16)   
This allows the convolution to be expressed as a simple power:

$$\overline{\rho _n(G)}=\left[\overline{\rho _1(G)}\right]^n.$$
(17)   
Note this equation holds even for $`n=0`$, since $`\overline{\rho _0(G)}=1`$. For any practical computation, the analytic Fourier transform can be approximated by a numerical Fast Fourier Transform (FFT) .

How does this help to determine $`\rho _{s+b}`$ and $`\rho _b`$? The probability density function for an experiment estimator with $`s`$ expected signal and $`b`$ expected background events is

$$\rho _{s+b}(F)=\underset{n=0}{\overset{\mathrm{}}{}}e^{(s+b)}\frac{(s+b)^n}{n!}\rho _n(F),$$
(18)   
where $`n`$ is the number of events observed in the experiment. Upon Fourier transformation, this becomes

$`\overline{\rho _{s+b}(G)}`$ $`=`$ $`{\displaystyle \underset{n=0}{\overset{\mathrm{}}{}}}e^{(s+b)}{\displaystyle \frac{(s+b)^n}{n!}}\overline{\rho _n(G)}`$ (19)   
$`=`$ $`{\displaystyle \underset{n=0}{\overset{\mathrm{}}{}}}e^{(s+b)}{\displaystyle \frac{(s+b)^n}{n!}}\left[\overline{\rho _1(G)}\right]^n`$ (20)   
$$\overline{)\overline{\rho _{s+b}(G)}=e^{(s+b)\left[\overline{\rho _1(G)}1\right]}}$$
(21)   
The function $`\rho _{s+b}(F)`$ may then be recovered by using the inverse transform. In general, this relation holds for any multiplicative estimator.

This final relation means that the probability density function for an arbitrary number of expected signal and background events can be calculated analytically once the probability density function of the estimator is known for a single event. This calculation is therefore just as fast for high background searches as for low background searches. In particular, it holds great promise for Higgs searches which, due to use of background subtraction and discriminant variables, are optimized to higher background levels than they have been in the past.

Two examples will provide practical proof of the principle. For the first, assume a hypothetical estimator results in a probability density function of simple Gaussian form

$$\rho _1(F)=\frac{1}{\sqrt{2\pi \sigma }}e^{\frac{(x\mu )^2}{2\sigma ^2}},$$
(22)   
where $`\sigma =0.2`$ and $`\mu =2.0`$. For an expected $`s+b=20.0`$, both the FFT method and the toy Monte Carlo method are used to evolute the event estimator probability density function to an experiment estimator probability density function. The agreement between the two methods (Fig. 1) is striking.

The higher precision of the FFT method is apparent, even when compared to 1 million toy Monte Carlo experiments. The periodic structure is due to the discontinuous Poisson distribution being convolved with a narrow event estimator probability function. In particular, the peak at $`\mathrm{ln}E=0`$ corresponds to the probability that exactly zero events be observed ($`e^{(s+b)}=2.1\times 10^9`$). The precision of the toy Monte Carlo method is limited by the number of Monte Carlo experiments, while the precision of the FFT method is limited only by computer precision.

For the second example, a more realistic estimator is calculated using a discriminant variable distribution from an imaginary $`\text{HZ}\text{H}\tau \tau `$ search. The variable used here is the reconstructed Higgs mass of the event. This estimator’s probability density function is then calculated for an experiment with $`s=5`$ and $`b=3`$ expected events (Fig. 2).

Again, the two methods agree well in regions where the toy Monte Carlo method is useful.

These examples support the mathematical proof of the FFT method described above. Because the final calculations $`c_{s+b}`$ and $`c_b`$ are simply integrals of the experiment estimator probability density function, any confidence levels calculated with the FFT method and the toy Monte Carlo method are identical. The examples also show the precision achievable with the FFT method, a precision that will be important when testing discovery hypotheses at the $`5\sigma =5\times 10^7`$ level.

## 4 Combining results from several searches

Given the multiplicative properties of the likelihood ratio estimator, the combination of several search channels proceeds intuitively. The estimator for any combination of events is simply the product of the individual event estimators. Consequently, construction of the estimator probability density function for the combination of channels parallels the construction of the estimator probability density function for the combination of events in a single channel. In particular, for a combination with $`N`$ search channels:

$`\overline{\rho _{s+b}(G)}`$ $`=`$ $`{\displaystyle \underset{j=1}{\overset{N}{}}}\overline{\rho _{s+b}^j(G)}`$ (23)   
$`=`$ $`e^{_{j=1}^N(s_j+b_j)\left[\overline{\rho _1^j(G)}1\right]}`$ (24)   
Due to the strictly multiplicative nature of the estimator, this combination method is internally consistent. No matter how subsets of the combinations are rearranged (i.e., combining channels in different orders, combining different subsets of data runs), the result of the combination does not change.

Once a results are obtained for $`\rho _{s+b}(F)`$ and $`\rho _b(F)`$, simple integration gives the confidence coefficients $`c_{s+b}`$ and $`c_b`$. From this point, confidence levels for the two search hypotheses may be calculated in a number of ways . Those straightforward calculations are outside the scope of this note.

## 5 Final remarks and conclusions

A few short remarks conclude this note and emphasize the advantages of calculations using the likelihood ratio with the Fast Fourier Transform (FFT) method.

1. The likelihood ratio estimator is an optimal ordering estimator for maximizing both discovery and exclusion potential. Such an estimator can only improve the discovery or exclusion potential of a search.
2. As a multiplicative estimator, the likelihood ratio estimator ensures internal consistency when results are combined. For example, if the dataset is split into several smaller pieces, the combined result always remains the same.
3. The probability density function of an ensemble estimator may be calculated analytically from the event estimator probability density function. Avoiding toy Monte Carlo generation brings revolutionary advances in speed and precision. For a $`\text{HZ}\text{4-jets}`$ search with 25 expected background events, a full confidence level calculation with $`2^{18}`$ toy MC experiments and 60 Higgs mass hypotheses takes approximately fifteen CPU hours. By contrast, the same calculation using the FFT method takes approximately two CPU minutes. This discrepancy only increases as the required confidence level precision and the number of toy MC experiments increase. For example, confidence level calculations for discovery at the $`5\sigma `$ level would require $`𝒪(10^8)`$ toy MC experiments. Given the approximately linear scaling of calculating time with number of toy experiments, such a calculation would take up almost a year in the 4-jet channel alone! The precision of the analytic FFT method is more than sufficient for a $`5\sigma `$ discovery.

A fast confidence level calculation makes possible studies that might have otherwise been too CPU-intensive with the toy MC method. These include studies of improvements in the event selections, of various working points, and of systematic errors and their effects, among others. A precise calculation makes possible rejection of null hypotheses at the level necessary for discovery.

The marriage of the likelihood ratio estimator and the FFT method seems well-suited for producing extremely fast and precise confidence level results, and the flexibility and ease of use of the clfft package should make this a powerful tool in interpreting searches for new particles.