# Interpretable Machine Learning Methods
Applied to Jet Background Subtraction in Heavy Ion Collisions

Tanner Mengel  Patrick Steffanic  Charles Hughes  Antonio Carlos Oliveira da
Silva  Christine Nattrass  University of Tennessee, Knoxville, TN, USA-37996.

###### Abstract

Jet measurements in heavy ion collisions can provide constraints on the
properties of the quark gluon plasma, but the kinematic reach is limited by a
large, fluctuating background. We present a novel application of symbolic
regression to extract a functional representation of a deep neural network
trained to subtract background from jets in heavy ion collisions. We show that
the deep neural network is approximately the same as a method using the
particle multiplicity in a jet. This demonstrates that interpretable machine
learning methods can provide insight into underlying physical processes.

##  I Introduction

The Quark Gluon Plasma (QGP) is a hot, dense, strongly interacting liquid of
quarks and gluons that is created briefly in high energy heavy ion collisions
[1, 2, 3, 4]. Measurements of jets produced by hard scatterings between
partons in heavy ion collisions can be used to investigate the properties of
the QGP [5]. Quantitative comparisons between jet measurements and physics
models can provide further constraints on these properties [6, 7]. However,
heavy ion events are dominated by a fluctuating background of soft particles
not due to hard scatterings. The details of these fluctuations are sensitive
to correlations from hydrodynamical flow and the shape of the single particle
spectra [8], and as such are unlikely to be exactly the same in data and
models. Mixed events are able to successfully describe the background in
measurements of hadron-jet correlations by the STAR collaboration [9] at the
Relativistic Heavy Ion Collider (RHIC). Studies of the background at the Large
Hadron Collider (LHC) by the ALICE Collaboration found that the distribution
of background energy density in random cones is well described by a random
background with correlations due to hydrodynamical flow and Poissonian
fluctuations [10]. A better understanding of this background will facilitate
more precise jet measurements for comparisons between data and models.

Measurement precision and kinematic range is limited by the ability to correct
for this background and its fluctuations. Background correction in jet
measurements requires subtraction of contributions from soft particles within
the jet, and suppression of fluctuations which have been reconstructed as
combinatorial jets. At low momenta, combinatorial jets limit the kinematic
reach of the measurement. Improved background subtraction methods would
increase measurements’ sensitivity to partonic energy loss. Measurements of
jet spectra which extend to low momenta primarily use the area method [11] for
background subtraction. This method was initially proposed to correct for the
underlying event in $p$+$p$ collisions in high pile-up conditions [11] and has
also been applied to heavy ion collisions [12, 13, 14, 15].

The complexity of jet background subtraction makes it an interesting
environment to apply machine learning techniques. However, application of
machine learning methods to background subtraction should be handled with care
since models are not able to fully reproduce background fluctuations in heavy
ion collisions [8]. Nuclear physics has prioritized the continued advancement
in machine learning analysis techniques with a focus on interpretable methods
that are robust, provide clear uncertainty quantification, and are explainable
[16]. Applications of non-interpretable machine learning methods are
insufficient when models available for training may be inaccurate, when it may
be necessary to understand the method to interpret the results, or when a
result is needed outside of the training space.

Application of a deep neural network, i.e. a neural network with multiple
hidden layers, to jet background subtraction in heavy ion collisions has
demonstrated significant improvements compared to the area method,
particularly at low jet momenta [17, 18]. Deep neural networks are susceptible
to model bias because their predictions risk being unreliable outside the
domain of their training space. These methods may break down when they are
extrapolated beyond this space, and due to their opaque nature, offer little
indication where and why this break down occurs. In addition, one cannot
validate the technique against data because we do not know the true jet
momenta in data.

Increased performance of machine learning methods over traditional methods is
an indication that there is information accessible to the machine learning
that accounts for this improvement. We present an interpretable machine
learning technique that allows us to understand why a deep neural network
improves the jet momentum resolution in heavy ion collisions. We empirically
derive an alternate method based on the background described in [19, 10], we
call the multiplicity method. We compare the widths of the fluctuations of the
jet momenta for the this method to the area and neural network methods and
estimate the impact of the methods on the kinematic range. We apply symbolic
regression to determine a functional form describing the mapping learned by
the neural network, which was trained using TennGen [8] for the background and
PYTHIA [20] for the signal. We compare this functional description of the
neural network to the form of the multiplicity method.

