# Segmentation of Cardiac Structures via Successive Subspace Learning with
Saab Transform from Cine MRI

Xiaofeng Liu1, Fangxu Xing1, Hanna K. Gaggin2, Weichung Wang3, C.-C. Jay Kuo4,
Georges El Fakhri1, Jonghye Woo1, 1X. Liu, F. Xing, G. El Fakhri, and J. Woo
are with the Gordon Center for Medical Imaging, Dept. of Radiology,
Massachusetts General Hospital and Harvard Medical School, Boston, MA, USA2H.
Gaggin is with Division of Cardiology, Corrigan Minehan Heart Center and Dept.
of Medicine, Massachusetts General Hospital and Harvard Medical School,
Boston, MA, USA3W. Wang is with Institute of Applied Mathematical Sciences,
National Taiwan University, Taipei, Taiwan4C.-C. J. Kuo is with Dept. of
Electrical and Computer Engineering, University of Southern California, Los
Angeles, CA, USA.

###### Abstract

Assessment of cardiovascular disease (CVD) with cine magnetic resonance
imaging (MRI) has been used to non-invasively evaluate detailed cardiac
structure and function. Accurate segmentation of cardiac structures from cine
MRI is a crucial step for early diagnosis and prognosis of CVD, and has been
greatly improved with convolutional neural networks (CNN). There, however, are
a number of limitations identified in CNN models, such as limited
interpretability and high complexity, thus limiting their use in clinical
practice. In this work, to address the limitations, we propose a lightweight
and interpretable machine learning model, successive subspace learning with
the subspace approximation with adjusted bias (Saab) transform, for accurate
and efficient segmentation from cine MRI. Specifically, our segmentation
framework is comprised of the following steps: (1) sequential expansion of
near-to-far neighborhood at different resolutions; (2) channel-wise subspace
approximation using the Saab transform for unsupervised dimension reduction;
(3) class-wise entropy guided feature selection for supervised dimension
reduction; (4) concatenation of features and pixel-wise classification with
gradient boost; and (5) conditional random field for post-processing.
Experimental results on the ACDC 2017 segmentation database, showed that our
framework performed better than state-of-the-art U-Net models with 200$\times$
fewer parameters in delineating the left ventricle, right ventricle, and
myocardium, thus showing its potential to be used in clinical practice.

Clinical relevance— Delineation of the left ventricular cavity, myocardium,
and right ventricle from cardiac MR images is a common clinical task to
establish diagnosis and prognosis of CVD.

##  I Introduction

Cardiovascular disease (CVD) continues to be the cause of the largest portion
of morbidity and mortality globally, accounting for over 18 million deaths
globally [1]. Assessment of CVD with cine magnetic resonance imaging (MRI) has
been shown to provide a non-invasive way to evaluate the detailed morphology
and function of the heart. In particular, cine MRI is considered to be the
most accurate imaging modality for assessing various quantitative parameters
with important prognostic implications.

Segmentation of the left ventricle (LV), right ventricle (RV), and myocardium
(MYO) from cardiac cine MR images plays an important role in characterizing
clinically important parameters [2], such as ejection fraction (EF), end
diastolic volume (EDV), end systolic volume (ESV), and myocardial mass. These
parameters, in turn, can be used to identify disease phenotypes, stratify
disease risks, and develop diagnostic and prognostic tools [3]. In clinical
practice, semi-automated segmentation is still predominantly used, partly due
to the lack of fully-automated and accurate segmentation tools [4], which is
time-consuming and suffers from inter-observer variability.

With the recent progress of deep learning [5], numerous convolutional neural
networks (CNN) models, e.g., U-Net [6], have been developed, demonstrating
their accuracy in many medical image analysis tasks [7]. While deep learning
has achieved impressive results for segmentation and classification, a number
of challenges arise in developing and deploying deep learning models for
clinical applications [7]. First, CNN models typically require a large number
of labeled training datasets [5]. Sparse and inaccurate labels caused by
privacy issues and the high cost of labeling, however, lead to difficulty in
collecting sufficient and high-quality training sample datasets [5]; with the
limited training datasets, an accurate model fitting at the training stage is
challenging. Recently, to address this, efforts have been made to generate
samples using data augmentation or adversarial training [8], which, however,
results in an unavoidable problem of appearance shift between real and
generated data. Second, importantly, many CNN models are seen as a “black-box”
model [9, 5]. Accordingly, CNN models remain largely elusive how a particular
CNN model makes a decision and when it can be trusted. Therefore, it is
crucial to develop an explainable model that works with a limited number of
datasets for clinical applications.

To address the aforementioned challenges, in this work, we propose to develop
a lightweight, interpretable, and fully-automated segmentation framework with
successive subspace learning (SSL) [10]. Specifically, our framework is
comprised of the following steps: (1) sequential expansion of near-to-far
neighborhood at different resolutions; (2) channel-wise subspace approximation
using the subspace approximation with adjusted bias (Saab) transform for
unsupervised dimension reduction; (3) a novel class-wise entropy guided
feature selection for supervised dimension reduction; (4) concatenation of
features and pixel-wise classification with gradient boost; and (5)
conditional random field for post-processing.

To the best of our knowledge, this is the first attempt at exploring the SSL
framework with the Saab transform for a segmentation task. Our framework is
lightweight and interpretable, yet achieving a superior segmentation
performance with 200$\times$ fewer parameters, compared with state-of-the-art
U-Net models.

##  II Methodology

###  II-A Fundamentals of SSL and Saab Transform

Inspired by the recent stacked design of CNN models, the SSL principle [10]
has been targeted for classifying 2D natural images (e.g., PixelHop [11, 12]),
3D MR images [13], and point clouds (e.g., PointHop [12]). In each layer of
SSL, the Saab transform [14], a variant of Principal Component Analysis (PCA),
is used as an alternative to nonlinear activation, thereby alleviating the
sign confusion problem [9]. Furthermore, the Saab transform is deemed more
interpretable than nonlinear activation functions in CNNs [14, 15], as the
model parameters are computed stage-by-stage in a feedforward manner, without
backpropagation. Accordingly, the training of our SSL-based method is more
efficient and interpretable than that of CNN models [11].

###  II-B Our Saab-based SSL segmentation Framework

In this work, we have a 2D MR image ${\bf x}\in\mathbb{R}^{H\times W\times 1}$
and its corresponding label ${\bf y}\in\mathbb{R}^{H\times W\times 4}$, where
$H$ and $W$ denote the horizontal and vertical dimensions, respectively. The
channel of the gray-value sample is 1, and the label of each pixel is encoded
as a four-dimensional one-hot vector for four tissue classes. The architecture
of our framework is illustrated in Fig. 1, as detailed below.

####  II-B1 Module 1: Unsupervised Feature Selection

We first construct $I$ cascade SSL units and $I-1$ max-pooling operations to
extract the attributes at different spatial scales in the unsupervised Module
1. Similar to PixelHop [11], in each SSL unit, we construct the neighboring
region on the $H\times W$ plane. For instance, in the first SSL unit, for the
single-channel data, we construct the $3\times 3$ region for each pixel
position. Each of them is then flattened to a 9-dimensional vector. With a
padding operation, x is transformed to a cubic with the size of $H\times
W\times 9$. Then, the Saab transform is used for unsupervised dimension
reduction in the channel direction. Each 9-dimensional vector is mapped to a
$F_{1}$-dimensional feature vector, where $F_{1}$ is a hyperparameter to
control the output dimension of the first PixelHop unit.

Specifically, the terms, direct current (DC) and alternating current (AC), are
adopted from the circuit theory. In the first Saab transform, we configure one
DC and $F_{1}-1$ AC anchor vectors with the size of $H\times W\times 1$. Then,
the $c$-th dimension of $f$ can be an affine transform of $x$, i.e.,

$f_{c}=a_{c}^{T}x+b_{c},~{}~{}~{}c=0,1,\cdots,F_{1}-1,$ (1)

and the Saab transform has a special design of the anchor vector
$a_{c}\in\mathbb{R}^{1\times(H\times W\times 1)}$ and the bias term
$b_{c}\in\mathbb{R}$ [14]. Similar to [14], we can set $b_{c}\equiv
d\sqrt{F_{1}},d\in\mathbb{R}$, and divide the anchor vector into two
categories:

$\displaystyle\bullet~{}{\text{DC anchor
vector}}~{}~{}a_{0}=\frac{1}{\sqrt{H\times W\times 1}}(1,\cdots,1)^{T},$ (2)
$\displaystyle\bullet~{}{\text{AC anchor
vector}}~{}~{}a_{c},~{}~{}c=1,\cdots,F_{1}-1.$

After computing $f_{1}\in\mathbb{R}^{H\times W\times F_{1}}$, we half its
spacial size with the max-pooling operation to
${\frac{H}{2}\times\frac{W}{2}\times F_{1}}$ and send to the next SSL unit.
With the multi-channel input, the neighborhood construction involves $3\times
3\times F_{1}$ region at each pixel position. Then, the neighborhood union is
flattened to a vector, which is further processed by the Saab transform for
dimension reduction. The detailed structure of our module 1 is provided in
Table I.

With the cascaded SSL units, the neighborhood union is correlated with more
pixels of $\bf x$ to extract global information. This process is similar to
CNN models in that a larger reception field is achieved in the deeper layers.

####  II-B2 Module 2: Supervised Feature Selection

In what follows, we resort to the supervised dimension reduction based on
class-wise entropy-guided feature selection to tailor the discriminative
feature for our segmentation task.

Figure 1: Illustration of our proposed framework using the Saab transform,
which consists of 3 modules.

Because of the resolution deduction in each unit, we have different spatial
size of the extracted features. The features in the later units correspond to
a larger reception field (i.e., more pixels) in ${\bf x}$ and ${\bf y}$. To
match the features in these units with the original pixels, we resize
$f_{2},\cdots,f_{I}$ to the size of $f_{1}$ and denote as
$\overline{f_{2}},\cdots,\overline{f_{I}}$. Therefore, we have
$\overline{f_{1}}\in\mathbb{R}^{H\times W\times C_{1}}$,
$\overline{f_{2}}\in\mathbb{R}^{H\times W\times C_{2}}$, $\cdots$,
$\overline{f_{I}}\in\mathbb{R}^{H\times W\times C_{I}}$.

Because of the disparate importance, depending on the different channels for
the segmentation decision, it is necessary to make supervised feature
selection. In related developments, PixelHop++ [16] proposes to classify each
channel with the size of $\mathbb{R}^{H\times W}$ and select the channels with
low cross-entropy score. However, it is not applicable to segmentation as a
channel selection, since the label in the segmentation task is pixel-wise and
the feature of a pixel in each channel is only a scalar, making it challenging
to be used as a feature for a classifier.

Instead, we propose to select the channel with the small entropy of each
class. Specifically, we would encourage the feature of a pixel in each channel
to be similar, if the label of the corresponding pixels is the same class. We
denote the feature of a pixel in each channel $p_{i}^{c}$ for the $i$-th pixel
of a class in the $c$-th channel. The entropy of a sample can be:

$\displaystyle
H=\sum_{j=1}^{4}H_{j},~{}~{}H_{j}=-\sum_{i}p_{i}^{c}{\text{log}}p_{i}^{c},$
(3)

where we use $j$ to index the four classes in our segmentation task. After
calculating the entropy of four classes for each channel, we rank the entropy
in descending order. Then, we select the top 80% channels for the subsequent
pixel-wise classification task.

TABLE I: The detailed structure of our 4 consecutive SSL units

Input Size | Type | Filter Shape
---|---|---
$[224\times 224\times 1]$ | Saab Trans |  $F_{1}$ kernels of $3\times 3$
$[224\times 224\times F_{1}]$ | MaxPool | (2$\times$2)-(1$\times$1)
$[112\times 112\times F_{1}]$ | Saab Trans |  $F_{2}$ kernels of $3\times 3$ for F1 channels
$[112\times 112\times F_{2}]$ | MaxPool | (2$\times$2)-(1$\times$1)
$[56\times 56\times F_{2}]$ | Saab Trans |  $F_{3}$ kernels of $3\times 3$ for F2 channels
$[56\times 56\times F_{3}]$ | MaxPool | (2$\times$2)-(1$\times$1)
$[28\times 28\times F_{3}]$ | Saab |  $F_{4}$ kernels of $3\times 3$ for F3 channels

####  II-B3 Module 3: Information fusion for segmentation and post-processing

With the extracted features
$\overline{f_{2}^{\prime}},\cdots,\overline{f_{I}^{\prime}}$ with both the
Saab transform and class-wise entropy guided selection, we concatenate them
along with the channel dimension to get the feature $f\in\mathbb{R}^{H\times
W\times C}$. The channel dimension $C$ is the sum of all channels in
$\overline{f_{2}^{\prime}},\cdots,\overline{f_{I}^{\prime}}$. Each feature
vector on the $H\times W$ plane of $f$ corresponding to an original pixel in
${\bf x}$ or ${\bf y}$. Then, we carry out the pixel-wise classification for
each of $C$ dimensional features with a classifier. We empirically choose the
extreme gradient boosting (XGBoost) [17], which is an optimized distributed
gradient boosting library designed to be highly efficient and flexible.
XGBoost is trained to learn the correlation of pixel-wise feature and ground
truth pixel class label in our training set.

We note that with a limited number of SSL units, it is challenging to support
the reception field to cover all of the pixels for global perception. In
contrast, too many SSL units will lead to very low resolution in the later
units, which is not sufficient to support the pixel-wise segmentation. In
addition, the channel size of the later units will be very large, leading to a
long and indiscriminative feature vector, which can distract the pixel-wise
classification.

