aainstitutetext: Département de Physique Nucléaire et Corpusculaire,
University of Geneva, 1211 Geneva, Switzerlandbbinstitutetext: Centre
Universitaire d’Informatique, University of Geneva, 1211 Geneva,
Switzerlandccinstitutetext: Institute for Theoretical Particle Physics and
Cosmology, RWTH Aachen University, 52074 Aachen, Germany

# Cluster Scanning: a novel approach to resonance searches

I. Oleksiyuk a J. A. Raine c M. Krämer b S. Voloshynovskiy a T. Golling
<EMAIL_ADDRESS><EMAIL_ADDRESS><EMAIL_ADDRESS><EMAIL_ADDRESS><EMAIL_ADDRESS>
###### Abstract

We propose a new model-independent method for new physics searches called
Cluster Scanning. It uses the k-means algorithm to perform clustering in the
space of low-level event or jet observables, and separates potentially
anomalous clusters to construct a signal-enriched region. The invariant mass
spectra in these two regions are then used to determine whether a resonant
signal is present. A pseudo-analysis on the LHC Olympics dataset with a
$Z^{\prime}$ resonance shows that Cluster Scanning outperforms the widely used
4-parameter functional background fitting procedures, reducing the number of
signal events needed to reach a $3\sigma$ significant access by a factor of
0.61. Emphasis is placed on the speed of the method, which allows the test
statistic to be calibrated on synthetic data.

##  1 Introduction

The Standard Model (SM) is the current apex of theoretical physics, describing
the electromagnetic, weak and strong interactions with unparalleled precision.
Unfortunately, it is still far from complete, as several phenomena remain
unexplained. In order to create a “theory of everything”, one would not only
need to combine the SM with general relativity, but also provide an
explanation for many other issues, including the existence of neutrino masses,
the origin of the matter-antimatter asymmetry, and most importantly, the
origin of dark matter. To solve these problems, researchers are collaborating
to formalise new theories, design, build and carry out new experiments, as
well as simulate and analyze research data.

One of the most renowned experimental facilities, the Large Hadron Collider
(LHC) was constructed with the purpose of testing the SM in the high energy
regime. The last elementary particle predicted by the SM, the Higgs boson, was
discovered in 2012 Aad_2012 ; Chatrchyan_2012 . Since then, LHC research has
shifted towards precision measurements and searches for beyond the Standard
Model (BSM) effects.

Many extensions of the SM imply the existence of as yet undiscovered massive
particles, often associated with proposed new symmetry groups. If a new
particle has a narrow decay width, the straightforward method is to search for
a resonant peak in the spectrum of a mass-like observable, such as the
invariant mass of a dijet event. However, such a bump hunt is not completely
free of assumptions. Often complex analytical functions need to be chosen to
model the background distribution, with the possibility to introduce spurious
signals and varying sensitivity under the assumption of different functional
forms. Furthermore, additional observables or fiducial cuts need to be chosen
and optimised to enhance sensitivity in the case where potential signal yields
are low, causing searches to become more model-specific.

Over the past decade, machine learning-based algorithms have become
increasingly popular for solving a multitude of problems. Deep learning, in
particular, has gained popularity for various tasks, with large neural
networks being utilised. For example, many methods were implemented to perform
anomaly detection (AD) tasks in various industries. Some of these AD methods
have been repurposed and extended to support BSM searches Collins:2019jip ;
Collins:2018epr ; DAgnolo:2018cun ; DAgnolo:2019vbw ; Farina:2018fyg ;
Heimel:2018mkt ; Roy:2019jae ; Cerri:2018anq ; Blance:2019ibf ; Hajer:2018kqm
; DeSimone:2018efk ; Mullin:2019mmh ; Dillon:2019cqt ; Aguilar-
Saavedra:2017rzt ; Romao:2019dvs ; Romao:2020ojy ; Amram:2020ykb ;
Cheng:2020dal ; Khosa:2020qrz ; Thaprasop:2020mzp ; Alexander:2020mbx ;
Mikuni:2020qds ; vanBeekveld:2020txa ; Park:2020pak ; Faroughy:2020gas ;
massiveissue ; Chakravarti:2021svb ; Batson:2021agz ; Blance:2021gcs ;
Bortolato:2021zic ; Collins:2021nxn ; Dillon:2021nxw ; Finke:2021sdf ;
Shih:2021kbt ; Atkinson:2021nlt ; Kahn:2021drv ; Dorigo:2021iyy ;
Caron:2021wmq ; Govorkova:2021hqu ; Kasieczka:2021tew ; Volkovich:2021txe ;
Govorkova:2021utb ; Ostdiek:2021bem ; Fraser:2021lxm ; Jawahar:2021vyu ;
Herrero-Garcia:2021goa ; Aguilar-Saavedra:2021utu ; Tombs:2021wae ;
Lester:2021aks ; Mikuni:2021nwn ; Chekanov:2021pus ; dAgnolo:2021aun ;
Canelli:2021aps ; Ngairangbam:2021yma ; Bradshaw:2022qev ; Aguilar-
Saavedra:2022ejy ; Buss:2022lxw ; Alvi:2022fkk ; Dillon:2022tmm ;
Birman:2022xzu ; Letizia:2022xbe ; Fanelli:2022xwl ; Finke:2022lsu ;
Verheyen:2022tov ; Dillon:2022mkq ; Caron:2022wrw ; Park:2022zov ;
Kamenik:2022qxs ; Kasieczka:2022naq ; Araz:2022zxk ; Schuhmacher:2023pro ;
Roche:2023int ; Vaslin:2023lig ; ATLAS:2023azi ; Chekanov:2023uot ;
CMSECAL:2023fvz ; Bickendorf:2023nej ; Freytsis:2023cjr ; Metodiev:2023izu
(see Refs. karagiorgi2022machine ; LHCO ; aarrestad2022darkMachines ;
Belis:2023mqs  for a comparison of various ML assisted BSM methods and Refs.
collins2021compare ; Golling:2023yjq  for a comparison of weakly supervised
and unsupervised approaches). The ATLAS collaboration produced the first
experimental results for such searches applied to experimental data using
weakly supervised methods atlas2020cwolaExp  and unsupervised ML anomaly
detection methods ATLAS:2023azi ; atlas2023twoBody . However, these efforts
have not observed any significant deviations from the SM expectation.

Many AD approaches rely on the assumption that any new signal would form a set
of outliers. However, in a bump hunt the assumption is instead that any new
signal would be localised in some feature space, in particular in an invariant
mass spectrum. Weakly supervised approaches, on the other hand, aim to enhance
the sensitivity by applying a cut on a classifier trained directly on the
data. However, in both instances the same bump hunt restrictions apply with
either functional forms or input observables impacting the sensitivity to a
model.

In this work we introduce a new data-driven method, cluster scanning (CS),
which builds on the foundations of the bump hunt but addresses several
limitations. By leveraging more information from the event CS is able to
enhance sensitivity to potential signals without enforcing any model specific
assumptions, and can also provide a direct estimate of the background
distribution. The proposed approach complements existing techniques and is
designed to be computationally efficient.

The paper is structured as follows. In Section 2, we briefly describe the LHCO
R&D dataset LHCOlympics , commonly used to benchmark the performance of
anomaly detection techniques, and introduce our data preprocessing steps.
Section 3 touches on the general topic of bump-hunting strategies in the
literature, introduces the novel CS method, and discusses similarities and
differences between them. In Section 4 we provide the results of applying CS
in an anomaly search. Finally, we draw conclusions in Section 5.

##  2 Dataset

The LHCO R&D dataset consists of one million background Standard Model dijet
events (also subsequently referred to as QCD) and 100 000 signal BSM
$Z^{\prime}\rightarrow XY$ events, where massive particles with $m_{X}=500$
GeV and $m_{Y}=100$ GeV decay into quark-antiquark pairs. The resonance itself
has a mass of $m_{Z^{\prime}}=3.5$ GeV. This anomaly model is discussed in
detail in Ref. Kim_2020 .

All the events were produced using Pythia 8.219 Sjostrand:2014zea  and Delphes
3.4.1 de_Favereau_2014 ; Mertens_2015 ; Selvaggi_2014  using default settings.
The jets were clustered using an anti-$k_{T}$ algorithm Cacciari_2008  with
$R=1$ using FastJet Cacciari_2012  with a python interface provided through
the pyjet library in Scikit-HEP Rodrigues:2020syo . Jets are required to have
$p_{T}>1.2$ TeV and fall within $|\eta|<2.5$.

###  2.1 Jet images

In addition to the di-jet invariant mass ($m_{jj}$) of the event, used in a
bump hunt, we extract additional information from the image representations of
the two jets. This allows for a more model agnostic approach than selecting
specific jet substructure observables. The jet images are processed following
a prescription similar to that used in Ref. Jet_images ; Kasieczka_2019 ;
Macaluso_2018 ; Heimel:2018mkt  from the $\eta$, $\phi$ and $p_{\mathrm{T}}$
of the jet constituents. Individual jet images are centred, rotated, and
flipped in order to provide a consistent input to a convolutional neural
network, reducing the number of symmetries the ML method would need to learn.

The jet images are cropped to $[-0.8,0.8]\times[-0.8,0.8]$ in $\eta-\phi$
space relative to the jet centre, binned with a $40\times 40$ pixel grid, and
normalised such that the sum of all pixels is equal to one. Fig. 1 shows the
average jet images for QCD background, and the separate averages of all
lighter (mostly $Y$) and heavier (mostly $X$) jets in each $Z^{\prime}$ event.

Figure 1: From left to right: average of all 2M available QCD jet images,
average image of all 100K lighter jets in a $Z^{\prime}$ event and average
image of all 100K heavier jets in a $Z^{\prime}$ event before smearing and
pixel scaling.

Despite being used in many applications, the jet image representation has two
main drawbacks, namely the sparsity of non-zero pixels (see app. B) and the
imbalance in the magnitudes of their intensities. This is particularly
problematic for approaches that depend on the $L_{2}$ (Euclidean) distance. We
address both of these problems with the solutions introduced in Refs.
2021AEanomalies ; Buss:2022lxw .

To take the soft constituents into account, which have intensities orders of
magnitudes lower than hard constituents, we apply a non-linear scaling to all
pixels of $I_{ij}\rightarrow I_{ij}^{\gamma}$. To address sparsity we convolve
(smear) the whole image with a two-dimensional Gaussian kernel with an
isotropic standard deviation $\sigma_{k}$. We find that using a value of
$\gamma=0.5$ for the pixel scaling alongside $\sigma_{k}=1$ for the Gaussian
kernel provides a adequate solution to both issues without excessive impact on
the structure of the jets.

##  3 Method

###  3.1 Bump hunt

