First measurement of the cross section for top quark pair production with additional charm jets using dileptonic final states in pp collisions at $\sqrt{s} = $ 13 TeV

The first measurement of the inclusive cross section for top quark pairs ($\mathrm{t\bar{t}}$) produced in association with two additional charm jets is presented. The analysis uses the dileptonic final states of $\mathrm{t\bar{t}}$ events produced in proton-proton collisions at a centre-of-mass energy of 13 TeV. The data correspond to an integrated luminosity of 41.5 fb$^{-1}$, recorded by the CMS experiment at the LHC. A new charm jet identification algorithm provides input to a neural network that is trained to distinguish among $\mathrm{t\bar{t}}$ events with two additional charm ($\mathrm{t\bar{t}c\bar{c}}$), bottom ($\mathrm{t\bar{t}b\bar{b}}$), and light-flavour or gluon ($\mathrm{t\bar{t}LL}$) jets. By means of a template fitting procedure, the inclusive $\mathrm{t\bar{t}c\bar{c}}$, $\mathrm{t\bar{t}b\bar{b}}$, and $\mathrm{t\bar{t}LL}$ cross sections are simultaneously measured, together with their ratios to the inclusive $\mathrm{t\bar{t}}$ + two jets cross section. This provides measurements of the $\mathrm{t\bar{t}c\bar{c}}$ and $\mathrm{t\bar{t}b\bar{b}}$ cross sections of 10.1 $\pm$ 1.2 (stat) $\pm$ 1.4 (syst) pb and 4.54 $\pm$ 0.35 (stat) $\pm$ 0.56 (syst) pb, respectively, in the full phase space. The results are compared and found to be consistent with predictions from two different matrix element generators with next-to-leading order accuracy in quantum chromodynamics, interfaced with a parton shower simulation.


Introduction
The modelling of top quark pair (tt) production in association with jets from the hadronization of bottom (b) or charm (c) quarks, referred to as b jets and c jets, respectively, in proton-proton (pp) collisions at the CERN LHC is challenging. Calculations of the production cross section for top quark pairs with additional pairs of b jets (ttbb) are available at next-to-leading order (NLO) in quantum chromodynamics (QCD) [1][2][3][4][5], but suffer from large uncertainties due to the choice of factorization (µ F ) and renormalization (µ R ) scales. The uncertainties arise from the different energy (or mass) scales in ttbb production, which range from the large scales associated with the top quark mass (m t ), to the relatively small scales associated with additional jets resulting mostly from gluon splitting into bb pairs. To better understand the ttbb process, the ATLAS and CMS experiments have conducted several measurements in pp collisions at centre-of-mass energies of 7, 8, and 13 TeV [6-13].
The production of a tt pair with an additional pair of c jets (ttcc) has received far less attention, both theoretically and experimentally. Whereas the experimental signature of a b jet looks quite different from that of a light-flavour (LF) or gluon jet, the differences are much less pronounced for c jets. This explains the challenge of simultaneously separating ttcc and ttbb events from a large background of tt events with additional LF or gluon jets (ttLL). With the development of a charm jet identification algorithm ("c tagger") [14], these signatures can now be more efficiently disentangled. We present the first measurement of the inclusive ttcc cross section and its ratio to the inclusive tt + two jets (ttjj) cross section. A fully consistent treatment of the different additional jet flavours is ensured using a technique that simultaneously extracts the ttcc, ttbb, and ttLL cross sections. The measurement is performed in the dileptonic decay channel of the tt events using a data sample of pp collisions at a centre-of-mass energy of 13 TeV, collected with the CMS detector in 2017, corresponding to an integrated luminosity of 41.5 fb −1 [15]. This data set benefits from the upgrade of the pixel tracking detector [16], which was installed in winter 2016-2017 and which has been shown to significantly improve the performance of heavy-flavour (HF) jet identification [17].
Although the additional b jets in ttbb events are predominantly produced via gluon splitting into bb pairs, they can also originate from the decay of a Higgs boson (H). Previous measurements of Higgs boson production in association with a top quark pair, where the Higgs boson decays into a pair of b quarks (ttH, H → bb) [18][19][20][21], suffer from a nonresonant background of gluon-induced ttbb events, and to a lesser extent also from ttcc events due to the misidentification of c jets as b jets. The techniques described here provide a basis for a simultaneous measurement of the ttbb and ttcc processes from data that can be adopted in future ttH analyses to significantly reduce the uncertainties related to these backgrounds.