To balance this conflict, we propose to validate the most effective $I$ and
adopt the well-established post-processing tool of conditional random field
(CRF) to further refine the segmentation results and get the final results of
our framework $\tilde{y}=CRF(\hat{y})$.

##  III Experiments

To demonstrate the performance of our Saab transform-based SSL framework, we
validated it on the Automated Cardiac Diagnosis Challenge (ACDC 2017)
database, which contains 100 subjects. The cine MRI short-axis slices were
acquired with 1.5T or 3.0T MRI scanners. The acquired cine MRI short-axis
slices covered the LV, RV, and MYO from the base (upper slice) to the apex
(lower slice), with 5–8 mm slice thickness, 5 or 10 mm inter-slice gap and the
spatial resolution of 1.37–1.68 $mm^{2}$.

Figure 2: Comparison of the segmentation results with 60 subjects for
training. TABLE II: Comparison of the Dice score with 50 training subjects

Methods | Parameters | RV | MYO | LV | Average
---|---|---|---|---|---
U-Net | 5.88 M | 81.45% | 78.54% | 90.07% | 83.35%
AttenUNet | 6.40 M | 81.02% | 78.40% | 89.32% | 82.91%
SSL | 0.03 M | 82.91% | 81.83% | 90.62% | 85.12%

TABLE III: Comparison of the Dice score with 60 training subjects

Methods | Parameters | RV | MYO | LV | Average
---|---|---|---|---|---
U-Net | 5.88 M | 85.54% | 78.81% | 92.15% | 85.50%
AttenUNet | 6.40 M | 85.64% | 77.37% | 91.29% | 84.76%
SSL | 0.03 M | 83.89% | 82.57% | 91.62% | 86.03%

For each patient, the delineations of the LV, RV, and MYO, were obtained by
two clinical experts. On average, each subject had about 27 labeled slices. We
reported the average Dice similarity score with 30 subjects for testing, 10
subjects for validation, and 50 or 60 subjects for training.

###  III-A Implementation Details

All the experiments were implemented using Python on a server with a Xeon E5
v4 CPU/Nvidia Tesla V100 GPU with 128GB memory. We also used the widely
adopted deep learning library, Pytorch, to implement U-Net [6] and AttnUNet
[18]. For a fair comparison, we resized all of the slices to $224\times
224\times 1$, which was consistent with the input of the U-Net models.

We empirically used four SSL units and set $F_{1}$=5, $F_{2}$=10, $F_{3}$=30,
and $F_{4}$=100. We note that the number of the Saab AC filters in the
unsupervised dimension reduction procedure controls the preserved energy
ratio.

###  III-B Experimental Results

Fig. 2 shows the segmentation results of U-Net with ResNet50 backbone and our
SSL framework. We can see that SSL is able to achieve comparable or even
better performance than the widely used U-Net models.

For quantitative evaluation, we compared the Dice similarity score in Tables
II and III, which used 50 or 60 subjects for training, respectively. Note that
the larger Dice similarity score indicates the better segmentation
performance. The best results are bolded. With 50 subjects for training, our
SSL framework outperformed U-Net [6] and attention-based U-Net [18] in all of
the three classes. We can observe that with relatively limited training
datasets, the performance of the CNN models is inferior to our framework. In
addition, the statistics of the network parameters are provided and compared
in Table II. We can see that the number of parameters of our SSL framework was
about 200 times fewer than the popular U-Net structures [6, 18]. The much
fewer parameters can largely alleviate the difficulty of a small number of
training datasets. In the case of using 60 subjects for training, our SSL
framework achieved better performance than the U-Net based methods in the
average Dice similarity score.

###  III-C Sensitivity Analysis and Ablation Study

Figure 3: Sensitivity analysis of the number of SSL units for the case of
using 60 subjects in training.

With four SSL units, we achieved a state-of-the-art Dice similarity score in
both 50 and 60 training subjects settings. The number of SSL units is
important for our segmentation framework to balance the efficiency and
perception area. The low resolution can be challenging to provide accurate
information for fine-grained pixel-wise classification. In Fig. 3, we have
shown the detailed sensitivity study using different SSL units. The standard
deviation was computed with five random choices of training and validation
splits. The class-wise entropy-guided feature selection was developed to
simplify the subsequent classification modules. In addition, CRF was applied
as a post-processing step. To demonstrate their effectiveness, we provide the
ablation study in Table IV and the effect of CRF in Fig. 4.

TABLE IV: Results of Ablation Studies with 60 training subjects

Methods | Average DICE
---|---
SSL | 86.03%
SSL without CRF | 84.76%
SSL without entropy-guided feature selection | 85.91%

##  IV Conclusion

In this work, we presented a lightweight, interpretable, and fully-automated
SSL framework with the Saab transform to segment the LV, RV, and MYO from cine
MRI. A novel class-wise entropy-guided feature selection was proposed to
achieve accurate segmentation. Our thorough experiments carried out using the
ACDC 2017 database with different number of training subjects demonstrated
that our framework achieved a superior performance, compared with the U-Net-
based approaches, with about 200$\times$ fewer parameters.

Figure 4: Comparison of with or without CRF post-processing.

## References

  * [1] G. A. Roth, G. A. Mensah, and V. Fuster, “The Global Burden of Cardiovascular Diseases and Risks: A Compass for Global Action,” _American College of Cardiology_ , pp. 2980–2981, 2020.
  * [2] H. K. Gaggin and J. L. Januzzi Jr, “Biomarkers and diagnostics in heart failure,” _Biochimica et Biophysica Acta (BBA)-Molecular Basis of Disease_ , vol. 1832, no. 12, pp. 2442–2450, 2013.
  * [3] A. Ammar, O. Bouattane, and M. Youssfi, “Automatic cardiac cine mri segmentation and heart disease classification,” _Computerized Medical Imaging and Graphics_ , vol. 88, p. 101864, 2021.
  * [4] O. Bernard, A. Lalande, C. Zotti, F. Cervenansky, X. Yang, P.-A. Heng, I. Cetin, K. Lekadir, O. Camara, M. A. G. Ballester _et al._ , “Deep learning techniques for automatic MRI cardiac multi-structures segmentation and diagnosis: is the problem solved?” _IEEE transactions on medical imaging_ , vol. 37, no. 11, pp. 2514–2525, 2018.
  * [5] I. Goodfellow, Y. Bengio, A. Courville, and Y. Bengio, _Deep learning_. MIT press Cambridge, 2016, vol. 1, no. 2.
  * [6] O. Ronneberger, P. Fischer, and T. Brox, “U-net: Convolutional networks for biomedical image segmentation,” in _International Conference on Medical image computing and computer-assisted intervention_. Springer, 2015, pp. 234–241.
  * [7] D. Shen, G. Wu, and H.-I. Suk, “Deep learning in medical image analysis,” _Annual review of biomedical engineering_ , vol. 19, pp. 221–248, 2017.
  * [8] X. Liu, B. V. Kumar, P. Jia, and J. You, “Hard negative generation for identity-disentangled facial expression recognition,” _Pattern Recognition_ , vol. 88, pp. 1–12, 2019.
  * [9] C.-C. J. Kuo, “Understanding convolutional neural networks with a mathematical model,” _Journal of Visual Communication and Image Representation_ , vol. 41, pp. 406–413, 2016.
  * [10] M. Rouhsedaghat, M. Monajatipoor, Z. Azizi, and C.-C. J. Kuo, “Successive subspace learning: An overview,” _arXiv_ , 2021.
  * [11] Y. Chen and C.-C. J. Kuo, “Pixelhop: A successive subspace learning (SSL) method for object recognition,” _Journal of Visual Communication and Image Representation_ , p. 102749, 2020.
  * [12] M. Zhang, H. You, P. Kadam, S. Liu, and C.-C. J. Kuo, “Pointhop: An explainable machine learning method for point cloud classification,” _IEEE Transactions on Multimedia_ , 2020.
  * [13] X. Liu, F. Xing, C. Yang, C.-C. J. Kuo, S. Babu, G. E. Fakhri, T. Jenkins, and J. Woo, “Voxelhop: Successive subspace learning for als disease classification using structural MRI,” _arXiv preprint arXiv:2101.05131_.
  * [14] C.-C. J. Kuo, M. Zhang, S. Li, J. Duan, and Y. Chen, “Interpretable convolutional neural networks via feedforward design,” _Journal of Visual Communication and Image Representation_ , vol. 60, pp. 346–359, 2019.
  * [15] F. Fan, J. Xiong, and G. Wang, “On interpretability of artificial neural networks,” _arXiv preprint arXiv:2001.02522_ , 2020.
  * [16] Y. Chen, M. Rouhsedaghat, S. You, R. Rao, and C.-C. J. Kuo, “Pixelhop++: A small successive-subspace-learning-based (SSL-based) model for image classification,” _arXiv preprint arXiv:2002.03141_ , 2020.
  * [17] T. Chen, T. He, M. Benesty, V. Khotilovich, Y. Tang, H. Cho _et al._ , “Xgboost: extreme gradient boosting,” _R package version 0.4-2_ , vol. 1, no. 4, 2015.
  * [18] J. Schlemper, O. Oktay, M. Schaap, M. Heinrich, B. Kainz, B. Glocker, and D. Rueckert, “Attention gated networks: Learning to leverage salient regions in medical images,” _Medical image analysis_ , vol. 53, pp. 197–207, 2019\.


EUROPEAN ORGANIZATION FOR NUCLEAR RESEARCH (CERN)

 CERN-EP-2023-141 LHCb-PAPER-2023-010 Augest 25, 2023

Measurement of $Z$ boson production cross-section in $pp$ collisions

at $\sqrt{s}$ $=5.02\text{\,Te\kern-2.07413ptV}$

LHCb collaboration†††Authors are listed at the end of this paper.

The first measurement of the $Z$ boson production cross-section at centre-of-
mass energy $\sqrt{s}=5.02\text{\,Te\kern-1.00006ptV}$ in the forward region
is reported, using $pp$ collision data collected by the LHCb experiment in
year 2017, corresponding to an integrated luminosity of $100\pm
2\text{\,pb}^{-1}$. The production cross-section is measured for final-state
muons in the pseudorapidity range $2.0<\eta<4.5$ with transverse momentum
$p_{\mathrm{T}}>20\text{\,Ge\kern-1.00006ptV\\!/}c$. The integrated cross-
section is determined to be

$\sigma_{{Z}\\!\rightarrow{\mu^{+}\mu^{-}}}=39.6\pm 0.7\text{\,(stat)}\pm
0.6\text{\,(syst)}\pm 0.8\text{\,(lumi)}\,\text{\,pb}$

for the di-muon invariant mass in the range
$60<M_{\mu\mu}<120\text{\,Ge\kern-1.00006ptV\\!/}c^{2}$. This result and the
differential cross-section results are in good agreement with theoretical
predictions at next-to-next-to-leading order in the strong coupling.

Based on a previous LHCb measurement of the $Z$ boson production cross-section
in $p$Pb collisions at $\sqrt{s_{NN}}=5.02\text{\,Te\kern-1.00006ptV}$, the
nuclear modification factor $R_{{p}\mathrm{Pb}}$ is measured for the first
time at this energy. The measured values are
$1.2^{+0.5}_{-0.3}\text{\,(stat)}\pm 0.1\text{\,(syst)}$ in the forward region
($1.53<y^{*}_{\mu}<4.03$) and $3.6^{+1.6}_{-0.9}\text{\,(stat)}\pm
0.2\text{\,(syst)}$ in the backward region ($-4.97<y^{*}_{\mu}<-2.47$), where
$y^{*}_{\mu}$ represents the muon rapidity in the centre-of-mass frame.

Submitted to JHEP

© 2024 CERN for the benefit of the LHCb collaboration. CC BY 4.0 licence.

##  1 Introduction

Measurements of the $pp\rightarrow{Z}\rightarrow\mu^{+}\mu^{-}$ process111The
production process should be interpreted as $pp\rightarrow
Z/\gamma^{*}\rightarrow{\mu^{+}}{\mu^{-}}$ in the strict sense. In this
article, the label $Z$ boson is defined to also include contributions from
virtual photons and the interference between the $Z$ boson and the virtual
photon. provide important tests of the quantum chromodynamics (QCD) and
electroweak (EW) sectors of the Standard Model. High precision measurements of
the $Z$ boson production cross-section at various experiments provide valuable
insights for testing theoretical predictions [1, 2, 3, 4, 5, 6, 7, 8, 9],
which are further validated through their comparison with precision
perturbative QCD calculations up to the order of $\alpha_{s}^{3}$ [8, 9].

High precision measurements of the $Z$ boson production cross-sections at
different rapidities at LHCb provide important constraints on parton
distribution functions (PDFs). Results from deep inelastic scattering and
hadronic collisions [10, 11, 12, 13, 14, 15, 16, 17, 18, 19, 20, 21, 22, 23,
24, 25], parameterized in terms of the Bjorken variable, $x$, indicating the
fraction of the proton momentum carried by a single parton, are used in the
global fits of the PDFs. However, these measurements provide limited
information on the PDFs in the very large ($x$ $\sim 0.8$) or very small ($x$
$\sim 10^{-4}$) Bjorken-$x$ regions. This leads to large uncertainties on the
PDFs, and on the theoretical predictions that make use of them. Forward
acceptance of the LHCb detector covers a unique region of phase space,
allowing the measurement of highly boosted $Z$ boson candidates. Measurements
within this region are sensitive to both large and small Bjorken-$x$ values.
Previous measurements of single $W$ and $Z$ production by the LHCb
collaboration [26, 27, 28, 29, 30, 31, 32] have been included in PDF
calculations [33, 34, 35, 36, 37, 38, 39, 40, 41]; the LHCb data contribute
significantly to the precision of these determinations of the quark PDFs at
large and small values of $x$.

