# Spatial-aware Transformer-GRU Framework for Enhanced Glaucoma Diagnosis from
3D OCT Imaging

Mona Ashtari-Majlan  Mohammad Mahdi Dehshibi Member, IEEE  David Masip  Senior
Member, IEEE This work received partial funding from the European Research
Council (ERC) under the Horizon 2020 research and innovation program (grant
agreement number 101002711). Additionally, support was provided by the Spanish
Ministry of Science and Innovation through the PID2022-138721NB-I00 grant,
funded by the FEDER initiative. M. A. and D. M. are with the Department of
Computer Science, Multimedia, and Telecommunications, Universitat Oberta de
Catalunya, Barcelona, Spain (e-mail: {mashtarimajlan, dmasipr}@uoc.edu)M. M.
D. is with the Department of Computer Science and Engineering, Universitat
Carlos III de Madrid, Madrid, Spain (e-mail: mohammad.dehshibi@yahoo.com).

###### Abstract

Glaucoma, a leading cause of irreversible blindness, necessitates early
detection for accurate and timely intervention to prevent irreversible vision
loss. In this study, we present a novel deep learning framework that leverages
the diagnostic value of 3D Optical Coherence Tomography (OCT) imaging for
automated glaucoma detection. In this framework, we integrate a pre-trained
Vision Transformer on retinal data for rich slice-wise feature extraction and
a bidirectional Gated Recurrent Unit for capturing inter-slice spatial
dependencies. This dual-component approach enables comprehensive analysis of
local nuances and global structural integrity, crucial for accurate glaucoma
diagnosis. Experimental results on a large dataset demonstrate the superior
performance of the proposed method over state-of-the-art ones, achieving an
F1-score of 93.58%, Matthews Correlation Coefficient (MCC) of 73.54%, and AUC
of 95.24%. The framework’s ability to leverage the valuable information in 3D
OCT data holds significant potential for enhancing clinical decision support
systems and improving patient outcomes in glaucoma management.

Glaucoma detection, 3D Optical Coherence Tomography, Gated Recurrent Units,
Vision Transformer, Spatial Coherence

##  1 Introduction

Glaucoma, a leading cause of irreversible blindness worldwide [30], is
particularly insidious due to its asymptomatic nature in the early stages.
This characteristic underscores the critical importance of timely diagnosis,
as early intervention is vital for preventing visual impairment and
irreversible vision loss [34]. Optical Coherence Tomography (OCT), a three-
dimensional (3D) non-invasive imaging modality, has revolutionized
ophthalmological diagnostics by enabling high-resolution visualization of the
eye’s intricate anatomy [12]. This technology has proven valuable for the
early detection and management of glaucoma.

Despite the widespread adoption of OCT imaging, a significant proportion of
studies on glaucoma diagnosis have primarily focused on analyzing two-
dimensional (2D) OCT scans [5, 6]. These approaches typically concentrate on
the middle slice of the OCT volume, centered around the Optic Nerve Head
(ONH), which has been demonstrated to be effective in characterizing
structural changes related to glaucoma [33]. However, this localized analysis
inherently overlooks the potential diagnostic value contained within the
comprehensive three-dimensional data provided by OCT imaging.

Glaucoma is a multifaceted disease that can manifest progressive structural
alterations throughout the retina, including notable thinning of the Retinal
Nerve Fiber Layer (RNFL) [13]. This widespread impact suggests that analyzing
the total B-scan slices within the OCT volume could reveal important
characteristics and patterns indicative of glaucoma. Consequently, a holistic
analysis that leverages the rich information within the entire ocular
structure is crucial for enhancing diagnostic accuracy and enabling early
intervention[11, 20, 26].

Herein lies the potential of Artificial Intelligence (AI)-based clinical
decision support systems, which can assist practitioners by leveraging
advanced algorithms for automated disease detection and management. To address
the challenges associated with the manual inspection of volumetric OCT data,
in this work, we propose a novel deep learning-based framework designed to
maximize the diagnostic potential of 3D OCT imaging for automated glaucoma
screening. In contrast to traditional analysis techniques that often overlook
valuable information, our approach systematically extracts and integrates
features from the entire OCT scan, revealing subtle glaucomatous indicators
distributed across the ocular structure.

The proposed end-to-end framework comprises two core components: (1) a
Transformer-based model for extracting features from individual slices within
the OCT volumes and (2) a Recurrent Neural Network (RNN) for integrating these
features into a comprehensive representation of the retinal structure. This
dual-component architecture (depicted in Fig. 1) enables the model to capture
both local nuances present in individual slices and the global structural
integrity of the retina while incorporating inter-slice dependencies essential
for a comprehensive analysis.

Figure 1: The schematic of a proposed framework for detecting glaucoma. The
framework uses a large Vision Transformer (ViT-large) feature extractor called
RETFound, which was pre-trained on a large OCT dataset. The symbol $\oplus$
represents the concatenation of spatial states at timestep $i$, combining
information from anterior and posterior spatial locations to capture long-
range spatial dependencies. The framework includes a Bidirectional Gated
Recurrent Unit (GRU), Adaptive Max Pooling (AMP), and a Fully Connected (FC)
layer.

We used a pre-trained large Vision Transformer (ViT-large) backbone developed
by Zhou et al. [35] for feature extraction from individual B-scan slices. This
foundational model, called RETFound, leverages self-supervised learning on 1.6
million unlabeled retinal images to acquire generalizable representations.
Notably, this model exhibits proficiency in identifying intricate patterns and
features relevant to ocular health, ensuring high-precision analysis of
individual slices within 3D OCT images. After feature extraction, the RNN
component, implemented using bidirectional Gated Recurrent Units (GRUs),
integrates these slice-based representations to account for the sequential and
interconnected nature of slices within the three-dimensional space. This
integration is crucial for capturing the dynamic spatial relationships and
dependencies between slices, thereby reflecting the comprehensive structural
complexity of the retina. Experimental results on a 3D OCT dataset [20]
validate the superior performance of our proposed framework in glaucoma
detection compared to state-of-the-art methods.

The following is the organization of the paper: Section 2 reviews the previous
studies. Section 3 describes the proposed method. In Section 4, we present and
discuss the experimental setup and results. Finally, in Section 5, we provide
the conclusion, summarizing our findings and suggesting future research
directions.

##  2 Literature review

Glaucoma is a medical condition characterized by structural changes in the
optic nerve head and thinning of the retinal nerve fiber layer in the retina
[22]. Numerous studies have demonstrated the capability of OCT imaging to
detect these structural alterations, facilitating early diagnosis and
monitoring of the disease [3, 14, 17, 21, 25]. A comprehensive literature
review covering the use of OCT imaging for glaucoma diagnosis is available in
[1].

Numerous previous studies have focused on examining two-dimensional (2D) OCT
scans, also known as B-scans, for glaucoma detection [5]. These studies have
explored various deep learning methods, integrating ophthalmological domain
knowledge to enhance the diagnostic efficacy of OCT imaging. For instance,
Raja et al. [25] developed a deep convolutional autoencoder network to segment
RNFL regions in the retina. The region-aware encoder component of their
proposed model was employed for glaucoma detection, achieving an accuracy of
94.91%, F1-score of 95.77%, and an area under the curve (AUC) of 98.71% on the
AFIO dataset [24]. Additionally, Garcia et al. [6] developed a hybrid neural
network incorporating a VGG16 backbone tailored for processing raw OCT B-scans
and an RNFL descriptor for extracting this region thickness information. Their
approach achieved an accuracy of 87.88% and an F1-score of 81.82% for
discriminating between healthy, early, and advanced glaucoma samples.

The availability of pre-trained deep learning models has facilitated the
transfer of features learned from large-scale 2D datasets, such as ImageNet,
for more accurate classification of 2D OCT scans. However, one of the primary
limitations of relying solely on 2D OCT scans for glaucoma detection is that a
single 2D slice can only capture a static view of the retina, failing to
account for the complex nature of RNFL thinning and damage distribution across
the retinal structure [33]. In an attempt to address this limitation, Wang et
al. [33] proposed selecting every three successive B-scan images from a 3D OCT
volume to form a three-channel input with volume-level labeling. They curated
the Hong Kong dataset, comprising 2,926 glaucomatous and 1,951 normal OCT
scans, and employed it to train their model. This approach achieved an
accuracy of 92.70%, an F1-score of 94.10%, and an AUC of 97.70% for diagnosing
glaucoma. However, their methodology did not fully consider inter-slice
dependencies, a critical aspect of spatial context that could potentially
improve the accuracy of glaucoma detection.

To better capture the intricacies of 3D OCT scans, several studies have
employed three-dimensional convolutional neural network (3D-CNN) models [8,
20, 23, 28]. For instance, Noury et al. [23] proposed a 3D-CNN model with
multiple layers of dense and 3D spatial convolutional blocks. They curated the
Stanford dataset, comprising 1,617 glaucomatous and 844 normal OCT scans, for
training their model. Subsequently, they evaluated their proposed model on
this dataset and three additional external datasets, achieving an AUC of
91.00% for the Stanford dataset, 80.00% for the Hong Kong dataset, 94.00% for
the India dataset, and 87.00% for the Nepal dataset. Furthermore, Ran et al.
[26] proposed training a 3D ResNet34 model from scratch on the Hong Kong
dataset for glaucoma diagnosis, yielding an AUC of 96.90% and accuracy of
91.00% on the primary validation set. George et al. [8] proposed a novel model
for glaucoma diagnosis and VF Index estimation. The model is an attention-
guided 3D-CNN with three streams, where one stream receives raw 3D-OCT scans
as inputs, while the inputs for the other two streams are determined during
training, guided by attention heatmaps from 3D Grad-CAM. The model was trained
on a dataset containing 3,355 glaucomatous and 427 normal OCT scans. The
results showed that the model achieved an AUC of 93.77%, an accuracy of
91.07%, and an F1-score of 94.88% for glaucoma diagnosis.

Training 3D-CNN architectures for analyzing 3D images is computationally
intensive, posing efficiency challenges and the risk of the curse of
dimensionality. To mitigate these issues, some studies in medical imaging have
proposed decomposing 3D images into a sequence of 2D slices [7, 10, 16, 19,
27, 32]. This strategy aims to balance computational burden and depth of
analysis while preserving the crucial inter-slice context for accurate image
classification. Specifically for diagnosing glaucoma, García et al. [7]
proposed a method that combines CNNs and RNNs with Long Short-Term Memory
(LSTM), using a pre-trained VGG16 model to extract discriminative features
from 2D slices of OCT images. They achieved an accuracy of 81.25%, an F1-score
of 78.57%, and an AUC of 80.79% on a dataset containing 144 glaucomatous and
176 normal OCT scans.

Other imaging modalities have also been used to further improve the accuracy
of glaucoma diagnosis. For instance, Song et al. [29] proposed a transformer-
based architecture that uses OCT and visual field (VF) results to diagnose
glaucoma. Their model incorporated deep reasoning to examine implicit pairwise
relations between these two complementary modalities. The experimental
results, conducted on a dataset comprising 697 glaucomatous and 698 normal
scans, demonstrated an accuracy of 88.30%, an F1-score of 88.90%, and an AUC
of 93.90% for glaucoma diagnosis.

Our approach deviates from previous studies in several key aspects, offering
enhanced diagnostic capabilities. Specifically, we employ a transformer-based
feature extractor trained on a large-scale OCT dataset, surpassing the
performance of models pre-trained on ImageNet and fine-tuned on relatively
small OCT datasets. Additionally, using a bidirectional GRU model for
analyzing OCT slices ensures a more comprehensive capture of anterior and
posterior spatial dependencies between slices. Furthermore, our model is
trained on a larger yet imbalanced dataset, better reflecting real-world
clinical scenarios. These distinctions underscore our model’s superior
generalizability and applicability in accurately diagnosing glaucoma through
OCT image analysis.

##  3 Proposed Method