##  II Method

###  II.1 Simulation

TennGen [21, 22] generates heavy ion collisions with $\pi^{\pm}$, K±, p and
$\bar{p}$ hadrons with yields [23], momentum distributions [24, 25], and
azimuthal anisotropies [26, 27] matched to published data. TennGen was updated
to simulate collision energies per nucleon of $\sqrt{s_{\mathrm{NN}}}$ = 200
GeV collisions as well as $\sqrt{s_{\mathrm{NN}}}$ = 2.76 TeV, including
multiplicity fluctuations, and improved computational efficiency. Proton-
proton collisions at $\sqrt{s}$ = 200 GeV were simulated with the PYTHIA 8.307
[20] Monash 2013 tune [28] in 25 $p_{T}^{hard}$ bins starting at 5 GeV, with 1
million $p$+$p$ events in each bin. Only final state charged particles from
PYTHIA are mixed with a TennGen background event. Charged particles from both
PYTHIA and TennGen are required to have a minimum $p_{T}$ of 150 MeV and be
within pseudo-rapidity $|\eta|<$ 0.9.

Jets are clustered using the anti-$k_{T}$ algorithm with FastJet 3.4.0 [29]
with jet resolution parameters $R$ = 0.2, 0.4, and 0.6. To determine the true
momentum, jets are reconstructed separately in both PYTHIA and the combined
event. Jets in the combined PYTHIA and TennGen event are geometrically matched
to a PYTHIA jet if $\Delta R=\sqrt{\Delta\eta^{2}+\Delta\phi^{2}}<0.1$ where
$\Delta\eta$ and $\Delta\phi$ are the differences in $\eta$ and $\phi$ between
the jets and there is a bijective match. Reconstructed jets are required to
have $p_{T}$ $>$ 5 GeV and be within pseudo-rapidity $|\eta_{jet}|<0.9-R$. The
momentum of the PYTHIA jet is taken as the truth momentum,
$p_{T,Jet}^{Truth}\equiv p_{T,Jet}^{PYTHIA}$.

The improvement in jet momentum resolution extends the kinematic range of the
measurement to lower jet momenta. Increased precision in jet momentum allows
for combinatorial jets to be reconstructed in jet momentum bins closer to
zero. The lower threshold for unfolding is typically set to be between 2-5
times the width of the jet momentum resolution to suppress effects of
combinatorial jets on the unfolded results [12, 18]. The suppression of
combinatorial jets for each method is demonstrated by taking the ratio of
unmatched reconstructed jet spectra including combinatorial background to
matched reconstructed PYTHIA jet spectra.

###  II.2 Area and multiplicity methods

For area-based background subtraction, the jet area [30] is estimating through
the use of “ghost” particles, jets are reconstructed using the $k_{T}$ jet
finder [31]. The corrected jet momentum is then estimated as

$p_{T,Jet}^{Corr,A}=p_{T,Jet}^{tot}-\rho A,$ (1)

where $A$ is the jet area, $\rho$ is the background momentum density per unit
area, and $p_{T,Jet}^{tot}$ is the total momentum in the jet. The $\rho$ in an
event is approximated as the median $p_{T,Jet}^{tot}/A$ for $k_{T}$ jets
because $k_{T}$ jets are dominated by background.

To a good approximation, the variance of the momentum residual $\delta
p_{T}=p_{T}^{Corr}-p_{T}^{Truth}$ with the area method is given by

$\sigma_{\delta
p_{T}}=\sqrt{N\sigma^{2}_{p_{T}}+(N+2N^{2}\sum_{n=1}^{\infty}v_{n}^{2})\langle
p_{T}\rangle^{2}},$ (2)

where $N$ is the number of background particles in the jet, $\sigma_{p_{T}}$
is the standard deviation of the single track momentum distribution, $v_{n}$
are the coefficients of the azimuthal anisotropies of the single particle
distributions, and $\langle p_{T}\rangle$ is the average momentum of
background particles [10]. This is derived by assuming each of the $N$
particles is drawn from a single track momentum distribution which is
approximately a Gamma distribution, giving rise to the first term [19]. The
second term is from Poissonian fluctuations in the number of background
particles and the third term is from fluctuations in the number of particles
due to hydrodynamical flow. Deviations of the single track momentum
distribution from a Gamma distribution and momentum dependence of the $v_{n}$
lead to slightly larger widths [8].

