Interaction networks for the identification of boosted $H\to b\overline{b}$ decays
Eric A. Moreno Thong Q. Nguyen Jean-Roch Vlimant Olmo Cerri Harvey B. Newman Avikar Periwal Maria Spiropulu Javier M. Duarte Maurizio Pierini


FERMILAB-PUB-19-492-CMS-E










California Institute of Technology, Pasadena, California 91125, USA

jduarte@ucsd.edu
University of California San Diego, La Jolla, California 92093, USA
Fermi National Accelerator Laboratory, Batavia, Illinois 60510, USA

European Organization for Nuclear Research, 1211 Geneva 23, Switzerland



  We develop an algorithm based on an interaction network to identify high-transverse-momentum Higgs bosons decaying to bottom quark-antiquark pairs and distinguish them from ordinary jets that reflect the configurations of quarks and gluons at short distances.
  The algorithm's inputs are features of the reconstructed charged particles in a jet and the secondary vertices associated with them.
  Describing the jet shower as a combination of particle-to-particle and particle-to-vertex interactions, the model is trained to learn a jet representation on which the classification problem is optimized.
  The algorithm is trained on simulated samples of realistic LHC collisions, released by the CMS Collaboration on the CERN Open Data Portal.
  The interaction network achieves a drastic improvement in the identification performance with respect to state-of-the-art algorithms.


  DOI: https://doi.org/10.1103/PhysRevD.102.01201010.1103/PhysRevD.102.012010
 



Interaction networks for the identification of boosted → decays
    Maurizio Pierini 0000-0003-1939-4268
    December 30, 2023
===============================================================





§ INTRODUCTION




Jets are collimated showers of hadrons that reflect the configurations of quarks and gluons produced at particle colliders.
Each shower, consisting of quarks and gluons emitted by the primary particle, results in an approximately cone-shaped spray of hadrons, which are then observed in particle detectors.
Jet identification, or tagging, algorithms are designed to identify the nature of the primary particle that initiates a shower by studying the collective features of the hadrons inside the jet.

Traditionally, jet tagging was limited to light-flavor quarks (), gluons (), or  quarks.
At the CERN Large Hadron Collider (LHC), jet tagging becomes a more complex task as new jet topologies are accessible (see Fig. <ref>).
Due to the large center-of-mass energy available in LHC collisions, heavy particles, such as , , and Higgs bosons () and top quarks (), may be produced with large transverse momentum ().
These particles can decay to all-quark final states.
Due to the large  of the original particle, these quarks are produced within a small solid angle.
The overlapping showers produced by these quarks may be reconstructed as a single massive jet.
As shown in Fig. <ref>, the presence of  quarks in the jet gives rise to unique experimental signatures.
In particular,  hadrons are characterized by a lifetime of approximately 1.5 ps, which results in a detectable displacement between the proton-collision point and the point where the  hadron decays.


The identification of jets from heavy resonances relies on jet substructure techniques, designed to quantify the number of clusters of energetic particles, or prongs, inside the jet.
The study of jet substructure was pioneered in the 1990s and early 2000s <cit.>, but interest skyrocketed after its proposed application to reconstruct Higgs bosons when produced in association with a vector boson <cit.>.
Extensive reviews of these techniques are provided in Refs. <cit.>.
Additional discrimination is provided by the reconstructed jet mass, usually computed after a jet grooming algorithm.
A review of the techniques used to reconstruct jets and their substructure at the LHC experiments can be found in Ref. <cit.>.
The jet mass plays a special role in physics analyses exploiting jet substructure, as described for instance in Ref. <cit.>.
The jet mass distribution is typically used to separate jets from boosted heavy particles, characterized by a peaking distribution, from the smoothly falling background, due to ordinary quark and gluon jets.
For certain applications, it is desirable to avoid any distortion of the jet mass distribution when applying a jet-tagging selection.




Due to its lifetime, the presence of a  hadron inside of a jet typically results in a reconstructed secondary vertex (SV) that is displaced from the primary vertex (PV).
Modern particle detectors are equipped with a vertex detector that can accurately determine SV positions and their separation from the PV, even in a dense environment like a high- jet.
This feature is particularly important for tagging a Higgs boson decaying to a bottom quark-antiquark pair (→) because all of the relevant jet constituents originate from two displaced vertices.

Recently, several approaches based on deep learning (DL) have been proposed to optimize jet tagging algorithms (see Sec <ref>), both using expert features with dense layers or raw data representations (e.g., images or lists of particle properties) with more complex architectures.
For instance, the LHC collaborations and other researchers have investigated the optimal way to combine substructure, tracking, and vertexing information to enhance the tagging efficiency for high- → decays <cit.>.
This is an important task in particle physics because measurements of high- → decays may help resolve the loop-induced and tree-level contributions to the gluon fusion process, providing an complementary approach to study the  Yukawa beyond the  process <cit.>.
These measurements are also sensitive probes for physics beyond the standard model <cit.>.
Finally, improving these measurements is important for measuring the Higgs boson self-coupling through the production of → <cit.>.