The CMS detector
The central feature of the CMS apparatus is a superconducting solenoid of 6 m internal diameter, providing a magnetic field of 3.8 T. 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. The silicon tracker measures charged particles within the pseudorapidity range |η| < 2.5. During the LHC running period when the data used for this analysis were recorded, the silicon tracker consisted of 1856 silicon pixel and 15 148 silicon strip detector modules. For nonisolated particles of 1 < p T < 10 GeV and |η| < 2.5, the track resolutions are typically 1.5% in p T and 20-75 µm in the transverse impact parameter [22]. Forward calorimeters extend the η coverage

Signal definition
A typical Feynman diagram describing the dileptonic ttcc process is shown in Fig. 1. To provide an unambiguous interpretation of the results, a generator-level definition of the event categories is needed, based on the flavours of the additional jets. The HF jets are identified at the generator level using the procedure of ghost matching [38], where in addition to the reconstructed final-state particles, the generated b and c hadrons are clustered into the jets. However, the modulus of the hadron four-momentum is set to a small number to prevent these generated hadrons from affecting the reconstructed jet momentum and to ensure that only their directional information is retained. Jets that contain both b and c hadrons are labelled as b jets, since the c hadrons most likely originate from a decay of a b hadron. The results of this measurement are reported in terms of the fiducial and full phase spaces. Fiducial phase space: In the definition of the fiducial phase space, all of the final-state particles (except for the neutrinos) resulting from the decay chain: pp → ttjj → + ν b − ν bjj are required to be within the region of the detector in which these objects can be properly reconstructed. The fiducial phase space is therefore defined by the presence of two oppositely charged electrons, muons, or τ leptons ( ) at the generator level with p T > 25 GeV and |η| < 2.4. Each lepton is required to originate from a decay of a W boson, which in turn results from a top quark decay. Particle-level jets are defined by clustering generated final-state particles with a mean lifetime greater than 30 ps, excluding neutrinos, using the anti-k T algorithm [39,40] with a distance parameter of 0.4. We demand two particle-level b jets from the top quark decays with p T > 20 GeV and |η| < 2.4. Besides these two b jets, at least two additional particlelevel jets must be present, with the same kinematic requirements imposed. Jets that lie within ∆R = √ (∆η) 2 + (∆φ) 2 < 0.4 of either one of the leptons from the W boson decays, where φ denotes the azimuthal angle, are excluded. The above requirements define the inclusive ttjj signal, which is then subdivided based on the flavour content of the additional particle-level jets (not from the top quark decays) into the following categories: ttbb: At least two additional b jets are present, each containing at least one b hadron. ttbL: Only one additional b jet is present, containing at least one b hadron. In addition, at least one additional LF or c jet is present. This category results from ttbb events in which one of the two additional b jets is outside the acceptance, or two b jets are merged into one.
ttcc: No additional b jets are present, but at least two additional c jets are found, each containing at least one c hadron.
ttcL: No additional b jets are present, and only one additional c jet is found, containing at least one c hadron. In addition, at least one additional LF jet is present. This category results from ttcc events in which one of the two additional c jets is outside the acceptance, or two c jets are merged into one.
ttLL: No additional b or c jets are present, but at least two additional LF jets are within the acceptance.
All other tt events that do not fit in any of the above categories, because they do not fulfill the acceptance requirements described in the definition of the fiducial phase space, are labelled as "tt+other". These could for example be events with semileptonic or fully hadronic tt decays that pass the event selection criteria outlined in Section 6, or dileptonic tt events for which the leptons or the particle-level jets as described above are not within the fiducial detector volume. These contributions are estimated from the same simulations as those used for the signal events.
Full phase space: The definition of the full phase space comprises dileptonic, semileptonic, and fully hadronic tt decays that contain in addition at least two particle-level jets with p T > 20 GeV and |η| < 2.4. These jets must not originate from the decays of the top quarks or W bosons. There are no requirements imposed on the generator-level leptons or on the particle-level jets that result from the top quark and W boson decays. The measurement in the fiducial phase space is extrapolated to the full phase space by applying an acceptance factor, estimated from simulations, to each signal category based on the additional jet flavours.
The definition of the fiducial phase space is much closer to the reconstructed phase space in which the measurement is performed, and therefore expected to suffers less from theoretical uncertainties that affect the extrapolation from the fiducial to the full phase space. However, the full phase space definition is better suited for comparison with theoretical calculations.