In addition, the $Z$ boson production cross-section is useful for constraining
nuclear PDFs (nPDFs), providing a clean probe of nuclear-matter effects in the
initial state. These effects are typically studied in terms of the nuclear
modification factor, $R_{{p}\mathrm{Pb}}$, defined as the ratio of the yield
observed in $p$Pb collisions to that in $pp$ collisions, scaled by the mean
number of nucleon-nucleon interactions. This quantity is used to study the
modification of particle production in heavy-ion collisions compared to $pp$
collisions. The LHCb experiment published the first inclusive $Z$ production
result in $p$Pb collisions at a nucleon-nucleon centre-of-mass energy of
$\sqrt{s_{NN}}=5.02\text{\,Te\kern-1.00006ptV}$ [42]. However, the nuclear
modification factors are not reported for this collision energy due to the
absence of a measurement for the cross-section in $pp$ collisions. With this
new measurement for the $Z$ boson production cross-section in $pp$ collisions
at $\sqrt{s}=5.02\text{\,Te\kern-1.00006ptV}$, the nuclear modification
factors in the forward region and backward regions are been reported for the
first time. The forward and backward regions, defined by the muon rapidity
$y^{*}_{\mu}$ in the centre-of-mass frame, are $1.53<y^{*}_{\mu}<4.03$ and
$-4.97<y^{*}_{\mu}<-2.47$, as the ${p}\mathrm{Pb}$ collision system
experiences an asymmetric distribution of beam energy.

In this article, the integrated and differential $Z$ boson production cross-
sections are measured at the Born level in QED, using $pp$ collision data
collected by the LHCb detector at a centre-of-mass energy of
$\sqrt{s}=5.02\text{\,Te\kern-1.00006ptV}$ in 2017, corresponding to an
integrated luminosity of $100\text{\,pb}^{-1}$ [43]. The production cross-
sections are measured in a fiducial region that closely matches the acceptance
of the LHCb detector, following the analysis strategy developed in Ref. [32].
The fiducial region is defined by requiring that both muons have a
pseudorapidity in the range of $2.0<\eta<4.5$ and transverse momentum
$p_{\mathrm{T}}>20\text{\,Ge\kern-1.00006ptV\\!/}c$, and that the di-muon
invariant mass is the interval
$60<M_{\mu\mu}<120\text{\,Ge\kern-1.00006ptV\\!/}c^{2}$. The differential
cross-section is measured as a function of the $Z$ boson rapidity ($y^{{Z}}$),
transverse momentum ($p_{\mathrm{T}}^{Z}$) and $\phi_{\eta}^{*}$ [44], which
is defined as

$\phi_{\eta}^{*}=\tan\left[(\pi-\Delta\phi^{\ell\ell})/2\right]\sin(\theta^{*}_{\eta}),$
(1)

where $\Delta\phi^{\ell\ell}$ represents the difference in the azimuthal angle
between the two muons in the laboratory frame. The variable
$\theta^{*}_{\eta}$ is defined by
$\cos(\theta^{*}_{\eta})=\tanh[(\eta^{-}-\eta^{+})/2]$, with $\eta^{-}$ and
$\eta^{+}$ denoting the pseudorapidities of the negatively and positively
charged muons in the laboratory frame. The observable $\phi_{\eta}^{*}$ probes
similar physics as the transverse momentum $p_{\mathrm{T}}^{Z}$, but is
measured with near-perfect resolution.

##  2 Detector and simulation

The LHCb detector [45, 46] is a single-arm forward spectrometer covering the
pseudorapidity range $2<\eta<5$, designed for the study of hadrons containing
$b$ or $c$ quarks. The detector includes a high-precision tracking system
consisting of a silicon-strip vertex detector surrounding the $pp$ interaction
region [47], a large-area silicon-strip detector (TT) [48], located upstream
of a dipole magnet with a bending power of about $4{\mathrm{\,Tm}}$, and three
stations of silicon-strip detectors and straw drift tubes [49] placed
downstream of the magnet. The tracking system provides a measurement of the
momentum, $p$, of charged particles with a relative uncertainty that varies
from 0.5 % at low momentum to 1.0 % at 200$\text{\,Ge\kern-1.00006ptV\\!/}c$.
The minimum distance of a track to a primary $pp$ collision vertex, the impact
parameter, is measured with a resolution of
$(15+29/p_{\mathrm{T}})\,\upmu\text{m}$, where $p_{\mathrm{T}}$ is the
component of the momentum transverse to the beam, in
$\text{\,Ge\kern-1.00006ptV\\!/}c$. Photons, electrons and hadrons are
identified by a calorimeter system consisting of scintillating-pad and
preshower detectors, an electromagnetic and a hadronic calorimeter. Muons are
identified by a system composed of alternating layers of iron and multiwire
proportional chambers [50]. The online event selection is performed by a
trigger [51], which consists of a hardware stage, based on information from
the calorimeter and muon systems, followed by a software stage, which applies
a full event reconstruction.

Simulation is used to model the effects of the detector acceptance and the
imposed selection requirements. In the simulation, $pp$ collisions are
generated using Pythia [52, 53] with a specific LHCb configuration [54]. Final
state radiation (FSR) is generated using Photos [55]. The interaction of the
generated particles with the detector, and its response, are implemented using
the Geant4 toolkit [56] as described in Ref. [57].

##  3 Event selection and background estimation

The muon triggers are responsible for the online event selection. At the
hardware trigger stage, a muon candidate with high $p_{\mathrm{T}}$ is
required. The muon candidate is required to have
$p_{\mathrm{T}}>6\text{\,Ge\kern-1.00006ptV\\!/}c$ and
$p>8\text{\,Ge\kern-1.00006ptV\\!/}c$, along with a good track fit quality in
the first software trigger stage. In the second software trigger stage, an
additional requirement of
$p_{\mathrm{T}}>12.5\text{\,Ge\kern-1.00006ptV\\!/}c$ is imposed on the muon
candidate. For a ${Z}\\!\rightarrow{\mu^{+}\mu^{-}}$ candidate, it is
necessary for at least one of the muons to pass both the hardware and software
trigger stages. These stages ensure the identification of the particle as a
muon with high transverse momentum.

A high purity ${Z}\\!\rightarrow{\mu^{+}\mu^{-}}$ sample is reconstructed from
a pair of opposite-signed tracks identified as muons. The invariant mass of
the di-muon is required to be within the range
$60<M_{\mu\mu}<120\text{\,Ge\kern-1.00006ptV\\!/}c^{2}$. For each muon track,
the fiducial requirements are
$p_{\mathrm{T}}>20\text{\,Ge\kern-1.00006ptV\\!/}c$ and pseudorapidity in the
range $2.0<\eta<4.5$. The muons are required to have momentum measurements
with relative uncertainties below 10%. In total, 3265
${Z}\\!\rightarrow{\mu^{+}\mu^{-}}$ candidates meet these selection criteria,
and the distribution of the di-muon invariant mass for the selected candidates
is shown in Fig. 1.

Figure 1: Mass distribution of the ${Z}\\!\rightarrow{\mu^{+}\mu^{-}}$ signal
candidates. The data are overlaid with model of the signal and background
models. The signal component is scaled such that the sum of the signal and
background matches the integral of the data.

The background contribution from decays of heavy flavour hadrons is estimated
using two control samples. Tracks from heavy flavour decays degrade the
primary vertex (PV) fit quality when included in the PV fit, since the heavy
hadron travels a certain distance before decaying. Hence, the first control
sample is obtained by applying a requirement that the selected candidate have
a PV with a low fit quality ($\chi^{2}>95$). Additionally, muons produced from
semileptonic decays of heavy flavour hadrons are less isolated. The variable
$I_{\mu}$ is defined as the ratio of the muon $p_{\mathrm{T}}$ to the vector
sum of the $p_{\mathrm{T}}$ of all charged particles in a cone around the
muon, where the cone is constrained using the radius parameter
$R=\sqrt{\Delta\eta^{2}+\Delta\phi^{2}}<0.5$. The second control sample is
selected by requiring that the two muons are not spatially isolated
($I_{\mu}<0.91$) from other rest of the event.

The yields of these two control samples are determined from a fit to the di-
muon invariant mass distributions, using an exponential function. The heavy
flavour background is concentrated at low mass
(${50<M_{\mu\mu}<70\text{\,Ge\kern-1.00006ptV\\!/}c^{2}}$), so to obtain a
larger sample and a more stable fit, the background yield is determined in the
region ${50<M_{\mu\mu}<110\text{\,Ge\kern-1.00006ptV\\!/}c^{2}}$. The
extrapolation corrects for the efficiency of the vertex quality and isolation
requirements. The efficiency of the muon isolation (vertex quality) selection
is calculated separately by applying the muon isolation requirement (vertex
quality requirement) to the first (second) control sample. The estimated
background contributions for these two subsamples are consistent, and the
average value is taken as the background contribution from the heavy flavour
decay process, which is determined to be $(1.98\pm 1.56)$% for the selected
${Z}\\!\rightarrow{\mu^{+}\mu^{-}}$ sample.

The contribution from the combinatorial background including misidentified
hadrons and $B$-$\bar{B}$ mixing is estimated using pairs of same-sign muons
in the data. In the same-sign events, a muon from heavy flavour decay combined
with a misidentified hadron are expected to make sizable contributions. After
removing the contribution from heavy flavour processes, the contribution from
misidentified hadrons is determined to be negligible. The electroweak
background contributions from the ${{t}{\overline{t}}}$, $W^{+}W^{-}$,
$W^{\pm}{Z}$, ${Z}{Z}$ and ${Z}\\!\rightarrow{\tau^{+}\tau^{-}}$ processes are
estimated from the simulation, with leading-order (LO) to next-to-next-to-
leading-order (NNLO) correction factors determined with the MCFM [58] package.

In summary, the total background contribution to the
${Z}\\!\rightarrow{\mu^{+}\mu^{-}}$ sample in the mass range
$60<M_{\mu\mu}<120\text{\,Ge\kern-1.00006ptV\\!/}c^{2}$ is determined to be
$(2.02\pm 1.56)\%$ and the background composition of the candidate sample is
summarised in Table 1. The data-driven backgrounds refer to processes
involving heavy flavor decays and hadron misidentification, which are
estimated based on the observed events in the data.

Table 1: Summary of the background composition in the ${Z}\\!\rightarrow{\mu^{+}\mu^{-}}$ data sample of candidates satisfying the signal selection.  Background | Estimation | Events | Fraction
---|---|---|---
Heavy flavour (${b}{\overline{b}}$, ${c}{\overline{c}}$) | data-driven | $65\pm 51$ | $2.0\times 10^{-2}$
Hadron misidentification | data-driven | $0.20\pm 0.10$ | $6.0\times 10^{-5}$
${Z}\\!\rightarrow{\tau^{+}\tau^{-}}$ | Simulation | $0.39\pm 0.12$ | $1.2\times 10^{-4}$
$WZ/ZZ$ | Simulation | $0.25\pm 0.04$ | $7.7\times 10^{-5}$
$WW$ | Simulation | $0.19\pm 0.02$ | $5.8\times 10^{-5}$
${t}{\overline{t}}$ | Simulation | $0.14\pm 0.10$ | $4.4\times 10^{-5}$
Total |  | $66\pm 51$ | $2.0\times 10^{-2}$

##  4 Cross-section determination

Only single-differential cross-section measurements are performed due to the
limited sample yields. The differential production cross-section is measured
as a function of $y^{{Z}}$, $p_{\mathrm{T}}^{Z}$ and $\phi_{\eta}^{*}$, and
the integrated cross-section is obtained by integrating over all bins, which
are chosen based on the detector resolution and sample size. Indicating with a
generic variable $a$, the differential cross-section is defined as

$\frac{d\sigma_{{Z}\\!\rightarrow{\mu^{+}\mu^{-}}}}{da}(i)=\frac{N_{Z}(i)\cdot
f_{FSR}^{Z}(i)}{\mathcal{L}\cdot\varepsilon^{Z}(i)\cdot\Delta a(i)},$ (2)

where the generic variable $a$ represents the observable $y^{{Z}}$,
$p_{\mathrm{T}}^{Z}$ or $\phi_{\eta}^{*}$, the index $i$ indicates the bin of
the variable under study, $N_{Z}(i)$ is the signal yield in bin $i$ after
background subtraction, $f_{FSR}^{Z}(i)$ is the FSR correction factor (as
discussed in Sec. 4.3), $\mathcal{L}$ is the integrated luminosity, $\Delta
a(i)$ is bin width for $i$-th bin (as presented in Tables 4 to 6 in Appendix
B), and $\varepsilon^{Z}(i)$ is the total efficiency in the $i$-th bin.

To account for detector misalignment effects, the $Z$ mass peak position and
resolution in simulated events are corrected to be compatible with the data,
using momentum scaling and smearing factors[59]. The impact from this
correction on the integrated cross-section measurement is found to be
negligible.

###  4.1 Efficiency

The selection efficiencies are determined for the muon tracking, muon
identification and trigger requirements. These are derived using the
${Z}\\!\rightarrow{\mu^{+}\mu^{-}}$ data and a tag-and-probe method [60].

