Measurement of the CKM angle γ using B± → DK± with D → KS0π+π−, KS0K+K− decays

A binned Dalitz plot analysis of B± → DK± decays, with D → KS0π+π− and D → KS0K+K−, is used to perform a measurement of the CP-violating observables x± and y±, which are sensitive to the Cabibbo-Kobayashi-Maskawa angle γ. The analysis is performed without assuming any D decay model, through the use of information on the strong-phase variation over the Dalitz plot from the CLEO collaboration. Using a sample of proton-proton collision data collected with the LHCb experiment in 2015 and 2016, and corresponding to an integrated luminosity of 2.0 fb−1, the values of the CP violation parameters are found to be x− = (9.0 ± 1.7 ± 0.7 ± 0.4) × 10−2, y− = (2.1 ± 2.2 ± 0.5 ± 1.1) × 10−2, x+ = (−7.7 ± 1.9 ± 0.7 ± 0.4) × 10−2, and y+ = (−1.0 ± 1.9 ± 0.4 ± 0.9) × 10−2. The first uncertainty is statistical, the second is systematic, and the third is due to the uncertainty (on the strong-phase measurements. These values are used to obtain γ = (87− 12+ 11)∘, rB = 0.086− 0.014+ 0.013, and δB = (101±11)°, where rB is the ratio between the suppressed and favoured B-decay amplitudes and δB is the corresponding strong-interaction phase difference. This measurement is combined with the result obtained using 2011 and 2012 data collected with the LHCb experiment, to give γ = (80− 9+ 10)∘, rB = 0.080 ± 0.011, and δB = (110 ± 10)°.


Introduction
The Standard Model (SM) description of CP violation [1,2] can be tested by overconstraining the angles of the Unitarity Triangle. The Cabibbo-Kobayashi-Maskawa (CKM) angle γ ≡ arg(−V ud V * ub /V cd V * cb ) is experimentally accessible through the interference between b →cus andb →ūcs transitions. It is the only CKM angle easily accessible in tree-level processes and it can be measured with negligible uncertainty from theory [3]. Hence, in the absence of new physics effects at tree level, a precision measurement of γ provides a SM benchmark that can be compared with other CKM-matrix observables more likely to be affected by physics beyond the SM. Such comparisons are currently limited by the uncertainty on direct measurements of γ, which is about 5 • [4] and is driven by the LHCb average.
The effects of interference betweenb →cus andb →ūcs transitions can be probed by studying CP -violating observables in B ± → DK ± decays, where D represents a D 0 or a D 0 meson reconstructed in a final state that is common to both [5][6][7]. This decay mode has been studied at LHCb with a wide range of D-meson final states to measure observables with sensitivity to γ [8 -11]. In addition to these studies, other B decays have also been used with a variety of techniques to determine γ [12][13][14][15].
This paper presents a model-independent study of the decay mode B ± → DK ± , using D → K 0 S π + π − and D → K 0 S K + K − decays (denoted D → K 0 S h + h − decays). The analysis utilises pp collision data accumulated with LHCb in 2015 and 2016 at a centre-of-mass energy of √ s = 13 TeV and corresponding to a total integrated luminosity of 2.0 fb −1 . The result is combined with the result obtained by LHCb with the same analysis technique, using data collected in 2011 and 2012 (Run 1) at centre-of-mass energies of √ s = 7 TeV and √ s = 8 TeV [9]. The sensitivity to γ is obtained by comparing the distributions in the Dalitz plots of D → K 0 S h + h − decays from reconstructed B + and B − mesons [6,7]. For this comparison, the variation of the strong-phase difference between D 0 and D 0 decay amplitudes within the Dalitz plot needs to be known. An attractive, model-independent, approach makes use of direct measurements of the strong-phase variation over bins of the Dalitz plot [6,16,17]. The strong phase can be directly accessed by exploiting the quantum correlation of D 0 D 0 pairs from ψ(3770) decays. Such measurements have been performed by the CLEO collaboration [18] and have been used by the LHCb [9] and Belle [19] collaborations to measure γ in B ± → DK ± decays, and have also been used to study B 0 → DK * 0 decays [20,21]. An alternative method relies on amplitude models of D → K 0 S h + h − decays, determined from flavour-tagged D → K 0 S h + h − decays, to predict the strong-phase variation over the Dalitz plot. This method has been used for a variety of B decays [22][23][24][25][26][27][28].
The separation of data into binned regions of the Dalitz plot leads to a loss of statistical sensitivity in comparison to using an amplitude model [16,17]. However, the advantage of using the direct strong-phase measurements resides in the model-independent nature of the systematic uncertainties. Where the direct strong-phase measurements are used, there is only a systematic uncertainty associated with the finite precision of such measurements. Conversely, systematic uncertainties associated with determining a phase from an amplitude model are difficult to evaluate, as common approaches to amplitude-model building break the optical theorem [29]. Therefore, the loss in statistical precision is compensated by reliability in the evaluation of the systematic uncertainty, which is increasingly important as the overall precision on the CKM angle γ improves.

Overview of the analysis
The amplitude of the decay B − → DK − , D → K 0 S h + h − can be written as a sum of the favoured B − → D 0 K − and suppressed B − → D 0 K − contributions as where m 2 − and m 2 + are the squared invariant masses of the K 0 S h − and K 0 S h + particle combinations, respectively, that define the position of the decay in the Dalitz plot, . The cosine of δ D (m 2 − , m 2 + ) weighted by the D-decay amplitude and averaged over bin i is written as c i [6], and is given by where the integrals are evaluated over the phase space of bin i. An analogous expression can be written for s i , which is the sine of the strong-phase difference weighted by the decay amplitude and averaged over the bin phase space. The values of c i and s i have been directly measured by the CLEO collaboration, exploiting quantum-correlated D 0 D 0 pairs produced at the ψ(3770) resonance [18]. The measurements of c i and s i are available in four different 2 × 8 binning schemes for the D → K 0 S π + π − decay. This analysis uses the 'optimal binning' scheme where the bins have been chosen to optimise the statistical sensitivity to γ, as described in Ref. [18]. The optimisation was performed assuming a strong-phase difference distribution as predicted by the BaBar model presented in Ref. [23]. For the K 0 S K + K − final state, three choices of binning schemes are available, containing 2 × 2, 2 × 3, and 2 × 4 bins. The guiding model used to determine the bin boundaries is taken from the BaBar study described in Ref. [24]. The 2 × 2 binning scheme is chosen, due to the low signal yields in the D → K 0 S K + K − mode. The same choice of bins was used in the LHCb Run 1 analysis [9]. The measurements of c i and s i are not biased by the use of a specific amplitude model in defining the bin boundaries. The choice of the model only affects this analysis to the extent that a poor model description of the underlying decay would result in a reduced statistical sensitivity of the γ measurement. The binning choices for the two decay modes are shown in Fig. 1.
The physics parameters of interest, r B , δ B , and γ, are translated into four CP observables [22] that are measured in this analysis. These observables are defined as x ± ≡ r B cos(δ B ± γ) and y ± ≡ r B sin(δ B ± γ). It follows from Eq. (1) that the expected numbers of B + and B − decays in bin i, N + i and N − i , are given by where F i are the fractions of decays in bin i of the D 0 → K 0 S h + h − Dalitz plot, and h B ± are normalisation factors, which can be different for B + and B − due to production, detection, and CP asymmetries. In this measurement, the integrated yields are not used to provide information on x ± and y ± , and so the analysis is insensitive to such effects. From Eq. (4) it is seen that studying the distribution of candidates over the D → K 0 S h + h − Dalitz plot gives access to the x ± and y ± observables. The detector and selection requirements placed on the data lead to a non-uniform efficiency over the Dalitz plot, which affects the F i parameters. The efficiency profile for the signal candidates is denoted as η(m 2 − , m 2 + ). The parameters F i can then be expressed as The values of F i are determined from the control decay mode where the D * − meson decays to D 0 π − and the D 0 meson decays to either the K 0 S π + π − or K 0 S K + K − final state. The symbol X indicates other particles which may be produced in the decay but are not reconstructed. Samples of simulated events are used to correct for the small differences in efficiency arising through unavoidable differences in selecting ( ) B → D * ± µ ∓ ( ) ν µ X and B ± → DK ± decays, as discussed further in Sect. 5. In addition to B ± → DK ± and ( ) B → D * ± µ ∓ ( ) ν µ X candidates, B ± → Dπ ± decays are selected. These provide an important control sample that is used to constrain the invariant-mass shape of the B ± → DK ± signal, as well as to determine the yield of B ± → Dπ ± decays misidentified as B ± → DK ± candidates. Note that this channel is not optimal for determining the values of F i as the small level of CP violation in the decay leads to a significant systematic uncertainty, as was reported in Ref. [30]. This uncertainty is eliminated when using the flavour-specific semileptonic decay, in favour of a smaller systematic uncertainty associated with efficiency differences.
The effect of D 0 -D 0 mixing was ignored in the above discussion. If the parameters F i are obtained from where the D * − decays to D 0 π − , D 0 -D 0 mixing has been shown to lead to a bias of approximately 0.2 • in the γ determination [31], which is negligible for the current analysis. The effects of CP violation in the neutral kaon system and of the different nuclear interaction cross-sections for K 0 and K 0 mesons are discussed in Sect. 7, where a systematic uncertainty is assigned.

Detector and simulation
The LHCb detector [32,33] is a single-arm forward spectrometer covering the pseudorapidity range 2 < η < 5, designed for the study of particles containing b or c quarks. The detector includes a high-precision tracking system consisting of a siliconstrip vertex detector surrounding the pp interaction region, a large-area silicon-strip detector located upstream of a dipole magnet with a bending power of about 4 Tm, and three stations of silicon-strip detectors and straw drift tubes placed downstream of the magnet. The polarity of the dipole magnet is reversed periodically throughout data-taking.
The tracking system provides a measurement of momentum, p, of charged particles with relative uncertainty that varies from 0.5% at low momentum to 1.0% at 200 GeV/c. The minimum distance of a track to a primary vertex (PV), the impact parameter (IP), is measured with a resolution of (15 + 29/p T ) µm, where p T is the component of the momentum transverse to the beam, in GeV/c. Different types of charged hadrons are distinguished using information from two ring-imaging Cherenkov detectors. Photons, electrons, and hadrons are identified by a calorimeter system consisting of scintillating-pad and preshower detectors, an electromagnetic calorimeter and a hadronic calorimeter. Muons are identified by a system composed of alternating layers of iron and multiwire proportional chambers.
The online event selection is performed by a trigger, which consists of a hardware stage based on information from the calorimeter and muon systems, followed by a software stage, which applies a full event reconstruction. At the hardware trigger stage, events are required to have a muon with high p T or a hadron, photon or electron with high transverse energy in the calorimeters. For hadrons, the transverse energy threshold is 3.5 GeV. The software trigger requires a two-, three-or four-track secondary vertex with a significant displacement from any primary pp interaction vertex. At least one charged particle must have transverse momentum p T > 1.6 GeV/c and be inconsistent with originating from a PV. A multivariate algorithm [34] is used for the identification of secondary vertices consistent with the decay of a b hadron. Small changes in the trigger thresholds were made throughout both years of data taking.
In the simulation, pp collisions are generated using Pythia 8 [35] with a specific LHCb configuration [36]. Decays of hadronic particles are described by EvtGen [37], in which final-state radiation is generated using Photos [38]. The interaction of the generated particles with the detector, and its response, are implemented using the Geant4 toolkit [39] as described in Ref. [40].
4 Event selection and fit to the invariant-mass spectrum for B ± → DK ± and B ± → Dπ ± decays Decays of K 0 S mesons to the π + π − final state are reconstructed in two categories, the first containing K 0 S mesons that decay early enough for the pions to be reconstructed in the vertex detector and the second containing K 0 S mesons that decay later such that track segments of the pions cannot be formed in the vertex detector. These categories are referred to as long and downstream, respectively. The candidates in the long category have better mass, momentum and vertex resolution than those in the downstream category.
Hereinafter, B ± candidates are denoted long or downstream depending on which category of K 0 S candidate they contain. For many of the quantities used in the selection and analysis of the data, a kinematic fit [41] is imposed on the full B ± decay chain. Depending on the quantity being calculated, the D and K 0 S candidates may be constrained to have their known masses [42], as described below. The fit also constrains the B ± candidate momentum vector to point towards the associated PV, defined as the PV for which the candidate has the smallest IP significance. These constraints improve the resolution of the calculated quantities, and thus help improve separation between signal and background decays. Furthermore, it improves the resolution on the Dalitz plot coordinates and ensures that all candidates lie within the kinematically allowed D → K 0 S h + h − phase space. The D (K 0 S ) candidates are required to be within 25 MeV/c 2 (15 MeV/c 2 ) of their known mass [42]. These requirements are placed on masses obtained using kinematic fits in which all constraints are applied except for that on the mass under consideration. Combinatorial background is primarily suppressed through the use of a boosted decision tree (BDT) multivariate classifier [43,44]. The BDT is trained on simulated signal events and background taken from the high B ± mass sideband (5800-7000 MeV/c 2 ). Separate BDTs are trained for the long and downstream categories.
Each BDT uses the same set of variables: the χ 2 of the kinematic fit of the whole decay chain; p and p T of the companion, D, and B ± after the kinematic refit (here and in the following, companion refers to the final state π ± or K ± meson produced in the B ± → Dh ± decay); the vertex quality of the K 0 S , D, and B ± candidates; the distance of closest approach between tracks forming the D and B ± vertices; the cosine of the angle between the momentum vector and the vector between the production and decay vertices of a given particle, for each of the K 0 S , D, and B ± candidates; the minimum and maximum values of the χ 2 IP of the pions from both the D and K 0 S decays, where χ 2 IP is defined as the difference in χ 2 of the PV fit with and without the considered particle; the χ 2 IP for the companion, K 0 S , D, and B ± candidates; the B ± flight-distance significance; the radial distance from the beamline to the D and B ± -candidate vertices; and a B ± isolation variable, which is designed to ensure the B ± candidate is well isolated from other tracks in the event. The B ± isolation variable is the asymmetry between the p T of the signal candidate and the sum of the p T of other tracks in the event that lie within a distance of 1.5 in η-φ space, where φ is the azimuthal angle measured in radians. Candidates in the data samples that have a BDT output value below a threshold are rejected. An optimal threshold value is determined for each of the two BDTs, using a series of pseudoexperiments to obtain the values that provide the best sensitivity to x ± and y ± . Across all B ± → DK ± channels this requirement is found to reject 99.1 % of the combinatorial background in the high B mass sideband that survives all other requirements, while having an efficiency of 92.4 % on simulated B ± → DK ± signal samples.
Particle identification (PID) requirements are placed on the companion to separate B ± → DK ± and B ± → Dπ ± candidates, and on the charged decay products of the D meson to remove cross-feed between different D → K 0 S h + h − decays. To ensure good control of the PID performance it is required that information from the RICH detectors is present. To remove background from D → π + π − π + π − or D → π + π − K + K − decays, long K 0 S candidates are required to have travelled a significant distance from the D vertex. This requirement is not necessary for downstream candidates. Similarly, the D decay vertex is required be significantly displaced from the B ± decay vertex in order to remove charmless B ± decays.
The Dalitz plots for B ± → DK ± candidates in a narrow region of ±25 MeV/c 2 around the B ± mass are shown in Fig. 2, for both D → K 0 S h + h − final states samples. Separate plots are shown for B + and B − decays. The Dalitz coordinates are calculated from the kinematic fit with all mass constraints applied.
In order to determine the parameterisation of the signal and background components that are used in the fit of partitioned regions of the Dalitz plot described in Sect. 6, an extended maximum likelihood fit to the invariant-mass distributions of the B ± candidates is performed, in which the B + and B − candidates in all of the Dalitz bins are combined. The invariant mass of each B ± candidate is calculated using the results of a kinematic fit in which the D and K 0 S masses are constrained to their known values. The sample is split into B ± → DK ± and B ± → Dπ ± candidates, by D decay mode and by K 0 S category. In order to allow sharing of some parameters, the fit is performed simultaneously for all of the above categories. The projections of the fit and the invariant-mass distributions of the selected B ± candidates are shown in Figs. 3 and 4 for D → K 0 S π + π − and D → K 0 S K + K − candidates, respectively. The fit range is between 5080 MeV/c 2 and 5800 MeV/c 2 in the B ± candidate invariant mass.
The peaks corresponding to actual B ± → DK ± and B ± → Dπ ± candidates are fitted with a sum of two Crystal Ball [45] functions, which are parameterised as where α > 0, and The sum is implemented such that the Crystal Ball functions have tails pointing in either direction. They share a common width, σ, and mean, µ. In practice, the signal probability density function (PDF) is defined as The tail parameters, n L,R and α L,R , are fixed from simulation, while the other parameters are left as free parameters in the fit. Separate tail parameters, f CB , and σ are used for long and downstream candidates. Different widths are used for the B ± → DK ± and B ± → Dπ ± channels, with their ratio r σ = σ DK /σ Dπ shared between all categories. The mean is shared among all categories. The yield of B ± → Dπ ± decays in each K 0 S and D-meson decay category, N cat (B ± → Dπ ± ), is determined in the fit. Instead of fitting the yield of B ± → DK ± decays directly in each category, it is determined from the B ± → Dπ ± yield in the corresponding category and the ratio which is a free parameter in the fit. The category-dependent PID efficiencies, ε cat PID (B ± → Dh ± ), are taken into account, so that a single R parameter can be shared between all categories in the fit. How these efficiencies are obtained is described below.
As the parameter R is efficiency corrected, it is equal to the ratio of branching fractions between the B ± → DK ± and B ± → Dπ ± decay modes. The measured ratio is found to be R = (7.66 ± 0.14) %, where the uncertainty is statistical only, and this is consistent with the expected value of (7.8 ± 0.4) % [42]. The background consists of random track combinations, partially reconstructed B decays, and B ± → Dh ± decays in which the companion has been misidentified. The random track combinations are modelled by an exponential PDF. The slopes of the exponentials are free parameters in the fit to the data. These slopes are independent for each of the B ± → Dπ ± categories, while they are shared for the B ± → DK ± categories to improve the stability of the fit. When these slopes are allowed to be independent, the fit returns results that are statistically compatible.
In the B ± → DK ± sample there is a clear contribution from B ± → Dπ ± decays in which the companion particle is misidentified as a kaon by the RICH system. The rate for B ± → DK ± decays to be misidentified and placed in the B ± → Dπ ± sample is much lower due to the reduced branching fraction. Nevertheless, this contribution is still accounted for in the fit. The yields of these backgrounds are fixed in the fit, using knowledge of misidentification efficiencies and the fitted yields of reconstructed decays with the correct particle hypothesis. The misidentification efficiencies are obtained from large samples of These decays are selected using only kinematic variables in order to provide pure samples of K ∓ and π ± that are unbiased in the PID variables. The PID efficiency is parameterised as a function of the companion momentum and pseudorapidity, and the charged-particle multiplicity in the event. The calibration sample is weighted so that the distribution of these variables matches that of the candidates in the signal region of the B ± sample, thereby ensuring that the measured PID performance is representative for the samples used in this measurement. The efficiency to identify a kaon correctly is found to be approximately 86 %, while that for a pion is approximately 97 %. The PDFs of the backgrounds due to misidentified companion particles are determined using data. As an example, consider the case of true B ± → Dπ ± decays misidentified as B ± → DK ± candidates. The sPlot method [46] is used on the B ± → Dπ ± sample in order to isolate the contribution from the signal decays. The B ± invariant mass is then calculated using the kaon mass hypothesis for the companion pion, and weighting by PID efficiencies in order to properly reproduce the kinematic properties of pions misidentified as kaons in the signal B ± → DK ± sample. The weighted distribution is fitted with a sum of two Crystal Ball shapes. The fitted parameters are subsequently fixed in the fit to the B ± invariant-mass spectrum, with the procedure applied separately for long and downstream candidates. An analogous approach is used to determine the shape of the misidentified B ± → DK ± contribution in the B ± → Dπ ± sample.
Partially reconstructed b-hadron decays contaminate the sample predominantly at invariant masses smaller than that of the signal peak. These decays contain an unreconstructed pion or photon, which predominantly comes from an intermediate resonance.
There are contributions from B 0 → D * ± h ∓ and B ± → D * 0 h ± decays in all channels (denoted as B → D * h ± decays), while B ± → Dρ ± and B ± → DK * ± decays contribute to the B ± → Dπ ± and B ± → DK ± channels, respectively. In the B ± → DK ± channels there is also a contribution from B 0 where the charged pion is not reconstructed. The invariant-mass distributions of these backgrounds depend on the spin and mass of the missing particle, as described in Ref. [47]. The shape of the background from B 0 s decays is based on the results of Ref. [48]. Additionally, each of the above backgrounds of B ± → Dπ ± decays can contribute in the B ± → DK ± channels if the pion is misidentified. The inverse contribution is negligible and is neglected. The shapes for the decays in which a pion is misidentified as a kaon are parameterised with semi-empirical PDFs formed from sums of Gaussian and error functions. The parameters of these backgrounds are fixed to the results of fits to data from two-body D decays [47], where they were obtained with a much larger data sample. However, the width of the resolution function and a shift along the B ± mass are allowed to differ in order to accommodate small differences between the D decay modes.
In each of the B ± → Dπ ± channels, the total yield of the partially reconstructed background is fitted independently. The relative amount of each B → D * π ∓ mode is fixed from efficiencies obtained from simulation and known branching fractions, while the fraction of B ± → Dρ ± decays is left free. In the B ± → DK ± channels, the yield of the B 0 s → D 0 π + K − background is fixed relative to the corresponding B ± → Dπ ± yield, using efficiencies from simulation and the known branching fraction. The total yield of the remaining partially reconstructed backgrounds is expressed via a single fraction, R low DK/Dπ , relative the B ± → Dπ ± yields. It is free in the fit, and common to all channels after taking into account the different particle-identification efficiencies. The relative amount of each B → D * K ∓ mode is fixed using efficiencies from simulation and known branching fractions, while the fraction of B ± → DK * ± decays is fixed using the results of Ref. [47]. The yields of the partially reconstructed modes with a companion pion misidentified as a kaon are fixed via the known PID efficiencies, based on the fitted yield of the partially reconstructed backgrounds in the corresponding B ± → Dπ ± channel.
In the B ± → Dπ ± channels, a total signal yield of approximately 56 100 (7750) is found in the signal region 5249-5319 MeV/c 2 of the D → K 0 S π + π − (D → K 0 S K + K − ) channel, 31 % (32 %) of which are in the long K 0 S category. The purity in the signal region is found to be 98.4 % (97.7 %), with the dominant background being combinatorial. In the B ± → DK ± channels, a total signal yield of approximately 3900 (530) is found in the signal region of the D → K 0 S π + π − (D → K 0 S K + K − ) channel, again finding 31 % (32 %) of the candidates in the long K 0 S category. The purity in the signal region is found to be 81 %. The dominant background is from misidentified B ± → Dπ ± decays, which accounts for 66 % of the background in the signal region. Equal amounts of combinatorial background and partially reconstructed decays, predominantly including a misidentified companion pion, make up the remaining background. (5), as the expected fractions of D 0 decays falling into the ith Dalitz plot bin, taking into account the efficiency profile of the signal decay. The semileptonic decay of the B meson and the strong-interaction decay of the D * ± meson allow the flavour of the D 0 meson to be determined from the charges of the muon and the soft pion from the D * ± decay. This particular decay chain, involving a flavour-tagged D 0 decay, is chosen due to its high yield, low background level, and low mistag probability. The selection requirements are chosen to minimise changes to the efficiency profile with respect to those associated with the B ± → DK ± channels.

Event selection and yield determination for
The selection is identical to that applied in Ref. [9], except for a tighter requirement on the significance of the D 0 flight distance that helps to suppress backgrounds from charmless B decays. To improve the resolution of the distribution of candidates across the Daltiz plot, the B-decay chain is refitted [41] with the D 0 and K 0 S candidates constrained to their known masses. An additional fit, in which only the K 0 S mass is constrained, is performed to improve the D 0 and D * ± mass resolution in the invariant-mass fit used to determine signal yields.
The invariant mass of the D 0 candidate, m(K 0 S h + h − ), and the invariant-mass difference, , are fitted simultaneously to determine the signal yields. This two-dimensional parameterisation allows the yield of selected candidates to be measured in three categories: true D * ± candidates (signal), candidates containing a true D 0 meson but random soft pion (RSP) and candidates formed from random track combinations that fall within the fit range (combinatorial background). Background contributions from real D * ± decays paired with a random µ are determined to be negligible by selecting pairs of D * ± mesons and µ ± with the same charge.
An example projection of m(K 0 S π + π − ) and ∆m is shown in Fig. 5. The result of a two-dimensional extended, unbinned, maximum likelihood fit is superimposed. The fit is performed simultaneously for the two D 0 final states and the two K 0 S categories with some parameters allowed to be independent between categories. Candidates selected from data recorded in 2015 and 2016 are fitted separately, in order to accommodate different trigger threshold settings that result in slightly different Dalitz plot efficiency profiles. The fit region is defined by 1830 < m(K 0 S h + h − ) < 1910 MeV/c 2 and 139.5 < ∆m < 153.0 MeV/c 2 . Within this m(K 0 S h + h − ) range, the ∆m resolution does not vary significantly. The signal is parameterised in ∆m with a sum of two Crystal Ball functions, as for the B ± → Dh ± signal. The mean, µ, is shared between all categories, while the other parameters are different for long and downstream candidates. The tail parameters are  fixed from simulation. The combinatorial and RSP backgrounds are both parameterised with an empirical model given by for ∆m−∆m 0 > 0 and f (∆m) = 0 otherwise, where ∆m 0 , x, p 1 , and p 2 are free parameters. The parameter ∆m 0 , which describes the kinematic threshold for a D * ± → ( ) D 0 π ± decay, is shared in all data categories and for both the combinatorial and RSP shapes. The remaining parameters are determined separately for D → K 0 S π + π − and D → K 0 S K + K − candidates.
In the m(K 0 S h + h − ) fit, all of the parameters in the signal and RSP PDFs are constrained to be the same as both describe a true D 0 candidate. These are also fitted with a sum of two Crystal Ball functions, with the tail parameters fixed from simulation. The parameters are fitted separately for the D → K 0 S π + π − and D → K 0 S K + K − shapes, due to the different phase space available in the D 0 decay. The combinatorial background is parameterised by an exponential function in m(K 0 S h + h − ). A total signal yield of approximately 113 000 (15 000) D → K 0 S π + π − (D → K 0 S K + K − ) decays is obtained. This is approximately 25 times larger than the B ± → DK ± yield. In the range surrounding the signal peaks, defined as 1840-1890 (1850-1880) MeV/c 2 in m(K 0 S π + π − ) (m(K 0 S K + K − )) and 143.9-146.9 MeV/c 2 in ∆m, the background components account for 2-5 % of the yield depending on the category.
The two-dimensional fit in m(K 0 S h + h − ) and ∆m of the ( ) B → D * ± µ ∓ ( ) ν µ X decay is repeated in each Dalitz plot bin with all of the PDF parameters fixed, resulting in a raw control-mode yield, R i , for each bin i. The measured R i are not equivalent to the F i fractions required to determine the CP parameters due to unavoidable differences from selection criteria in the efficiency profiles of the signal and control modes. Examples of the efficiency profiles from simulation of the downstream candidates in 2016 data are shown in Fig. 6. For each Dalitz plot bin i a correction factor ξ i is determined to account for B → D * ± µ ∓ ( ) ν µ X decays in the simulation. The top (bottom) plots are for D → K 0 S π + π − (D → K 0 S K + K − ) decays. These plots refer to downstream K 0 S candidates under 2016 data taking conditions. The normalisation is chosen so that the average over the Dalitz plot is unity. these efficiency differences, defined as where η(m 2 − , m 2 + ) Dπ and η(m 2 − , m 2 + ) D * µ are the efficiency profiles of the B ± → Dπ ± and ( ) B → D * ± µ ∓ ( ) ν µ X decays, respectively, and are determined from simulation. The B ± → Dπ ± decay mode is used rather than B ± → DK ± as the simulation is more easily compared to the data, due to the larger decay rate and the smaller interference between B ± → D 0 π ± and B ± → D 0 π ± decays, compared to in the B ± → DK ± decay mode. It is verified using simulation that the efficiency profiles of the B ± → Dπ ± and B ± → DK ± decays are the same. The simulated events are generated with a flat distribution across the D → K 0 S h + h − phase space; hence the distribution of simulated events after triggering, reconstruction and selection is directly proportional to the efficiency profile. The amplitude models used to determine the Dalitz plot intensity for the correction factor are those from Ref.
[24] for the D → K 0 S π + π − and D → K 0 S K + K − decays, respectively. The amplitude models provide a description of the intensity distribution over the Dalitz plot and introduce no significant model dependence into the analysis. The F i values can be determined via the relation F i = h ξ i R i , where h is a normalisation factor such that the sum of all F i is unity. The F i values are determined separately for each year of data taking and K 0 S category and are then combined in the fractions observed in the B ± → Dπ ± signal region in data. This method of determining the F i parameters is preferable to using solely the amplitude models and B ± → Dπ ± simulated events, since the method is data-driven. The amplitude models and simulation data enter the correction factor as a ratio, and thus imperfections in the simulation and the model cancel at first order. The average correction factor over all bins is approximately 2 % from unity and the largest correction factor is within 7 %. Uncertainties on these correction factors are driven by the size of the simulation samples and are of a similar size as the corrections themselves.
6 Dalitz plot fit to determine the CP -violating parameters x ± and y ± The Dalitz plot fit is used to measure the CP -violating parameters x ± and y ± , as introduced in Sect. 2. Following Eq. (4), these parameters are determined from the populations of the B + and B − Dalitz plot bins, given the external information of the c i and s i parameters from CLEO-c data and the values of F i from the semileptonic control decay modes. Although the absolute numbers of B + and B − decays integrated over the D Dalitz plot have some dependence on x ± and y ± , the sensitivity gained compared to using just the relations in Eq. (4) is negligible [49] given the available sample size. Consequently, as stated previously, the integrated yields are not used to provide information on x ± and y ± and the analysis is insensitive to B meson production and detection asymmetries. A simultaneous fit is performed on the B ± → Dh ± data, split into the two B charges, the two K 0 S categories, the B ± → DK ± and B ± → Dπ ± candidates, and the two D → K 0 S h + h − final states. The invariant mass of each B ± candidate is calculated using the results of a kinematic fit in which the D and K 0 S masses are constrained to their known values. Each category is then divided into the Dalitz plot bins shown in Fig. 1, where there are 16 bins for D → K 0 S π + π − and 4 bins for D → K 0 S K + K − . The B ± → DK ± and B ± → Dπ ± samples are fitted simultaneously because the yield of B ± → Dπ ± signal in each Dalitz plot bin is used to determine the yield of misidentified candidates in the corresponding B ± → DK ± Dalitz plot bin. The PDF parameters for both the signal and background invariant-mass distributions are fixed to the values determined in the invariant-mass fit described in Sect. 4. The B ± mass range is reduced to 5150-5800 GeV/c 2 to avoid the need of a detailed description of the shape of the partially reconstructed background. The yields of signal candidates for each bin in the B ± → Dπ ± sample are free parameters. In each of the B ± → DK ± channels, the total yield integrated over the Dalitz plot is a free parameter. The fractional yields in each bin are defined using the expressions for the Dalitz plot distribution in terms of x ± , y ± , F i , c i , and s i in Eq. (4), where the x ± and y ± parameters are free and the values of F i are Gaussian-constrained within their uncertainties. The values of c i and s i are fixed to their central values, which is taken into account as a source of systematic uncertainty. The yields of the component due to B ± → Dπ ± decays, where the companion has been misidentified as a kaon, are fixed in each B ± → DK ± bin, relative to the yield in the corresponding B ± → Dπ ± bin, using the known PID efficiencies. A component for misidentified B ± → DK ± decays in the B ± → Dπ ± channels is not included, as it is found to contribute less than 0.5% of the yield in the signal region in the global fit described in Sect. 4. The total yield of the partially reconstructed B ± and B 0 backgrounds is fitted in each bin, using the same shape in all bins, with the fractions of each component taken from the global fit. The total yield of the B 0 s → D 0 π + K − (B 0 s → D 0 π − K + ) background is fixed in each channel, using the results of the global fit. The yield in each bin is then fixed from the F i parameters, using the known Dalitz distribution of D 0 (D 0 ) → K 0 S h + h − decays. The separate treatment of the partially reconstructed background from B 0 s decays is necessary due to the significantly different Dalitz distribution, arising because only a D 0 meson is produced along with a K − meson, while for the remaining modes, the D meson is either a D 0 meson or an admixture where the D 0 component is r B -suppressed. The yield of the combinatorial background in each bin is a free parameter. In bins in which an auxiliary fit determines the yield of the partially reconstructed or combinatorial background to be negligible, the corresponding yields are set to zero to facilitate the calculation of the covariance matrix [50,51].
A large ensemble of pseudoexperiments is performed to validate the fit procedure. In each pseudoexperiment the numbers and distributions of signal and background candidates are generated according to the expected distribution in data, and the full fit procedure is then executed. The input values for x ± and y ± correspond to γ = 70 • , r B = 0.1, and δ B = 130 • . The uncertainties determined by the fit to data are consistent with the size of the uncertainties determined by the pseudoexperiments. Small biases are observed in the central values and are due to the low event yields in some of the bins. These biases are observed to decrease in simulated experiments of larger size. The central values are corrected for the biases and a systematic uncertainty is assigned, as described in Sect. 7.
The total B ± → DK ± yields in the signal region, where the invariant mass of the B candidate is in the interval 5249-5319 MeV/c 2 , are shown in Table 2.
The measured values of (x ± , y ± ) from the fit to data are displayed in Fig. 7, along with their likelihood contours, corresponding to statistical uncertainties only. The systematic uncertainties are discussed in the next section. The two vectors defined by the coordinates (x − , y − ) and (x + , y + ) are not consistent with zero magnitude and they have a non-zero opening angle. Therefore the data sample exhibits the expected features of CP violation. The opening angle is equal to 2γ, as illustrated in Fig. 7.
In order to assess the goodness of fit, and to demonstrate that the equations in (x ± , y ± ) provide a good description of data, an alternative fit is performed where the B ± → DK ± yields are measured independently in each bin. In Fig. 8 (left) the obtained yields are compared with the yields predicted from the values of (x ± , y ± ) obtained in the default fit. The yields from the direct fit agree with the prediction with a p-value of 0.33. In Fig. 8 (right) the difference N i B + − N −i B − in each bin is calculated using the results of the direct fit of the B ± → DK ± yields. This distribution is compared to that predicted by the central (x ± , y ± ) values. The measured yield differences are compatible with the prediction with a p-value of 0.58. In addition, data are fitted with the assumption of no CP violation by enforcing x + = x − ≡ x 0 and y + = y − ≡ y 0 . The obtained x 0 and y 0 values are used to determine the predicted values of N i B + − N −i B − , which are also shown in Fig. 8 (right). This prediction is not zero because the B meson production and various detection effects can induce a global asymmetry in the measured yields. The comparison of the data to this hypothesis yields a p-value of 1 × 10 −6 , which strongly disfavours the CP -conserving hypothesis.

Systematic uncertainties
Systematic uncertainties on the measurements of the x ± and y ± parameters are evaluated and are presented in Table 3. The source of each systematic uncertainty is described below. The systematic uncertainties are generally determined from an ensemble of pseudoexperiments where the simulated data are generated in an alternative configuration and fitted with the default method. The mean shifts in the fitted values of x ± and y ± in comparison to their input values are taken as the systematic uncertainty.
The limited precision on (c i , s i ) coming from the CLEO measurement induces uncertainties on x ± and y ± [18]. These uncertainties are evaluated by fitting the data multiple times, each with different (c i , s i ) values sampled according to their experimental uncertainties and correlations. The resulting widths in the distributions of x ± and y ± values are assigned as the systematic uncertainties. Values of (0.4-1.1) × 10 −2 are found for the fit to the full sample. The uncertainties are similar to, but different from, those reported in Ref. [9]. This is as expected since it is found from simulation studies that the (c i , s i )-related uncertainty depends on the particular sample under study. It is found that the uncertainties do become constant when simulated samples with very high signal yields are studied. The uncertainties arising from the CLEO measurements are kept separate from the other experimental uncertainties.
A systematic uncertainty arises from imperfect modelling in the simulation used to derive the efficiency correction for the determination of the F i parameters. As the simulation enters the correction in a ratio, it is expected that imperfections cancel to first order. To determine the residual systematic uncertainty associated with this correction, an additional set of correction factors is calculated and used to evaluate an alternative set of F i parameters. To determine this additional factor, a new rectangular binning scheme is used, which is shown in Fig. 9. The bin-to-bin efficiency variation in this rectangular scheme is significantly larger than for the default partitioning and is more sensitive to imperfections in the simulated data efficiency profile. The yields of the B ± → Dπ ± and ( ) B → D * ± µ ∓ ( ) ν µ X decays in each bin of the rectangular scheme are compared to the predictions from the amplitude model and the simulated data efficiency profile. The usage of the rectangular binning also helps to dilute the small level of CP violation in B ± → Dπ ± such that differences from this comparison will come primarily from efficiency effects. The alternative correction factors ξ alt i are calculated as where the C(m 2 − , m 2 + ) terms are the ratios between the predicted and observed data yields in the rectangular bins. Many pseudoexperiments are performed, in which the data are generated according to the alternative F i parameters and then fitted with the default F i parameters. The overall shift in the fitted values of the CP parameters in comparison to their input values is taken as the systematic uncertainty, yielding 0.6 × 10 −2 for x ± and 0.1(0.2) × 10 −2 for y + (y − ).
Various effects are considered to assign an uncertainty for the imperfections in the description of the invariant-mass spectrum. For the PDF used to fit the signal, the parameters of the PDF used in the binned fit are varied according to the uncertainties obtained in the global fit. An alternative shape is also tested. The global fit is repeated with the mean and width of the shape used to describe the background due to misidentified companions allowed to vary freely. The results are used to generate data sets with an alternative PDF, and fit them using the default setup. The description of the partially reconstructed background is changed to a shape obtained from a fit of the PDF to simulated decays. The slope of the exponential used to fit the combinatorial background is also Two systematic uncertainties associated with the misidentified B ± → Dπ ± background in the B ± → DK ± sample are considered. First, the uncertainties on the particle misidentification probabilities are found to have a negligible effect on the measured values of x ± and y ± . Second, it is possible that the invariant-mass distribution of the misidentified background (the mis-ID shape) is not uniform over the Dalitz plot, as assumed in the fit. This can occur through kinematic correlations between the reconstruction efficiency across the Dalitz plot of the D decay and the momentum of the companion pion from the B ± decay. Alternative mass shapes are constructed by repeating the procedure used to obtain the default shape for each Dalitz bin individually. The alternative shapes are used when generating data sets for pseudoexperiments, and the fits then performed assuming a single shape, as in the fit to data. The resulting uncertainty is at most 0.2 × 10 −2 for all CP parameters.
In the fit to the data, the relative contributions of the partially reconstructed B ± and B 0 backgrounds are kept the same in each Dalitz bin. This is a simplification as some partially reconstructed backgrounds will be distributed as D 0 (D 0 ) → K 0 S h + h − for reconstructed B − (B + ) candidates, while partially reconstructed B ± → D ( * ) K ( * )± decays will be distributed as a D 0 − D 0 admixture depending on the relevant CP violation parameters. Pseudoexperiments are generated, where the D-decay Dalitz plot distribution for B ± → D * K + is based on the CP parameters reported in Ref. [52] and those for B ± → DK * + are taken from Ref. [53]. The generated samples are fitted with the standard method. The resulting uncertainty is at most 0.2 × 10 −2 for all CP parameters.
The total yield of the B 0 s → D 0 π + K − background in the B ± → DK ± channels is fixed relative to the corresponding B ± → Dπ ± yield. The systematic uncertainty due to the uncertainty on the relative rate is estimated via pseudoexperiments, where data sets are generated with the rate varied by ±1σ and fitted using the default value. The maximal mean bias for each parameter is taken as the uncertainty. The resulting uncertainty is 0.1 × 10 −2 for all CP parameters.
An uncertainty is assigned to each CP parameter to accompany the correction that is applied for the small bias observed in the fit procedure. These uncertainties are determined by performing sets of pseudoexperiments, each generated with different values of x ± and y ± throughout a range around the values predicted by the world averages. The spread in observed bias is combined in quadrature with the uncertainty in the precision of the pseudoexperiments. This is taken as the systematic uncertainty and is 0.1 × 10 −2 for all CP parameters.
The systematic uncertainty from the effect of candidates being assigned the wrong Dalitz bin number is considered. The resolution in m 2 + and m 2 − is approximately 0.006 GeV 2 /c 4 for candidates with long K 0 S decays and 0.007 GeV 2 /c 4 for candidates with downstream K 0 S decays. While this is small compared to the typical width of a bin, net migration can occur in regions where the presence of resonances cause the density to change rapidly. To first order this effect is accounted for by use of the control channel. However, differences in the distributions of the Dalitz plots due to efficiency differences or the nonzero value of r B in the signal decay may cause residual effects. The uncertainty from this is determined via pseudoexperiments, in which different input F i values are used to reflect the residual migration. The size of any possible bias is found to be 0.1 × 10 −2 for all CP parameters.
There is a systematic uncertainty related to CP violation in the neutral kaon system due to the fact that the K 0 S state is not an exact CP eigenstate and, separately, due to different nuclear interaction cross-sections of the K 0 and K 0 mesons. The measurement is insensitive to global asymmetries, but is affected by the different Dalitz distributions of D → K 0 S h − h + and D → K 0 L h − h + decays, as well as any correlations between Dalitz coordinates and the net material interaction. The potential bias on x ± and y ± is assessed using a series of pseudoexperiments, where data are generated taking the effects into account and fitted using the default fit. The D → K 0 L h − h + Dalitz distribution is estimated by transforming an amplitude model of D → K 0 S h − h + [22], following arguments and assumptions laid out in Ref. [18]. The effect of material interaction is treated using the formalism described in Ref. [54]. The size of the potential bias is found to be ≤ 0.2 × 10 −2 for all CP parameters, corresponding to a bias on γ of approximately 0.8 • , which is within expected limits [55].
The nonuniform efficiency profile over the Dalitz plot means that the values of (c i , s i ) appropriate for this analysis can differ from those measured by the CLEO collaboration, which correspond to the constant-efficiency case. Amplitude models are used to calculate the values of c i and s i both with and without the efficiency profiles determined from simulation. The models are taken from Ref.
[24] for D → K 0 S K + K − decays. The difference is taken as an estimate of the size of this effect. Pseudoexperiments are generated in which the values have been shifted by this difference, and then fitted with the default (c i , s i ) values. The resulting bias on x ± and y ± is found to be negligible.
The effect that a detection asymmetry between hadrons of opposite charge can have on the symmetry of the efficiency across the Dalitz plot is found to be negligible. Changes in the mass model used to describe the semileptonic control sample are also found to have a negligible effect on the F i values.
Finally, several checks are conducted to assess the stability of the results. These include repeating the fits separately for both K 0 S categories, for each data-taking year, and by splitting the candidates depending on whether the hardware trigger decision was due to particles in the signal-candidate decay chain or other particles produced in the pp collision. No anomalies are found and no additional systematic uncertainties are assigned.
In total the systematic uncertainties are less than half of the corresponding statistical uncertainties. The correlation matrix obtained for the combined effect of the sources of experimental and strong-phase related systematic uncertainties is given in Table 4.

Results and interpretation
The CP observables are measured to be x + = (−7.7 ± 1.9 ± 0.7 ± 0.4) × 10 −2 , y + = (−1.0 ± 1.9 ± 0.4 ± 0.9) × 10 −2 , where the first uncertainty is statistical, the second is the total experimental systematic uncertainty and the third is that arising from the precision of the CLEO measurements.
The signature for CP violation is that (x + , y + ) = (x − , y − ). The distance between (x + , y + ) and (x − , y − ) is calculated, taking all uncertainties and correlations into account, and found to be |(x + , y + ) − (x − , y − )| = (17.0 ± 2.7) × 10 −2 , which is different from zero by 6.4 standard deviations. This constitutes the first observation of CP violation in B ± → DK ± decays for the D → K 0 S h + h − final states. These results are compared to the expected central values of x ± and y ± that can be computed from r B , δ B , and γ as determined in the LHCb combination in Ref. [52], and the results are shown in Fig. 10 (the later LHCb combination in Ref. [56] includes the results of this measurement and is therefore unsuitable for comparison). The two sets of (x + , y + ) are in agreement within 1.6 standard deviations when the uncertainties and correlations of both the LHCb combination and this measurement are taken into account. There is a 2.7 standard deviation tension between the measured values of (x − , y − ) and the values calculated from the LHCb combination. This tension will be investigated further when this measurement and the LHCb combination are updated using data taken in 2017 and 2018.
The results for x ± and y ± are interpreted in terms of the underlying physics parameters γ, r B and δ B . The interpretation is done via a maximum likelihood fit using a frequentist treatment as described in Ref. [57]. The solution for the physics parameters has a two-fold The values for γ and r B are consistent with those presented in Ref. [52]. This is the most precise measurement of γ from a single analysis. The value of δ B shows some disagreement with Ref. [52], where the angle is determined to be 139.9 +4.8 The values of x ± , y ± measured in this analysis can be combined with those from the corresponding analysis of Run 1 data [9]. This procedure is done via a maximum likelihood fit, as implemented in the gammacombo package [57]. The previous measurements are identified by the index I, and the results within this paper are identified by the index II. When combining the two results, the fit determines the (x ± ,ŷ ± ) parameters that maximize the multivariate Gaussian likelihood function where z = (x I ± , y I ± , x II ± , y II ± ) T andẑ = (x ± ,ŷ ± ,x ± ,ŷ ± ) T are 8 × 1 vectors and Σ is the 8 × 8 covariance matrix The covariance matrix is expressed in terms of the covariance matrices obtained for the individual measurements, Σ I and Σ II , and the cross-covariance matrix Σ I-II describing  correlations between the measurements. The covariance matrix for this measurement, Σ II , is calculated using the total statistical and systematic uncertainties, and the correlation matrices in Tables 1 and 4. The covariance matrix for the Run 1 measurement, Σ I , is taken from Ref. [9], where it was calculated taking strong-phase-related correlations into account, but treating the experimental systematic uncertainties as uncorrelated. The impact of using the correlation matrix in Table 4 for these instead is found to be negligible. The dominant uncertainty in both measurements is the statistical uncertainty. As the measurements use independent data sets, the statistical uncertainties are uncorrelated. The cross-correlations of the systematic errors between measurements due to the strong phase inputs are obtained from the results of a series of fits to the two data sets in which the strong phases are varied identically. This mirrors the procedure used to evaluate the uncertainties within a single data set. The obtained cross-correlations between the fit results are given in Table 5. The elements on the diagonal do not have unit value because the obtained correlations depend on the specific data sets for the two measurements.
The combination is performed assuming full correlation between the non-strong-phase related experimental systematic uncertainties in Run 1 and this measurement. The correlation matrix for the experimental uncertainties of this analysis is used as the crossrun correlation of the experimental systematic uncertainties. The complete correlation matrix for the experimental and strong-phase-related systematic uncertainties is given in Table 6. The impact on the combination due to different assumptions on the crosscorrelations of the systematic uncertainties is found to be negligible. This is unsurprising as both measurements remain limited in precision by their statistical uncertainties.  Figure 11: Two-dimensional 68.3 % and 95.5 % confidence regions for (γ, r B , δ B ) for the x ± , y ± parameters obtained in the fit to 2015 and 2016 data, the fit to Run 1 data, and their combinations.
The interpretation in terms of the underlying physics parameters is performed on the combined values of x ± and y ± and the central values and their 68% (95%) confidence intervals are The results of the interpretation for both the combined and individual data sets are shown in Fig. 11, where the projections of the three-dimensional surfaces bounding the one and two standard deviation volumes on the (γ, r B ) and (γ, δ B ) planes are shown. The uncertainty on γ is inversely proportional to r B . Therefore the lower central value of r B in the combined results lead to a larger than naively expected uncertainty on γ when both data sets are used. The contribution of each source of uncertainty are estimated by performing the combination while taking only subsets of the uncertainties into account. It is found that the statistical uncertainty on γ is 8.5 • , the uncertainty due to strong-phase inputs is 4 • , and the uncertainty due to experimental systematic effects is 2 • .

Conclusions
Approximately 4100 (560) B ± → DK ± decays with the D meson decaying to K 0 S π + π − (K 0 S K + K − ) are selected from data corresponding to an integrated luminosity of 2.0 fb −1 collected with the LHCb detector in 2015 and 2016. These samples are analysed to determine the CP -violating parameters x ± ≡ r B cos(δ B ± γ) and y ± ≡ r B sin(δ B ± γ), where r B is the ratio of the absolute values of the B + → D 0 K − and B + → D 0 K − amplitudes, δ B is their strong-phase differences, and γ is an angle of the Unitarity Triangle. The analysis is performed in bins of the D-decay Dalitz plot and existing measurements performed by the CLEO collaboration [18] are used to provide input on the D-decay strong-phase parameters (c i , s i ). Such an approach allows the analysis to be free from model-dependent assumptions on the strong-phase variation across the Dalitz plot. This paper also gives the combination with the results obtained with an earlier data set, thereby allowing further improvements in the precision on γ. Considering only the data collected in 2015 and 2016 and choosing the solution that satisfies 0 < γ < 180 • yields r B = 0.086 +0.013 −0.014 , δ B = (101±11) • , and γ = (87 +11 −12 ) • . The values of r B and γ are consistent with world averages, while there is some tension in the determined value of δ B . This could be resolved by future analyses of the B → DK mode in a variety of D decays, including those analysed here, utilising the data set that is being collected with LHCb in 2017 and 2018. The measurement reported in this paper represents the most precise determination of γ from a single analysis.