Object reconstruction
The global event reconstruction is based on the particle-flow algorithm [41], which aims to reconstruct and identify each individual particle in an event by combining 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 zero-suppression effects and 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 these reconstructed particles using the infrared and collinear safe anti-k T algorithm 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 simulations to be, on average, within 5 to 10% of the true momentum over the whole p T spectrum and detector acceptance. Pileup interactions can contribute additional tracks and calorimetric energy depositions, increasing the apparent jet momentum. To mitigate this effect, tracks identified as originating from pileup vertices are discarded and an offset correction is applied to correct for remaining contributions. Jet energy corrections are derived from simulation studies so that the average measured response of jets becomes identical to that of particle-level jets. In situ measurements of the momentum balance in dijet, photon+jets, Z + jets, and multijet events are used to determine any residual differences between the jet energy scale in data and simulations [42], and appropriate corrections are made. Additional selection criteria are applied to each jet to remove jets potentially dominated by instrumental effects or reconstruction failures.
The missing transverse momentum vector p miss T is computed as the negative of the vector sum of the p T of all the particle-flow candidates in an event, and its magnitude is denoted as p miss T [43]. The p miss T is modified to account for corrections to the energy scale of the reconstructed jets in the event.
The candidate vertex with the largest value of summed physics-object p 2 T is taken to be the primary pp interaction vertex. The physics objects are the jets and the associated p miss T . Jets originating from the hadronization of b and c quarks are identified using the deep combined secondary vertex (DeepCSV) algorithm [14], which uses information on the decay vertices of long-lived mesons and the impact parameters of the charged particle tracks as input to a deep neural network (NN) classifier. For the identification of b jets, a medium working point is chosen, corresponding to a ≈70% efficiency for correctly selecting b jets and a misidentification probability for LF (c) jets of ≈1 (12)%, derived from simulated tt events. The same algorithm also provides discriminators to distinguish c jets from LF and b jets, which are collectively referred to as c tagging discriminators. The information from the distributions of these observables plays a key role in this analysis, and the calibration of the c tagger is further discussed in Section 8.

Event selection
An event selection has been employed to select a subset of events that consists almost exclusively (more than 95% as evaluated from simulations) of dileptonic tt events with at least two additional jets. Exactly two reconstructed, oppositely charged leptons (either electrons or muons) are required to be present. This procedure also selects τ leptons that decay into an electron or a muon. The electrons and muons are required to have p T > 25 GeV and |η| < 2.4, and should be isolated from other objects in the event. At least four jets are required in the event, all of which must be spatially separated from the isolated leptons by imposing ∆R( , jet) > 0.5. Only jets with p T > 30 GeV and |η| < 2.4 are considered. An assignment of the jets to the expected partons is made to identify the b jets from top quark decays and jets originating from additional radiation (described in detail in Section 7). The two jets assigned to the b quarks from the top quark decays are required to be b tagged. Neutrinos from leptonically decaying W bosons are not detected, but instead contribute to the p miss T , which is required to exceed 30 GeV in events with two electrons (ee) or two muons (µµ), in order to reduce contributions from DY events. In events with one electron and one muon (eµ), no requirement is imposed on p miss T . In order to further reduce the contribution from DY production in ee and µµ events, the invariant mass of the two leptons (m ) is required to be outside of the Z boson mass window, m ∈ [m Z − 15, m Z + 15] GeV, with m Z = 91.2 GeV [44]. For all events, it is required that m > 12 GeV in order to minimize contributions from low-mass resonances.

