Adversarially Learned Anomaly Detection on CMS Open Data: re-discovering the top quark

We apply an Adversarially Learned Anomaly Detection (ALAD) algorithm to the problem of detecting new physics processes in proton-proton collisions at the Large Hadron Collider. Anomaly detection based on ALAD matches performances reached by Variational Autoencoders, with a substantial improvement in some cases. Training the ALAD algorithm on 4.4 fb-1 of 8 TeV CMS Open Data, we show how a data-driven anomaly detection and characterization would work in real life, re-discovering the top quark by identifying the main features of the t-tbar experimental signature at the LHC.


Introduction
CERN's Large Hadron Collider (LHC) delivers proton-proton collisions in unique experimental conditions. Not only it accelerates protons to an unprecedented energy (6.5 TeV for each proton beam). It also operates at the highest collision frequency, producing one proton-beam crossing (an "event") every 25 nsec. By recording the information of every sensor, the two LHC multipurpose detectors, ATLAS and CMS, generate O(1) MB of data for each interesting event. The LHC big data problem consists in the incapability of storing each event, that would require writing O(10) TB/sec. In order to deal with this limitation, a typical LHC detector is equipped with a real-time event selection system, the trigger.
The need of a trigger system for data taking has important consequences downstream. In particular, it naturally shapes the data analysis strategy for many searches for new physics phenomena (new particles, new forces, etc.) into a hypothesis test [1]: one specifies a signal hypothesis upfront (in fact, as upfront as the trigger system) and tests the presence of the predicted kind of events (the experimental signature) on top of a background from known physics processes, against the alternative hypothesis (i.e., known physics processes alone). From a data-science perspective, this corresponds to a supervised strategy. This procedure was very successful so far, thanks to well established signal hypotheses to test (e.g., the existence of the Higgs boson). On the other hand, following this paradigm didn't produce so far any evidence for new-physics signals. While this is teaching us a lot about our Universe, 1 it also raises questions on the applied methodology. trained against each other in a saddle-point problem: min where p X (x) is the distribution over the data space X and p Z (z) is the distribution over the latent space Z. The solution to this problem will have the property p X = p G , where p G is the distribution induced by the generator [19]. The training typically involves alternating gradient descent on the parameters of G and D x to maximize for D x (treating G as fixed) and to minimize for G (treating D x as fixed).
Deep learning for anomaly detection [20] is usually discussed in the context of (variational) autoencoders [21,22]. With autoencoders (cf. Fig. 1), one projects the input x to a point z of a latent-space through an encoder network E : X → Z. An approximation D(z) = D(E(x)) of the input information is then reconstructed through the decoder network, D : Z → X . The intuition is that the decoder D can only reconstruct the input from the latent space representation z if x ∼ p X . Therefore, the reconstruction for an anomalous sample, which belongs to a different distribution, would typically have a higher reconstruction loss. One can then use a metric D R defining the output-to-input distance (e.g., the one used in the reconstruction loss function) to derive an anomaly-score A: (2) While this is not directly possible with GANs, since a generated G(z) doesn't correspond to a specific x, several GANbased solutions have been proposed that would be suitable for anomaly detection, as for instance in Refs. [18,20,[27][28][29].
In this work, we focus on the ALAD method [18], built upon the use of bidirectional-GANs (BiGAN) [30]. As shown in Fig. 1, a BiGAN model adds an encoder E : X → Z to the GAN construction. This encoder is trained simultaneously to the generator. The saddle point problem in Eq.(1) is then extended as follows: min where D xz is a modified discriminator, taking inputs from both the X and Z. Provided there is convergence to the global minimum, the solution has the distribution matching property [30]. To help reaching full convergence, the ALAD model is equipped with two additional discriminators: D xx and D zz . The former discriminator together with the value function enforces the cycle-consistency condition G(E(x)) ≈ x. The latter is added to further regularize the latent space through a similar value function: enforcing the cycle condition E(G(z)) ≈ z. The ALAD training objective consists in solving: min Having multiple outputs at hand, one can associate the ALAD algorithm to several anomaly-score definitions. Following Ref. [18], we consider the following four anomaly scores: • A "Logits" score, defined as: A L (x) = log(D xx (x, G(E(x))).
• A "Features" score, defined as: are the activation values in the last hidden layer of D xx . • The L 1 distance between an input x and its reconstructed output G(E(x)): A L1 (x) = ||x − G(E(x))|| 1 .
• The L 2 distance between an input x and its reconstructed output G(E(x)): We first apply this model to the problem described in Ref [3], in order to obtain a direct comparison with VAE-based anomaly detection. Then, we apply this model to real LHC data (2012 CMS Open Data), showing how anomaly detection could guide physicists to discover and characterize new processes.

ALAD Performance Benchmark
We consider a sample of simulated LHC collisions, pre-filtered by requiring the presence of a muon with large transverse momentum (p T ) 2 and isolated from other particles. Proton-proton collision events at a center-of-mass energy √ s = 13 TeV are generated with the PYTHIA8 event-generation library [31]. The generated events are further processed with the DELPHES library [32] to model the detector response. Subsequently the DELPHES particle-flow (PF) algorithm is applied to obtain a list of reconstructed particles for each event, the so-called PF candidates.
Events are filtered requiring p T > 23 GeV and isolation 3 Iso < 0.45. Each collision event is represented by 21 physics-motivated high-level features (see Ref. [3]). These input features are pre-processed before being fed to the ALAD algorithm. The discrete quantities 4 (q l , IsEle, N µ and N e ) are represented through one-hot encoding. The other features are standardized to a zero median and unit variance. The resulting vector, containing the one-hot encoded and continuous features, has a dimension of 39 and is given as input to the ALAD algorithm.
A SM cocktail is assembled from a weighted mixture of those four processes, with weights given by the production cross section. This cocktail's composition is given in Table 1. Table 1: Composition of the cocktail of SM processes. The first column gives the production cross section for each process. The trigger efficiency is the fraction of events passing a loose selection on the muon p T , corresponding to an online selection in the trigger system. The acceptance is the fraction of events passing our selection criteria (p T > 23 GeV and Iso < 0.45). Each fraction is computed with respect to the previous step. The last column gives the composition of the SM-cocktail.

Process
Cross We train our ALAD model on this SM cocktail and subsequently apply it to a test dataset, containing a mixture of SM events and events of physics beyond the Standard Model (BSM). In particular, we consider the following BSM datasets, also available on Zenodo: • A leptoquark with mass 80 GeV, decaying to a b quark and a τ lepton: LQ → bτ [38].
As a starting point, we consider the ALAD architecture [18] used for the KDD99 dataset, which has similar dimensionality as our input feature vector. In this configuration, both the D xx and D zz discriminators take as input the concatenation of the two input vectors, which is processed by the network up to the single output node, activated by a sigmoid function. The D xz discriminator has one dense layer for each of the two inputs. The two intermediate representations are concatenated and passed to another dense layer and then to a single output node with sigmoid activation, as for the other discriminators. The hidden nodes of the generator are activated by ReLU functions [42], while Leaky ReLU [43] are used for all the other nodes. The slope parameter of the Leaky ReLU function is fixed to 0.2. The network is optimized using the Adam [44] minimizer and minibatches of 50 events each. The training is regularized using dropout layers in the three discriminators.
Starting from this baseline architecture, we adjust the architecture hyperparameters one by one, repeating the training while maximizing a figure of merit for anomaly detection efficiency. We perform this exercise using as anomalies the benchmark models described in Ref. [3] and looking for a configuration that performs well on all of them. To quantify performance, we consider both the area under the receiver operating characteristic (ROC) curve and the positive likelihood ratio LR + . We define the LR + as the ratio between the BSM signal efficiency, i.e., the true positive rate (TPR), and the SM background efficiency, i.e., the false positive rate (FPR). The training is performed on half of the available SM events (3.4M events), leaving the other half of the SM events and the BSM samples for validation. From the resulting anomaly scores, we compute the ROC curve and compare it to the results of the VAE in Ref. [3]. We further quantify the algorithm performance considering the LR + values corresponding to an FPR of 10 −5 .
The optimized architecture, adapted from Ref. [18], is summarized in Table 2. This architecture is used for all subsequent studies. We consider as hyperparameters the number of hidden layers in the five networks, the number of nodes in each hidden layer, and the dimensionality of the latent space, represented in the table by the size of the E output layer. Having trained the ALAD on the training dataset, we compute the anomaly scores for the validation samples as well as for the four BSM samples, where each BSM process has O(0.5M) samples. Figure 2 shows the ROC curves of each BSM benchmark process, for the four considered anomaly scores. The best VAE result from Ref. [3] is also shown for comparison. In the rest of this paper, we use the L 1 score as the anomaly score. Similar results would have been obtained using any of the other three anomaly scores. Figure 3 compares the A L1 distribution for each BSM process with the SM cocktail. One can clearly see that all BSM processes have an increased probability in the high-score regime The VAE curve corresponds to the best result of Ref. [3], which is shown here for comparison. The other four lines correspond to the different anomaly score models of the ALAD.
compared to the SM cocktail. We further verified that the anomaly score distributions obtained on the SM-cocktail training and validation sets are consistent. This test excludes the occurrence of over-training issues.
The ALAD algorithm outperforms the VAE by a substantial margin on the A → 4 sample, providing similar performance overall, and in particular for FPR ∼ 10 −5 , the working point chosen as a reference in Ref. [3]. We verified that the uncertainty on the TPR at fixed FPR, computed with the Agresti-Coull interval [45], is negligible when compared to the observed differences between ALAD and VAE ROC curves, i.e., the difference is statistically significant.
The left plot in Fig. 4 provides a comparison across different BSM models. As for the VAE, ALAD performs better on A → 4 and h ± → τ ν than for the other two BSM processes. The right plot in Figure 4 shows the LR + values as a function of the FPR ones. The LR + peaks at a SM efficiency of O(10 −5 ) for all four BSM processes and is basically constant for smaller SM-efficiency values.

Re-discovering the top quark with ALAD
In order to test the performance of ALAD on real data, and in general to show how an anomaly detection technique could guide physicists to a discovery in a data-driven manner, we consider a scenario in which collision data from the LHC are available, but no previous knowledge about the existence of the top quark is at hand. The list of known SM  processes contributing to a dataset with one isolated high-p T lepton would then include W production, Z/γ * production and QCD multijet events, neglecting more rare processes such as diboson production. Top-quark pair production, the "unknown anomalous process", represents ∼ 0.1% of the total dataset.
We consider a fraction of LHC events collected by the CMS experiment in 2012 and corresponding to an integrated luminosity of about 4.4 fb −1 at a center-of-mass energy √ s = 8 TeV. For each collision event in the dataset, we compute a vector of physics-motivated high-level features, which we give as input to ALAD for training and inference. Details on the dataset, event content, and data processing can be found in Appendix A.
We select events applying the following requirements: • At least one lepton with p T > 23 GeV and PF isolation Iso < 0.1 within |η| < 1.4.
• At least two jets with p T > 30 within |η| < 2.4. This selection is tuned to reduce the expected QCD multijet contamination to a negligible level, which avoids problems with the small size of the available QCD simulated dataset. This selection should not be seen as a limiting factor for the generality of the study: in real life, one would apply multiple sets of selection requirements on different triggers, in order to define multiple datasets on which different anomaly detection analyses would run.
Our goal is to employ ALAD to isolate tt events as due to a rare (and pretended to be) unknown process in the selected dataset. Unlike the case discussed in Ref [3] and Section 3, we are not necessarily interested in a pre-defined algorithm to run on a trigger system. Given an input dataset, we want to isolate its outlier events without explicitly specifying which tail of which kinematic feature one should look at. Because of this, we don't split the data into a training and a validation dataset. Instead, we run the training on the full dataset. The training is performed fixing the ALAD architecture to the values shown in Table 2. In our procedure, we implicitly rely on the assumption that the anomalous events are rare, so that their modeling is not accurately learned by the ALAD algorithm.
In order to evaluate the ALAD performance, we show in Fig. 5 the ROC and LR + curves on labelled Monte Carlo (MC) simulated data, considering the tt sample as the signal anomaly and the SM W and Z production as the background. An event is classified as anomalous whenever the L 1 anomaly score is above a given threshold. The threshold is set such that the fraction of selected events is about 10 −3 . The anomaly selection results in a factor-20 enhancement of the tt contribution over the SM background, for the anomaly-defining FPR threshold. Figure 6 shows the distributions of a subset of the input quantities before and after anomaly selection (see Appendix A): H T , M J , N J , and N b . These are the quantities that will become relevant in the rest of the discussion. The corresponding distributions for MC simulated events are shown in Figure 7, where the tt contribution to the selected anomalous data is highlighted. At this stage, we don't attempt a direct comparison of simulated distributions to data, since we couldn't apply the many energy scale and resolution corrections, normally applied for CMS studies. Anyhow, such a comparison is not required for our data-driven definition of anomalous events.
In order to quantify the impact of the applied anomaly selection on a given quantity, we consider the bin-by-bin fraction of accepted events. We identify this differential quantity as anomaly selection transfer function (ASTF). When compared to the expectation from simulated MC data, the dependence of ASTF values on certain quantities allows one to characterize the nature of the anomaly, indicating in which range of which quantity the anomalies cluster together. Figure 8 shows the ASTF for the data and the simulated SM processes (W and Z production) defining the background data. The comparison between the data and the background ASTF suggests the presence of an unforeseen class of events, clustering at large number of jets. An excess of anomalies is observed at large jet multiplicity, which also induces an excess of anomalies at large values of H T and M J . Notably, a large fraction of anomalous events has jets originating from b quarks. This is the first time that a MC simulation enters our analysis. We stress the fact that this comparison between data and simulation is qualitative. At this stage, we don't need the MC-predicted ASTF values to agree with data in absence of a signal, since we are not attempting a background estimate like those performed in data-driven searches at the LHC. For us, it is sufficient to observe qualitative differences in the dependence of ASTFs on specific features. Nevertheless, a qualitative difference like those we observe in Fig. 8 could still be induced by systematic differences between data and simulation. That would still be an anomaly, but not of the kind that would lead to a discovery. In our case, we do know that the observed discrepancy is too large to be explained by a wrong modeling of W and Z events, given the level of accuracy reached by the CMS simulation software in Run I. In a real-life situation, one would have to run further checks to exclude this possibility.
Starting from the ASTF plots of Fig. 8, we define a post-processing selection (PPS), aiming to isolate a subset of anomalous events in which the residual SM background could be suppressed. In particular, we require N J ≥ 6 and N b ≥ 2. Figure 9 shows the distributions of some of the considered features after the PPS. According to MC expectations, almost all background events should be rejected. The same should apply to background events in data. Instead, a much larger number of events is observed. This discrepancy points to the presence of an additional process, characterized by many jets (particularly coming from b jets). As a closure test, we verified that the agreement between the observed distributions in data and the expectation from MC simulation is restored once the tt contribution is taken into account.
In summary, the full procedure consists of the following steps: 1. Define an input dataset with an event representation which is as generic as possible.
2. Train the ALAD (or any other anomaly detection algorithm) on it.
3. Isolate a fraction α of the events, by applying a threshold on the anomaly score.
4. Study the ASTF on as many quantities as possible, in order to define a post-processing selection that allows one to isolate the anomaly.
5. (Optionally) once a pure sample of anomalies is isolated, a visual inspection of the events could also guide the investigation, as already suggested in Ref. [3]. clustering on specific run periods would most likely be due to transient detector malfunctioning or experimental problems of other nature (e.g., a bug in the reconstruction software). If instead the significant ASTF ratios point to events with a lepton, large jet multiplicity, with an excess of b-jets, one might have discovered a new heavy particle decaying to leptons and b-jets. In a world with no previous knowledge of the top quark, some bright theorist would explain the anomaly proposing the existence of a third up-type quark with unprecedented heavy mass.

Conclusions
We presented an application of ALAD to the search for new physics at the LHC. Following the study presented in Ref. [3], we show how this algorithm matches (and in some cases improves) the performance obtained with variational autoencoders. The ALAD architecture also offers practical advantages with respect to the VAE of Ref. [3]. Besides providing an improved performance in some cases, it offers an easier training procedure. Good performance can be achieved using a standard MSE loss, unlike the VAE model in Ref. [3], for which good performance was obtained only after a heavy customization of the loss function.
We train the ALAD algorithm on a sample of real data, obtained by processing part of the 2012 CMS Open data. On these events, we show how one could detect the presence of tt events, a < 0.1% population pretended to originate from a previously unknown process. Using the anomaly score of the trained ALAD algorithm, we define a post-selection procedure that let us isolate an almost pure subset of anomalous events. Furthermore, we present a strategy based on ASTF distributions to characterize the nature of an observed excess and we show its effectiveness on the specific example at hand. Further studies should be carried on to demonstrate the robustness of this strategy for more rare processes.
For the first time, our study shows with real LHC data that anomaly detection techniques can highlight the presence of rare phenomena in a data-driven manner. This result could help promoting a new class of data-driven studies in the next run of the LHC, possibly offering additional guidance in the search for new physics. As already stressed in Ref [3], the promotion of these algorithms to the trigger system of the LHC experiments could allow to mitigate for the limited bandwidth of the experiments. Unlike Ref [3], the strategy presented in this paper would also require a substantial fraction of normal events to be saved (through pre-scaled triggers), in order to compute the ASTF ratios with sufficient precision. In this respect, there would be a cost in terms of trigger throughput, that could be compensated by foreseeing dedicated scouting streams [46][47][48] with reduced event content, tailored to the computation of the ASTF ratios in offline analyses. Putting in place this strategy for the LHC Run III could be an interesting way to extend the physics reach of the LHC experiments.

A CMS Single Muon Open Data
The re-discovery of the top quark, described in Section 4, is performed on a sample of real LHC collisions, collected by the CMS experiment and released on the CERN Open Data portal [24]. The collision data correspond to the SingleMu dataset in the Run2012B period [49]. This dataset results from the logic OR of many triggers requiring a reconstructed muon. Most of the events were collected by an inclusive isolated-muon trigger, with requirements very similar to those applied as a pre-selection in Section 3 [3].
Once collected, the raw data recorded by the CMS detector were processed by a global event reconstruction algorithm, based on Particle Flow (PF) [50]. This processing step gives as output a list of so-called PF candidates (electrons, muons, photons, charged hadrons and neutral hadrons), reconstructed combining measurements from different detectors to achieve maximal accuracy.
While the ALAD training is performed directly on data, samples from MC simulation are also used in the study, in order to validate the model and its findings on a labelled sample. To this purpose, we considered samples of W +jets [51][52][53], Z/γ * +jets [54][55][56][57], and tt [58] events. No sample of QCD multijet events was considered, since its contribution was found to be negligible after the baseline requirements on the muon and on two additional jets (see Section 4). These MC samples were generated by the CMS collaboration with different libraries and processed by a full Geant4 simulation [59]. The same reconstruction software used on data was applied to the output of the simulation, so that the same lists of PF candidates are available in this case.
Following a procedure similar to that of Ref. [3], we take as input the lists of PF candidates and compute a set of physics motivated features, which are used as input to train the ALAD algorithm: We consider 14 event-related quantities: • H T -The scalar sum of the transverse momenta (p T ) of all jets having p T > 30 GeV and |η| < 2.4.
• M J -The invariant mass of all jets entering the H T sum.
• N J -The number of jets entering the H T sum.
• N B -The number of jets identified as originating from a b quark. We further consider 10 quantities, specific to the highest-p T lepton in the event: • p T -The lepton p T .
• q -The lepton charge (either −1 or +1) • Iso ch -The lepton isolation, defined as the ratio between the scalar sum of the p T of the other charged PF candidates with angular distance ∆R = ∆η 2 + ∆φ 2 < 0.3 from the lepton, and the lepton p T .
• Iso neu -Same as Iso ch but with the sum going over all neutral hadrons.
• Iso γ -Same as Iso ch but with the sum going over all photons.
• Iso -The total isolation given by Iso = Iso ch + Iso neu + Iso γ .
• M T -The combined transverse mass of the lepton and the E miss T system, which is given by M T = 2p T E miss T (1 − cos ∆φ). .

• p miss
T, -The parallel component of p miss T with respect to the lepton.
• p miss T,⊥ -The orthogonal component of p miss T with respect to the lepton.
• IsEle -A flag set to 1 if the lepton is an electron, 0 if it is a muon.