Extending the Bump Hunt with Machine Learning

Introduction

Searching for new resonances as bumps in the invariant mass spectrum of the new particle decay products is one of the oldest and most robust techniques in particle physics, from the $\rho$ meson discovery PhysRev.126.1858 and earlier up through the recent Higgs boson discovery Aad:2012tfa; Chatrchyan:2012xdj. This technique is very powerful because sharp structures in invariant mass spectra are not common in background processes, which tend to produce smooth distributions. As a result, the background can be estimated directly from data by fitting a shape in a region away from the resonance (sideband) and then extrapolating to the signal region. It is often the case that the potential resonance mass is not known a priori and a technique like the BumpHunter Choudalakis:2011qn is used to scan the invariant mass distribution for a resonance. In some cases, the objects used to construct the invariant mass (e.g. jet substructure) and their surroundings (e.g. presence of additional forward jets) have properties that can be used to increase the signal purity. Both ATLAS and CMS The references here are the Run 2 results; Run 1 results can be found within the cited papers. The techniques described in this paper also apply to leptonic or photonic final states, but jets are used as a prototypical example due to their inherent complex structure. have conducted extensive searches for resonances decaying into jets originating from generic quarks and gluons Aaboud:2017yvp; Sirunyan:2016iap; Khachatryan:2015dcf, from boosted $W$ Sirunyan:2017acf; Aaboud:2017eta, $Z$ Sirunyan:2016wqt, $Z^{\prime}$ Sirunyan:2017dnz; Sirunyan:2017nvi; Aaboud:2018zba or Higgs bosons Sirunyan:2017dgc; Sirunyan:2017isc; Aaboud:2017ahz, from $b$-quarks Aaboud:2016nbq, as well as from boosted top quarks Sirunyan:2017uhk; Aaboud:2018juj. There is some overlapping sensitivity in these searches, but in general the sensitivity is greatly diminished away from the target process (see e.g. Aguilar-Saavedra:2018xpl; Aguilar-Saavedra:2017zuc; boosted_diboson for examples). It is not feasible to perform a dedicated analysis for every possible topology and so some signals may be missed. Global searches for new physics have been performed by the the LHC experiments and their predecesors, but only utilize simple objects and rely heavily on simulation for background estimation ATLAS-CONF-2017-001; CMS-PAS-EXO-10-021; ATLAS-CONF-2012-107; ATLAS-CONF-2014-006; Aktas:2004pz; Aaron:2008aa; Abbott:2000fb; Abbott:2000gx; Aaltonen:2007dg; Aaltonen:2008vt; sleuth; Knuteson:2004nj.

The tagging techniques used to isolate different jet types have increased in sophistication with the advent of modern machine learning classifiers Larkoski:2017jix; Cogan:2014oua; Almeida:2015jua; deOliveira:2015xxd; Baldi:2016fql; Barnard:2016qma; Kasieczka:2017nvn; Butter:2017cot; Komiske:2016rsd; Louppe:2017ipp; ATLAS-CONF-2017-064; ATL-PHYS-PUB-2017-013; ATL-PHYS-PUB-2017-004; CMS-DP-2017-005; CMS-DP-2017-013; ATL-PHYS-PUB-2017-003; Pearkes:2017hku; Datta:2017rhs; Datta:2017lxt; ATL-PHYS-PUB-2017-017; CMS-DP-2017-027; Fraser:2018ieu; Andreassen:2018apy; Macaluso:2018tck. These new algorithms can use all of the available information to achieve optimal classification performance and could significantly improve the power of hadronic resonance searches. Deep learning techniques are able to outperform traditional methods by exploiting subtle correlations in the radiation pattern inside jets. These correlations are not well-modeled in general Barnard:2016qma which renders classifiers sub-optimal when training on simulation and testing on data. This is already apparent for existing multivariate classifiers where post-hoc mis-modeling corrections can be large Aad:2015ydr; Chatrchyan:2012jua; Aad:2014gea; CMS:2013kfa; CMS-DP-2016-070; Aad:2015rpa; Khachatryan:2014vla; Aad:2016pux; CMS:2014fya. Ideally, one would learn directly from data (if possible) and/or combine with other approaches to mitigate potential mis-modeling effects during training (e.g. with adversaries Louppe:2016ylz).

We propose a new method that combines resonance searches with recently proposed techniques for learning directly from data Dery:2017fap; Metodiev:2017vrx; Komiske:2018oaa; Cohen:2017exh. Simply stated, the new algorithm trains a fully supervised classifier to distinguish a signal region from a mass sideband using auxiliary observables which are decorrelated from the resonance variable under the background-only hypothesis. A bump hunt is then performed on the mass distribution after applying a threshold on the classifier output. This is Classification Without Labels (CWoLa) Metodiev:2017vrx where the two mixed samples are the signal region and sideband and the signal is a potential new resonance and the background is the Standard Model continuum. The algorithm naturally inherits the property of CWoLa that it is fully based on data and thus is insensitive to simulation mis-modeling The algorithm also inherits the assumptions of the CWoLa method. In this context, the main assumption will be that signal region and sideband region can only be distinguished with the mass. More details on this are in the next sections. . The key difference with respect to Ref. Metodiev:2017vrx; Komiske:2018oaa is that the signal process need not be known a priori. Therefore, we can become sensitive to new signatures for which we did not think to construct dedicated searches.