In the determination of the tracking efficiency, a particle reconstructed in
all the tracking subdetectors, and fulfilling the muon trigger and muon
identification requirements, is used as the tag. An object reconstructed by
combining hits in the muon stations and the TT downstream tracking stations,
denoted as a MuonTT track, then acts as the probe. As described in Ref. [60],
the tracking efficiency is calculated as the fraction of probe candidates
matched with a reconstructed track. However, the precision of the measured
tracking efficiency is limited by the low number of $Z$-boson candidates in
the data sample. In this analysis, we therefore use the tracking efficiency
$\varepsilon^{\rm{MC},5.02}_{\rm{Tracking}}$ found from the 5.02 TeV
${Z}\\!\rightarrow{\mu^{+}\mu^{-}}$ simulation, and apply a correction to this
using a scale factor determined using the tag-and-probe method in the 13 TeV
analysis [32], as

$\varepsilon^{\rm{Data},5.02}_{\rm{Tracking}}=\varepsilon^{\rm{MC},5.02}_{\rm{Tracking}}\times\frac{\varepsilon^{\rm{\rm{Data},13}}_{\rm{Tracking}}}{\varepsilon^{\rm{MC},13}_{\rm{Tracking}}}.$
(3)

For the muon identification efficiency determination, the yield in the data
sample is also limited which leads to statistical fluctuations in the results.
Similar to the correction method employed for muon tracking efficiency, with
the data and simulated ID efficiency already available from the published
results at 13 TeV [32], a scale factor is applied to the 5.02 TeV simulation
efficiency.

The efficiency of the muon trigger is determined with the tag-and-probe
method, where a tag particle is chosen from a particle reconstructed in all
tracking subdetectors, which must be identified and triggered as a muon.
Meanwhile, the probe particle must pass all selection requirements used in the
analysis apart from the trigger requirements. The invariant mass of tag and
probe particles is further required to be within the range
[60,120]$\text{\,Ge\kern-1.00006ptV\\!/}c^{2}$, and the azimuthal separation,
$|\Delta\phi|$, greater than 2.7 radians. The efficiency is computed as the
ratio of the number of probes within the selected sample meeting the muon
trigger criteria to the total number of probes.

After correcting the muon tracking efficiency using Eq. 3 and applying this
correction method to the muon identification efficiency in a similar way, the
determined muon tracking reconstruction efficiency
$\varepsilon_{\rm{Track}}^{\mu^{\pm}}$, identification efficiency
$\varepsilon_{\rm{ID}}^{\mu^{\pm}}$ and trigger efficiency
$\varepsilon_{\rm{Trig}}^{\mu^{\pm}}$ vary between 94% and 98%, 90% and 97%,
62% and 85% respectively. These determined efficiencies as a function of muon
pseudorapidity are summarized in Appendix A.

This measurement follows the efficiency correction method employed in the 13
TeV analysis [32], which involved a one-dimensional efficiency correction
based solely on the muon pseudorapidity. The validity of this correction
method is ensured by conducting the so-called closure tests using simulated
samples. The differences between this correction method and the true
information from the simulation is considered as a source of systematic
uncertainty (as discussed in Sec. 5).

With the previous muon efficiencies $\varepsilon_{\rm{Track}}^{\mu^{\pm}}$,
$\varepsilon_{\rm{ID}}^{\mu^{\pm}}$, and $\varepsilon_{\rm{Trig}}^{\mu^{\pm}}$
determined, the total efficiency $\varepsilon^{Z}$ depends on the
pseudorapidities of the two final-state muons and can be calculated as

$\varepsilon^{Z}=\left(\varepsilon_{\rm{Track}}^{\mu^{+}}\cdot\varepsilon_{\rm{Track}}^{\mu^{-}}\right)\cdot\left(\varepsilon_{\rm{ID}}^{\mu^{+}}\cdot\varepsilon_{\rm{ID}}^{\mu^{-}}\right)\cdot\left(\varepsilon_{\rm{Trig}}^{\mu^{+}}+\varepsilon_{\rm{Trig}}^{\mu^{-}}-\varepsilon_{\rm{Trig}}^{\mu^{+}}\cdot\varepsilon_{\rm{Trig}}^{\mu^{-}}\right).$
(4)

###  4.2 Unfolding

The detector resolution causes migration between kinematic bins. Due to the
good angular resolution of the LHCb detector, negligible migration effects are
observed among $y^{{Z}}$ and $\phi_{\eta}^{*}$ bins.

To assess the necessity of using an unfolding method to correct the
$p_{\mathrm{T}}^{Z}$ measurement at the reconstruction level, the ratio of the
$p_{\mathrm{T}}^{Z}$ distribution before and after unfolding is calculated
from data, and then compared to the ratio of the reconstructed
$p_{\mathrm{T}}^{Z}$ distribution and the generated $p_{\mathrm{T}}^{Z}$
distribution, obtained from the simulated sample. Since notable migration is
observed, the unfolding procedure is performed for the $p_{\mathrm{T}}^{Z}$
distribution, using the Bayesian unfolding method [61, 62].

###  4.3 Final state radiation correction

The measured cross-section is corrected to the Born level in QED, so that it
can be directly compared with theoretical predictions. The final-state
radiation (FSR) correction is developed and applied to the measurements, by
comparing the ResBos [63] predictions with and without the implementation of
Photos [55], which corrects the quantities of muons after final state
radiation to the Born level. The FSR corrections in bins of $y^{{Z}}$,
$p_{\mathrm{T}}^{Z}$ and $\phi_{\eta}^{*}$ are shown in Fig. 2. The
corrections for single-differential cross-section measurements are also
presented in Appendix B.

Figure 2: Final state radiation correction estimated for the (top-left)
$y^{{Z}}$, (top-right) $p_{\mathrm{T}}^{Z}$, and (bottom) $\phi_{\eta}^{*}$
differential cross-section measurements. The error bars represent the total
(statistical and systematic) uncertainties.

##  5 Systematic uncertainties

To account for potential sources of systematic uncertainties in the cross-
section measurements, several sources are considered. These uncertainties
include the background estimation, the calibration of the momentum scale, the
efficiency determination, the unfolding procedure, the results of efficiency
closure tests, the FSR correction, and the measurement of the integrated
luminosity.

The determination of the heavy flavour background uses the averaged yield
between two methods as the background contribution. The associated uncertainty
is calculated as the discrepancy between the background yields obtained from
the two control samples. Additionally, a systematic uncertainty is introduced
by varying the mass region and selection requirements of the control samples.
The hadron misidentification and other backgrounds are estimated from the data
and simulation, with the statistical uncertainties being treated as systematic
uncertainties that depend on the limited size of the same-sign data and
simulation samples. The systematic uncertainties on the ${t}{\overline{t}}$,
$W^{+}W^{-}$, $W^{\pm}{Z}$, ${Z}{Z}$, and
${Z}\\!\rightarrow{\tau^{+}\tau^{-}}$ components derive from the statistical
uncertainties and theoretical uncertainties on the LO to NNLO correction
factors. To estimate the uncertainty due to inadequate calibration of the
detector, the momentum scaling and smearing correction for the simulation are
studied in order to improve the modeling of the
${Z}\\!\rightarrow{\mu^{+}\mu^{-}}$ data. Only a small fraction, less than
0.01% of events, exhibits changes when comparing the simulation before and
after studying the momentum scaling and smearing. This value is conservatively
assigned as the systematic uncertainty related to the detector alignment. The
impact from momentum calibration uncertainty on the integrated cross-section
measurement is found to be negligible. However, the alignment uncertainty is
considered and added for the $p_{\mathrm{T}}^{Z}$ distribution.

The track reconstruction and identification efficiencies for high
$p_{\mathrm{T}}$ muons are calculated using simulation and data taking at
$\sqrt{s}=13\text{\,Te\kern-1.00006ptV}$. The systematic uncertainties are
determined by the size of the 5.02 TeV simulation samples and the
uncertainties are propagated from the 13 TeV results. For the trigger
efficiency, which is directly measured using control samples in the 5.02 TeV
data, a systematic uncertainty is assigned for variations due to the limited
size of the samples studied. As studied in Ref. [32], an additional systematic
uncertainty is evaluated, based on the method used to determine the tracking
efficiency. This is evaluated to be 0.47% which is already considered in the
systematic uncertainty of the muon tracking efficiency.

To examine whether it is sufficient to perform efficiency corrections only as
a function of the muon pseudorapidity variable or if the efficiency also
depends on other variables such as muon $p_{\mathrm{T}}$, a closure test is
performed. The reconstructed event yields in simulation are adjusted using the
efficiencies determined solely based on the muon pseudorapidity obtained from
the simulation samples and then compared with the generated yield. The
differences, which do not exhibit any systematic pattern across the
differential cross-section measurement regions, are attributed as an
additional source of uncertainty.

To estimate the uncertainty attributed to the unfolding process, the
$p_{\mathrm{T}}^{Z}$ distribution is unfolded using the so-called Invert
Approach [64], which employs a simple inversion of the response matrix without
regularisation. The deviation of the results from the Bayesian Approach [61,
62] is taken as a systematic uncertainty on the $p_{\mathrm{T}}^{Z}$
differential cross-section. The systematic uncertainty arising from the FSR
correction is evaluated by comparing the default correction with one
determined using the Powheg generator showered using Pythia. The differences
in FSR corrections between ResBos with Photos and Powheg with Pythia are then
taken into consideration as a systematic uncertainty. Regarding the data
sample used, the luminosity is determined with a precision of 2.0% [43], which
is quoted separately to the other sources of systematic uncertainty. The
statistical and systematic uncertainties in the integrated cross-section
measurement are listed in Table 2. The different sources of systematic
uncertainty for each bin of the differential cross-sections are summarised in
Tables 13 to 15 in Appendix D.

Table 2: The uncertainties for the integrated ${Z}\\!\rightarrow{\mu^{+}\mu^{-}}$ cross-section measurement. Source |  $\Delta\sigma\,$[pb] |  $\Delta\sigma/\sigma$ [%]
---|---|---
Luminosity | 0.79 | 2.00
Statistical | 0.70 | 1.77
Tracking | 0.40 | 1.01
Efficiency Closure | 0.24 | 0.61
Trigger | 0.21 | 0.54
Background | 0.19 | 0.48
Identification | 0.10 | 0.25
FSR | 0.07 | 0.18
Calibration | $<4.0\times 10^{-3}$ | $<0.01$
Total Systematic (excl. lumi.) | 0.56 | 1.42

##  6 Results

###  6.1 Differential cross-section results

The single differential cross-sections in regions of $y^{{Z}}$,
$p_{\mathrm{T}}^{Z}$ and $\phi_{\eta}^{*}$ are shown in Fig. 3, together with
the ratios of theoretical predictions to data. The numerical results of single
differential cross-section within each bin are summarised in Tables 16–18 in
Appendix E.

Figure 3: (Left) Measured single differential cross-section as a function of
$y^{{Z}}$, $p_{\mathrm{T}}^{Z}$ and $\phi_{\eta}^{*}$ compared with different
theoretical predictions. (Right) Ratio of theoretical predictions to measured
values, with the horizontal bars showing the uncertainty from the PDFs. The
green band, centered at unity, shows the uncertainty of the measurement.

The ResBos [63] program performs resummation of large logarithms, up to next-
to-next-leading-logarithm accuracy, within the Collins-Soper-Sterman
resummation formalism [65, 66, 67] and matches to next-to-leading-order (NLO)
fixed order calculations. The enhanced QCD prediction is obtained from the
simulated ${Z}\\!\rightarrow{\mu^{+}\mu^{-}}$ sample, generated with ResBos
[63] using the CT18 PDFs [36] for all measurements. Powheg-BOX [68, 69, 70,
71] is an NLO generator, and can be interfaced with Pythia for QCD and EW
showering. Although Pythia is a LO generator, it approximates higher order
effects in initial and final states via a parton showering algorithm [72]. In
this analysis, Powheg-BOX is used to generate
${Z}\\!\rightarrow{\mu^{+}\mu^{-}}$ events, followed by hadronization using
Pythia for all measurements. The MCFM package [58], a fixed-order NNLO
generator, is used here to estimate the $Z$ boson production cross-section as
a function of $y^{{Z}}$ in the acceptance of the LHCb detector.

For these three generators, different PDF sets NNPDF3.1 [37], NNPDF4.0 [38],
MSHT20 [39], and CT18 [36] are employed to provide theoretical predictions.
Measurements are in reasonable agreement with the different theoretical
predictions. In both the lower and higher $p_{\mathrm{T}}^{Z}$ and
$\phi_{\eta}^{*}$ region, the measurements agree with predictions from both
ResBos and Powheg. For the measurement of $y^{{Z}}$, the theoretical
predictions from the three generators are also compatible with the
experimental results within the uncertainty range.

###  6.2 Correlation matrices

The event migration between bins causes statistical correlations, which are
determined using simulation. Large correlations are found in the
low-$p_{\mathrm{T}}^{Z}$ region and small correlations in the
high-$p_{\mathrm{T}}^{Z}$ region, while the statistical correlations are
negligible for both the $y^{{Z}}$ and $\phi_{\eta}^{*}$ distributions.

For the differential cross-section measurements, background, alignment,
efficiency closure test, and FSR uncertainties are assumed to be 50%
correlated between different bins, while the luminosity uncertainty is
considered to be 100% correlated. The calculated correlation matrices for the
efficiencies are presented in Appendix C. Large correlations between different
bins are present in the $p_{\mathrm{T}}^{Z}$ differential cross-section
measurement, but small correlations are also present between most bins for the
$y^{{Z}}$ and $\phi_{\eta}^{*}$ measurements.

###  6.3 Integrated cross-section results

Using the LHCb 2017 $pp$ collision data at
$\sqrt{s}=5.02\text{\,Te\kern-1.00006ptV}$, the integrated Born-level $Z$
boson production cross-section, with two muons in the final state and within
the LHCb acceptance is