The proposed model is designed for end-to-end glaucoma classification. It
integrates two core components: Feature Extraction and Sequential Processing.
The model employs a visual transformer and an RNN with the gated recurrent
unit to enhance the efficiency and accuracy of glaucoma diagnosis. By
leveraging these techniques, the model aims to harness the diagnostic
potential of 3D OCT scans for reliable glaucoma detection.

###  3.1 Feature Extraction

The Feature Extraction component uses the ViT-large encoder from the RETFound
model [35], a self-supervised autoencoder framework pre-trained on large-scale
OCT data. The RETFound encoder employs a ViT-large architecture [4] with 24
Transformer blocks and a 1024-dimensional embedding vector. We discard the
decoder layers for our 3D OCT scans and leverage the pre-trained ViT-large
encoder to extract rich feature representations.

Mathematically, let $V$ denote a 3D OCT volume, which is decomposed into $D$
individual slices, where the $i$-th slice is represented as
$S_{i}\in\mathbb{R}^{H\times W}$, for $i=1,2,\dots,D$. Each slice $S_{i}$ is
independently processed by the ViT-large model to derive a feature vector
$f_{i}\in\mathbb{R}^{1024}$, thereby extracting rich, slice-specific feature
representations from the 3D OCT data.

###  3.2 Sequential Processing

Following the feature extraction, to encompass all inter-slice correlations
and spatial dependencies intrinsic to OCT imaging, we propose using a
bidirectional GRU architecture for glaucoma classification. The GRU, a variant
of RNNs tailored for sequence modeling, offers an effective framework for
analyzing sequential data. In our methodology, each feature vector $f_{i}$
extracted from slice $S_{i}$ of the OCT scan constitutes a distinct state at
timestep $i$ within the sequence. This bidirectional approach enables the GRU
layers to adeptly capture both anterior and posterior spatial dependencies
among the sequential OCT slices. Such a comprehensive analysis facilitates an
in-depth understanding of the spatial relationships and pathological
indicators critical for accurate glaucoma diagnosis. Moreover, the GRU
architecture efficiently mitigates the vanishing gradient problem and enables
the learning of long-term dependencies with fewer parameters than conventional
long-term memory models, thereby enhancing the model’s efficiency and
effectiveness in glaucoma classification.

To further delineate the sequential processing component of our proposed
framework, the architecture employs two consecutive bidirectional GRU layers,
each with a hidden size of 256. Given a sequence of feature vectors
$\\{f_{1},f_{2},\cdots,f_{D}\\}$ where $D=64$, the bidirectional GRU processes
each $f_{i}$ using Eq. 1.

$\displaystyle h_{i}^{\mathrm{(fwd)}}$
$\displaystyle=\mathrm{GRU}^{\mathrm{(fwd)}}(f_{i},h_{i-1}^{\mathrm{(fwd)}})$
(1) $\displaystyle h_{i}^{\mathrm{(bwd)}}$
$\displaystyle=\mathrm{GRU}^{\mathrm{(bwd)}}(f_{i},h_{i+1}^{\mathrm{(bwd)}})$

where $h_{i}^{\mathrm{(fwd)}}$ and $h_{i}^{\mathrm{(bwd)}}$ denote the forward
and backward hidden states at timestep $i$, respectively. The final hidden
states are concatenated to form a comprehensive sequence representation, which
is expressed by Eq. 2.

$H_{G}=\left[h_{1}^{\mathrm{(fwd)}}\oplus
h_{1}^{\mathrm{(bwd)}},h_{2}^{\mathrm{(fwd)}}\oplus
h_{2}^{\mathrm{(bwd)}},\ldots,h_{D}^{\mathrm{(fwd)}}\oplus
h_{D}^{\mathrm{(bwd)}}\right]$

(2)

where $\oplus$ denotes the concatenation operation. After a Dropout layer,
this sequence $H_{G}$ is then aggregated along the spatial dimension to form a
unified feature representation $h_{\mathrm{pooled}}$ through an Adaptive Max
Pooling (AMP) layer. This layer effectively extracts the essential features
from the sequence, capturing the most important information across all
timesteps. Finally, a Fully Connected (FC) layer maps the aggregated feature
representation to the probability distribution over the two classes, i.e.,
glaucoma and normal, represented in Eq. 3.

$p=\sigma\left(W^{\prime}\cdot h_{\mathrm{pooled}}+b^{\prime}\right)$ (3)

where $p$ denotes the probability vector associated with the classification of
a specific volume $V$, $\sigma$ denotes the Sigmoid activation function,
$W^{\prime}$ and $b^{\prime}$ are the weight and bias of the FC layer,
respectively.

We used the Focal Loss function to train our proposed model, which is a
variant of the binary cross entropy loss that addresses the class imbalance
issue by adjusting the loss contribution from each sample [18]. This method
focuses more on hard-to-classify examples, ultimately reducing the impact of
class imbalance on the model’s learning process. The Focal Loss is defined in
Eq. 4.

$\mathrm{FL}(p_{t})=-\alpha_{t}(1-p_{t})^{\gamma}\log(p_{t})$ (4)

where $p_{t}$ is the model’s estimated probability for the class with label
$t$, $\alpha_{t}$ is a weighting factor to counteract class imbalance, and
$\gamma$ controls the rate at which easy examples are down-weighted. Note that
the details for identifying our model’s optimal set of hyper-parameters are
presented in Section 4.2.

##  4 Experiments

###  4.1 Dataset and Pre-processing

In this study, we used a dataset introduced by Maetschke et al.
[20]111https://doi.org/10.5281/zenodo.1481223, comprising 3D OCT scans
centered on the optic nerve head. The dataset includes scans from 624 patients
obtained using a Cirrus SD-OCT Scanner (Zeiss, Dublin, CA, USA), with physical
dimensions of 6$\times$6$\times$2 mm and a resolution of
64$\times$64$\times$128 voxels per volume. Only scans with a signal strength
of 7 or higher were included, resulting in a total of 1110 scans. These scans
maintain their original laterality, avoiding any left to right eye
alterations. Among these, 263 scans were classified as healthy, and 847 were
diagnosed with primary open-angle glaucoma, based on glaucomatous visual field
defects confirmed by at least two consecutive abnormal test results. The
demographic details like gender and race distribution, along with mean values
and standard deviations for patient’s age, Intraocular Pressure (IOP), Mean
Field Defects (MD), and Glaucoma Hemifield Test (GHT) [2] results are
presented in Table 1. It should be noted that demographic data was incomplete
for some patients, and as a result, the aggregate numbers may not align
exactly with the overall dataset size. For a detailed description of the
dataset, you can refer to the Maestack et al [20].

Table 1: Demographic details of the dataset | Healthy | Glaucoma
---|---|---
Female | 88 | 217
Male | 49 | 215
White | 101 | 318
Black | 30 | 154
Asian | 5 | 12
Age | 54.1 ± 15.3 [22.1-88.9] | 64.3 ± 12.5 [25.2-93.8]
IOP | 13.5 ± 2.4 [9-23] | 16.7 ± 5.8 [2-51]
MD | -0.8 ± 1.7 [-9.9-2.8] | -6.8 ± 8.1 [-32.9-2.17]
GHT | 1.6 ± 1.0 [1-6] | 2.4 ± 0.9 [1-6]

In our study, we employed a 5-fold cross-validation approach to ensure our
model’s robustness and generalizability across different subsets of the data.
The rationale behind using 5-fold cross-validation is to provide a more
reliable estimate of the model’s performance on unseen data by averaging
results across multiple folds. The average number of samples allocated for
training, validation, and testing in each fold is detailed in Table 2.

Table 2: Dataset summary | Total | 5-fold cross validation
---|---|---
Train | Validation | Test
Glaucoma | 843 | 611 | 67 | 169
Normal | 263 | 189 | 21 | 53
Total | 1110 | 800 | 88 | 222

We performed a two-step pre-processing to standardize the input data for
optimal integration with our proposed model. Initially, the scans underwent
normalization using the mean and standard deviation of ImageNet. This step is
crucial for adjusting pixel intensity values to a common scale, thereby
enhancing the model’s ability to learn meaningful patterns from the data.
Subsequently, the scans were resized to a uniform dimension of
64$\times$128$\times$128, ensuring consistency across the dataset and aligning
with the input specifications of the pre-trained components of our model.

###  4.2 Implementation details and Evaluation metrics

The proposed framework is implemented using Python based on the PyTorch
library version 1.8.1 on a computer powered by Intel(R) Xeon(R) Silver 4210R
CPU with 16G of RAM and the NVIDIA GeForce RTX 2080 GPU with 12G of RAM. We
trained the proposed framework using Adam optimizer [15], configured with the
first momentum of 0.9, the second momentum of 0.999, and the stability
constant of $10^{-8}$. To ensure efficient convergence, we explored the impact
of varying combinations of learning rates and batch sizes on model performance
across the $\\{10^{-5},5\times 10^{-5},10^{-4},5\times 10^{-4},10^{-3}\\}$ set
for learning rates and $\\{8,16,32,64,128\\}$ set for batch sizes. The optimal
performance on the validation set was achieved with an initial learning rate
of $10^{-4}$ and a batch size of 32. We employed a learning rate scheduler to
adjust the learning rate by a factor of 0.9 every 5 epochs. To address the
challenge of data imbalance in our dataset, we adopted a balanced batch
training strategy to ensure fair class representation and mitigate bias. The
model was trained for a maximum of 100 epochs and incorporated early stopping
triggered by validation loss stagnation after 6 epochs, which improved
computational efficiency and prevented overfitting. We conducted an empirical
exploration to find the optimal combination of five other hyper-parameters for
improving the performance of the proposed glaucoma classification framework.
The hyper-parameters include the hidden size of two GRU layers, the model
dropout rate, and the $\alpha$ and $\gamma$ parameters in the Focal Loss
function.

One important parameter for determining the capacity of the GRU architecture
to learn and represent data is the hidden size of the GRU layers. Given that
the architecture includes two GRU layers, we evaluated different
configurations for the hidden sizes, ranging from 512, 256, to 128 on the
validation set. The summary of the results is presented in Table 3. Notably,
as the size of the hidden layers decreases from 512 to 128, there is a general
trend toward reduced F1-scores, indicating that larger hidden layer sizes tend
to capture more complex patterns and dependencies in the data. However, the
optimal configuration is achieved with both GRU layers set to a hidden size of
256, resulting in the highest F1-score of 93.51%. This suggests a threshold
for increasing layer size, highlighting the importance of balancing model
capacity and the risk of overfitting.

Table 3: F1-score (%) with respect to the hidden sizes of GRU layers GRU-1 | GRU-2
---|---
512 | 256 | 128 |
512 | 92.11 | 93.29 | 89.34 |
256 | - | 93.51 | 92.91 |
128 | - | - | 91.24 |

Another parameter that helps to improve the GRU architecture generalizability
is the dropout rate. As illustrated in Fig. 2(b), we investigated the effect
of different dropout rates, ranging from 0.1 to 0.5, and identified 0.3 as the
optimal rate. Dropout rates below this optimal threshold are associated with
overfitting, where the model learns the noise in the training data too well,
reducing its ability to generalize to unseen data. Conversely, dropout rates
above 0.3 hinder effective learning, leading to an inadequate representation
of the underlying distribution of the data.

(a)

(b)

Figure 2: F1-score with respect to the (a) $\alpha$ and $\gamma$ parameters in
Focal Loss and (b) Dropout rate