In addition to CWoLa, the extended bump hunt shares some features with the sPlot technique Pivk:2004ty. Our proposed extension to the bump hunt makes use of auxiliary features to enhance the presence of signal events over background events in a target distribution, where the signal is expected to be resonant. Similarly, sPlot provides a procedure for using auxiliary features (‘discriminating variables’ in the language of Ref. Pivk:2004ty) to extract the distribution of signal and background events in a target distribution (‘control variable’ in Ref. Pivk:2004ty). In both cases, the auxiliary features must be uncorrelated with the target feature. One main difference between the methods is that the extended bump hunt uses machine learning to identify regions of phase space that are signal-like. A second key distinction between methods is that sPlot takes the distribution of the auxiliary features as input, whereas this information is not required for the extended bump hunt.

This paper is organized as follows. Section 2 formally introduces the CWoLa hunting approach and briefly discusses how auxiliary information can be useful for bump hunting. Then, Sec. 3 uses a simplified example to show how a neural network can be used to identify new physics structures from pseudodata. A complete procedure for applying the CWoLa hunting approach is given in Sec. 4. Finally, a realistic example based on a hadronic resonance search is presented in Sec. 5. Conclusions and future outlook are presented in Sec. 6.

Bump Hunting using Classification Without Labels

In a typical resonance search, events have at least two objects whose four-vectors are used to construct an invariant mass spectrum. The structure of these objects as well as other information in the event may be useful for distinguishing signal from background even though there may be no other resonance structures. Let $m_{\text{res}}$ be a random variable that represents the invariant mass. The distribution of $m_{\text{res}}$ given background is smooth while $m_{\text{res}}$ given signal is expected to be localized near some $m_{0}$. Let $Y$ be another random variable that represents all other information available in the events of interest. Define two sets of events:

where $\epsilon>\delta$. The value of $\delta$ is chosen such that $M_{1}$ should have much more signal than $M_{2}$ and the value of $\epsilon$ is chosen such that the distribution of $Y$ is nearly the same between $M_{1}$ and $M_{2}$. CWoLa hunting entails training a classifier to distinguish $M_{1}$ from $M_{2}$ using $Y$ and then performing a usual bump hunt on $m_{\text{res}}$ after placing a threshold on the classifier output. This procedure is then repeated for all mass hypotheses $m_{0}$. Note that nothing is assumed about the distribution of $Y$ other than that it should be nearly the same for $M_{1}$ and $M_{2}$ under the background-only hypothesis.

Ideally, $Y$ incorporates as much information as possible about the properties of the objects used to construct the invariant mass and their surroundings. The subsequent sections will show how this can be achieved with neural networks. To build intuition for the power of auxiliary information, the rest of this section provides analytic scaling results for a simplified bump hunt with the most basic case: $Y\in\{0,1\}$.

Suppose that we have two mass bins $M_{1}$ and $M_{2}$ and the number of expected events in each mass bin is $N_{b}$. Further suppose that the signal is in at most one of the $M_{i}$ (not required in general) and the expected number of signal events is $N_{s}$. A version of the bump hunt would be to compare the number of events in $M_{1}$ and $M_{2}$ to see if they are significantly different. As a Bernoulli random variable, $Y$ is uniquely specified by $\Pr(Y=1)$. Define $\Pr(Y=1|\text{background})=p$ and $\Pr(Y=1|\text{signal})=q$. The purpose of CWoLa hunting is to incorporate the information about $Y$ into the bump hunt. By only considering events with $Y=1$, the significance of the signal scales as $qN_{s}/\sqrt{N_{b}p}$. Therefore, the information about $Y$ is useful when $q>\sqrt{p}$.

More quantitatively, suppose that we declare discovery of new physics when the number of events with $Y=1$ in $M_{1}$ exceeds the number of events with $Y=1$ in $M_{2}$ by some amount. Under the background-only case, for $N_{b}\gg 1$, the difference between the number of events in $M_{1}$ and $M_{2}$ with $Y=1$ is approximately normally distributed with mean $0$ and variance $2N_{b}p$. If we want the probability for a false positive to be less than 5%, then the threshold value is simply $\sqrt{2N_{b}p}\times\Phi^{-1}(0.95)$, where $\Phi$ is the cumulative distribution function of a standard normal distribution. Ideally, we would like to reject the SM often when there is BSM, $N_{s}>0$. Figure 1 shows the probability to reject the SM for a one-bin search using $N_{b}=1000$ and $N_{s}=20$ for different values of $p$ as a function of $q$. The case $p=q=1$ corresponds to the standard search that does not gain from having additional information. However, away from this case, there can be a significant gain from using $Y$, especially when $p$ is small and $q$ is close to $1$. In the case where $Y$ is a truth bit, i.e. $p=1-q=0$, the SM is rejected as long as a single BSM event is observed. By construction, when $q\rightarrow 0$ (for $p>0$), the rejection probability is 0.05. Note that when $q<p$, only considering events with $Y=1$ is sub-optimal - this is a feature that is corrected in the full CWoLa hunting approach.

While the model used here is simple, it captures the key promise of CWoLa hunting that will be expanded upon in more detail in the next sections. In particular, the main questions to address are: how to find $Y$ and how to use the information about $Y$ once it is identified.