The area method is usually used instead of iterative background subtraction
methods [32, 33, 34] for measurements of jets at lower momenta. Iterative
methods may suppress the fluctuations described in eq. 2 by estimating the
local background and suppress combinatorial jets by requiring high momentum or
energy constituents. At low momenta, these requirements may impose a bias on
the surviving jets. Fluctuations and the contribution from combinatorial jets
are generally higher with the area method, but with less bias.

We propose a multiplicity-based method as an alternative to the area method

$p_{T,Jet}^{Corr,N}=p_{T,Jet}^{tot}-\rho_{Mult}(N_{tot}-N_{signal}),$ (3)

where $N_{tot}$ is the total number of particles in the jet, $N_{signal}$ is
the number of particles in the signal, and $N=N_{tot}-N_{signal}$. This
leverages the fact that the natural variable in the variance is the number of
background particles, largely eliminating the second and third terms in eq. 2.
The $\rho_{Mult}$ in an event is the mean transverse momentum per background
particle, which is approximated as the median $p_{T,Jet}^{tot}/N_{tot}$ for
$k_{T}$ jets. $N_{signal}$ is roughly described by models [35] and therefore
can be estimated. This method would introduce an additional systematic
uncertainty, which can be estimated using measurements of fragmentation
functions to be around one additional particle in $N_{signal}$ [36]. The
variance of the momentum residual in eq. 2 is then generally smaller for the
multiplicity method than the area method for 0–80% central collisions.

###  II.3 Machine learning methods

A sufficiently complex neural network can interpolate any function, at the
cost of transparency to the user. This poses an obstacle to application of
deep neural networks in physics where understanding predictions and
identifying their potential biases is crucial. Our approach to addressing this
challenge is through symbolic regression, one example of interpretable machine
learning, to extract mathematical expressions from trained deep neural
networks. The resulting equations provide an effective description of the
neural network’s mapping between the input and output. By constraining the
types of operations available, we can impose complexity and smoothness
requirements.

We train a deep neural network to predict the corrected jet momentum from the
following input features: the uncorrected jet momentum, jet area, jet
angularity, number of jet constituents, and seven leading constituent momenta.
The architecture and input features of the network are motivated by previous
application of neural networks to proton-proton jets with a thermal background
[17]. The deep neural network is implemented with TensorFlow 2.10.0 [37]. The
deep neural network has three hidden layers consisting of 100, 100 and 50
nodes, each activated by a rectified linear unit (ReLU) [38] function. The
model is optimized with ADAM [38] and the loss function is a modified mean
squared error

$\mathcal{L}=\langle||p_{T,Jet}^{Truth}-p_{T,Jet}^{DNN}||^{2}\rangle-\lambda\sum_{l=1}^{L}||\mathbf{W}_{l}||^{2},$
(4)

where $p_{T,Jet}^{DNN}$ is the predicted jet momentum, $p_{T,Jet}^{Truth}$ is
the truth momentum, the last term is an $L^{2}(\lambda)$ regularization where
$\lambda=0.001$, $\mathbf{W}_{l}$ is the weight matrix of layer $l$, and the
sum is over the $L$ layers. The regularization term penalizes redundancy and
encourages sparsity in the final trained network. The network is trained using
50% of the simulated jets while the remaining 50% are reserved for testing.

Once the neural network is trained, it represents an approximate mapping
between the input jet features and the truth jet momentum. We apply a genetic
algorithm to symbolically regress a functional form which describes this
mapping using the PySR 0.11.11 [39] package. The PySR model samples the phase
space of analytic expressions defined by operators, input features, and
constants for minimization through genetic programming. The input features are
comparable to those of the neural network, and the pool of operations are
arithmetic, exponential, trigonometric, and exponentiation. The model mutates
over 50 generations of 20 different population samples, with each population
containing 33 individuals. The loss function for the PySR model

$\mathcal{L}=\langle||p_{T,Jet}^{DNN}-p_{T,Jet}^{PySR}||^{2}\rangle,$ (5)

is the mean squared error between the prediction from PySR $p_{T,Jet}^{PySR}$
and the corrected jet momentum predicted by the neural network. PySR evaluates
expressions based on a score $S$ that rewards minimizing the loss function
$\mathcal{L}$ and penalizes equation complexity $C$

$S=-\frac{\delta\ln\mathcal{L}}{\delta C},$ (6)