We further investigated how different $\alpha$ and $\gamma$ parameters in the
Focal Loss function affected the model performance for imbalanced
classification tasks on the validation set. As demonstrated in Fig. 2(a),
varying levels of $\alpha$ and $\gamma$ show distinct impacts on the F1-score.
Results show that the optimal performance is achieved with $\alpha$ at 0.3 and
$\gamma$ at 2, indicating a balanced approach to addressing class imbalance
and focusing on hard examples. With $\gamma=0$, where the loss is equivalent
to Cross-Entropy, the importance given to the majority class ($\alpha$) is
best at 0.2, which is close to the ratio of the minority (Normal) class in our
dataset. Interestingly, increasing $\gamma$ to 5 does not consistently improve
model performance, contrary to expectations that higher $\gamma$ values would
more aggressively penalize easy examples and thus benefit model learning from
hard examples. This trend suggests a threshold beyond which increasing
$\gamma$ may not yield additional benefits and could potentially lead to
overfitting on hard examples or under-representing easy examples. Furthermore,
the performance improvement observed with increasing $\alpha$ values up to
0.3, followed by a slight decline, highlights the importance of $\alpha$ in
modulating the loss contribution of different classes. It suggests an optimal
balance at $\alpha=0.3$, where the model sufficiently addresses class
imbalance without disproportionately penalizing the majority class or
overlooking the minority class.

To evaluate the performance of the proposed framework, we use the metrics as
in Eq. 5 and 6.

$\begin{split}\mathrm{ACC}=&~{}\frac{TP+TN}{TP+TN+FP+FN}.\\\
\mathrm{SPE}=&~{}\frac{TN}{TN+FP}.\\\ \mathrm{SEN}=&~{}\frac{TP}{TP+FN},\\\
\mathrm{PRC}=&~{}\frac{TP}{TP+FP},\\\
\mathrm{F1-score}=&~{}\frac{2\times\mathrm{SEN}\times\mathrm{PRC}}{\mathrm{SEN}+\mathrm{PRC}}.\\\
\end{split}$ (5)

$\mathrm{MCC}=\frac{TP\times TN-FP\times
FN}{\sqrt{(TP+FP)(TP+FN)(TN+FP)(TN+FN)}}$

(6)

where $TP$, $TN$, $FP$, and $FN$ denote true positive, true negative, false
positive, and false negative, respectively. We also report area under receiver
operating characteristic curve (AUC).

###  4.3 Experimental Results

In this section, we present the experimental results of our proposed end-to-
end framework for glaucoma classification using 3D OCT images. We compared our
method against two approaches: a 3D Convolutional Neural Network (3D-CNN)
specifically designed for 3D OCT images and the RETFound model encoder
extended with two FC layers for individual 2D slice classification.

Table 4: Summary of glaucoma detection results averaged from a 5-fold cross-validation with 95% confidence intervals (±). Method | ACC | AUC | SEN | SPE | PRC | F1-score | MCC
---|---|---|---|---|---|---|---
3D-CNN [20] | 75.95 (± 10.43) | 91.66 (± 1.98) | 84.36 | 49.25 | 84.06 | 84.21 | 33.74
RETFound [35] | 84.77 (± 3.66) | 88.99 (± 4.02) | 88.90 | 69.71 | 91.47 | 90.17 | 56.61
Proposed (ViT-large + GRU) | 90.27 (± 1.44) | 95.24 (± 1.69) | 94.25 | 78.18 | 92.92 | 93.58 | 73.54

To evaluate the performance of our proposed framework, we first froze all the
pre-trained layers within our model’s backbone feature extractors that were
trained on OCT images. This strategy was employed to leverage the rich,
learned representations specific to the OCT image domain, thus ensuring the
extraction of precise and relevant features for glaucoma classification. We
then trained our proposed framework using a 5-fold cross-validation method.
The performance metrics averaged across all folds are meticulously presented
in Table 4, offering a comprehensive overview of the model’s effectiveness in
classifying glaucoma.

The 3D-CNN architecture used for comparative analysis is identical to that
proposed by Maetschke et al. [20]. Additionally, we extended the RETFound
model encoder by integrating two FC layers with 256 and 64 nodes,
respectively. This model was designed to process 2D slices from the 3D OCT
scans, each annotated with subject-level labels, thereby aligning with the
dimensional requirements of our study. To ensure a fair and rigorous
comparison, both of these architectures underwent the same 5-fold cross-
validation training approach as our proposed framework. This consistency
ensured that the same instances were used across each fold, thereby
facilitating an equitable and direct evaluation against our framework.

As summarized in Table 4, our proposed framework significantly outperforms the
other methods across all evaluation metrics. Specifically, our model achieved
an accuracy of 90.27%, which is substantially higher than the 75.95% accuracy
of the 3D-CNN model and the 84.77% accuracy of the standalone RETFound model.
Moreover, the 3D-CNN model and standalone RETFound model exhibit wide accuracy
confidence intervals of ±10.43% and ±3.66%, respectively, derived from the
5-fold cross-validation process. This indicates a lower level of reliability
and variable performance across different data subsets, underscoring
uncertainty in their predictive capabilities. The superior AUC of 95.24%
further emphasizes the proposed model’s outstanding ability to distinguish
between glaucomatous and non-glaucomatous scans.

Importantly, despite the challenges posed by an imbalanced dataset, the
proposed model demonstrated a remarkable capability to accurately identify
true positives and true negatives, achieving sensitivity and specificity
scores of 94.25% and 78.18%, respectively. This performance is particularly
noteworthy when compared to the 3D-CNN model, which showed a clear bias
towards the majority class, as evidenced by its significantly lower
specificity score of 49.25%. The precision, F1-score, and Matthews Correlation
Coefficient (MCC) further substantiate the robustness of the proposed model,
with scores that notably outperform those of the 3D-CNN and standalone
RETFound models.

These results underscore our model’s nuanced capability to capture both the
local and global structural integrity of the retina through the integration of
slice-based feature extraction and sequential analysis, thereby ensuring a
comprehensive understanding of the intricate patterns associated with
glaucoma.

###  4.4 Ablation Study

We evaluated the contribution of individual components of the proposed
framework in glaucoma detection. This experiment involves selectively
replacing components and assessing the resulting impact on model
effectiveness, thereby identifying the most crucial parts for achieving
optimal outcomes. Initially, we replaced the feature extraction component,
which is based on the pre-trained ViT-large model, with a pre-trained ResNet34
[9] model. While our ViT-large backbone feature extractor is a specialized
model pre-trained on OCT images, ResNet34 has been broadly pre-trained on the
ImageNet dataset, leading to an inferior performance as depicted in Table 5.
This outcome highlights the importance of task-specific model training.

Table 5: Summary of ablation study results averaged from a 5-fold cross-validation with 95% confidence intervals (±). Method | ACC | AUC | SEN | SPE | PRC | F1-score | MCC
---|---|---|---|---|---|---|---
ResNet34 + GRU | 87.39 (± 1.73) | 92.99 (± 1.75) | 92.03 | 72.86 | 91.38 | 91.71 | 65.40
ViT-large + SVM (Slice 1) | 77.47 (± 3.39) | 83.79 (± 1.25) | 91.29 | 51.68 | 77.92 | 84.08 | 48.15
ViT-large + SVM (Slice 31) | 80.09 (± 1.43) | 88.08 (± 1.17) | 92.18 | 55.71 | 80.76 | 86.09 | 53.02
ViT-large + SVM (Slice 32) | 82.44 (± 2.69) | 89.34 (± 1.46) | 92.89 | 59.71 | 83.35 | 87.87 | 57.49
ViT-large + SVM (Slice 37) | 82.52 (± 3.58) | 89.34 (± 2.23) | 92.79 | 59.94 | 83.59 | 87.95 | 57.49
ViT-large + SVM (Slice 52) | 78.11 (± 3.71) | 86.25 (± 1.98) | 92.42 | 52.51 | 77.69 | 84.41 | 50.67
ViT-large + SVM (Majority voting) | 85.05 (± 2.23) | 92.62 (± 1.33) | 94.16 | 64.31 | 85.71 | 89.74 | 63.33
Proposed (ViT-large + GRU) | 90.27 (± 1.44) | 95.24 (± 1.69) | 94.25 | 78.18 | 92.92 | 93.58 | 73.54

(a)

(b)

(c)

Figure 3: t-SNE visualization of embedding spaces for (a) ViT-large output
features for slice 32, (b) ResNet34 output features for slice 32, (c) ViT-
large + GRU output features

In the subsequent experiment, we assessed the impact of replacing the GRU
model in the proposed framework with a voting ensemble model. To address the
challenge of incorporating all available slices in 3D OCT images, we initially
reduced the number of slices. Therefore, from the original set of 64 OCT
slices, we chose 5 slices based on their entropy values. Entropy, in this
context, measures the information content for identifying the most relevant
slices for glaucoma detection. We computed the entropy of all slices in the
OCT images. Considering the significance of the central slice in glaucoma
diagnosis [33], we selected four slices with entropy values close to the
$32^{\mathrm{nd}}$ slice and the central slice itself. This method ensures a
balance between data-driven selection and clinical relevance.

After extracting features from the five slices, we proceeded with feature
selection to reduce the feature vector’s dimensionality. We used the gain
ratio method to reduce each slice’s feature dimensions from 1024 to 128
individually. Each slice, now represented by the reduced feature set, was
independently processed through a Support Vector Machine (SVM) to classify it.
Additionally, we employed a majority voting ensemble technique to integrate
the insights obtained from all five slices into a coherent diagnostic
evaluation. The results of these experiments are detailed in Table 5.

We used t-distributed Stochastic Neighbor Embedding (t-SNE) [31] to explore
the model outputs further. Specifically, for slice 32, we used t-SNE to
visualize the output feature spaces generated by pre-trained ViT-large and
ResNet34 feature extractors, as shown in Fig. 3(a) and Fig. 3(b),
respectively. By comparing the spatial distributions of features for the
normal and glaucoma classes, we observed that the normal class samples were
more densely clustered within the ViT-large feature space than in the ResNet34
feature space. This finding suggests that ViT-large has superior
discriminative capabilities and is more effective in identifying and
distinguishing the intricate patterns characteristic of glaucoma in OCT
images.

We have visualized the feature space in our proposed framework, which is the
last layer before the FC layer. This is shown in Fig. 3(c). The visualization
clearly differentiates between the normal and glaucoma classes, which is a
testament to the effectiveness of our framework. Our framework is highly
capable of extracting clinically relevant features from OCT image slices and
distinguishing sequential patterns to improve glaucoma detection accuracy.

##  5 Conclusion

In this study, we introduced a novel deep learning framework that addresses
the critical challenge of accurate and early glaucoma detection from 3D OCT
imaging. By leveraging the complementary strengths of a pre-trained Vision
Transformer and a bidirectional GRU, the proposed approach enables rich
features to be extracted from individual B-scan slices and capture inter-slice
spatial dependencies associated with the ocular structure. The holistic
analysis of local retinal characteristics and global structural integrity
captures the intricate patterns associated with glaucoma.

The experimental results validate the framework’s superior performance,
achieving an AUC of 95.24%, F1-score of 93.58%, and MCC of 73.54%. These
metrics outperform state-of-the-art methods, including the 3D-CNN [20] and
RETFound model [35] extended for slice-based classification. These results
underscore the potential of the proposed approach to address challenges posed
by class imbalance and subtle glaucomatous patterns distributed across the
ocular anatomy, enhancing clinical decision support systems and improving
patient outcomes in glaucoma management.

In our opinion, there are several opportunities for future research. Exploring
the integration of complementary clinical data modalities, such as visual
field tests or patient demographics, could further enhance the diagnostic
capabilities of the proposed framework. Moreover, investigating alternative
sequence processing models and incorporating attention mechanisms may yield
further performance improvements in glaucoma detection and screening tasks.
Additionally, extending the application scope of our framework coupled with
interpretability to other ocular pathologies could significantly contribute to
the advancement of early diagnosis and treatment strategies within the broader
field of ophthalmology. Furthermore, developing robust techniques to mitigate
the impact of data imbalance could substantially enhance the diagnostic
accuracy and generalizability of the framework across diverse datasets and
imaging modalities. Such strategies would enable a comprehensive assessment of
the framework’s robustness and facilitate its broader deployment in clinical
ophthalmological practice.

## References