The bump hunt approach is a standard method used to search for excesses over a
non-resonant background in HEP (high-energy physics) data. This method usually
follows four main steps that we briefly discuss below. Each of these steps is
a complex topic in itself with several different approaches in the literature,
thus for our study, we choose only simplified and basic approaches.

####  3.1.1 Signal enrichment

Signal enrichment, in general, refers to the selection of a subset of
experimental data in such a manner that the fraction of signal events in it is
increased compared to the initial sample. Most often, this is done by cutting
out a region of the observable space where the signal is expected to be
abundant compared to the rest of the space, typically using a theoretical
model of the signal of interest.

These approaches, despite being sensitive to specific signal processes, make
the search less model-agnostic and are ill-suited for general anomaly
detection searches. Alternatively, one can hope to define a signal-rich region
of the experimental data using a plethora of unsupervised ML (machine
learning) techniques, which are expected to provide enhanced sensitivity over
a wider range of potential signal processes.

In our particular example of LHCO data, we choose to explore a wide, smoothly
falling region of the spectrum of dijet events with invariant dijet mass
$m_{jj}$ from 3000 GeV to 4600 GeV. We choose this lower bound to avoid the
turn on curve of the mass distribution, resulting from the jet trigger, and
the upper bound is selected to remain in a region with relatively high
statistics, so that we work in the region where the fit functions from
Subsection 3.1.2 are applicable. This interval contains, in total, around
380,000 QCD events and nearly all $Z^{\prime}$ events. We divide this region
into 16 non-intersecting bins with 100 GeV width each, as in Ref. cathode ;
curtains .

####  3.1.2 Background estimation

To perform a hypothesis test, one must first postulate a null hypothesis,
which in counting experiments takes form of the expected background coming
from the Standard Model processes. Often the background prediction relies on a
theoretical basis to calculate the cross sections of the hard process and a
simulation to account for detector response and measurement uncertainties.
Still there are a number of searches where theory and simulation cannot
provide a reliable background estimate. In these cases the background has to
be estimated from the data itself in an empiric manner, using some general
assumptions.

In dijet-like searches a background is often estimated by fitting a function
of the form

$f(x)=p_{1}(1-x)^{p_{2}}x^{p_{3}+p_{4}\ln(x)+p_{5}\ln(x^{2})}$ (1)

to a smoothly falling part of the dijet mass distribution CMS:2022eud ;
4par_ATLAS_09_04_2015 ; 3par_ATLAS_03_03_2016 ; 4par_CMS_2017 ;
Khachatryan_2016 ; Khachatryan_2016_2 ; 4par_ATLAS_20_04_2011 ; 4par_CMS_2015
; PhysRevD.87.114015 ; 4par_CMS_2013 ; 4par_ATLAS_2012 ; 4par_ATLAS_2013 ;
4par_CMS_2011 ; 4par_CDF_2009 ; Aad_2011 ; PhysRevLett.105.211801 ;
5-parameter , where $x=m_{jj}/\sqrt{s}$. This function is referred to as the
“$n$-parameter dijet fit function”, where $n$ is the number of nonzero free
parameters $p_{i}$ used in the function. Despite being a good fit to the
simulated data, this functional form is still an empirical assumption and thus
is subjects to a systematic error. Furthermore, after applying some selection
criteria on the events which could be correlated with $m_{jj}$, this function
may no longer well describe the resulting distribution.

More advanced methods of fitting, such as the Sliding Window Fit (SWIFT) SWIFT
and the ABCD method used in Alison:2229034 ; PhysRevD.103.035021  are other
methods that reduce the assumption of a functional form but introduce their
own assumptions instead. However, due to the simplicity and wide use of the
$n$-parameter fit function, we choose to use global 3-parameter and
4-parameter function fits as the benchmark analysis strategy. Further details
of the (pseudo-)analysis on the LHCO R&D data performed using these background
estimates are given in Appendix A. To access the upper bound on the
performance of all background estimation methods, we use the underlying
background distribution as an idealised fit, i.e. a fit with no systematic
error. The (pseudo-)analysis using this is also described in Appendix A.

####  3.1.3 Test Statistic definition and calibration

There are several ways to calculate a global test statistic for two spectra.
In HEP one of the more popular tests in model agnostic searches, called
BumpHunter vaslin2023pybumphunter , relies on the maximal local significance
(MLS) as the test statistic, where it is computed using a range of different
windows over the spectrum. One of the benefits of the MLS test statistic is
its simplicity and that it is well suited for signals that give rise to
narrow, localised resonances. Here the MLS is applied to the binned $m_{jj}$
distributions of the data. Given a set $\mathcal{B}=\\{b_{1},...,b_{n_{\rm
bins}}\\}$ of non-intersecting bins with $N_{\text{sig},b}$ events or jets
from the signal-rich (experimental) distribution and $N_{\text{bkg},b}$ events
or jets from the background estimation, the MLS can be written as

$\begin{split}\text{MLS}=\underset{b\in\mathcal{B}}{\text{max}}Z_{b}=\underset{b\in\mathcal{B}}{\text{max}}(\text{CDF}^{-1}_{\mathcal{N}(0,1)}(\text{CDF}_{Poisson(N_{\text{bkg},b})}(N_{\text{sig},b})))\,,\end{split}$
(2)

where CDF is the cumulative density function of the respective distribution.
In equation 2 only overdensities are taken into account, i.e. $Z_{b}>0$ only
for $N_{\text{sig},b}>N_{\text{bkg},b}$ as we are searching for a resonance.

For bins with $N_{\text{sig}/\text{bkg},b}\gg 1$ one can approximate the
Poisson distribution with a normal distribution
$\mathcal{N}(N_{\text{sig}/\text{bkg},b},\sqrt{N_{\text{sig}/\text{bkg},b}})$.
Equation 2 then reduces to a much simpler form

$\begin{split}\text{MLS}=\underset{b\in\mathcal{B}}{\text{max}}Z_{b}=\underset{b\in\mathcal{B}}{\text{max}}\frac{N_{\text{bkg},b}-N_{\text{sig},b}}{\sqrt{N_{\text{bkg},b}}}\,.\end{split}$
(3)

Although some test statistics, like $\chi^{2}$, have well-known distributions,
other more unusual test statistics, like the BumpHunter test statistic,
require calibration. This is commonly done by modelling its distribution using
Monte Carlo simulation.

Moreover, as systematic uncertainties arise from the definition of a signal
region selection and the background estimate, this calibration should be
performed even in the case where the distribution is known a priori. The
calibration for the BumpHunter test statistic is performed in Ref.
vaslin2023pybumphunter  by running pseudo-experiments in which the counts in
each bin are varied according to Poisson’s law. This can be extended to higher
dimensions by resampling the background events with bootstrapping. By
calculating the test statistic for each of our bootstrapped background-only
pseudo-experiments, we obtain the distribution of the test statistic in the
background-only hypothesis. To ensure good modelling of the tail of the test
statistic distribution, which corresponds to large significance values in the
presence of signal, a large number of pseudo-experiments is required.

####  3.1.4 Significance evaluation

To obtain a calibrated $p$-value for a given value of the test statistic $t$,
one counts the number of background only pseudo-experiments exceeding this
value $N_{>t}$ and divides it by the total number of pseudo-experiments done,
$N_{tot}$.

The (one-sided) significance is computed using the inverse cumulative density
function of the normal distribution
$Z=\mathrm{CDF}^{-1}_{\mathcal{N}(0,1)}(1-p\text{-value})$.

In the case of $N_{>t}=0$ arising from the limited number of pseudo-
experiments, we instead set a lower bound:

$p<\frac{1}{N_{tot}},\;Z>\mathrm{CDF}^{-1}_{\mathcal{N}(0,1)}(1-\frac{1}{N_{tot}})\,.$
(4)

For every experiment with added signal events, we still bootstrap the
background (for consistency) and combine it with a given number of signal
events chosen at random from 100,000 signal events (around 5% of events fall
outside of the evaluation region). Due to statistical fluctuations we also
perform several pseudo-experiments in the signal enriched case in order to
obtain a robust estimate of the significance for each level of signal doping.

###  3.2 Cluster scanning

In this section we present a novel approach called Cluster Scanning, which
follows the same bump hunting scheme, but relies on a distinct set of
assumptions than the commonly employed methods and thus has several favorable
characteristics. Our approach can be divided into several key steps given
below, with the hyperparameters chosen in order to search for narrow
resonances in the $m_{jj}$ spectrum of the LHCO R&D data. The motivation for
these hyperparameters in each step and the argumentation on how to choose them
for a different application case is given in App. C.

##### Training region selection:

We select a narrow $m_{jj}$ window [3000, 3100] GeV for training of the
k-means algorithm. This window contains 56,486 original background events. In
this publication, we focus on relatively small signal injections that include
only 5% or less of the total number of Z’ signals available. Therefore the
training region is expected to contain 89 signal events or less, which can be
regarded as negligible. Despite observing a qualitative improvement in
performance in case the training region matches the resonant peak and thus has
a larger portion of signals events involved in clustering, in an actual
analysis the position of the peak will be unknown, thus we choose to discuss a
more representative case given here, when the training region happens to be in
the tail of the signal peak and thus has a negligible number of signal events.

##### K-means Clustering:

We apply a mini-batch $k$-means clustering algorithm with $k=50$ implemented
in the scikit-learn scikit-learn  library, with a batch size of 2048 on the
set containing jet images of the leading two jets from each event in this
$m_{jj}$ window. The mini-batch implementation is chosen due to its
computational speed. The seeding of the cluster centroids is performed using
the k-means++ prescription described and motivated in Ref. kmeans .

##### Cluster Spectra:

After performing the fit of $k$-centroids to the data in the training region,
we fix the centroid positions and evaluate how many jet images from each of
the 16 $m_{jj}$ bins of the evaluation region, defined in Subsection 3.1.1,
fall into each of the $k$ clusters $N_{i,b}$ where
$i\in\\{1...k\\},\,b\in\\{1...n_{bins}\\}$. Fig. 2 shows the resulting 50
normalised cluster spectra $N_{i,b}/\sum_{b}(N_{i,b})$ for one pseudo-
experiment with signal injection.

Figure 2: The $m_{jj}$ distributions for the jets in each of the 50 clusters,
each normalised to unity. Here, 5,000 signal events have been injected into
the evaluation dataset, which corresponds to 5% of the total available signal
events.

##### Per bin standardisation:

We note that in each bin the normalised cluster spectra follow an
approximately normal distribution with several outliers from the anomalous
clusters (see discussion in App. D). Therefore we standardise the normalised
cluster spectra in each bin using outlier robust estimators (described in App.
E) for mean and standard deviation with an outlier factor of 0.2. Here we make
the assumption that the majority of the signal is located in a small number of
clusters, and the rest of the clusters are signal depleted. Figure 3 shows the
cluster spectra from Fig. 2 after normalising with the outlier robust
estimator.