Illustrative Example: Learning to Find Auxiliary Information

The model setup described above and used for the rest of this section is depicted in Fig. 2. The numerical examples presented below use $N_{b}=10,000$, $N_{s}=300$, and $w=0.2$. Without using $Y$, these values correspond to $N_{s}/\sqrt{N_{b}}=3\sigma$. The ideal tagger (one that is optimal by the Neyman-Pearson lemma Neyman289) should reject all events outside of the square in the $(x,y)$ plane centered at zero with side length $w$. For the $N_{s}$ and $N_{b}$ used here, the expected significance of the ideal tagger is $15\sigma$. The goal of this section is to show that without using any truth information, the CWoLa approach can recover much of the discriminating power from a neural network trained in the $(x,y)$ plane. Note that optimal classifier is simply given by thresholding the likelihood ratio Neyman289 $p_{s}(Y)/p_{b}(Y)$; in this two-dimensional case it is possible to provide an accurate approximation to this classifier without neural networks. However, these approximations often do not scale well with the dimensionality and will thus be less useful for the realistic example presented in the next section. This is illustrated in the context of the CWoLa hunting in Fig. 3.

To perform CWoLa hunting, a neural network is trained on $(x,y)$ values to distinguish events in the mass sidebands from the signal region. Due to the simple nature of the example, it is also possible to easily visualize what the network is learning. A fully-connected feed-forward network is trained using the Python deep learning library Keras keras with a Tensorflow tensorflow backend. The network has three hidden layers with (256, 256, 64) nodes. The network was trained with the categorical cross-entropy loss function using the Adam algorithm adam with a learning rate of 0.003 and a batch size of 1024. The data are split into three equal sets, one used for training, one for validation, and one for testing. The training is terminated based on the efficiency of the signal region cut on the validation data at a fixed false-positive-rate of $2\%$ for the sideband data. If it fails to improve for 60 epochs, the training is halted and the network reverts to the last epoch for which there was a training improvement. This simple scheme is robust against enhancing statistical fluctuations but reduces the number of events used for the final search by a factor of three as only the classifier output on the test set is used for the bump hunt. In the physical example described later, a more complicated scheme maximizes the statistical power of the available data.

Visualizations of the neural network trained as described above are presented in Fig. 4. In the top two examples, the network finds the signal region and correctly estimates the magnitude of the likelihood ratio. In both these cases, the network also overtrains on a (real) fluctuation in the training data, despite the validation procedure. Such regions will tend to decrease the effectiveness of the classifier, since a given cut threshold will admit more background in the test data. In the bottom left example of Fig. 4, the network finds a function approximately monotonic to $h(x,y)$ but with different normalization – while the cost function would have preferred to optimize this network to reach $h(x,y)$, the validation procedure cut off the optimization when the correct shape to isolate the signal region had been found. Due to the nature of the cuts, there is no performance loss for this network, since crucially it has found the correct shape near the signal region. The last network fails to converge to the signal region, and instead focuses its attention on the fluctuation in the training data. The variation in the network performance illustrates the importance of training multiple classifiers and using schemes to mitigate the impact of statistical fluctuations in the training dataset.

Figure 5 shows the mass distribution in the three bins after applying successfully tighter threshold on the neural network output. Since $Y$ is not a truth bit, the data are reduced in both the signal region and the mass sidebands. For each threshold, the background expectation $\hat{n}_{b}$ assuming a uniform distribution is estimated by fitting a straight line to the mass sidebands. Then, the significance is estimated from the number of observed events in the signal region, $n_{o}$, via $\mathcal{S}\approx(n_{o}-\hat{n}_{b})/\sqrt{\hat{n}_{b}}$. Of the threshold presented, the maximum significance corresponds to the 5% efficiency with $\mathcal{S}\approx 10.8\sigma$. Even though the ideal significance is $15\sigma$, for the particular pseudodataset shown in Fig. 5, the ideal classifier significance is $13.9\sigma$.

We can study the behavior of our NN classifiers by looking at the significance generated by ensembles of models trained on signals of different strength, as shown in Fig. 6. The top histogram shows the significance for an ensemble of models trained on the example signal (blue) and on a control dataset with no signal (green). The control ensemble appears to be normally distributed around $s/\sqrt{b}=0$, while the example signal ensemble is approximately normally distributed around $12\sigma$ (compared to $13.9\sigma$ for the ideal cut), along with a small $O(5\%)$ population of networks that fail to find the signal. The middle histogram shows the effect of decreasing the size of the signal region $w_{s}$ while modifying $N_{s}$ to maintain an expected significance of $15\sigma$ with ideal cuts. When $w_{s}$ is decreased, the training procedure appears to have a harder time picking up the signal, possibly due to our choice of an operating point of $2\%$ false-positive rate for the sideband validation. For $w_{s}=0.1$ (green), about $50\%$ of the networks effectively find the signal. For $w_{s}=0.05$ (red), only about $5\%$ find the signal. The bottom plot shows the effect of increasing $w_{s}$ while keeping $N_{s}$ fixed, so that the strength of the signal decreases. When the size of the signal region is doubled to $w_{s}=0.4$ (green), giving a expected signifance of $7.5\sigma$, the network performs similarly to the $w_{s}=0.2$ example (blue). When the signal distribution is identical to the background distribution ($w_{s}=1.0$, red), there is on average a small decrease in performance compared to simply not using a classifier.