## References

  * [1] Mona Ashtari-Majlan, Mohammad Mahdi Dehshibi, and David Masip.  Deep learning and computer vision for glaucoma detection: A review.  arXiv preprint arXiv:2307.16528, 2023.
  * [2] Peter Åsman and Anders Heijl.  Glaucoma Hemifield Test: Automated Visual Field Evaluation.  Archives of Ophthalmology, 110(6):812–819, 1992.
  * [3] Zailiang Chen, Xianxian Zheng, Hailan Shen, Ziyang Zeng, Qing Liu, and Zhuo Li.  Combination of enhanced depth imaging optical coherence tomography and fundus images for glaucoma screening.  Journal of Medical Systems, 43:1–12, 2019.
  * [4] Alexey Dosovitskiy, Lucas Beyer, Alexander Kolesnikov, Dirk Weissenborn, Xiaohua Zhai, Thomas Unterthiner, Mostafa Dehghani, Matthias Minderer, Georg Heigold, Sylvain Gelly, Jakob Uszkoreit, and Neil Houlsby.  An image is worth 16x16 words: Transformers for image recognition at scale.  In International Conference on Learning Representations, 2021.
  * [5] Gabriel García, Rocío del Amor, Adrián Colomer, and Valery Naranjo.  Glaucoma detection from raw circumpapillary oct images using fully convolutional neural networks.  In 2020 IEEE International Conference on Image Processing (ICIP), pages 2526–2530, 2020.
  * [6] Gabriel García, Rocío Del Amor, Adrián Colomer, Rafael Verdú-Monedero, Juan Morales-Sánchez, and Valery Naranjo.  Circumpapillary oct-focused hybrid learning for glaucoma grading using tailored prototypical neural networks.  Artificial Intelligence in Medicine, 118:102132, 2021.
  * [7] Gabriel García, Adrián Colomer, and Valery Naranjo.  Glaucoma detection from raw sd-oct volumes: A novel approach focused on spatial dependencies.  Computer Methods and Programs in Biomedicine, 200:105855, 2021.
  * [8] Yasmeen George, Bhavna J. Antony, Hiroshi Ishikawa, Gadi Wollstein, Joel S. Schuman, and Rahil Garnavi.  Attention-guided 3d-cnn framework for glaucoma detection and structural-functional association using volumetric images.  IEEE Journal of Biomedical and Health Informatics, 24(12):3421–3430, 2020.
  * [9] Kaiming He, Xiangyu Zhang, Shaoqing Ren, and Jian Sun.  Deep residual learning for image recognition.  In 2016 IEEE Conference on Computer Vision and Pattern Recognition (CVPR), pages 770–778, 2016.
  * [10] Sheng He, Randy L. Gollub, Shawn N. Murphy, Juan David Perez, Sanjay Prabhu, Rudolph Pienaar, Richard L. Robertson, P. Ellen Grant, and Yangming Ou.  Brain age estimation using lstm on children’s brain mri.  In 2020 IEEE 17th International Symposium on Biomedical Imaging (ISBI), pages 1–4, 2020.
  * [11] Ruben Hemelings, Bart Elen, João Barbosa-Breda, Matthew B Blaschko, Patrick De Boever, and Ingeborg Stalmans.  Deep learning on fundus images detects glaucoma beyond the optic disc.  Scientific Reports, 11(1):20313, 2021.
  * [12] David Huang, Eric A Swanson, Charles P Lin, Joel S Schuman, William G Stinson, Warren Chang, Michael R Hee, Thomas Flotte, Kenton Gregory, Carmen A Puliafito, and James G. Fujimoto.  Optical coherence tomography.  Science, 254(5035):1178–1181, 1991.
  * [13] Jost B. Jonas, Tin Aung, Rupert R. Bourne, Alain M. Bron, Robert Ritch, and Songhomitra Panda-Jonas.  Glaucoma.  The Lancet, 390(10108):2183–2193, 2017.
  * [14] Mamta Juneja, Janmejai Singh Minhas, Naveen Singla, Sarthak Thakur, Niharika Thakur, and Prashant Jindal.  Fused framework for glaucoma diagnosis using optical coherence tomography (oct) images.  Expert Systems with Applications, 201:117202, 2022.
  * [15] Diederik P Kingma and Jimmy Ba.  Adam: A method for stochastic optimization.  In Third International Conference on Learning Representations, ICLR, pages 1––15, 2015.
  * [16] Thomas Kurmann, Pablo Márquez-Neila, Siqing Yu, Marion Munk, Sebastian Wolf, and Raphael Sznitman.  Fused detection of retinal biomarkers in oct volumes.  In Dinggang Shen, Tianming Liu, Terry M. Peters, Lawrence H. Staib, Caroline Essert, Sean Zhou, Pew-Thian Yap, and Ali Khan, editors, Medical Image Computing and Computer Assisted Intervention – MICCAI 2019, pages 255–263. Springer International Publishing, 2019.
  * [17] Jinho Lee, Young Kook Kim, Ki Ho Park, and Jin Wook Jeoung.  Diagnosing glaucoma with spectral-domain optical coherence tomography using deep learning classifier.  Journal of glaucoma, 29(4):287–294, 2020.
  * [18] Tsung-Yi Lin, Priya Goyal, Ross Girshick, Kaiming He, and Piotr Dollar.  Focal loss for dense object detection.  In Proceedings of the IEEE International Conference on Computer Vision (ICCV), Oct 2017.
  * [19] Manhua Liu, Danni Cheng, Weiwu Yan, and Alzheimer’s Disease Neuroimaging Initiative .  Classification of alzheimer’s disease by combination of convolutional and recurrent neural networks using fdg-pet images.  Frontiers in Neuroinformatics, 12, 2018.
  * [20] Stefan Maetschke, Bhavna Antony, Hiroshi Ishikawa, Gadi Wollstein, Joel Schuman, and Rahil Garnavi.  A feature agnostic approach for glaucoma detection in oct volumes.  PLOS ONE, 14(7):1–12, 2019.
  * [21] Felipe A. Medeiros, Alessandro A. Jammal, and Atalie C. Thompson.  From machine to machine: An oct-trained deep learning algorithm for objective quantification of glaucomatous damage in fundus photographs.  Ophthalmology, 126(4):513–521, 2019.
  * [22] Kouros Nouri-Mahdavi, Vahid Mohammadzadeh, Alessandro Rabiolo, Kiumars Edalati, Joseph Caprioli, and Siamak Yousefi.  Prediction of visual field progression from oct structural measures in moderate to advanced glaucoma.  American Journal of Ophthalmology, 226:172–181, 2021.
  * [23] Erfan Noury, Suria S. Mannil, Robert T. Chang, An Ran Ran, Carol Y. Cheung, Suman S. Thapa, Harsha L. Rao, Srilakshmi Dasari, Mohammed Riyazuddin, Dolly Chang, Sriharsha Nagaraj, Clement C. Tham, and Reza Zadeh.  Deep Learning for Glaucoma Detection and Identification of Novel Diagnostic Areas in Diverse Real-World Datasets.  Translational Vision Science & Technology, 11(5):11–11, 2022.
  * [24] Hina Raja, M. Usman Akram, Sajid Gul Khawaja, Muhammad Arslan, Aneeqa Ramzan, and Noman Nazir.  Data on oct and fundus images for the detection of glaucoma.  Data in Brief, 29:105342, 2020.
  * [25] Hina Raja, Taimur Hassan, Muhammad Usman Akram, and Naoufel Werghi.  Clinically verified hybrid deep learning system for retinal ganglion cells aware grading of glaucomatous progression.  IEEE Transactions on Biomedical Engineering, 68(7):2140–2151, 2021.
  * [26] An Ran Ran, Carol Y Cheung, Xi Wang, Hao Chen, Lu yang Luo, Poemen P Chan, Mandy O M Wong, Robert T Chang, Suria S Mannil, Alvin L Young, Hon wah Yung, Chi Pui Pang, Pheng-Ann Heng, and Clement C Tham.  Detection of glaucomatous optic neuropathy with spectral-domain optical coherence tomography: a retrospective training and validation deep-learning analysis.  The Lancet Digital Health, 1(4):e172–e182, 2019.
  * [27] Ashraf Haroon Rashid, Aditya Gupta, Jhalak Gupta, and Mohammad Tanveer.  Biceph-net: A robust and lightweight framework for the diagnosis of alzheimer’s disease using 2d-mri scans and deep similarity learning.  IEEE Journal of Biomedical and Health Informatics, 27(3):1205–1213, 2023.
  * [28] Daniel B. Russakoff, Suria S. Mannil, Jonathan D. Oakley, An Ran Ran, Carol Y. Cheung, Srilakshmi Dasari, Mohammed Riyazzuddin, Sriharsha Nagaraj, Harsha L. Rao, Dolly Chang, and Robert T. Chang.  A 3D Deep Learning System for Detecting Referable Glaucoma Using Full OCT Macular Cube Scans.  Translational Vision Science & Technology, 9(2):12–12, 2020.
  * [29] Diping Song, Bin Fu, Fei Li, Jian Xiong, Junjun He, Xiulan Zhang, and Yu Qiao.  Deep Relation Transformer for Diagnosing Glaucoma With Optical Coherence Tomography and Visual Field Function.  IEEE Transactions on Medical Imaging, 40(9):2392–2402, 2021.
  * [30] Jaimie D. Steinmetz et al.  Causes of blindness and vision impairment in 2020 and trends over 30 years, and prevalence of avoidable blindness in relation to VISION 2020: the Right to Sight: an analysis for the Global Burden of Disease Study.  The Lancet Global Health, 9(2):e144–e160, 2021.
  * [31] Laurens van der Maaten and Geoffrey Hinton.  Visualizing data using t-sne.  Journal of Machine Learning Research, 9(86):2579–2605, 2008.
  * [32] Ramdas Vankdothu, Mohd Abdul Hameed, and Husnah Fatima.  A brain tumor identification and classification using deep learning based on cnn-lstm method.  Computers and Electrical Engineering, 101:107960, 2022.
  * [33] Xi Wang, Hao Chen, An-Ran Ran, Luyang Luo, Poemen P. Chan, Clement C. Tham, Robert T. Chang, Suria S. Mannil, Carol Y. Cheung, and Pheng-Ann Heng.  Towards multi-center glaucoma oct image screening with semi-supervised joint structure and function multi-task learning.  Medical Image Analysis, 63:101695, 2020.
  * [34] Robert N. Weinreb, Tin Aung, and Felipe A. Medeiros.  The Pathophysiology and Treatment of Glaucoma: A Review.  JAMA, 311(18):1901–1911, 05 2014.
  * [35] Yukun Zhou et al.  A foundation model for generalizable disease detection from retinal images.  Nature, 622(7981):156–163, 2023.


HIG-21-008

HIG-21-008

[cern]CMS Collaboration

# Search for Higgs boson decay to a charm quark-antiquark pair in proton-
proton collisions at $\sqrt{s}=13\TeV$

###### Abstract

A search for the standard model Higgs boson decaying to a charm quark-
antiquark pair, ${\PH\to\ccbar}$, produced in association with a leptonically
decaying (or ) boson is presented. The search is performed with proton-proton
collisions at ${\sqrt{s}=13\TeV}$ collected by the CMS experiment,
corresponding to an integrated luminosity of 138. Novel charm jet
identification and analysis methods using machine learning techniques are
employed. The analysis is validated by searching for ${\PZ\to\ccbar}$ in
$\PV\PZ$ events, leading to its first observation at a hadron collider with a
significance of 5.7 standard deviations. The observed (expected) upper limit
on
$\sigma_{\text{}}\left(\PV\PH\right)\mathcal{B}_{\text{}}\left(\PH\to\PQc\PAQc\right)$
is 0.94 ($0.50^{+0.22}_{-0.15}$)pbat $95\%$ confidence level (), corresponding
to 14 ($7.6^{+3.4}_{-2.3}$) times the standard model prediction. For the
Higgs-charm Yukawa coupling modifier, $\kappa_{\PQc}$, the observed (expected)
$95\%$ interval is ${1.1<\abs{\kappa_{\PQc}}<5.5}$
(${\abs{\kappa_{\PQc}}<3.4}$), the most stringent constraint to date.

