B2G-21-004

B2G-21-004

[cern]The CMS Collaboration

# Search for pair-produced vector-like leptons in final states with third-
generation leptons and at least three quark jets in proton-proton collisions
at $\sqrt{s}=13\TeV$

###### Abstract

The first search is presented for vector-like leptons (VLLs) in the context of
the “4321 model”, an ultraviolet-complete model with the potential to explain
existing physics measurements that are in tension with standard model
predictions. The analyzed data, corresponding to an integrated luminosity of
96.5, were recorded in 2017 and 2018 with the CMS detector at the LHC in
proton-proton collisions at $\sqrt{s}=13\TeV$. Final states with ${\geq}3$
-tagged jets and two third-generation leptons ($\PGt\PGt$, $\PGt\PGnGt$, or
$\PGnGt\PGnGt$) are considered. Upper limits are derived on the VLL production
cross section in the VLL mass range 500–1050. The maximum likelihood fit
prefers the presence of signal at the level of 2.8 standard deviations, for a
representative VLL mass point of 600. As a consequence, the observed upper
limits are approximately double the expected limits.

##  0.1 Introduction

We present a search for vector-like leptons (VLLs), considering the mass range
500–1050, in the context of the “4321 model” [1, 2, 3]. The search uses LHC
proton-proton ($\Pp\Pp$) collision data collected by the CMS experiment in
2017 and 2018 corresponding to a total integrated luminosity of 96.5 [4, 5,
6].

The 4321 model is an ultraviolet-complete model that extends the standard
model (SM) gauge groups to a larger
$\text{SU}(4)\times\text{SU}(3)^{\prime}\times\text{SU}(2)_{L}\times\text{U}(1)^{\prime}$
group. It is motivated by recent measurements of hadron decays that are in
tension with the SM. This tension particularly concerns $R({\PD})$,
$R({\PD}^{*})$, $R({\PK})$, and $R({\PK^{*}})$ measurements [7, 8, 9, 10, 11,
12, 13, 14, 15, 16] that provide evidence for lepton flavour nonuniversality.
The results of the search described here are the first search results for VLLs
in the context of the 4321 model and are also expected to hold for other
models with a Pati–Salam leptoquark explaining the B anomalies [17, 18]. Final
states are probed that are not explored by existing searches for VLLs in
models without leptoquarks [19].

The 4321 model gives a possible explanation for the aforementioned flavour-
nonuniversal results, while simultaneously respecting many other measurements
that are in agreement with SM expectations and lepton flavour universality
[20, 21, 22, 23]. Here, we search for pair production of the lightest new
particles in this model, the vector-like leptons. These leptons are referred
to as vector-like because they are non-chiral, their left- and right-handed
components have the same charges.

The VLLs come in electroweak doublets with one charged VLL, $\mathrm{E}$, and
one neutral VLL, $\mathrm{N}$, whose masses are taken to be equal. For the
model to explain the physics anomalies while also remaining consistent with
other measurements, the mass of the VLLs cannot be too large. In particular,
requiring compatibility with the measured $R({\PD})$and
$R({\PD}^{*})$anomalies and with measurements of - mixing suggests that the
VLLmass cannot be more than a few [24].

The VLLs may be produced via electroweak production through their couplings to
the SM and / bosons or via interactions with a new boson introduced by the
4321 model. In this paper, we consider only electroweak production, and ignore
potential contributions from the boson; the expected size of the contribution
is dependent on the unknown mass and coupling strength. Despite not explicitly
considering this contribution, the results are also expected to hold when
including production. Examples of Feynman diagrams showing the electroweak
pair production of VLLs, as well as diagrams of the VLL decays, are shown in
Fig. 1.

Figure 1: Left side: example Feynman diagrams showing electroweak production
of VLLpairs through $s$-channel bosons, as expected at the LHC. In these
diagrams, $\mathrm{L}$ represents either the neutral VLL, $\mathrm{N}$, or the
charged VLL, $\mathrm{E}$. Right side: vector-like lepton decays are mediated
by a vector leptoquark, $\mathrm{U}$. These decays are primarily to third-
generation leptons and quarks.

The VLLs decay, via a virtual leptoquark, $\mathrm{U}$, to two quarks and one
lepton. In order to explain the anomalies, the leptoquark is expected to
couple most strongly to the third generation and therefore its decays are
expected to be almost entirely to third-generation fermions. For each second-
generation fermion, a one to two order of magnitude suppression in the
branching fraction is expected, and even larger suppressions are expected for
any first-generation fermion. For this analysis, we have considered only the
dominant decays, decays to third-generation fermions.

The analysis selection is driven by the flavour content of the VLLdecay
products. Given the expectation of two third-generation quarks in every
VLLdecay, we search for pairs of VLLs by selecting events with a high jet
multiplicity. These events are further categorized by the number of leptons.
For each multiplicity, dedicated selections are made to divide the category
into a signal-enriched “signal region” (SR) and one or multiple background-
enriched “control regions” (CRs). Table 1 shows the multiplicity categories
and the decay modes of the different VLLpairs that contribute to each
category. While topologies with electrons or muons in the final state (coming
from top quark or tau decays) are possible, we focus only on all-hadronic
final states in this analysis.

Table 1:  Illustrative contributions from different VLLproduction and decay modes to the 0-, 1-, and 2-signal regions. The decay products in parentheses represent the objects coming from the intermediate vector leptoquark, $\mathrm{U}$, in the decay. In the table, no distinction is made between particles and antiparticles, the multiplicities of each decay mode are not shown, and the impacts of object misidentification are not considered. The charged and neutral VLLs are represented by $\mathrm{E}$ and $\mathrm{N}$; $\mathrm{j}$ represents any quark other than or .  Tau | VLLproduction | Final
---|---|---
multiplicity | \+ decay mode | state
0 $\PGt$ | $\mathrm{E}\mathrm{E}\to\PQb(\PQt\PGnGt)\PQb(\PQt\PGnGt)$ | $4\PQb+4\mathrm{j}+2\PGnGt$
$\mathrm{E}\mathrm{N}\to\PQb(\PQt\PGnGt)\PQt(\PQt\PGnGt)$ | $4\PQb+6\mathrm{j}+2\PGnGt$
$\mathrm{N}\mathrm{N}\to\PQt(\PQt\PGnGt)\PQt(\PQt\PGnGt)$ | $4\PQb+8\mathrm{j}+2\PGnGt$
1 $\PGt$ | $\mathrm{E}\mathrm{E}\to\PQb(\PQb\PGt)\PQb(\PQt\PGnGt)$ | $4\PQb+2\mathrm{j}+\PGt+\PGnGt$
$\mathrm{E}\mathrm{N}\to\PQb(\PQt\PGnGt)\PQt(\PQb\PGt)$ | $4\PQb+4\mathrm{j}+\PGt+\PGnGt$
$\mathrm{E}\mathrm{N}\to\PQb(\PQb\PGt)\PQt(\PQt\PGnGt)$ | $4\PQb+4\mathrm{j}+\PGt+\PGnGt$
$\mathrm{N}\mathrm{N}\to\PQt(\PQb\PGt)\PQt(\PQt\PGnGt)$ | $4\PQb+6\mathrm{j}+\PGt+\PGnGt$
2 $\PGt$ | $\mathrm{E}\mathrm{E}\to\PQb(\PQb\PGt)\PQb(\PQb\PGt)$ | $4\PQb+2\PGt$
$\mathrm{E}\mathrm{N}\to\PQb(\PQb\PGt)\PQt(\PQb\PGt)$ | $4\PQb+2\mathrm{j}+2\PGt$
$\mathrm{N}\mathrm{N}\to\PQt(\PQb\PGt)\PQt(\PQb\PGt)$ | $4\PQb+4\mathrm{j}+2\PGt$