Figure 3: Spectra in Fig. 2 standardised over clusters in each bin.
Potentially signal-rich cluster spectra are shown in red.

##### Selecting anomalous clusters:

Utilising the assumption that the signal is localised in $m_{jj}$, we select
potentially signal-rich cluster spectra as those with a deviation of more than
a threshold value of $\theta=3$ standard deviations from the robust mean in
the positive direction as we are only interested in a resonance leading to
excess of events. The rest of the clusters are labelled as signal-depleted.
The threshold and the selected signal-rich clusters are shown in Fig. 3 in
red.

##### Signal-rich and signal-depleted regions:

After the selection, we combine the non-normalised distributions corresponding
to our selected signal-rich clusters. This results in a signal-rich spectrum
$N_{\text{sig},b}$ with an example shown in red in Fig. 4.

The remaining cluster spectra are combined to form a signal-depleted spectrum
$N_{\text{poor},b}$. The estimate of the background is then constructed by
normalising it to the same total entries as in signal-rich spectrum, namely
$N_{\text{bkg},b}=N_{\text{poor},b}\frac{\sum_{b}N_{\text{sig},b}}{\sum_{b}N_{\text{poor},b}}$.
It is shown in blue in Fig. 4.

Figure 4: Curves corresponding to the sum of signal-rich and signal-poor
spectra in Fig. 3. The blue signal-poor curve is rescaled to have the same
total jet number as the signal-rich curve. The coloured region around the blue
curve is a $\sigma_{\text{bkg},b}=\sqrt{N_{\text{bkg},b}}$ Poisson deviation
after recaling used to compute MLS.

##### Test statistic:

As previously discussed, to test the significance of an observed excess we use
the simple maximum local significance, as defined in Equation 2. It may occur
that no cluster is selected as anomalous. In this case we assign a default
value of 0, in order to show good agreement with the null-hypothesis
expectation. This is similar in motivation to setting the value for an
observed deficit in events to zero. Following the discussion in Subsection 3.1
for the calibration process, we construct 3,900 pseudo-experiments using
bootstrap resampling on 1 million background events. The distribution of the
test statistic is discussed in App. F.

##### Ensembling:

Different initialisations lead to a broader distribution over the final test
statistic obtained with cluster scanning. In order to obtain a final value for
the test statistic, the cluster scanning method is performed 15 times with
independent initialisations. The mean of the test statistic from all the runs
forms the final ensembled test statistic. The distribution of this statistic
is presented in app. F.

###  3.3 Discussion

As we can see, CS follows the general bump hunt strategy, but introduces novel
approaches for the first two steps of this strategy. First of all, CS selects
the most anomalous looking clusters to define the signal-enriched region, and
constructs a background estimate from the rest of the clusters. Notably
though, this selection is completely data-driven and does not target a
specific family of signal models. However, CS relies on a set of assumptions
that fundamentally differ from those commonly used in other anomaly detection
approaches.

##### Search for overdensity instead of outliers:

Most anomaly search methods like Autoencoders 2021AEanomalies  and SVDDs
aarrestad2022darkMachines  rely on outlier detection, namely, identifying the
data instances that lie in a region of very low probability density or outside
the support of the “normal” distribution. Notably, while all normal events
share similar characteristics and exhibit easily recognisable trends,
anomalous data, such as defects or fraud, can differ in numerous ways and are
thus given a wide prior. Although model-agnostic searches should accommodate a
wide range of possible anomaly models, it is usually assumed that a signal is
produced by only one or a few unknown BSM process. Thus, all anomalous events
have many features in common and exhibit some similarity to SM events, as any
new particle must radiate and decay into SM particles to be detectable.

Therefore we use the localisation of anomalies in both low-level (e.g. jet
images) and high-level variable (e.g. $m_{jj}$) space as the first main
assumption of the CS method. Localisation of anomalies in low-level variable
space means that only a few out of all clusters contain a fraction of
anomalies much higher than the rest of the clusters. This way clustering plays
a role of data-driven binning in low-level variable space. Localisation in
$m_{jj}$ gives us a possibility to distinguish these anomaly rich clusters
from the rest, namely, by searching for an overdensity in $m_{jj}$ in one
cluster spectrum compared to all others. Thus, CS is able to select a signal-
rich region of events by leveraging the assumption of signal being localised
rather than consisting of outliers.

Although semi-supervised methods based on CWoLa (see Refs. salad ; cathode ;
curtains ; feta ; curtainsf4f ; buhmann2023phase ; Sengupta:2023vtm ) and
density estimation methods are also sensitive to overdensities, they usually
require construction of a background template, which until recent developments
buhmann2023phase ; Sengupta:2023vtm  was preferably constructed for a smooth
distribution of low dimensionality, typically using a few high-level
observables. In this publication, we show that CS is able to draw significant
improvement from a high-dimensional distribution of low-level jet observables.
In this way, it can be considered less signal-model dependent than the methods
that rely on hand-crafted high-level observables.

##### Assume cluster mass independence instead of smoothness:

CS proposes a solution to the second step of the analysis, namely, it
estimates the form of the background by combining the signal-depleted
clusters. In this way, we do not rely on any assumptions on smoothness or on a
particular functional form of the background-only spectrum in $m_{jj}$, which
are heavily relied upon by most other bump hunt methods, such as the global
functional fit mentioned in Subsection 3.1, SWIFT SWIFT  and even Gaussian
processes frate2017modeling .

Instead, the second main assumption of CS is that in the background-only case,
the assigned cluster centroid is approximately independent from $m_{jj}$.
Ideally, we would want the distribution of background events over $m_{jj}$ in
each cluster to be identical within statistical uncertainty, such that the
probability of a jet belonging to a cluster and having a specific mass
factorises, $p(i,m_{jj})=p(i)p(m_{jj})$, or at least that the correlation is
weak. This would minimise the rate of incorrectly identified signal-enriched
clusters. In practice, although Fig. 2 shows that the distributions all follow
a similar trend, there are still some systematic deviations. These are a
result of the finite width in $m_{jj}$ of the training window and slight
correlations between the distribution of the jet constituents and $m_{jj}$
arising from the transverse momenta of the two non-resonant jets depending on
$m_{jj}$. Therefore, for the selection of the clusters, we estimate the full
uncertainty (including Poisson fluctuations and uncertainty from mass
dependence) separately for each experiment and for bin based on the sample of
our $k$ cluster spectra values.

Unlike in methods with sliding window approach SWIFT ; curtains , in CS the
fit only needs to be performed once. Moreover $k$-means clustering is a simple
classical algorithm that typically requires less training than deep learning
approaches, making CS a relatively fast analysis method. This is important in
the context of the ensembling and calibration, which both require a large
number of analysis iterations, and are thus a notable obstacle to
incorporating deep learning in HEP analysis under the constraint in computing
resources. Fast analysis is also advantageous for testing its efficiency for
simulated BSM events in order to produce the exclusion limits (see the RECAST
Cranmer_2011  framework). Moreover, CS avoids other disadvantages inherent to
sliding window approaches, such as limited search range due to the definition
of the sidebands and the need to optimise sideband and signal window widths.

###  3.4 Idealised CS

Despite choosing a narrow $m_{jj}$ window to reduce mass dependence
systematics, the variables that we use for clustering are in general not
independent of $m_{jj}$. Thus, we observe the background-only spectra of some
clusters do not just statistically fluctuate around the expected shape of the
background, but exhibit some degree of smooth mass sculpting. This affects the
performance of the method by introducing false positives at the cluster
selection stage and increasing the discrepancy between signal and background
spectra when evaluating the test statistic. This may be partially remedied by
a more sophisticated method of selecting anomalous clusters or a better
background estimate, both of which rely on more assumptions. These studies are
outside the scope of this publication. However, to give an upper bound on the
performance one may achieve with such improvements we propose an idealised
version of clusters scanning.

Idealised CS version requires us to modify the distribution of the jets
between the clusters. First, we count the numbers of jets that fall into each
cluster in the first $m_{jj}$ bin. If no mass dependence were present, the
fractions of QCD jets in each cluster
$x_{\text{QCD},i,b}=N_{\text{QCD},i,b}/\sum_{i}N_{\text{QCD},i,b}$ should be
independent of bin number $b$ within statistical uncertainties. To simulate
this case in all the consecutive bins except the first we distribute the QCD
jets in these bins among clusters using a multinomial distribution with
weights equal to the fractions obtained in the first bin
$x_{\text{QCD},i,b}=x_{\text{QCD},i,1}$, thus generating cluster spectra that
follow the original background spectrum with statistical fluctuations, i.e.
the case with no mass dependence. The signal jets are distributed as before
according to which cluster they belong to, such that the fractions of
$Z^{\prime}$ jets may differ between different bins. This is done because we
assume that only the background is distributed roughly proportionally between
clusters, which is equivalent to assumption 2, but not the signal.

This distribution of jets creates idealised cluster spectra for each
clustering, and the rest of the algorithm remains unchanged.

##  4 Results

As a proof of concept we perform an analysis applying CS and global fit based
bump-hunting with the above mentioned hyperparameters to the LHCO R&D dataset
with different amounts of signal injection, given in figures either as an
absolute number of injected events $\epsilon$ or as a signal to background
ratio $S/B$ of events in the considered [3000, 4600] Gev $m_{jj}$ region. For
each pseudo-experiment with signal injection we calculate the significance $Z$
as discussed in Subsection 3.1.4 using the calibration test statistic
distribution. For each signal injection level we run 100 pseudo-experiments
with bootstrapped background data and randomly sampled signal events. As a
reference for the significance and its statistical variation for each
contamination level, we report the median significance of these pseudo-
experiments and 0.25 and 0.75 quantiles. We define the ratio between the
median significance provided by CS to the significance of a baseline method as
the significance improvement (SI). We also quote the relation between the
number of events needed to obtain a $3\sigma$ evidence in each analysis
strategy.

Figure 5(a) shows how the global significance of CS and the parametric fit-
based methods depends on the signal contamination in the pseudo-experiment. It
characterises the performance of these realistic analysis strategies, which do
not use any truth information for the evaluation of the test statistic, thus
including all the systematical uncertainties coming from partially fulfilled
assumptions needed for the respective method.

(a)

(b)

