aainstitutetext: Theory Center, IPNS, KEK, 1-1 Oho, Tsukuba, Ibaraki 305-0801,
Japanbbinstitutetext: The Graduate University of Advanced Studies (Sokendai),
1-1 Oho, Tsukuba, Ibaraki 305-0801, Japanccinstitutetext: Kavli IPMU (WPI),
University of Tokyo, 5-1-5 Kashiwanoha, Kashiwa, Chiba 277-8583, Japan

# Interpretable Deep Learning for Two-Prong Jet Classification with Jet
Spectra

Amit Chakraborty a Sung Hak Lim a,b,c and Mihoko M. Nojiri<EMAIL_ADDRESS><EMAIL_ADDRESS><EMAIL_ADDRESS>
###### Abstract

Classification of jets with deep learning has gained significant attention in
recent times. However, the performance of deep neural networks is often
achieved at the cost of interpretability. Here we propose an interpretable
network trained on the jet spectrum $S_{2}(R)$ which is a two-point
correlation function of the jet constituents. The spectrum can be derived from
a functional Taylor series of an arbitrary jet classifier function of energy
flows. An interpretable network can be obtained by truncating the series. The
intermediate feature of the network is an infrared and collinear safe
C-correlator which allows us to estimate the importance of an $S_{2}(R)$
deposit at an angular scale $R$ in the classification. The performance of the
architecture is comparable to that of a convolutional neural network (CNN)
trained on jet images, although the number of inputs and complexity of the
architecture is significantly simpler than the CNN classifier. We consider two
examples: one is the classification of two-prong jets which differ in color
charge of the mother particle, and the other is a comparison between Pythia 8
and Herwig 7 generated jets.

###### Keywords:

Jets, QCD Phenomenology

††preprint: KEK-TH-2117††arxiv: 1904.02092

##  1 Introduction

Deep learning is gaining significant interest recently in the field of
collider data analysis. One of the primary motivations is to extract the
maximum information from the complex collision events. The deep learning in
collider physics takes advantage of a large influx of data from experiments,
more precise theoretical predictions, significant improvement in computing
power, and ongoing progress in the field of machine learning itself. Such
techniques offer advances in areas ranging from event selection to particle
identification.

The large center-of-mass energy at the Large Hadron Collider (LHC) enables the
production of boosted particles whose decay products are highly collimated.
These collimated objects are reconstructed as a jet, and it is often
misidentified as a QCD jet originated from light quarks or gluons. Many jet
substructure techniques using the information of subjets Butterworth:2008iy ;
Thaler:2008ju ; Kaplan:2008ie ; Plehn:2009rk ; Plehn:2010st ; Soper:2011cr ;
Soper:2012pb ; Dasgupta:2013ihk ; Soper:2014rya ; Larkoski:2014wba  and the
distribution of jet constituents Gallicchio:2010sw ; Thaler:2010tr ;
Gallicchio:2011xq ; Chien:2013kca ; Larkoski:2013eya ; Larkoski:2014gra ;
Moult:2016cvt  have been developed in order to improve the sensitivity of
tagging and to classify these boosted particle jets. The deep learning methods
Almeida:2015jua ; deOliveira:2015xxd ; Komiske:2016rsd ; Kasieczka:2017nvn ;
Louppe:2017ipp ; Komiske:2017ubm ; Butter:2017cot ; Cheng:2017rdo ;
Andreassen:2018apy ; Choi:2018dag ; Lim:2018toa ; Dreyer:2018nbf ; Lin:2018cin
; Komiske:2018cqr ; Martinez:2018fwc ; Kasieczka:2018lwf ; Qu:2019gqs  have
provided useful insight into the internal structure of the jets and, thereby,
shown better performances than those jet substructure techniques.111For a
review on the recent theoretical and machine learning developments in jet
substructure techniques at the LHC, we refer Larkoski:2017jix ;
Asquith:2018igt . The flexibility of deep learning also enables us to solve
problems beyond supervised classifications, such as weakly supervised learning
Dery:2017fap ; Cohen:2017exh ; Metodiev:2017vrx , adversarial learning to
suppress learning from unwanted information Louppe:2016ylz ; Shimmin:2017mfk ,
and unsupervised learning for finding anomalous signatures Chakraborty:2017mbz
; Hajer:2018kqm ; Heimel:2018mkt ; Farina:2018fyg ; Cerri:2018anq ;
Roy:2019jae . The neural network can also be useful to new physics searches
with deep learning at the LHC Roxlo:2018adx ; Brehmer:2018kdj ;
Brehmer:2018eca ; Guo:2018hbv ; Collins:2018epr ; DAgnolo:2018cun ;
DeSimone:2018efk ; Englert:2019xhk ; Collins:2019jip .

The output of a neural network is, in general, a highly non-linear function of
the inputs. A neural network classifier often acts like a “black box.” One may
consider architectures with post-hoc interpretability
DBLP:journals/corr/Lipton16a , which allows us to extract information other
than its prediction from the learned model after training. A simple strategy
is using a predefined functional form to restrict the representation power of
the neural network Komiske:2018cqr ; Datta:2019ndh . Then the network is
interpreted in terms of the functional form. The aim of this paper is also to
construct an interpretable neural network architecture that allows us not only
to interpret the predictions of the network but also to visualize it in terms
of trained weights connected to physical variables.

In Lim:2018toa , a multilayer perceptron (MLP) trained on two-point
correlation functions $S_{2}$ and $S_{2,\mathrm{trim}}$ of angular scale $R$
was introduced. The $S_{2}(R)$ and $S_{2,\mathrm{trim}}(R)$ spectra are
constructed from the constituents of a jet before and after the trimming
Krohn:2009th  respectively. The angular scale $R$ is an important parameter
for describing the kinematics of a decaying particle and parton shower (PS);
hence, these spectra efficiently encode the radiation pattern inside a jet.
The MLP trained on these inputs learns relevant features for the
classification among the Higgs boson jet (Higgs jet) and QCD jet.

In this paper, we connect the spectra to energy flow functionals
$P_{T}(\vec{R})$ Tkachov:1995kk , i.e., we consider transverse energy of a jet
constituent as particle-specific information at $\vec{R}$ in the $\eta-\phi$
plane 7974879 . The spectra are basis vectors of infrared and collinear (IRC)
safe variables called bilinear $C$-correlators Tkachov:1995kk  whose angular
weighting function depends only on the relative distance between two
constituents. Those correlators naturally appear in the functional Taylor
series of a classifier of $P_{T}(\vec{R})$, and the MLP can be considered as a
subseries of the Taylor series. We show that the performance of the MLP and
neural networks trained on jet images Cogan:2014oua ; Almeida:2015jua ;
deOliveira:2015xxd ; Kasieczka:2017nvn  are comparable. This strongly suggests
that $S_{2}$ and $S_{2,\mathrm{trim}}$ contain sufficient information for jet
classification. Encouraged by this feature, we construct an interpretable
architecture by truncating the series. Namely, $\int
dR\,S_{2}(R)w(R;\vec{x}_{\mathrm{kin}})$ can be implemented in a classifier
after proper discretization in $R$, where $\vec{x}_{\mathrm{kin}}$ is a set of
kinematic variables of the jet and $w$ is a smooth function. By reading the
weights $w(R;\vec{x}_{\mathrm{kin}})$, we could quantify important features
for the given classification problem.

Jet substructure studies often suffer from systematic uncertainties of soft
activities. The soft radiations generated by a Monte Carlo program are
strongly model dependent. While this mismodeling could be corrected by using
real data, it is certainly useful to use input variables with less systematic
uncertainties. When hard substructures are important for solving the problem,
we may use jet grooming techniques Butterworth:2008iy ; Krohn:2009th ;
Ellis:2009su ; Ellis:2009me ; Larkoski:2014wba  to remove the soft activity.
Instead of throwing this soft activity away, we encode it in
$S_{2,\mathrm{soft}}(R)$, which is $S_{2}(R)-S_{2,\mathrm{trim}}(R)$. Then,
the inputs $S_{2,\mathrm{trim}}$ and $S_{2,\mathrm{soft}}$ include hard and
soft substructure information, respectively. The interpretable architecture
trained on $S_{2,\mathrm{trim}}$ and $S_{2,\mathrm{soft}}$ is able to quantify
these features. We study two classification problems: one is a classification
of two-prong jets to understand their hard substructures and color coherence,
and the other is a comparison of Pythia 8 Sjostrand:2014zea  and Herwig 7
Bellm:2015jjp ; Bahr:2008pv  events to quantify the differences.

The paper is organized as follows. In section 2, we review $S_{2}$ and
$S_{2,\mathrm{trim}}$ and show its relation to energy flow and
$C$-correlators. We also show $S_{2}$ and $S_{2,\mathrm{trim}}$ distributions
of typical Higgs jet and QCD jet. A hypothetical color octet scalar particle,
sgluon, decaying to $b\bar{b}$ is considered to study the color connection in
two-prong jets. In section 3, we first discuss the capability of $S_{2}$ and
$S_{2,\mathrm{trim}}$ for the classification of two-prong jets and show the
result of an MLP trained on those inputs. The results are then compared with
that of a CNN trained on jet images. In section 4, we introduce a two-level
architecture consists of a softmax classifier and an MLP trained on $S_{2}$
and $S_{2,\mathrm{trim}}$. The intermediate feature of this architecture is
the bilinear $C$-correlator whose basis vectors are $S_{2,\mathrm{trim}}$ and
$S_{2,\mathrm{soft}}$, and the MLP generates its components. We visualize and
interpret the weights of the given classification problem. Finally, the
summary and outlook are given in section 5.

##  2 Two-Point Correlation Spectrum and Two-Prong Jets

###  2.1 Jet Spectra

In Lim:2018toa , we introduced a two-point correlation spectral function
$S_{2}(R)$ which maps a jet to a function of angular scale $R$,

$\displaystyle S_{2}(R)$ $\displaystyle=$ $\displaystyle\int
d\vec{R}_{1}\,d\vec{R}_{2}\,P_{T}(\vec{R}_{1})\,P_{T}(\vec{R}_{2})\,\delta(R-R_{12}),$
(1) $\displaystyle P_{T}(\vec{R})$ $\displaystyle=$
$\displaystyle\sum_{\begin{subarray}{c}i\in\mathbf{J}\end{subarray}}p_{T,i}\,\delta(\vec{R}-\vec{R}_{i}),$
(2)

where $\mathbf{J}$ is a set of jet constituents,
$\vec{R}_{i}=(\eta_{i},\phi_{i})$ is the position of the $i$-th jet
constituent in the pseudorapidity-azimuth plane,
$R_{ij}=\sqrt{(\eta_{i}-\eta_{j})^{2}+(\phi_{i}-\phi_{j})^{2}}$ is the angular
distance between the two jet constituents $i$ and $j$, and $P_{T}(\vec{R})$ is
an energy flow functional Tkachov:1995kk  of $\mathbf{J}$. For practical
purpose, $S_{2}(R)$ is discretized as below,

$\displaystyle S_{2}(R;\Delta R)$ $\displaystyle=$
$\displaystyle\frac{1}{\Delta R}\int_{R}^{R+\Delta R}dR\,S_{2}(R)$ (3)
$\displaystyle=$ $\displaystyle\frac{1}{\Delta
R}\sum_{i,j\in\mathbf{J}}p_{T,i}\,p_{T,j}\,I_{[R,R+\Delta R)}(R_{ij}),$

where $I_{A}(R_{ij})$ is an indicator function of the angular distance
$R_{ij}$ of the domain $A$,

$\displaystyle I_{A}(x)$ $\displaystyle=$
$\displaystyle\begin{cases}1&\mathrm{if}\;x\in A,\\\ 0&\mathrm{if}\;x\notin
A.\end{cases}$

The spectral function $S_{2}(R;\Delta R)$ is, therefore, the sum of the
product of $p_{T}$’s of the two jet constituents with an angular distance
$R_{ij}$ lying between $R$ and $R+\Delta R$.

We obtain IRC safe quantities by multiplying smooth functions222 Continuous
functions are sufficient for the convergence and IRC safety Tkachov:1995kk ,
but we further restrict $w$’s to smooth functions for perturbative
calculations.  $w(\vec{R})$ and $P_{T}(\vec{R})$ (or $S_{2}(R)$), and
integrating over $\vec{R}$. To understand the IRC safety of $P_{T}(\vec{R})$,
let us consider splitting of a given constituent $i_{0}$ in $\mathbf{J}$ into
two constituents, $i_{0}\rightarrow i_{1}i_{2}$. The inner product of
$w(\vec{R})$ and the difference of the energy flow before and after the
splitting, $\delta P_{T}(\vec{R})$, is given as follows,

$\displaystyle\int d\vec{R}\,\delta P_{T}(\vec{R})w(\vec{R})$ $\displaystyle=$
$\displaystyle
p_{T,i_{1}}w(\vec{R}_{i_{1}})+p_{T,i_{2}}w(\vec{R}_{i_{2}})-p_{T,i_{0}}w(\vec{R}_{i_{0}})$
(4) $\displaystyle=$ $\displaystyle\left[\delta
p_{T,i_{0}}-p_{T,i_{1}}(\delta\vec{R}_{i_{1}}\cdot\nabla_{\vec{R}})-p_{T,i_{2}}(\delta\vec{R}_{i_{2}}\cdot\nabla_{\vec{R}})+\cdots\right]w(\vec{R}_{i_{0}}),$

where $\delta p_{T,i_{0}}=p_{T,i_{1}}+p_{T,i_{2}}-p_{T,i_{0}}$, and
$\delta\vec{R}_{i_{1}(i_{2})}=\vec{R}_{i_{1}(i_{2})}-\vec{R}_{i_{0}}$. The
soft limit, where $i_{2}$ carries a small momentum, corresponds to $\delta
p_{T,i_{0}},\,\delta\vec{R}_{i_{1}},\;p_{T,i_{2}}\rightarrow 0$, while $\delta
p_{T,i_{0}},\,\delta\vec{R}_{i_{1}},\;\delta\vec{R}_{i_{2}}\rightarrow 0$ in
the collinear limit. The integral vanishes in these limits, namely the energy
flow after parton splitting converges weakly Tkachov:1995kk  to the one before
splitting.

The spectrum $S_{2}(R)$ inherits the same property. The inner product of the
smooth function $w(R)$ and the difference of the spectrum, $\delta S_{2}(R)$,
before and after the splitting $i_{0}\rightarrow i_{1}i_{2}$ is given as
follows,

$\int dR\,\delta S_{2}(R)w(R)=2\sum_{j\in\mathbf{J}}\left[\delta
p_{T,i_{0}}+p_{T,i_{1}}(\delta\vec{R}_{i_{1}}\cdot\nabla_{\vec{R}})+p_{T,i_{2}}(\delta\vec{R}_{i_{2}}\cdot\nabla_{\vec{R}})+\cdots\right]p_{T,j}\,w(R_{i_{0}j}).$
(5)