A maximum likelihood fit is performed simultaneously across all multiplicity
categories, including both SRs and CRs, to determine the signal strength for
each VLLmass hypothesis. Within each region, differential distributions are
used as input to the fit. The distributions are chosen to provide additional
separation between the signal and backgrounds. In the fit, the data are also
separated by the year in which it was collected (either 2017 or 2018), to
ensure a good description of the data-taking and detector conditions in each
case.

Machine learning techniques are used to build two deep neural network (DNN)
classifiers [25, 26, 27] to separate signal events from background events. One
of the classifiers, $\mathrm{DNN}_{\mathrm{QCD}}$, is trained to distinguish
signal from quantum chromodynamic (QCD) multijet events and is used to define
the signal region in the 0- category. The other DNN, $\mathrm{DNN}_{\ttbar}$,
is trained to discriminate between signal and events and its output
distribution is used in the fit for the 1- and 2- categories.

Tabulated results for this analysis are provided in HEPData [28].

##  0.2 The CMS detector and object reconstruction

The central feature of the CMS apparatus is a superconducting solenoid of 6m
internal diameter, providing a magnetic field of 3.8T. Within the solenoid
volume are a silicon pixel and strip tracker, a lead tungstate crystal
electromagnetic calorimeter (ECAL), and a brass and scintillator hadron
calorimeter (HCAL), each composed of a barrel and two endcap sections. Forward
calorimeters extend the pseudorapidity ($\eta$) coverage provided by the
barrel and endcap detectors. Muons are measured in gas-ionization detectors
embedded in the steel flux-return yoke outside the solenoid. A more detailed
description of the CMS detector, together with a definition of the coordinate
system used and the relevant kinematic variables, can be found in Ref. [29].

Events of interest are selected using a two-tiered trigger system. The first
level, composed of custom hardware processors, uses information from the
calorimeters and muon detectors to select events at a rate of around 100kHz
within a fixed latency of about 4 [30]. The second level, known as the high-
level trigger, consists of a farm of processors running a version of the full
event reconstruction software optimized for fast processing, and reduces the
event rate to around 1kHz before data storage [31]. The primary vertex is
taken to be the vertex corresponding to the hardest scattering in the event,
evaluated using tracking information alone, as described in Section 9.4.1 of
Ref. [32].

Reconstruction is performed using the particle-flow algorithm [33], which aims
to reconstruct and identify each individual particle in an event with an
optimized combination of information from the various elements of the CMS
detector. The energy of photons is obtained from the ECAL measurement. The
energy of electrons is determined from a combination of the electron momentum
at the primary interaction vertex as determined by the tracker, the energy of
the corresponding ECAL cluster, and the energy sum of all bremsstrahlung
photons spatially compatible with originating from the electron track. The
energy of muons is obtained from the curvature of the corresponding track. The
energy of charged hadrons is determined from a combination of their momentum
measured in the tracker and the matching ECAL and HCAL energy deposits,
corrected for the response function of the calorimeters to hadronic showers.
Finally, the energy of neutral hadrons is obtained from the corresponding
corrected ECAL and HCAL energies.

For each event, hadronic jets are clustered from the reconstructed particles
using the infrared and collinear safe anti-algorithm [34, 35] with a distance
parameter of 0.4. Jet momentum is determined as the vectorial sum of all
particle momenta in the jet, and is found from simulation to be, on average,
within 5 to 10% of the true momentum over the whole transverse momentum ()
spectrum and detector acceptance. Additional $\Pp\Pp$ interactions within the
same or nearby bunch crossings can contribute additional tracks and
calorimetric energy depositions, increasing the apparent jet momentum. To
mitigate this effect, tracks identified to be originating from pileup vertices
are discarded and an offset correction is applied to correct for remaining
contributions [36]. Jet energy corrections are derived from simulation studies
so that the average measured energy of jets becomes identical to that of
particle level jets. In situ measurements of the momentum balance in dijet,
photon+jet, Z+jet, and multijet events are used to determine any residual
differences between the jet energy scale in data and in simulation, and
appropriate corrections are made [37]. Additional selection criteria are
applied to each jet to remove jets potentially dominated by instrumental
effects or reconstruction failures. All jets considered in the analysis are
required to have $\pt>40\GeV$ and $\abs{\eta}\Xspace<2.5$.

To identify jets originating from quarks, -tagging algorithms are applied to
all reconstructed jets in an event. To identify jets at the level of the
trigger selections, the CSV v2 algorithm is used in the 2017 data sample and
the DeepCSV algorithm is used in the 2018 sample [38]. In the offline
selection and analysis, the DeepJet algorithm is used to identify jets [39,
40] for both data-taking years.

Each of the tagging algorithms outputs a continuous discriminator value. By
requiring this value to be above a given threshold (“working point”), a jet
can be identified as a -tagged jet. The tagging efficiency for jets is 79
(82)% for a 1% misidentification probability of light flavour quark and gluon
jets in 2017 (2018); this is referred to as the medium working point (WP).
These values, taken from simulation, are further corrected towards the data
using dedicated measurements that correct the identification and
misidentification probabilities as functions of the kinematical properties of
the jets. In this analysis, the tagging discriminant score is used as input to
the analysis-specific machine learning discriminators, and therefore
corrections are applied to the full simulated discriminant score distribution,
rather than just to the efficiencies of the individual -tagging WPs.

Hadronically decaying $\PGt$ leptons ($\tauh$) are reconstructed from jets,
using the “hadrons-plus-strips” algorithm described in Ref. [41], which
combines 1–3 tracks with energy deposits in the calorimeters, to identify the
decay modes. Modes with two reconstructed tracks target cases in which one of
the tracks from a $\PGt$ decay with three charged particles in the final state
is not reconstructed. To take into account the effects of pair conversion and
bremsstrahlung in the material in front of the ECAL, neutral pions are
reconstructed by combining electron and photon candidates within a defined
region in $\eta$–$\phi$ space, where $\phi$ is the azimuthal angle. The
dimensions of this region are in the range 0.05–0.15 in $\eta$ and 0.05–0.30
in $\phi$ , depending on the of the electron and photon candidates. All
candidates in the analysis are required to have $\pt>20\GeV$ and
$\abs{\eta}\Xspace<2.1$.