Matching jets to partons
The distinction among the ttcc, ttbb, and ttLL categories relies on the correct identification of the additional jets not coming from the decay of the top quarks. Assuming a 100% branching fraction for the decay t → bW, and focusing on the dileptonic decay channel, two b jets are expected from the top quark decays and at least two additional jets are required through the other event selection criteria. In practice, not all b jets from top quark decays will be reconstructed within the acceptance of the detector and additional jets will also not necessarily pass the reconstruction and selection criteria. In this section, a matching procedure is described to achieve the most accurate correspondence between the final-state jets and the expected partons. This is done by considering all possible permutations of four jets in the collection of jets passing the selection criteria described in Section 6 and training a NN to identify the correct jet-parton assignment.
Whether a given permutation corresponds to a correct jet-parton assignment can be inferred from different quantities, such as the jet kinematic variables, b and c tagging discriminators, and angular separations and invariant masses between pairs of jets (and between jets and leptons). The NN processes, for each event, all possible jet-parton assignments and is trained on the aforementioned observables to give the highest possible score to the correct permutation. For the ttLL, ttcc, and ttcL categories, the b and c tagging discriminators dominate the decision made by the NN. For the ttbb category all four partons have the same flavour, so the correct assignment of the four b jets relies mainly on the angular separations and invariant masses between jets or leptons. The best assignment for the ttbL category benefits from a combination of all these observables.
The main objective is to identify additional HF jets in the event. In the assignment of b jets from the tt decays, it does not matter which b jet is matched to the decay of the top quark or antiquark. Any permutation for which these b jets are reversed can still be considered appropriate for this measurement. If at least one additional b or c jet is present, a correct permutation has to identify these as the first or second additional jets. The NN is also trained to choose a permutation in which the first additional jet has a larger b tagging discriminator value than the second additional jet. With these considerations in mind, three output classes are defined. One of the NN outputs denotes the probability for a given permutation to correspond to the correct jet-parton assignment (P + ). Another output class refers to those permutations for which the additional jets are correctly matched, but the b jets from the tt decays are reversed as explained above (P × ). The third output represents all permutations for which the matching is wrong (P − ), meaning that either at least one of the b jets from the tt decays is not correctly matched, or an additional HF jet is found but is not identified as either the first or the second additional jet. The best jet-parton assignment is then identified by selecting the permutation with the highest value of: The NN is trained using the KERAS deep learning library [45], interfaced to TENSORFLOW [46] as a back end. Its architecture is composed of two fully connected hidden layers with 50 neurons each, and with a rectified linear unit activation function. The training is performed using an independent data set of simulated tt events that pass all the event selections outlined in Section 6, except for the b tagging requirement. Only events for which the two generatorlevel b quarks from top quark decays lie within ∆R < 0.3 of a reconstructed b jet are used for the training. These constitute ≈76% of the simulated tt events with at least two additional jets. Of those, the NN correctly identifies the two additional c (b) jets in 50 (30)% of the cases for ttcc (ttbb) events. For ttbb events, the matching is more challenging because the HF tagging information cannot help in separating additional b jets from those originating from top quark decays. A comparison between data and simulations of the NN score for the best permutation in each event is shown in Fig. 2, indicating good agreement between the two.   Figure 2: Comparison between data (points) and simulated predictions (histograms) for the distribution of the NN score for the best permutations of jet-parton assignments found in each event. Underflow is included in the first bin. The distributions are found after the event selections outlined in Section 6, but before fitting the predicted signal yields to the data. The lower panel shows the ratio of the yields in data to those predicted in simulations. The vertical bars represent the statistical uncertainties in data, while the hatched bands show the statistical uncertainty in the simulated predictions.

Charm jet identification and calibration
The DeepCSV HF tagging algorithm has a multiclass output structure that predicts the probabilities for each jet to contain a single b hadron (P(b)), two b hadrons (P(bb)), one or more c hadrons (P(c)), or no b or c hadrons (P(udsg)). The b tagger used throughout this analysis is constructed from the sum P(b) + P(bb). However, different combinations of output values provide different types of discrimination. Since the displacements of tracks and secondary vertices for c jets are on average smaller than those for b jets and larger than those for LF jets, a charm jet identification algorithm uses a combination of two discriminators. The first is used to distinguish c jets from LF jets (CvsL) and the second separates c jets from b jets (CvsB). The CvsL and CvsB discriminators are defined from the multiclass output structure of the DeepCSV algorithm as: The c tagging discriminators require calibration with the data, given that these algorithms are trained on simulated events and are therefore prone to mismodelling effects in the input variables. The observed initial discrepancies between the data and the simulated predictions can be as large as 50%, as can be seen from the distributions before calibration in Appendix A. In order to use these c tagger discriminators for the fit of the template distributions (as discussed in Section 9), the shape of the two-dimensional (CvsL, CvsB) distribution is corrected to reproduce the distribution observed in data. To this end, a novel calibration technique is employed that uses three control regions selecting for semileptonic tt, W+ c, and DY + jets events, which are enriched in b, c, and LF jets, respectively [47]. By means of an iterative fit in these three control regions, a set of scale factors for each jet flavour is derived, as a function of both the CvsL and CvsB discriminator values of a given jet. The effectiveness of this calibration for different jet flavours has been validated in the three corresponding control regions outlined above, showing that the calibrated distributions in simulations indeed match those in data within the associated uncertainties. Additionally, the method shows no bias when applied to pseudo-data constructed from simulated events that have artificial scale factors applied to them.
After applying this calibration to the simulated events, the distributions of the CvsL and CvsB discriminators provide a good description of the data, as shown in Fig. 3 for the first (upper row) and second (lower row) additional jet. The uncertainties related to this calibration are further discussed in Section 10.