Again, this integral vanishes in the IRC limits. Note that the binned spectrum
$S_{2}(R;\Delta R)$ is not completely IRC safe because of the discontinuity of
the indicator function at the bin boundaries. Nevertheless, when the domain is
discretized into small sections $[R_{i},R_{i}+\Delta R_{i})$, the IRC unsafe
terms cancel in the sum, $\sum_{i}S_{2}(R_{i};\Delta R_{i})\,w(R_{i})$, and it
is approximately IRC safe up to binning errors.

The resulting IRC safe observables belong to $C$-correlators Tkachov:1995kk ,
which are multilinear forms of the energy flow. An $n$-linear $C$-correlator
is expressed as follows,

$\int d\vec{R}_{1}\cdots d\vec{R}_{n}\,P_{T}(\vec{R}_{1})\cdots
P_{T}(\vec{R}_{n})\,w(\vec{R}_{1},\cdots,\vec{R}_{n}),$ (6)

where $w$ is a continuous function of $\vec{R}_{1},\cdots,\vec{R}_{n}$. For
example, an inner product of $P_{T}(\vec{R})$ and $w(\vec{R})$ is a linear
$C$-correlator, and an inner product of $S_{2}(R)$ and $w(R)$ is a bilinear
$C$-correlator with $w$ depending only on the relative distance $R_{12}$,

$\int dR\,S_{2}(R)\,w(R)=\int
d\vec{R}_{1}d\vec{R}_{2}\,P_{T}(\vec{R}_{1})P_{T}(\vec{R}_{2})\,w(R_{12}).$
(7)

Many well-known jet observables belong to the $C$-correlator, for example, a
jet transverse momentum $p_{T,\mathbf{J}}$ is a linear $C$-correlator with
$w(\vec{R}_{1})\approx 1$, a jet mass $m_{\mathbf{J}}$ is a bilinear
$C$-correlator with $w(\vec{R}_{1},\vec{R}_{2})\approx R_{12}^{2}/2$.

The $S_{2}(R)$ spectra use all the jet constituents, but it is useful to
separate the correlations of constituents of the hard subjets; we consider jet
trimming for this purpose. We recluster the constituents of a jet of a radius
parameter $R_{\mathbf{J}}$ to subjets with a smaller radius parameter
$R_{\mathrm{trim}}$. A subjet $\mathbf{J}_{a}$ is discarded if
$p_{T,\mathbf{J}_{a}}<f_{\mathrm{trim}}\,p_{T,\mathbf{J}}$, where
$p_{T,\mathbf{J}}$ and $p_{T,\mathbf{J}_{a}}$ are the transverse momenta of
the jet and $a$-th subjet respectively. The trimmed jet
$\mathbf{J}_{\mathrm{trim}}$ is defined as a union of the remaining subjets,

$\mathbf{J}_{\mathrm{trim}}=\bigcup_{\begin{subarray}{c}a\\\
\frac{p_{T,\mathbf{J}_{a}}}{p_{T,\mathbf{J}}}\geq
f_{\mathrm{trim}}\end{subarray}}\mathbf{J}_{a}~{}.$ (8)

The jet trimming is beneficial because it does not introduce additional
angular scale parameters other than $R_{\mathrm{trim}}$. The trimmed spectrum
is then calculated using the constituents of the trimmed jet. We denote it as
$S_{2,\mathrm{trim}}(R)$ and its binned version $S_{2,\mathrm{trim}}(R;\Delta
R)$, which are defined as follows:

$\displaystyle S_{2,\mathrm{trim}}(R)$ $\displaystyle=$ $\displaystyle\int
d\vec{R}_{1}\,d\vec{R}_{2}\,P_{T,\mathrm{trim}}(\vec{R}_{1})\,P_{T,\mathrm{trim}}(\vec{R}_{2})\cdot\delta(R-R_{12}),$
(9) $\displaystyle P_{T,\mathrm{trim}}(\vec{R})$ $\displaystyle=$
$\displaystyle\sum_{\begin{subarray}{c}i\in\mathbf{J}_{\mathrm{trim}}\end{subarray}}p_{T,i}\,\delta(\vec{R}-\vec{R}_{i}),$
(10) $\displaystyle S_{2,\mathrm{trim}}(R;\Delta R)$ $\displaystyle=$
$\displaystyle\frac{1}{\Delta
R}\sum_{i,j\in\mathbf{J}_{\mathrm{trim}}}p_{T,i}~{}p_{T,j}\cdot I_{[R,R+\Delta
R)}(R_{ij}),$ (11)

where $P_{T,\mathrm{trim}}(\vec{R})$ is the energy flow of
$\mathbf{J}_{\mathrm{trim}}$.

In the limit of the constituents of each subjet $\mathbf{J}_{a}$ are
localized, the energy flow and the jet spectrum can be approximated in terms
of the subjet momenta. The energy flow of such a jet is decomposed into a sum
of energy flows of all the subjets,

$P_{T}(\vec{R})=\sum_{a}P_{T,a}(\vec{R}),\quad
P_{T,a}(\vec{R})=\sum_{i\in\mathbf{J}_{a}}p_{T,i}\delta(\vec{R}-\vec{R}_{i}).$
(12)

The energy flow of each subjet converges weakly to
$p_{T,\mathbf{J}_{a}}\delta(\vec{R}-\vec{R}_{\mathbf{J}_{a}})$. The $S_{2}(R)$
spectrum can be approximated by the momenta of the subjets, i.e.,

$\displaystyle S_{2}(R;\Delta R)$ $\displaystyle\approx$
$\displaystyle\sum_{\begin{subarray}{c}a,b\\\
\mathbf{J}_{a},\mathbf{J}_{b}\subset\mathbf{J}\end{subarray}}p_{T,\mathbf{J}_{a}}\,p_{T,\mathbf{J}_{b}}\cdot
I_{[R,R+\Delta R)}(R_{ab}).$ (13)

The jet trimming also introduces a $p_{T}$ scale hierarchy among the subjets,
and so their pairwise contributions to $S_{2}(R;\Delta R)$ can be classified
by the scale. We define a quantity $S_{2,\mathrm{soft}}(R;\Delta R)$ where

$S_{2,\mathrm{soft}}(R;\Delta R)=S_{2}(R;\Delta
R)-S_{2,\mathrm{trim}}(R;\Delta R).$ (14)

In the r.h.s. of the above equation, the correlations among the constituents
of the hard subjets are canceled, and we have

$\displaystyle S_{2,\mathrm{trim}}(R;\Delta R)$ $\displaystyle=$
$\displaystyle p_{T,\mathbf{J}}^{2}\cdot\mathcal{O}\left[1\right],$ (15)
$\displaystyle S_{2,\mathrm{soft}}(R;\Delta R)$ $\displaystyle=$
$\displaystyle
p_{T,\mathbf{J}}^{2}\cdot\left(\mathcal{O}\left[f_{\mathrm{trim}}\right]+\mathcal{O}\left[f_{\mathrm{trim}}^{2}\right]\right).$
(16)

The dominant contributions to $S_{2,\mathrm{soft}}(R;\Delta R)$ (i.e., the
$\mathcal{O}\left[f_{\mathrm{trim}}\right]$ terms) come from the correlations
between a constituent in $\mathbf{J}_{\mathrm{trim}}$ and a constituent in
$\mathbf{J}-\mathbf{J}_{\mathrm{trim}}$. The subleading
$\mathcal{O}\left[f_{\mathrm{trim}}^{2}\right]$ terms denote the correlations
among the constituents in $\mathbf{J}-\mathbf{J}_{\mathrm{trim}}$.

###  2.2 Derivation of Classifiers based on Energy Flows and Jet Spectra

We discuss the relation between $S_{2}(R)$ and neural network classifiers
trained on the energy flow $P_{T}(\vec{R})$. A general softmax classifier that
solves $K$-class jet classification problem can be expressed as a functional
$\hat{\Psi}_{i}$ which maps the energy flow to real numbers $h_{i}$, i.e.,

$\displaystyle h_{i}$ $\displaystyle=$ $\displaystyle\hat{\Psi}_{i}[P_{T}]\,$
(17) $\displaystyle\hat{y}$ $\displaystyle=$
$\displaystyle\varphi_{\mathrm{softmax}}\left(\vec{z}\right),\quad
z_{k}=w_{ki}^{(\mathrm{out})}h_{i}+b_{k}^{(\mathrm{out})},\;\;k\in\\{1,\cdots,K\\},$
(18)

where $w^{\mathrm{(out)}}_{ki}$ and $b^{\mathrm{(out)}}_{k}$ are the weights
and biases of the output layer, and $\hat{y}$ is the prediction of the
classifier. Here the $\varphi_{\mathrm{softmax}}$ is the softmax function
whose $k$-th component is expressed as follows,

$\varphi_{\mathrm{softmax},k}(\vec{z})=\frac{e^{z_{k}}}{\sum^{K}_{k=1}e^{z_{k}}}.$
(19)

Many jet classifiers can be expressed in the form of eq. (17). For example, in
the cut-based analysis, $\hat{\Psi}_{i}$ is a jet substructure variable, such
as a ratio of $n$-subjettiness Thaler:2010tr , a ratio of energy correlation
functions Larkoski:2014gra ; Moult:2016cvt , etc. The deep neural network
classifiers, such as artificial neural network tagger Almeida:2015jua ,
convolutional neural network using pixelated jet images deOliveira:2015xxd ,
energy flow network Komiske:2018cqr , etc., are also described by eq. (17).
The neural networks that are introduced in section 3 and section 4 also belong
to this category.

The jet spectra $S_{2}$ and $S_{2,\mathrm{trim}}$ can be derived from eq. (17)
using a functional Taylor expansion. The energy flow is decomposed by trimming
as follows,

$P_{T,a}(\vec{R})=\begin{cases}P_{T,\mathrm{trim}}(\vec{R})&a=1,\\\
P_{T}(\vec{R})-P_{T,\mathrm{trim}}(\vec{R})&a=2.\end{cases}$ (20)

One can express $\hat{\Psi}_{i}[P_{T,a}]$ as a functional series at a
reference point $P_{T,a}(\vec{R})=0$,

$\displaystyle h_{i}$ $\displaystyle=$ $\displaystyle w^{(0)}_{i}+\int
d\vec{R}\;P_{T,a}(\vec{R})w^{(1)}_{i,a}(\vec{R})+\frac{1}{2!}\int
d\vec{R}_{1}d\vec{R}_{2}\;P_{T,a}(\vec{R}_{1})P_{T,b}(\vec{R}_{2})w^{(2)}_{i,ab}(\vec{R}_{1},\vec{R}_{2})+\cdots,\phantom{0000}$
(21)

where $w^{(n)}_{i,a_{1}\cdots a_{n}}(\vec{R}_{1},\cdots,\vec{R}_{n})$ is the
coefficient of $n$-th correlation function. The first order coefficient
$w_{i,a}^{(1)}$ can be chosen as a constant if we are not interested in
features depending on reference vectors, for example, jet axes, beam
directions, etc. The linear term in $P_{T}(\vec{R})$ of eq. (21) is related to
the jet momentum $p_{T,\mathbf{J}}$ and trimmed jet momentum
$p_{T,\mathbf{J},\mathrm{trim}}$ as follows,

$\int d\vec{R}\,P_{T,1}(\vec{R})\simeq
p_{T,\mathbf{J},\mathrm{trim}},\quad\int d\vec{R}\,P_{T,2}(\vec{R})\simeq
p_{T,\mathbf{J}}-p_{T,\mathbf{J},\mathrm{trim}}.$ (22)

The second order coefficient $w^{(2)}_{i,ab}$, the first non-trivial term of
the series expansion, is a function of the relative distance of $\vec{R}_{1}$
and $\vec{R}_{2}$. The basis vectors of $w^{(2)}_{i,ab}$ are two-point
correlation functions $S_{2,ab}(R)$,

$\displaystyle h_{i}$ $\displaystyle=$ $\displaystyle w_{i}^{(0)}+\int
d\vec{R}\,P_{T,a}(\vec{R})\,w_{i,a}^{(1)}+\frac{1}{2!}\int
dR\,S_{2,ab}(R)\,w_{i,ab}^{(2)}(R)+\cdots$ (23) $\displaystyle S_{2,ab}(R)$
$\displaystyle=$ $\displaystyle\int
d\vec{R}_{1}d\vec{R}_{2}\,P_{T,a}(\vec{R}_{1})P_{T,b}(\vec{R}_{2})\,\delta(R-R_{12}).$
(24)

The spectra $S_{2}$ and $S_{2,\mathrm{trim}}$ are expressed in terms of
$S_{2,ab}$ as follows,

$S_{2}(R)=\sum_{a,b}S_{2,ab}(R),\quad S_{2,\mathrm{trim}}(R)=S_{2,11}(R).$
(25)

Instead of the energy flows, we consider a classifier of $S_{2,A}$ ($A$ =
$\mathrm{trim}$, $\mathrm{soft}$),

$h_{i}=\Psi_{i}[S_{2,A};\vec{x}_{\mathrm{kin}}],$ (26)

where $\vec{x}_{\mathrm{kin}}$ is a set of additional inputs to the classifier
based on the kinematics of the jet. Similar to eq. (23), we expand eq. (26)
around $S_{2,A}(R)=0$ as

$\displaystyle h_{i}$ $\displaystyle=$ $\displaystyle
w_{i}^{(0)}(\vec{x}_{\mathrm{kin}})+\int
dR\,S_{2,A}(R)\,\frac{w_{i,A}^{(2)}(R;\vec{x}_{\mathrm{kin}})}{2}$ (27)
$\displaystyle\phantom{00}+\frac{1}{2}\int
dR_{1}dR_{2}\,S_{2,A_{1}}(R_{1})S_{2,A_{2}}(R_{2})\,\frac{w_{i,A_{1}A_{2}}^{(4)}(R_{1},R_{2};\vec{x}_{\mathrm{kin}})}{12}+\cdots,$

where $w_{i,A_{1}\cdots A_{\frac{n}{2}}}^{(n)}$ is the weight function
corresponding to $w_{i}^{(n)}$ in eq. (26). One may further truncate the
series to get a linear form,

$h_{i}=\frac{1}{2}\int dR\,S_{2,A}(R)w_{i,A}^{(2)}(R;\vec{x}_{\mathrm{kin}}).$
(28)

The above-mentioned linear setup has an advantage on the interpretability and
visualization of the network predictions; we discuss more on this network in
section 4.

###  2.3 Relation between Two-Point Correlation Spectra and Energy Flow
Polynomials

Both the two-point correlation spectra and the energy flow polynomials
Komiske:2017aww  with two vertices span the set of bilinear $C$-correlators;
therefore, there is a transformation rule between them. We first extend the
definition of the energy flow polynomials to compare them to $S_{2,ab}$. Since
$S_{2,ab}$ is a multivariate function of energy flows, we introduce a
multivariate energy flow polynomial with two labeled vertices,