$\sigma_{{Z}\\!\rightarrow{\mu^{+}\mu^{-}}}=39.6\pm 0.7\text{\,(stat)}\pm
0.6\text{\,(syst)}\pm 0.8\text{\,(lumi)}\,\text{\,pb}$ (5)

where the uncertainties are due to statistical effects, systematic effects,
and the luminosity measurement, respectively.

In this article, the $Z$ boson is defined to also include contributions from
virtual photons, and the interference between them since they cannot be
distinguished experimentally. The measured results obtained for the total
cross-section of the $pp\rightarrow Z\rightarrow{\mu^{+}}{\mu^{-}}$ process at
5.02 TeV in the fiducial region of the LHCb detector have been compared to the
predictions obtained using MCFM with CT18NNLO and other theoretical models,
including Powheg-Box with NNPDF3.1, NNPDF4.0, MSHT20, CT18, and ResBos with
CT18, all of which account for both statistical and PDF uncertainties. These
theoretical predictions are compared to the results in Fig. 4, which
demonstrate a reasonable agreement.

By comparing the theoretical predictions at different center-of-mass energies
provided by MCFM, with the experimental measurements previously obtained by
LHCb at $\sqrt{s}$ = 7, 8 and $13\text{\,Te\kern-1.00006ptV}$, as illustrated
in Fig. 5, a good level of consistency can be observed overall.

Figure 4: Comparison of the integrated cross-section,
$\sigma_{{Z}\\!\rightarrow{\mu^{+}\mu^{-}}}$, between data and theoretical
predictions. The bands correspond to the data, with the inner band
corresponding to the statistical uncertainty and the outer bands corresponding
to the systematic uncertainty and total uncertainty.  Figure 5: Measured
$\sigma_{{Z}\\!\rightarrow{\mu^{+}\mu^{-}}}$ for $pp$ collisions, as a
function of $\sqrt{s}$. The uncertainty on the data is too small to be
visible. The data are overlaid with a curve showing the theoretical
prediction.

###  6.4 Nuclear modification factors

The nuclear modification factors are determined based on the present
measurements and those in Ref. [42]. Here, the statistical and systematic
uncertainties (including uncertainties due to integrated luminosity) of the
cross-sections measured between $pp$ and $p$Pb collisions are assumed to be
fully uncorrelated when being propagated to the nuclear modification factors.
Efficiency and background modeling uncertainties are also treated as fully
uncorrelated, owing to the very different running conditions between $pp$ and
$p$Pb collisions such as the heavier detector contamination from the much
higher charged hadron multiplicities in case of $p$Pb. All other sources of
systematic uncertainty are considered to be uncorrelated, an assumption that
has negligible effect to the reported results.

The unique forward geometry coverage allows the LHCb detector to probe nPDFs
at very small Bjorken-$x$ ($10^{-4}<x<10^{-3}$). The nuclear modification
factors can be calculated using the $pp$ cross-section measured here together
with the $Z$ boson production cross-section in $p$Pb collisions measured in
Ref. [42]. The nuclear modification factors are defined as the ratio of the
cross-sections between ${p}\mathrm{Pb}$ and $pp$ collisions under the same
muon rapidity acceptance. Due to the asymmetric beam energy in the
${p}\mathrm{Pb}$ collision system, the muon rapidity acceptance
($2.0<\eta<4.5$) becomes $1.53<y^{*}_{\mu}<4.03$ in case of the forward
collisions, and $-4.97<y^{*}_{\mu}<-2.47$ in case of backward collisions,
where $y^{*}_{\mu}$ represents the muon rapidity in the centre-of-mass frame.
The two rapidity quantities are identical in $pp$ collisions,
$\eta_{\mu}\equiv y^{*}_{\mu}$, due to the symmetric beam energy. Forward
collisions refer to the proton beam entering the LHCb detector along the
positive direction of the $z$ axis, while backward collisions correspond to
the proton beam going in the opposite direction.

The cross-sections of the ${p}\mathrm{Pb}$ collisions measured under different
rapidity acceptances need to be corrected before being used to calculate the
nuclear modification factors. The $k_{{p}\mathrm{Pb}}$ factor corrects for the
different $\eta$ acceptance between ${p}\mathrm{Pb}$ and $pp$ collisions and
can be calculated using Powheg with the proton PDF set CTEQ6.1 [73],

$k_{{p}\mathrm{Pb}}^{\rm
F}=\dfrac{\sigma_{(pp,~{}2.0<\eta<4.5)}}{\sigma_{(pp,~{}1.53<\eta<4.03)}}=0.706\pm
0.002,$ (6)

and

$k_{{p}\mathrm{Pb}}^{\rm
B}=\dfrac{\sigma_{(pp,~{}2.0<\eta<4.5)}}{\sigma_{(pp,~{}-4.97<\eta<-2.47)}}=1.518\pm
0.003,$ (7)

for forward and backward collisions, respectively. The integrated nuclear
modification factors can then be defined as follows

$R_{{p}\mathrm{Pb}}^{\rm F}=k_{{p}\mathrm{Pb}}^{\rm
F}\cdot\dfrac{\sigma_{({p}\mathrm{Pb},~{}1.53<y^{*}_{\mu}<4.03)}}{208\cdot\sigma_{(pp,~{}2.0<\eta<4.5)}},$
(8)

and

$R_{{p}\mathrm{Pb}}^{\rm B}=k_{{p}\mathrm{Pb}}^{\rm
B}\cdot\dfrac{\sigma_{({p}\mathrm{Pb},~{}-4.97<y^{*}_{\mu}<-2.47)}}{208\cdot\sigma_{(pp,~{}2.0<\eta<4.5)}},$
(9)

for forward and backwards ${p}\mathrm{Pb}$ collisions, respectively, with 208
being the number of binary nucleon-nucleon collisions in ${p}\mathrm{Pb}$
collisions.

The nuclear modification factors are calculated to be

$R_{pPb}^{\rm F}=1.2^{+0.5}_{-0.3}\text{\,(stat)}\pm 0.1\text{\,(syst)}$

in the forward collisions, and

$R_{pPb}^{\rm B}=3.6^{+1.6}_{-0.9}\text{\,(stat)}\pm 0.2\text{\,(syst)}$

in the backward collisions. The large statistical uncertainties are due to the
small size of the $p$Pb data sample.

Based on the theoretical prediction of the cross-sections from Ref. [42]
derived using FEWZ [74] at NNLO with EPS09 nPDFs [40] and MSTW08 PDF [41], the
predicted $R_{{p}\mathrm{Pb}}$ can be derived as $R_{{p}\mathrm{Pb}}^{\rm
F,\,\rm theo.}=0.906^{+0.002}_{-0.007}$ and $R_{pPb}^{\rm B,\,\rm
theo.}=0.929^{+0.011}_{-0.028}$ for the forward and backward collisions,
respectively. The measurement and theoretical prediction agree within the
uncertainties in the forward region, whereas there is some tension between the
results in the backward region. However, due to the limited size of the
${p}\mathrm{Pb}$ sample, the tension could be caused by a statistical
fluctuation.

##  7 Conclusion

This paper reports the first measurement of $Z$ boson production cross-section
at the centre-of-mass energy $\sqrt{s}=5.02\text{\,Te\kern-1.00006ptV}$, using
the LHCb 2017 $pp$ collision dataset. The techniques employed in this analysis
closely follow the techniques and strategies established in a previous LHCb
Run 2 analysis [32]. The results show reasonable agreement with different
Standard Model predictions. Combining with the first inclusive $Z$ production
result in $p$Pb collisions at $\sqrt{s_{NN}}=5.02\text{\,Te\kern-1.00006ptV}$
[42] the nuclear modification factors in the forward and backward regions have
been measured for the first time.

## Appendices

##  Appendix A Efficiency

Results of muon efficiencies used in the total $Z$-boson efficiency
calculation are presented in Fig. 6 and Table 3.

Figure 6: Muon tracking, identification and trigger efficiency as a function of pseudorapidity, estimated from the data and simulation at 5.02 TeV and 13 TeV.  Table 3: The muon tracking, muon identification and trigger efficiency in each pseudorapidity bin. $\eta$ | $\varepsilon^{\mu}_{\rm{Track}}\,[\%]$ | $\varepsilon^{\mu}_{\rm{ID}}\,[\%]$ | $\varepsilon^{\mu}_{\rm{Trig}}\,[\%]$
---|---|---|---
$[\,$2.00, 2.25$\,]$ | 93.77 | 96.26 | 84.23
$[\,$2.25, 2.50$\,]$ | 96.82 | 96.30 | 83.95
$[\,$2.50, 2.75$\,]$ | 97.51 | 95.85 | 85.20
$[\,$2.75, 3.00$\,]$ | 95.14 | 96.55 | 80.08
$[\,$3.00, 3.25$\,]$ | 96.26 | 96.73 | 79.41
$[\,$3.25, 3.50$\,]$ | 98.32 | 96.92 | 81.50
$[\,$3.50, 3.75$\,]$ | 97.64 | 95.02 | 75.90
$[\,$3.75, 4.00$\,]$ | 96.00 | 97.08 | 69.23
$[\,$4.00, 4.25$\,]$ | 96.41 | 97.28 | 69.29
$[\,$4.25, 4.50$\,]$ | 93.86 | 90.10 | 61.54

##  Appendix B Final state radiation corrections

Tabulated results of final state radiation corrections used in the single
differential cross-section measurements are presented in Tables 4 to 6.

Table 4: Final state radiation correction used in the $y^{{Z}}$ cross-section measurement. The first uncertainty is statistical and the second is systematic. $y^{{Z}}$ | Correction
---|---
$[$ 2.000, 2.125 $]$ | 1.018 $\pm$ 0.004 $\pm$ 0.092
$[$ 2.125, 2.250 $]$ | 1.018 $\pm$ 0.002 $\pm$ 0.000
$[$ 2.250, 2.375 $]$ | 1.019 $\pm$ 0.002 $\pm$ 0.006
$[$ 2.375, 2.500 $]$ | 1.021 $\pm$ 0.002 $\pm$ 0.057
$[$ 2.500, 2.625 $]$ | 1.021 $\pm$ 0.001 $\pm$ 0.046
$[$ 2.625, 2.750 $]$ | 1.022 $\pm$ 0.001 $\pm$ 0.122
$[$ 2.750, 2.875 $]$ | 1.024 $\pm$ 0.002 $\pm$ 0.009
$[$ 2.875, 3.000 $]$ | 1.026 $\pm$ 0.002 $\pm$ 0.162
$[$ 3.000, 3.125 $]$ | 1.026 $\pm$ 0.002 $\pm$ 0.070
$[$ 3.125, 3.250 $]$ | 1.027 $\pm$ 0.002 $\pm$ 0.128
$[$ 3.250, 3.375 $]$ | 1.027 $\pm$ 0.002 $\pm$ 0.011
$[$ 3.375, 3.625 $]$ | 1.025 $\pm$ 0.002 $\pm$ 0.023
$[$ 3.625, 4.000 $]$ | 1.020 $\pm$ 0.004 $\pm$ 0.001
Table 5: Final state radiation correction used in the $p_{\mathrm{T}}^{Z}$ cross-section measurement. The first uncertainty is statistical and the second is systematic. $p_{\mathrm{T}}^{Z}$ [$\text{\,Ge\kern-1.00006ptV\\!/}c$ ] | Correction
---|---
$[$ 0.0, 2.2$\,]$ | 1.092 $\pm$ 0.002 $\pm$ 0.020
$[$ 2.2, 3.4$\,]$ | 1.080 $\pm$ 0.002 $\pm$ 0.018
$[$ 3.4, 4.6$\,]$ | 1.063 $\pm$ 0.002 $\pm$ 0.003
$[$ 4.6, 5.8$\,]$ | 1.044 $\pm$ 0.002 $\pm$ 0.015
$[$ 5.8, 7.2$\,]$ | 1.027 $\pm$ 0.002 $\pm$ 0.029
$[$ 7.2, 8.7$\,]$ | 1.012 $\pm$ 0.002 $\pm$ 0.025
$[$ 8.7, 10.5$\,]$ | 1.001 $\pm$ 0.002 $\pm$ 0.003
$[$ 10.5, 12.8$\,]$ | 0.987 $\pm$ 0.002 $\pm$ 0.006
$[$ 12.8, 15.4$\,]$ | 0.977 $\pm$ 0.002 $\pm$ 0.010
$[$ 15.4, 19.0$\,]$ | 0.968 $\pm$ 0.002 $\pm$ 0.005
$[$ 19.0, 34.0$\,]$ | 0.985 $\pm$ 0.001 $\pm$ 0.000
$[$ 34.0, 120.0$\,]$ | 1.038 $\pm$ 0.002 $\pm$ 0.000
Table 6: Final state radiation correction used in the $\phi_{\eta}^{*}$ cross-section measurement. The first uncertainty is statistical and the second is systematic. $\phi_{\eta}^{*}$ | Correction
---|---
$[$ 0.00, 0.01 $]$ | 1.035 $\pm$ 0.002 $\pm$ 3.802
$[$ 0.01, 0.02 $]$ | 1.034 $\pm$ 0.002 $\pm$ 0.919
$[$ 0.02, 0.03 $]$ | 1.032 $\pm$ 0.002 $\pm$ 1.384
$[$ 0.03, 0.05 $]$ | 1.027 $\pm$ 0.001 $\pm$ 0.878
$[$ 0.05, 0.07 $]$ | 1.021 $\pm$ 0.002 $\pm$ 1.060
$[$ 0.07, 0.10 $]$ | 1.017 $\pm$ 0.002 $\pm$ 0.378
$[$ 0.10, 0.15 $]$ | 1.013 $\pm$ 0.002 $\pm$ 0.220
$[$ 0.15, 0.20 $]$ | 1.011 $\pm$ 0.002 $\pm$ 0.288
$[$ 0.20, 0.30 $]$ | 1.010 $\pm$ 0.002 $\pm$ 0.236
$[$ 0.30, 0.60 $]$ | 1.011 $\pm$ 0.002 $\pm$ 0.039
$[$ 0.60, 1.20 $]$ | 1.017 $\pm$ 0.006 $\pm$ 0.001