To distinguish genuine $\tauh$ decays from jets originating from the
hadronization of quarks or gluons, and from electrons or muons, the DeepTau
algorithm is used [42]. Information from all particles reconstructed near the
$\tauh$ axis is combined with properties of the $\tauh$ candidate and the
event.

Several different WPs of the DeepTau algorithm are used: the very loose WP (in
the rest of the paper referred to simply as the “loose WP“) , medium WP, and
tight WP. The rate at which jets are misidentified as a is given in Table 0.2,
as estimated from simulated samples of boson production in association with
jets.

Estimated identification efficiencies and jet misidentification rates for the
DeepTau algorithm used in the analysis.  Working point Efficiency Jet
misidentification rate Loose 90% 10% Medium 70% 0.4% Tight 60% 0.2%

The missing transverse momentum vector is computed as the negative vector sum
of the transverse momenta of all the particle-flow candidates in an event, and
its magnitude is denoted as [43]. The is modified to account for corrections
to the energy scale of the reconstructed jets in the event.

##  0.3 Data and simulated samples

The analysis is performed on a sample of $\Pp\Pp$ collision data with
$\sqrt{s}=13\TeV$, collected with the CMS detector in 2017–2018 and
corresponding to a total integrated luminosity of 96.5 [4, 5, 6]. The 2017 and
2018 data sets benefit from the upgrade of the pixel tracking detector
installed in the winter of 2016–2017 [44]. The upgraded tracker significantly
improves the jet identification performance [45], which is crucial to the
analysis and online event selection. The 2016 data set, collected prior to
this upgrade, is therefore not considered.

The signal is modelled at leading order (LO) with [46] v2.6.5, using the
universal FeynRules outputinterface [47] to implement the 4321 model. The
matrix element calculations are interfaced to 8.2 [48], with the description
of the underlying event is provided by the CP5 tune [49]. The NNPDF 3.1 next-
to-next-to-leading order (NNLO) global fit parton distribution function (PDF)
set [50] is used.

Signal production is restricted to use only SM electroweak couplings. We
studied the impact of contributions from ’ production and found that they do
not alter our conclusions. In the simulated samples, the couplings of the
leptoquark are set to 0 for all first- and second-generation SM fermions. The
mass of the leptoquark is always taken to be 3.5and the masses of the two VLLs
are always taken to be equal to each other. The simulated model parameters,
including the leptoquark and ’ masses, are taken from benchmark scenarios
proposed to explain the anomalies [24]. However, due to the much lighter mass
scale of the VLLs, the results are mostly insensitive to the values of the
boson masses.

We generate signal samples for 12 different VLL mass hypotheses, ranging from
500 to 1050, in 50increments. The widths of the VLLs are evaluated for each
mass. The decay of the VLLs is performed within , however, the decays of any
top quarks produced in the VLL decays are simulated with .

Simulated samples are also used for background estimation. Samples of events
are simulated at next-to-leading order using powheg box [51, 52, 53]
interfaced to the NNPDF3.1 NNLO PDF set [50]. These samples are generated
assuming a 5-flavour scheme, where quarks are considered massless. The top
quark mass is set to 172.5, and the hdamp parameter, which regulates emissions
in , is set to 237.8775. The inclusive production cross section is taken to be
831.76pb, as calculated at NNLO in QCD, including resummation of next-to-next-
to-leading logarithmic soft gluon terms with Top++ 2.0 [54, 55, 56, 57, 58,
59, 60].

Because the signal events are expected to produce four quarks, the component
of the simulation plays a particularly important role in the analysis. This
component is taken from the simulation described above. Events are taken as
part of the component if there is at least one jet containing a hadron not
coming from a top quark decay. Otherwise the event is categorized as jj.
Recent measurements from both the ATLAS and CMS Collaborations indicate that
the cross section of the component is underpredicted by the simulation [61,
62, 63, 64]. In order to account for this discrepancy, an overall scale factor
of 1.2 is applied to the component.

Smaller backgrounds coming from production in association with a boson,
denoted X, from vector boson production with additional associated jets,
denoted +jets, and from multiboson production, denoted (), are also
considered. The contributions from these backgrounds are taken from
simulation. These processes are all simulated using for the matrix element
calculation, which is interfaced with 8 for showering and hadronization.

Samples of QCD multijet events are generated at LO using , merging two-,
three-, and four-jet matrix elements, and interfaced to 8 for showering and
hadronization. These samples are only used for training the
$\mathrm{DNN}_{\mathrm{QCD}}$ classifier; in the analysis, background
contributions from QCD multijet processes are estimated using control samples
in data.

For all simulated data, additional minimum bias $\Pp\Pp$ collisions are
simulated and superimposed on the primary simulated interaction. Simulation of
the passage of particles through the detector material is performed by [65].

##  0.4 Event selection and categorization

All events considered in the analysis must satisfy a trigger that requires
three -tagged jets and places additional requirements on the of the jets in
the event. A further selection, after the trigger, is used in the analysis to
avoid the region near the trigger thresholds, which may be sensitive to
mismodelling. In order to pass this basic selection, an event is required to
contain at least three jets passing the medium tagging WP and at least four
jets with greater than 80, 65, 50, and 50. In addition, the scalar sum of the
of all jets in the event () must be greater than 400, and no isolated
electrons ($\pt>15\GeV,\abs{\eta}\Xspace<2.5$) or muons
($\pt>30\GeV,\abs{\eta}\Xspace<2.5$) must be present.

These events are then further categorized based on the candidate multiplicity.
Three categories are considered, with either 0, 1, or $\geq 2{\tauh}$
candidates. Within each of these categories, SRs and CRs are defined. The 1-
and 2- categories also make use of several determination regions, which are
not themselves included in the signal extraction fit, but are used for the
background estimation described in Section 0.6. A diagram summarizing all
these regions is shown in Fig. 2.

Figure 2: Diagram of the event categorization and the various signal and
control regions used in the analysis. The regions within the solid box are all
used in the signal extraction fit. The regions in the dashed box are used to
determine transfer factors from the control regions to the signal region for
events with jets misidentified as ’s. The selections are mutually exclusive
between all regions, so that each event only enters a single region. For
brevity, not all selection criteria are shown; the detailed selection criteria
are described in the text.

The 0- region is defined by requiring that there are no candidates passing the
loose, medium, or tight WP. No additional event selections beyond the basic
ones described above are required. The primary background in this region is
from QCD multijet events. This background in the SR is estimated using a
technique that makes use of three background-enriched regions in data. The SR
is defined by requiring $\ptmiss>160\GeV$, as expected from the two neutrinos,
as well as by making a selection on the $\mathrm{DNN}_{\mathrm{QCD}}$
(described in Section 0.5). The jet multiplicity distribution,
$\mathrm{N}_{\text{jets}}$, is used for the final fit.