While existing DL approaches have been successfully applied to jet tagging, particle jets involve multiple entities with complex interactions that are not easily encoded as images or lists.
Graphs provide a natural representation for such relational information.
Traditional machine learning methods use feature engineering and preprocessing to learn from these graphs, which can be time consuming and costly, and may miss important features present in the data.
Graph representation learning, including graph convolution networks <cit.> and graph generative models <cit.>, leverages DL to learn directly from graph-structured data.
In contrast to other DL methods, graph representation learning can (1) handle irregular grids with non-Euclidean geometry <cit.>, (2) encode physics knowledge via graph construction <cit.>, and (3) introduce relational inductive bias into data-driven learning systems <cit.>.
For example, while convolutional neural networks (CNNs) are powerful classifiers that work extremely well for data represented on a grid <cit.>, geometric DL algorithms, such as graph neural networks (GNNs) <cit.>, are applicable even without an underlying grid structure.
Because the data in many scientific domains are not Euclidean, GNNs emerge as a more natural choice.

In this work, we propose to identify → jets with an interaction network (IN), a type of graph network.
In Ref. <cit.>, INs were introduced to describe complex physical systems and predict their evolution after a certain amount of time.
This was achieved by constructing graph networks to learn the interactions between the physical objects, represented as the nodes of the graph.

Just as noted jet substructure variables like N^β=1_2 and D^β=1_2 compute 2-point energy correlation functions between jet constituents to quantify the number of prongs in a jet <cit.>, we posit that the ability of INs to learn complex pairwise relationships aids in identifying the patterns present in → decays.
Moreover, Ref. <cit.> showed that the learned features of an IN correlate with known jet substructure variables.
It was further demonstrated that the IN architecture outperformed other deep neural networks (DNNs), such as dense, convolutional, and recurrent networks, for a jet-substructure classification task.
However, this study was limited because the simulation considered was not fully realistic.

In this paper, we demonstrate that an interaction network with an extended feature representation outperforms state of the art methods for → tagging with -based <cit.> realistic simulation, while relying on less parameters.
In particular, we investigate the use of INs to learn a collective representation of the tracking, vertexing, and substructure properties of the jet and employ this optimized representation to enhance the tagging efficiency.
By placing charged particles and secondary vertices on a graph, the network can learn a representation of each particle-to-particle and particle-to-vertex interaction, and exploit this information to categorize a given jet as signal (→) or background (QCD).

The study is carried out using a sample of fully simulated LHC collision events, released by the CMS Collaboration on the CERN Open Data portal <cit.>.
Previously, many machine learning studies were limited to studies based on generator-level physics with simple detector emulation.
The released CMS full-simulation samples allow for a more in depth and realistic study of the efficacy of machine learning methods on high-energy physics experiments.
We compare the performance to several different algorithms that we trained with open simulation for → tagging based on the architecture of the deep double- (DDB) tagger created by the CMS Collaboration <cit.>.

The IN and DDB taggers only rely on information related to charged particles, which (unlike neutral particles) can be traced back to their point of origin: the PV of the high- collision, any SV generated in the collision, or additional PVs originating from simultaneous proton-proton interactions (pileup).
This choice makes the algorithm particularly robust against the large pileup contamination expected in future LHC runs since this contamination can be removed via so-called charged hadron subtraction (CHS) <cit.>.
For the IN tagger, we consider an extended representation of each charged particle (secondary vertex), with 22 (12) additional features with respect to the nominal DDB tagger (as discussed in Sec. <ref>).
To enable a fair comparison between network architectures, we also report results for an extended variant of the DDB tagger, the deep double-+ (DDB+) tagger, which consumes the same information as the IN tagger.


This paper is structured as follows: we discuss related work in Sec. <ref>.
Section <ref> gives a brief description of the datasets used.
Sections <ref> and <ref> describe the IN architecture and the algorithms used to decorrelate its score from the jet mass distribution.
Section <ref> describes our reconstruction and training of the DDB and DDB+ algorithms.
Results are presented in Sec. <ref> and conclusions are given in Sec. <ref>.



§ RELATED WORK



The use of DNNs has recently found a great deal of success in particle physics <cit.>, especially jet tagging.

Driving this innovation are increasingly complex architectures that are tailored to particular domains, including CNNs <cit.>, which are well suited to computer vision, and recurrent neural networks (RNNs) <cit.> like long short-term memory units (LSTMs) <cit.> and gated recurrent units (GRUs) <cit.>, which are appropriate for natural language processing.