Figure 5: Upper figure: Median and the quartile bounds of global significance
of the signal contaminated pseudo-experiments as a function of the number of
signal events $\epsilon$ injected shown for the (a) realistic and (b)
idealised analysis methods. The dotted lines mark lower bounds, as there was
not enough statistics to access higher significance levels. Bottom figure:
Significance improvement of the CS method compared to the 4-parameter fit for
the (a) and the significance improvement of idealised CS compared to idealised
fit for the (b) as a function of the signal-to-background ratio $S/B$.

We observe that although 3- and 4-parameter fits give approximately the same
results, CS outperforms them by a significant margin in the region from 1500
to 4000 signal events. Beyond 3000 signal events, the significance yield from
CS is limited by the number of bootstrap pseudo-experiments in the calibration
set, but the lower bound on its significance still remains substantially
higher than the significance of the parametric fits. This is the most
interesting region as there the transition between non-significant signal
(below $1\sigma$) and new physics evidence (above $3\sigma$) takes place.
Looking at the lower subplot in Fig. 5(a) we see that CS gives us an SI of 2
and higher on the majority of regions of interest. We can also see that CS
produces a 3$\sigma$ evidence for only 61% of the events needed to obtain this
evidence with the parametric fit. This shows that although both suffer from
fit and assumption induced systematic uncertainties, CS has a clear advantage
over parametric fitting procedures.

Above we have also described the idealised version of the cluster scanning
method and in Appendix A the analysis with an idealised background fit. Both
methods rely on event labels to remove systematic uncertainties introduced by
the limitations of the assumptions of our methods and to make the background
estimate in both cases close to the true background, with only statistical
fluctuation taken in consideration. This is done to separate the influence of
additional information, namely the low-level observables used in the analysis,
from the systematic uncertainties introduced by the fits, and to construct the
upper bound on the performance of our methods.

Figure 5(b) shows how the global significance of both idealised methods depend
on the signal contamination in the pseudo-experiment. It can be seen that one
needs substantially less signal for it to be significant in the idealised
methods compared to the realistic methods, as it removes false-positives
induced by systematic uncertainties in the fits. Still, we see that the
idealised CS outperforms the idealised functional fit on the majority of the
interval between 1000 and 2200 signal events. Looking at the lower subplot in
Fig. 5(b) we see that by using CS in the region of interest we gain a
significance improvement factor of up to 1.5. We can also see that CS produces
a $3\sigma$ evidence with only 69% of the events needed to gain this evidence
with the idealised fit. This shows that in the case of negligible systematic
uncertainties, CS gives an improvement over any smooth fit as, in addition to
just using information from $m_{jj}$, it also makes use of the low-level event
information. From the difference between idealised and non-idealised CS we can
see that there is some room for improvement of CS to reduce the false positive
rate, and improve the analysis efficiency.

##  5 Conclusions and outlook

This paper is a first proof of concept for the cluster scanning anomaly search
method, which is designed to search for resonant overdensities on the
distribution of an observable using clustering techniques in auxiliary
observables.

We found that it outperforms the widely used bump-hunting method, which relies
on the functional background fits, in several metrics relevant to an analysis.
In the transition region, where the benchmark algorithm achieves 1$\sigma$ to
3$\sigma$ significance, CS improves the result by a factor of 2 or more for
the realistic case, or by a factor of 1.5 for the idealised comparison. This
reduces the number of signal events required to produce a $3\sigma$
significance by a factor of 0.61 in the realistic case and by a factor of 0.69
in the idealised case. The former factor of improvement should be expected in
a real application. We also discuss the comparison of cluster scanning with
other anomaly detection algorithms in Subsection 3.3, outlining its advantages
and limitations.

The CS method should not be seen as a direct competitor to background fitting
methods, but rather as a complementary approach that relies on a different set
of assumptions about the nature of the anomaly and the background
distributions, which are not well known.

There remains a large unexplored field of potential extensions and
improvements to this method or synergies with other methods. Straightforward
follow-up studies can explore the use of clustering methods other than
$k$-means. One can look for other ways of selecting the anomalous clusters,
alternatives to the one proposed in Subsection 3.2, that would rely on
different assumptions. For example, one can require that all anomalous
clusters are neighbors in the space of clustered inputs. One can also unify
the assumptions of CS and functional fits to produce separate background
estimates for each of the clusters separately, greatly reducing the $m_{jj}$
dependent systematic uncertainties.

CS could benefit from using features developed by other algorithms that have
already been optimised for other tasks, such as flavour tagging, or even using
unsupervised learning for feature extraction.

Furthermore, since many other ML approaches to improve sensitivity in model-
independent searches rely on a bump hunt for the final statistical analysis,
CS could also be used to further enhance sensitivity. This could be of
particular interest when the background distribution is no longer well
described by simple, smoothly decreasing functional forms.

###### Acknowledgements.

The authors would like to acknowledge funding through the SNSF Sinergia grant
CRSII5_193716 “Robust Deep Density Models for High-Energy Particle Physics and
Solar Flare Analysis (RODEM)” and the SNSF project grant 200020_212127 “At the
two upgrade frontiers: machine learning and the ITk Pixel detector”. The
research of MK is supported by the DFG under grant 396021762 - TRR 257:
Particle Physics Phenomenology after the Higgs Discovery.

## Code availability

The code used to produce all results presented in this paper is available at
https://github.com/IvanOleksiyuk/jet_cluster_scanning .

## Appendix

##  Appendix A Idealised fit and n-parameter fit pseudo-analysis

In general, if a distribution $\mathcal{H}_{b}$ of a background events is
perfectly known a-priory, given a binning for this distribution, one can
calculate the expected number of events in each bin. Being provided directly
from the true underlying $\mathcal{H}_{b}$, this background estimate will on
average provide the most efficient tests to discriminate samples drawn from
$\mathcal{H}_{b}$ from samples drawn from an alternative hypothesis
$\mathcal{H}_{b+s}$ with signal, compared to any other estimate of the
expected background for these samples. Hence we call the expectation from
$\mathcal{H}_{b}$ an ”idealised fit”.

As discussed extensively in section 3.1, estimation of the expected background
is only one step of the analysis. To create a benchmark, we do pseudo-analysis
on LHCO R&D dataset represented by spectrum $N_{\text{bkg},\text{orig},b}$
using the other choices defined in section 3.1. Namely, we generate pseudo-
experiments by bootstrap resampling the events from
$N_{\text{bkg},\text{orig},b}$ and add a number of signal events if needed.
The ”idealised” background estimation for every pseudo-experiment is equal to
$N_{\text{bkg},\text{orig},b}$ itself (as the samples were generated with
these expected values). Following the discussion in subsection 3.1.3 we use
MLS test statistic between this estimate and the generated pseudo-experiments,
to generate null-hypothesis test statistic distribution and its value for
signal contaminations and thereafter estimate the significance. Depending on
the number of doped signal events, the median and quartile region significance
given by this test is provided in the main text in Fig. 5(b).

Unfortunately the background model is usually unknown, so for each
experimental sample the background should be estimated in some less precise
way relying on weaker assumptions.

As a realistic benchmark to our method we explore how sensitive the analysis
is using global n-parameter functional to the kind of signal presented in LHCO
R&D dataset. We use the binning with 16 bins defined in subsection 3.1.2 and
count the number of background events in each bin to get an original
background spectrum $N_{\text{bkg},\text{orig},b}$.

For all the fits in this studies, we use Trust Region Reflective nonlinear
least squares fitting algorithm implementation from Scipy python package
2020SciPy-NMeth . The chosen bins generally contain more than 5000 counts, so
the Poisson distributions of these counts can be well approximated by a
Gaussian distributions with the variances equal to the bin counts. Using
variances to scale the summands in the least squares objective we make it
equivalent to the maximum likelihood objective for this setup.

First, we fit our 3- and 4-parameter functions to the spectrum to see if the
fit is valid. Resulting fits with 13 and 12 degrees of freedom score
$\frac{\chi^{2}_{3-par}}{n_{dof}}\approx 1.201$ and
$\frac{\chi^{2}_{4-par}}{n_{dof}}\approx 1.338$ that correspond to p-values of
$0.275$ and $0.182$ which signify validity of these fits.

Unlike the CS method that doesn’t generally rely on the smoothness of the
background, global n-parameter takes it as the main assumption, so as
$N_{\text{bkg},\text{orig},b}$ already has some statistical fluctuations a
distribution resampled from it will have even larger statistical fluctuations
than the ones expected for Poisson distribution. To simulate the proper scale
Poisson fluctuations in the chosen region for our pseudo-experiments we
resample events not from $N_{\text{bkg},\text{orig},b}$ but from the best
possible fit. This also negates the systematic error from null-hypothesis not
corresponding to the empirical functional form, so these experiments can be
viewed as semi-idealised. In a more realistic cases, the space of functions
given by all possible parameter values, does not contain the true form of
null-hypothesis distribution and can only yield an approximation of it with
limited precision. It is usual for fit functions with a small number of
parameters, but with increasing number of parameters the function fit problem
becomes over-defined and the function can fit the signal bump as well.
Experimentally we have observed only insignificant increase in performance
when comparing sampling from $N_{\text{bkg},\text{orig},b}$ or from the best
fit distributions. On top of the resampled background events we add a number
of signal events from signal’s original distribution when needed.

The initial parameters of the fit in each experiment are chosen to be equal to
the optimal parameters of the initial fit discussed above, so that one gets an
”idealised” background fit if no optimisation is done. However, because of the
statistical fluctuations and/or added signal contamination, the maximisation
of likelihood results in a different set of parameters for this functional
form. This error of background mismodeling under its statistical fluctuations
and addition of the signal is exactly the type of error we want to demonstrate
with this pseudo-analysis.

The results of such analysis for different signal contamination is given in
Fig. 5(a). We can see that the 3-parameter fit provides a slightly better
result than 4-parameter fit as the latter has more flexibility to overfit the
signal and the statistical fluctuations. This is so because the samples are
drawn from 3- and 4-parameter functions with fixed parameters themselves. If
we were to sample from other distribution the error coming from mismatch in
true end expected functional forms may switch this ordering but it will reduce
both performances. Therefore, the curves shown in Fig. 5(a) are upper limits
of these realistic n-parameter fit analyses achievable only when the true
distribution is described by one of the functions in the chosen parameterised
space.

##  Appendix B Sparsity of the jet images

Fig. 6 show that the jet images are very sparsely populated ususally having
less tha 100 non-zero pixels per 1600 pixels total.

Figure 6: Distribution of images in QCD and top datasets from top-tagging task
vs the number of non-0 pixels in them.

##  Appendix C Hyperparameter selection and motivation

In this appendix we give motivation for every not yet discussed choice of
hyperparameter in our pseudo-analysis. All the hyperparameter suggestions are
done in an unsupervised way coming from general assumptions about signal and
background and are not optimised using the truth information from LHCO R&D
data. As such the levels of significance improvements may be further increased
by performing a dedicated parameter scan for a specific application, however,
we recommend to follow the same reasoning when applying CS in other analyses.