##  Appendix C Correlation matrices

The calculated statistical correlation matrices are shown in Fig. 7 and
presented in Tables 7 to 9. The correlation matrices for the efficiency
uncertainty are shown in Fig. 8 for single differential cross-section
measurements, and presented in Tables 10 to 12.

Figure 7: Statistical correlation matrices for the differential cross-section
measurements as functions of (top-left) $y^{{Z}}$, (top-right)
$p_{\mathrm{T}}^{Z}$ and (bottom) $\phi_{\eta}^{*}$.

Table 7: Statistical correlation matrix for the one-dimensional $y^{{Z}}$ measurement. Bin | 1 | 2 | 3 | 4 | 5 | 6 | 7 | 8 | 9 | 10 | 11 | 12 | 13
---|---|---|---|---|---|---|---|---|---|---|---|---|---
1 | 1.00 |  |  |  |  |  |  |  |  |  |  |  |
2 | 0.00 | 1.00 |  |  |  |  |  |  |  |  |  |  |
3 | 0.00 | 0.00 | 1.00 |  |  |  |  |  |  |  |  |  |
4 | 0.00 | 0.00 | 0.00 | 1.00 |  |  |  |  |  |  |  |  |
5 | 0.00 | 0.00 | 0.00 | 0.01 | 1.00 |  |  |  |  |  |  |  |
6 | 0.00 | 0.00 | 0.00 | 0.00 | 0.01 | 1.00 |  |  |  |  |  |  |
7 | 0.00 | 0.00 | 0.00 | 0.00 | 0.00 | 0.01 | 1.00 |  |  |  |  |  |
8 | 0.00 | 0.00 | 0.00 | 0.00 | 0.00 | 0.00 | 0.01 | 1.00 |  |  |  |  |
9 | 0.00 | 0.00 | 0.00 | 0.00 | 0.00 | 0.00 | 0.00 | 0.01 | 1.00 |  |  |  |
10 | 0.00 | 0.00 | 0.00 | 0.00 | 0.00 | 0.00 | 0.00 | 0.00 | 0.01 | 1.00 |  |  |
11 | 0.00 | 0.00 | 0.00 | 0.00 | 0.00 | 0.00 | 0.00 | 0.00 | 0.00 | 0.01 | 1.00 |  |
12 | 0.00 | 0.00 | 0.00 | 0.00 | 0.00 | 0.00 | 0.00 | 0.00 | 0.00 | 0.00 | 0.02 | 1.00 |
13 | 0.00 | 0.00 | 0.00 | 0.00 | 0.00 | 0.00 | 0.00 | 0.00 | 0.00 | 0.00 | 0.00 | 0.00 | 1.00

Table 8: Statistical correlation matrix for the one-dimensional $p_{\mathrm{T}}^{Z}$ measurement. Bin | 1 | 2 | 3 | 4 | 5 | 6 | 7 | 8 | 9 | 10 | 11 | 12
---|---|---|---|---|---|---|---|---|---|---|---|---
1 | 1.00 |  |  |  |  |  |  |  |  |  |  |
2 | 0.16 | 1.00 |  |  |  |  |  |  |  |  |  |
3 | 0.00 | 0.18 | 1.00 |  |  |  |  |  |  |  |  |
4 | 0.00 | 0.01 | 0.18 | 1.00 |  |  |  |  |  |  |  |
5 | 0.00 | 0.00 | 0.01 | 0.18 | 1.00 |  |  |  |  |  |  |
6 | 0.00 | 0.00 | 0.00 | 0.01 | 0.14 | 1.00 |  |  |  |  |  |
7 | 0.00 | 0.00 | 0.00 | 0.00 | 0.01 | 0.12 | 1.00 |  |  |  |  |
8 | 0.00 | 0.00 | 0.00 | 0.00 | 0.00 | 0.00 | 0.10 | 1.00 |  |  |  |
9 | 0.00 | 0.00 | 0.00 | 0.00 | 0.00 | 0.00 | 0.00 | 0.10 | 1.00 |  |  |
10 | 0.00 | 0.00 | 0.00 | 0.00 | 0.00 | 0.00 | 0.00 | 0.00 | 0.05 | 1.00 |  |
11 | 0.00 | 0.00 | 0.00 | 0.00 | 0.00 | 0.00 | 0.00 | 0.00 | 0.00 | 0.04 | 1.00 |
12 | 0.00 | 0.00 | 0.00 | 0.00 | 0.00 | 0.01 | 0.00 | 0.00 | 0.00 | 0.01 | 0.01 | 1.00

Table 9: Statistical correlation matrix for the one-dimensional $\phi_{\eta}^{*}$ measurement. Bin | 1 | 2 | 3 | 4 | 5 | 6 | 7 | 8 | 9 | 10 | 11
---|---|---|---|---|---|---|---|---|---|---|---
1 | 1.00 |  |  |  |  |  |  |  |  |  |
2 | 0.01 | 1.00 |  |  |  |  |  |  |  |  |
3 | 0.00 | 0.01 | 1.00 |  |  |  |  |  |  |  |
4 | 0.00 | 0.00 | 0.01 | 1.00 |  |  |  |  |  |  |
5 | 0.00 | 0.00 | 0.00 | 0.01 | 1.00 |  |  |  |  |  |
6 | 0.00 | 0.00 | 0.00 | 0.00 | 0.00 | 1.00 |  |  |  |  |
7 | 0.00 | 0.00 | 0.00 | 0.00 | 0.00 | 0.00 | 1.00 |  |  |  |
8 | 0.00 | 0.00 | 0.00 | 0.00 | 0.00 | 0.00 | 0.00 | 1.00 |  |  |
9 | 0.00 | 0.00 | 0.00 | 0.00 | 0.00 | 0.00 | 0.00 | 0.00 | 1.00 |  |
10 | 0.00 | 0.00 | 0.00 | 0.00 | 0.00 | 0.00 | 0.00 | 0.00 | 0.00 | 1.00 |
11 | 0.00 | 0.00 | 0.00 | 0.00 | 0.00 | 0.00 | 0.00 | 0.00 | 0.00 | 0.00 | 1.00

Figure 8: Correlation matrices for the efficiency uncertainties as functions
of (top-left) $y^{{Z}}$, (top-right) $p_{\mathrm{T}}^{Z}$ and (bottom)
$\phi_{\eta}^{*}$.

Table 10: Correlation matrix for the efficiency uncertainty of the one-dimensional $y^{{Z}}$ measurement. Bin | 1 | 2 | 3 | 4 | 5 | 6 | 7 | 8 | 9 | 10 | 11 | 12 | 13
---|---|---|---|---|---|---|---|---|---|---|---|---|---
1 | 1.00 |  |  |  |  |  |  |  |  |  |  |  |
2 | 0.92 | 1.00 |  |  |  |  |  |  |  |  |  |  |
3 | 0.68 | 0.89 | 1.00 |  |  |  |  |  |  |  |  |  |
4 | 0.56 | 0.74 | 0.93 | 1.00 |  |  |  |  |  |  |  |  |
5 | 0.49 | 0.62 | 0.77 | 0.92 | 1.00 |  |  |  |  |  |  |  |
6 | 0.34 | 0.47 | 0.65 | 0.81 | 0.95 | 1.00 |  |  |  |  |  |  |
7 | 0.31 | 0.41 | 0.54 | 0.67 | 0.82 | 0.94 | 1.00 |  |  |  |  |  |
8 | 0.19 | 0.28 | 0.40 | 0.53 | 0.69 | 0.81 | 0.93 | 1.00 |  |  |  |  |
9 | 0.14 | 0.21 | 0.29 | 0.38 | 0.51 | 0.63 | 0.78 | 0.94 | 1.00 |  |  |  |
10 | 0.13 | 0.19 | 0.27 | 0.35 | 0.44 | 0.54 | 0.69 | 0.87 | 0.99 | 1.00 |  |  |
11 | 0.13 | 0.18 | 0.25 | 0.33 | 0.42 | 0.52 | 0.68 | 0.87 | 0.98 | 1.00 | 1.00 |  |
12 | 0.16 | 0.22 | 0.29 | 0.37 | 0.46 | 0.57 | 0.72 | 0.89 | 0.99 | 1.00 | 1.00 | 1.00 |
13 | 0.16 | 0.21 | 0.27 | 0.34 | 0.42 | 0.52 | 0.68 | 0.87 | 0.98 | 1.00 | 1.00 | 1.00 | 1.00

Table 11: Correlation matrix for the efficiency uncertainty of the one-dimensional $p_{\mathrm{T}}^{Z}$ measurement. Bin | 1 | 2 | 3 | 4 | 5 | 6 | 7 | 8 | 9 | 10 | 11 | 12
---|---|---|---|---|---|---|---|---|---|---|---|---
1 | 1.00 |  |  |  |  |  |  |  |  |  |  |
2 | 1.00 | 1.00 |  |  |  |  |  |  |  |  |  |
3 | 1.00 | 0.99 | 1.00 |  |  |  |  |  |  |  |  |
4 | 1.00 | 1.00 | 1.00 | 1.00 |  |  |  |  |  |  |  |
5 | 1.00 | 1.00 | 1.00 | 1.00 | 1.00 |  |  |  |  |  |  |
6 | 0.99 | 0.98 | 0.99 | 0.99 | 0.98 | 1.00 |  |  |  |  |  |
7 | 1.00 | 1.00 | 1.00 | 1.00 | 1.00 | 0.98 | 1.00 |  |  |  |  |
8 | 0.99 | 0.99 | 0.99 | 1.00 | 0.99 | 0.99 | 0.99 | 1.00 |  |  |  |
9 | 0.99 | 0.99 | 0.98 | 0.99 | 0.99 | 0.98 | 0.98 | 0.99 | 1.00 |  |  |
10 | 0.98 | 0.98 | 0.98 | 0.99 | 0.98 | 0.99 | 0.97 | 0.99 | 0.98 | 1.00 |  |
11 | 0.96 | 0.94 | 0.95 | 0.96 | 0.95 | 0.97 | 0.93 | 0.97 | 0.97 | 0.98 | 1.00 |
12 | 0.94 | 0.93 | 0.93 | 0.95 | 0.93 | 0.96 | 0.92 | 0.97 | 0.97 | 0.98 | 0.99 | 1.00

Table 12: Correlation matrix for the efficiency uncertainty of the one-dimensional $\phi_{\eta}^{*}$ measurement. Bin | 1 | 2 | 3 | 4 | 5 | 6 | 7 | 8 | 9 | 10 | 11
---|---|---|---|---|---|---|---|---|---|---|---
1 | 1.00 |  |  |  |  |  |  |  |  |  |
2 | 0.99 | 1.00 |  |  |  |  |  |  |  |  |
3 | 0.99 | 1.00 | 1.00 |  |  |  |  |  |  |  |
4 | 1.00 | 1.00 | 1.00 | 1.00 |  |  |  |  |  |  |
5 | 1.00 | 0.99 | 0.99 | 1.00 | 1.00 |  |  |  |  |  |
6 | 0.99 | 0.98 | 0.99 | 0.99 | 0.99 | 1.00 |  |  |  |  |
7 | 0.99 | 0.98 | 0.99 | 0.99 | 0.99 | 0.99 | 1.00 |  |  |  |
8 | 0.99 | 0.98 | 0.98 | 0.99 | 0.99 | 1.00 | 0.98 | 1.00 |  |  |
9 | 0.97 | 0.96 | 0.97 | 0.97 | 0.98 | 0.99 | 0.98 | 0.98 | 1.00 |  |
10 | 0.96 | 0.95 | 0.96 | 0.96 | 0.97 | 0.99 | 0.98 | 0.98 | 0.99 | 1.00 |
11 | 0.94 | 0.95 | 0.96 | 0.95 | 0.95 | 0.96 | 0.95 | 0.95 | 0.95 | 0.97 | 1.00

##  Appendix D Summary of systematic uncertainties

The summarised systematic uncertainties for single differential cross-sections
are shown in Tables 13 to 15.