The discovery of a Higgs boson () with the LHC Run 1 data by the ATLAS [1] and
CMS [2, 3] experiments in 2012 was a major advance in the understanding of the
electroweak (EW) symmetry breaking mechanism. The measured Higgs boson mass is
$125.38\pm 0.14$ [4]. The observed interactions with gauge bosons and third-
generation fermions [5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15], and all measured
properties [16, 17, 18, 19, 20, 21, 22, 23, 24, 4], are compatible with
standard model (SM) predictions. Recently, the CMS Collaboration reported the
first evidence of Higgs boson decays to muons, , second-generation leptons
[25]. An important next milestone is the observation of its coupling to
second-generation quarks. In this Letter, we focus on a search for Higgs boson
decay to $\ccbar$, a charm quark-antiquark pair. The corresponding Yukawa
coupling, $y_{\PQc}$, can be significantly modified in the presence of physics
beyond the SM [26, 27, 28, 29]. However, the small branching ratio predicted
by the SM, ubiquitous production of quark and gluon jets at the LHC, and the
difficulty of identifying charm quark jets in a hadronic environment,
including distinguishing them from bottom quark jets, make this a challenging
measurement. Searches for ${\PH\to\ccbar}$ reported in Refs. [30, 31, 32] by
the ATLAS and CMS Collaborations target the associated production of a Higgs
boson with a (or ) boson. Using 139 of data at 13, the most recent search by
the ATLAS Collaboration obtains an observed (expected) upper limit on the
product of the production cross section $\sigma_{\text{}}\left(\PV\PH\right)$
and branching fraction $\mathcal{B}_{\text{}}\left(\PH\to\PQc\PAQc\right)$ of
26 (31) times the SM prediction at 95% confidence level () [32].

This Letter presents a search for ${\PH\to\ccbar}$ in $\PV\PH$ production
using proton-proton ($\Pp\Pp$) collision data at $\sqrt{s}=13\TeV$, collected
by the CMS detector in 2016–2018, and corresponding to an integrated
luminosity of 138 [33, 34, 35]. Building upon a previous analysis by the CMS
Collaboration [31], advanced jet reconstruction and identification algorithms
and sophisticated analysis techniques using machine learning (ML) have been
further developed to significantly improve sensitivity.

The CMS apparatus [36] is a multipurpose, nearly hermetic detector, designed
to trigger on [37, 38] and identify electrons, muons, photons, and (charged
and neutral) hadrons [39, 40, 41]. A global reconstruction “particle-flow”
(PF) algorithm [42] combines the information provided by the all-silicon inner
tracker and by the crystal electromagnetic and brass-scintillator hadron
calorimeters, operating inside a 3.8T superconducting solenoid, with data from
gas-ionization muon detectors embedded in the solenoid flux-return yoke, to
reconstruct charged and neutral hadrons, electrons, photons, and muons (PF
candidates) and build leptons, jets, missing transverse momentum, and other
physics objects [43, 44, 45].

Two collections of jets, formed by PF candidates clustered using the anti-
algorithm [46, 47], are used in the search. The first uses a distance
parameter $R=0.4$, and will be referred to as “small-$R$” jets. The second
uses $R=1.5$ and contains what are referred to as “large-$R$” jets. The impact
of particles from additional $\Pp\Pp$ interactions within the same or nearby
bunch crossings (pileup) is mitigated via the charged hadron subtraction
algorithm [42] for small-$R$ jets, and via the PUPPI [48, 49] algorithm for
large-$R$ jets. A regression algorithm [50] is developed to improve large-$R$
jet mass reconstruction that exploits properties of the PF candidates and
secondary vertices associated to the jet using the ParticleNet graph neural
network [51]. Mass resolution is improved by about 50% over traditional jet
grooming algorithms [52, 53]. The small-$R$ (large-$R$) jets are required to
have transverse momentum () above 25 (200)and to be within the tracker
acceptance.

Signal and background processes are simulated using various Monte Carlo event
generators. The detector response is modeled with [54]. The $\PW\PH$ and
quark-induced $\PZ\PH$ signal processes are generated at next-to-leading order
(NLO) accuracy in quantum chromodynamics (QCD) using the v2 [55, 56, 57] event
generator extended with the multi-scale improved NLO (MiNLO) procedure [58,
59], while the gluon-induced $\PZ\PH$ process is generated at leading order
(LO) accuracy with v2. The Higgs boson mass is set to 125for all simulations.
The production cross sections of the signal processes [60] are corrected as a
function of $\pt(\PV)$ to next-to-next-to-leading order (NNLO) QCD + NLO EW
accuracy combining the vhnnlo [61, 62, 63, 64], vh@nnlo [65, 66], and hawk
v2.0 [67] generators, as described in Ref. [60].

The +jets background samples are generated with v2.6.0 [68] at NLO with up to
two additional partons. The top quark pair () [69] and single top quark
production processes [70, 71, 72] are generated to NLO accuracy with v2. The
production cross sections for the samples are scaled to the NNLO prediction
with the next-to-next-to-leading logarithmic resummation result obtained from
v2.0 [73], and the differential cross sections as a function of top quark are
corrected to the NNLO QCD + NLO EW prediction [74]. Diboson backgrounds are
generated at NLO with v2 ( v2.4.2) for the $\PW\PW$ [75] ($\PW\PZ$ and
$\PZ\PZ$) process. Production cross sections of the diboson processes are
reweighted as a function of the subleading vector boson to NNLO QCD + NLO EW
accuracy [76].

The NLO NNPDF3.0 [77] (NNLO NNPDF3.1 [78]) parton distribution function (PDF)
set is used for the 2016 (2017 and 2018) simulations. For parton showering and
hadronization, including the ${\PH\to\ccbar}$ decay, matrix element generators
are interfaced with v8.230 [79] with the CUETP8M1 [80] (CP5 [81]) underlying
event tune for 2016 (2017 and 2018) samples. The matching of jets from matrix
element calculations and those from parton showers is done with the FxFx [82]
(MLM [83]) prescription for NLO (LO) samples. For all samples, pileup
interactions are simulated with and added to the hard-scattering process.
Events are then reweighted to match the pileup profile observed in data.

The analysis is carried out in mutually exclusive channels targeting leptonic
decays of the vector bosons: $\PZ\to\nu\nu$, $\PW\to\ell\nu$, and
$\PZ\to\ell\ell$, where $\ell$ is an electron or a muon, and referred to as
the 0L, 1L, and 2L channels. Events are collected using triggers based on
large missing transverse momentum, or the presence of one, or two well-
identified and isolated leptons. The event selection criteria are detailed in
Ref. [31].

As for the previous search [31], the analysis reconstructs the Higgs boson
candidate ($\PH_{\text{cand}}$) assuming either a “merged-jet” topology, in
which the hadronization products of the two charm quarks are reconstructed as
a single large-$R$ jet, or a “resolved-jet” topology, in which the
$\PH_{\text{cand}}$ is reconstructed from two well-separated and individually
resolved small-$R$ jets. These two topologies can have significant overlap and
are made distinct in reference to whether a given $\PH_{\text{cand}}$,
identified through a large-$R$ jet in the event, has above or below a
threshold of 300for the definition of the merged-jet and the resolved-jet
topology, respectively, where the dividing line is chosen to maximize
sensitivity to ${\PV\PH(\PH\to\ccbar)}$.

On average, the $\PH_{\text{cand}}$ large-$R$ jets in the signal process have
larger than those from the +jets and backgrounds. Thus, the high-regime
explored in the merged-jet topology, although amounting to less than 5% of the
signal cross section, provides significant sensitivity to the search. Because
the highly boosted Higgs boson decay produces a more narrowly concentrated
collection of final-state particles, the use of a single large-$R$ jet
enhances the probability of correctly capturing the hadronization products of
both charm quarks, and facilitates the accounting of final-state radiation
(FSR) emitted by the quarks. A detailed discussion of the advantages can be
found in Refs. [53, 84, 31].

Figure 1:  Performance of ParticleNet (blue lines) for identifying a $\ccbar$
pair for large-$R$ jets with $\pt>300$. The solid (dashed) line shows the
efficiency to correctly identify ${\PH\to\ccbar}$ vs. the efficiency of
misidentifying quarks or gluons from the +jets process (vs. ${\PH\to\bbbar}$).
The red crosses represent the three working points used in the merged-jet
analysis. The performance of DeepAK15 (yellow lines) used in Ref. [31] is
shown for comparison.

State-of-the-art performance in the reconstruction and identification of the
pair of quarks from the Higgs boson decay is achieved with ParticleNet [51], a
novel jet identification algorithm. Using PF candidates and secondary vertices
associated to large-$R$ jets as inputs, ParticleNet simultaneously exploits
information related to jet substructure, flavor, and pileup with a graph
neural network [85], yielding substantial gains over other approaches [86,
87]. Decorrelation of the algorithm’s response with the jet mass is achieved
by training it with a dedicated set of simulations produced with the same jet
mass distributions for the signal and background processes [86]. Figure 1
shows the performance of the $\ccbar$ discriminant in identifying a pair of
quarks from Higgs boson decay for large-$R$ jets with $\pt>300$. ParticleNet
is compared to the previous state-of-the-art $\ccbar$ discriminant “DeepAK15”
[31, 88], yielding an improvement by a factor of 4 to 7 in the rejection of
other jet flavors. Three working points are defined on the $\ccbar$
discriminant distribution with approximately 58, 40, and 16% efficiencies for
identifying a $\ccbar$ pair. The corresponding misidentification rates of
light quark and gluon jets ($\bbbar$ jets) are 2 (9), 0.7 (5), and 0.08 (1)%.
These working points are used to separate events into three mutually exclusive
categories with different $\ccbar$ purity to improve the sensitivity of the
analysis. The $\ccbar$ identification efficiency in data is measured using a
sample of events containing a gluon splitting to $\ccbar$. To increase the
similarity to ${\PH\to\ccbar}$ decay, a dedicated BDT classifier is developed
to enrich jets where a large fraction of momentum is carried by the quark pair
from gluon splitting rather than by additionally radiated gluons [89]. The
-dependent data-to-simulation efficiency ratios (used as corrective scale
factors) are typically 0.9–1.3 with corresponding uncertainties of 20–30%.

The main backgrounds, and +jets, are suppressed by a separate boosted decision
tree (BDT) classifier for each channel, using kinematical variables that are
not correlated with the $\PH_{\text{cand}}$ mass,
$m\left(\PH_{\text{cand}}\right)$, or the $\ccbar$ discriminant as inputs. The
BDT design relies on previous developments [31] with improvements in variable
selection and training procedure, leading to ${\approx}$15% enhancement of the
sensitivity of the analysis. The BDT discriminants are used to define 2 (1)
signal regions (SRs) in the 1L (0L and 2L) channel. Events in the SRs are
further subdivided into the three $\ccbar$ discriminant categories mentioned
above. The $m\left(\PH_{\text{cand}}\right)$ distributions are used to
separate signal and background contributions in each SR, as both the BDTs and
the $\ccbar$ discriminant are designed to be largely independent of
$m\left(\PH_{\text{cand}}\right)$.

More than 95% of the $\PV\PH$ events have a Higgs boson with
$\pt(\PH)<300\GeV$, corresponding to the phase space region where the Higgs
boson decay products generally give rise to two distinctly reconstructed
small-$R$ jets. The resolved-jet topology exploits a large fraction of this
phase space, which, however, contains higher background contamination than
that used in the merged-jet analysis. The $\PH_{\text{cand}}$ is reconstructed
via two distinct small-$R$ jets. The identification of jets relies on the ML-
based DeepJet algorithm [90, 91]. The discrimination between jets and light
quark or gluon jets (jets) is achieved via the ratio of the corresponding
probabilities, defined as CvsL (CvsB) [92]. Potential differences in the
discriminant shapes between data and simulation are taken into account by
flavor-dependent simulation-to-data efficiency scale factors [92], which are
typically in the 0.9–1.0 range. The corresponding uncertainties range from
${\approx}$2% for light quark and gluon jets or jets, to ${\approx}$5% for
jets.