Several classification algorithms have been studied in the context of jet tagging at the LHC using CNNs <cit.> and physics-inspired DNN models <cit.>.
Recurrent and recursive layers have been used to define jet classifiers starting from a list of reconstructed particle momenta <cit.>.
Recently, several different approaches, applied to the specific case of  jet identification have been compared <cit.> on a public  jet tagging dataset <cit.>.
This study found ParticleNet <cit.>, a GNN based on the dynamic graph CNN <cit.> to be the best performing for that task.
In Ref. <cit.>, it was shown that the area under the receiver operating characteristic (ROC) curve (AUC), accuracy, and background rejection at a 30% true positive rate (TPR) of a simple IN architecture trained with the same dataset is within 1%, 0.5%, and 40% of those of ParticleNet, while using 70% fewer parameters.
Unsupervised, semisupervised, and weakly supervised methods have also been proposed, mainly to tag  jets or jets coming from postulated new particles <cit.>.
Finally, others have also explored the CMS open data and simulation to study jet properties and jet classification algorithms in a realistic setting <cit.>.

For the task of identifying → specifically, several machine learning approaches have been applied.
In generator-level studies, Ref. <cit.> uses images, representing both the  candidate jet and the full event, as inputs to a CNN.
In conditions more closely resembling real data, the CMS Collaboration created a boosted decision tree based on expert chosen features to identify the presence of two  hadrons within a single anti- <cit.> R=0.8 jet (AK8 jet) <cit.>.
This approach was extended using a deep neural network and additional particle-level and vertex level information, the DDB tagger <cit.>.
Other more generic CMS algorithms, also based on deep neural networks and known as the boosted event shape tagger (BEST) and the DeepAK8 tagger, were created to classify the decays of multiple heavy resonances, including , , , and  <cit.>.
The ATLAS collaboration has also designed an algorithm to identify two  hadrons within an anti- R=1 jet using  tagging of track-based subjets <cit.>.
For the task of → identification, the CMS DDB tagger, DeepAK8 algorithm, and the ATLAS tagger achieve similar state-of-the-art performance.

Graph networks <cit.> and the related particle flow networks <cit.> have recently been used for other kinds of jet tagging, matching or exceeding the performances of other DL approaches, for event classification <cit.>, for charged particle tracking in a silicon detector <cit.>, for mitigation of the effects pileup <cit.>, and for particle reconstruction in irregular calorimeters <cit.> and the IceCube experiment <cit.>.


While applying GNNs is natural for particle physics data, one issue we confront in this paper is how to deal with heterogeneous hierarchical data, i.e. data composed of different sets of elements with different numbers and types of features.
The primary original contributions of this paper are (1) designing an IN with data comprising a heterogeneous graph with two types of graph nodes: particles and SVs), (2) demonstrating that an IN achieves competetive performance on public, realistic simulation for the task of → tagging with fewer trainable parameters in a way that is robust to the effects of pileup, and (3) comparing and evaluating mass decorrelation methods.



§ DATA SAMPLES



The CMS open data and simulation are available from the CERN Open Data Portal <cit.>, including releases of 2010, 2011, and 2012 CMS collision data as well as 2011, 2012, and 2016 CMS simulated data.

Samples of → jets are available from simulated events containing Randall-Sundrum gravitons <cit.> decaying to two Higgs bosons, which subsequently decay to  pairs.
The event generation was done by the CMS Collaboration with  2.2.2 at leading order, with graviton masses ranging between 0.6 and 4.5.
Generation of this process enables better sampling of events with large Higgs boson .
The main source of background originates from multijet events.
The background dataset was generated with  8.205 <cit.> in different bins of the average  of the final-state partons (p̂_T).
The parton showering and hadronization was performed with  8.205 <cit.>, using the CMS underlying event tune CUETP8M1 <cit.> and the NNPDF 2.3 <cit.> parton distribution functions.
Pileup interactions are modeled by overlaying each simulated event with additional minimum bias collisions, also generated with  8.205.
The CMS detector response is modeled by  <cit.>.

The outcome of the default CMS reconstruction workflow is provided in the open simulation <cit.>.
In particular, particle candidates are reconstructed using the particle-flow (PF) algorithm <cit.>.
Charged particles from pileup interactions are removed using the CHS algorithm.
Jets are clustered from the remaining reconstructed particles using the anti- algorithm <cit.> with a jet-size parameter R=0.8.
The standard CMS jet energy corrections are applied to the jets.
In order to remove soft, wide-angle radiation from the jet, the soft-drop (SD) algorithm <cit.> is applied, with  angular exponent β = 0, soft cutoff threshold z_cut < 0.1, and characteristic radius R_0 = 0.8 <cit.>.
The SD mass () is then computed from the four-momenta of the remaining constituents.

A signal → jet is defined as a jet geometrically matched to the generator-level Higgs boson and both  quark daughters.
Jets from QCD multijet events are used to define a sample of fake → candidates.