Full Method

The previous sections uses key elements of the full extended bump hunt but do not include all components, including the full background estimation and statistical analysis. This section gives a concrete prescription for applying the CWoLa hunting method in practice, which will be used in an explicit example in Sec. 5. The setup is as in the previous sections: there is feature $m_{\text{res}}$ where the signal is expected to be resonant and then a set of other features $Y$ that are uncorrelated with $m_{\text{res}}$, but potentially useful for distinguishing signal from background. It is important to state that while a detailed model of $Y$ is not required to perform the CWoLa hunting procedure that is described in the rest of the section, a limited model of $Y$ is required to ensure the correlations with $m_{\text{res}}$ are minimal. Such a model could come from simulation, from theory, or directly from a sufficiently signal-devoid data sample.

While in the presence of signal, the CWoLa hunting method would ideally learn systematic correlations between $m_{\text{res}}$ and $Y$, instead, it may focus on statistical fluctuations in the background distributions. A naive application of CWoLa directly on the data may produce bumps in $m_{\text{res}}$ by seeking local statistical excesses in the background distribution. This corresponds to a large look-elsewhere effect over the space of observables $Y$ – the classifier may search this entire space and find the selection with the largest statistical fluctuation. In Sec. 3, we took the approach of splitting the dataset into training, validation and test samples which eliminates this affect, since the statistical fluctuations in the three samples will be uncorrelated. However, applying this approach in practice would reduce the effective luminosity available for the search and thus degrade sensitivity. We therefore apply a cross-validation technique which allows all data to be used for testing while ensuring that event subsamples are never selected using a classifier that was trained on them. We split the events randomly, bin-by-bin, into five event samples of equal size. The first sample is set aside, and the first classifier is trained on the signal- and sideband-region events of the remaining four samples. This classifier may learn the statistical fluctuations in these event samples, but those will be uncorrelated with the fluctuations of the first sample. Applying the classifier to the set-aside event sample will then correspond to only one statistical test, eliminating the look elsewhere effect. By repeating this procedure five times – each time setting aside one $k$-fold for testing and four for training and validation, all the data can be used for the bump hunt by adding up the selected events from each $k$-fold.

Further technical details about the statistical methods can be found in Appendix A. Asymptotic formulae can be used to determine the local $p$-value of an excess, but such formulae must be validated using more computationally expensive methods for each application of CWoLa hunting, as is demonstrated in the appendix.

The main result following the application of the method from Sec. 4 is the local $p$-value. To determine the compatibility of the entire mass range with the no-resonance hypothesis, it is desirable to be able to compute a global $p$-value. In the result presented here, the mass bins were fixed ahead of time and were also non-overlapping. Therefore, it is relatively simple to estimate a global $p$-value using e.g. a Bonferroni correction. However, this is not ideal (over-conservative) when the mass bin width is scanned as part of the procedure. It is still possible to determine a global $p$-value, in the same spirit as the full bumphunter statistic Choudalakis:2011qn. This would require a significant computational overhead as a large number of neural networks would need to be trained for each of many pseudo-experiments. An additional trials factor would be associated with scanning the threshold fraction on the neural network output. In the simplest approach, a small number of well-separated working points would be chosen, such as 10%, 1%, and 0.1%. These should be sufficiently different that the three local $p$-values could be treated as independent. However, a finer scan would require a proper assessment of the global $p$-value using pseudo-experiments. It may be possible to significantly reduce the computational cost by estimating the correlation between mass windows and threshold fractions in order to properly account for the look-elsewhere-effect Gross:2010qma; VITELLS2011230.

One final remark is about how one would use CWoLa hunting to set limits. In the form described above, the CWoLa hunting approach is designed to find new signals in data without any model assumptions. However, it is also possible to recast the lack of an excess as setting limits on particular BSM models. Given a simulated sample for a particular model, it would be possible to set limits on this model by mixing the simulation with the data and training a series of classifiers as above and running toy experiments, re-estimating the background each time. This is similar to the usual bump hunt, except that there is more computational overhead because the background distribution is determined in part by the neural networks, and the distribution in expected signal efficiencies cannot be determined except by these toy experiments This complicates the legacy utility of the results, but it would be possible to tweak procedures like those advocated by RECAST Cranmer:2010hk in which neutral networks would be automatically trained for a new signal model.. In the absence of an excess, it is also possible to directly recast the results by taking the classifier trained on data with no significant signal. However, without a real excess, the classifier will have nothing to learn. Such a classifier will likely not be useful for any particular signal model. Therefore, while it is technically possible to do a standard re-interpretation of the results, the most powerful limit setting requires access to the data to retrain the neural networks for an injected signal.

Physical Example

This section uses a dijet resonance search at the LHC to show the potential of CWoLa hunting in a realistic setting. As discussed in Sec. 1, both ATLAS and CMS have a broad program targeting resonance decays into a variety of SM particles. Due to significance advances in jet substructure-based tagging Larkoski:2017jix, searches involving hadronic decays of the SM particles can be just as if not more powerful than their leptonic counterparts. The usual strategy for these searches is to develop dedicated single-jet classifiers, including These are the latest $\sqrt{s}=13$ TeV results - see references within to find the complete history. $W/Z$- CMS-DP-2015-043; ATLAS-CONF-2017-064, $H$- CMS-DP-2015-038; ATLAS-CONF-2016-039, top- CMS-DP-2015-043; CMS-PAS-JME-15-002; ATLAS-CONF-2017-064, $b$- CMS-DP-2017-012; ATL-PHYS-PUB-2017-013, and quark-jet taggers CMS-DP-2016-070; ATL-PHYS-PUB-2017-009. Simulated events with per-instance labels are used for training and then these classifiers are deployed in data. However, the best classifier in data may not be the best classifier in simulation. This problem is alleviated when learning directly from data.