$\displaystyle{\mathrm{EFP}}_{2,ab}^{(n)}$ $\displaystyle=$ $\displaystyle\int
d\vec{R}_{1}d\vec{R}_{2}\,P_{T,a}(\vec{R}_{1})P_{T,b}(\vec{R}_{1})\,R_{12}^{n}=\sum_{i\in{\mathbf{J}_{a}}}\sum_{j\in{\mathbf{J}_{b}}}p_{T,i}p_{T,j}R_{ij}^{n}.$
(29)

This expression suggests that $R^{n}$ can be considered as an angular
weighting function $w_{i,ab}^{(2)}(R)$ in eq. (23).

The resulting transformation from $S_{2,ab}(R)$ to
${\mathrm{EFP}}_{2,ab}^{(n)}$ is the Mellin transformation,

$\displaystyle{\mathrm{EFP}}_{2,ab}^{(n-1)}$ $\displaystyle=$
$\displaystyle\int_{0}^{\infty}dR\,R^{n-1}\cdot S_{2,ab}(R).$ (30)

The integral on the right-hand side is finite because $S_{2,ab}(R)$ vanishes
on $R\gg 2R_{\mathbf{J}}$. The inverse transform is also well-defined if we
allow the exponent $n$ in the angular weighting function of
${\mathrm{EFP}}_{2,ab}^{(n)}$ to be a complex number.

###  2.4 Spectra of Two-Prong Jets

The $S_{2}(R;\Delta R)$ spectrum is useful to identify the substructures of
the jet and also to characterize the jet. Typically, $S_{2}(R;\Delta R)$ of
QCD jet has a peak at $R=0$ with a long tail towards large $R$. The peak
originates from the autocorrelation term $\sum_{i}p^{2}_{T,i}$ in eq. (3). On
the other hand, if a jet originates from a Higgs boson decaying into
$b\bar{b}$, the $b$-partons create two isolated cores inside the jet. The
spectrum of the Higgs jet has a peak at the angular scale equal to the angle
between the two clusters. In addition, $S_{2}(R;\Delta R)$ encodes the
fragmentation pattern of $b$-partons.

At the LHC, boosted heavy objects such as top quark, gauge bosons and Higgs
boson decaying into quarks can be studied by identifying jet substructures.
Usually, these substructures are characterized by parameters such as $D_{2}$
defined as,

$\displaystyle D_{2}^{\beta}$ $\displaystyle=$ $\displaystyle
e_{3}^{\beta}/(e^{\beta}_{2})^{3},$ (31) $\displaystyle e_{2}^{\beta}$
$\displaystyle=$
$\displaystyle\frac{1}{p^{2}_{T,\mathbf{J}}}\sum_{i,j\in\mathbf{J},i<j}p_{T,i}p_{T,j}R^{\beta}_{ij},$
(32) $\displaystyle e_{3}^{\beta}$ $\displaystyle=$
$\displaystyle\frac{1}{p^{3}_{T,\mathbf{J}}}\sum_{i,j,k\in\mathbf{J},i<j<k}p_{T,i}p_{T,j}p_{T,k}R^{\beta}_{ij}R^{\beta}_{jk}R^{\beta}_{ki},$
(33)

where $\beta$ is the angular exponent. If a jet has a two-prong substructure,
$D_{2}$ is much less than one. The jet spectrum $S_{2}(R;\Delta R)$ contains
more information than $D_{2}$, and therefore, the analysis with
$S_{2}(R;\Delta R)$ goes beyond the one using $D_{2}$. It was shown that a
neural network trained on $S_{2}(R;\Delta R)$ distinguishes Higgs jet from QCD
jet better than the one trained on $D_{2}$ Lim:2018toa .

To study the fragmentation pattern of the $b$-partons and their color
connection to the mother particle, we introduce a color-octet scalar, sgluon
($\sigma$). We assume that the Higgs boson ($h$) and $\sigma$ decay into
$b\bar{b}$ through the interaction,

$\displaystyle\mathcal{L_{\mathrm{SM}}}$ $\displaystyle\ni$ $\displaystyle
y_{hb\bar{b}}\,h\,\bar{b}b+\mathrm{h.c.}$ (34)
$\displaystyle\mathcal{L_{\mathrm{Sgluon}}}$ $\displaystyle\ni$ $\displaystyle
y_{\sigma b\bar{b}}\,\sigma^{a}\,\bar{b}T^{a}b+\mathrm{h.c.}.$ (35)

The Higgs boson is a color singlet particle, and the decay $h\rightarrow
b\bar{b}$ is isolated in color flows. Therefore, $S_{2}(R;\Delta R)$ beyond
the angle between the b-partons, i.e., $R_{b\bar{b}}$, is suppressed due to
the color coherence. No such constraint on the angular scale exists for sgluon
and QCD jets. Meanwhile, the Higgs jet and sgluon jet have the same two-prong
substructure, unlike the QCD jet, as both are originating from a particle
decaying into $b\bar{b}$ final states.

To study the spectra of two-prong jets, we simulate events as follows. We use
Madgraph5 2.6.1 Alwall:2014hca  to generate the events of $pp\rightarrow Zh$,
$pp\rightarrow Z\sigma$, and $pp\rightarrow Zj$ processes with a collision
energy of 13 TeV and the $Z$ boson decaying to a pair of neutrinos. These
events are then passed to Pythia 8.226 Sjostrand:2014zea  for the parton
shower and hadronizations. To study the impact of the parton shower and
hadronization schemes, we also pass those parton level events to Herwig 7.1.3
Bellm:2015jjp ; Bahr:2008pv . A color octet scalar UFO model Degrande:2014sta
; NLOModels  generated by Feynrule 2.0 Alloul:2013bka ; Degrande:2011ua  is
used to simulate $pp\rightarrow Z\sigma$ process. The masses and widths of
Higgs boson and sgluon are $m_{h}=m_{\sigma}=$ 125 GeV and
$\Gamma_{\sigma}=\Gamma_{h}=6.4$ MeV. The detector response is simulated by
Delphes 3.4.1 deFavereau:2013fsa  with the default ATLAS detector
configuration. We use FastJet 3.3.0 Cacciari:2011ma ; Cacciari:2005hq  to
reconstruct jets from the calorimeter towers using anti-$k_{T}$ algorithm
Cacciari:2008gp  with the radius parameter $R_{\mathbf{J}}=1.0$. For jet
trimming, we use $R_{\mathrm{trim}}=0.2$ and $f_{\mathrm{trim}}=0.05$. We
select the events with the leading jet transverse momentum
$p_{T,\mathbf{J}}\in[300,400]$ GeV and its mass $m_{\mathbf{J}}\in[100,150]$
GeV. For Higgs jet and sgluon jet, we additionally require that the two
$b$-partons originating from their decay are located within $R_{\mathbf{J}}$
from the jet axis. More details on our simulations are described in appendix
A.

Figure 1:  A jet image (left) and the corresponding $S_{2}$ (blue) and
$S_{2,\mathrm{trim}}$ (red) spectra (right) of the leading jet of a
$pp\rightarrow Zj$ event. In the jet image, a red triangle is a position of a
parton level quark in the jet, and a green “+” shows a leading gluon emitted
from the quark.

Figure 2:  A jet image (left) and the corresponding $S_{2}$ (blue) and
$S_{2,\mathrm{trim}}$ (red) spectra (right) of the leading jet of a
$pp\rightarrow Zh$ event. The blue triangles in the jet image are positions of
the parton level bottom quarks from the Higgs decay, and a green “+” shows the
position of a leading gluon emitted from a bottom quark.

Figure 3:  A jet image (left) and the corresponding $S_{2}$ (blue) and
$S_{2,\mathrm{trim}}$ (red) spectra (right) of the leading jet of a
$pp\rightarrow Z\sigma$ event. The blue triangles in the jet image are
positions of the parton level bottom quarks from the sgluon decay, and green
“+” show the position of a leading gluon emitted from a bottom quark.

In figure 1, we show the pixelated jet image (left panel) and $S_{2}$ and
$S_{2,\mathrm{trim}}$ spectra (right panel) of a QCD jet. There are high
energy deposits in the jet image near the jet center along with a wide spray
of soft activity. It also has a moderate amount of radiation at $(-0.4,0.0)$.
As a result, $S_{2}(R;\Delta R)$ spectra has a long tail starting from
$R=0.4$. The jet trimming eliminates a significant amount of soft particles
and, therefore, the tail does not appear in $S_{2,\mathrm{trim}}(R;\Delta R)$.
The remaining cross-correlations contributing to $S_{2,\mathrm{trim}}(R;\Delta
R)$ are the ones between high and moderate energy deposits. Most of the energy
deposits are concentrated at the center, and the peak intensity at $R=0.4$ is
much lower than the intensity from autocorrelations at $R=0$.

In figure 2, we show $S_{2}(R;\Delta R)$ and $S_{2,\mathrm{trim}}(R;\Delta R)$
distributions of a Higgs jet. For this particular event, $S_{2}(R;\Delta R)$
distribution is similar to $S_{2,\mathrm{trim}}(R;\Delta R)$ distribution, and
their difference $S_{2,\mathrm{soft}}(R;\Delta R)$ is hard to be seen. No
significant activity has been observed beyond the peak at $R\sim 0.8$, mostly
because the Higgs jet is very compact compared to the QCD jet.
Correspondingly, there are two prominent subjets in the jet image, while most
of the cells have no jets.

Finally, we show the $S_{2}(R;\Delta R)$ and $S_{2,\mathrm{trim}}(R;\Delta R)$
distributions of a sgluon jet in figure 3. The $S_{2}(R;\Delta R)$
distribution has a large peak at $R=0.6$ which is as significant as the one at
$R=0$. This spectrum is qualitatively similar to the Higgs jet in figure 2.
However, the $S_{2}(R;\Delta R)$ spectrum has a long tail beyond
$R_{\mathbf{J}}$ as compared with that of a Higgs jet. The tail disappears
after jet trimming, like the QCD jet in figure 1, that makes the
$S_{2,\mathrm{trim}}(R;\Delta R)$ distribution more compact. From figure 1-3,
we observe that $S_{2,\mathrm{trim}}$ and $S_{2,\mathrm{soft}}$ include useful
complementary information.

In Lim:2018toa , it was shown that a neural network classifier trained on
$S_{2}(R;\Delta R)$ and $S_{2,\mathrm{trim}}(R;\Delta R)$ spectra performs
better than one without $S_{2,\mathrm{trim}}(R;\Delta R)$. The reason is that
the hard and soft correlation terms in $S_{2}(R)$, i.e., $\mathcal{O}[1]$
terms in eq. (15) and
$\mathcal{O}[f_{\mathrm{trim}}]+\mathcal{O}[f_{\mathrm{trim}}^{2}]$ in eq.
(16) respectively, can be resolved by the jet trimming. Therefore, we use the
orthogonal combinations, namely $S_{2,\mathrm{trim}}$ and
$S_{2,\mathrm{soft}}$, throughout this paper.

The $S_{2,\mathrm{trim}}$ and $S_{2,\mathrm{soft}}$ spectra encode the
important features of the parton shower and fragmentation, and, thus, may be
regarded as a well-motivated prototype. The hard partons evolve by the parton
splittings $i\rightarrow i_{1}i_{2}$, which are parameterized by the angle
$R_{i_{1}i_{2}}$ and momentum fraction $z$ with $p_{T,i_{1}}=zp_{T,i}$ and
$p_{T,i_{2}}=(1-z)p_{T,i}$. The splitting generates two-point correlation
$z(1-z)p_{T,i}^{2}$ at $R_{i_{1}i_{2}}$. Therefore, $S_{2}$ spectra encode the
parton splitting at any angular scale.

##  3 Classifying Higgs jet, sgluon jet, and QCD jet

In this section, we introduce a neural network trained on the jet spectra for
classifying Higgs jet, sgluon jet, and QCD jet. We first discuss the basic
kinematic features of these jets and then outline their dependence on the
parton shower simulators. Afterward, we show the details of the neural network
and then present our results in terms of the receiver operating characteristic
(ROC) curves.

###  3.1 Basic Kinematics

In figure 4, we show $p_{T,\mathbf{J}}$ and $m_{\mathbf{J}}$ distributions for
the Higgs boson, sgluon, and QCD jets. The solid and dashed lines correspond
to Pythia 8 (PY8) and Herwig 7 (HW7) generated jets, respectively. The mild
differences in the $p_{T}$ distribution are due to the difference in their
matrix elements. The Higgs jet is produced via $s$-channel process, while the
sgluon and QCD jet are produced via $t$-channel and $u$-channel processes;
hence, $p_{T,\mathbf{J}}$ scalings are different. Not much difference is
observed between the $p_{T,\mathbf{J}}$ distributions of PY8 and HW7 samples.
This is because $p_{T,\mathbf{J}}$ is mostly determined by the matrix level
$p_{T}$ of the leading parton and the jet algorithm with large radius
parameter clusters most of the radiations from this parton into a single jet.
However, the difference between $m_{\mathbf{J}}$ distributions is large. The
peak at $m_{\sigma}$ of sgluon jet is significantly broader than that of Higgs
jet because radiations of the b-partons from the Higgs boson decay are mostly
confined due to the color coherence, but those of the sgluon are not. As a
consequence, PY8 and HW7 generate different $m_{\mathbf{J}}$ distributions.

Figure 4: The distribution of $p_{T,\mathbf{J}}$ (left) and $m_{\mathbf{J}}$
(right) of the leading jet. The red, green, and blue solid (dashed) lines
correspond to the Higgs jet, sgluon jet, and QCD jet of PY8 (HW7) samples,
respectively.

We assume that both Higgs boson and sgluon have narrow-widths although sgluon
width can be large. An increase in the width will broaden $R_{b\bar{b}}$
distribution of $\sigma\rightarrow b\bar{b}$ that has a peak at the
characteristic angular scale,

$\hat{R}_{b\bar{b}}=\frac{2m_{h}}{p_{T,\mathbf{J}}}=\frac{2m_{\sigma}}{p_{T,\mathbf{J}}}\simeq\frac{2m_{\mathbf{J}}}{p_{T,{\mathbf{J}}}}.$
(36)

For example, the variation of $R_{b\bar{b}}$ is only about 0.07 for
$p_{T,\mathbf{J}}=300$ GeV, $\Gamma_{\sigma}=10$ GeV and 0.05 for
$p_{T,\mathbf{J}}=400$ GeV, $\Gamma_{\sigma}=10$ GeV. Those variations are
close to the calorimeter angular resolution $\sim$ 0.1 and do not affect the
calorimeter level analysis.

We first make a quantitative estimate of the radiation pattern inside the jet.
To do so, we define two quantities comparing $S_{2}(R)$ spectra around
$\hat{R}_{b\bar{b}}$,