Fit to an event-based neural network discriminator
The extraction of the ttcc, ttbb, and ttLL cross sections proceeds by means of a template fit to an observable that can distinguish among the different flavour categories. Since the differentiation relies mainly on the additional-jet flavour, the c tagging discriminators of the first and second additional jets provide a natural choice for separating the different signals. The CvsL discriminator distinguishes the ttcc and ttbb from the ttLL events, whereas the CvsB discriminator provides additional information that can be used to distinguish between the ttcc and ttbb events. The flavour tagging information for the second additional jet allows the ttcL and ttbL processes to be identified. Additional information is extracted from two kinematic variables. The first is the angular separation ∆R between the two additional jets. In the ttcc and ttbb processes, the additional jets arise predominantly from gluon splitting into cc and bb pairs, respectively, and are therefore expected to be more collimated. The second is the NN score for the best permutation (shown in Fig. 2), which is expected to be larger on average for ttLL events because the additional jets are well distinguished from the b jets from top quark decays. This observable indirectly incorporates information on the event kinematic features through its input variables.
Using the six aforementioned observables, a NN is trained. Given the relatively small number of inputs, an architecture for the NN is chosen with one hidden layer that comprises 30 neurons with a rectified linear unit activation function. This NN predicts output probabilities for five output classes: P(ttcc), P(ttcL), P(ttbb), P(ttbL), and P(ttLL). To obtain the distributions that are used in the fit, these probabilities are projected onto a two-dimensional phase space spanned by two derived discriminators: These discriminators can be interpreted as topology-specific c tagger discriminators that augment the information on the jet flavour of the two additional jets with additional event kinematic features to optimally distinguish different signal categories. The two-dimensional distributions of these discriminators, normalized to unit area, are shown in Fig. 4 for simulated Single t Rare CvsL discriminator first add. jet Single t Rare CvsB discriminator first add. jet Single t Rare CvsL discriminator second add. jet Single t Rare CvsB discriminator second add. jet : Comparison between data (points) and simulated predictions (histograms) for the CvsL (left column) and CvsB (right column) c tagging discriminator distributions of the first (upper row) and second (lower row) additional jet after applying the c tagging calibration. The distributions are found after the event selections outlined in Section 6, but before fitting the predicted signal yields to the data. The lower panels show the ratio of the yields in data to those predicted in simulations. The vertical bars represent the statistical uncertainties in data, while the hatched bands show the systematic uncertainty from the c tagging calibration only. dileptonic tt events. The different signal categories occupy different parts of this phase space, demonstrating that a fit of templates derived from these distributions to the data can be used to extract the ttcc, ttbb, and ttLL cross sections. It can be seen that the lower right corner of this phase space is almost exclusively populated with ttbb events, whereas the upper right corner is dominated by ttcc events.
Templates are constructed separately for each dilepton channel (ee, µµ, and eµ) and are fitted simultaneously to data by means of a binned maximum likelihood fit assuming Poisson statistics. The negative logarithm of the likelihood is minimized using the "Combine" framework developed for the combined Higgs boson measurements performed by ATLAS and CMS [48,49]. Uncertainties are included as nuisance parameters in the definition of the likelihood.
A first fit is performed to extract the absolute cross section (σ) of the ttcc, ttbb, and ttLL events in the fiducial phase space, where the expected yield in each bin is parametrized using the Here, L int denotes the total integrated luminosity and f norm i represents a given bin (with index i) of the simulated two-dimensional template, normalized to unit area. The measurements in the reconstructed phase space are corrected to the fiducial phase space through an efficiency ( ). To extract the result in the full phase space, an additional acceptance factor, A, is applied to the extracted fiducial cross section to account for the difference in acceptance between the fiducial and full phase spaces. The efficiency and acceptance factors are summarized in Table 1, and are calculated from simulations. The largest contribution to the acceptance can be attributed to the extrapolation from the dileptonic to the fully inclusive decays of the tt pairs. The remaining contributions are due to changes in kinematic requirements on the generatorlevel objects. As discussed in Section 4, the ttcL and ttbL categories result from ttcc and ttbb events, respectively, where one of the additional HF jets is outside the acceptance or both are merged into one jet. Therefore, these components are scaled with the same factors as the ttcc and ttbb templates, respectively, such that the relative yield of ttcL (ttbL) with respect to ttcc (ttbb) events is fixed to that predicted in the simulations. The predicted yields from the simulations are denoted by N MC k , where k denotes the signal process. The tt+other component is scaled with the same factor as the ttLL component, motivated by their similar LF content. Uncertainties in the ratios of simulated yields are taken into account in the fit. The background processes are summed together into one template and their yield is fixed to the predictions from simulations, with uncertainties taken into account as discussed in Section 10. The sum of the cross sections for the background processes (obtained from simulations) is denoted as σ bkg . A second binned maximum likelihood fit is performed assuming Poisson statistics to extract the ratios of the ttcc and ttbb cross sections to the overall inclusive ttjj cross section (denoted R c and R b , respectively) in the fiducial phase space, with the expected yield in bin i calculated using the function: with The parameters R cL and R bL are used to denote, respectively, the ratios of the ttcL and ttbL cross sections to the inclusive ttjj cross section, and are defined as a function of R c and R b in Eq. (6).