The dataset is reduced by requiring the AK8 jets to have 300 <  < 2400, |η| < 2.4, and 40 <  < 200.
After this reduction, the dataset consists of 3.9 million → jets and 1.9 million inclusive QCD jets.
Charged particles are required to have > 0.95 and reconstructed secondary vertices (SVs) are associated with the AK8 jet using Δ R = √(Δϕ^2+Δη^2) < 0.8.
The dataset is divided into blocks of features, referring to different objects.
Different blocks are used as input by the models described in the rest of the paper.

The IN uses 30 features related to charged particles (see Table <ref> in Appendix <ref>).
The IN also uses 14 SV features listed in Table <ref>.
The DDB tagger <cit.> uses a subset of the above features (8 features for each particle and 2 features for each SV), chosen to minimize the correlation with the jet mass.
In addition, the DDB tagger uses 27 high-level features (HLF) listed in Table <ref> and first used in a previous version of the algorithm, described in Ref. <cit.>.
To isolate the effects of the different architecture, the DDB+ tagger uses the same inputs as the IN tagger, while retaining the architecture of the DDB tagger.
The charged particles (SVs) are sorted in descending order of the 2D impact parameter significance (2D flight distance significance) and only the first 60 (5) are considered.




§ THE INTERACTION NETWORK MODEL



The IN is based on two input collections comprising N_p particles, each represented by a feature vector of length P, and N_v vertices, each represented by a feature vector of length S.
Although kinematic features of neutral particles could also be taken into account with an additional input graph, we verified that doing so does not significantly improve the performance for this task as shown in Sec. <ref>.
Further, excluding neutral particles has the benefit of improved robustness to pileup.
For a single jet, the input consists of an X and a Y matrix, with sizes P × N_p and S × N_v, respectively.
The X matrix contains the input features (columns) of the charged particles (rows), while the Y matrix contains the input features of the SVs.

A particle graph 𝒢_p is constructed by connecting each particle to every other particle through N_pp = N_p(N_p - 1) directed edges.
Similarly, a particle-vertex graph 𝒢_pv is constructed by connecting each vertex to each particle  through N_pv = N_pN_v directed edges.
As described below, we only consider those edges that are received by particles because the final aggregation is performed over the particles.
These graphs are pictorially represented in Fig. <ref> for the case of three particles and two vertices.
As shown in the figure, the graph nodes and edges are arbitrarily enumerated.
The result of the graph processing is independent of the labeling order, as described below.




For the graph 𝒢_p, a receiving matrix (R_R) and a sending matrix (R_S) are defined, both of size N_p× N_pp.
The element (R_R)_ij is set to 1 when the ith particle receives the jth edge and is 0 otherwise.
Similarly, the element (R_S)_ij is set to 1 when the ith particle sends the jth edge and is 0 otherwise.
For the second graph, the corresponding adjacency matrices R_K (of size N_p× N_vp) and R_V (of size N_v× N_vp) are defined.
In the example of Fig. <ref>, the R_R, R_S, R_K, and R_V matrices would be written as:


    R_R    =
             (pp)_1     (pp)_2     (pp)_3     (pp)_4     (pp)_5     (pp)_6
                      p_1     1     1     0     0     0     0
                      p_2     0     0     1     1     0     0
                      p_3     0     0     0     0     1     1,
    
         R_S    =
             (pp)_1     (pp)_2     (pp)_3     (pp)_4     (pp)_5     (pp)_6
                      p_1     0     0     1     0     1     0
                      p_2     1     0     0     0     0     1
                      p_3     0     1     0     1     0     0,
    
        R_K    =
             (vp)_1     (vp)_2     (vp)_3     (vp)_4     (vp)_5     (vp)_6
                      p_1     1     1     0     0     0     0
                      p_2     0     0     1     1     0     0
                      p_3     0     0     0     0     1     1,
    
         R_V    =
             (vp)_1     (vp)_2     (vp)_3     (vp)_4     (vp)_5     (vp)_6
                      v_1     1     0     1     0     1     0
                      v_2     0     1     0     1     0     1.




    R_R    = 
             (pp)_1     (pp)_2     (pp)_3     (pp)_4     (pp)_5     (pp)_6
                      p_1     1     1     0     0     0     0
                      p_2     0     0     1     1     0     0
                      p_3     0     0     0     0     1     1,
    
         R_S    =
             (pp)_1     (pp)_2     (pp)_3     (pp)_4     (pp)_5     (pp)_6
                      p_1     0     0     1     0     1     0
                      p_2     1     0     0     0     0     1
                      p_3     0     1     0     1     0     0,
    
        R_K    = 
             (vp)_1     (vp)_2     (vp)_3     (vp)_4     (vp)_5     (vp)_6
                      p_1     1     1     0     0     0     0
                      p_2     0     0     1     1     0     0
                      p_3     0     0     0     0     1     1,
    
         R_V    = 
             (vp)_1     (vp)_2     (vp)_3     (vp)_4     (vp)_5     (vp)_6
                      v_1     1     0     1     0     1     0
                      v_2     0     1     0     1     0     1.