The two jets with the highest CvsL discriminants in each event are selected to
reconstruct the $\PH_{\text{cand}}$ momentum four-vector. To improve sample
purity, criteria are imposed on CvsL and CvsB for the leading jet that
correspond to ${\approx}40\%$ jet efficiency, for ${\approx}4$
(${\approx}16$)% light quark or gluon jet (jet) misidentification rate.
Multiple corrections are applied to improve the
$m\left(\PH_{\text{cand}}\right)$ reconstruction. To account for a potential
underestimation of the jet energy, due to the presence of undetected neutrinos
in hadron decays, an ML-based jet energy regression algorithm [93], originally
developed for jets, is adapted for jets. In addition, small-$R$ jets
reconstructed in the vicinity of $\PH_{\text{cand}}$ jets, often stemming from
FSR, are included in the $\PH_{\text{cand}}$ reconstruction [31]. The
$m\left(\PH_{\text{cand}}\right)$ resolution in the 2L channel is further
improved via a kinematic fit by balancing the momenta of the two small-$R$
jets and the lepton pair within experimental uncertainties [94]. These steps
improve the $m\left(\PH_{\text{cand}}\right)$ resolution up to 20%. Finally, a
BDT classifier is developed to maximize the discrimination power between
signal and background processes in each channel, using event-level kinematical
variables, jet identification discriminants, and properties of
$\PH_{\text{cand}}$, including $m\left(\PH_{\text{cand}}\right)$, as inputs
[31].

The signal strength modifier $\mu$, defined as
${\left(\sigma\mathcal{B}\right)_{\text{obs}}/\left(\sigma\mathcal{B}\right)_{\text{SM}}}$
where $\sigma$ is the signal production cross section and $\mathcal{B}$ is the
branching fraction, is measured via a binned maximum likelihood fit to data.
The best-fit value of $\mu$ and an approximate 68% confidence interval are
extracted following the procedure in Ref. [17]. The fitted variable is
$m\left(\PH_{\text{cand}}\right)$ in the merged-jet analysis, and the BDT
discriminant in the resolved-jet analysis. The normalizations of the main
backgrounds, namely +jets and , are estimated by including dedicated control
regions in the fit following the strategy detailed in Ref. [31]. Contributions
from single top, diboson and ${\PV\PH(\PH\to\bbbar)}$ processes are estimated
from simulation assuming SM production rates. Because of improvements in jet
identification, the difficult ${\PV\PH(\PH\to\bbbar)}$ background contribution
in the high $\ccbar$ purity SRs is reduced from about twenty times to only
about twice that expected for ${\PV\PH(\PH\to\ccbar)}$, and the
${\PV\PZ(\PZ\to\bbbar)}$ yield is reduced from more than 100% to about 10% of
that expected for ${\PV\PZ(\PZ\to\ccbar)}$. Contributions from
${\PH\to\PGt\PGt}$ decays are negligible after the above-mentioned selection
criteria, and are not considered.

Systematic uncertainties affecting normalizations and shapes of fitted
variables are taken into account via nuisance parameters. The relative
contribution of each uncertainty source to the total uncertainty in the fitted
$\mu$ is summarized in Table Search for Higgs boson decay to a charm quark-
antiquark pair in proton-proton collisions at $\sqrt{s}=13\TeV$. The leading
uncertainty is statistical because of the limited number of events in the SRs
as well as the control regions used to extract background normalizations. The
main experimental systematic uncertainties are associated with limited
simulation sample sizes, particularly due to large fractions of negatively
weighted events in the NLO +jets samples, and the jet identification
efficiencies, representing ${\approx}37$ and ${\approx}23\%$ of the total,
respectively. Theoretical uncertainties in the cross sections, spectra, PDFs,
renormalization and factorization scales, represent ${\approx}22\%$ of the
total uncertainty in $\mu$.

The relative contributions to the total uncertainty in the signal strength
modifier $\mu$ for the ${\PV\PH(\PH\to\ccbar)}$ process, where the best fit is
${\mu_{{\PV\PH(\PH\to\ccbar)}}=7.7^{+3.8}_{-3.5}}$. lr Uncertainty source &
$\Delta\mu/\left(\Delta\mu\right)_{\text{tot}}$

Statistical 85%
Background normalizations 37%

Experimental 48%
Sizes of the simulated samples 37%
jet identification efficiencies 23%
Jet energy scale and resolution 15%
Simulation modeling 11%
Integrated luminosity 6%
Lepton identification efficiencies 4%

Theory 22%
Backgrounds 17%
Signal 15%

A search for the analogous SM process ${\PV\PZ(\PZ\to\ccbar)}$ is performed to
validate the analysis strategy. The BDTs in the resolved-jet topology are
modified by training them with ${\PV\PZ(\PZ\to\ccbar)}$ as signal. No
modification is needed for the BDTs in the merged-jet topology as they are
independent of $m\left(\PH_{\text{cand}}\right)$. The best fit $\mu$ of this
process is ${\mu_{{\PV\PZ(\PZ\to\ccbar)}}=1.01_{-0.21}^{+0.23}}$, in agreement
with the SM expectation. Figure 2 (left) shows the distribution of events in
all channels, sorted into bins of similar signal-to-background ratios. The
observed data shows a visible excess over the expected backgrounds. The
significance of the excess is computed using the asymptotic distribution of a
test statistic based on the profile likelihood ratio [95, 96]. The observed
(expected) significance is 4.4 (4.7) standard deviations for the merged-jet
analysis, 3.1 (3.3) standard deviations for the resolved-jet analysis, and 5.7
(5.9) standard deviations for their combination. This is the first observation
of ${\PZ\to\ccbar}$ at a hadron collider.

Figure 2: Distribution of events as a function of $S/B$ in the
${\PV\PZ(\PZ\to\ccbar)}$ (left) and ${\PV\PH(\PH\to\ccbar)}$ (right) searches,
where $S$ and $B$ are the postfit signal and background yields, respectively,
in each bin of the fitted $m\left(\PH_{\text{cand}}\right)$ or BDT
discriminant distributions. The bottom panel shows the ratio of data to the
total background, with the uncertainty in background indicated by gray
hatching. The red line represents background plus SM signal divided by
background.

Figure 3: Combined $m\left(\PH_{\text{cand}}\right)$ distribution in all
channels of the merged-jet analysis. The fitted
$m\left(\PH_{\text{cand}}\right)$ distribution in each SR is weighted by
$S/(S+B)$, where $S$ and $B$ are the postfit ${\PV\PH(\PH\to\ccbar)}$ signal
and total background yields. The lower panel shows data (points) and the
fitted ${\PV\PH(\PH\to\ccbar)}$ (red) and ${\PV\PZ(\PZ\to\ccbar)}$ (grey)
distributions after subtracting all other processes. Error bars represent pre-
subtraction statistical uncertainties in data, while the gray hatching
indicates the total uncertainty in the signal and all background processes.

Figure 2 (right) compares the observed data to the SM prediction in the search
of ${\PV\PH(\PH\to\ccbar)}$, where the best fit is
${\mu_{{\PV\PH(\PH\to\ccbar)}}=7.7^{+3.8}_{-3.5}}$. The fitted
$m\left(\PH_{\text{cand}}\right)$ distribution in the merged-jet topology is
displayed in Fig. 3. No significant excess over the background-only hypothesis
is observed. An upper limit on $\mu_{{\PV\PH(\PH\to\ccbar)}}$ is extracted
using the criterion [97, 98]. The test statistic is the profile likelihood
ratio modified for upper limits [95], and the asymptotic approximation [96] is
used in the limit setting procedure. The observed (expected) 95% upper limit
on $\mu_{{\PV\PH(\PH\to\ccbar)}}$ is 14 ($7.6^{+3.4}_{-2.3}$), which is
equivalent to an observed (expected) upper limit on
$\sigma_{\text{}}\left(\PV\PH\right)\mathcal{B}_{\text{}}\left(\PH\to\PQc\PAQc\right)$
of 0.94 ($0.50^{+0.22}_{-0.15}$)pb. Contributions from the individual channels
are summarized in Fig. 4. Tabulated results are provided in the HEPData record
for this analysis [99].

The result is interpreted in the $\kappa$-framework [100, 60] by
reparameterizing $\mu_{{\PV\PH(\PH\to\ccbar)}}$ in terms of the Higgs-charm
Yukawa coupling modifier $\kappa_{\PQc}$, assuming only the Higgs boson decay
widths are altered:

$\mu_{{\PV\PH(\PH\to\ccbar)}}=\frac{\kappa_{\PQc}^{2}}{1+\mathcal{B}_{\text{SM}}\left({\PH\to\ccbar}\right)\,(\kappa_{\PQc}^{2}-1)}.$
(1)

The observed 95% interval is ${1.1<\abs{\kappa_{\PQc}}<5.5}$, and the
corresponding expected constraint is ${\abs{\kappa_{\PQc}}<3.4}$.

Figure 4: The 95% upper limits on $\mu_{{\PV\PH(\PH\to\ccbar)}}$. Green and
yellow bands indicate the 68 and 95% intervals on the expected limits,
respectively. The vertical red line indicates the SM value
$\mu_{{\PV\PH(\PH\to\ccbar)}}=1$.

In summary, a search for the SM Higgs boson decaying to a pair of charm quarks
in the CMS experiment is presented. Novel jet reconstruction and
identification tools, and analysis techniques are developed for this analysis,
which is validated by measuring the ${\PV\PZ(\PZ\to\ccbar)}$ process. The
observed boson signal relative to the SM prediction is
${\mu_{{\PV\PZ(\PZ\to\ccbar)}}=1.01_{-0.21}^{+0.23}}$, with an observed
(expected) significance of 5.7 (5.9) standard deviations above the background-
only hypothesis. This is the first observation of ${\PZ\to\ccbar}$ at a
hadronic collider.

The observed (expected) upper limit on
$\sigma_{\text{}}\left(\PV\PH\right)\mathcal{B}_{\text{}}\left(\PH\to\PQc\PAQc\right)$
is 0.94 ($0.50^{+0.22}_{-0.15}$)pb, corresponding to 14 ($7.6^{+3.4}_{-2.3}$)
times the theoretical prediction for an SM Higgs boson mass of 125.38. The
observed (expected) 95% interval on the modifier, $\kappa_{\PQc}$, for the
Yukawa coupling of the Higgs boson to the charm quark is
${1.1<\abs{\kappa_{\PQc}}<5.5}$ (${\abs{\kappa_{\PQc}}<3.4}$). This is the
most stringent constraint on $\kappa_{\PQc}$ to date.