where the equation complexity $C$ is defined as the total number of
operations, variables, and constants used in an equation [40]. The simulated
jets, designated for testing, are used to sufficiently sample the neural
network outputs throughout the possible input feature space. The highest
scoring PySR expression is a functional representation of the mapping from
input jet features to corrected jet momentum learned by the deep neural
network.

##  III Results

Figure 1 shows the width of the jet momentum residual distributions as a
function of jet momentum for each background subtraction method in both Au+Au
collisions at $\sqrt{s_{\mathrm{NN}}}$ = 200 GeV and Pb+Pb collisions at
$\sqrt{s_{\mathrm{NN}}}$ = 2.76 TeV. The $\sigma_{\delta p_{T}}$ increases
with increasing jet resolution parameter, as expected because there is more
background when the jet is larger. The $\sigma_{\delta p_{T}}$ also increases
with $\sqrt{s_{\mathrm{NN}}}$ because the particle multiplicity increases. As
seen in [17], the deep neural network reconstructs the momentum considerably
more accurately than the area method. The performance of the multiplicity
method is comparable to that of the deep neural network in Au+Au collisions
and small jet resolution parameter.

Figure 1: Comparisons of jet $p_{T}$ residual width for each background
subtraction method as a function of reconstructed jet momentum for Au+Au
collisions at $\sqrt{s_{\mathrm{NN}}}$ = 200 GeV and Pb+Pb collisions at
$\sqrt{s_{\mathrm{NN}}}$ = 2.76 TeV for jet resolution parameters $R=0.2,0.4,$
and $0.6$.

The ability of each method to sufficiently suppress contributions from
combinatorial jets at low $p_{T}$ is demonstrated with the ratios of the
reconstructed jet spectra to the true jet spectra, shown in Fig. 2. The
contributions from combinatorial jets decreases with increasing jet momentum
for all methods, with all jet resolution parameters, and for both collision
energies. The ratios for the deep neural network and multiplicity methods are
both lower than those of the area method.

Figure 2: Ratio of the reconstructed jet spectrum over the truth spectrum for
Au+Au collisions at $\sqrt{s_{\mathrm{NN}}}$ = 200 GeV and Pb+Pb collisions at
$\sqrt{s_{\mathrm{NN}}}$ = 2.76 TeV for jet resolution parameters $R=0.2,0.4,$
and $0.6$. Low momentum points for the area method for LHC energies at $R$ =
0.4 and $R$ = 0.6 are off scale.

For all jet resolution parameters and both collision energies, the symbolic
regression found that the best description of the deep neural network has the
functional form

$p_{T,Jet}^{Corr.PySR}=p_{T,Jet}^{tot}-C_{1}(N_{tot}-C_{2}),$ (7)

where the two parameters, $C_{1}$ and $C_{2}$, are optimization constants
defined by PySR. These parameters are plotted in Fig. 3 and compared to the
average value of the parameters used in the multiplicity method. We find that
the symbolic regression parameters $C_{1}$ and $C_{2}$ are comparable to the
averages of those for the multiplicity method, $\langle\rho_{Mult}\rangle$ and
$\langle N_{signal}\rangle$, respectively, with greater deviations at LHC
energies and larger $R$. This indicates that the deep neural network is using
a relationship similar to the multiplicity method to predict jet momenta.

Figure 3: PySR optimization constants compared to average value of
multiplicity method parameters versus jet resolution parameter for Au+Au
collisions at $\sqrt{s_{\mathrm{NN}}}$ = 200 GeV and Pb+Pb collisions at
$\sqrt{s_{\mathrm{NN}}}$ = 2.76 TeV for jet resolution parameters $R=0.2,0.4,$
and $0.6$.

This approach to machine learning enables use of domain knowledge. The
optimization parameters from PySR would otherwise not have a clear physical
interpretation. Since these parameters are understood in the multiplicity
method, it is possible to assign a physically motivated uncertainty to them.
Assumptions inherent in the method can then be understood.

##  IV Conclusions