The 1- region is defined by requiring that there is exactly one candidate. In
the SR, the candidate must pass the tight WP and there must be no other
candidates passing the medium WP. In the CR, the candidate must pass the loose
WP and there must be no candidates passing the medium WP. Additionally, in
both the SR and the CR, the candidates must have $\pt>60\GeV$, and the event
must have $\ptmiss>60\GeV$. In the 1- region, significant backgrounds arise
from both and QCD multijetevents. Multijet events can pass the selection
criteria when a jet has been misidentified as a . The background has
components both with and without genuine ’s. The
$\mathrm{DNN}_{\ttbar}$classifier distribution, described in Section 0.5, is
used for the final fit.

The 2- region is defined by requiring that there are at least two candidates.
In our data set, no events with more than two candidates are found. The SR
requires that there are at least two candidates passing the medium WP; the CR
requires that there are at least two candidates passing the loose WP and no
candidates passing the medium WP. Additionally, for all events in 2- region,
the candidates must have opposite charge, and be separated by
$\Delta{R}\Xspace=\sqrt{\smash[b]{(\Delta\eta)^{2}+(\Delta\phi)^{2}}}$ of at
least 0.5. The event is required to have $\ptmiss>40\GeV$. In this region, is
the dominant background and the QCD multijet contribution is highly
suppressed. The same $\mathrm{DNN}_{\ttbar}$used in the 1- region is also used
in the final fit of the 2- region.

##  0.5 Machine learning classifiers

Two DNN classifiers are trained to separate signal and background events. One
of them, $\mathrm{DNN}_{\mathrm{QCD}}$, is specifically trained to
discriminate between signal and QCD multijetevents. The other,
$\mathrm{DNN}_{\ttbar}$, is trained to distinguish signal from events.

Both DNNs use the same underlying architecture, a graph neural network called
ABCNet [26]. In this implementation, individual final state objects are passed
as inputs to the DNN, and are represented as nodes in the graph. The graph
connects sets of $\mathrm{k}$-nearest neighbour objects with edges, so that
relationships between the objects are taken into account. In all cases, we use
$\mathrm{k}=8$ and a Euclidean metric in the $\eta$–$\phi$ coordinates.

Both candidates and jets are passed to ABCNet as inputs, with up to ten
objects per event used as input. If more than ten objects are present in the
event, all of the candidates are passed to the network first, followed by the
jets, ordered by , until the total object count has reached ten.

For each object, several kinematic features are passed to ABCNet. They are
shown in Table 0.5.

List of input features used during training of the ABCNet model for
VLLclassification.  Variable Description $\eta$ Pseudorapidity $\phi$
Azimuthal angle $\log(\frac{\pt}{\GeV})$ Logarithm of the
$\log(\frac{m}{\GeV})$ Logarithm of the mass $Q$ Charge DeepJet score tagging
discriminant

For jets, the charge entry is always set to zero. For $\tauh$ candidates, the
DeepJet score is always set to zero, which corresponds to an extremely low
probability of being a jet. With this information, the model is able to
distinguish $\PGt$ leptons from jets during training. Using the same model,
two different trainings are performed to distinguish VLLevents from QCD and
events. Signal samples with different mass hypotheses are used in the
training. The simulated events used for training the
$\mathrm{DNN}_{\mathrm{QCD}}$ classifier are required to pass the 0- selection
criteria, whereas those used for training the $\mathrm{DNN}_{\ttbar}$
classifier must have at least one candidate passing the medium WP.

##  0.6 Background estimation

The primary SM processes that pass the event selection and form backgrounds in
the search are and QCD multijetproduction. Additional contributions from X,
+jets, and () production are also considered, but play smaller roles. The
estimated contributions of these minor background processes are taken directly
from simulation, and the +jets and () contributions are found to be
negligible. On the other hand, the QCD multijet contribution is estimated
using control samples in data.

The expected contributions from events, including the component, are estimated
mostly from simulation, as described in Section 0.3. The number of events with
misidentified candidates, where at least one identified candidate does not
come from a true hadronic $\PGt$ lepton decay, is estimated from a combination
of simulation and control samples in data. The estimation is performed by
means of a “fake factor” method. The method uses a CR with a looser WP to
estimate the number of events in the SR with a tighter WP. The ratio of the
number of events in the two regions, called the fake factor (FF), is first
measured in a dedicated determination region (DR) which is orthogonal to the
primary analysis region. Then in the primary analysis region, the estimated
number of events in the control region is multiplied by the fake factor to
estimate the number of events in the SR:

FF
$\displaystyle\equiv\frac{N^{\text{misid}}\Xspace(\text{SR}_{\mathrm{DR}})}{N^{\text{misid}}\Xspace(\text{CR}_{\mathrm{DR}})},$
(1) $\displaystyle N^{\text{misid}}\Xspace(\text{SR})$
$\displaystyle=N^{\text{misid}}\Xspace(\text{CR})\text{FF},$

where $N^{\text{misid}}\Xspace(X)$ is the number of events in region $X$ with
at least one misidentified candidate. The FF and $N^{\text{misid}}\Xspace$ are
estimated in bins of the $\mathrm{DNN}_{\ttbar}$ score.

To estimate $N^{\text{misid}}\Xspace$ in the DRs, the number of events in data
with correctly identified objects is estimated from simulation and subtracted
from the total; the remaining component is taken to be
$N^{\text{misid}}\Xspace$. In the primary analysis region,
$N^{\text{misid}}\Xspace(\text{CR})$ is taken directly from simulation. An
additional uncertainty is applied to this component to account for possible
mismodelling.

The DRs are designed to be enriched in events. After passing a single-muon
trigger, these events are required to pass the same selection as the main
analysis, except that they must have at least two (not three) -tagged jets.
Events are categorized based on the number of candidates as in the main
analysis region. However, unlike in the main analysis regions, no additional
selection criteria (on di-$\PGt$ charge, $\tauh$ , or ) are applied when
performing the categorization in the DRs. These differences in selection are
all designed to increase the number of events in these regions, as the primary
limitation of our FF estimation is the statistical uncertainty. An additional
correction and uncertainty, to account for these differences in the phase
space between the determination and primary analysis region, is estimated from
simulation.

The fake factor method is applied independently for the 1- and 2- regions,
using the same methodology in both cases.

The QCD multijet background in the 1- and 2- regions is also determined using
a fake factor method. Unlike the approach adopted for , where
$N^{\text{misid}}\Xspace(\text{CR})$ is taken from simulation, this component
is also taken from data for the QCD multijet background. The number of
expected events for all other background processes is subtracted from the data
to estimate $N^{\text{misid}}\Xspace(\text{CR})$. The signal contribution is
ignored in the subtraction, since it is never expected to contribute more than
$\mathcal{O}(1\%)$ of the QCD multijetyield in the control region. The FF DR
is constructed by requiring $\ptmiss<40\GeV$. As for the fake factor, a
correction is applied to account for the difference between the FFs in the DR
and the primary analysis region. However, in the 2- region, no simulated QCD
multijet events pass the signal region selection criteria and therefore no
correction is applied.