## References

  * [1] ATLAS Collaboration, “Observation of a new particle in the search for the standard model Higgs boson with the ATLAS detector at the LHC”,  Phys. Lett. B 716 (2012) 1, 10.1016/j.physletb.2012.08.020, arXiv:1207.7214.
  * [2] CMS Collaboration, “Observation of a new boson at a mass of 125 GeV with the CMS experiment at the LHC”,  Phys. Lett. B 716 (2012) 30, 10.1016/j.physletb.2012.08.021, arXiv:1207.7235.
  * [3] CMS Collaboration, “Observation of a new boson with mass near 125 GeV in pp collisions at $\sqrt{s}$ = 7 and 8 TeV”,  JHEP 06 (2013) 081, 10.1007/JHEP06(2013)081, arXiv:1303.4571.
  * [4] CMS Collaboration, “A measurement of the Higgs boson mass in the diphoton decay channel”,  Phys. Lett. B 805 (2020) 135425, 10.1016/j.physletb.2020.135425, arXiv:2002.06398.
  * [5] ATLAS Collaboration, “Measurement of Higgs boson production in the diphoton decay channel in pp collisions at center-of-mass energies of 7 and 8 TeV with the ATLAS detector”,  Phys. Rev. D 90 (2014) 112015, 10.1103/PhysRevD.90.112015, arXiv:1408.7084.
  * [6] CMS Collaboration, “Observation of the diphoton decay of the Higgs boson and measurement of its properties”,  Eur. Phys. J. C 74 (2014) 3076, 10.1140/epjc/s10052-014-3076-z, arXiv:1407.0558.
  * [7] ATLAS Collaboration, “Measurements of Higgs boson production and couplings in the four-lepton channel in pp collisions at center-of-mass energies of 7 and 8 TeV with the ATLAS detector”,  Phys. Rev. D 91 (2015) 012006, 10.1103/PhysRevD.91.012006, arXiv:1408.5191.
  * [8] CMS Collaboration, “Measurement of the properties of a Higgs boson in the four-lepton final state”,  Phys. Rev. D 89 (2014) 092007, 10.1103/PhysRevD.89.092007, arXiv:1312.5353.
  * [9] ATLAS Collaboration, “Observation and measurement of Higgs boson decays to WW∗ with the ATLAS detector”,  Phys. Rev. D 92 (2015) 012006, 10.1103/PhysRevD.92.012006, arXiv:1412.2641.
  * [10] ATLAS Collaboration, “Study of (W/Z)H production and Higgs boson couplings using H $\rightarrow$ WW∗ decays with the ATLAS detector”,  JHEP 08 (2015) 137, 10.1007/JHEP08(2015)137, arXiv:1506.06641.
  * [11] CMS Collaboration, “Measurement of Higgs boson production and properties in the WW decay channel with leptonic final states”,  JHEP 01 (2014) 096, 10.1007/JHEP01(2014)096, arXiv:1312.1129.
  * [12] ATLAS Collaboration, “Evidence for the Higgs-boson Yukawa coupling to tau leptons with the ATLAS detector”,  JHEP 04 (2015) 117, 10.1007/JHEP04(2015)117, arXiv:1501.04943.
  * [13] CMS Collaboration, “Evidence for the 125 GeV Higgs boson decaying to a pair of $\tau$ leptons”,  JHEP 05 (2014) 104, 10.1007/JHEP05(2014)104, arXiv:1401.5041.
  * [14] CMS Collaboration, “Observation of the Higgs boson decay to a pair of tau leptons”,  Phys. Lett. B 779 (2017) 283, 10.1016/j.physletb.2018.02.004, arXiv:1708.00373.
  * [15] CMS Collaboration, “Measurements of properties of the Higgs boson decaying to a W boson pair in pp collisions at $\sqrt{s}=$ 13 TeV”,  Phys. Lett. B 791 (2019) 96, 10.1016/j.physletb.2018.12.073, arXiv:1806.05246.
  * [16] ATLAS Collaboration, “Measurements of the Higgs boson production and decay rates and coupling strengths using pp collision data at $\sqrt{s}=7$ and 8 TeV in the ATLAS experiment”,  Eur. Phys. J. C 76 (2016) 6, 10.1140/epjc/s10052-015-3769-y, arXiv:1507.04548.
  * [17] CMS Collaboration, “Precise determination of the mass of the Higgs boson and tests of compatibility of its couplings with the standard model predictions using proton collisions at 7 and 8 TeV”,  Eur. Phys. J. C 75 (2015) 212, 10.1140/epjc/s10052-015-3351-7, arXiv:1412.8662.
  * [18] CMS Collaboration, “Study of the mass and spin-parity of the Higgs boson candidate via its decays to Z boson pairs”,  Phys. Rev. Lett. 110 (2013) 081803, 10.1103/PhysRevLett.110.081803, arXiv:1212.6639.
  * [19] ATLAS Collaboration, “Evidence for the spin-0 nature of the Higgs boson using ATLAS data”,  Phys. Lett. B 726 (2013) 120, 10.1016/j.physletb.2013.08.026, arXiv:1307.1432.
  * [20] ATLAS and CMS Collaborations, “Measurements of the Higgs boson production and decay rates and constraints on its couplings from a combined ATLAS and CMS analysis of the LHC pp collision data at $\sqrt{s}=7$ and 8 TeV”,  JHEP 08 (2016) 045, 10.1007/JHEP08(2016)045, arXiv:1606.02266.
  * [21] ATLAS and CMS Collaborations, “Combined measurement of the Higgs boson mass in pp collisions at $\sqrt{s}=7$ and 8 TeV with the ATLAS and CMS experiments”,  Phys. Rev. Lett. 114 (2015) 191803, 10.1103/PhysRevLett.114.191803, arXiv:1503.07589.
  * [22] CMS Collaboration, “Measurements of properties of the Higgs boson decaying into the four-lepton final state in pp collisions at $\sqrt{s}=13$ TeV”,  JHEP 11 (2017) 047, 10.1007/JHEP11(2017)047, arXiv:1706.09936.
  * [23] ATLAS Collaboration, “Combined measurement of differential and total cross sections in the $\PH\rightarrow\PGg\PGg$ and the $\PH\rightarrow\PZ\PZ^{*}\rightarrow 4\ell$ decay channels at $\sqrt{s}=13$ TeV with the ATLAS detector”,  Phys. Lett. B 786 (2018) 114, 10.1016/j.physletb.2018.09.019, arXiv:1805.10197.
  * [24] ATLAS Collaboration, “Measurement of the Higgs boson mass in the $\PH\rightarrow\PZ\PZ^{*}\rightarrow 4\ell$ and $\PH\rightarrow\PGg\PGg$ channels with $\sqrt{s}=13$ TeV pp collisions using the ATLAS detector”,  Phys. Lett. B 784 (2018) 345, 10.1016/j.physletb.2018.07.050, arXiv:1806.00242.
  * [25] CMS Collaboration, “Search for the Higgs boson decaying to two muons in proton-proton collisions at $\sqrt{s}=$ 13 TeV”,  Phys. Rev. Lett. 122 (2019) 021801, 10.1103/PhysRevLett.122.021801, arXiv:1807.06325.
  * [26] D. Ghosh, R. S. Gupta, and G. Perez, “Is the Higgs mechanism of fermion mass generation a fact? A Yukawa-less first-two-generation model”,  Phys. Lett. B 755 (2016) 504, 10.1016/j.physletb.2016.02.059, arXiv:1508.01501.
  * [27] F. J. Botella, G. C. Branco, M. N. Rebelo, and J. I. Silva-Marcos, “What if the masses of the first two quark families are not generated by the standard model Higgs boson?”,  Phys. Rev. D 94 (2016) 115031, 10.1103/PhysRevD.94.115031, arXiv:1602.08011.
  * [28] R. Harnik, J. Kopp, and J. Zupan, “Flavor violating Higgs decays”,  JHEP 03 (2013) 026, 10.1007/JHEP03(2013)026, arXiv:1209.1397.
  * [29] W. Altmannshofer et al., “Collider signatures of flavorful Higgs bosons”,  Phys. Rev. D 94 (2016) 115032, 10.1103/PhysRevD.94.115032, arXiv:1610.02398.
  * [30] ATLAS Collaboration, “Search for the decay of the Higgs boson to charm quarks with the ATLAS experiment”,  Phys. Rev. Lett. 120 (2018) 211802, 10.1103/PhysRevLett.120.211802, arXiv:1802.04329.
  * [31] CMS Collaboration, “A search for the standard model Higgs boson decaying to charm quarks”,  JHEP 03 (2020) 131, 10.1007/JHEP03(2020)131, arXiv:1912.01662.
  * [32] ATLAS Collaboration, “Direct constraint on the Higgs-charm coupling from a search for Higgs boson decays into charm quarks with the ATLAS detector”, 2022. arXiv:2201.11428. Submitted to EPJC.
  * [33] CMS Collaboration, “Precision luminosity measurement in proton-proton collisions at $\sqrt{s}=$ 13 TeV in 2015 and 2016 at CMS”,  Eur. Phys. J. C 81 (2021) 800, 10.1140/epjc/s10052-021-09538-2, arXiv:2104.01927.
  * [34] CMS Collaboration, “CMS luminosity measurement for the 2017 data-taking period at $\sqrt{s}$ = 13 TeV”, CMS Physics Analysis Summary CMS-PAS-LUM-17-004, 2018.
  * [35] CMS Collaboration, “CMS luminosity measurement for the 2018 data-taking period at $\sqrt{s}$ = 13 TeV”, CMS Physics Analysis Summary CMS-PAS-LUM-18-002, 2019.
  * [36] CMS Collaboration, “The CMS experiment at the CERN LHC”,  JINST 3 (2008) S08004, 10.1088/1748-0221/3/08/S08004.
  * [37] CMS Collaboration, “Performance of the CMS Level-1 trigger in proton-proton collisions at $\sqrt{s}=13$ TeV”,  JINST 15 (2020) P10017, 10.1088/1748-0221/15/10/P10017, arXiv:2006.10165.
  * [38] CMS Collaboration, “The CMS trigger system”,  JINST 12 (2017) P01020, 10.1088/1748-0221/12/01/P01020, arXiv:1609.02366.
  * [39] CMS Collaboration, “Electron and photon reconstruction and identification with the CMS experiment at the CERN LHC”,  JINST 16 (2021) P05014, 10.1088/1748-0221/16/05/P05014, arXiv:2012.06888.
  * [40] CMS Collaboration, “Performance of the CMS muon detector and muon reconstruction with proton-proton collisions at $\sqrt{s}=$ 13 TeV”,  JINST 13 (2018) P06015, 10.1088/1748-0221/13/06/P06015, arXiv:1804.04528.
  * [41] CMS Collaboration, “Description and performance of track and primary-vertex reconstruction with the CMS tracker”,  JINST 9 (2014) P10009, 10.1088/1748-0221/9/10/P10009, arXiv:1405.6569.
  * [42] CMS Collaboration, “Particle-flow reconstruction and global event description with the CMS detector”,  JINST 12 (2017) P10003, 10.1088/1748-0221/12/10/P10003, arXiv:1706.04965.
  * [43] CMS Collaboration, “Performance of reconstruction and identification of $\tau$ leptons decaying to hadrons and $\nu_{\tau}$ in pp collisions at $\sqrt{s}=$ 13 TeV”,  JINST 13 (2018) P10005, 10.1088/1748-0221/13/10/P10005, arXiv:1809.02816.
  * [44] CMS Collaboration, “Jet energy scale and resolution in the CMS experiment in pp collisions at 8 TeV”,  JINST 12 (2017) P02014, 10.1088/1748-0221/12/02/P02014, arXiv:1607.03663.
  * [45] CMS Collaboration, “Performance of missing transverse momentum reconstruction in proton-proton collisions at $\sqrt{s}=13$ TeV using the CMS detector”,  JINST 14 (2019) P07004, 10.1088/1748-0221/14/07/P07004, arXiv:1903.06078.
  * [46] M. Cacciari, G. P. Salam, and G. Soyez, “The anti-jet clustering algorithm”,  JHEP 04 (2008) 063, 10.1088/1126-6708/2008/04/063, arXiv:0802.1189.
  * [47] M. Cacciari, G. P. Salam, and G. Soyez, “FastJet user manual”,  Eur. Phys. J. C 72 (2012) 1896, 10.1140/epjc/s10052-012-1896-2, arXiv:1111.6097.
  * [48] CMS Collaboration, “Pileup mitigation at CMS in 13 TeV data”,  JINST 15 (2020) P09018, 10.1088/1748-0221/15/09/P09018, arXiv:2003.00503.
  * [49] D. Bertolini, P. Harris, M. Low, and N. Tran, “Pileup per particle identification”,  JHEP 10 (2014) 059, 10.1007/JHEP10(2014)059, arXiv:1407.6013.
  * [50] CMS Collaboration, “Mass regression of highly-boosted jets using graph neural networks”, CMS Detector Performance Note CMS-DP-2021-017, 2021.
  * [51] H. Qu and L. Gouskos, “Jet tagging via particle clouds”,  Phys. Rev. D 101 (2020) 056019, 10.1103/PhysRevD.101.056019, arXiv:1902.08570.
  * [52] M. Dasgupta, A. Fregoso, S. Marzani, and G. P. Salam, “Towards an understanding of jet substructure”,  JHEP 09 (2013) 029, 10.1007/JHEP09(2013)029, arXiv:1307.0007.
  * [53] 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, 10.1103/PhysRevLett.100.242001, arXiv:0802.2470.
  * [54] GEANT4 Collaboration, “ — A simulation toolkit”,  Nucl. Instrum. Meth. A 506 (2003) 250, 10.1016/S0168-9002(03)01368-8.
  * [55] P. Nason, “A new method for combining NLO QCD with shower Monte Carlo algorithms”,  JHEP 11 (2004) 040, 10.1088/1126-6708/2004/11/040, arXiv:hep-ph/0409146.
  * [56] S. Frixione, P. Nason, and C. Oleari, “Matching NLO QCD computations with parton shower simulations: the POWHEG method”,  JHEP 11 (2007) 070, 10.1088/1126-6708/2007/11/070, arXiv:0709.2092.
  * [57] S. Alioli, P. Nason, C. Oleari, and E. Re, “A general framework for implementing NLO calculations in shower Monte Carlo programs: the POWHEG BOX”,  JHEP 06 (2010) 043, 10.1007/JHEP06(2010)043, arXiv:1002.2581.
  * [58] K. Hamilton, P. Nason, and G. Zanderighi, “MINLO: Multi-scale improved NLO”,  JHEP 10 (2012) 155, 10.1007/JHEP10(2012)155, arXiv:1206.3572.
  * [59] G. Luisoni, P. Nason, C. Oleari, and F. Tramontano, “HW±/HZ + 0 and 1 jet at NLO with the POWHEG BOX interfaced to GoSam and their merging within MiNLO”,  JHEP 10 (2013) 083, 10.1007/JHEP10(2013)083, arXiv:1306.2542.
  * [60] LHC Higgs Cross Section Working Group, “Handbook of LHC Higgs cross sections: 4. deciphering the nature of the Higgs sector”,  CERN (2016) 10.23731/CYRM-2017-002, arXiv:1610.07922.
  * [61] G. Ferrera, M. Grazzini, and F. Tramontano, “Higher-order QCD effects for associated WH production and decay at the LHC”,  JHEP 04 (2014) 039, 10.1007/JHEP04(2014)039, arXiv:1312.1669.
  * [62] G. Ferrera, M. Grazzini, and F. Tramontano, “Associated ZH production at hadron colliders: the fully differential NNLO QCD calculation”,  Phys. Lett. B 740 (2015) 51, 10.1016/j.physletb.2014.11.040, arXiv:1407.4747.
  * [63] G. Ferrera, M. Grazzini, and F. Tramontano, “Associated WH production at hadron colliders: a fully exclusive QCD calculation at NNLO”,  Phys. Rev. Lett. 107 (2011) 152003, 10.1103/PhysRevLett.107.152003, arXiv:1107.1164.
  * [64] G. Ferrera, G. Somogyi, and F. Tramontano, “Associated production of a Higgs boson decaying into bottom quarks at the LHC in full NNLO QCD”,  Phys. Lett. B 780 (2018) 346, 10.1016/j.physletb.2018.03.021, arXiv:1705.10304.
  * [65] O. Brein, R. V. Harlander, and T. J. E. Zirke, “vh@nnlo — Higgs Strahlung at hadron colliders”,  Comput. Phys. Commun. 184 (2013) 998, 10.1016/j.cpc.2012.11.002, arXiv:1210.5347.
  * [66] R. V. Harlander, S. Liebler, and T. Zirke, “Higgs Strahlung at the Large Hadron Collider in the 2-Higgs-doublet model”,  JHEP 02 (2014) 023, 10.1007/JHEP02(2014)023, arXiv:1307.8122.
  * [67] A. Denner, S. Dittmaier, S. Kallweit, and A. Mück, “HAWK 2.0: A Monte Carlo program for Higgs production in vector-boson fusion and Higgs strahlung at hadron colliders”,  Comput. Phys. Commun. 195 (2015) 161, 10.1016/j.cpc.2015.04.021, arXiv:1412.5390.
  * [68] J. Alwall 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, 10.1007/JHEP07(2014)079, arXiv:1405.0301.
  * [69] S. Frixione, P. Nason, and G. Ridolfi, “A positive-weight next-to-leading-order Monte Carlo for heavy flavour hadroproduction”,  JHEP 09 (2007) 126, 10.1088/1126-6708/2007/09/126, arXiv:0707.3088.
  * [70] E. Re, “Single-top Wt-channel production matched with parton showers using the POWHEG method”,  Eur. Phys. J. C 71 (2011) 1547, 10.1140/epjc/s10052-011-1547-z, arXiv:1009.2450.
  * [71] R. Frederix, E. Re, and P. Torrielli, “Single-top $t$-channel hadroproduction in the four-flavour scheme with POWHEG and aMC@NLO”,  JHEP 09 (2012) 130, 10.1007/JHEP09(2012)130, arXiv:1207.5391.
  * [72] S. Alioli, P. Nason, C. Oleari, and E. Re, “NLO single-top production matched with shower in POWHEG: $s$\- and $t$-channel contributions”,  JHEP 09 (2009) 111, 10.1088/1126-6708/2009/09/111, arXiv:0907.4076. [Erratum: 10.1007/JHEP02(2010)011].
  * [73] M. Czakon and A. Mitov, “Top++: A program for the calculation of the top-pair cross-section at hadron colliders”,  Comput. Phys. Commun. 185 (2014) 2930, 10.1016/j.cpc.2014.06.021, arXiv:1112.5675.
  * [74] M. Czakon et al., “Top-pair production at the LHC through NNLO QCD and NLO EW”,  JHEP 10 (2017) 186, 10.1007/JHEP10(2017)186, arXiv:1705.04105.
  * [75] T. Melia, P. Nason, R. Rontsch, and G. Zanderighi, “W+W-, WZ and ZZ production in the POWHEG BOX”,  JHEP 11 (2011) 078, 10.1007/JHEP11(2011)078, arXiv:1107.5051.
  * [76] M. Grazzini et al., “NNLO QCD + NLO EW with Matrix+OpenLoops: precise predictions for vector-boson pair production”,  JHEP 02 (2020) 087, 10.1007/JHEP02(2020)087, arXiv:1912.00068.
  * [77] NNPDF Collaboration, “Parton distributions for the LHC Run II”,  JHEP 04 (2015) 040, 10.1007/JHEP04(2015)040, arXiv:1410.8849.
  * [78] NNPDF Collaboration, “Parton distributions from high-precision collider data”,  Eur. Phys. J. C 77 (2017) 663, 10.1140/epjc/s10052-017-5199-5, arXiv:1706.00428.
  * [79] T. Sjöstrand et al., “An introduction to PYTHIA 8.2”,  Comput. Phys. Commun. 191 (2015) 159, 10.1016/j.cpc.2015.01.024, arXiv:1410.3012.
  * [80] CMS Collaboration, “Event generator tunes obtained from underlying event and multiparton scattering measurements”,  Eur. Phys. J. C 76 (2016) 155, 10.1140/epjc/s10052-016-3988-x, arXiv:1512.00815.
  * [81] CMS Collaboration, “Extraction and validation of a new set of CMS PYTHIA8 tunes from underlying-event measurements”,  Eur. Phys. J. C 80 (2020) 4, 10.1140/epjc/s10052-019-7499-4, arXiv:1903.12179.
  * [82] R. Frederix and S. Frixione, “Merging meets matching in MC@NLO”,  JHEP 12 (2012) 061, 10.1007/JHEP12(2012)061, arXiv:1209.6215.
  * [83] J. Alwall et al., “Comparative study of various algorithms for the merging of parton showers and matrix elements in hadronic collisions”,  Eur. Phys. J. C 53 (2008) 473, 10.1140/epjc/s10052-007-0490-5, arXiv:0706.2569.
  * [84] T. Plehn, G. P. Salam, and M. Spannowsky, “Fat jets for a light Higgs”,  Phys. Rev. Lett. 104 (2010) 111801, 10.1103/PhysRevLett.104.111801, arXiv:0910.5472.
  * [85] Y. Wang et al., “Dynamic Graph CNN for learning on point clouds”,  ACM Trans. Graph. 38 (2019) 146, 10.1145/3326362, arXiv:1801.07829.
  * [86] CMS Collaboration, “Identification of highly Lorentz-boosted heavy particles using graph neural networks and new mass decorrelation techniques”, CMS Detector Performance Note CMS-DP-2020-002, 2020.
  * [87] A. Butter et al., “The machine learning landscape of top taggers”,  SciPost Phys. 7 (2019) 014, 10.21468/SciPostPhys.7.1.014, arXiv:1902.09914.
  * [88] CMS Collaboration, “Identification of heavy, energetic, hadronically decaying particles using machine-learning techniques”,  JINST 15 (2020) P06005, 10.1088/1748-0221/15/06/P06005, arXiv:2004.08262.
  * [89] CMS Collaboration, “Calibration of the mass-decorrelated ParticleNet tagger for boosted $\bbbar$ and $\ccbar$ jets using LHC Run 2 data”, CMS Detector Performance Note CMS-DP-2022-005, 2022.
  * [90] E. Bols et al., “Jet flavour classification using DeepJet”,  JINST 15 (2020) P12012, 10.1088/1748-0221/15/12/P12012, arXiv:2008.10519.
  * [91] CMS Collaboration, “Performance of the DeepJet b tagging algorithm using 41.9/fb of data from proton-proton collisions at 13 TeV with Phase 1 CMS detector”, CMS Detector Performance Note CMS-DP-2018-058, 2018.
  * [92] CMS Collaboration, “A new calibration method for charm jet identification validated with proton-proton collision events at $\sqrt{s}=13\TeV$”,  JINST 17 (2022) P03014, 10.1088/1748-0221/17/03/P03014, arXiv:2111.03027.
  * [93] CMS Collaboration, “A deep neural network for simultaneous estimation of b jet energy and resolution”,  Comput. Softw. Big Sci. 4 (2020) 10, 10.1007/s41781-020-00041-z, arXiv:1912.06046.
  * [94] CMS Collaboration, “Observation of Higgs boson decay to bottom quarks”,  Phys. Rev. Lett. 121 (2018) 121801, 10.1103/PhysRevLett.121.121801, arXiv:1808.08242.
  * [95] ATLAS and CMS Collaborations, and the LHC Higgs Combination Group, “Procedure for the LHC Higgs boson search combination in Summer 2011”, Technical Report CMS-NOTE-2011-005, ATL-PHYS-PUB-2011-11, 2011.
  * [96] G. Cowan, K. Cranmer, E. Gross, and O. Vitells, “Asymptotic formulae for likelihood-based tests of new physics”,  Eur. Phys. J. C 71 (2011) 1554, 10.1140/epjc/s10052-011-1554-0, arXiv:1007.1727. [Erratum: 10.1140/epjc/s10052-013-2501-z].
  * [97] T. Junk, “Confidence level computation for combining searches with small statistics”,  Nucl. Instrum. Meth. A 434 (1999) 435, 10.1016/S0168-9002(99)00498-2, arXiv:hep-ex/9902006.
  * [98] A. L. Read, “Presentation of search results: The technique”,  J. Phys. G 28 (2002) 2693, 10.1088/0954-3899/28/10/313.
  * [99] “HEPData record for this analysis”, 2022.  10.17182/hepdata.128675.
  * [100] LHC Higgs Cross Section Working Group, “Handbook of LHC Higgs cross sections: 3. Higgs properties”, CERN Yellow Report CERN-2013-004, 2013.  10.5170/CERN-2013-004, arXiv:1307.1347.