$\displaystyle R_{\mathrm{sym}}$ $\displaystyle=$
$\displaystyle\frac{\int_{a\hat{R}_{b\bar{b}}}^{\min[a^{\prime}\hat{R}_{b\bar{b}},R_{\mathbf{J}}]}dR\,S_{2}(R)}{\int_{0}^{a\hat{R}_{b\bar{b}}}dR\,S_{2}(R)+\int_{\min[a^{\prime}\hat{R}_{b\bar{b}},R_{\mathbf{J}}]}^{\infty}dR\,S_{2}(R)},$
(37) $\displaystyle R_{\mathrm{rad}}$ $\displaystyle=$
$\displaystyle\frac{C\int_{\min[a^{\prime}\hat{R}_{b\bar{b}},R_{\mathbf{J}}]}^{\infty}dR\,S_{2}(R)}{\int_{0}^{\min[a^{\prime}\hat{R}_{b\bar{b}},R_{\mathbf{J}}]}dR\,S_{2}(R)+C\int_{\min[a^{\prime}\hat{R}_{b\bar{b}},R_{\mathbf{J}}]}^{\infty}dR\,S_{2}(R)}$
(38)

with $a=0.75$, $a^{\prime}=1.25$ and $C=40$. The ratio $R_{\mathrm{sym}}$
compares energy deposits around $\hat{R}_{b\bar{b}}$ and in its surrounding
angular scales Lim:2018toa . The ratio is sensitive to the correlation between
the two hard substructures of the Higgs jet. On the other hand, The
$R_{\mathrm{rad}}$ is sensitive to the color of mother particle as it compares
energy deposits in the large angular scales.

Figure 5: The distribution of $R_{\mathrm{prong,sym}}$ (left) and
$R_{\mathrm{rad}}$ (right) for Higgs jet (red), sgluon jet (green), and QCD
jet (blue). The solid (dashed) lines correspond to the jets of PY8 (HW7)
samples.  Figure 6:  The distribution of $R_{\mathrm{rad}}$ for Higgs jet
(red), sgluon jet (green), and QCD jet (blue) with an additional selection of
$R_{\mathrm{sym}}>0.85$. The solid (dashed) lines correspond to the jets of
PY8 (HW7) samples.

We show the $R_{\mathrm{sym}}$ distributions in the left panel of figure 5.
The distributions of the Higgs jet and sgluon jet are similar because both of
the $S_{2}(R)$ spectra peak at ${R}_{b\bar{b}}$. Meanwhile, the two-point
correlations for the QCD jet are not localized around the ${R}_{b\bar{b}}$
scale, so the $R_{\mathrm{sym}}$ is smaller than that of a Higgs jet and a
sgluon jet. In the right panel of figure 5, we show the $R_{\mathrm{rad}}$
distributions. The $R_{\mathrm{rad}}$ of the sgluon jet and QCD jet are large
on average, while $R_{\mathrm{rad}}$ is smaller for Higgs jet because large
angle radiations are suppressed.

The difference in $R_{\mathrm{sym}}$ and $R_{\mathrm{rad}}$ distributions
between PY8 and HW7 samples is small; however, there is an appreciable
difference in the restricted phase space. In figure 6, we plot
$R_{\mathrm{rad}}$ distributions after the selection, $R_{\mathrm{sym}}>0.85$,
so that the jets always contain two hard subjets with similar transverse
momenta. The PY8 (solid line) and HW7 (dashed line) samples have significantly
different $R_{\mathrm{rad}}$ distributions for the Higgs jet. Such a
difference is not observed for the QCD/sgluon jets. The observed deviation for
the Higgs jets could be originating from the difference of the parton shower
scheme. The angular-ordered shower is adopted in HW7. On the other hand, the
$p_{T}$-ordered shower is the default shower algorithm for PY8 where angular
ordering is enforced by hand. An artificial veto in $p_{T}$-ordered shower
introduces the mismatch to the angular-ordered shower at double-leading log
level Webber:2010vz ; Bhattacherjee:2015psa .

###  3.2 Multilayer Perceptron of Spectra

We introduce a neural network trained on the kinematic and spectrum
($S_{2,\mathrm{trim}}$ and $S_{2,\mathrm{soft}}$) variables to classify the
jets. A schematic diagram of the architecture of the classifier is shown in
figure 7. The following set of inputs is used,

$\displaystyle\vec{x}=\\{p_{T,\mathbf{J}},m_{\mathbf{J}},p_{T,\mathbf{J},\mathrm{trim}},m_{\mathbf{J},\mathrm{trim}}\\}\cup\left\\{S_{2,\mathrm{trim}}^{k},S_{2,\mathrm{soft}}^{k}\,|\,k\in\\{0,\cdots,19\\}\right\\},$
(39)

where $p_{T,\mathbf{J},\mathrm{trim}}$ and $m_{\mathbf{J},\mathrm{trim}}$ are
the transverse momentum and mass of the trimmed jet, respectively. The
discretized spectra $S_{2,\mathrm{trim}}^{k}$ and $S_{2,\mathrm{soft}}^{k}$
are used to analyze the radiation pattern of the jet,

$\displaystyle S_{2,\mathrm{trim}}^{k}$ $\displaystyle=$ $\displaystyle
S_{2,\mathrm{trim}}(0.1\,k;0.1),$ (40) $\displaystyle S_{2,\mathrm{soft}}^{k}$
$\displaystyle=$ $\displaystyle
S_{2}(0.1\,k;0.1)-S_{2,\mathrm{trim}}(0.1\,k;0.1).$ (41)

Here we take the bin width $\Delta R=0.1$, which is approximately the angular
resolution of hadronic calorimeter of the ATLAS detector. Note that the
maximum separation between any two constituents of the jet is
$2R_{\mathbf{J}}$.

jetspectralfunctions$S_{2,\mathrm{trim}}$$S_{2,\mathrm{soft}}$$p_{T,\mathbf{J}}$$m_{\mathbf{J}}$$p_{T,\mathbf{J},\mathrm{trim}}$$m_{\mathbf{J},\mathrm{trim}}$kinematicsMLPsoftmaxclassifier$\hat{y}_{1}$$\hat{y}_{2}$
Figure 7:  Schematic diagram of the classifier, including the multilayer
perceptron (MLP). The double bordered boxes represent trainable modules.

A multilayer perceptron (MLP) with $L$ layers is used to map the inputs to the
class prediction. The following first-order recurrence relation between the
layers describes an MLP,

$\displaystyle h^{(\ell)}_{i}$ $\displaystyle=$
$\displaystyle\varphi^{(\ell)}\left(w^{(\ell)}_{ij}h^{(\ell-1)}_{j}+b^{(\ell)}_{i}\right),\quad\vec{h}^{(0)}=\vec{x},$
(42)

where $w^{(\ell)}_{ij}$ and $b^{(\ell)}_{i}$ are the weight and bias of the
$\ell$-th layer. The activation function of the $\ell$-th layer,
$\varphi^{(\ell)}:\mathbb{R}\rightarrow\mathbb{R}$, is a monotonic and
nonlinear function. We use four hidden layers with 1000, 800, 400, and 200
nodes, respectively, with a rectified linear unit (ReLU),
$\varphi_{\mathrm{ReLU}}(x)=\max(0,x)$, as the activation function. This MLP
will identify important features of inputs for the classification after
training. To make a class prediction, we provide the outputs of the MLP to a
softmax classifier in eq. (18). The whole network architecture is illustrated
in figure 8.

layer:inputdensedensedensedensedense(output)activation:$\mathrm{ReLU}$$\mathrm{ReLU}$$\mathrm{ReLU}$$\mathrm{ReLU}$softmaxbias:000$\vdots$000$\vdots$000$\vdots$000$\vdots$000$\vdots$11111000$\hat{y}_{3}$$\hat{y}_{2}$$\hat{y}_{1}$$\vec{x}$
Figure 8: Schematic diagram of the multilayer perceptron.

The MLP is trained by minimizing a loss function including categorical cross-
entropy and $L_{2}$ weight regularization NIPS1991_563 ,

$\mathcal{L}=\frac{1}{N_{\mathrm{events}}}\sum_{\mathrm{events}}^{N_{\mathrm{events}}}\sum_{i}y_{i}\log\hat{y}_{i}\,+\,\lambda\sum_{\ell=1}^{L}\sum_{i,j}|w_{ij}^{(\ell)}|^{2},$
(43)

where $N_{\mathrm{events}}$ is the total number of events in the training data
set, $\lambda$ is a weight decay constant associated to the $L_{2}$
regularization. We choose $\lambda=0.01$. The $y_{i}$ ($\hat{y}_{i}$) denotes
the components of the truth (predicted) label vector $\vec{y}$ ($\hat{y}$).
The $L_{2}$ weight regularization reduces the over-fitting on the training
data and also allows smooth extrapolation to the phase space that is not
covered by the training sample. The minimization is done with ADAM optimizer
DBLP:journals/corr/KingmaB14 . We stop training when the validation loss has
stopped improving for 50 epochs. After the minimization of the loss function,
the softmax layer provides scores of the classes of a given event. The truth
label vectors are defined as follows,

$\displaystyle\vec{y}$ $\displaystyle=$
$\displaystyle\begin{cases}(1,0,0)&\textrm{Higgs jet},\\\
(0,1,0)&\textrm{sgluon jet},\\\ (0,0,1)&\textrm{QCD jet}.\end{cases}$ (44)

The unnecessary symmetries in the neural network are broken by using the
Glorot uniform initialization method pmlr-v9-glorot10a . The weights in the
hidden layers are initialized by assigning random numbers between
$[-\sqrt{6/(N_{\mathrm{in}}+N_{\mathrm{out}})},\sqrt{6/(N_{\mathrm{in}}+N_{\mathrm{out}})}]$,
where $N_{\mathrm{in}}$ and $N_{\mathrm{out}}$ are numbers of inputs and
outputs of a layer, respectively. The biases are initialized to zero. All the
inputs are standardized before training. The architecture is implemented in
Keras chollet2015keras  with backend TensorFlow tensorflow2015-whitepaper .

Figure 9:  Ternary plots of the predicted label vector $\hat{y}$ of the MLP
for the Higgs jet (left), sgluon jet (center), and QCD jet (right).

In figure 9, we show ternary plots of the predicted label vector $\hat{y}$.
The three sides of the triangle (starting from the base of the triangle and
then counterclockwise) are $\hat{y}_{1}$, $\hat{y}_{2}$, $\hat{y}_{3}$ axis;
we denote them as $\hat{y}_{h}$, $\hat{y}_{\sigma}$ and $\hat{y}_{\rm QCD}$,
respectively. The $\hat{y}$ distributions of the Higgs jet and QCD jet have
high-density spots that do not overlap with each other. It means that the
network has found the exclusive features of those two kinds of jets. The two-
prong substructure of a Higgs jet and the one-prong structure of a QCD jet are
the exclusive features. However, the two-prong substructure of a sgluon jet is
more radiative and less exclusive, and therefore, there are no high-density
spots in the $\hat{y}$ distribution of the sgluon jet.

Next, we show ROC curves of binary classifications in figure 10 with the red
dotted lines. The following signal-background classifications are considered:
Higgs-QCD, sgluon-QCD, and Higgs-sgluon. We assign the truth label vectors
$\vec{y}=(1,0)$ for the signal and $\vec{y}=(0,1)$ for the background. The QCD
jet mistag rates are comparable for both Higgs-QCD and sgluon-QCD
classifications; however, the separation between the Higgs jet and sgluon jet
is weaker.

We now compare these ROC curves with that of CNN trained on jet images.333The
CNN setup is explained in detail in appendix C. The CNN classifier takes
$20\times 20$ inputs of the jet images, while $2\times 20$ inputs of
$S_{2,\mathrm{trim}}$ and $S_{2,\mathrm{soft}}$ spectra are used for the MLP.
The solid blue lines in figure 10 denote the ROC curves of the CNN. Some
improvement in the background mistag rates is observed compared with the MLP
classifier. Quantitatively, it is only 0.2% ($=2.5\%-2.3\%$) at the signal
acceptance of 20% for Higgs-QCD classification.

Figure 10:  The ROC curves of the binary classifiers: the MLP trained on
$S_{2,\mathrm{trim}}$ and $S_{2,\mathrm{soft}}$ (red dashed), the CNN trained
on jet images (blue solid), and the two-level architecture (see section 4)
trained on $S_{2,\mathrm{trim}}$ and $S_{2,\mathrm{soft}}$ (green dotted) with
PY8 samples. The dashed gray lines represent the ROC curves of the random
guess. We show the results of Higgs jet vs. QCD jet (left), sgluon jet vs. QCD
jet (center), and Higgs jet vs. sgluon jet (right) classifications.

###  3.3 Event Generator Dependence

The classifier introduced in the previous subsection uses not only the
information of hard subjets encoded in $S_{2,\mathrm{trim}}$ but also the soft
activities captured in $S_{2,\mathrm{soft}}$ as well. This leads to concerns
about the accuracy of the models of soft physics. Specifically, the
performance of the classifier could be sensitive to the soft activities in the
jet while the simulated soft activities may be significantly different from
the truth.

In figure 11, we compare the ROC curves of the MLP trained with PY8 and HW7
samples. As these two event generators are based on different modeling of
parton shower and hadronization scheme, the comparison would give us a
reasonable estimate of the systematic uncertainty originating from the
generator choice.

Figure 11:  The ROC curves of the MLP trained on $S_{2,\mathrm{trim}}$ and
$S_{2,\mathrm{soft}}$. The solid and dashed lines correspond to the classifier
trained with PY8 and HW7, respectively. The red and blue lines correspond to
the classifier tested with PY8 and HW7, respectively. The dashed gray lines
represent the ROC curves of the random guess. We show the results of
discriminating Higgs jet vs. QCD jet (left), sgluon jet vs. QCD jet (center),
and Higgs jet vs. sgluon jet (right).

In the left panel of figure 11, we compare the ROC curves of the Higgs jet vs.
QCD jet classification for different generator choices. By doing this
exercise, we estimate a systematic uncertainty in the predictions of the
classifier by comparing ROC(PY8, PY8) and ROC(HW7, HW7) curves, where the
first and second entries in the parenthesis correspond to the generators used
to simulate the training and test samples, respectively. On the other hand,
ROC(HW7, PY8) and ROC(PY8, HW7) show the degradation of the performance of
classifier trained on the “wrong sample” to analyze “real events.”