##### Training region:

Training on the full spectrum would likely result in each cluster
corresponding to a specific mass region, thus the background spectrum for each
cluster would not be close to the original mass spectrum. Therefore we perform
clustering in a narrow mass window [3000, 3100] GeV. We choose this window as
it lies in the studied region defined in subsection 3.1.1 and has the largest
statistic of all other 100 GeV windows.

##### Number of clusters:

The most important parameter we had to choose is the number of clusters $k$.
Two factors play the key role in this choice. On one hand, the number of
clusters has to be as large as possible to better narrow down the anomaly-rich
region. On the other hand, for a given number of events in the evaluation
region and the binning of this region one has to take the number of clusters
sufficiently small so that the least populated clusters in the smallest
$m_{jj}$ bin $\underset{i,m_{jj}}{min}(N_{i}(m_{jj}))$ still has enough
statistics for a meaningful statistical analysis. We assume that
$\underset{k,m_{jj}}{min}(N_{i}(m_{jj}))=O(50)$ should be sufficient. Using a
coarse search, we determine, that for our choice of binning and overall
statistic at hand choosing $k=50$ gives a good trade-off as it has a median of
55 events in smallest cluster-bin and it goes below 20 only 1 time in 1000
pseudo-experiment runs, as can be seen in Fig. 7.

Figure 7: Distribution of the number of jet image counts in the least
populated bin of the least populated cluster in each of the 1000 random
background only runs of the CS algorithm on backround only data.

##### Batch size:

Scikit-learn scikit-learn  documentation states that the parallelisation is
performed on all available $N_{cores}$ computing cores if the batch size is
$N_{cores}\cdot 256$ or larger. We performed all computations with 8 core
parallelisation, thus the natural choice of a batch size was 2048. It is also
important to maintain the batch size much larger than the number of clusters
to ensure faster convergence.

##### Outlier fraction:

To quantify the performance of our method we introduce the signal fraction
improvement score (SFI) that characterises a subset $\mathcal{S}$ of the
events in evaluation set $\mathcal{E}$ by the relative increase in the signal
to background ratio

$SFI(\mathcal{S})=\frac{N_{\text{sig}}(\mathcal{S})}{N_{\text{QCD}}(\mathcal{S})}\frac{N_{\text{QCD}}(\mathcal{E})}{N_{\text{sig}}(\mathcal{E})}.$
(5)

Our main assumption is that the signal is distributed in clusters unevenly and
there only several clusters have a significantly large SFI. To put it in
numbers, we assume that not more than 20% of clusters have SFI of 2 or more.
Following this assumption we choose the outlier fraction of 0.2 for outlier
robust estimators. This is an ad. hoc prior assumption about the data at hand,
and it has to be made prior to analysis and has no way to be validated without
knowing the truth lables. Still we can show that this assumption is satisfied
in our case with a margin for the pseudo-experiment shown on all the figures
of section 3. 5000 signal events were giving an overdensity on the original
spectrum that was not identifiable as a deviation from smooth background by
human eye (without knowing the background truth), but in Fig. 2 one can easily
notice two spectra with a significant bump around 3.5 TeV that stand out of
the crowd of other spectra. Unsurprisingly these two spectra have SFIs of 9.1
and 8.9. Three more clusters also have a visible overdensity at this position
possessing SFIs of 6.3, 5.6, 4.4. In total, exactly 8 clusters have $SFI>2$.
Still as we will see later only 3 of these clusters have a signal significant
enough to be selected as anomalous, showing that our assumption is quite
conservative in its limit and either the threshold $SFI$ can be increased or
the percentage of clusters to path the threshold reduced for it to still
remain a valid assumption. Runs of the analysis on other (pseudo-)experiments
behave in the similar manner.

##### Cluster selection threshold $\theta$:

First of all, we use the threshold only for positive deviations as we only
search for excesses of events. Apart from the signal-rich outlier clusters the
threshold can be passed by signal poor clusters, but only with an expected
false positive rate of
$1-(1-\text{p-value}_{\mathcal{N}(0,1)}(\theta))^{n_{bins}}$. Then for large
enough thresholds the average number of false positives can be estimated as
$k\cdot n_{bins}\cdot\text{p-value}_{\mathcal{N}(0,1)}$. Higher thresholds
result in lower false positive and lower true positive rates. To retain the
sensitivity for statistically small signal we choose to use $\theta=3$ that
will result in approximately $50\cdot 16\cdot 0.00135=1.08$ signal poor
cluster being assigned a false positive label on average. Fig. 3 shows 4
clusters being chosen using this threshold. Three of them have an overdensity
at 3.5 TeV and one does not, implying that it is a likely false posive.

##### Ensemble size:

We recommend to take the ensemble size as high as possible, for given
computation resource constrains to reduce the width of the test statistic
distribution (see appendix F).

Parameter | value | motivation
---|---|---
k | 50 |  $\underset{i,b}{min}(N_{i,b})=O(50)$ with binning below
mini-batch | 2048 |  $N_{cores}\cdot 256$, must be $\gg k$
Training region | [3, 3.1] TeV | narrow mass window with high statistic
Evaluation region | [3, 4.6] TeV | the n-parameter fit is applicable
|  | excluding low statistic regions
Bin width | 100 Gev | broad enough to have sufficient statistics in each bin
outlier | 0.2 | consistent with assumption
fraction $f$ |  | on the maximum number of signal clusters
Cluster selection | 3 $\sigma$ | low enough to let trough many true positives
threshold $\theta$ |  | but high enough to filter most false positives
Test statistic | MLS | simple and specialised for local excesses
Default TS | 0 | minimal test statistic possible
Ensemble size | 15 | As large as possible realistic compute limitations
Table 1: Summary of the hyperparameters used in cluster scanning.

##  Appendix D Gaussianity of cluster scanning bin entries

Assumption on the Gaussianity of cluster spectra in each bin can be shown to
be valid by standardising the background-only cluster counts in each bin and
checking if these distributions match $\mathcal{N}(0,1)$. 50 samples are
usually not enough to determine whether the distribution is Gaussian or not,
but by marginalising over 16 bins we get 800 samples in total. Fig. 8 shows
the said distribution for background only spectra that fits the
$\mathcal{N}(0,1)$ distribution well visually and by a consensus of 3
normality tests Shapiro-Wilk, Kolmogorov-Smirnov and Jarque-Bera (note that
the p-values for all the tests are high). The same distribution with added
signal has some clear outliers, which reflects in lower p-values for the
normality tests, but apart from these it still can be well approximated by a
unit Gaussian.

Figure 8: Histogram of deviations from the average of normalised cluster
spectra measured in standard deviations of each bin. Blue histogram gives it
for signal-free case and orange histogram for case with 5000 injected signal
events. Red line show the expected normal distribution.

##  Appendix E Outlier robust estimators

While searching for outliers, it is preferred to use outlier robust estimators
for standard deviation (SD) and mean. We define them as follows: given a
sample of observations $S=\\{x_{1},x_{2},...\,x_{n}\\}$ we find a median
$med(S)$ (which is itself an outlier robust estimator) of this sample and take
a subsample $\tilde{S_{f}}$ that is constructed from $S$ by discarding a
fraction $0<f<1$ of all samples that have largest absolute distance to this
median. In this way we have discarded the outliers. After that we construct
estimators $\tilde{\mu}_{f}=mean(\tilde{S}_{f})$ and
$\tilde{\sigma}_{f}=SD(\tilde{S}_{f})\cdot g(f)$. If $S$ is a sample from
$\mathcal{N}(\mu,\sigma)$ it is obvious that with
$\underset{n\rightarrow\infty}{lim}\tilde{\mu}_{f}=\underset{n\rightarrow\infty}{lim}mean(S)=\mu$.
If one takes $S$ from $\mathcal{N}(0,1)$ and rescales $x_{i}\rightarrow\sigma
x_{i}$, then both estimators transform as
$\tilde{\sigma}_{f}\rightarrow\sigma\tilde{\sigma}_{f}$ and
$SD(S)\rightarrow\sigma SD(S)$ by definition, so both estimators
$\tilde{\sigma}_{f}$ and $SD(S)$ are proportional to a true $\sigma$ of the
Gaussian distribution.

Both $\tilde{\sigma}_{f}$ and $SD(S)$ are independent of $\mu$ and there are
no other parameters of the normal distribution for estimators to depend on,
therefore for a family of Gaussian distribution estimators
$\tilde{\sigma}_{f}$ and $SD(S)$ are proportional to each other by some
constant factor $g(f)$ in the limit of infinite sample. In other words,
adjusting numerically
$g(f)=\frac{SD(\mathcal{N}(0,1))}{\tilde{\sigma}_{f}(\mathcal{N}(0,1))}=\frac{1}{\tilde{\sigma}_{f}(\mathcal{N}(0,1))}$
is sufficient to sattisfy
$\underset{n\rightarrow\infty}{lim}\tilde{\sigma}_{f}=\underset{n\rightarrow\infty}{lim}SD(S)=\sigma$.
So $\tilde{\mu}_{f}$ and $\tilde{\sigma}_{f}$ are unbiased estimators of $\mu$
and $\sigma$ of a normal distribution, although depending on $f$ they are less
efficient than usual non-robust mean and SD.

Fig. 9 shows us the cluster spectra from Fig. 2 with subtracted normalised
original spectrum (which is only needed for better visualisation as this step
has no effect on the standardisation). Fig. 9 also shows the conventional and
the outlier robust estimations of mean and SD of the cluster spectra values in
each bin. As expected for lower $m_{jj}$ the SD is higher as these deviations
is partially caused by the Poisson fluctuations which are proportional to
$\sqrt{N_{i,b}}$. We can also see the conventional estimators have a bump
around 3.5 TeV that is induced by our outlier signal-rich clusters, while the
robust estimators are unaffected by the outliers.

Figure 9:  Normalised spectra with subtracted normalised original $m_{jj}$
spectrum. Amount of signal is 5000. The selection of the anomalous clusters is
taken from Fig. 4.

##  Appendix F Calibration distributions

The distribution of the test statistics given by CS without ensembling for all
background only pseudo-experiments is shown in Fig. 10(a) as a histogram. We
see that around 300 of those were assigned test statistic of 0 as they had no
clusters selected as anomalous. Other cases where one or more anomalous
clusters were selected form a smooth continuous distribution.

The median CS test statistic for 100 signal contaminated pseudo-experiments is
represented in Fig. 10(a) by a vertical line, and the vertical band represent
the region between the quartiles of such a test statistic sample. For each
signal-doped pseudo-experiment we calculate significance as it is described in
subsection 3.1.4. The median significance is quoted in the legend of the
figure.