Systematic uncertainties
This section summarizes the systematic uncertainties related to the extraction of the ttcc, ttbb, and ttLL cross sections (and the ratios R c and R b ), as well as corrections applied to the simulated events to account for differences with respect to data. Systematic uncertainties are treated as nuisance parameters in the template fit to data that can affect both the shapes of the templates and the yields of the signal and background processes. A smoothing procedure [50] is applied to the templates that describe the uncertainty variations affecting the template shape. The sources of systematic uncertainties are subdivided into experimental and theoretical components and are discussed below.
Experimental uncertainties: These uncertainties affect both the shape and normalization of the templates. The jet energy resolution is known to be worse in data than in simulations, and a corresponding additional smearing is applied to the simulated jet energies [42]. Systematic uncertainties are estimated by varying the smearing of the jet energy within its uncertainties in the calculation of the cross sections. Similarly, we take into account corrections and uncertainties from observed differences in the jet energy scale. These corrections are evaluated and applied in different regions of jet p T and |η|. Observed differences in electron and muon identification, isolation, reconstruction, and trigger efficiencies between data and simulations are taken into account through p T -and η-dependent scale factors, with the corresponding uncertainties accounted for. The distribution of the number of pileup collisions in simulated events is reweighted to match the distribution observed in data, using an inelastic pp cross section of 69.2 mb [51]. An uncertainty related to this correction is applied by varying this inelastic cross section by ± 4.6%. An uncertainty of 2.3% [15] in the total integrated luminosity is also taken into account. The scale factors extracted from the c tagging calibration are applied to the simulated events, and corresponding uncertainties are considered. Uncertainties related to this calibration are found to be dominated by the choice of µ R and µ F in the W+ c and DY + jets control regions, affecting scale factors for c and LF jets, respectively, and by PS uncertainties in semileptonic tt events that affect scale factors for b jets, together with a significant contribution from statistical uncertainties for all jet flavours. Most of the theoretical and experimental sources of uncertainty are in common between the control regions in which the c tagging calibration is derived and the tt dileptonic signal region considered in this analysis. In such cases, the common uncertainties are considered fully correlated and evaluated simultaneously.
Theoretical uncertainties: In the ME calculation, the choice of µ R and µ F can have an impact on the kinematic distributions of the final-state objects. Uncertainties in these scales are taken into account by rescaling µ F and µ R up or down by a factor of two at the ME level [52,53]. The choice was made not to include the difference between PYTHIA and alternative PS simulations as a systematic uncertainty in this analysis. Instead a variety of parameters in PYTHIA, sensitive to the PS and hadronization, are consistently varied to assess the uncertainty in an unambiguous way. In the PS, the uncertainty in the value of the strong coupling constant (α S ) evaluated at m Z is taken into account by varying the renormalization scale of QCD emissions in the initial-and final-state radiation up and down by a factor of two. Similarly, uncertainties in the momentum transfer from b quarks to b hadrons (b fragmentation) in the PS have been included. The b fragmentation in PYTHIA is parametrized using a Bowler-Lund model [54][55][56], with uncertainties calculated by tuning the internal parameters of this model to measurements from the ALEPH [57], DELPHI [58], OPAL [59] and SLD [60] experiments. In practice, it was found that this parametrization can be varied within its uncertainties by a reweighting at the generator level of the so-called "transfer function", x b = p T (b hadron)/p T (b jet). No such detailed assessment of the c fragmentation uncertainties has yet been performed, but a similar level of variation is observed when comparing the available experimental measurements [61] with the default PYTHIA tune used in this analysis [62]. We verified that variations of the c jet transfer function (x c = p T (c hadron)/p T (c jet)) induce similar variations of the c jet p T and c tagging discriminator distributions as the analogous variations in b fragmentation for b jets, and found that these effects are comfortably covered by an uncertainty a factor of two larger than that for b fragmentation. This uncertainty is modelled as an uncertainty in the ttcc and ttcL yields in the fit in the fiducial phase space, and an additional uncertainty in the acceptance correction from the full to the fiducial phase space. Uncertainties associated with the PDF, as well as with the value of α S in the PDF of the proton are considered, following the PDF4LHC prescription [63]. For all of the aforementioned theoretical sources of uncertainty (except for the c fragmentation uncertainty), the effects of these variations on the shape of the fitted templates are taken into account in the fit, whereas their impact on the signal yields is considered as an uncertainty in the theoretical prediction to which the measurement is compared in Section 11. The residual theoretical uncertainty that enters the measured cross sections through the i terms in Eqs. (4) and (5) is accounted for by a separate theoretical uncertainty in the effi-ciency. The experimental uncertainties in the efficiency are already taken into account through the uncertainties in the normalization of the templates. The matching between the ME and PS is governed by a parameter called h damp . The value of this parameter is varied according to h damp = 1.379 +0.926 −0.505 m t [31] in a separate simulation. Since the size of this simulated data set is insufficient to reliably estimate the effect on the shapes of the templates, this uncertainty is conservatively estimated through its effect on the overall yield. The remnants of the pp collisions that do not take part in the hard scattering are referred to as the underlying event. Their kinematic distributions are tuned in the generators to match those observed in the data [31]. The resulting parametrization of the CP5 tune, used in this analysis, is varied within its uncertainties in separate simulations. Here too the effect of the uncertainty in the overall yield, rather than in the template shapes, is propagated to the measured cross sections. The effects of the theoretical uncertainties listed above on the fixed ratios of ttbL to ttbb, ttcL to ttcc, and tt+other to ttLL yields, taken from simulations in Eqs. (4) and (5), are also included in the fit. An uncertainty of 25% is assigned to the total cross section of all background processes, based on the precision of recent measurements of the dominant background processes [64,65]. The statistical uncertainty due to the finite size of the simulated samples is taken into account in the fit. In extrapolating the results from the fiducial to the full phase space, some of these theoretical uncertainties also affect the acceptance for the different signal categories. The individual and combined impacts of different sources of uncertainty in the acceptance are summarized in Table 2. The individual impacts from each source of uncertainty in the cross sections (and ratios) for the fiducial phase space after the fit are summarized in Table 3. The fitted nuisance parameters do not deviate significantly from their initial values and are not significantly constrained.
No strong correlations between any of the nuisance parameters and the fitted cross sections or ratios are observed. The dominant systematic uncertainties are related to the c tagging calibration, followed by jet energy scale and fragmentation uncertainties, as well as uncertainties related to the matching between ME and PS, and the choice of µ R and µ F scales in the ME calculation.