The performance of the classifier improves as we vary generator combinations
in the following order: ROC(PY8, HW7), ROC(HW7, HW7), ROC(HW7, PY8), and
ROC(PY8, PY8). We find that the classification performance is significantly
better for PY8 test samples than that of HW7 samples. On the other hand, the
classification performance for the same test samples hardly depends on the
classifiers, namely ROC(PY8, HW7) $\sim$ ROC(HW7, HW7) and ROC(HW7, PY8)
$\sim$ ROC(PY8, PY8). For the Higgs jet vs. QCD jet classification, the
classifier mostly concentrates on the core substructures within the jet, and
here both PY8 and HW7 provide similar kinematics and radiation spectra.
Therefore, we do not observe any significant change in the ROC curves by
varying training samples while keeping the test samples the same.

In the middle panel of figure 11, we compare the classifier performance for
the sgluon jet vs. QCD jet classification. It improves in the following order:
ROC(PY8, HW7), ROC(HW7, PY8), ROC(HW7, HW7), and ROC(PY8, PY8). The
classifiers are indeed sensitive to the choice of generators. The network
trained with PY8 (HW7) samples has failed to capture the features of HW7 (PY8)
test samples. The networks have focused on different portions of the
distribution of the fragmentation functions. In the right panel of figure 11,
the ROC curves for the Higgs jet vs. sgluon jet classification show similar
behavior.

##  4 Interpretable Two-level Architecture

A quantitative understanding of a neural network is not straightforward
because the parameters and intermediate outputs of the neural network are less
readable. In this section, we propose an architecture constructed from the
truncated series in eq. (28) and try to explain quantitatively how this
network classifies events. In the case of binary classifications, the
discretized architecture is defined as follows,

$\displaystyle h$ $\displaystyle=$
$\displaystyle\sum_{k}S^{k}_{2,\mathrm{trim}}\,w^{k}_{\mathrm{trim}}+\sum_{k}S^{k}_{2,\mathrm{soft}}\,w^{k}_{\mathrm{soft}},\quad
w_{A}^{k}=\frac{1}{2}\int_{R_{k}}^{R_{k}+\Delta
R_{k}}dR\,w_{A}^{(2)}(R)\;\;(A=\mathrm{trim},\,\mathrm{soft}),$ (45)
$\displaystyle\hat{y}_{1}$ $\displaystyle=$
$\displaystyle\frac{e^{h}}{e^{h}+1},\quad\hat{y}_{2}=1-\hat{y}_{1}$ (46)

where $w^{k}_{A}$ is a trainable weight. We change the activation of the
output layer to a sigmoid activation since the softmax function for binary
classification is essentially a sigmoid with a scale factor on its argument.
The loss function is the categorical cross-entropy as defined in eq. (43).
This setup is effectively a logistic classifier on $S_{2,\mathrm{trim}}^{k}$
and $S_{2,\mathrm{soft}}^{k}$. After the training, the magnitude of
$S_{2,\mathrm{trim}}^{k}w^{k}_{\mathrm{trim}}$ or
$S_{2,\mathrm{soft}}^{k}w^{k}_{\mathrm{soft}}$ is high when the corresponding
$S_{2,\mathrm{trim}}^{k}$ or $S_{2,\mathrm{soft}}^{k}$ is useful for the
classification.

jetradiation module: jet spectrumspectral functions
$S_{2,\mathrm{trim}}$$S_{2,\mathrm{soft}}$kinetic module: weight
generator$p_{T,\mathbf{J}}$$m_{\mathbf{J}}$$w_{\mathrm{soft}}$$w_{\mathrm{trim}}$kinematicsMLPinnerproduct$\displaystyle{\sum}$$\times$$\times$sigmoid$\hat{y}_{1}$
Figure 12:  A schematic diagram of a two-level architecture for binary
classification. An MLP trained on $p_{T,\mathbf{J}}$ and $m_{\mathbf{J}}$
generates weights $w_{\mathrm{trim}}$ and $w_{\mathrm{soft}}$ for analyzing
radiation patterns encoded in $S_{2,\mathrm{trim}}$ and $S_{2,\mathrm{soft}}$
spectra. Double bordered boxes represent trainable modules. The $\hat{y}_{2}$
is given by the normalization, i.e., $\hat{y}_{2}=1-\hat{y}_{1}$.

The logistic classifier does not take into account the $p_{T,\mathbf{J}}$
dependence of $\hat{R}_{b\bar{b}}$; therefore, we introduce a two-level
architecture, which is a variant of the logistic classifier. The weights
$w_{A}^{k}$ are calculated by a kinetic module
$\Phi^{k}_{A}(\vec{x}_{\mathrm{kin}})$ of an MLP trained on
$\vec{x}_{\mathrm{kin}}=(p_{T,\mathbf{J}},m_{\mathbf{J}})$,

$\displaystyle w^{k}_{A}$ $\displaystyle=$
$\displaystyle\Phi^{k}_{A}(\vec{x}_{\mathrm{kin}}).$ (47)

A schematic diagram of this setup is shown in figure 12. The inputs
$\vec{x}_{\mathrm{kin}}$ are standardized before training, and
$S_{2,\mathrm{trim}}^{k}$ and $S_{2,\mathrm{soft}}^{k}$ are divided by their
maximum value of the training sample because standardizing the spectra
reintroduce the zeroth order term of eq. (27). This architecture is similar to
the self-explaining neural network NIPS2018_8003 . The $\Phi^{k}_{A}$ is
modeled with an MLP of two hidden layers with exponential linear unit (ELU)
activations DBLP:journals/corr/ClevertUH15 ,

$\varphi_{\mathrm{ELU}}(x)=\begin{cases}x&x>0,\\\ e^{x}-1&x<0.\end{cases}$
(48)

The nodes of the successive layers are configured as 400 $\mathrm{ELU}$, 200
$\mathrm{ELU}$, 2 $\times$ 20 linear respectively. We do not use
$\mathrm{ReLU}$ in modeling $\Phi^{k}_{A}$ because dead $\mathrm{ReLU}$ nodes
with a zero gradient kill $\vec{x}_{\mathrm{kin}}$ dependency of the weights.
This is known as the dying $\mathrm{ReLU}$ problem. In this case, the
architecture is reduced to the logistic classifier.

The vanishing gradient problem arises when the momentum range of the training
sample is too small, which generates constant weights. The characteristic
scale of $\hat{R}_{b\bar{b}}$ is $[0.625,0.833]$ for the $p_{T,\mathbf{J}}$
range [300, 400] GeV. The variation in $\hat{R}_{b\bar{b}}$ is about 0.2,
which is not significantly large compared with the calorimeter resolution of
0.1. Therefore, we extend the $p_{T,\mathbf{J}}$ range of all the samples to
$[300,\,600]$ GeV. In addition, we avoid vanishing gradient problem by using
He uniform initializer He_2015_ICCV . The weights and biases are initialized
by uniform random numbers in
$[-\sqrt{6/N_{\mathrm{in}}},\sqrt{6/N_{\mathrm{in}}}]$ where $N_{\mathrm{in}}$
is the number of inputs to the layer. The advantage of using He initializer
over Glorot initializer is that it generates random numbers in a wider range
so that the neural network can start up from wider initial weights and
gradient. The weight decay parameter $\lambda$ of the $L_{2}$ weight
regularizer is set to 0.001 so that the weights do not vanish too early.

After the successful training, the performance of the two-level architecture
is close to that of the MLP in section 3. The green dotted lines in figure 10
are the ROC curves of the classifier. The difference is smaller than the
systematic uncertainty shown in figure 11. This makes a good reason to believe
that the weights in eq. (47) capture the essential features of the MLP and CNN
in section 3. The correlation between the output of the two-level architecture
and the CNN model is shown in figure 13. We can see a positive correlation
between them, but the correlation is slightly tilted towards the lower
triangle (upper triangle) for small (large) $\hat{y}_{h}$ values because the
CNN performs better than the two-level architecture.

Figure 13:  Two-dimensional histograms of two $\hat{y}_{h}$’s from the two-
level architecture and the jet image CNN in Higgs jet vs. QCD jet
classification. The left panel is the histogram of Higgs jets, and the right
panel is that of QCD jets.

Figure 14:  The weights $w_{\mathrm{trim}}$ (left) and $w_{\mathrm{soft}}$
(right) of Higgs jet vs. QCD jet classification. We show the weights for
$p_{T,\mathbf{J}}$ values in the range [300, 550] GeV while fixing the
$m_{\mathbf{J}}$ to 115 GeV. The weights at an angular scale smaller than 1.0
are magnified by 10. The value of $w_{\mathrm{soft}}$ in the bin [0,0.1),
i.e., $w_{\mathrm{soft}}^{0}$, is approximately 32 for all the values of
$p_{T,\mathbf{J}}$. In the lower panels, we show the statistical uncertainty
of the weights from the training dataset at $p_{T,\mathbf{J}}=350$ GeV.

In figure 14, we show the weight functions $w_{\mathrm{trim}}$ (left) and
$w_{\mathrm{soft}}$ (right) of Higgs jet vs. QCD jet classifier trained with
MG5+PY8+Delphes samples. Note that the weights are outputs of the neural
network for the given $\vec{x}_{\mathrm{kin}}$ and not the output with the
sample at the indicated $p_{T,\mathbf{J}}$ and $m_{\mathbf{J}}$. The weights
in $R>R_{\mathbf{J}}$ and $w_{\mathrm{soft}}$ in $R<R_{\mathrm{trim}}$ are
large compared to the weights in other angular scales. The
$S_{2,\mathrm{trim}}$ and $S_{2,\mathrm{soft}}$ on these angular scales are
typically smaller than that in other scales. Therefore, their weights become
large to compensate for the energy difference when the corresponding value of
the spectrum is useful for jet classification. The dotted lines in figure 14
denote the values of $\hat{R}_{b\bar{b}}$, defined in eq. (36), for different
values of $p_{T,\mathbf{J}}$.

The $w_{\mathrm{trim}}$ around $\hat{R}_{b\bar{b}}$ is positive because the
correlation at the scale is a characteristic feature of Higgs jet. If a Higgs
boson is decaying to a pair of bottom quarks perpendicular to the boosted
direction in its rest frame, the relative angular separation of the decay
products is $\hat{R}_{b\bar{b}}$ in the lab frame. Due to the phase space of
the decay, most of the events are distributed near $\theta=\pi/2$, where
$\theta$ is Higgs decay angle relative to the boost direction. As
$p_{T,\mathbf{J}}$ increases, $\hat{R}_{b\bar{b}}$ decreases, and the lower
edge of $w_{\mathrm{trim}}>0$ moves toward smaller values. The region with
$w_{\mathrm{trim}}>0$ also shifts towards smaller values of $R$. The weight
$w_{\mathrm{trim}}$ on $R>\hat{R}_{b\bar{b}}$ is positive for capturing a
Higgs jet whose $\theta$ is smaller than $\pi/2$. These events have $p_{T}$
asymmetric subjets, and the cross-correlation terms in
$S_{2,\mathrm{trim}}(R)$ are smaller than that of $p_{T}$ symmetric case.
These correlations are still useful for the classification because $S_{2}$
spectrum of QCD jet reduces much faster than that of Higgs jet. As a result,
the $w_{\mathrm{trim}}$ is an increasing function in this region to compensate
for the $S_{2,\mathrm{trim}}$ reduction.

For $R\gtrsim R_{\mathbf{J}}$, weight $w_{\mathrm{trim}}$ is negative. The
score $\hat{y}_{h}$ decreases whenever there are any energy deposits at
$R\gtrsim R_{\mathbf{J}}$. The crossover point from $w_{\mathrm{trim}}>0$ to
$w_{\mathrm{trim}}<0$ shifts towards smaller values of $R$ with an increase of
$p_{T,\mathbf{J}}$ because the Higgs decay products become more asymmetric
with respect to the boost direction. In such a case, one of the subjets tends
to be soft so that the two-point correlations are included in
$S_{2,\mathrm{soft}}$, instead of $S_{2,\mathrm{trim}}$. These contributions
to $S_{2,\mathrm{soft}}$ do not affect $w_{\mathrm{soft}}$, because
$S_{2,\mathrm{soft}}$ spectra of QCD jets are overwhelming at large $R$.

The $S_{2,\mathrm{soft}}$ on $R>R_{\mathrm{trim}}$ always reduces
$\hat{y}_{h}$, and there is no prominent structure around
$\hat{R}_{b\bar{b}}$. Moreover, $|w_{\mathrm{soft}}|$ decreases as
$p_{T,\mathbf{J}}$ increases. The reduction of $w_{\mathrm{soft}}$ compensates
the increase of $S_{2,\mathrm{soft}}$, and the prediction is more or less
$p_{T,\mathbf{J}}$ independent. On $R>R_{\mathbf{J}}$, $|w_{\mathrm{soft}}|$
increases with $R$ because activity in this region is a sign of QCD jet even
though corresponding $S_{2,\mathrm{soft}}$ decreases due to suppressed large
angle radiations.

The $w_{\mathrm{soft}}$ on $R\lesssim R_{\mathrm{trim}}$ is positive and
$w_{\mathrm{soft}}$ has a break at $R\sim R_{\mathrm{trim}}$. Correlations
between the constituents in a soft subjet contributes to the
$S_{2,\mathrm{soft}}$ on $R<R_{\mathrm{trim}}$, i.e.,
$S_{2,\mathrm{soft}}(0;R_{\mathrm{trim}})\propto f^{2}_{\mathrm{trim}}$. Let
us assume $\mathbf{J}_{a}$ is a single soft subjet, then
$S_{2,\mathrm{soft}}(0;R_{\mathrm{trim}})\sim
p_{T,\mathbf{J}_{a}}^{2}\sim(p_{T,\mathbf{J}}-p_{T,\mathbf{J},\mathrm{trim}})^{2}$.
If there are multiple soft subjets, then
$S_{2,\mathrm{soft}}(0;R_{\mathrm{trim}})\sim\sum_{a}p_{T,\mathbf{J}_{a}}^{2}<(\sum_{a}p_{T,\mathbf{J}_{a}})^{2}\sim(p_{T,\mathbf{J}}-p_{T,\mathbf{J},\mathrm{trim}})^{2}$.
This triangular inequality suggests that the magnitude of
$S_{2,\mathrm{soft}}(0;R_{\mathrm{trim}})$ is small for a jet with a given
$p_{T,\mathbf{J}}-p_{T,\mathbf{J},\mathrm{trim}}$ when there are multiple soft
jets. The positive $w_{2}$ on $R<R_{\mathrm{trim}}$ means that Higgs jet has
less soft subjets than QCD jet. The $S_{2,\mathrm{soft}}^{0}$ consists of the
autocorrelation of soft subjets, which has different energy scaling behavior
compared with the other $S_{2,\mathrm{soft}}^{k}$. The $S_{2,\mathrm{soft}}$
on $R<R_{\mathrm{trim}}$ is $S_{2,22}\sim\mathcal{O}[f_{\mathrm{trim}}^{2}]$
in eq. (24). On the other hand, $S_{2,\mathrm{soft}}^{k}$ $(k\geq 1)$ is
dominated by $S_{2,12}$. The $S_{2,12}$ on $R<R_{\mathrm{trim}}$ does not
contribute to $S_{2,\mathrm{soft}}^{0}$ because it vanishes. Therefore, we may
rewrite $h$ as follows,