Learning directly from data has another advantage - the decay products of a new heavy resonance may themselves be beyond the SM. If the massive resonance decays into new light states such as BSM Higgs bosons or dark sector particles that decay hadronically, then no dedicated SM tagger will be optimal Aguilar-Saavedra:2018xpl; Aguilar-Saavedra:2017zuc. A tagger trained to directly find non-generic-jet structure could find these new intermediate particles and thus also find the heavy resonance. This was the approach taken in Aguilar-Saavedra:2017rzt, but that method is fully supervised and so suffers the usual theory prior bias and potential sources of mismodelling. Here we will illustrate how the CWoLa hunting approach could be used instead to find such a signal. The next section (Sec. 5.1) describes the benchmark model in more detail, as well as the simulation details for both signal and background.

For a benchmark signal, we consider the process $pp\to W^{\prime}\to WX,X\to WW$, where $W^{\prime}$ and $X$ are a new vector and scalar particle respectively. This process is predicted, for example, in the warped extra dimensional construction of Agashe:2016rle; Agashe:2017wss; boosted_diboson. The typical opening angle between the two $W$ bosons resulting from the $X$ decay is given by $\Delta R(W,W)\simeq 4\,m_{X}/m_{W^{\prime}}$ for $2m_{W}\ll m_{X}\ll m_{W^{\prime}}$, and so the $X$ particle will give rise to a single large-radius jet in the hadronic channel when $m_{X}\lesssim m_{W^{\prime}}/4$. Taking the mass choices $m_{W^{\prime}}=3\;\text{TeV}$ and $m_{X}=400\;\text{GeV}$, the signal in the fully hadronic channel is a pair of large-radius jets $J$ with $m_{JJ}\simeq 3\;\text{TeV}$, one of which has a jet mass $m_{J}\simeq 80\;\text{GeV}$ and a two-pronged substructure, and the other has mass $m_{J}\simeq 400\;\text{GeV}$ with a four-prong substructure which often is arranged as a pair of two-pronged subjets.

Events are generated with Madgraph5_aMC@NLO Alwall:2014hca v2.5.5 to generate $10^{4}$ signal events, using a model file implementing the tensor couplings of Agashe:2017wss and selecting only the fully hadronic decays of the three $W$ bosons. The events are showered using Pythia 8.226 Sjostrand:2007gs, and are passed through the fast detector simulator Delphes 3.4.1 deFavereau:2013fsa. Jets are clustered from energy-flow tracks and towers using the FastJet Cacciari:2011ma implementation of the anti-$k_{t}$ algorithm Cacciari:2008gp with radius parameter $\Delta R=1.2$. We require events to have at least two ungroomed large-radius jets with $p_{T}>400\;\text{GeV}$ and $|\eta|<2.5$. The selected jets are groomed using the soft drop algorithm Larkoski:2014wba in grooming mode, with $\beta=0.5$ and $z_{\text{cut}}=0.02$. The two hardest groomed jets are selected as a dijet candidate, and a suite of substructure variables are recorded for these two jets. With the same simulation setup, $4.45\times 10^{6}$ Quantum Chromodynamic (QCD) dijet events are generated with parton level cuts $p_{T,\,j}>300\;\text{GeV}$, $|\eta_{j}|<2.5$, $m_{jj}>1400\;\text{GeV}$.

In order to study the behaviour of the CWoLa hunting procedure both in the presence and absence of a signal, we produce samples both with and without an injected signal. The events are binned uniformly in $\log(m_{JJ})$, with 15 bins in the range $2001\;\text{GeV}<m_{JJ}<4350\;\text{GeV}$.

2 Training a Classifier

In order to test for a signal with mass hypothesis $m_{JJ}\simeq m_{\text{res}}$, we construct a ‘signal region’ consisting of all the events in the three bins centered around $m_{\text{res}}$. We also construct a low- and a high-mass sideband consisting of the events in the two bins below and above the signal region, respectively. The mass hypothesis will be scanned over the range $2278\;\text{GeV}\leq m_{\text{res}}\leq 3823\;\text{GeV}$, to avoid the first and last bins that can not have a reliable background fit without constraints on both sides of the signal region. The signal region width is motivated by the width of the $m_{JJ}$ peak for the benchmark signal process described earlier. Because all particles in the process are very narrow, this width corresponds to the resolution allowed by the jet reconstruction and detector smearing effects and will be relevant for other narrow signal processes also. For processes giving rise to wider bumps, the width of the signal hypothesis could be scanned over just as we scan over the mass hypothesis. We will then train a classifier to distinguish the events in the signal region from those in the sideband on the basis of their substructure. The objective in constructing the training framework is that the classifier should be very poor (equal efficiency in signal region and sideband for any threshold) in the case that no signal is present in the signal region, but if a signal is present with unusual jet substructure then the classifier should be able to locate the signal and provide discrimination power between signal and SM dijet events.