Results
The binning of the two-dimensional ∆ , and then analogously for the remainder of the bins. This histogram is shown in Fig. 5 after normalizing the simulated templates according to the fitted cross sections. The factors by which the templates of the ttcc, ttbb, and ttLL processes (using the POWHEG ME generator) are scaled to match the data are denoted µ tt cc , µ tt bb , and µ tt LL , respectively. Their measured values from the fit are displayed in the top panel of Fig. 5, together with their combined statistical and systematic uncertainties. An improved agreement between data and simulated predictions is also observed for the c tagging discriminators of the first and second additional jets after normalizing the simulated templates according to the fitted cross sections. This is demonstrated in Appendix B and can be directly compared to the agreement before the fit in Fig. 3. The measured cross sections in the fiducial and full phase spaces, together with their statistical  in the simulations (histograms) and in data (points), after normalizing the simulated templates according to the fitted cross sections. The lower panel shows the ratio of the yields in data to those predicted in the simulations. The brown and grey uncertainty bands denote, respectively, the statistical and total uncertainties from the fit. The factors (µ) by which the templates of the different processes (using the POWHEG ME generator) are scaled, are also displayed, together with their combined statistical and systematic uncertainties. and systematic uncertainties, are summarized in Table 4 and compared to predictions using the POWHEG and MADGRAPH5 aMC@NLO ME generators. Uncertainties in the measured values of the cross sections and ratios are determined from the points where the decrease in the logarithm of the profiled likelihood from the best-fit value intersects with 0.5. The inclusive ttcc cross section and the ratio R c are measured here for the first time. They are in agreement with the predictions from both ME generators, within the uncertainties.
In addition, two-dimensional likelihood scans are performed over different combinations of cross sections or ratios. These are shown in Fig. 6 for the measurements in the fiducial phase space. The 68 and 95% confidence level contours are shown, along with the predictions from the POWHEG and MADGRAPH5 aMC@NLO ME generators. Agreement is observed at the level of one to two standard deviations between the measured values and simulated predictions for the ttcc, ttbb, and ttLL processes. The most significant tension is observed in the ratio R b , at the level of 2.5 standard deviations. Results for σ tt bb and R b are consistent with previous measurements targeting specifically this signature [6-13] and show the same tendency to be slightly above the predictions from simulations.