Fig. 10(b) shows the distribution that is analogous to the one in Fig. 10(a),
but with an ensemble of 15 runs of CS algorithm for each pseudo-experiment. We
notice that the distribution in Fig. 10(b) is significantly narrower than in
10(a) which reduces the frequency of background only experiment having large
test statistic, thus increasing the sensitivity to signal injection. An
additional benefit is that the uncertainty region (between two quarterlies)
for each signal doping have significantly decreased which is important for
lower uncertainty in the analysis on experimental data on the excess
significance or on the exclusion limits.

This motivates, that in general the ensemble size should be taken as large as
reasonably possible. Our choice of ensemble size 15 together with the number
of pseudo-experiments 3900 were dictated by the computing time and storage
memory limits as the amount of full CS algorithm iterations is the product of
those numbers (excluding the pseudo-experiments with signal injection)

(a)

(b)

(c)

Figure 10: Histogram of the CS test statistic for pseudo-experiments with
bootstrapped background only samples. Vertical lines and vertical bands show
median and region between lower and higher quartiles of test statistics for
pseudo-experiments with signal injection. Several signal injection levels are
represented by different colours. Panel (a) shows a case with only 1
initialisation of clusters in CS per pseudo-experiment, panel (b) shows a case
for ensembling 15 runs of CS with different intialisations per pseudo-
experiment and panel (c) shows a case for ensembling 15 runs of idealised CS
with different initialisations per pseudo-experiment.

Finally, Fig. 10(c) shows that for idealised CS without systematics introduced
by mass correlations the MLS between our signal spectrum and background
estimate is lower. Moreover, as expected, it improved the sensitivity of the
method to the signal. Obviously this technique cannot be utilised in an actual
analysis as jet labels are needed to distribute signal and background jets in
a different manner. Moreover, as expected, it improved the sensitivity of the
method to the signal. Obviously this technique cannot be utilised in an actual
analysis as jet labels are needed to distribute signal and background jets in
a different manner.