We have shown that interpretable machine learning methods can learn an
underlying physical correlation, such as the multiplicity dependence for jet
background, that was previously overlooked. We previously showed that when we
used a random forest to classify jets as either combinatorial or signal, the
optimal selection was on the leading hadron momentum [41], already used as a
standard technique [12, 13, 14, 15]. We argue that applying machine learning
to scientific problems requires methods that are interpretable. The definition
of interpretability is often ambiguous or under-specified, but [42] presents
several definitions of interpretability to guide our selection of machine
learning methods. We argue that for a machine learning method to be
interpretable (1) it should be applicable equivalently to data and simulation,
(2) the method’s output can be understood outside the range of the training
set, and (3) a measurement uncertainty can be calculated. We argue that an
uncertainty on the method is not a proxy for a measurement uncertainty. These
stricter criteria are consistent with those outlined in [16]. Symbolic
regression satisfies these requirements because the output is a formula. The
convergence of the empirically-based multiplicity method and the formula
produced through symbolic regression is a clear indication of the usefulness
of an interpretable method. Machine learning should be used to gain knowledge
about the underlying physical processes that drive the relationships in our
data. We must interpret the details of any method in terms of these underlying
physical processes.

###### Acknowledgements.

We are grateful to Friederike Bock, Hannah Bossi, Adrian Del Maestro, Jamie
Nagle, Ken Read, and Austin Schmier for useful discussions and feedback on the
manuscript. This work was supported in part by funding from the Division of
Nuclear Physics of the U.S. Department of Energy under Grant No. DE-
FG02-96ER40982. This work was performed on the computational resources at the
Infrastructure for Scientific Applications and Advanced Computing (ISAAC)
supported by the University of Tennessee.