$\displaystyle h=\int
dR\,S_{2,\mathrm{trim}}(R)w_{\mathrm{trim}}(R)+\int_{0}^{R_{\mathrm{trim}}}dR\,S_{2,\mathrm{soft}}(R)w_{\mathrm{soft}}^{\prime}(R)+\int
dR\,S_{2,\mathrm{soft}}(R)w_{\mathrm{soft}}^{\prime\prime}(R)$ (49)

where $w_{\mathrm{soft}}^{\prime}$ and $w_{\mathrm{soft}}^{\prime\prime}$ are
continuous functions with $w_{\mathrm{soft}}^{\prime}(R_{\mathrm{trim}})=0$
and
$w_{\mathrm{soft}}(R)=w_{\mathrm{soft}}^{\prime}(R)+w_{\mathrm{soft}}^{\prime\prime}(R)$.
The second term is essentially the same as $\int
dR\,S_{2,22}(R)w_{2,22}^{(2)}(R)$ in eq. (23), and the last term is $\int
dR\,S_{2,12}(R)w_{2,12}^{(2)}(R)+\int dR\,S_{2,21}(R)w_{2,21}^{(2)}(R)$.

The sudden changes of $w_{\mathrm{trim}}$ and $w_{\mathrm{soft}}$ at $R\simeq
1.8$ are due to the statistical fluctuations of the training sample. The
$S_{2,\mathrm{trim}}$ and $S_{2,\mathrm{soft}}$ may have a non-zero value at
large R if the jet has multiple large angle radiations with the opposite
direction from the jet axis; however, the probability of such a radiation
pattern is small. As a result, the number of events used for training the
weights at large $R$ is not sufficient. The large weights also do not
contribute much to the classification (see figure 16).

To estimate this statistical uncertainty, we use both training and test
datasets. The merged dataset is divided into ten subsets and we train a
network for each of them. This will decrease the number of events in each
subset by a factor 5. We estimate the uncertainty of the fit by calculating
the mean and variance of the $w_{\mathrm{trim}}^{k}$ and
$w_{\mathrm{soft}}^{k}$ from the ten subsets. The green band in the bottom
panels of figure 14 represents the estimated uncertainty. Note that
$w_{\mathrm{trim}}^{k}$ and $w_{\mathrm{soft}}^{k}$ are not sensitive to the
network initialization because the two-level network is effectively a logistic
regression for a fixed $p_{T,\mathbf{J}}$ and $m_{\mathbf{J}}$ and the loss
function is a convex function of them.

Figure 15:  The weights $w_{\mathrm{trim}}$ (left) and $w_{\mathrm{soft}}$
(right) at $p_{T,\mathbf{J}}=350$ GeV and $m_{\mathbf{J}}=115$ GeV for the
sgluon jet vs. QCD jet (top) and the Higgs jet vs. sgluon jet (bottom)
classification. We show the weights with statistical uncertainty from the
training dataset.

In the top panels of figure 15, we show the weights $w_{\mathrm{trim}}$ and
$w_{\mathrm{soft}}$ for the sgluon jet vs. QCD jet classification with
$p_{T,\mathbf{J}}=350$ GeV and $m_{\mathbf{J}}=115$ GeV. The
$w_{\mathrm{trim}}$ distribution is similar to that of the Higgs jet vs. QCD
jet classification. However, the $|w_{\mathrm{soft}}|$ is much smaller. This
comes from the fact that $S_{2,\mathrm{soft}}$ of sgluon jet is similar to
that of QCD jet and it is less important in the classification. No peak of
$S_{2,\mathrm{soft}}$ around $R\lesssim R_{\mathrm{trim}}$ also indicates that
the soft substructures of sgluon jet are as radiative as QCD jet.
Additionally, there is no color coherence restriction of soft radiations for
the sgluon jet. This leads to small $|w_{\mathrm{soft}}|$ for
$R>R_{\mathbf{J}}$. In the bottom panels of figure 15, we show the weights for
the Higgs jet vs. sgluon jet classification. The peak of $w_{\mathrm{trim}}$
around $R=\hat{R}_{b\bar{b}}$ is small as the hard substructures of Higgs jet
and sgluon jet are (almost) the same. However, a sgluon is more radiative than
a Higgs boson, and $w_{\mathrm{soft}}$ is negative in the entire region of
$R<1.5$.

As described above, weights $w_{\mathrm{trim}}$ and $w_{\mathrm{soft}}$ may
take large values, but it does not necessarily mean that the corresponding
$S_{2,\mathrm{trim}}$ and $S_{2,\mathrm{soft}}$ contribute dominantly in the
jet classification. The energy scaling factors on the $S_{2,\mathrm{trim}}$
($S_{2,\mathrm{soft}}$) and its weight $w_{\mathrm{trim}}$
($w_{\mathrm{soft}}$) cancel out in the quantity of our interest
$h=\sum_{k}(S_{2,\mathrm{trim}}^{k}w_{\mathrm{trim}}^{k}+S_{2,\mathrm{soft}}^{k}w_{\mathrm{soft}}^{k})$.
For example, $\mathcal{O}[1]$ terms in $S_{2,\mathrm{trim}}$ and
$\mathcal{O}[f_{\mathrm{trim}}]$ terms in $S_{2,\mathrm{soft}}$ contribute
equally to the classifier if $w_{\mathrm{soft}}$ is around
$f_{\mathrm{trim}}\,w_{\mathrm{trim}}$. In the left panel of figure 16, we
draw the mean values $\langle
S_{2,\mathrm{trim}}^{k}w^{k}_{\mathrm{trim}}\rangle$ and $\langle
S_{2,\mathrm{soft}}^{k}w^{k}_{\mathrm{soft}}\rangle$ of Higgs jet vs. QCD jet
classification, which are more directly related to the jet classification. The
solid and dashed red lines correspond to the distributions of the Higgs jet,
while the solid and dashed blue lines are for the QCD jet. The regions where
Higgs jet and QCD jet distributions differ significantly are important for the
network predictions. We find $\langle
S_{2,\mathrm{trim}}^{k}w^{k}_{\mathrm{trim}}\rangle$ around
$R\sim\hat{R}_{b\bar{b}}$ and $\langle
S_{2,\mathrm{soft}}^{k}w^{k}_{\mathrm{soft}}\rangle$ in the region $R<1.2$
mostly contribute to the jet classification.

Figure 16:  The distribution of the mean value $\langle
S_{2,\mathrm{trim}}(R)\,w_{\mathrm{trim}}(R)\rangle$ (solid) and $\langle
S_{2,\mathrm{soft}}(R)\,w_{\mathrm{soft}}(R)\rangle$ (dashed). We show the
$\langle S_{2,A}(R)\,w_{A}(R)\rangle$ for Higgs jet vs. QCD jet (top), sgluon
jet vs. QCD jet (center), and Higgs jet vs. sgluon jet (bottom)
classifications. In the right figure, we additionally demand that
$\hat{y}_{h}$ of the Higgs jet and $\hat{y}_{\mathrm{QCD}}$ of the QCD jet are
larger than their 95th percentile respectively. We show their statistical
uncertainty from the training samples as colored bands.

The average distribution may not illustrate all the features of the classifier
performance. The energy deposits in each bin fluctuate, and the bins with hits
higher than the average value contribute more to the network decisions. For
example, soft emissions outside the angle between the two hardest subjets are
rare in the Higgs jet. Once there is large angle radiation outside the cone of
hard subjets, the network is likely to identify the jet as a QCD jet. In the
right panel of figure 16, we plot the $\langle
S_{2,\mathrm{trim}}^{k}w_{\mathrm{trim}}^{k}\rangle$ and $\langle
S_{2,\mathrm{soft}}^{k}w_{\mathrm{soft}}^{k}\rangle$ distributions of the
Higgs jet (QCD jet) with $\hat{y}_{h}$ ($\hat{y}_{\mathrm{QCD}}$) higher than
the 95th percentile. The distributions indicate that the selected Higgs jets
are mostly classified because of the $S_{2,\mathrm{trim}}$ excess at
$\hat{R}\sim R_{b\bar{b}}$, while the QCD jets are classified using
$S_{2,\mathrm{soft}}$ excess above $R>0.2$.

Figure 17:  Top figures are the weights $w_{\mathrm{trim}}$ (left) and
$w_{\mathrm{soft}}$ (right) at $p_{T,\mathbf{J}}=350$ GeV and
$m_{\mathbf{J}}=115$ GeV for classifying Higgs jet of PY8 and HW7 events.
Bottom figures are $\langle S_{2,A}(R)\,w_{A}(R)\rangle$. In the right bottom
figure, we additionally demand that $\hat{y}_{1}$ of the PY8 generated jets
and $\hat{y}_{2}$ of the HW7 generated jets are larger than their 95th
percentile respectively. We show their statistical uncertainty from the
training samples as colored bands.

Figure 18:  Top figures are the weights $w_{\mathrm{trim}}$ (left) and
$w_{\mathrm{soft}}$ (right) at $p_{T,\mathbf{J}}=350$ GeV and
$m_{\mathbf{J}}=115$ GeV for classifying sgluon jet of PY8 and HW7 events.
Bottom figures are $\langle S_{2,A}(R)\,w_{A}(R)\rangle$. In the right bottom
figure, we additionally demand that $\hat{y}_{1}$ of the PY8 generated jets
and $\hat{y}_{2}$ of the HW7 generated jets are larger than their 95th
percentile respectively. We show their statistical uncertainty from the
training samples as colored bands.

Figure 19:  Top figures are the weights $w_{\mathrm{trim}}$ (left) and
$w_{\mathrm{soft}}$ (right) at $p_{T,\mathbf{J}}=350$ GeV and
$m_{\mathbf{J}}=115$ GeV for classifying QCD jet of PY8 and HW7 events. Bottom
figures are $\langle S_{2,A}(R)\,w_{A}(R)\rangle$. In the right bottom figure,
we additionally demand that $\hat{y}_{1}$ of the PY8 generated jets and
$\hat{y}_{2}$ of the HW7 generated jets are larger than their 95th percentile
respectively. We show their statistical uncertainty from the training samples
as colored bands.

We now use the two-level architecture to compare PY8 and HW7. As we have
already shown in section 3, the performance of the classifier depends on the
event generators significantly. We show the weights of the classifiers for
$p_{T,\mathbf{J}}=350$ GeV and $m_{\mathbf{J}}=115$ GeV trained with Higgs
jets in figure 17, sgluon jets in figure 18, and QCD jets in figure 19. For
each plot, the signals are PY8 events, and the backgrounds are HW7 events. The
$w_{\mathrm{trim}}$ in $R\lesssim R_{\mathbf{J}}$ is close to zero everywhere,
representing $S_{2,\mathrm{trim}}$ spectra of PY8 and HW7 events are similar.
It is not surprising because both PY8 and HW7 events from identical hard
partons and these partons create the trimmed subjets inside the jet. On the
other hand, the correlation involving constituents of the soft activities,
$S_{2,\mathrm{soft}}$, is manifestly different and so $w_{\mathrm{soft}}$ is
nonzero. For Higgs jet and sgluon jet, the $w_{\mathrm{soft}}$ distribution of
PY8 events is significantly large (and positive) for $R\sim R_{\mathrm{trim}}$
and it decreases as $R$ increases. The weight $w_{\mathrm{soft}}$ is negative
for $R>R_{\mathbf{J}}$, which means that HW7 events have more soft activity in
the region $R\gg\hat{R}_{b\bar{b}}$.

For the case of QCD jet, the distribution of $w_{\mathrm{soft}}$ is always
positive and flat for $R<R_{\mathbf{J}}$ and negative for $R>1.5$. It would be
interesting to evaluate the weights for the classifiers trained with the
experimental data and compare with the simulated results to tune the
parameters of the event generators further.

##  5 Summary and Outlook

The classification of jets with deep learning has gained significant attention
in recent times. Majority of these analyses take advantage of the significant
development in computing power. These deep learning architectures utilize the
complete event information in terms of low-level observables. These deep
learning based strategies can be compared with the previous approaches for
tagging jets, for example, mass drop tagger, $n$-subjettiness, energy
correlation function where each of them has solid physics motivation.

In this paper, we introduce neural networks trained on “jet spectrum”
$S_{2}(R)$, which are essentially two-point correlation functions of jet
constituents. We also introduce $S_{2,\mathrm{trim}}(R)$, which is $S_{2}(R)$
calculated from the trimmed jet, to encode hard substructure of the jet. The
difference $S_{2,\mathrm{soft}}(R)=S_{2}(R)-S_{2,\mathrm{trim}}(R)$ encodes
the remaining correlations with soft radiations and is less affected by the
correlations among the hard constituents. Our neural networks are trained on
$S_{2,\mathrm{trim}}$ and $S_{2,\mathrm{soft}}$ integrated over certain bins.
If the $S_{2,\mathrm{trim}}$ and $S_{2,\mathrm{soft}}$ spectra are multiplied
by smooth functions and integrated over $R$, it forms an IRC safe
C-correlator. This feature assures that the classifiers trained on binned
$S_{2,\mathrm{trim}}$ and $S_{2,\mathrm{soft}}$ are approximately IRC safe.

The performance of MLP trained on $S_{2,\mathrm{trim}}$ and
$S_{2,\mathrm{soft}}$ is compared to that of CNN trained on jet images. The
CNN shows better performance than the MLP, but the difference is small. The
key reason is efficient preprocessing of parton shower effects with less
number of free parameters. Parton shower is the multiple splittings of the
partons where each splitting is parametrized by the angular scale and momentum
fraction of the partons. The binned $S_{2}(R)$ spectra collect the information
of the parton splitting successfully. The spectra provide comparable jet
classification performance with a fewer number of inputs. Furthermore, the MLP
is computationally economical than the CNN because the MLP has smaller
complexity than the CNN and takes only $\mathcal{O}(40)$ inputs.

The $S_{2,\mathrm{trim}}(R)$ and $S_{2,\mathrm{soft}}(R)$ spectra can be
obtained from a functional Taylor series of an arbitrary classifier in energy
flows. The spectra are basis vectors of the second order term in the
expansion. In this context, the MLP trained on $S_{2,\mathrm{trim}}(R)$ and
$S_{2,\mathrm{soft}}(R)$ can be considered as a sum of $2n$-linear
$C$-correlators which can be reduced to products of the bilinear
$C$-correlators in $S_{2,\mathrm{trim}}$ and $S_{2,\mathrm{soft}}$. The (mild)
difference in the performance of the CNN and MLP comes from the remaining
irreducible $n$-linear $C$-correlators.