Summary
The production of a top quark pair (tt) in association with additional bottom or charm jets at the LHC provides challenges both in the theoretical modelling and experimental measurement of this process. Whereas tt production with two additional bottom jets (ttbb) has been measured by the ATLAS and CMS Collaborations at different centre-of-mass energies [6-13], this analysis presents the first measurement of the cross section for tt production with two additional charm jets (ttcc). The analysis is conducted using data from proton-proton collisions recorded by the CMS experiment at a centre-of-mass energy of 13 TeV, corresponding to an integrated luminosity of 41.5 fb −1 . The measurement is performed in the dileptonic channel of Table 4: Measured parameter values in the fiducial (upper rows) and full (lower rows) phase spaces with their statistical and systematic uncertainties listed in that order. The last two columns display the expectations from the simulated tt samples using the POWHEG or MAD-GRAPH5 aMC@NLO ME generators. The uncertainties quoted for these predictions include the contributions from the theoretical uncertainties listed in the lower rows of Table 3, as well as the uncertainty in the tt cross section. 1.51 ± 0.11 ± 0.16 1.03 ± 0.08 1.03 ± 0.09 the tt decays and relies on the use of recently developed charm jet identification algorithms (c tagging). A template fitting method is used, based on the outputs of a neural network classifier trained to identify the signal categories defined by the flavour of the additional jets. This allows the simultaneous extraction of the cross section for the ttcc, ttbb, and tt with two additional light-flavour or gluon jets (ttLL) processes. A novel multidimensional calibration of the shape of the c tagging discriminator distributions is employed, such that this information can be reliably used in the neural network classifier.

Result
The ttcc cross section is measured for the first time to be 0.207 ± 0.025 (stat) ± 0.027 (syst) pb in the fiducial phase space (matching closely the sensitive region of the detector) and 10.1 ± 1.2 (stat) ± 1.4 (syst) pb in the full phase space. The ratio of the ttcc to the inclusive tt + two jets cross sections is found to be (3.01 ± 0.34 (stat) ± 0.31 (syst))% in the fiducial phase space and (3.36 ± 0.38 (stat) ± 0.34 (syst))% in the full phase space. These results are compared with predictions from two different matrix element generators with next-to-leading order accuracy in quantum chromodynamics, in which an inclusive description of the tt process, with up to two additional radiated hard gluons at the ME level, is interfaced with a parton shower simulation to generate the additional radiation. Agreement is observed at the level of one to two standard deviations between the measured values and simulated predictions for the ttcc, ttbb, and ttLL processes. The observed ratio of the ttbb to the inclusive tt + two jets cross sections exceeds the predictions by about 2.5 standard deviations, consistent with the tendency seen in previous measurements.  [8] ATLAS Collaboration, "Measurements of fiducial cross sections for tt production with one or two additional b jets in pp collisions at √ s = 8 TeV using the ATLAS detector", Eur. Phys. J. C 76 (2016) 11, doi:10.1140/epjc/s10052-015-3852-4, arXiv:1508.06868.

A Charm tagging discriminators before calibration
The c tagging algorithms are trained on simulated events and are therefore prone to mismodelling effects in the input variables. The initial CvsL and CvsB c tagging discriminator distributions of the first and second additional jet are shown in Fig. 7. Discrepancies between the data and the simulated predictions of up to 50% are observed, demonstrating the need for the c tagging calibration described in Section 8.

B Charm tagging discriminators after scaling to the fitted cross sections.
After scaling the simulated templates of the c tagging discriminators to their fitted yields, the agreement between data and simulations is shown in Fig. 8. This can be compared to the agreement before the fit was performed in Fig. 3, and indeed shows an improved agreement especially in those bins which have a relatively large contribution from ttbb and ttcc events.   : Comparison between data (points) and simulated predictions (histograms) for the CvsL (left column) and CvsB (right column) c tagging discriminator distributions of the first (upper row) and second (lower row) additional jet, after normalizing the simulated templates according to the fitted cross sections. The lower panels show the ratio of the yields in data to those predicted in the simulations. The brown and grey uncertainty bands denote, respectively, the statistical and total uncertainties from the fit. The factors (µ) by which the templates of the different processes (using the POWHEG ME generator) are scaled, are also displayed, together with their combined statistical and systematic uncertainties.