Each column of an adjacency matrix corresponds to a directional connection from one particle to another, (pp)_i, or from a vertex and to a particle, (vp)_j.
Column entries that are 1 in a given row in the receiving matrix R_R indicate that the corresponding particle receives that connection.
Likewise, if a column entry is 1 in a given row in the sending matrix R_S, the corresponding particle is the sender for that connection.
Because the fully connected particle graph we consider has no self-connections, i.e. no particle sends and receives the same connection, the rows of R_R and R_S do not share any of the same nonzero column entries.
For the R_R and R_V adjacency matrices, we only consider those connections that are sent to particles because the final aggregation is performed over the particles.
We tested a version of the IN architecture in which we considered connections that are sent to vertices as well and aggregated separately before being processed by the final network, but found no significant improvement.


The data flow of the IN model is pictorially represented in Fig. <ref>.
The input processing starts by creating the 2P× N_pp particle-particle interaction matrix B_pp and the (P+S)× N_vp particle-vertex interaction matrix B_vp defined as:

    B_pp   = [ X · R_R; X · R_S ],
    
      B_vp   = [ X · R_K; Y · R_V ],

where · indicates the ordinary matrix product.
Each column of B_pp consists of the 2P features of the sending and receiving nodes of each particle-particle interaction, while each column of B_vp consists of the P+S features of each particle-vertex one.

Processing each column of B_pp by the function f_R^pp, one builds an internal representation of the particle-particle interaction with a function f_R^pp: ℝ^2P↦ℝ^D_E, where D_E is the size of the internal representation.
This results in an effect matrix E_pp with dimensions D_E× N_pp.
We similarly build the E_vp matrix, with dimensions D_E× N_vp, using a function f_R^vp:ℝ^P+S↦ℝ^D_E.




We then propagate the particle-particle interactions back to the particles receiving them, by building  E_pp = E_pp R_R^⊤ with dimension D_E × N_p.
We also build E_vp = E_vp R_V^⊤ with dimension D_E × N_p, which collects the information of the particle-vertex interactions for each particle and across all of the vertices.

The next step consists of building the C matrix, with dimensions (P + 2D_E)× N_p, by combining the input information for each particle (X) with the learned representation of the particle-particle (E_pp) and particle-vertex (E_vp) interactions:

    C    = [    X; E_pp; E_vp ].

The final aggregator combines the input and interaction information to build the postinteraction representation of the graph, summarized by the matrix O, with dimensions D_O × N_p.
The aggregator consists of a function f_O: ℝ^P + 2D_E↦ℝ^D_O, which computes the elements of the O matrix
The elements of the O matrix are computed by a function f_O: ℝ^P + 2D_E↦ℝ^D_O, which returns the postinteraction representation for each of the input nodes.
As is done for f_R^pp and f_R^vp, f_O is applied to each column of C.

We stress the fact that the by-column processing applied by the f_R^pp,  f_R^vp, and f_O functions and the sum across interactions by defining the E_pp and E_vp matrices are essential ingredients to make the outcome of the IN tagger independent of the order used to label the N_p input particles and N_v input vertices.
In other words, while the representations of the R_R, R_S, R_K, and R_V matrices depend on the adopted labeling convention, the final representation of each particle does not.

The learned representation of the post-interaction graph, given by the elements of the O matrix, can be used to solve the specific task at hand.
Depending on the task, the final function that computes the classifier output may be chosen to preserve the permutation invariance of the input particles and vertices.
In this case, we first sum along each row (corresponding to a sum over particles) of O to produce a feature vector O with length D_O for the jet as a whole.
This is passed to a function ϕ_C: ℝ^D_O↦ℝ^N, which produces the output of the classifier.


The training of the IN is performed with the CMS open simulation with 2016 conditions.
The input dataset is split into training, validation, and test samples with percentages of 80%, 10%, and 10%, respectively.

We use PyTorch <cit.> to implement and train the classifier on one NVIDIA GeForce GTX 1080 GPU.
We also convert the interaction network into a TensorFlow model, as discussed in Appendix <ref>.
The model is implemented with each of f^pp_R and f^vp_R expressed as a sequence of 3 dense layers of sizes (60, 30, 20) with a rectified linear unit (ReLU) activation function after each layer.
The function f_O is a similar sequence of dense layers of sizes (60, 30, 24) with ReLU activations.
We use up to N_p = 60 charged particles and N_v = 5 secondary vertices as inputs to the IN tagger.
Given the size of these layers, the total number of trainable parameters is 18,144.
We train the model using the Adam optimizer <cit.> with an initial learning rate of 10^-4 and a batch size of 128 for up to 200 epochs, enforcing early stopping <cit.> on the validation loss with a patience of 5 epochs.
The size of the batch is constrained by the required memory utilization of the GPU.
The training takes approximately 25 minutes per epoch on the GPU and stopped after 110 epochs.