The background is estimated by fitting the regions outside of the signal region to a smoothly falling distribution. In practice, this requires that the auxiliary information $Y$ is nearly independent of $m_{JJ}$; otherwise, the distribution could be sculpted. To illustrate the problem, consider a classifier trained to distinguish the sideband and signal regions using the observables $m_{J}$ and the N-subjettiness variable $\tau_{1}^{(2)}$ Thaler:2011gf. The ratio $m_{J}/\sqrt{\tau_{1}^{(2)}}$ is approximately the jet $p_{\text{T}}$, which is highly correlated with $m_{JJ}$ for the background. While it is often possible to find ways to decorrelate substructure observables Dolen:2016kst; Shimmin:2017mfk; Aguilar-Saavedra:2017rzt; Moult:2017okx, we take a simpler approach and instead select a basis of substructure variables which have no strong correlations with $m_{JJ}$. We will use the following set of 12 observables which does not provide learnable correlations with $m_{JJ}$ sufficient to create signal-like bumps in our simulated background dijet event samples, as we shall demonstrate later in this section

In our study, the classifiers used are dense neural networks built and trained using Keras with a TensorFlow backend. We use four hidden layers consisting of a first layer of 64 nodes with a leaky Rectified Linear Unit (ReLU) activation (using an inactive gradient of 0.1), and second through fourth layers of 32, 16, 4 nodes respectively with Exponential Linear Unit (ELU) activation clevert2015fast. The output node has a sigmoid activation. The first three hidden layers are regulated with dropout layers with 20% dropout rate JMLR:v15:srivastava14a. The neural networks are trained to minimize binary cross-entropy loss using the Adam optimizer with learning rate of 0.001, batch size of 20000, first and second moment decay rates of 0.8 and 0.99, respectively, and learning rate decay of $5\times 10^{-4}$. The training data is reweighted such that the low sideband has equal total weight to the high sideband, the signal region has the same total weight as the sum of the sidebands, and the sum of all events weights in the training data is equal to the total number of training events. This ensures that the NN output will be peaked around 0.5 in the absence of any signal, and ensures that low and high sideband regions contribute equally to the training in spite of their disparity in event rates.

3 Results

The fact that there is no significant bump in the left plot of Fig. 8 is an important method closure test. When deploying the CWoLa hunting approach in practice, we advocate to test the method in simulation in order to validate that there are no bump-catalyzing correlations in the selected classification features. A residual concern may be that there are correlations in the data which are not present in simulation. Residual correlations may come in two forms: process and kinematic. Process correlations occur when $Y$ depends on the production channel (e.g. $pp\rightarrow qq$ or $pp\rightarrow gg$) and $m_{JJ}$ also depends on the production channel; kinematic correlations are the case when $m_{JJ}$ is correlated with $Y$ given the process. Residual process correlations do not cause bumps because the $m_{JJ}$ distribution of each process type (aside from signal) is smoothly falling. Thus, even if the classifier can exactly pick out one process, no bumps will be artificially sculpted. Residual kinematic correlations could cause artificially bumps in the $m_{JJ}$ distribution. Physically, kinematic correlations occur because $Y_{i}$ is correlated with $p_{\text{T,i}}$. One way to show in data that residual kinematic correlations are negligible is to use a mixed sample in which pairs of jets from different events are combined. As long as the potential signal fraction is small, this mixed sample will have no resonance peak. While the features chosen in this section were designed to be uncorrelated with $m_{JJ}$ and not sculpt bumps, it may be possible to utilize correlated features in a modified CWoLa hunting procedure that includes systematic uncertainties for strong residual correlations. We leave studies of this possibility to future work.

The ability of the CWoLa approach to discriminate signal from background depends on the number of signal and background events in the signal and sideband regions. In Fig. 11, we keep the number of background events fixed but vary the size of the signal, and plot truth-label ROC curves for each example. This allows us to directly asses the performance of the taggers for the signal. For varying thresholds, the $x$-axis corresponds to the efficiency on true signal events in the signal region, $\epsilon_{S}$, while the $y$-axis represents the inverse efficiency on true QCD events in the signal region, $\epsilon_{B}$. The gray dashed lines labelled 1 to 32 indicate the significance improvement, $\epsilon_{S}/\sqrt{\epsilon_{B}}$, which quantifies the gain in statistical significance compared to the raw $m_{JJ}$ distribution with no cuts applied. In solid black we show the performance of a dedicated tagger trained with labelled signal and background events using a fully supervised approach. This gives a measure of the maximum achievable performance for this signal using the selected variables. A true dedicated tagger which could be used in a realistic dedicated search would be unlikely to reach this performance, since this would require careful calibration over 12 substructure variables with only simulated data available for the signal. While the CWoLa-based taggers do not reach the supervised performance in these examples, we find that performance does gradually improve with increasing statistics.