In the 0- region, the QCD multijet contribution is estimated using an “ABCD”
method. This method makes use of two observables, over which the background
distribution is uncorrelated. Each observable is divided into two subregions,
to form four regions in total. Then, from the fact that the two observables
are uncorrelated, the background contribution in a signal-enriched region is
estimated from the background yields in three signal-depleted regions:

$N^{\mathrm{QCD}}(\text{SR})=\frac{N^{\mathrm{QCD}}(\text{CR1})N^{\mathrm{QCD}}(\text{CR3})}{N^{\mathrm{QCD}}(\text{CR2})},$
(2)

where $N^{\mathrm{QCD}}(X)$ represents the number of QCD multijetevents in
region $X$ and is determined independently for each bin of the
$\mathrm{N}_{\text{jets}}$distribution used in the fit. We use the and the
$\mathrm{DNN}_{\mathrm{QCD}}$ classifier score as the two observables for the
ABCD method. The regions are divided at
$\mathrm{DNN}_{\mathrm{QCD}}\Xspace=0.6$ and $\ptmiss=60\GeV$, with the region
$\mathrm{DNN}_{\mathrm{QCD}}\Xspace>0.6,\ptmiss>60\GeV$ being the SR. The
training of the DNN is designed to ensure that its output and the are
uncorrelated [27], and this assumption was further tested both in simulated
samples and in a dedicated validation region.

##  0.7 Uncertainties

The dominant uncertainty comes from the limited number of events in the 1- and
2- SRs. A number of other experimental and theoretical uncertainties are
considered and included in the fit as nuisance parameters. In each case, the
nuisance parameters alter one or more of the processes included in the fit.
Except where otherwise noted, the nuisance parameters alter both the shape and
normalization of the affected processes.

Uncertainties in the measurement of the energies and resolutions of jets and
candidates are estimated from the calibrations described in Section 0.2.

Uncertainties in the -tagging efficiencies and misidentification rates are
taken from dedicated measurements [38]. The tagging rates and their
uncertainties are determined by looking separately at events enriched in
genuine jets, using a selection, and events enriched in light flavour quark
and gluon jets, using $\PZ\to\Pe^{+}\Pe^{-}$ and $\PZ\to\PGm^{+}\PGm^{-}$
events with additional jets. An iterative fit to these two regions is
performed and the uncertainties are characterized by nine separate components,
which account for different uncertainty sources such as the limited sample
size and uncertainties in the flavour compositions of the regions.

Uncertainties in the identification and misidentification rates of candidates
are also considered. The identification rate uncertainties are taken from
measurements using a set of $\PZ\to\PGt\PGt$ enriched events, where one of the
$\PGt$ leptons decays into neutrinos and a muon [41]. The muonically decaying
$\PGt$ lepton is used as a probe to measure the identification rate for the
hadronically decaying $\PGt$ lepton. For each identification WP, there is one
nuisance parameter that alters the identification efficiency.

Three additional sets of uncertainties related to misidentification are
included. A 20% normalization uncertainty is applied to the contributions
arising from misidentified candidates in the CRs, where the loose WP is used.
To account for the misidentification uncertainties at tighter WPs, the
uncertainties from the FF derivation are propagated to background components
with misidentified candidates in the signal regions, where the medium or tight
WP is used. Systematic uncertainties in the simulated samples used in the
derivation of the FFs are small compared with the associated statistical
uncertainties, and are neglected. Each bin in which FFs are derived is
therefore assigned an independent uncertainty corresponding to the statistical
uncertainty in the FF derivation. There are 20 such uncertainties for the QCD
multijet FFs and nine for the FFs. Finally, a set of uncertainties accounts
for differences between the FFs in the determination regions and the main fit
region. The uncertainty is taken to be the relative difference between the two
FFs, as found in simulation. For the 2- QCD multijet FFs, where no simulated
events pass the selection criteria, an 80% uncertainty is applied. These
uncertainties are fully correlated within each region and across years.

A simplified Barlow–Beeston approach [66], with a single nuisance parameter
per bin of the distribution, is used to treat the statistical uncertainties
due to the limited size of the simulated samples.

Other experimental uncertainties that have a smaller impact on the analysis
include: uncertainties in the total delivered luminosity (normalization only),
the uncertainty in the number of inelastic $\Pp\Pp$ collisions in a single
bunch crossing, uncertainties related to detector conditions, and
uncertainties in the determination of the trigger efficiencies.

Theoretical uncertainties arising from the PDFs are assessed separately for
the , X, and signal processes. For each process, a single uncertainty is
constructed by adding in quadrature the contributions from all of the Hessian
eigenvectors of the process PDF set and its variations in the strong coupling
() [67]. Independent variations of the renormalization and factorization
scales by factors of two are used to assess the uncertainty from the
approximations used in the matrix element calculations; an envelope that
describes the largest impacts of these variations is formed and encoded as a
single nuisance parameter. Variations where the two scales differ by a factor
of four are not included. Uncertainties in the parton shower simulation are
assessed by changing the scale at which is evaluated in the shower by factors
of two up and down. This is done independently for the initial state and final
state radiation, creating two independent nuisance parameters.

An additional 20% normalization uncertainty is applied to the component. This
accounts for the fact that the component is not as well determined as the
overall cross section.

##  0.8 Fit setup and validation

The results of the analysis are extracted by performing a simultaneous binned
maximum likelihood fit over both years and all signal and control regions. The
SRs defined in Fig. 2 provide the main sensitivity to any potential signal,
while the CRs are used to provide estimates of some of the important
background processes from data. For the 0- region, the distribution of the
number of jets in the event ($\mathrm{N}_{\text{jets}}$) is used in the fit,
whereas in the 1- and 2- regions the distributions of the
$\mathrm{DNN}_{\ttbar}$ scores are fitted.

The fit is performed separately for each mass hypothesis to extract any
dependence of the result on the VLLmass. In each case, a single signal
strength parameter that multiplies the VLLproduction cross section is used.
Systematic uncertainties are treated as nuisance parameters in the fit.
Nuisance parameters that affect the overall rate of a process are described by
log-normal probability density functions, whereas nuisance parameters that
affect the shapes of distributions are described by Gaussian probability
density functions. Upper limits at the 95% confidence level () on the signal
strength are determined using the criterion [68, 69] with the profile
likelihood ratio as the test statistic. In this procedure, we make use of the
asymptotic approximation [70]. The signal significance is calculated from the
likelihoods of the background-only and signal-plus-background models using
Wilks’ theorem [71].