For the baseline algorithm, we minimize the categorical cross-entropy loss function for this classification task L_C and let the network exploit all of the discriminating information in the dataset.

To determine the impact of neutral particles, we also train an augmented all-particle IN model, which consumes an additional input set with 10 kinematic features for up to 100 charged or neutral particles, listed in Table <ref>.
This additional input set is processed by the model in a similar way to the SV input set: the set of all particles is fully connected to the set of charged particles. 
The effect matrix for these interactions is computed by an independent neural network and then appended to an enlarged C matrix, now of size (P+3D_E)× N_p, before being processed by the network f_O.
The remaining steps of the model proceed as described above.
The total number of trainable parameters for this model is 24,254.



§ DECORRELATION WITH THE JET MASS


Many possible applications of a jet tagging algorithm would require the final score to be uncorrelated from the jet mass, so that a selection based on the tagger score does not change the jet mass distribution. This is particularly relevant for the background distribution, but is required to some extent also for the signal one.
Several techniques exist to deliver a tagger with minimal effects on the jet mass distribution.
For taggers based on high-level features, one could remove those features more correlated to the jet mass or divide those correlated features by the jet mass.
For taggers based on a more raw representation of the jet (as in this case), one could perform an adversarial training <cit.>.
One could also reweight or remove background events such that the background  distribution is indistinguishable from the signal  distribution <cit.>.
Finally, one could also define a mass-dependent threshold based on simulation as in the “designing decorrelated taggers” (DDT) procedure proposed in Ref. <cit.>.

We test and compare all three methods in App. <ref>.
We found the DDT method to be the most robust and performant deocorrelation procedure.
As such, we use it as the nominal decorrelation method in the following results.



 §.§ Designing decorrelated taggers



Following the DDT procedure <cit.>, the tagger threshold for a given false positive rate (FPR) or “working point” is determined as a function of .
By creating a -dependent tagger threshold, the background jet  distribution for events passing and failing this threshold can be made identical.
In practice, this is done by considering the distribution of the network score versus the jet  for the training dataset.
A quantile regression was used to find the threshold on the network score as a function of  distribution that would correspond to a fixed quantile (the chosen 1-FPR value).
By construction, this procedure results in near-perfect mass decorrelation.

In this case, a gradient boosted regressor <cit.> with the following parameters was used:

    
  * α-quantile of 1 - FPR,
    
  * number of estimators of 500,
    
  * minimum number of samples at a leaf node of 50,
    
  * minimum number of samples to split an internal node of 2500,
    
  * maximum depth of 5,
    
  * validation set of 20%,
    
  * early stopping with tolerance of 10.




§ DEEP DOUBLE-B TAGGER MODELS




The DDB tagger is a convolutional and recurrent neural network model developed by CMS <cit.> to identify boosted → jets.
We reconstruct this model based on publicly available information from the CMS Collaboration as follows.
The model takes as input 27 HLFs used in Ref. <cit.>, as well as 8 particle-specific features of up to 60 charged particles, and 2 properties of up to 5 SVs associated with the jet (see Appendix <ref>).
Each block of inputs is treated as a one-dimensional list, with batch normalization <cit.> applied directly to the input layers.
For each collection of charged particles and SVs, separate 1D convolutional layers <cit.>, with a kernel size of 1, are applied: 2 hidden layers with 32 filters each and ReLU <cit.> activation.
The outputs are then separately fed into two gated recurrent units (GRUs) with 50 output nodes each and ReLU activations.
Finally, the GRU outputs are concatenated with the HLFs and processed by a dense layer with 100 nodes and ReLU activation, and another final dense layer with 2 output nodes with softmax activation.
Dropout <cit.> (with a rate of 10%) is used in each layer to prevent overfitting.
The nominal DDB tagger model has 40,344 trainable parameters, 32% of which are found in the fully connected layers.

We define a variant of this model, the DDB+ model, which takes as input all 30 features of charged particles and all 14 features of the SVs.
In this variant, we do not consider the HLFs.
Thus, the final dense layer only receives the GRU outputs from processing the low-level charged particle and SV information.
This extended DDB+ tagger algorithm has 38,746 trainable parameters.
The number of parameters is less overall because the increase in the size of the convolutional and recurrent layers is compensated by the decrease in the size of the fully connected layers.







We train the DDB and DDB+ models using the CMS open simulation dataset with Keras <cit.> over up to 200 epochs with an early stopping patience of 5 epochs and a batch size of 4096 using the Adam optimizer with an initial learning rate of 10^-3.
For both models, one training epoch takes about 3 minutes and training stops after approximately 50 epochs.
In this case, the larger batch size is possible due to the smaller GPU memory utilization of the model during training.
We find consistent performance for different batch size choices with no evidence of overfitting with larger batch sizes.