Table 13: Systematic uncertainties on the single differential cross-sections in bins of $y^{{Z}}$, presented in percentage.  $y^{{Z}}$ | Efficiency | Background | FSR | Closure
---|---|---|---|---
$[\,$2.000, 2.125$\,]$ | 1.42 | 0.04 | 2.10 | 0.45
$[\,$2.125, 2.250$\,]$ | 1.37 | 0.77 | 0.00 | 0.31
$[\,$2.250, 2.375$\,]$ | 1.32 | 0.62 | 0.02 | 0.34
$[\,$2.375, 2.500$\,]$ | 1.31 | 1.47 | 0.21 | 0.67
$[\,$2.500, 2.625$\,]$ | 1.32 | 1.42 | 0.13 | 1.17
$[\,$2.625, 2.750$\,]$ | 1.33 | 0.28 | 0.30 | 0.92
$[\,$2.750, 2.875$\,]$ | 1.36 | 0.52 | 0.02 | 0.89
$[\,$2.875, 3.000$\,]$ | 1.40 | 0.39 | 0.42 | 1.13
$[\,$3.000, 3.125$\,]$ | 1.46 | 0.51 | 0.23 | 0.68
$[\,$3.125, 3.250$\,]$ | 1.56 | 0.00 | 0.50 | 0.58
$[\,$3.250, 3.375$\,]$ | 1.63 | 0.00 | 0.07 | 1.01
$[\,$3.375, 3.625$\,]$ | 1.85 | 1.14 | 0.34 | 0.85
$[\,$3.625, 4.000$\,]$ | 2.23 | 0.00 | 0.21 | 0.80
Table 14: Systematic uncertainties on the single differential cross-sections in bins of $p_{\mathrm{T}}^{Z}$, presented in percentage. $p_{\mathrm{T}}^{Z}$ [$\text{\,Ge\kern-1.00006ptV\\!/}c$ ] | Efficiency | Background | FSR | Closure | Calibration | Unfolding
---|---|---|---|---|---|---
$[\,$0.0, 2.2$\,]$ | 1.42 | 0.56 | 1.15 | 0.43 | 0.13 | 0.54
$[\,$2.2, 3.4$\,]$ | 1.41 | 0.50 | 0.64 | 0.37 | 0.40 | 1.73
$[\,$3.4, 4.6$\,]$ | 1.43 | 0.09 | 0.09 | 0.50 | 0.28 | 1.66
$[\,$4.6, 5.8$\,]$ | 1.40 | 0.04 | 0.55 | 0.05 | 0.74 | 1.18
$[\,$5.8, 7.2$\,]$ | 1.43 | 0.25 | 1.27 | 0.12 | 0.91 | 0.36
$[\,$7.2, 8.7$\,]$ | 1.38 | 0.60 | 1.19 | 2.11 | 0.59 | 0.02
$[\,$8.7, 10.5$\,]$ | 1.44 | 0.20 | 0.14 | 0.50 | 0.54 | 0.04
$[\,$10.5, 12.8$\,]$ | 1.39 | 1.34 | 0.48 | 0.97 | 0.31 | 0.03
$[\,$12.8, 15.4$\,]$ | 1.39 | 2.00 | 1.11 | 0.18 | 0.36 | 0.01
$[\,$15.4, 19.0$\,]$ | 1.39 | 1.10 | 0.71 | 1.16 | 0.17 | 0.00
$[\,$19.0, 34.0$\,]$ | 1.37 | 1.09 | 0.09 | 1.97 | 0.10 | 0.00
$[\,$34.0, 120.0$\,]$ | 1.36 | 1.03 | 0.31 | 2.18 | 0.16 | 0.00
Table 15: Systematic uncertainties on the single differential cross-sections in bins of $\phi_{\eta}^{*}$, presented in percentage. $\phi_{\eta}^{*}$ | Efficiency | Background | FSR | Closure
---|---|---|---|---
$[\,$0.00, 0.01$\,]$ | 1.41 | 1.10 | 0.81 | 0.01
$[\,$0.01, 0.02$\,]$ | 1.42 | 0.43 | 0.21 | 0.00
$[\,$0.02, 0.03$\,]$ | 1.41 | 0.14 | 0.36 | 0.12
$[\,$0.03, 0.05$\,]$ | 1.42 | 0.52 | 0.28 | 0.45
$[\,$0.05, 0.07$\,]$ | 1.41 | 0.01 | 0.50 | 0.44
$[\,$0.07, 0.10$\,]$ | 1.39 | 0.29 | 0.24 | 2.22
$[\,$0.10, 0.15$\,]$ | 1.39 | 1.45 | 0.25 | 0.81
$[\,$0.15, 0.20$\,]$ | 1.39 | 0.82 | 0.63 | 1.35
$[\,$0.20, 0.30$\,]$ | 1.38 | 0.77 | 0.97 | 1.00
$[\,$0.30, 0.60$\,]$ | 1.37 | 1.45 | 0.61 | 2.72
$[\,$0.60, 1.20$\,]$ | 1.38 | 0.09 | 0.13 | 2.13

##  Appendix E Numerical results of single differential cross-sections

The measured single differential cross-sections in bins of $y^{{Z}}$,
$p_{\mathrm{T}}^{Z}$ and $\phi_{\eta}^{*}$ are presented in Tables 16 to 18.

Table 16: Measured single differential cross-sections in bins of $y^{{Z}}$. The first uncertainty is statistical, the second systematic, and the third is from the uncertainty on the integrated luminosity. $y^{{Z}}$ |  $d\sigma({Z}\\!\rightarrow{\mu^{+}\mu^{-}})/d{y^{{Z}}}$ [ pb]
---|---
$[\,$2.000, 2.125$\,]$ | 4.4 | $\pm$ | 0.7 | $\pm$ | 0.1 | $\pm$ | 0.1
$[\,$2.125, 2.250$\,]$ | 14.6 | $\pm$ | 1.2 | $\pm$ | 0.2 | $\pm$ | 0.3
$[\,$2.250, 2.375$\,]$ | 24.3 | $\pm$ | 1.5 | $\pm$ | 0.4 | $\pm$ | 0.5
$[\,$2.375, 2.500$\,]$ | 27.5 | $\pm$ | 1.6 | $\pm$ | 0.6 | $\pm$ | 0.6
$[\,$2.500, 2.625$\,]$ | 36.1 | $\pm$ | 1.8 | $\pm$ | 0.8 | $\pm$ | 0.7
$[\,$2.625, 2.750$\,]$ | 39.9 | $\pm$ | 2.0 | $\pm$ | 0.7 | $\pm$ | 0.8
$[\,$2.750, 2.875$\,]$ | 42.6 | $\pm$ | 2.0 | $\pm$ | 0.7 | $\pm$ | 0.9
$[\,$2.875, 3.000$\,]$ | 38.8 | $\pm$ | 1.9 | $\pm$ | 0.7 | $\pm$ | 0.8
$[\,$3.000, 3.125$\,]$ | 30.8 | $\pm$ | 1.7 | $\pm$ | 0.5 | $\pm$ | 0.6
$[\,$3.125, 3.250$\,]$ | 25.5 | $\pm$ | 1.6 | $\pm$ | 0.4 | $\pm$ | 0.5
$[\,$3.250, 3.375$\,]$ | 16.0 | $\pm$ | 1.3 | $\pm$ | 0.3 | $\pm$ | 0.3
$[\,$3.375, 3.625$\,]$ | 6.8 | $\pm$ | 0.6 | $\pm$ | 0.2 | $\pm$ | 0.1
$[\,$3.625, 4.000$\,]$ | 0.7 | $\pm$ | 0.2 | $\pm$ | 0.0 | $\pm$ | 0.0
Table 17: Measured single differential cross-sections in bins of $p_{\mathrm{T}}^{Z}$. The first uncertainty is statistical, the second systematic, and the third is due to the luminosity.  $p_{\mathrm{T}}^{Z}$ [$\text{\,Ge\kern-1.00006ptV\\!/}c$ ] |  $d\sigma({Z}\\!\rightarrow{\mu^{+}\mu^{-}})/d{p_{\mathrm{T}}^{Z}}$ [ pb/($\text{\,Ge\kern-1.00006ptV\\!/}c$)]
---|---
$[\,$0.0, 2.2$\,]$ | 1.769 | $\pm$ | 0.105 | $\pm$ | 0.036 | $\pm$ | 0.035
$[\,$2.2, 3.4$\,]$ | 2.77 | $\pm$ | 0.17 | $\pm$ | 0.07 | $\pm$ | 0.06
$[\,$3.4, 4.6$\,]$ | 3.23 | $\pm$ | 0.19 | $\pm$ | 0.07 | $\pm$ | 0.06
$[\,$4.6, 5.8$\,]$ | 2.71 | $\pm$ | 0.17 | $\pm$ | 0.06 | $\pm$ | 0.05
$[\,$5.8, 7.2$\,]$ | 2.27 | $\pm$ | 0.14 | $\pm$ | 0.05 | $\pm$ | 0.05
$[\,$7.2, 8.7$\,]$ | 2.09 | $\pm$ | 0.13 | $\pm$ | 0.06 | $\pm$ | 0.04
$[\,$8.7, 10.5$\,]$ | 1.806 | $\pm$ | 0.110 | $\pm$ | 0.030 | $\pm$ | 0.036
$[\,$10.5, 12.8$\,]$ | 1.234 | $\pm$ | 0.076 | $\pm$ | 0.028 | $\pm$ | 0.025
$[\,$12.8, 15.4$\,]$ | 0.911 | $\pm$ | 0.063 | $\pm$ | 0.025 | $\pm$ | 0.018
$[\,$15.4, 19.0$\,]$ | 0.759 | $\pm$ | 0.050 | $\pm$ | 0.017 | $\pm$ | 0.015
$[\,$19.0, 34.0$\,]$ | 0.328 | $\pm$ | 0.016 | $\pm$ | 0.009 | $\pm$ | 0.007
$[\,$34.0, 120.0$\,]$ | 0.0325 | $\pm$ | 0.0021 | $\pm$ | 0.0009 | $\pm$ | 0.0007
Table 18: Measured single differential cross-sections in bins of $\phi_{\eta}^{*}$. The first uncertainty is statistical, the second systematic, and the third is due to the luminosity. $\phi_{\eta}^{*}$ |  $d\sigma({Z}\\!\rightarrow{\mu^{+}\mu^{-}})/d{\phi_{\eta}^{*}}$ [ pb]
---|---
$[\,$0.00, 0.01$\,]$ | 468.5 | $\pm$ | 23.8 | $\pm$ | 9.2 | $\pm$ | 9.4
$[\,$0.01, 0.02$\,]$ | 442.4 | $\pm$ | 23.4 | $\pm$ | 6.6 | $\pm$ | 8.8
$[\,$0.02, 0.03$\,]$ | 384.6 | $\pm$ | 21.9 | $\pm$ | 5.6 | $\pm$ | 7.7
$[\,$0.03, 0.05$\,]$ | 311.6 | $\pm$ | 13.8 | $\pm$ | 5.0 | $\pm$ | 6.2
$[\,$0.05, 0.07$\,]$ | 212.1 | $\pm$ | 11.5 | $\pm$ | 3.3 | $\pm$ | 4.2
$[\,$0.07, 0.10$\,]$ | 160.4 | $\pm$ | 8.0 | $\pm$ | 4.2 | $\pm$ | 3.2
$[\,$0.10, 0.15$\,]$ | 87.63 | $\pm$ | 4.51 | $\pm$ | 1.91 | $\pm$ | 1.75
$[\,$0.15, 0.20$\,]$ | 45.59 | $\pm$ | 3.31 | $\pm$ | 1.00 | $\pm$ | 0.91
$[\,$0.20, 0.30$\,]$ | 24.21 | $\pm$ | 1.70 | $\pm$ | 0.51 | $\pm$ | 0.48
$[\,$0.30, 0.60$\,]$ | 6.30 | $\pm$ | 0.49 | $\pm$ | 0.22 | $\pm$ | 0.13
$[\,$0.60, 1.20$\,]$ | 0.571 | $\pm$ | 0.108 | $\pm$ | 0.015 | $\pm$ | 0.011