The analysis procedure and optimization is performed “blinded” (without
looking at the data) and a number of validation tests are also performed on
the fit setup without using the observed data in the main fit region. Fits are
performed on pseudo-data generated based on the expected background model, and
then applying Poisson-distributed smearing to the values in each bin. This
test is performed in four configurations: without signal, and also including
signal simulation with three different input cross sections. For each
configuration, 500 pseudo-experiments are performed. The signal strength
extracted from the fits is consistent with it being an unbiased estimator of
the true signal strength, and the uncertainties quoted on the signal strengths
are consistent with the expected frequentist coverage.

Fits are also performed on collision data in validation regions. These regions
are chosen to be orthogonal to the primary analysis regions, and with almost
no signal contamination. Separate validation regions are chosen for each of
the 0-, 1-, and 2- categories. In every case, the validation region is chosen
to have exactly two -tagged jets (as compared to $\geq 3$ in the main analysis
region). In the 1- and 2- categories, the validation regions also require that
the event score from the DNN classifier separating signal from is less than
0.8, in order to ensure that there is no signal contamination. The results of
these fits are consistent with the background-only hypothesis.

##  0.9 Results

Postfit distributions for the background-only and signal-plus-background fits
are shown in Fig. 3. For the fits including signal, the 600 VLL mass point is
shown as an example. The signal fit shows consistency across the sensitive
channels.

The observed data show excesses in the highest $\mathrm{DNN}_{\ttbar}$ bins
for both the 1- and 2- channels, for both 2017 and 2018. The 0- channel does
not show any excess, which is expected given its lower sensitivity.

Figure 3: Postfit plots for 2017 (top two rows) and 2018 (bottom two rows)
showing the fitted distributions in the signal region for the different
multiplicity channels (from left to right: 0-, 1-, 2-). The jet multiplicity
is fit for the 0- region, whereas the $\mathrm{DNN}_{\ttbar}$ score is fit for
the 1- and 2- regions. The first and third row show the background-only fits
and the second and fourth rows show the fit including the signal.

The expected and observed limits are shown in Fig. 4. The
$\mathrm{DNN}_{\ttbar}$ output, used in the 1- and 2- signal regions, is not
very sensitive to the signal mass, and therefore the results are highly
correlated across all mass points. The higher than expected limits are
consistent with the existence of a signal. As a result, none of the mass range
tested here (500–1050) is excluded at the 95% , and limits are set between 10
and 30fb, depending on the VLLmass hypothesis. At the representative 600mass
point, the signal-plus-background model is preferred over the background-only
model at the level of 2.8 standard deviations ($\sigma$). The highest
$\mathrm{DNN}_{\ttbar}$ bins of each region fit the data better when the
signal hypothesis is included, especially for the 2018 2- region, where the
background-dominated part of the fit is also substantially improved. Nuisance
parameters related to the fake factor estimates and the statistical
uncertainties of the simulated events in high $\mathrm{DNN}_{\ttbar}$ bins are
also pulled to lower likelihood values in the background-only fit,
contributing to the significance.

A number of post-unblinding checks were performed to test the robustness of
the result, reliability of the model, and consistency of the data with the
signal hypothesis. As indicated by the significance in the combined fit, the
signal interpretation is consistent across the channels. No excess is observed
in the significantly less sensitive 0- regions, whereas the excesses observed
in the 1- and 2- regions are consistent with each other. The observed limits
for each of the 1- and 2- channels, when fit individually, are above the
expected limits. Taking the 600mass point as an example, the significances for
these channels when fit individually range from 0.8 to 2.0$\sigma$ and are all
within 1$\sigma$ of their expected significances for a signal strength
corresponding to the best combined fit signal strength.

Several alternative analysis procedures were also tested and found to give
similar results to those presented here. These included rejecting candidates
with only two reconstructed tracks, applying a reweighting to account for
possible mismodelling of the spectrum, increasing the size of certain
experimental uncertainties, and including an additional uncertainty related to
the difference between data and simulation in the FFs. All of these were found
to give very similar results to the main result presented here.

These results have been obtained considering only electroweak production of
VLL pairs. However, tests showed that the results hold when including
production of VLLpairs as well.

Figure 4: Expected and observed 95% upper limits on the product of the VLLpair
production cross section and the branching fraction to third generation quarks
and leptons, combining the 2017 and 2018 data and all multiplicity channels.
The theoretical prediction in the 4321 model for electroweak production of
VLLs is also shown.

##  0.10 Summary

The first search for vector-like leptons (VLLs) in the context of the 4321
model has been presented, using proton-proton collision data collected with
the CMS detector at $\sqrt{s}=13\TeV$, corresponding to an integrated
luminosity of 96.5. The probed model consists of an extension of the standard
model with an
$\text{SU}(4)\times\text{SU}(3)^{\prime}\times\text{SU}(2)_{L}\times\text{U}(1)^{\prime}$
gauge sector that can provide a combined explanation for multiple anomalies
observed in hadron decays, which point to lepton flavour nonuniversality. In
the model, a leptoquark is predicted as the primary source of lepton flavour
nonuniversality, while the ultraviolet-completion predicts additional vector-
like fermion families. In particular, VLLs are investigated via their
couplings to standard model fermions through leptoquark interactions,
resulting in third-generation fermion signatures. Final states containing at
least three -tagged jets and varying $\PGt$ lepton multiplicities are
considered. To improve the search sensitivity, graph neural networks are
trained to discriminate between classes of events using the kinematic
relationships between particles in an environment with high jet multiplicity.

An excess of 2.8 standard deviations over the standard model background-only
hypothesis was observed in the data at the representative VLLmass point of
600. As a result, no VLLmasses are excluded at the 95% confidence level and
limits on the product of the VLL pair production cross section and their
branching ratio to third generation fermions are set between 10 and 30fb,
depending on the VLLmass hypothesis.