The terms linear in $S_{2,\mathrm{trim}}(R)$ and $S_{2,\mathrm{soft}}(R)$
provide an opportunity to visualize and interpret the network predictions;
therefore, we study a novel two-level architecture that involves an
interpretable layer of a single node in the form of a $C$-correlator. The
output is the sum of the product of the trained weights ($w_{\mathrm{trim}}$
and $w_{\mathrm{soft}}$) and jet spectra ($S_{2,\mathrm{trim}}$ and
$S_{2,\mathrm{soft}}$). The absolute values of the weights signify the impact
of the corresponding $S_{2,\mathrm{trim}}$ and $S_{2,\mathrm{soft}}$ bin
values on the jet classification. In the context of classification between
Higgs jet and QCD jet, the distribution of $w_{\mathrm{trim}}$ shows that
$S_{2,\mathrm{trim}}$ spectrum around
$\hat{R}_{b\bar{b}}=2m_{h}/p_{T,\mathbf{J}}$ increases the output, and the
classifier regards the jet as a Higgs jet. We have also shown that the
dependence of $w_{\mathrm{trim}}$ on jet $p_{T}$ can be qualitatively
understood (at the parton level) from the decay of a boosted Higgs boson. In
short, the network is using $S_{2,\mathrm{trim}}$ inputs to obtain the core
substructure information inside the jet.

The soft activity is also useful for Higgs jet vs. QCD jet classification. The
probability for assigning the given jet as a QCD jet increases with increase
in $S_{2,\mathrm{soft}}$ on $R>R_{\mathrm{trim}}$. To study the impact of soft
physics in jet classification, we also introduce sgluon, a hypothetical color
octet scalar, and compare the classifier performance among Higgs jet, sgluon
jet, and QCD jet. The network predictions for sgluon jet vs. QCD jet
classification are primarily determined by the core substructure information
as expected. However, the network uses the difference in the
$S_{2,\mathrm{soft}}$ spectra arising from the different color structure of
the decaying particle for Higgs jet vs. sgluon jet classification.

The non-trivial role of soft radiations in the predictions of the classifiers
implies the results are highly sensitive to the choice of event generators.
The weights associated with the $S_{2,\mathrm{trim}}$ are almost insensitive
to the choice; however, the weights of the $S_{2,\mathrm{soft}}$ are strongly
affected. This behavior is expected as modeling of soft physics is quite
different in Pythia 8 and Herwig 7.

The two-point correlation spectra and the architectures introduced in this
paper can be applied for solving other interesting problems, thanks to
flexibility on designing neural network. For jets with more complex
substructures, e.g., top jet, the higher order terms in the energy flow series
expansion may be included. It would be worthwhile to study the classifier
performances when the network is trained with the experimental data and
compare with the predictions of event generators to tune their parameters to
reduce the uncertainty in modeling the soft physics. It is also interesting to
use this interpretable architecture as a model-agnostic interpreter for black
box architectures 2016arXiv160605386T . We leave these possibilities for
future works.

## Acknowledgements

We thank the organizers of _Beyond the BSM_ and _Machine Learning for Jet
Physics 2018_ workshops. This work is supported by the Grant-in-Aid for
Scientific Research on Scientific Research B (No.16H03991, 17H02878) and
Innovative Areas (16H06492), and by World Premier International Research
Center Initiative (WPI Initiative), MEXT, Japan.

##  Appendix A Event generation and reconstruction

The parton level event samples, namely $pp\rightarrow Zj$, $pp\rightarrow Zh$
and $pp\rightarrow Z\sigma$ events, are generated at the leading order in QCD
using MadGraph5_aMC@NLO 2.6.1 Alwall:2014hca . We force the Higgs boson ($h$)
and the sgluon ($\sigma$) to decay to a pair of bottom quarks, while $Z$ boson
to decay invisibly. For sgluon, we use a UFO model in Degrande:2014sta ;
NLOModels  with the following interaction term for the decay,

$\displaystyle\mathcal{L_{\mathrm{sgluon}}}$ $\displaystyle\ni$ $\displaystyle
y_{\sigma b\bar{b}}\,\sigma^{a}\,\bar{b}T^{a}b+\mathrm{h.c.}$ (50)

The parton distribution function (PDF) set NNPDF 2.3 LO at
$\alpha_{S}(m_{Z})=0.130$ Ball:2012cx  is used. To generate Higgs jets and
Sgluon jets, we impose a parton level selection criterion on the $Z$ boson
transverse momentum, $p_{T,Z}>250$ GeV. We simulate approximately 3 million
events of $pp\rightarrow Zh$ and $pp\rightarrow Z\sigma$ processes and 18
million events of $pp\rightarrow Zj$.

We use two parton shower and hadronization simulators to compare the results.
Namely, we use Pythia 8.226 Sjostrand:2014zea  with Monash tune Skands:2014pea
and Herwig 7.1.3 Bellm:2015jjp ; Bahr:2008pv  with default tune herwigtune ;
Gieseke:2012ft . The shower starting scale is $H_{T}/2$ for $pp\rightarrow Zh$
and $pp\rightarrow Z\sigma$ processes and $p_{T,j}$ for $pp\rightarrow Zj$
process, where $H_{T}$ is the transverse energy sum of the produced partons.
The effects of underlying events and multi-parton interactions are taken into
account, but we neglect the contaminations coming from the pile-ups. The PDF
set for simulating all these effects are the same as that in the parton level
simulation.

We use Delphes 3.4.1 deFavereau:2013fsa  with its default ATLAS configuration
for fast detector simulations. Jets are reconstructed from the calorimeter
towers using FastJet 3.3.0 Cacciari:2011ma ; Cacciari:2005hq  with
anti-$k_{T}$ algorithm Cacciari:2008gp  and jet radius parameter
$R_{\mathbf{J}}=1$. The leading jet of each event with
$p_{T,\mathbf{J}}\in[300,400]\,\mathrm{GeV}$ and
$m_{\mathbf{J}}\in[100,150]\,\mathrm{GeV}$ is selected. Since the scale
$H_{T}/2$ for $pp\rightarrow Zh$ and $pp\rightarrow Z\sigma$ is higher than
$p_{T,h}$ and $p_{T,\sigma}$, respectively, there is a chance that the leading
jet is from the initial state radiation rather than from the decay of Higgs
boson or sgluon. To filter out such jets from the Higgs jet and sgluon jet
samples, we require that $b$-partons produced from the decay are within
$R_{\mathbf{J}}$ from the leading jet axis.

##  Appendix B Oversampling and $p_{T,\mathbf{J}}$-bias removal

Figure 20: Validation of the $p_{T}$-bias removal on the training samples.
Each row is a conditional probability density on a given $p_{T,\mathbf{J}}$
range, i.e., the sum of each row is 1. Upper (lower) panel displays the
conditional probability density histograms of the predicted class vector
$\hat{y}_{h}$ ($\hat{y}_{\sigma}$) for the Higgs jet (left), sgluon jet
(center), and QCD jet (right).

The neural networks may learn the inherent $p_{T,\mathbf{J}}$ difference among
Higgs jet, sgluon jet, and QCD jet in figure 4 to classify them instead of
learning the difference in their substructures. To penalize the learning from
the $p_{T,\mathbf{J}}$ distribution, we augment the training and validation
samples by oversampling as follows,

  1. 1.

The samples (of each class) are binned in $p_{T,\mathbf{J}}$ with bin-width 1
GeV with $b_{i}$ entries in the $i$-th bin and $b_{\max}=\max\\{b_{i}\\}$.

  2. 2.

For each bin, oversample the events so that the number of the bin contents
becomes a certain value $n_{\max}=c_{d}\,b_{\max}$. This oversampling is
identical to repeating events in bin $i$ sequentially for

$M_{i}=\mathrm{ceil}\left(\frac{c_{d}\cdot b_{\max}}{b_{i}}\right)$ (51)

times and stop the oversampling when the number of events in the bin reaches
$n_{\max}$. It may introduce a small bias due to unequal oversampling among
different bins. We choose $c_{d}=2$ and ignore the small residual bias.

In the upper (lower) panel of figure 20, we show the conditional probability
density of the predicted class vector $\hat{y}_{h}$ ($\hat{y}_{\sigma}$) for a
given $p_{T}$. The probability density has a mild dependence on
$p_{T,\mathbf{J}}$ which is originating from the interplay of the phase-space
selection and the jet radius parameter.

##  Appendix C Jet Image and Convolution Neural Network

jetjet imaging
CNN$p_{T,\mathbf{J}}$$m_{\mathbf{J}}$$p_{T,\mathbf{J},\mathrm{trim}}$$m_{\mathbf{J},\mathrm{trim}}$kinematicsMLPsoftmaxclassifier$\hat{y}_{1}$$\hat{y}_{2}$
Figure 21:  A schematic diagram of a convolutional neural network trained on
jet image. The double bordered boxes represent trainable modules.

We obtain and pre-process the jet image as follows,444This setup is similar to
Kasieczka:2017nvn .

  1. 1.

Recluster the jet constituents by $k_{T}$ algorithm with a jet radius
parameter $R_{\mathbf{J}}=0.2$.

  2. 2.

Set the center of $(\eta,\phi)$ coordinate to the leading (in $p_{T}$) subjet.

  3. 3.

If a second leading subjet is found, rotate the jet constituents on
$(\eta,\phi)$ plane about the jet center so that the sub-leading jet is on the
positive $y$-axis.

  4. 4.

If a third leading subjet is found, flip the image about $y$-axis when $x$
coordinate of the subjet is negative.

  5. 5.

Select jet constituents within $[-1.5,1.5]\otimes[-1.5,1.5]$.

  6. 6.

Finally, pixelate the jet constituents with pixel size $0.1\times 0.1$. The
$(k,l)$-th pixel intensity $P_{T}^{k,l}$ is determined by the total transverse
energy of the jet constituents present in a given pixel, i.e.,

$P_{T}^{k,l}=\sum_{i\in\mathbf{J}}p_{T,i}I_{\mathrm{bin}_{k,l}}(\vec{R}_{i})=\int_{\mathrm{bin}_{k,l}}d\vec{R}\,P_{T}(\vec{R})$
(52)

where $\mathrm{bin}_{k,l}$ is the region of $(k,l)$-th bin.

  7. 7.

Standardize all the $P_{T}^{k,l}$.

This jet image is analyzed by a CNN which consists of two-dimensional
convolutional layers (CONV) and max-pooling layers (figure 21). In particular,
we use the following CNN setup,

  * •

Layer 1: Convolutional layer with 64 filters with kernel size $3\times 3$, and
$\mathrm{ReLU}$ activation,

  * •

Layer 2: Max-pooling layer with pool size $2\times 2$,

  * •

Layer 3: Convolutional layer with 32 filters with kernel size $4\times 4$, and
$\mathrm{ReLU}$ activation,

  * •

Layer 4: Max-pooling with pool size $2\times 2$.

The first convolutional layer deals with angular scale up to 0.3 to treat
collinear radiations while the second convolutional layer operates up to 0.8.
The outputs of Layer 4 are flattened into a one-dimensional array and
concatenated with a set of kinematic inputs,
$\\{p_{T,\mathbf{J}},m_{\mathbf{J}},p_{T,\mathbf{J},\mathrm{trim}},m_{\mathbf{J},\mathrm{trim}}\\}$.
The flattened output array is fed into an MLP with two hidden layers with 300,
100 filters respectively, and $\mathrm{ReLU}$ activation function. The outputs
of the MLP are fed into a softmax layer to make a prediction. The training
setup is the same as in section 3.