In order to decorrelate the tagger output from the jet mass, we use the same DDT procedure described in Sec. <ref> applied to both the DDB and DDB+ taggers.






§ RESULTS



In Fig. <ref> the performance of the IN, all-particle IN, DDB, and DDB+ algorithms are quantified in a ROC curve.
The axes are the TPR, or → tagging efficiency and the false positive rate, or QCD mistagging rate.
As shown in Fig. <ref>, the IN provides an improved performance with respect to the DDB and DDB+ taggers.
At a 1% FPR, the IN tagger outperforms the DDB and DDB+ taggers by 37% and 2% in TPR, respectively.
 Likewise, at a 50% TPR, the IN tagger yields a factor of 6 or 1.2 better background rejection (1/FPR) than the DDB or DDB+ tagger, respectively.
Thus, while the additional inputs provide a significant improvement for the DDB+ model, the IN architecture is also important to achieve a better performance with significantly less parameters than the DDB+ model.

We verified that one could match the performance obtained by the IN with a DDB-inspired architecture and expanding the model size.
With 150,786 trainable parameters, a DDB architecture achieves the same performance as the IN at the cost of 8 times more parameters.

Because of this the IN model holds an advantage in terms of memory usage during inference over this alternative model.

Figure <ref> also shows that there is only a modest improvement in the AUC and accuracy by including information in the IN model from neutral particles. 
For this reason and to preserve robustness to increased pileup, in the following results, we consider the original IN model that excludes neutral particles.

Figure <ref> shows an illustration of how the signal and background jet mass distributions change after applying a threshold on the different baseline and DDT-decorrelated tagger scores.
Following Ref. <cit.>, we quantify the impact of these algorithms on the mass decorrelation by computing the Jensen-Shannon (JS) divergence:

    D_JS(P ∥ Q)    = 1/2D_KL(P ∥ M) + 1/2D_KL(Q∥ M),

where M= 1/2(P+Q) is the average of the normalized  distributions of the background jets passing (P) and failing (Q) a given tagger score and D_KL(P ∥ Q  ) = ∑_i P_i log (P_i/Q_i) is the Kullback-Leibler (KL) divergence.
Larger values of the metric 1/D_JS correspond to a better decorrelation.



After applying the mass decorrelation techniques, the performance of each of the taggers worsens slightly but the IN algorithm still significantly outperforms the DDB and DDB+ taggers.
Figure <ref> displays the trade-off between the background rejection and 1/D_JS at different TPRs for the baseline and DDT-decorrelated algorithms.
At a 50% TPR, the decorrelated IN algorithm achieves a significantly better 1/D_JS by a factor of about 2,200 while the background rejection decreases by a factor of about 3.3 compared to the baseline IN algorithm.
At a 1% FPR, the DDT-decorrelated IN tagger has a TPR of 75.6% compared to the DDT-decorrelated DDB (DDB+) tagger with a 55.9% (72.9%) TPR, corresponding to an improvement of 35% (4%).
Table <ref> summarizes different performance metrics for the three considered models and their decorrelated versions.
For the DDT models, the corresponding accuracy is listed for the tagger after the decorrelation is performed for a FPR of 50%.


To quantify the dependence on the number of pileup interactions, Fig. <ref> shows the performance of the different algorithms as a function of the number of primary vertices in the event, which scales linearly with the number of pileup collisions.
Using only charged particles and secondary vertices as input, the IN tagger is robust against an increasing number of pileup interactions, exhibiting behavior similar to the DDB and DDB+ taggers.






§ CONCLUSIONS



We presented a novel technique using a graph representation of the jet's constituents and secondary vertices based on an interaction network to identify Higgs bosons decaying to bottom quark-antiquark pairs (→) in LHC collisions.
This model can operate on a variable number of jet constituents and secondary vertices and does not depend on the ordering schemes of these objects.
The interaction network was trained on an open simulation dataset released by the CMS Collaboration in the CERN Open Data Portal.
A significant improvement in performance is observed with respect to two alternative taggers based on the deep double- tagger created by the CMS Collaboration.
By design, the interaction network uses extended low-level input features for particles and vertices, offers a more flexible representation of jet data, and is robust against the noise generated by pileup collisions.
Even when trained with the same set of input features, the interaction network architecture outperforms the deep double- architecture.
Thus, while part of the improvement is due to the extended input representation, additional improvement comes from the interaction network architecture, despite using on half as many parameters.
The interaction network algorithm implementation and its training code are available at Ref. <cit.>.

Together with the best-performing models, we presented additional models, obtained by applying different decorrelation techniques between the network score and the jet-mass distribution.
This was done to minimize the selection bias of the classifier output towards any values of the jet mass, which would make the algorithms suitable for physics analyses relying on the jet mass as a discriminating variable.
As expected, the decorrelation procedure results in a reduction of the → identification performance.
Nevertheless, the decorrelated interaction network model outperforms the decorrelated deep double- models.