## References

  * (1) ATLAS Collaboration, _Observation of a new particle in the search for the standard model higgs boson with the ATLAS detector at the LHC_ , _Physics Letters B_ 716 (2012) 1–29 [1207.7214].
  * (2) CMS Collaboration, _Observation of a new boson at a mass of 125 GeV with the CMS experiment at the LHC_ , _Physics Letters B_ 716 (2012) 30–61 [1207.7235].
  * (3) J.H. Collins, K. Howe and B. Nachman, _Extending the search for new resonances with machine learning_ , _Phys. Rev. D_ 99 (2019) 014038 [1902.02634].
  * (4) J.H. Collins, K. Howe and B. Nachman, _Anomaly Detection for Resonant New Physics with Machine Learning_ , _Phys. Rev. Lett._ 121 (2018) 241803 [1805.02664].
  * (5) R.T. D’Agnolo and A. Wulzer, _Learning New Physics from a Machine_ , _Phys. Rev. D_ 99 (2019) 015014 [1806.02350].
  * (6) R.T. D’Agnolo, G. Grosso, M. Pierini, A. Wulzer and M. Zanetti, _Learning multivariate new physics_ , _Eur. Phys. J. C_ 81 (2021) 89 [1912.12155].
  * (7) M. Farina, Y. Nakai and D. Shih, _Searching for New Physics with Deep Autoencoders_ , _Phys. Rev. D_ 101 (2020) 075021 [1808.08992].
  * (8) T. Heimel, G. Kasieczka, T. Plehn and J.M. Thompson, _QCD or What?_ , _SciPost Phys._ 6 (2019) 030 [1808.08979].
  * (9) T.S. Roy and A.H. Vijay, _A robust anomaly finder based on autoencoders_ , 1903.02032.
  * (10) O. Cerri, T.Q. Nguyen, M. Pierini, M. Spiropulu and J.-R. Vlimant, _Variational Autoencoders for New Physics Mining at the Large Hadron Collider_ , _JHEP_ 05 (2019) 036 [1811.10276].
  * (11) A. Blance, M. Spannowsky and P. Waite, _Adversarially-trained autoencoders for robust unsupervised new physics searches_ , _JHEP_ 10 (2019) 047 [1905.10384].
  * (12) J. Hajer, Y.-Y. Li, T. Liu and H. Wang, _Novelty Detection Meets Collider Physics_ , _Phys. Rev. D_ 101 (2020) 076015 [1807.10261].
  * (13) A. De Simone and T. Jacques, _Guiding New Physics Searches with Unsupervised Learning_ , _Eur. Phys. J. C_ 79 (2019) 289 [1807.06038].
  * (14) A. Mullin, S. Nicholls, H. Pacey, M. Parker, M. White and S. Williams, _Does SUSY have friends? A new approach for LHC event analysis_ , _JHEP_ 02 (2021) 160 [1912.10625].
  * (15) B.M. Dillon, D.A. Faroughy and J.F. Kamenik, _Uncovering latent jet substructure_ , _Phys. Rev. D_ 100 (2019) 056002 [1904.04200].
  * (16) J.A. Aguilar-Saavedra, J.H. Collins and R.K. Mishra, _A generic anti-QCD jet tagger_ , _JHEP_ 11 (2017) 163 [1709.01087].
  * (17) M. Romão Crispim, N.F. Castro, R. Pedro and T. Vale, _Transferability of Deep Learning Models in Searches for New Physics at Colliders_ , _Phys. Rev. D_ 101 (2020) 035042 [1912.04220].
  * (18) M. Crispim Romão, N.F. Castro, J.G. Milhano, R. Pedro and T. Vale, _Use of a generalized energy Mover’s distance in the search for rare phenomena at colliders_ , _Eur. Phys. J. C_ 81 (2021) 192 [2004.09360].
  * (19) O. Amram and C.M. Suarez, _Tag N’ Train: a technique to train improved classifiers on unlabeled data_ , _JHEP_ 01 (2021) 153 [2002.12376].
  * (20) T. Cheng, J.-F. Arguin, J. Leissner-Martin, J. Pilette and T. Golling, _Variational autoencoders for anomalous jet tagging_ , _Phys. Rev. D_ 107 (2023) 016002 [2007.01850].
  * (21) C.K. Khosa and V. Sanz, _Anomaly Awareness_ , _SciPost Phys._ 15 (2023) 053 [2007.14462].
  * (22) P. Thaprasop, K. Zhou, J. Steinheimer and C. Herold, _Unsupervised Outlier Detection in Heavy-Ion Collisions_ , _Phys. Scripta_ 96 (2021) 064003 [2007.15830].
  * (23) S. Alexander, S. Gleyzer, H. Parul, P. Reddy, M.W. Toomey, E. Usai et al., _Decoding Dark Matter Substructure without Supervision_ , 2008.12731.
  * (24) V. Mikuni and F. Canelli, _Unsupervised clustering for collider physics_ , _Phys. Rev. D_ 103 (2021) 092007 [2010.07106].
  * (25) M. van Beekveld, S. Caron, L. Hendriks, P. Jackson, A. Leinweber, S. Otten et al., _Combining outlier analysis algorithms to identify new physics at the LHC_ , _JHEP_ 09 (2021) 024 [2010.07940].
  * (26) S.E. Park, D. Rankin, S.-M. Udrescu, M. Yunus and P. Harris, _Quasi Anomalous Knowledge: Searching for new physics with embedded knowledge_ , _JHEP_ 21 (2020) 030 [2011.03550].
  * (27) D.A. Faroughy, _Uncovering hidden new physics patterns in collider events using Bayesian probabilistic models_ , _PoS_ ICHEP2020 (2021) 238 [2012.08579].
  * (28) T. Golling et al., _The Mass-ive Issue: Anomaly Detection in Jet Physics_ , 2303.14134.
  * (29) P. Chakravarti, M. Kuusela, J. Lei and L. Wasserman, _Model-Independent Detection of New Physics Signals Using Interpretable Semi-Supervised Classifier Tests_ , 2102.07679.
  * (30) J. Batson, C.G. Haaf, Y. Kahn and D.A. Roberts, _Topological Obstructions to Autoencoding_ , _JHEP_ 04 (2021) 280 [2102.08380].
  * (31) A. Blance and M. Spannowsky, _Unsupervised event classification with graphs on classical and photonic quantum computers_ , _JHEP_ 21 (2020) 170 [2103.03897].
  * (32) B. Bortolato, A. Smolkovič, B.M. Dillon and J.F. Kamenik, _Bump hunting in latent space_ , _Phys. Rev. D_ 105 (2022) 115009 [2103.06595].
  * (33) J.H. Collins, P. Martín-Ramiro, B. Nachman and D. Shih, _Comparing weak- and unsupervised methods for resonant anomaly detection_ , _Eur. Phys. J. C_ 81 (2021) 617 [2104.02092].
  * (34) B.M. Dillon, T. Plehn, C. Sauer and P. Sorrenson, _Better Latent Spaces for Better Autoencoders_ , _SciPost Phys._ 11 (2021) 061 [2104.08291].
  * (35) T. Finke, M. Krämer, A. Morandini, A. Mück and I. Oleksiyuk, _Autoencoders for unsupervised anomaly detection in high energy physics_ , _JHEP_ 06 (2021) 161 [2104.09051].
  * (36) D. Shih, M.R. Buckley, L. Necib and J. Tamanas, _via machinae: Searching for stellar streams using unsupervised machine learning_ , _Mon. Not. Roy. Astron. Soc._ 509 (2021) 5992 [2104.12789].
  * (37) O. Atkinson, A. Bhardwaj, C. Englert, V.S. Ngairangbam and M. Spannowsky, _Anomaly detection with convolutional Graph Neural Networks_ , _JHEP_ 08 (2021) 080 [2105.07988].
  * (38) A. Kahn, J. Gonski, I. Ochoa, D. Williams and G. Brooijmans, _Anomalous jet identification via sequence modeling_ , _JINST_ 16 (2021) P08012 [2105.09274].
  * (39) T. Dorigo, M. Fumanelli, C. Maccani, M. Mojsovska, G.C. Strong and B. Scarpa, _RanBox: anomaly detection in the copula space_ , _JHEP_ 01 (2023) 008 [2106.05747].
  * (40) S. Caron, L. Hendriks and R. Verheyen, _Rare and Different: Anomaly Scores from a combination of likelihood and out-of-distribution models to detect new physics at the LHC_ , _SciPost Phys._ 12 (2022) 077 [2106.10164].
  * (41) E. Govorkova, E. Puljak, T. Aarrestad, M. Pierini, K.A. Woźniak and J. Ngadiuba, _LHC physics dataset for unsupervised New Physics detection at 40 MHz_ , _Sci. Data_ 9 (2022) 118 [2107.02157].
  * (42) G. Kasieczka, B. Nachman and D. Shih, _New Methods and Datasets for Group Anomaly Detection From Fundamental Physics_ , in _Conference on Knowledge Discovery and Data Mining_ , 7, 2021 [2107.02821].
  * (43) S. Volkovich, F. De Vito Halevy and S. Bressler, _A data-directed paradigm for BSM searches: the bump-hunting example_ , _Eur. Phys. J. C_ 82 (2022) 265 [2107.11573].
  * (44) E. Govorkova et al., _Autoencoders on field-programmable gate arrays for real-time, unsupervised new physics detection at 40 MHz at the Large Hadron Collider_ , _Nature Mach. Intell._ 4 (2022) 154 [2108.03986].
  * (45) B. Ostdiek, _Deep Set Auto Encoders for Anomaly Detection in Particle Physics_ , _SciPost Phys._ 12 (2022) 045 [2109.01695].
  * (46) K. Fraser, S. Homiller, R.K. Mishra, B. Ostdiek and M.D. Schwartz, _Challenges for unsupervised anomaly detection in particle physics_ , _JHEP_ 03 (2022) 066 [2110.06948].
  * (47) P. Jawahar, T. Aarrestad, N. Chernyavskaya, M. Pierini, K.A. Wozniak, J. Ngadiuba et al., _Improving Variational Autoencoders for New Physics Detection at the LHC With Normalizing Flows_ , _Front. Big Data_ 5 (2022) 803685 [2110.08508].
  * (48) J. Herrero-Garcia, R. Patrick and A. Scaffidi, _A semi-supervised approach to dark matter searches in direct detection data with machine learning_ , _JCAP_ 02 (2022) 039 [2110.12248].
  * (49) J.A. Aguilar-Saavedra, _Anomaly detection from mass unspecific jet tagging_ , _Eur. Phys. J. C_ 82 (2022) 130 [2111.02647].
  * (50) R. Tombs and C.G. Lester, _A method to challenge symmetries in data with self-supervised learning_ , _JINST_ 17 (2022) P08024 [2111.05442].
  * (51) C.G. Lester and R. Tombs, _Using unsupervised learning to detect broken symmetries, with relevance to searches for parity violation in nature. (Previously: ”Stressed GANs snag desserts”)_ , 2111.00616.
  * (52) V. Mikuni, B. Nachman and D. Shih, _Online-compatible unsupervised nonresonant anomaly detection_ , _Phys. Rev. D_ 105 (2022) 055006 [2111.06417].
  * (53) S. Chekanov and W. Hopkins, _Event-Based Anomaly Detection for Searches for New Physics_ , _Universe_ 8 (2022) 494 [2111.12119].
  * (54) R.T. d’Agnolo, G. Grosso, M. Pierini, A. Wulzer and M. Zanetti, _Learning new physics from an imperfect machine_ , _Eur. Phys. J. C_ 82 (2022) 275 [2111.13633].
  * (55) F. Canelli, A. de Cosa, L.L. Pottier, J. Niedziela, K. Pedro and M. Pierini, _Autoencoders for semivisible jet detection_ , _JHEP_ 02 (2022) 074 [2112.02864].
  * (56) V.S. Ngairangbam, M. Spannowsky and M. Takeuchi, _Anomaly detection in high-energy physics using a quantum autoencoder_ , _Phys. Rev. D_ 105 (2022) 095004 [2112.04958].
  * (57) L. Bradshaw, S. Chang and B. Ostdiek, _Creating simple, interpretable anomaly detectors for new physics in jet substructure_ , _Phys. Rev. D_ 106 (2022) 035014 [2203.01343].
  * (58) J.A. Aguilar-Saavedra, _Taming modeling uncertainties with mass unspecific supervised tagging_ , _Eur. Phys. J. C_ 82 (2022) 270 [2201.11143].
  * (59) T. Buss, B.M. Dillon, T. Finke, M. Krämer, A. Morandini, A. Mück et al., _What’s anomalous in LHC jets?_ , _SciPost Phys._ 15 (2023) 168 [2202.00686].
  * (60) S. Alvi, C.W. Bauer and B. Nachman, _Quantum anomaly detection for collider physics_ , _JHEP_ 02 (2023) 220 [2206.08391].
  * (61) B.M. Dillon, R. Mastandrea and B. Nachman, _Self-supervised anomaly detection for new physics_ , _Phys. Rev. D_ 106 (2022) 056005 [2205.10380].
  * (62) M. Birman, B. Nachman, R. Sebbah, G. Sela, O. Turetz and S. Bressler, _Data-directed search for new physics based on symmetries of the SM_ , _Eur. Phys. J. C_ 82 (2022) 508 [2203.07529].
  * (63) M. Letizia, G. Losapio, M. Rando, G. Grosso, A. Wulzer, M. Pierini et al., _Learning new physics efficiently with nonparametric methods_ , _Eur. Phys. J. C_ 82 (2022) 879 [2204.02317].
  * (64) C. Fanelli, J. Giroux and Z. Papandreou, _‘Flux+Mutability’: a conditional generative approach to one-class classification and anomaly detection_ , _Mach. Learn. Sci. Tech._ 3 (2022) 045012 [2204.08609].
  * (65) T. Finke, M. Krämer, M. Lipp and A. Mück, _Boosting mono-jet searches with model-agnostic machine learning_ , _JHEP_ 08 (2022) 015 [2204.11889].
  * (66) R. Verheyen, _Event Generation and Density Estimation with Surjective Normalizing Flows_ , _SciPost Phys._ 13 (2022) 047 [2205.01697].
  * (67) B.M. Dillon, L. Favaro, T. Plehn, P. Sorrenson and M. Krämer, _A Normalized Autoencoder for LHC Triggers_ , 2206.14225.
  * (68) S. Caron, R.R. de Austri and Z. Zhang, _Mixture-of-Theories training: can we find new physics and anomalies better by mixing physical theories?_ , _JHEP_ 03 (2023) 004 [2207.07631].
  * (69) S.E. Park, P. Harris and B. Ostdiek, _Neural embedding: learning the embedding of the manifold of physics data_ , _JHEP_ 07 (2023) 108 [2208.05484].
  * (70) J.F. Kamenik and M. Szewc, _Null hypothesis test for anomaly detection_ , _Phys. Lett. B_ 840 (2023) 137836 [2210.02226].
  * (71) G. Kasieczka, R. Mastandrea, V. Mikuni, B. Nachman, M. Pettee and D. Shih, _Anomaly detection under coordinate transformations_ , _Phys. Rev. D_ 107 (2023) 015009 [2209.06225].
  * (72) J.Y. Araz and M. Spannowsky, _Quantum-probabilistic Hamiltonian learning for generative modelling & anomaly detection_, 2211.03803.
  * (73) J. Schuhmacher, L. Boggia, V. Belis, E. Puljak, M. Grossi, M. Pierini et al., _Unravelling physics beyond the standard model with classical and quantum anomaly detection_ , _Mach. Learn. Sci. Tech._ 4 (2023) 045031 [2301.10787].
  * (74) S. Roche, Q. Bayer, B. Carlson, W. Ouligian, P. Serhiayenka, J. Stelzer et al., _Nanosecond anomaly detection with decision trees for high energy physics and real-time application to exotic Higgs decays_ , 2304.03836.
  * (75) L. Vaslin, V. Barra and J. Donini, _GAN-AE: an anomaly detection algorithm for New Physics search in LHC data_ , _Eur. Phys. J. C_ 83 (2023) 1008 [2305.15179].
  * (76) ATLAS Collaboration, _Anomaly detection search for new resonances decaying into a Higgs boson and a generic new particle $X$ in hadronic final states using $\sqrt{s}=13$ TeV $pp$ collisions with the ATLAS detector_, _Phys. Rev. D_ 108 (2023) 052009 [2306.03637].
  * (77) S.V. Chekanov and R. Zhang, _Enhancing the hunt for new phenomena in dijet final-states using anomaly detection filter at High-Luminosity Large Hadron Collider_ , 2308.02671.
  * (78) CMS ECal Collaboration, _Autoencoder-based Anomaly Detection System for Online Data Quality Monitoring of the CMS Electromagnetic Calorimeter_ , 2309.10157.
  * (79) G. Bickendorf, M. Drees, G. Kasieczka, C. Krause and D. Shih, _Combining Resonant and Tail-based Anomaly Detection_ , 2309.12918.
  * (80) M. Freytsis, M. Perelstein and Y.C. San, _Anomaly Detection in Presence of Irrelevant Features_ , 2310.13057.
  * (81) E.M. Metodiev, J. Thaler and R. Wynne, _Anomaly Detection in Collider Physics via Factorized Observables_ , 2312.00119.
  * (82) G. Karagiorgi, G. Kasieczka, S. Kravitz, B. Nachman and D. Shih, _Machine learning in the search for new fundamental physics_ , _Nature Reviews Physics_ 4 (2022) 399 [2112.03769].
  * (83) G. Kasieczka, B. Nachman, D. Shih, O. Amram, A. Andreassen, K. Benkendorfer et al., _The LHC olympics 2020 a community challenge for anomaly detection in high energy physics_ , _Reports on Progress in Physics_ 84 (2021) 124201 [2101.08320].
  * (84) T. Aarrestad, M. van Beekveld, M. Bona, A. Boveia, S. Caron, J. Davies et al., _The Dark Machines Anomaly Score Challenge: Benchmark Data and Model Independent Event Classification for the Large Hadron Collider_ , _SciPost Phys._ 12 (2022) 043 [2105.14027].
  * (85) V. Belis, P. Odagiu and T.K. Aarrestad, _Machine learning for anomaly detection in particle physics_ , _Rev. Phys._ 12 (2024) 100091 [2312.14190].
  * (86) J.H. Collins, P. Martín-Ramiro, B. Nachman and D. Shih, _Comparing weak- and unsupervised methods for resonant anomaly detection_ , _The European Physical Journal C_ 81 (2021)  [2104.02092].
  * (87) T. Golling, G. Kasieczka, C. Krause, R. Mastandrea, B. Nachman, J.A. Raine et al., _The Interplay of Machine Learning–based Resonant Anomaly Detection Methods_ , 2307.11157.
  * (88) ATLAS Collaboration, _Dijet resonance search with weak supervision using $\sqrt{s}=13$ TeV $pp$ collisions in the ATLAS detector_, _Physical Review Letters_ 125 (2020) 131801 [2005.02983].
  * (89) ATLAS Collaboration, _Search for new phenomena in two-body invariant mass distributions using unsupervised machine learning for anomaly detection at $\sqrt{s}=13$ TeV with the ATLAS detector_, 2023.  All figures including auxiliary figures are available at https://atlas.web.cern.ch/Atlas/GROUPS/PHYSICS/CONFNOTES/ATLAS-CONF-2023-022.
  * (90) G. Kasieczka, B. Nachman and D. Shih, _Official Datasets for LHC Olympics 2020 Anomaly Detection Challenge (Version v6)_ , 2019.  10.5281/zenodo.4536624.
  * (91) J.H. Kim, K. Kong, B. Nachman and D. Whiteson, _The motivation and status of two-body resonance decays after the LHC run 2 and beyond_ , _Journal of High Energy Physics_ 2020 (2020)  [1907.06659].
  * (92) T. Sjöstrand, S. Ask, J.R. Christiansen, R. Corke, N. Desai, P. Ilten et al., _An introduction to PYTHIA 8.2_ , _Comput. Phys. Commun._ 191 (2015) 159 [1410.3012].
  * (93) J. de Favereau, C. Delaere, P. Demin, A. Giammanco, V. Lemaître, A. Mertens et al., _DELPHES 3: a modular framework for fast simulation of a generic collider experiment_ , _Journal of High Energy Physics_ 2014 (2014) 026.
  * (94) A. Mertens, _New features in delphes 3_ , _Journal of Physics: Conference Series_ 608 (2015) 012045.
  * (95) M. Selvaggi, _Delphes 3: A modular framework for fast-simulation of generic collider experiments_ , _Journal of Physics: Conference Series_ 523 (2014) 012033.
  * (96) M. Cacciari, G.P. Salam and G. Soyez, _The anti-ktjet clustering algorithm_ , _Journal of High Energy Physics_ 2008 (2008) 063–063 [0802.1189].
  * (97) M. Cacciari, G.P. Salam and G. Soyez, _Fastjet user manual_ , _The European Physical Journal C_ 72 (2012)  [1111.6097].
  * (98) E. Rodrigues et al., _The Scikit HEP Project – overview and prospects_ , _EPJ Web Conf._ 245 (2020) 06028 [2007.03577].
  * (99) L. de Oliveira, M. Kagan, L. Mackey, B. Nachman and A. Schwartzman, _Jet-images — deep learning edition_ , _Journal of High Energy Physics_ 2016 (2016)  [1511.05190].
  * (100) G. Kasieczka, T. Plehn, A. Butter, K. Cranmer, D. Debnath, B.M. Dillon et al., _The machine learning landscape of top taggers_ , _SciPost Physics_ 7 (2019)  [1902.09914].
  * (101) S. Macaluso and D. Shih, _Pulling out all the tops with computer vision and deep learning_ , _Journal of High Energy Physics_ 2018 (2018)  [1803.00107].
  * (102) T. Finke, M. Krämer, A. Morandini, A. Mück and I. Oleksiyuk, _Autoencoders for unsupervised anomaly detection in high energy physics_ , _JHEP_ 2021 (2021)  [2104.09051].
  * (103) A. Hallin, J. Isaacson, G. Kasieczka, C. Krause, B. Nachman, T. Quadfasel et al., _Classifying anomalies through outer density estimation_ , _Phys. Rev. D_ 106 (2022) 055006 [2109.00546].
  * (104) J.A. Raine, S. Klein, D. Sengupta and T. Golling, _CURTAINs for your sliding window: Constructing unobserved regions by transforming adjacent intervals_ , _Front. Big Data_ 6 (2023) 899345 [2203.09470].
  * (105) CMS Collaboration, _Search for narrow resonances in the b-tagged dijet mass spectrum in proton-proton collisions at s=13 TeV_ , _Phys. Rev. D_ 108 (2023) 012009 [2205.01835].
  * (106) ATLAS Collaboration, _Search for new phenomena in the dijet mass distribution using $pp$ collision data at $\sqrt{s}=8\text{ }\text{ }\mathrm{TeV}$ with the atlas detector_, _Phys. Rev. D_ 91 (2015) 052007 [1407.1376].
  * (107) ATLAS Collaboration, _Search for new phenomena in dijet mass and angular distributions from pp collisions at $\sqrt{s}=$ 13 tev with the ATLAS detector_, _Physics Letters B_ 754 (2016) 302–322 [1512.01530].
  * (108) CMS Collaboration, _Search for dijet resonances in proton–proton collisions at $\sqrt{s}=13$ tev and constraints on dark matter and other models_, _Physics Letters B_ 769 (2017) 520–542.
  * (109) CMS Collaboration, _Search for narrow resonances decaying to dijets in proton-proton collisions at $\sqrt{s}=13$ tev_, _Physical Review Letters_ 116 (2016)  [1512.01224].
  * (110) CMS Collaboration, _Search for narrow resonances in dijet final states at $\sqrt{s}=8$ tev with the novel cms technique of data scouting_, _Physical Review Letters_ 117 (2016) .
  * (111) ATLAS Collaboration, _A search for new physics in dijet mass and angular distributions in pp collisions at $\sqrt{s}=7$ tev measured with the ATLAS detector_, _New Journal of Physics_ 13 (2011) 053044 [1512.01530].
  * (112) CMS Collaboration, _Search for resonances and quantum black holes using dijet mass spectra in proton-proton collisions at $\sqrt{s}=8\text{ }\text{ }\mathrm{TeV}$_, _Phys. Rev. D_ 91 (2015) 052009 [1501.04198].
  * (113) CMS Collaboration, _Search for narrow resonances using the dijet mass spectrum in $pp$ collisions at $\sqrt{s}\mathbf{=}8\text{ }\text{ }\mathrm{TeV}$_, _Phys. Rev. D_ 87 (2013) 114015 [1302.4794].
  * (114) CMS Colaboration, _Search for narrow resonances and quantum black holes in inclusive and b-tagged dijet mass spectra from pp collisions at $\sqrt{s}=7$ tev_, _Journal of High Energy Physics_ 2013 (2013)  [1210.2387].
  * (115) ATLAS Collaboration, _Search for new physics in the dijet mass distribution using 1 $\text{fb}^{-1}$ of pp collision data at $\sqrt{s}=7$ tev collected by the atlas detector_, _Physics Letters B_ 708 (2012) 37 [1108.6311].
  * (116) ATLAS Collaboration, _Atlas search for new phenomena in dijet mass and angular distributions using pp collisions at $\sqrt{s}=7$tev_, _Journal of High Energy Physics_ 2013 (2013)  [1210.1718].
  * (117) CMS Collaboration, _Search for resonances in the dijet mass spectrum from 7 tev pp collisions at cms_ , _Physics Letters B_ 704 (2011) 123–142 [1107.4771].
  * (118) CDF Collaboration, _Search for new particles decaying into dijets in proton-antiproton collisions at $\sqrt{s}=1.96$ tev_, _Physical Review D_ 79 (2009)  [0812.4036].
  * (119) ATLAS Collaboration, _A search for new physics in dijet mass and angular distributions inppcollisions at $\sqrt{s}=7$ tev measured with the atlas detector_, _New Journal of Physics_ 13 (2011) 053044 [1103.3864].
  * (120) CMS Collaboration, _Search for dijet resonances in 7 tev $pp$ collisions at cms_, _Phys. Rev. Lett._ 105 (2010) 211801 [1010.0203].
  * (121) ATLAS Collaboration, _Search for low-mass dijet resonances using trigger-level jets with the atlas detector in $pp$ collisions at $\sqrt{s}=13\text{ }\text{ }\mathrm{TeV}$_, _Physical Review Letters_ 121 (2018)  [1804.03496].
  * (122) K. Sekhon, R.C. Edgar and D. Amidei, _SWiFt: Sliding Window Fit Method for Resonance Searches_ , Tech. Rep. ATL-COM-PHYS-2018-161, CERN, Geneva (2018).
  * (123) J. Alison, F.A. Di Bello, S. Chekanov, A. Coccaro, T. Golling, G. Facini et al., _Search for resonances in the di-jet mass distribution with one or two jets identified as $b$-jets in proton-proton collisions at $\sqrt{s}$=13 TeV with the ATLAS detector_, Tech. Rep. ATL-COM-PHYS-2016-1561, CERN, Geneva (2016).
  * (124) G. Kasieczka, B. Nachman, M.D. Schwartz and D. Shih, _Automating the abcd method with machine learning_ , _Phys. Rev. D_ 103 (2021) 035021 [2007.14400].
  * (125) L. Vaslin, S. Calvet, V. Barra and J. Donini, _pybumphunter: A model independent bump hunting tool in python for high energy physics analyses_ , 2208.14760.
  * (126) F. Pedregosa et al., _Scikit-learn: Machine Learning in Python_ , _J. Machine Learning Res._ 12 (2011) 2825 [1201.0490].
  * (127) D. Arthur and S. Vassilvitskii, _K-means++: The advantages of careful seeding_ , in _Proceedings of the Eighteenth Annual ACM-SIAM Symposium on Discrete Algorithms_ , SODA ’07, (USA), p. 1027–1035, Society for Industrial and Applied Mathematics, 2007.
  * (128) A. Andreassen, B. Nachman and D. Shih, _Simulation Assisted Likelihood-free Anomaly Detection_ , _Phys. Rev. D_ 101 (2020) 095004 [2001.05001].
  * (129) T. Golling, S. Klein, R. Mastandrea and B. Nachman, _Flow-enhanced transportation for anomaly detection_ , _Phys. Rev. D_ 107 (2023) 096025 [2212.11285].
  * (130) D. Sengupta, S. Klein, J.A. Raine and T. Golling, _CURTAINs Flows For Flows: Constructing Unobserved Regions with Maximum Likelihood Estimation_ , 2305.04646.
  * (131) E. Buhmann, C. Ewen, G. Kasieczka, V. Mikuni, B. Nachman and D. Shih, _Full phase space resonant anomaly detection_ , 2023.
  * (132) D. Sengupta, M. Leigh, J.A. Raine, S. Klein and T. Golling, _Improving new physics searches with diffusion models for event observables and jet constituents_ , 2312.10130.
  * (133) M. Frate, K. Cranmer, S. Kalia, A. Vandenberg-Rodes and D. Whiteson, _Modeling Smooth Backgrounds and Generic Localized Signals with Gaussian Processes_ , 1709.05681.
  * (134) K. Cranmer and I. Yavin, _RECAST: Extending the Impact of Existing Analyses_ , _JHEP_ 04 (2011) 038 [1010.2506].
  * (135) P. Virtanen, R. Gommers, T.E. Oliphant, M. Haberland, T. Reddy, D. Cournapeau et al., _SciPy 1.0: Fundamental Algorithms for Scientific Computing in Python_ , _Nature Methods_ 17 (2020) 261.