## References

  * [1] P. J. Rijken and W. L. van Neerven, _Order $\alpha_{s}^{2}$ contributions to the Drell-Yan cross-section at fixed target energies_, Phys. Rev. D51 (1995) 44, arXiv:hep-ph/9408366
  * [2] R. Hamberg, W. L. van Neerven, and T. Matsuura, _A complete calculation of the order $\alpha_{s}^{2}$ correction to the Drell-Yan $K$ factor_, Nucl. Phys. B359 (1991) 343, [Erratum: Nucl. Phys. B644, 403–404 (2002)]
  * [3] R. V. Harlander and W. B. Kilgore, _Next-to-next-to-leading order Higgs production at hadron colliders_ , Phys. Rev. Lett. 88 (2002) 201801, arXiv:hep-ph/0201206
  * [4] W. L. van Neerven and E. B. Zijlstra, _The $O(\alpha_{s}^{2})$ corrected Drell-Yan $K$ factor in the DIS and MS scheme_, Nucl. Phys. B382 (1992) 11, [Erratum: Nucl. Phys. B680, 513–514 (2004)]
  * [5] C. Anastasiou, L. J. Dixon, K. Melnikov, and F. Petriello, _High precision QCD at hadron colliders: Electroweak gauge boson rapidity distributions at NNLO_ , Phys. Rev. D69 (2004) 094008, arXiv:hep-ph/0312266
  * [6] I. Scimemi and A. Vladimirov, _Non-perturbative structure of semi-inclusive deep-inelastic and Drell-Yan scattering at small transverse momentum_ , JHEP 06 (2020) 137, arXiv:1912.06532
  * [7] M. Bury et al., _PDF bias and flavor dependence in TMD distributions_ , JHEP 10 (2022) 118, arXiv:2201.07114
  * [8] S. Camarda, L. Cieri, and G. Ferrera, _Drell–Yan lepton-pair production: qT resummation at N3LL accuracy and fiducial cross sections at N3LO_ , Phys. Rev. D104 (2021) L111503, arXiv:2103.04974
  * [9] C. Duhr and B. Mistlberger, _Lepton-pair production at hadron colliders at N 3LO in QCD_, JHEP 03 (2022) 116, arXiv:2111.10379
  * [10] H1, ZEUS Collaboration, F. D. Aaron et al., _Combined measurement and QCD analysis of the inclusive $e^{+}p$ scattering cross sections at HERA_, JHEP 01 (2010) 109, arXiv:0911.0884
  * [11] H1, ZEUS Collaboration, H. Abramowicz et al., _Combination of measurements of inclusive deep inelastic ${e^{\pm}p}$ scattering cross sections and QCD analysis of HERA data_, Eur. Phys. J. C75 (2015) 580, arXiv:1506.06042
  * [12] CDF collaboration, F. Abe et al., _Forward-backward charge asymmetry of electron pairs above the $Z^{0}$ pole_, Phys. Rev. Lett. 77 (1996) 2616
  * [13] CDF collaboration, D. Acosta et al., _Measurement of the forward-backward charge asymmetry from $W\rightarrow e\nu$ production in $p\bar{p}$ collisions at $\sqrt{s}=1.96$ TeV_, Phys. Rev. D71 (2005) 051104, arXiv:hep-ex/0501023
  * [14] CDF collaboration, T. A. Aaltonen et al., _Measurement of $d\sigma/dy$ of Drell-Yan $e^{+}e^{-}$ pairs in the $Z$ Mass Region from $p\bar{p}$ Collisions at $\sqrt{s}=1.96$ TeV_, Phys. Lett. B692 (2010) 232, arXiv:0908.3914
  * [15] CDF collaboration, T. Aaltonen et al., _Measurement of the inclusive jet cross section at the Fermilab Tevatron $p\bar{p}$ collider using a cone-based jet algorithm_, Phys. Rev. D78 (2008) 052006, arXiv:0807.2204, [Erratum: Phys. Rev. D79, 119902 (2009)]
  * [16] D0 collaboration, V. M. Abazov et al., _Measurement of the muon charge asymmetry from $W$ boson decays_, Phys. Rev. D77 (2008) 011106, arXiv:0709.4254
  * [17] D0 collaboration, V. M. Abazov et al., _Measurement of the ratios of the $Z/\gamma^{*}+\geq$ n jet production cross sections to the total inclusive $Z/\gamma^{*}$ cross section in $p\bar{p}$ collisions at $\sqrt{s}=1.96$ TeV_, Phys. Lett. B658 (2008) 112, arXiv:hep-ex/0608052
  * [18] D0 collaboration, V. M. Abazov et al., _Measurement of the electron charge asymmetry in $p\bar{p}\rightarrow W+X\rightarrow e\nu+X$ decays in $p\bar{p}$ collisions at $\sqrt{s}=1.96$ TeV_, Phys. Rev. D91 (2015) 032007, arXiv:1412.2862, [Erratum: Phys. Rev. D91, 079901 (2015)]
  * [19] D0 collaboration, V. M. Abazov et al., _Measurement of the inclusive jet cross-section in $p\bar{p}$ collisions at $\sqrt{s}=1.96$ TeV_, Phys. Rev. Lett. 101 (2008) 062001, arXiv:0802.2400
  * [20] ATLAS collaboration, G. Aad et al., _Measurement of the inclusive $W^{\pm}$ and $Z/\gamma^{*}$ cross sections in the electron and muon decay channels in $pp$ collisions at $\sqrt{s}=7$ TeV with the ATLAS detector_, Phys. Rev. D85 (2012) 072004, arXiv:1109.5141
  * [21] ATLAS collaboration, G. Aad et al., _Measurement of the transverse momentum and $\phi^{*}_{\eta}$ distributions of Drell–Yan lepton pairs in proton–proton collisions at $\sqrt{s}=8$ TeV with the ATLAS detector_, Eur. Phys. J. C76 (2016) 291, arXiv:1512.02192
  * [22] CMS collaboration, S. Chatrchyan et al., _Measurement of the muon charge asymmetry in inclusive $pp\rightarrow W+X$ production at $\sqrt{s}=7$ TeV and an improved determination of light parton distribution functions_, Phys. Rev. D90 (2014) 032004, arXiv:1312.6283
  * [23] CMS collaboration, S. Chatrchyan et al., _Measurement of the electron charge asymmetry in inclusive $W$ production in $pp$ collisions at $\sqrt{s}=7$ TeV_, Phys. Rev. Lett. 109 (2012) 111806, arXiv:1206.2598
  * [24] ATLAS collaboration, G. Aad et al., _Measurement of the transverse momentum distribution of Drell–Yan lepton pairs in proton–proton collisions at $\sqrt{s}=13$ TeV with the ATLAS detector_, Eur. Phys. J. C80 (2020) 616, arXiv:1912.02844
  * [25] CMS collaboration, A. M. Sirunyan et al., _Measurements of differential Z boson production cross sections in proton-proton collisions at $\sqrt{s}$ = 13 TeV_, JHEP 12 (2019) 061, arXiv:1909.04133
  * [26] LHCb collaboration, R. Aaij et al., _Measurement of the forward $W$ boson production cross-section in $p$ $p$ collisions at $\sqrt{s}=$7 TeV_, JHEP 12 (2014) 079, arXiv:1408.4354
  * [27] LHCb collaboration, R. Aaij et al., _Measurement of the forward $Z$ boson cross-section in $p$ $p$ collisions at $\sqrt{s}=$7 TeV_, JHEP 08 (2015) 039, arXiv:1505.07024
  * [28] LHCb collaboration, R. Aaij et al., _Measurement of ${Z}\\!\rightarrow{e^{+}e^{-}}$ production at $\sqrt{s}=$8 TeV_, JHEP 05 (2015) 109, arXiv:1503.00963
  * [29] LHCb collaboration, R. Aaij et al., _Measurement of forward $W$ and $Z$ boson production in $p$ $p$ collisions at $\sqrt{s}=$8 TeV_, JHEP 01 (2016) 155, arXiv:1511.08039
  * [30] LHCb collaboration, R. Aaij et al., _Measurement of forward ${W}\\!\rightarrow{e}{\nu}$ production in $p$ $p$ collisions at $\sqrt{s}=$8 TeV_, JHEP 10 (2016) 030, arXiv:1608.01484
  * [31] LHCb collaboration, R. Aaij et al., _Measurement of the forward $Z$ boson production cross-section in $p$ $p$ collisions at $\sqrt{s}=$13 TeV_, JHEP 09 (2016) 136, arXiv:1607.06495
  * [32] LHCb collaboration, R. Aaij et al., _Precision measurement of forward $Z$ boson production in proton-proton collisions at $\sqrt{s}=13$ TeV_, JHEP 07 (2022) 026, arXiv:2112.07458
  * [33] S. Dulat et al., _New parton distribution functions from a global analysis of quantum chromodynamics_ , Phys. Rev. D93 (2016) 033006, arXiv:1506.07443
  * [34] L. A. Harland-Lang, A. D. Martin, P. Motylinski, and R. S. Thorne, _Parton distributions in the LHC era: MMHT 2014 PDFs_ , Eur. Phys. J. C75 (2015) 204, arXiv:1412.3989
  * [35] NNPDF collaboration, R. D. Ball et al., _Parton distributions for the LHC Run II_ , JHEP 04 (2015) 040, arXiv:1410.8849
  * [36] T.-J. Hou et al., _New CTEQ global analysis of quantum chromodynamics with high-precision data from the LHC_ , Phys. Rev. D103 (2021) 014013, arXiv:1912.10053
  * [37] NNPDF collaboration, R. D. Ball et al., _Parton distributions from high-precision collider data_ , Eur. Phys. J. C77 (2017) 663, arXiv:1706.00428
  * [38] NNPDF collaboration, R. D. Ball et al., _The path to proton structure at 1% accuracy_ , Eur. Phys. J. C82 (2022) 428, arXiv:2109.02653
  * [39] S. Bailey et al., _Parton distributions from LHC, HERA, Tevatron and fixed target data: MSHT20 PDFs_ , Eur. Phys. J. C81 (2021) 341, arXiv:2012.04684
  * [40] K. J. Eskola, H. Paukkunen, and C. A. Salgado, _EPS09: A new generation of NLO and LO nuclear parton distribution functions_ , JHEP 04 (2009) 065, arXiv:0902.4154
  * [41] A. D. Martin, W. J. Stirling, R. S. Thorne, and G. Watt, _Parton distributions for the LHC_ , Eur. Phys. J. C63 (2009) 189, arXiv:0901.0002
  * [42] LHCb collaboration, R. Aaij et al., _Observation of $Z$ production in proton-lead collisions at LHCb_, JHEP 09 (2014) 030, arXiv:1406.2885
  * [43] LHCb collaboration, R. Aaij et al., _Precision luminosity measurements at LHCb_ , JINST 9 (2014) P12005, arXiv:1410.0149
  * [44] A. Banfi et al., _Optimisation of variables for studying dilepton transverse momentum distributions at hadron colliders_ , Eur. Phys. J. C71 (2011) 1600, arXiv:1009.1580
  * [45] LHCb collaboration, A. A. Alves Jr. et al., _The LHCb detector at the LHC_, JINST 3 (2008) S08005
  * [46] LHCb collaboration, R. Aaij et al., _LHCb detector performance_ , Int. J. Mod. Phys. A30 (2015) 1530022, arXiv:1412.6352
  * [47] R. Aaij et al., _Performance of the LHCb Vertex Locator_ , JINST 9 (2014) P09007, arXiv:1405.7808
  * [48] LHCb collaboration, _LHCb reoptimized detector design and performance: Technical Design Report_ , CERN-LHCC-2003-030, 2003
  * [49] P. d’Argent et al., _Improved performance of the LHCb Outer Tracker in LHC Run 2_ , JINST 12 (2017) P11016, arXiv:1708.00819
  * [50] A. A. Alves Jr. et al., _Performance of the LHCb muon system_ , JINST 8 (2013) P02022, arXiv:1211.1346
  * [51] R. Aaij et al., _The LHCb trigger and its performance in 2011_, JINST 8 (2013) P04022, arXiv:1211.3055
  * [52] T. Sjöstrand, S. Mrenna, and P. Skands, _PYTHIA 6.4 physics and manual_ , JHEP 05 (2006) 026, arXiv:hep-ph/0603175
  * [53] C. Bierlich et al., _A comprehensive guide to the physics and usage of PYTHIA 8.3_ , doi: 10.21468/SciPostPhysCodeb.8 arXiv:2203.11601
  * [54] I. Belyaev et al., _Handling of the generation of primary events in Gauss, the LHCb simulation framework_ , J. Phys. Conf. Ser. 331 (2011) 032047
  * [55] N. Davidson, T. Przedzinski, and Z. Was, _PHOTOS interface in C++: Technical and physics documentation_ , Comp. Phys. Comm. 199 (2016) 86, arXiv:1011.0937
  * [56] Geant4 collaboration, J. Allison et al., _Geant4 developments and applications_ , IEEE Trans. Nucl. Sci. 53 (2006) 270
  * [57] M. Clemencic et al., _The LHCb simulation application, Gauss: Design, evolution and experience_, J. Phys. Conf. Ser. 331 (2011) 032023
  * [58] J. Campbell and T. Neumann, _Precision Phenomenology with MCFM_ , JHEP 12 (2019) 034, arXiv:1909.09117
  * [59] LHCb collaboration, R. Aaij et al., _Measurement of the $W$ boson mass_, JHEP 01 (2022) 036, arXiv:2109.01113
  * [60] LHCb collaboration, R. Aaij et al., _Measurement of the track reconstruction efficiency at LHCb_ , JINST 10 (2015) P02007, arXiv:1408.1251
  * [61] H. B. Prosper and L. Lyons, eds., Proceedings, PHYSTAT 2011 Workshop on Statistical Issues Related to Discovery Claims in Search Experiments and Unfolding, CERN,Geneva, Switzerland 17-20 January 2011, CERN Yellow Reports: Conference Proceedings, (Geneva), CERN, 2011.  doi: 10.5170/CERN-2011-006
  * [62] G. D’Agostini, _A multidimensional unfolding method based on Bayes’ theorem_ , Nucl. Instrum. Meth. A362 (1995) 487
  * [63] C. Balazs and C. P. Yuan, _Soft gluon effects on lepton pairs at hadron colliders_ , Phys. Rev. D56 (1997) 5558, arXiv:hep-ph/9704258
  * [64] T. Adye, _Unfolding algorithms and tests using RooUnfold_ , CERN-2011-006 (2011) 313, arXiv:1105.1160
  * [65] J. C. Collins, D. E. Soper, and G. F. Sterman, _Transverse momentum distribution in Drell-Yan pair and $W$ and $Z$ boson production_, Nucl. Phys. B250 (1985) 199
  * [66] J. C. Collins and D. E. Soper, _Back-to-back jets in QCD_ , Nucl. Phys. B193 (1981) 381, [Erratum: Nucl. Phys. B 213, 545 (1983)]
  * [67] J. C. Collins and D. E. Soper, _Back-to-back jets: fourier transform from B to K-transverse_ , Nucl. Phys. B197 (1982) 446
  * [68] P. Nason, _A New method for combining NLO QCD with shower Monte Carlo algorithms_ , JHEP 11 (2004) 040, arXiv:hep-ph/0409146
  * [69] S. Frixione, P. Nason, and C. Oleari, _Matching NLO QCD computations with Parton Shower simulations: the POWHEG method_ , JHEP 11 (2007) 070, arXiv:0709.2092
  * [70] S. Alioli, P. Nason, C. Oleari, and E. Re, _NLO vector-boson production matched with shower in POWHEG_ , JHEP 07 (2008) 060, arXiv:0805.4802
  * [71] 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, arXiv:1002.2581
  * [72] T. Sjöstrand, S. Mrenna, and P. Skands, _A brief introduction to PYTHIA 8.1_ , Comput. Phys. Commun. 178 (2008) 852, arXiv:0710.3820
  * [73] J. Pumplin et al., _New generation of parton distributions with uncertainties from global QCD analysis_ , JHEP 07 (2002) 012, arXiv:hep-ph/0201195
  * [74] Y. Li and F. Petriello, _Combining QCD and electroweak corrections to dilepton production in FEWZ_ , Phys. Rev. D86 (2012) 094034, arXiv:1208.5967