## References

  * Adcox  _et al._ [2005] K. Adcox _et al._ (PHENIX), Formation of dense partonic matter in relativistic nucleus nucleus collisions at RHIC: Experimental evaluation by the PHENIX collaboration, Nucl. Phys. A757, 184 (2005), arXiv:nucl-ex/0410003 .
  * Adams  _et al._ [2005] J. Adams _et al._ (STAR), Experimental and theoretical challenges in the search for the quark gluon plasma: The STAR collaboration’s critical assessment of the evidence from RHIC collisions, Nucl. Phys. A757, 102 (2005), arXiv:nucl-ex/0501009 .
  * Back  _et al._ [2005] B. B. Back _et al._ , The PHOBOS perspective on discoveries at RHIC, Nucl. Phys. A757, 28 (2005), arXiv:nucl-ex/0410022 .
  * Arsene  _et al._ [2005] I. Arsene _et al._ (BRAHMS), Quark Gluon Plasma an Color Glass Condensate at RHIC? The perspective from the BRAHMS experiment, Nucl. Phys. A757, 1 (2005), arXiv:nucl-ex/0410020 .
  * Connors  _et al._ [2018] M. Connors, C. Nattrass, R. Reed, and S. Salur, Jet measurements in heavy ion physics, Rev. Mod. Phys. 90, 025005 (2018), arXiv:1705.01974 [nucl-ex] .
  * Burke  _et al._ [2014] K. M. Burke _et al._ (JET), Extracting the jet transport coefficient from jet quenching in high-energy heavy-ion collisions, Phys.Rev. C90, 014909 (2014), arXiv:1312.5003 [nucl-th] .
  * Cao  _et al._ [2021] S. Cao _et al._ (JETSCAPE), Determining the jet transport coefficient q from inclusive hadron suppression measurements using Bayesian parameter estimation, Phys. Rev. C 104, 024905 (2021), arXiv:2102.11337 [nucl-th] .
  * Hughes  _et al._ [2022] C. Hughes, A. C. Oliveira Da Silva, and C. Nattrass, Model studies of fluctuations in the background for jets in heavy ion collisions, Phys. Rev. C 106, 044915 (2022), arXiv:2005.02320 [hep-ph] .
  * Adamczyk  _et al._ [2014] L. Adamczyk _et al._ (STAR), Jet-Hadron Correlations in $\sqrt{s_{NN}}=200$ GeV $p+p$ and Central $Au+Au$ Collisions, Phys. Rev. Lett. 112, 122301 (2014), arXiv:1302.6184 [nucl-ex] .
  * Abelev  _et al._ [2012] B. Abelev _et al._ (ALICE), Measurement of Event Background Fluctuations for Charged Particle Jet Reconstruction in Pb-Pb collisions at $\sqrt{s_{NN}}=2.76$ TeV, JHEP 03, 053, arXiv:1201.2423 [hep-ex] .
  * Soyez [2009] G. Soyez, Jet areas as a tool for background subtraction, in  _25th Winter Workshop on Nuclear Dynamics_ (2009) arXiv:0905.2851 [hep-ph] .
  * Abelev  _et al._ [2014] B. Abelev _et al._ (ALICE), Measurement of charged jet suppression in Pb-Pb collisions at $\sqrt{s_{NN}}$ = 2.76 TeV, JHEP 03, 013, arXiv:1311.0633 [nucl-ex] .
  * Adam  _et al._ [2015] J. Adam _et al._ (ALICE), Measurement of jet suppression in central Pb-Pb collisions at $\sqrt{s_{\rm NN}}$ = 2.76 TeV, Phys. Lett. B 746, 1 (2015), arXiv:1502.01689 [nucl-ex] .
  * Acharya  _et al._ [2020] S. Acharya _et al._ (ALICE), Measurements of inclusive jet spectra in pp and central Pb-Pb collisions at $\sqrt{s_{\rm{NN}}}$ = 5.02 TeV, Phys. Rev. C 101, 034911 (2020), arXiv:1909.09718 [nucl-ex] .
  * Adam  _et al._ [2020] J. Adam _et al._ (STAR), Measurement of inclusive charged-particle jet production in Au + Au collisions at $\sqrt{s_{NN}}=$200 GeV, Phys. Rev. C 102, 054913 (2020), arXiv:2006.00582 [nucl-ex] .
  * Achenbach  _et al._ [2023] P. Achenbach _et al._ , The Present and Future of QCD (2023) arXiv:2303.02579 [hep-ph] .
  * Haake and Loizides [2019] R. Haake and C. Loizides, Machine Learning based jet momentum reconstruction in heavy-ion collisions, Phys. Rev. C 99, 064904 (2019), arXiv:1810.06324 [nucl-ex] .
  * ALI [2023] Measurement of the radius dependence of charged-particle jet suppression in Pb-Pb collisions at $\sqrt{s_{\rm NN}}$ = 5.02 TeV, (2023), arXiv:2303.00592 [nucl-ex] .
  * Tannenbaum [2001] M. Tannenbaum, The distribution function of the event-by-event average pt for statistically independent emission, Physics Letters B 498, 29 (2001).
  * Sjostrand  _et al._ [2008] T. Sjostrand, S. Mrenna, and P. Z. Skands, A Brief Introduction to PYTHIA 8.1, Comput. Phys. Commun. 178, 852 (2008), arXiv:0710.3820 [hep-ph] .
  * Hughes and Mengel [2019] C. Hughes and T. Mengel, Github: TennGen (2019).
  * Mengel and Hughes [2022] T. Mengel and C. Hughes, Github: tenngen200 (2022).
  * Aamodt  _et al._ [2010] K. Aamodt _et al._ (ALICE Collaboration), Centrality dependence of the charged-particle multiplicity density at mid-rapidity in Pb-Pb collisions at $\sqrt{s_{NN}}$ = 2.76 TeV, Phys. Rev. Lett. 106, 032301. 14 p (2010).
  * Adare  _et al._ [2013] A. Adare _et al._ (PHENIX), Spectra and ratios of identified particles in Au+Au and $d$+Au collisions at $\sqrt{s_{NN}}=200$ GeV, Phys. Rev. C 88, 024906 (2013), arXiv:1304.3410 [nucl-ex] .
  * Abelev  _et al._ [2013] B. Abelev _et al._ (ALICE), Centrality dependence of $\pi$, K, p production in Pb-Pb collisions at $\sqrt{s_{NN}}$ = 2.76 TeV, Phys. Rev. C88, 044910 (2013), arXiv:1303.0737 [hep-ex] .
  * Adare  _et al._ [2016] A. Adare _et al._ (PHENIX), Measurement of the higher-order anisotropic flow coefficients for identified hadrons in Au$+$Au collisions at $\sqrt{s_{{}_{NN}}}$ = 200 GeV, Phys. Rev. C 93, 051902 (2016), arXiv:1412.1038 [nucl-ex] .
  * Adam  _et al._ [2016] J. Adam _et al._ (ALICE), Higher harmonic flow coefficients of identified hadrons in Pb-Pb collisions at $\sqrt{s_{\rm NN}}$ = 2.76 TeV, JHEP 09, 164, arXiv:1606.06057 [nucl-ex] .
  * Skands  _et al._ [2014] P. Skands, S. Carrazza, and J. Rojo, Tuning PYTHIA 8.1: the Monash 2013 Tune, Eur. Phys. J. C 74, 3024 (2014), arXiv:1404.5630 [hep-ph] .
  * Cacciari  _et al._ [2012] M. Cacciari, G. P. Salam, and G. Soyez, FastJet User Manual, Eur. Phys. J. C 72, 1896 (2012), arXiv:1111.6097 [hep-ph] .
  * Cacciari  _et al._ [2008] M. Cacciari, G. P. Salam, and G. Soyez, The Catchment Area of Jets, JHEP 04, 005, arXiv:0802.1188 [hep-ph] .
  * Ellis and Soper [1993] S. D. Ellis and D. E. Soper, Successive combination jet algorithm for hadron collisions, Phys. Rev. D 48, 3160 (1993), arXiv:hep-ph/9305266 .
  * Hanks  _et al._ [2012] J. A. Hanks _et al._ , Method for separating jets and the underlying event in heavy ion collisions at the BNL Relativistic Heavy Ion Collider, Phys. Rev. C 86, 024908 (2012), arXiv:1203.1353 [nucl-ex] .
  * Khachatryan  _et al._ [2017] V. Khachatryan _et al._ (CMS), Measurement of inclusive jet cross sections in $pp$ and PbPb collisions at $\sqrt{s_{NN}}=$ 2.76 TeV, Phys. Rev. C 96, 015202 (2017), arXiv:1609.05383 [nucl-ex] .
  * Aad  _et al._ [2013] G. Aad _et al._ (ATLAS), Measurement of the jet radius and transverse momentum dependence of inclusive jet suppression in lead-lead collisions at $\sqrt{s_{NN}}$= 2.76 TeV with the ATLAS detector, Phys. Lett. B 719, 220 (2013), arXiv:1208.1967 [hep-ex] .
  * Abelev  _et al._ [2015] B. B. Abelev _et al._ (ALICE), Charged jet cross sections and properties in proton-proton collisions at $\sqrt{s}=7$ TeV, Phys. Rev. D 91, 112012 (2015), arXiv:1411.4969 [nucl-ex] .
  * Chatrchyan  _et al._ [2014] S. Chatrchyan _et al._ (CMS), Measurement of Jet Fragmentation in PbPb and pp Collisions at $\sqrt{s_{NN}}=2.76$ TeV, Phys. Rev. C 90, 024908 (2014), arXiv:1406.0932 [nucl-ex] .
  * Abadi  _et al._ [2015] M. Abadi, A. Agarwal, P. Barham, E. Brevdo, Z. Chen, C. Citro, G. S. Corrado, A. Davis, J. Dean, M. Devin, S. Ghemawat, I. Goodfellow, A. Harp, G. Irving, M. Isard, Y. Jia, R. Jozefowicz, L. Kaiser, M. Kudlur, J. Levenberg, D. Mané, R. Monga, S. Moore, D. Murray, C. Olah, M. Schuster, J. Shlens, B. Steiner, I. Sutskever, K. Talwar, P. Tucker, V. Vanhoucke, V. Vasudevan, F. Viégas, O. Vinyals, P. Warden, M. Wattenberg, M. Wicke, Y. Yu, and X. Zheng, TensorFlow: Large-scale machine learning on heterogeneous systems (2015), software available from tensorflow.org.
  * Goodfellow  _et al._ [2016] I. J. Goodfellow, Y. Bengio, and A. Courville,  _Deep Learning_ (MIT Press, Cambridge, MA, USA, 2016) http://www.deeplearningbook.org.
  * Cranmer [2020] M. Cranmer, Pysr: Fast & parallelized symbolic regression in python/julia (2020).
  * Cranmer  _et al._ [2020] M. Cranmer, A. Sanchez-Gonzalez, P. Battaglia, R. Xu, K. Cranmer, D. Spergel, and S. Ho, Discovering symbolic models from deep learning with inductive biases, NeurIPS 2020 (2020), arXiv:2006.11287 [cs.LG] .
  * Steffanic  _et al._ [2023] P. Steffanic, C. Hughes, and C. Nattrass, Separating signal from combinatorial jets in a high background environment, (2023), arXiv:2301.09148 [hep-ph] .
  * Lipton [2018] Z. C. Lipton, The mythos of model interpretability: In machine learning, the concept of interpretability is both important and slippery., Queue 16, 31–57 (2018).