## References

  * [1] L. Di Luzio, A. Greljo, and M. Nardecchia, “Gauge leptoquark as the origin of B-physics anomalies”,  Phys. Rev. D 96 (2017) 115011, 10.1103/physrevd.96.115011, arXiv:1708.08450.
  * [2] L. Di Luzio et al., “Maximal flavour violation: a Cabibbo mechanism for leptoquarks”,  JHEP 11 (2018) 081, 10.1007/jhep11(2018)081, arXiv:1808.00942.
  * [3] A. Greljo and B. A. Stefanek, “Third family quark–lepton unification at the TeV scale”,  Phys. Lett. B 782 (2018) 131, 10.1016/j.physletb.2018.05.033, arXiv:1802.04274.
  * [4] 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.
  * [5] CMS Collaboration, “CMS luminosity measurement for the 2017 data-taking period at $\sqrt{s}=13~{}\mathrm{TeV}$”, cms physics analysis summary, CERN, 2018.
  * [6] CMS Collaboration, “CMS luminosity measurement for the 2018 data-taking period at $\sqrt{s}=13~{}\mathrm{TeV}$”, cms physics analysis summary, CERN, 2019.
  * [7] BaBar Collaboration, “Evidence for an excess of $\bar{B}\to D^{(*)}\tau^{-}\bar{\nu}_{\tau}$ decays”,  Phys. Rev. Lett. 109 (2012) 101802, 10.1103/PhysRevLett.109.101802, arXiv:1205.5442.
  * [8] BaBar Collaboration, “Measurement of an excess of $\bar{B}\to D^{(*)}\tau^{-}\bar{\nu}_{\tau}$ decays and implications for charged Higgs bosons”,  Phys. Rev. D 88 (2013) 072012, 10.1103/PhysRevD.88.072012, arXiv:1303.0571.
  * [9] LHCb Collaboration, “Measurement of the ratio of branching fractions $\mathcal{B}(\bar{B}^{0}\to D^{*+}\tau^{-}\bar{\nu}_{\tau})/\mathcal{B}(\bar{B}^{0}\to D^{*+}\mu^{-}\bar{\nu}_{\mu})$”,  Phys. Rev. Lett. 115 (2015) 111803, 10.1103/PhysRevLett.115.111803, arXiv:1506.08614. [Erratum: https://doi.org/10.1103/PhysRevLett.115.159901].
  * [10] Belle Collaboration, “Measurement of the branching ratio of $\bar{B}\to D^{(\ast)}\tau^{-}\bar{\nu}_{\tau}$ relative to $\bar{B}\to D^{(\ast)}\ell^{-}\bar{\nu}_{\ell}$ decays with hadronic tagging at Belle”,  Phys. Rev. D 92 (2015) 072014, 10.1103/PhysRevD.92.072014, arXiv:1507.03233.
  * [11] Belle Collaboration, “Measurement of the $\tau$ lepton polarization and $R(D^{*})$ in the decay $\bar{B}\to D^{*}\tau^{-}\bar{\nu}_{\tau}$”,  Phys. Rev. Lett. 118 (2017) 211801, 10.1103/PhysRevLett.118.211801, arXiv:1612.00529.
  * [12] LHCb Collaboration, “Test of lepton flavor universality by the measurement of the $B^{0}\to D^{*-}\tau^{+}\nu_{\tau}$ branching fraction using three-prong $\tau$ decays”,  Phys. Rev. D 97 (2018) 072013, 10.1103/PhysRevD.97.072013, arXiv:1711.02505.
  * [13] Belle Collaboration, “Measurement of the $\tau$ lepton polarization and $R(D^{*})$ in the decay $\bar{B}\rightarrow D^{*}\tau^{-}\bar{\nu}_{\tau}$ with one-prong hadronic $\tau$ decays at Belle”,  Phys. Rev. D 97 (2018) 012004, 10.1103/PhysRevD.97.012004, arXiv:1709.00129.
  * [14] LHCb Collaboration, “Measurement of the ratio of the $B^{0}\to D^{*-}\tau^{+}\nu_{\tau}$ and $B^{0}\to D^{*-}\mu^{+}\nu_{\mu}$ branching fractions using three-prong $\tau$-lepton decays”,  Phys. Rev. Lett. 120 (2018) 171802, 10.1103/PhysRevLett.120.171802, arXiv:1708.08856.
  * [15] Belle Collaboration, “Measurement of $\mathcal{R}(D)$ and $\mathcal{R}(D^{*})$ with a semileptonic tagging method”,  Phys. Rev. Lett. 124 (2020) 161803, 10.1103/PhysRevLett.124.161803, arXiv:1910.05864.
  * [16] LHCb Collaboration, “Test of lepton universality in beauty-quark decays”,  Nature Phys. 18 (2022) 277, 10.1038/s41567-021-01478-8, arXiv:2103.11769.
  * [17] M. Blanke and A. Crivellin, “$B$ meson anomalies in a Pati-Salam model within the Randall-Sundrum background”,  Phys. Rev. Lett. 121 (2018) 011801, 10.1103/PhysRevLett.121.011801, arXiv:1801.07256.
  * [18] L. Calibbi, A. Crivellin, and T. Li, “Model of vector leptoquarks in view of the $B$-physics anomalies”,  Phys. Rev. D 98 (2018) 115002, 10.1103/PhysRevD.98.115002, arXiv:1709.00692.
  * [19] CMS Collaboration, “Search for vector-like leptons in multilepton final states in proton-proton collisions at $\sqrt{s}$ = 13 TeV”,  Phys. Rev. D 100 (2019) 052003, 10.1103/PhysRevD.100.052003, arXiv:1905.10853.
  * [20] F. Abe, D. Amidei, G. Apollinari, and M. Atac, “Measurement of the ratio B(W$\rightarrow$$\tau$$\nu$)/B(W$\rightarrow$e$\nu$) in $p\bar{p}$ collisions at $\surd$s =1.8 TeV”,  Phys. Rev. Lett. 68 (1992) 3398, 10.1103/PhysRevLett.68.3398.
  * [21] D0 Collaboration, “A measurement of the $W\to\tau\nu$ production cross section in $p\bar{p}$ collisions at $\sqrt{s}=1.8$ TeV”,  Phys. Rev. Lett. 84 (2000) 5710, 10.1103/PhysRevLett.84.5710, arXiv:hep-ex/9912065.
  * [22] LHCb Collaboration, “Measurement of forward $W\to e\nu$ production in $pp$ collisions at $\sqrt{s}=8$ TeV”,  JHEP 10 (2016) 030, 10.1007/JHEP10(2016)030, arXiv:1608.01484.
  * [23] ATLAS Collaboration, “Test of the universality of $\tau$ and $\mu$ lepton couplings in $W$-boson decays with the ATLAS detector”,  Nature Phys. 17 (2021) 813, 10.1038/s41567-021-01236-w, arXiv:2007.14040.
  * [24] C. Cornella et al., “Reading the footprints of the B-meson flavor anomalies”,  JHEP 08 (2021) 050, 10.1007/jhep08(2021)050, arXiv:2103.16558.
  * [25] T. Hastie, R. Tibshirani, and J. Friedman, “The elements of statistical learning”.  Springer-Verlag, second edition, 2013.  10.1007/978-0-387-84858-7, ISBN 978-0-387-84858-7.
  * [26] V. Mikuni and F. Canelli, “ABCNet: An attention-based method for particle tagging”,  Eur. Phys. J. Plus 135 (2020) 463, 10.1140/epjp/s13360-020-00497-3, arXiv:2001.05311.
  * [27] G. Kasieczka, B. Nachman, M. D. Schwartz, and D. Shih, “Automating the ABCD method with machine learning”,  Phys. Rev. D 103 (2021) 035021, 10.1103/PhysRevD.103.035021, arXiv:2007.14400.
  * [28] HEPData record for this analysis, 2022.  10.17182/hepdata.132482.
  * [29] CMS Collaboration, “The CMS experiment at the CERN LHC”,  JINST 3 (2008) S08004, 10.1088/1748-0221/3/08/S08004.
  * [30] 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.
  * [31] CMS Collaboration, “The CMS trigger system”,  JINST 12 (2017) P01020, 10.1088/1748-0221/12/01/P01020, arXiv:1609.02366.
  * [32] CMS Collaboration, “Technical proposal for the Phase-II upgrade of the Compact Muon Solenoid”, CMS Technical Proposal CERN-LHCC-2015-010, CMS-TDR-15-02, 2015.
  * [33] 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.
  * [34] M. Cacciari, G. P. Salam, and G. Soyez, “The anti-$\kt$ jet clustering algorithm”,  JHEP 04 (2008) 063, 10.1088/1126-6708/2008/04/063, arXiv:0802.1189.
  * [35] 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.
  * [36] CMS Collaboration, “Pileup mitigation at CMS in 13 TeV data”,  JINST 15 (2020) P09018, 10.1088/1748-0221/15/09/P09018, arXiv:2003.00503.
  * [37] 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.
  * [38] CMS Collaboration, “Identification of heavy-flavour jets with the CMS detector in pp collisions at 13 TeV”,  JINST 13 (2018) P05011, 10.1088/1748-0221/13/05/P05011, arXiv:1712.07158.
  * [39] E. Bols et al., “Jet flavour classification using DeepJet”,  JINST 15 (2020) P12012, 10.1088/1748-0221/15/12/p12012, arXiv:2008.10519.
  * [40] 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.
  * [41] 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.
  * [42] CMS Collaboration, “Identification of hadronic tau lepton decays using a deep neural network”,  JINST 17 (2022) P07023, 10.1088/1748-0221/17/07/P07023, arXiv:2201.08458.
  * [43] 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.
  * [44] CMS Collaboration, “CMS technical design report for the pixel detector upgrade”, CMS Technical Design Report CERN-LHCC-2012-016, CMS-TDR-011, 2012.
  * [45] CMS Collaboration, “Performance of b tagging algorithms in proton-proton collisions at 13 TeV with Phase 1 CMS detector”, CMS Detector Performance Note CMS-DP-2018-033, 2018.
  * [46] 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.
  * [47] A. Alloul et al., “FeynRules 2.0 – a complete toolbox for tree-level phenomenology”,  Comput. Phys. Commun. 185 (2014) 2250, 10.1016/j.cpc.2014.04.012, arXiv:1310.1921.
  * [48] 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.
  * [49] 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.
  * [50] R. D. Ball et al., “Parton distributions from high-precision collider data”,  Eur. Phys. J. C 77 (2017) 663, 10.1140/epjc/s10052-017-5199-5, arXiv:1706.00428.
  * [51] 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.
  * [52] 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.
  * [53] 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.
  * [54] M. Beneke, P. Falgari, S. Klein, and C. Schwinn, “Hadronic top-quark pair production with NNLL threshold resummation”,  Nucl. Phys. B 855 (2012) 695, 10.1016/j.nuclphysb.2011.10.021, arXiv:1109.1536v1.
  * [55] M. Cacciari et al., “Top-pair production at hadron colliders with next-to-next-to-leading logarithmic soft-gluon resummation”,  Phys. Lett. B 710 (2012) 612, 10.1016/j.physletb.2012.03.013, arXiv:1111.5869.
  * [56] P. Bärnreuther, M. Czakon, and A. Mitov, “Percent-level-precision physics at the Tevatron: Next-to-next-to-leading order QCD corrections to $q\bar{q}\to t\bar{t}+X$”,  Phys. Rev. Lett. 109 (2012) 132001, 10.1103/physrevlett.109.132001, arXiv:1204.5201.
  * [57] M. Czakon and A. Mitov, “NNLO corrections to top-pair production at hadron colliders: the all-fermionic scattering channels”,  JHEP 12 (2012) 054, 10.1007/jhep12(2012)054, arXiv:1207.0236.
  * [58] M. Czakon and A. Mitov, “NNLO corrections to top pair production at hadron colliders: the quark-gluon reaction”,  JHEP 01 (2013) 080, 10.1007/jhep01(2013)080, arXiv:1210.6832.
  * [59] M. Czakon, P. Fiedler, and A. Mitov, “Total top-quark pair-production cross section at hadron colliders through O($\alpha_{S}^{4}$)”,  Phys. Rev. Lett. 110 (2013) 252004, 10.1103/physrevlett.110.252004, arXiv:1303.6254.
  * [60] 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.
  * [61] CMS Collaboration, “Measurements of $\mathrm{t\bar{t}}$ cross sections in association with b jets and inclusive jets and their ratio using dilepton final states in pp collisions at $\sqrt{s}=$ 13 TeV”,  Phys. Lett. B 776 (2018) 355, 10.1016/j.physletb.2017.11.043, arXiv:1705.10141.
  * [62] ATLAS Collaboration, “Measurements of inclusive and differential fiducial cross-sections of production with additional heavy-flavour jets in proton-proton collisions at $\sqrt{s}$ = 13 TeV with the ATLAS detector”,  JHEP 04 (2019) 046, 10.1007/jhep04(2019)046, arXiv:1811.12113.
  * [63] CMS Collaboration, “Measurement of the $\mathrm{t\bar{t}}\mathrm{b\bar{b}}$ production cross section in the all-jet final state in pp collisions at $\sqrt{s}=$ 13 TeV”,  Phys. Lett. B 803 (2020) 135285, 10.1016/j.physletb.2020.135285, arXiv:1909.05306.
  * [64] CMS Collaboration, “Measurement of the cross section for $\mathrm{t\bar{t}}$ production with additional jets and b jets in pp collisions at $\sqrt{s}=$ 13 TeV”,  JHEP 07 (2020) 125, 10.1007/JHEP07(2020)125, arXiv:2003.06467.
  * [65] GEANT4 Collaboration, “—a simulation toolkit”,  Nucl. Instrum. Meth. A 506 (2003) 250, 10.1016/S0168-9002(03)01368-8.
  * [66] R. J. Barlow and C. Beeston, “Fitting using finite Monte Carlo samples”,  Comput. Phys. Commun. 77 (1993) 219, 10.1016/0010-4655(93)90005-W.
  * [67] J. Butterworth et al., “PDF4LHC recommendations for LHC Run II”,  J. Phys. G 43 (2016) 023001, 10.1088/0954-3899/43/2/023001, arXiv:1510.03865.
  * [68] 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.
  * [69] A. L. Read, “Presentation of search results: The CL${}_{\text{s}}$ technique”,  J. Phys. G 28 (2002) 2693, 10.1088/0954-3899/28/10/313.
  * [70] 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].
  * [71] S. S. Wilks, “The large-sample distribution of the likelihood ratio for testing composite hypotheses”,  Ann. Math. Statist. 9 (1938) 60, 10.1214/aoms/1177732360.