Once applied to a full data analysis, this graph-based tagging algorithm could contribute a substantial improvement to the experimental precision of → measurements, including those sensitive to beyond the standard model physics and the Higgs boson self-coupling.
These results motivate further exploration of applications based on interaction networks (and graph neural networks in general) for object tagging and other similar tasks in experimental high energy physics.

 

This work was possible thanks to the commitment of the CMS Collaboration to release its simulation data through the CERN Open Data Portal.
We would like to thank our CMS colleagues and the CERN Open Data team for their effort to promote open access to science.
In particular, we thank Kati Lassila-Perini for her precious help.

We are grateful to Caltech and the Kavli Foundation for their support of undergraduate student research in cross-cutting areas of machine learning and domain sciences.
This work was conducted at “iBanks,” the AI GPU cluster at Caltech.
We acknowledge NVIDIA, SuperMicro and the Kavli Foundation for their support of iBanks.
This project is partially supported by the United States Department of Energy, Office of High Energy Physics Research under Caltech Contract No. DE-SC0011925.
M. P. is supported by the European Research Council (ERC) under the European Union's Horizon 2020 research and innovation program (Grant Agreement No. 772369).
J. M. D. is supported by Fermi Research Alliance, LLC under Contract No. DE-AC02-07CH11359 with the U.S. Department of Energy, Office of Science, Office of High Energy Physics.
E. A. M. is supported by the Taylor W. Lawrence Research Fellowship and Mellon Mays Research Fellowship.

 



§ APPENDIX A: ADDITIONAL MASS DECORRELATION METHODS



In this appendix, we describe and compare two additional mass decorrelation methods to the DDT procedure described in Section <ref>.
In one method, we train two neural networks simultaneously: the original classifier and an additional network intended to regress the jet mass, known as the adversary.
The original classifer is trained to maximally confuse the adversary.
After training, the effect is the classifier is not able to discriminate the jet mass.
In the other method, we train the classifier with sample weights, such that the QCD background jet mass distribution is reweighted to be identical to that of the → signal.
We then compare these procedures to the DDT method.



 §.§ Adversarial training

The secondary adversary network is constructed that consists of three hidden layers each with 64 nodes.
The adversary is trained simultaneously with the classifier (interaction network) using the summed post-interaction feature vector 𝐎 as its input.
From this input, the adversary is trained to predict a one-hot encoding of the pivot feature , which we aim to decorrelate from the classifier output.
The chosen one-hot encoding corresponds to 40  bins from 40 to 200.
The training begins by initializing the weights from the best classifier training.
The adversary is then pre-trained for 10 epochs using the Adam algorithm with an initial learning rate of 10^-4.
During each epoch, the classifier is first trained by minimizing the total loss

    L    = L_C - λ L_adversary.

Subsequently, the adversary is trained by minimizing L_adversary using only the background QCD samples.
To balance tagging performance and  correlation, λ=10 was chosen.



 §.§ Sample reweighting

While adversarial training requires a complicated tuning process, sample reweighting is a simpler way to achieve the same goal.
Individual QCD events are weighted in the loss function based on their mass bin as to match the signal jet mass distribution of the training sample.
Given a background event in certain mass bin, with the number of background and signal events in that bin denoted as N^bin_b and N^bin_s, respectively, the event is weighted by w_bin=N^bin_s/N^bin_b.









 §.§ Results


Figure <ref> shows a comparison of the ROC curves for the baseline IN algorithm and the versions that were decorrelated using the DDT procedure, adversarial training, and sample reweighting.
Table <ref> summarizes a variety of performance metrics for the decorrelated algorithms including 1/D_JS, which quantifies the success of the decorrelation procedure for a given FPR.
As shown in Fig. <ref> and Table <ref>, the DDT procedure provides the best decorrelation in terms of 1/D_JS, and comparable to the best accuracy, AUC, background rejection, and tagging efficiency.



§ APPENDIX B: MODEL IMPLEMENTED IN ONNX AND TENSORFLOW



In order to integrate the IN algorithm into experimental workflows, it is often necessary to provide the algorithm in other formats.
For example, the CMS experimental software framework CMSSW <cit.> currently only supports ONNX <cit.>, TensorFlow <cit.>, and MXNet <cit.> models.
To perform this conversion, we first translate the PyTorch model into an ONNX representation using the built-in exporter.
Then the conversion to TensorFlow is performed with the dedicated TensorFlow backend for ONNX <cit.>.
The trained model in all three formats is available at Ref. <cit.>.



§ APPENDIX C: DATASET FEATURES



The charged particle features used by the IN and DDB taggers are listed in Table <ref>.
The SV features used by both taggers are listed in Table <ref>, and the high-level features used only by the reconstructed DDB tagger are shown in Table <ref>.
Finally, additional features of charged or neutral particles are listed in Table <ref> to demonstrate the change in the performance of the IN model by including neutral particles.