We also show in the dashed black curve the performance of $W$/$Z$ tagger in identifying this signal for which the tagger is not designed. This tagger is trained on a sample of $pp\to W^{\prime}\to WZ$ events in the fully hadronic channel. In this case, the tagger is trained on the individual $W$/$Z$ jets themselves rather than the dijet event, as is typical in the current ATLAS and CMS searches. In producing the ROC curve, dijet events are considered to pass the tagging requirement if both large-radius jets pass a threshold cut on the output of the $W$/$Z$-tagger. We see that for $\epsilon_{B}\sim 10^{-4}$, which is a typical background rejection rate for recent hadronic diboson searches, the signal rate is negligible since the $X$-jet rarely passes the cuts. This illustrates that CWoLa hunting may find unexpected signals which are not targeted by existing dedicated searches is $S/B$ if high enough. If $S/B$ is too low, then the CWoLa hunting approach is not able to identify the signal and it underperforms compared with the search that is targeting a different signal model.

The datasets and code used for the case study can be found at Refs. cwola_hunting_dataset; cwola_hunting_code.

Conclusions

We have presented a new anomaly detection technique for finding BSM physics signals directly from data. The central assumption is that the signal is localized as a bump in one variable in which the background is smooth, and that other features are available for additional discrimination power. This allows us to identify potential signal-enhanced and signal-depleted event samples with almost identical background characteristics on which a classifier can be trained using the Classification Without Labels approach. In the case that a distinctive signal is present, the trained classifier output becomes an effective discriminant between signal events and background events, while in the case that no signal is present the classifier output shows no clear pattern. An event selection based on a threshold cut on the classifier output produces a smooth distribution if no signal is present and produces a bump if a signal is present, and so standard bump hunting techniques can be used on the selected distribution.

The prototypical example used here is the dijet resonance search in which the dijet mass is the one-dimensional feature where the signal is localized. Related quantities could also be used, such as the single jet mass for boosted resonance searches Sirunyan:2017dnz; Sirunyan:2017nvi; Aaboud:2018zba or the average mass of pair produced objects Aaboud:2017nmi; CMS:2018sek; ATLAS:2012ds; Aad:2016kww; Chatrchyan:2013izb; Khachatryan:2014lpa. Jet substructure information was used to augment information from just the dijet mass and a CWoLa classifier was trained using a deep neural network to discriminate signal region events from sideband events based on their substructure distributions. Additional local information such as the number of leptons inside the jets, the number of displaced vertices, etc. could be used in the future to ensure sensitivity to a wide variety of models. Furthermore, event-level information such as the number of jets or the magnitude of the missing transverse momentum could be added to an extended CWoLa hunt.

The CWoLa hunting strategy is generalizable beyond this single case study. To summarize, the essential requirements are:

There is one bump-variable $m_{\text{res}}$ in which the background forms a smooth distribution, for which there is a background model such a parametric function, and a signal can be expected to be localized as a bump. This was the variable $m_{JJ}$ in the dijet case study.

There are additional features $Y$ in the events which may potentially provide discriminating power between signal and background, but the detailed topology of the of the signal in these variables is not known in advance. This was the set of substructure variables in the dijet study.

The background distribution in $Y$ should not have strong correlations with $m_{\text{res}}$ over the resonance width of the signal. In the case that such correlations exist, it may be possible to find a transformation of the variables that removes these correlations before being fed into the classifier, or alternatively to train the classifier in such a way that penalizes shaping of the $m_{\text{res}}$ distribution outside of the signal region. Closure tests in simulation or with mixed samples in data can be used to confirm that $Y$ is not strongly correlated with $m_{\text{res}}$.

By harnessing the power of modern machine learning, CWoLa hunting and other weakly supervised strategies may provide the key to uncovering BSM physics lurking in the unique datasets already collected by the LHC experiments.

Appendix A Statistical Analysis

The significance of the bump is evaluated in the following way. First, after selecting the signal-like events in each cross-validation sample using its corresponding classifier, we merge the selected events from the $k$ samples into a single selected dataset. This dataset is binned in $m_{JJ}$, and we estimate the background by fitting a smooth, parametric function to the dataset with the signal region masked out. We use the following three-parameter function (also used in the ATLAS Aaboud:2017eta and CMS Sirunyan:2016cao searches for fully hadronic diboson resonances More complex procedures for fitting the background such as Gaussian processes are also possible Frate:2017mai but their use is beyond our scope.):

which is fitted using a least-squares fit. The number of events in the signal region is predicted by summing the predictions in each of the three signal region bins. The systematic uncertainty in this fit is estimated by propagating linearly the uncertainties on the fit parameters onto an uncertainty in the signal region prediction. The fits and fit uncertainties are indicated in the left plot in Fig. 8 by the red dashed lines and gray bands. We tested the goodness of fit of this functional form in background-only simulations using Kolmogorov–Smirnov tests, and in any real search we would advocate similar tests in simulation. In the case that simulation is not completely reliable, it is possible to define data validation regions using non-signal selections in order to provide a cross-check of the fit function, as is done in e.g. Ref. Aaboud:2017eta. In the case of CWoLa hunting, this would entail selecting events in non-signal windows of the classifier output. For example, if using a 1% selection for the signal search, one could use other percentile windows of the NN output to define non-signal selections with similar statistics which should be well fitted by the fit function under both the null and alternate hypotheses.

Since the shape of the signal in the signal region is a-priori unknown we base our hypothesis test on the total number of events in the signal region. We form the profile likelihood ratio