## References

  * (1) J. M. Butterworth, A. R. Davison, M. Rubin and G. P. Salam, _Jet substructure as a new Higgs search channel at the LHC_ , _Phys. Rev. Lett._ 100 (2008) 242001, [0802.2470].
  * (2) J. Thaler and L.-T. Wang, _Strategies to Identify Boosted Tops_ , _JHEP_ 07 (2008) 092, [0806.0023].
  * (3) D. E. Kaplan, K. Rehermann, M. D. Schwartz and B. Tweedie, _Top Tagging: A Method for Identifying Boosted Hadronically Decaying Top Quarks_ , _Phys. Rev. Lett._ 101 (2008) 142001, [0806.0848].
  * (4) T. Plehn, G. P. Salam and M. Spannowsky, _Fat Jets for a Light Higgs_ , _Phys. Rev. Lett._ 104 (2010) 111801, [0910.5472].
  * (5) T. Plehn, M. Spannowsky, M. Takeuchi and D. Zerwas, _Stop Reconstruction with Tagged Tops_ , _JHEP_ 10 (2010) 078, [1006.2833].
  * (6) D. E. Soper and M. Spannowsky, _Finding physics signals with shower deconstruction_ , _Phys. Rev._ D84 (2011) 074002, [1102.3480].
  * (7) D. E. Soper and M. Spannowsky, _Finding top quarks with shower deconstruction_ , _Phys. Rev._ D87 (2013) 054012, [1211.3140].
  * (8) M. Dasgupta, A. Fregoso, S. Marzani and G. P. Salam, _Towards an understanding of jet substructure_ , _JHEP_ 09 (2013) 029, [1307.0007].
  * (9) D. E. Soper and M. Spannowsky, _Finding physics signals with event deconstruction_ , _Phys. Rev._ D89 (2014) 094005, [1402.1189].
  * (10) A. J. Larkoski, S. Marzani, G. Soyez and J. Thaler, _Soft Drop_ , _JHEP_ 05 (2014) 146, [1402.2657].
  * (11) J. Gallicchio and M. D. Schwartz, _Seeing in Color: Jet Superstructure_ , _Phys. Rev. Lett._ 105 (2010) 022001, [1001.5027].
  * (12) J. Thaler and K. Van Tilburg, _Identifying Boosted Objects with N-subjettiness_ , _JHEP_ 03 (2011) 015, [1011.2268].
  * (13) J. Gallicchio and M. D. Schwartz, _Quark and Gluon Tagging at the LHC_ , _Phys. Rev. Lett._ 107 (2011) 172001, [1106.3076].
  * (14) Y.-T. Chien, _Telescoping jets: Probing hadronic event structure with multiple R ’s_ , _Phys. Rev._ D90 (2014) 054008, [1304.5240].
  * (15) A. J. Larkoski, G. P. Salam and J. Thaler, _Energy Correlation Functions for Jet Substructure_ , _JHEP_ 06 (2013) 108, [1305.0007].
  * (16) A. J. Larkoski, I. Moult and D. Neill, _Power Counting to Better Jet Observables_ , _JHEP_ 12 (2014) 009, [1409.6298].
  * (17) I. Moult, L. Necib and J. Thaler, _New Angles on Energy Correlation Functions_ , _JHEP_ 12 (2016) 153, [1609.07483].
  * (18) L. G. Almeida, M. Backović, M. Cliche, S. J. Lee and M. Perelstein, _Playing Tag with ANN: Boosted Top Identification with Pattern Recognition_ , _JHEP_ 07 (2015) 086, [1501.05968].
  * (19) L. de Oliveira, M. Kagan, L. Mackey, B. Nachman and A. Schwartzman, _Jet-images — deep learning edition_ , _JHEP_ 07 (2016) 069, [1511.05190].
  * (20) P. T. Komiske, E. M. Metodiev and M. D. Schwartz, _Deep learning in color: towards automated quark/gluon jet discrimination_ , _JHEP_ 01 (2017) 110, [1612.01551].
  * (21) G. Kasieczka, T. Plehn, M. Russell and T. Schell, _Deep-learning Top Taggers or The End of QCD?_ , _JHEP_ 05 (2017) 006, [1701.08784].
  * (22) G. Louppe, K. Cho, C. Becot and K. Cranmer, _QCD-Aware Recursive Neural Networks for Jet Physics_ , 1702.00748.
  * (23) P. T. Komiske, E. M. Metodiev, B. Nachman and M. D. Schwartz, _Pileup Mitigation with Machine Learning (PUMML)_ , _JHEP_ 12 (2017) 051, [1707.08600].
  * (24) A. Butter, G. Kasieczka, T. Plehn and M. Russell, _Deep-learned Top Tagging with a Lorentz Layer_ , _SciPost Phys._ 5 (2018) 028, [1707.08966].
  * (25) T. Cheng, _Recursive Neural Networks in Quark/Gluon Tagging_ , _Comput. Softw. Big Sci._ 2 (2018) 3, [1711.02633].
  * (26) A. Andreassen, I. Feige, C. Frye and M. D. Schwartz, _JUNIPR: a Framework for Unsupervised Machine Learning in Particle Physics_ , _Eur. Phys. J._ C79 (2019) 102, [1804.09720].
  * (27) S. Choi, S. J. Lee and M. Perelstein, _Infrared Safety of a Neural-Net Top Tagging Algorithm_ , 1806.01263.
  * (28) S. H. Lim and M. M. Nojiri, _Spectral Analysis of Jet Substructure with Neural Networks: Boosted Higgs Case_ , _JHEP_ 10 (2018) 181, [1807.03312].
  * (29) F. A. Dreyer, G. P. Salam and G. Soyez, _The Lund Jet Plane_ , _JHEP_ 12 (2018) 064, [1807.04758].
  * (30) J. Lin, M. Freytsis, I. Moult and B. Nachman, _Boosting $H\to b\bar{b}$ with Machine Learning_, _JHEP_ 10 (2018) 101, [1807.10768].
  * (31) P. T. Komiske, E. M. Metodiev and J. Thaler, _Energy Flow Networks: Deep Sets for Particle Jets_ , _JHEP_ 01 (2019) 121, [1810.05165].
  * (32) J. Arjona Martínez, O. Cerri, M. Pierini, M. Spiropulu and J.-R. Vlimant, _Pileup mitigation at the Large Hadron Collider with Graph Neural Networks_ , 1810.07988.
  * (33) G. Kasieczka, N. Kiefer, T. Plehn and J. M. Thompson, _Quark-Gluon Tagging: Machine Learning meets Reality_ , 1812.09223.
  * (34) H. Qu and L. Gouskos, _ParticleNet: Jet Tagging via Particle Clouds_ , 1902.08570.
  * (35) A. J. Larkoski, I. Moult and B. Nachman, _Jet Substructure at the Large Hadron Collider: A Review of Recent Advances in Theory and Machine Learning_ , 1709.04464.
  * (36) L. Asquith et al., _Jet Substructure at the Large Hadron Collider : Experimental Review_ , 1803.06991.
  * (37) L. M. Dery, B. Nachman, F. Rubbo and A. Schwartzman, _Weakly Supervised Classification in High Energy Physics_ , _JHEP_ 05 (2017) 145, [1702.00414].
  * (38) T. Cohen, M. Freytsis and B. Ostdiek, _(Machine) Learning to Do More with Less_ , _JHEP_ 02 (2018) 034, [1706.09451].
  * (39) E. M. Metodiev, B. Nachman and J. Thaler, _Classification without labels: Learning from mixed samples in high energy physics_ , _JHEP_ 10 (2017) 174, [1708.02949].
  * (40) G. Louppe, M. Kagan and K. Cranmer, _Learning to Pivot with Adversarial Networks_ , 1611.01046.
  * (41) C. Shimmin, P. Sadowski, P. Baldi, E. Weik, D. Whiteson, E. Goul et al., _Decorrelated Jet Substructure Tagging using Adversarial Neural Networks_ , _Phys. Rev._ D96 (2017) 074034, [1703.03507].
  * (42) A. Chakraborty, A. M. Iyer and T. S. Roy, _A Framework for Finding Anomalous Objects at the LHC_ , _Nucl. Phys._ B932 (2018) 439–470, [1707.07084].
  * (43) J. Hajer, Y.-Y. Li, T. Liu and H. Wang, _Novelty Detection Meets Collider Physics_ , 1807.10261.
  * (44) T. Heimel, G. Kasieczka, T. Plehn and J. M. Thompson, _QCD or What?_ , _SciPost Phys._ 6 (2019) 030, [1808.08979].
  * (45) M. Farina, Y. Nakai and D. Shih, _Searching for New Physics with Deep Autoencoders_ , 1808.08992.
  * (46) O. Cerri, T. Q. Nguyen, M. Pierini, M. Spiropulu and J.-R. Vlimant, _Variational Autoencoders for New Physics Mining at the Large Hadron Collider_ , 1811.10276.
  * (47) T. S. Roy and A. H. Vijay, _A robust anomaly finder based on autoencoder_ , 1903.02032.
  * (48) T. Roxlo and M. Reece, _Opening the black box of neural nets: case studies in stop/top discrimination_ , 1804.09278.
  * (49) J. Brehmer, K. Cranmer, G. Louppe and J. Pavez, _Constraining Effective Field Theories with Machine Learning_ , _Phys. Rev. Lett._ 121 (2018) 111801, [1805.00013].
  * (50) J. Brehmer, K. Cranmer, G. Louppe and J. Pavez, _A Guide to Constraining Effective Field Theories with Machine Learning_ , _Phys. Rev._ D98 (2018) 052004, [1805.00020].
  * (51) J. Guo, J. Li, T. Li, F. Xu and W. Zhang, _Deep learning for $R$-parity violating supersymmetry searches at the LHC_, _Phys. Rev._ D98 (2018) 076017, [1805.10730].
  * (52) 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].
  * (53) R. T. D’Agnolo and A. Wulzer, _Learning New Physics from a Machine_ , _Phys. Rev._ D99 (2019) 015014, [1806.02350].
  * (54) A. De Simone and T. Jacques, _Guiding New Physics Searches with Unsupervised Learning_ , 1807.06038.
  * (55) C. Englert, P. Galler, A. Pilkington and M. Spannowsky, _Approaching robust EFT limits for CP-violation in the Higgs sector_ , 1901.05982.
  * (56) J. H. Collins, K. Howe and B. Nachman, _Extending the search for new resonances with machine learning_ , _Phys. Rev._ D99 (2019) 014038, [1902.02634].
  * (57) Z. C. Lipton, _The mythos of model interpretability_ , in _2016 ICML Workshop on Human Interpretability in Machine Learning (WHI 2016)_ , 2016.  1606.03490.
  * (58) K. Datta, A. Larkoski and B. Nachman, _Automating the Construction of Jet Observables with Machine Learning_ , 1902.07180.
  * (59) D. Krohn, J. Thaler and L.-T. Wang, _Jet Trimming_ , _JHEP_ 02 (2010) 084, [0912.1342].
  * (60) F. V. Tkachov, _Measuring multi - jet structure of hadronic energy flow or What is a jet?_ , _Int. J. Mod. Phys._ A12 (1997) 5411–5529, [hep-ph/9601308].
  * (61) M. M. Bronstein, J. Bruna, Y. LeCun, A. Szlam and P. Vandergheynst, _Geometric deep learning: Going beyond euclidean data_ , _IEEE Signal Processing Magazine_ 34 (July, 2017) 18–42, [1611.08097].
  * (62) J. Cogan, M. Kagan, E. Strauss and A. Schwarztman, _Jet-Images: Computer Vision Inspired Techniques for Jet Tagging_ , _JHEP_ 02 (2015) 118, [1407.5675].
  * (63) S. D. Ellis, C. K. Vermilion and J. R. Walsh, _Techniques for improved heavy particle searches with jet substructure_ , _Phys. Rev._ D80 (2009) 051501, [0903.5081].
  * (64) S. D. Ellis, C. K. Vermilion and J. R. Walsh, _Recombination Algorithms and Jet Substructure: Pruning as a Tool for Heavy Particle Searches_ , _Phys. Rev._ D81 (2010) 094023, [0912.0033].
  * (65) 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–177, [1410.3012].
  * (66) J. Bellm et al., _Herwig 7.0/Herwig++ 3.0 release note_ , _Eur. Phys. J._ C76 (2016) 196, [1512.01178].
  * (67) M. Bahr et al., _Herwig++ Physics and Manual_ , _Eur. Phys. J._ C58 (2008) 639–707, [0803.0883].
  * (68) P. T. Komiske, E. M. Metodiev and J. Thaler, _Energy flow polynomials: A complete linear basis for jet substructure_ , _JHEP_ 04 (2018) 013, [1712.07124].
  * (69) J. Alwall, R. Frederix, S. Frixione, V. Hirschi, F. Maltoni, O. Mattelaer et al., _The automated computation of tree-level and next-to-leading order differential cross sections, and their matching to parton shower simulations_ , _JHEP_ 07 (2014) 079, [1405.0301].
  * (70) C. Degrande, B. Fuks, V. Hirschi, J. Proudom and H.-S. Shao, _Automated next-to-leading order predictions for new physics at the LHC: the case of colored scalar pair production_ , _Phys. Rev._ D91 (2015) 094005, [1412.5589].
  * (71) “FeynRules models to be used for NLO calculations with<EMAIL_ADDRESS>https://feynrules.irmp.ucl.ac.be/wiki/NLOModels.
  * (72) A. Alloul, N. D. Christensen, C. Degrande, C. Duhr and B. Fuks, _FeynRules 2.0 - A complete toolbox for tree-level phenomenology_ , _Comput. Phys. Commun._ 185 (2014) 2250–2300, [1310.1921].
  * (73) C. Degrande, C. Duhr, B. Fuks, D. Grellscheid, O. Mattelaer and T. Reiter, _UFO - The Universal FeynRules Output_ , _Comput. Phys. Commun._ 183 (2012) 1201–1214, [1108.2040].
  * (74) DELPHES 3 collaboration, 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_ , _JHEP_ 02 (2014) 057, [1307.6346].
  * (75) M. Cacciari, G. P. Salam and G. Soyez, _FastJet User Manual_ , _Eur. Phys. J._ C72 (2012) 1896, [1111.6097].
  * (76) M. Cacciari and G. P. Salam, _Dispelling the $N^{3}$ myth for the $k_{t}$ jet-finder_, _Phys. Lett._ B641 (2006) 57–61, [hep-ph/0512210].
  * (77) M. Cacciari, G. P. Salam and G. Soyez, _The Anti-k(t) jet clustering algorithm_ , _JHEP_ 04 (2008) 063, [0802.1189].
  * (78) B. R. Webber, _QCD Jets and Parton Showers_ , in _Quantum chromodynamics and beyond: Gribov-80 memorial volume. Proceedings, Memorial Workshop devoted to the 80th birthday of V.N. Gribov, Trieste, Italy, May 26-28, 2010_ , pp. 82–92, 2011.  1009.5871.  DOI.
  * (79) B. Bhattacherjee, S. Mukhopadhyay, M. M. Nojiri, Y. Sakaki and B. R. Webber, _Associated jet and subjet rates in light-quark and gluon jet discrimination_ , _JHEP_ 04 (2015) 131, [1501.04794].
  * (80) A. Krogh and J. A. Hertz, _A simple weight decay can improve generalization_ , in _Advances in Neural Information Processing Systems 4_ (J. E. Moody, S. J. Hanson and R. P. Lippmann, eds.), pp. 950–957.  Morgan-Kaufmann, 1992.
  * (81) D. P. Kingma and J. Ba, _Adam: A method for stochastic optimization_ , in _The International Conference on Learning Representations (ICLR)_ , 2015\.  1412.6980.
  * (82) X. Glorot and Y. Bengio, _Understanding the difficulty of training deep feedforward neural networks_ , in _Proceedings of the Thirteenth International Conference on Artificial Intelligence and Statistics_ (Y. W. Teh and M. Titterington, eds.), vol. 9 of _Proceedings of Machine Learning Research_ , (Chia Laguna Resort, Sardinia, Italy), pp. 249–256, PMLR, 13–15 May, 2010.
  * (83) F. Chollet et al., “Keras.” https://keras.io, 2015.
  * (84) M. Abadi, A. Agarwal, P. Barham, E. Brevdo, Z. Chen, C. Citro et al., “TensorFlow: Large-scale machine learning on heterogeneous systems.” https://www.tensorflow.org/, 2015.
  * (85) D. Alvarez Melis and T. Jaakkola, _Towards robust interpretability with self-explaining neural networks_ , in _Advances in Neural Information Processing Systems 31_ (S. Bengio, H. Wallach, H. Larochelle, K. Grauman, N. Cesa-Bianchi and R. Garnett, eds.), pp. 7786–7795.  Curran Associates, Inc., 2018.  1806.07538.
  * (86) D. Clevert, T. Unterthiner and S. Hochreiter, _Fast and Accurate Deep Network Learning by Exponential Linear Units (ELUs)_ , in _The International Conference on Learning Representations (ICLR)_ , 2016.  1511.07289.
  * (87) K. He, X. Zhang, S. Ren and J. Sun, _Delving deep into rectifiers: Surpassing human-level performance on imagenet classification_ , in _The IEEE International Conference on Computer Vision (ICCV)_ , December, 2015.  1502.01852.
  * (88) M. Tulio Ribeiro, S. Singh and C. Guestrin, _Model-Agnostic Interpretability of Machine Learning_ , in _2016 ICML Workshop on Human Interpretability in Machine Learning (WHI 2016)_ , 2016.  1606.05386.
  * (89) R. D. Ball et al., _Parton distributions with LHC data_ , _Nucl. Phys._ B867 (2013) 244–289, [1207.1303].
  * (90) P. Skands, S. Carrazza and J. Rojo, _Tuning PYTHIA 8.1: the Monash 2013 Tune_ , _Eur. Phys. J._ C74 (2014) 3024, [1404.5630].
  * (91) J. Bellm et al., “Minimum-bias and underlying-event tunes.” https://herwig.hepforge.org/tutorials/mpi/tunes.html, 2015.
  * (92) S. Gieseke, C. Rohr and A. Siodmok, _Colour reconnections in Herwig++_ , _Eur. Phys. J._ C72 (2012) 2225, [1206.0041].