where $\mu$ indicates the signal rate and $\theta$ is the nuisance parameter associated with the systematic uncertainties on the background prediction. In the numerator, $\hat{\vphantom{\rule{1.0pt}{5.71527pt}}\smash{\hat{\theta}}}$ represents the best fit value for the nuisance parameter in the background-only hypothesis $\mu=0$, while in the denominator $\hat{\mu}$ and $\hat{\theta}$ represent the combined best fit for $\mu$ and $\theta$. The likelihood is formed from a product of a Poisson factor for the number of events in the signal region, and a Gaussian constraint for the background nuisance parameter

where $n$ is the observed number of events in the signal region, $b$ is the number of background events predicted by the sideband fit, $\theta$ is the nuisance parameter associated with the systematic uncertainty for the background prediction, and $\sigma$ is the uncertainty on that nuisance parameter.

Using asymptotic formulae Cowan:2010js, gives a significance $Z=\sqrt{q_{0}}$ and $p_{0}=1-\Phi(Z)$, where $\Phi$ is the cumulative distribution function of the normal distribution.

The null hypothesis is that the dijet invariant mass distribution after selection by the classifiers is well described by the smooth functional form of Equation (4). This requires that prior to any classification, the spectrum is smooth (already assumed by ATLAS and CMS) and that the classifiers are not able to generate localized features in the mass distribution following the CWoLa hunting procedure. In order to use the asymptotic formulae from Ref. Cowan:2010js, the bin counts in the selected, merged datasets must be Poissonian. The rest of the section investigates the validity of these approximations.

Let $f(x)$ represent the function described by Eq. 4 prior to any classification and consider a dataset with $N_{i}^{\text{(uncut)}}$ events in $m_{JJ}$ bin $i$ (bin center $m_{\text{res},\,i}$) with $N_{i}^{\text{(uncut)}}\sim\text{Poiss}(f(m_{\text{res},\,i}))$. Let $Y$ be a set of auxiliary observables whose probability distributions are independent of $m_{\text{res}}$. The goal is to demonstrate that the $p$-values reported from the statistical procedure described above are accurate. To begin, the dataset with $N_{i}^{\text{(uncut)}}$ events is partitioned into $k$ samples with equal probability for an event to be assigned to one of the samples. Next, a classifier is trained to discriminate signal region events from sideband events using all subsamples except the $j$’th. The classifier is then used to select a fraction $\epsilon$ of events in the held out $j$’th sample, using all other bins to determine $\epsilon$ The number of events used to determine $\epsilon$ is sufficiently large that the uncertainty on the value of the NN used to achieve $\epsilon$ efficiency is negligible.. This means that

In the case that the cross-validated selected event rates in the $k$ samples are uncorrelated, then it would follow that after merging these datasets the total selected event rate distribution would be given by

However, because the events in one sample are used to train a classifier applied to the other samples, it cannot necessarily be assumed that the event rates are uncorrelated between samples. If strong correlations are expected between selected samples then in order to calculate reliable $p$-values the test statistic would need to be calibrated by running many toy experiments on either new simulated event samples or on bootstrapped samples, with the NNs trained fresh each time. Since this is a computationally expensive procedure, it is preferable if a simpler alternative is available.

In order to check that the simpler approach (assuming no correlations between cross-validated samples) is valid, we have performed an empirical test of this effect in the following way. We generated $10^{3}$ toy datasets with binned event counts drawn from Poisson distributions with means determined by the distribution of Equation (4), with parameters obtained by a fit to the uncut dijet dataset used in Section. 5. Each event has 12 auxialliary variables, as in Section 5, but with these variables drawn from a random uniform distribution in the range $$. NNs were trained using a cross validation procedure exactly as described in Section 5, except for the following modifications that were required to reduce the computational time required. We used 4-fold cross validation (rather than 5-fold), trained only four NNs per iteration from which the best was selected (rather than 20), and the NNs were trained with a patience of 100 epochs of no improvement in validation performance before stopping (instead of 300 epochs). The trained NNs were used to select the 1% most ‘signal-like’ events for each toy. For each toy we then calculated the test statistic for rejection of the null hypothesis, and the distribution of these test statistics is shown by the black markers in Figure 12.

Additionally, we generated $10^{5}$ toy datasets with $m_{\text{res}}$ drawn in the same way. Instead of training NNs to select events, we randomly selected 1% of events. For each toy we then calculated the test statistic for rejection of the null hypothesis, and the distribution of these test statistics is shown by the orange histogram in Figure 12. Finally, we show the expected asymptotic distribution with the blue line.

The key feature of Figure 12 is that the test statistic distribution for the NN toys shows no apparent deviation from that for the simple toys or from the asymptotic form. We therefore find no evidence of any distortion caused by correlations between cross-validation samples in this toy experiment.

Finally, it is worth remarking that the $p$-value computed with the above procedure is only local. If a local $p$-value is below some threshold, a followup, dedicated analysis using an orthogonal dataset should target the identified region of phase space with no trials factor penalty. One could also estimate a global $p$-value in a standard way using e.g. a Bonferroni correction. Other methods like the full bumphunter statistic could be used Choudalakis:2011qn but that is not the standard practice in the current ATLAS and CMS diboson resonance searches.

Appendix B Dijet Mass Scans

In Figs. 13, 14 we plot the dijet invariant mass distributions before and after applying tagger cuts over the full range of the mass scan described in Sec. 5. The $p$-values calculated from the top four distributions in these plots are displayed in Fig. 8 (right).

References