Measurement of the CKM angle $\gamma$ in $B^\pm\to D K^\pm$ and $B^\pm \to D \pi^\pm$ decays with $D \to K_\mathrm S^0 h^+ h^-$

A measurement of $CP$-violating observables is performed using the decays $B^\pm\to D K^\pm$ and $B^\pm\to D \pi^\pm$, where the $D$ meson is reconstructed in one of the self-conjugate three-body final states $K_{\mathrm S}\pi^+\pi^-$ and $K_{\mathrm S}K^+K^-$ (commonly denoted $K_{\mathrm S} h^+h^-$). The decays are analysed in bins of the $D$-decay phase space, leading to a measurement that is independent of the modelling of the $D$-decay amplitude. The observables are interpreted in terms of the CKM angle $\gamma$. Using a data sample corresponding to an integrated luminosity of $9\,\text{fb}^{-1}$ collected in proton-proton collisions at centre-of-mass energies of $7$, $8$, and $13\,\text{TeV}$ with the LHCb experiment, $\gamma$ is measured to be $\left(68.7^{+5.2}_{-5.1}\right)^\circ$. The hadronic parameters $r_B^{DK}$, $r_B^{D\pi}$, $\delta_B^{DK}$, and $\delta_B^{D\pi}$, which are the ratios and strong-phase differences of the suppressed and favoured $B^\pm$ decays, are also reported.


Introduction
In the framework of the Standard Model, CP violation can be described by the angles and lengths of the Unitarity Triangle constructed from elements of the CKM matrix [1,2]. The angle γ ≡ arg (−V ud V * ub /V cd V * cb ), has particularly interesting features. It is the only CKM angle that can be measured in decays including only tree-level processes, and is experimentally accessible through the interference of b → cus and b → ucs (and CPconjugate) decay amplitudes. In addition, there are negligible theoretical uncertainties when interpreting the measured observables in terms of γ [3]. Hence, in the absence of unknown physics effects at tree level, a precision measurement of γ provides a Standard Model benchmark that can be compared with indirect determinations from other CKMmatrix observables more likely to be affected by physics beyond the Standard Model [4]. Such comparisons are currently limited by the precision of direct measurements of γ, which is about 5 • [5,6] dominated by LHCb results.
Decays such as B ± → DK ± , where D represents a superposition of D 0 and D 0 states, are used to observe the effects of interference between b → cus and b → ucs (and CPconjugate) decay amplitudes. The interference arises when the decay channel of the D meson is common to both D 0 and D 0 mesons. The B ± → DK ± decay has been studied extensively with a wide range of D-meson final states [7][8][9][10][11]. The exact choice of observables from each of these analyses is dependent on the method that is most appropriate for the D decay used [12][13][14][15][16][17][18][19][20]. The methods can be extended to a variety of different B-decay modes [8,[21][22][23][24][25].
This paper presents a model-independent study of the decay modes B ± → DK ± and B ± → Dπ ± where the chosen D decays are the self-conjugate decays D → K 0 S π + π − and D → K 0 The analysis of the B ± → DK ± , D → K 0 S h + h − decay chain is powerful due to the rich resonance structure of the D-decay modes, as has been described in Refs. [17][18][19]. The data used in this analysis were accumulated with the LHCb detector over the period 2011-2018 in pp collisions at energies of √ s =7, 8,13 TeV, corresponding to a total integrated luminosity of approximately 9 fb −1 . The presence of interference leads to differences in the phase-space distributions of D decays from reconstructed B + and B − decays. In order to interpret any observed difference in the context of the angle γ, knowledge of the strong phase of the D 0 decay amplitude, and how it varies over phase space, is required. An attractive model-independent approach makes use of direct measurements of the strong-phase difference between D 0 and D 0 decays, averaged over regions of the phase space [17,26,27]. Quantum correlated pairs of D mesons produced in decays of ψ(3770) give direct access to the strong-phase differences. These have been measured by the CLEO collaboration [28], and more recently the BESIII collaboration [29][30][31]. Measurements using the inputs in Ref. [28] have been used by the LHCb [10,21,32] and Belle [33,34] collaborations. An alternate method is to use an amplitude model of the D decay to determine the strong-phase variation [35][36][37]. 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. 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 violate the optical theorem [38]. 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. The analysis approach is laid out in Sect. 2, while Sect. 3 describes the LHCb detector used to collect the data sample, and Sect. 4 summarises the selection criteria. The measurement is based on a two-step fit procedure covered in Sect. 5, where the fit to the invariant-mass distribution is detailed, and Sect. 6, which describes how the CP observables are determined. The systematic uncertainties are reported in Sect. 7, and the results are interpreted to determine the value of γ in Section 8. Finally, the conclusions are presented in Sect. 9.

Analysis Overview
The sum of the favoured and suppressed contributions to the B − → DK − amplitude can be written as where The bins for which m 2 − > m 2 + are defined to have positive values of i 1 . The strong-phase difference between the D 0 -and D 0 -decay amplitudes at a given point on the Dalitz plot is denoted as δ D (m 2 − , m 2 + ). The cosine of δ D (m 2 − , m 2 + ) weighted by the D-decay amplitude and averaged over bin i is written as c i [17], and is given by where the integrals are evaluated over 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 expected yield of B − decays in bin i is found by integrating the square of the amplitude given in Eq. (1) over the region of phase space defined by the ith bin. The effects of charm mixing and CP violation are ignored, as is the presence of CP violation and matter regeneration in the neutral K 0 decays. These effects are expected to have a small impact [39,40] on the distribution of events on the Dalitz plot. Selection requirements and reconstruction effects lead to a non-uniform efficiency over phase space, denoted by η(m 2 − , m 2 + ). At LHCb the typical efficiency variation over phase space for a D → K 0 S h + h − decay from a region of high efficiency to low efficiency is approximately 60% [21]. The fractional yield of pure D 0 decays in bin i in the presence of this efficiency profile is denoted F i , given by where the sum in the denominator is over all Dalitz plot bins, indexed by j. Neglecting CP violation in these charm decays, the charge-conjugate amplitudes satisfy the relation where F i is the fractional yield of D 0 decays to bin i. The physics parameters of interest, r DK B , δ DK B , and γ, are translated into four CP -violating observables [41] that are measured in this analysis and are the real and imaginary parts of the ratio of the suppressed and favoured B decay amplitudes, Using the relations c i = c −i and s i = −s −i the B + (B − ) yields, N + (N − ), in bin i and −i are given by where h B + and h B − are normalisation constants. The value of r DK B is allowed to be different for each charge and is constructed from either (r DK B ) 2 = x DK The normalisation constants can be written as a function of γ, analogous to the global asymmetries studied in decays where the D meson decays to a CP eigenstate [8]. However, not only is this global asymmetry expected to be small since the CP -even content of the D → K 0 S π + π − and D → K 0 S K + K − decay modes is close to 0.5, it is also expected to be heavily biased due to the effects of K 0 S CP violation [40] on total yields. Therefore the global asymmetry is ignored and the loss of information is minimal. An advantage of this approach is that the normalisation constants h B + and h B − are independent of each other, and will implicitly contain the effects of the production asymmetry of B ± mesons in pp collisions and the detection asymmetries of the charged kaon from the B decay. This leads to a CP -violation measurement that is free of systematic uncertainties associated to these effects.
The system of equations provides 4N observables and 4+2N unknowns, assuming that the available measurements of c i and s i are used. This is solvable for N ≥ 2, but in practice the simultaneous fit of the F i , x DK ± , and y DK ± parameters leads to large uncertainties on the CP observables, and hence some external knowledge of the F i parameters is desirable. The F i parameters could be computed from simulation and an amplitude model, but the systematic uncertainties associated with the LHCb simulation would be significant. Recent analyses [10,32] have used the semileptonic decay B → D * µν, where the flavour-tagged yields of D 0 mesons are corrected for the differences in selection between the semileptonic channel and the signal mode. However, with the increased signal yields, the uncertainty due to this necessary correction will be approximately half the statistical uncertainty on the measurement presented in this paper, and therefore a different method is adopted.
The B ± → Dπ ± decay mode is expected to have F i parameters that are the same as those for B ± → DK ± if a similar selection is applied due to the common topology and the ability to use same signatures in the detector to select the candidates. The B ± → Dπ ± decay is expected to exhibit CP violation through the interference of b → cud and b → ucd transitions, analogous to the B ± → DK ± decay but suppressed by one order of magnitude [42]. Further effects from K 0 S CP violation and matter regeneration have been recently shown to have only a small impact on the distribution over the Dalitz plot [40], in contrast to their impact on the global asymmetry. Therefore the B ± → Dπ ± channel can be used to determine the F i parameters if the small level of CP violation in the B ± decay is accounted for.
Pseudoexperiments are performed in which the two B-decay modes are fit together assuming common F i parameters. Independent x ± and y ± observables are required for the two B decay modes due to different values of the hadronic parameters, r B and δ B . The value of r B in B ± → DK ± is approximately 0.1, and it is expected that it will be a factor 20 smaller in B ± → Dπ ± decays [42]. The yields of B ± → Dπ ± are described by a set of equations analogous to Eq. (5), with the substitutions x DK ± → x Dπ ± and y DK ± → y Dπ ± . An analysis that simultaneously measures the F i , x DK ± , y DK ± , x Dπ ± , and y Dπ ± parameters is found to be stable only if r Dπ B > 0.03. At the expected value r Dπ B = 0.005 the fit is unstable due to high correlations between the F i and x Dπ ± and y Dπ ± . Therefore an alternate parameterisation [43,44] is introduced, which utilises the fact that γ is a common parameter, and that the CP violation in B ± → Dπ ± decays can therefore be described by the addition of a single complex variable and in terms of x Dπ ξ ≡ Re(ξ Dπ ) and y Dπ ξ ≡ Im(ξ Dπ ), the x Dπ ± , y Dπ ± parameters are given by With this parameterisation, the simultaneous fit to x DK ± , y DK ± , x Dπ ξ , y Dπ ξ (the CP observables) and F i parameters is stable for all values of r Dπ B . The simultaneous fit of B ± → Dπ ± and B ± → DK ± candidates has two advantages. Firstly, the extraction of F i in this manner is expected to have negligible associated systematic uncertainty, and reduces significantly the reliance on simulation. Secondly, the CP -violating observables in B ± → Dπ ± using other D-decay modes [8,9] are not routinely included in the γ combination of all results because they allow for two solutions of r Dπ B , δ Dπ B , which makes the statistical interpretation of the full B ± → Dh ± combination problematic [45]. The measurement in the B ± → Dπ ± , D → K 0 S h + h − decays has the potential to resolve this redundancy, and allow for a more straightforward inclusion of all B ± → Dπ ± results in the combination. A small disadvantage is that the measurement of γ will incorporate information from both B ± → DK ± and B ± → Dπ ± decay modes and the contribution of each cannot be disentangled. However, since the size of contribution from the B ± → Dπ ± decay to the precision is expected to be negligible in comparison to that from the B ± → DK ± decay, this is considered an acceptable compromise. 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 scheme called the optimal binning, where the bins have been chosen to optimise the statistical sensitivity to γ, as described in Ref. [28]. The optimisation was performed assuming a strong-phase difference distribution as predicted by the BaBar model presented in Ref. [46]. 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. [47]. The D → K 0 S K + K − decay mode is dominated by the intermediate K 0 S φ and K 0 S a(980) states which are CP -odd and CP -even, respectively, and the narrow K 0 S φ resonance is encapsulated within the second bin of the 2 × 2 scheme. Therefore, most of the sensitivity is encompassed by this scheme, and the additional small gains from the more detailed schemes are offset by low yields and fit instabilities that arise when these bins are used. Therefore, the 2 × 2 bin is used for the analysis of the D → K 0 S K + K − decay mode. 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.
Measurements of the c i and s i parameters in the optimal binning scheme for the D → K 0 S π + π − decay and in the 2 × 2 binning scheme for the D → K 0 S K + K − decay are available from both the CLEO and BESIII collaborations. A combination of results from both collaborations is presented in Ref. [30] and Ref. [31] for the D → K 0 S π + π − and D → K 0 S K + K − decays, respectively. The combinations are used within this analysis.

LHCb Detector
The LHCb detector [48,49] 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 tracking system provides a measurement of the momentum, p, of charged particles with a relative uncertainty that varies from 0.5% at low momentum to 1.0% at 200 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 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. The events considered in the analysis are triggered at the hardware level when either one of the final-state tracks of the signal decay deposits enough energy in the calorimeter system, or when one of the other particles in the event, not reconstructed as part of the signal candidate, fulfils any trigger requirement. At the software stage, it is required that at least one particle should have high p T and high χ 2 IP , where χ 2 IP is defined as the difference in the primary vertex fit χ 2 with and without the inclusion of that particle. A multivariate algorithm [50] is used to identify secondary vertices consistent with being a two-, three-, or four-track b-hadron decay. The PVs are fitted with and without the B candidate tracks, and the PV that gives the smallest χ 2 IP is associated with the B candidate. Simulation is required to model the invariant-mass distributions of the signal and background contributions and determine the selection efficiencies of the background relative to the signal decay modes. It is also used to provide an approximation for the efficiency variations over the phase space of the D decay for systematic studies. In the simulation, pp collisions are generated using Pythia [51] with a specific LHCb configuration [52]. Decays of unstable particles are described by EvtGen [53], in which final-state radiation is generated using Photos [54]. The decays D → K 0 S π + π − and D → K 0 S K + K − are generated uniformly over phase space. The interaction of the generated particles with the detector, and its response, are implemented using the Geant4 toolkit [55] as described in Ref. [56]. With the exception of the signal decay, the simulated event is reused multiple times [57]. Some subdominant backgrounds are generated with a fast simulation [58] that can mimic the geometric acceptance and tracking efficiency of the LHCb detector as well as the dynamics of the decay.

Selection
The selection closely follows that of Ref. [10]. Decays of K 0 S → π + π − are reconstructed in two different ways: the first involving 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 that decay later such that track segments of the pions cannot be formed in the vertex detector. The first and second types of reconstructed K 0 S decays are referred to as long and downstream candidates, respectively. The long candidates have the best mass, momentum and vertex resolution, but approximately two-thirds of the signal candidates belong to the downstream category.
The D meson candidates are built by combining a K 0 S candidate with two tracks assigned either the pion or kaon hypothesis. A B candidate is then formed by combining the D meson candidate with a further track. At each stage of combination, selection requirements are placed to ensure good quality vertices, and K 0 S and D candidate invariantmasses are required to be close to their nominal mass [59]. Mutually exclusive particle identification (PID) requirements are placed on the companion track from the B decay to separate B ± → DK ± and B ± → Dπ ± candidates, where the companion refers to the final state π ± or K ± meson produced in the B ± → Dh ± decay. PID requirements are also placed on the charged decay products of the D meson to reduce combinatorial background. A series of selection requirements are placed on the candidates to remove background from other B meson decays. A background from B ± → Dh ± decays where the D meson decays to either π + π − π + π − or K + K − π + π − is rejected by requiring that the long K 0 S candidates decay a significant distance from the D vertex. Similarly, the D meson is required to have travelled a significant distance from the B vertex to suppress B decays with the same final state, but where there is no intermediate D meson decay. Semileptonic decays of the type D 0 → K * − l + ν, where charge-conjugate decays are implied, can be reconstructed as D → K 0 S h + h − with expected contamination rates of the order of a percent. To suppress electron to pion misidentification, a veto is placed on the pion from the D decay that has the opposite charge with respect to the companion particle, if the PID response suggests it is an electron. To suppress the similar muonic background, it is required that the charged track from the D decay has no corresponding activity in the muon detector. This veto also suppresses signal decays where the pion or kaon meson decays before reaching the muon detector. Therefore, it is applied on both charged tracks from the D decay, as these events have a worse resolution on the Dalitz plot, which is undesirable. Finally, the same requirement is placed on the companion track to suppress B → Dµν decays.
The large remaining combinatorial background is suppressed through the use of a boosted decision tree (BDT) [60,61] multivariate classifier. The BDT is trained on simulated signal events. The background training sample is obtained from the far upper sideband of the m(Dh ± ) mass distribution between 5800-7000 MeV/c 2 , in order to provide a sample independent from the data which will be used in the fit to determine the CP observables. A separate BDT is trained for B decays containing long or downstream K 0 S candidates. The input variables given to each BDT include momenta of the B, D, and companion particles, the absolute and relative positions of decay vertices, as well as parameters that quantify the fit quality in the reconstruction; the parameter set is identical to the one used in the previous LHCb measurement and listed in detail in Ref. [10]. The BDT has been proven not to bias the m(Dh ± ) distribution. A series of pseudoexperiments are run to find the threshold values for the two BDTs which provide the best sensitivity to γ. This requirement rejects approximately 98% of the combinatorial background that survives all other selection requirements, while having an efficiency of approximately 93% in simulated B ± → DK ± decays. The selection applied to B ± → DK ± and B ± → Dπ ± candidates is identical between the two decay modes with the exception of the PID requirement on the companion track.
A signal region is defined as within 30 MeV/c 2 of the B-meson mass [59]. The phasespace distributions for candidates in this range are shown in the Dalitz plots of Fig. 2 for B ± → DK ± candidates. The data are split by the final state of the D decay and by the charge of the B meson. Small differences between the phase-space distributions in B + → DK + and B − → DK − decays are visible in the K 0 S π + π − final state.

The DK and Dπ invariant-mass spectra
The analysis uses a two-stage strategy to determine the CP observables. First, an extended maximum-likelihood fit to the invariant-mass spectrum of all selected B ± candidates in LHCb LHCb LHCb LHCb The top (bottom) plots show data where the K 0 S candidate is long (downstream). A particle within square brackets in the legend denotes the particle that has not been reconstructed. the mass range 5080 to 5800 MeV/c 2 is performed, with no partition of the D phase space. This fit is referred to as the global fit. The global fit is used to determine the signal and background component parameterisations, which are subsequently used in a second stage where the data are split by B charge and partitioned into the Dalitz plot bins to determine the CP observables.
The invariant mass distributions of the selected B ± candidates are shown for Figs. 3 and 4, respectively, together with the results of the global fit superimposed. The invariant mass is kinematically constrained through a fit imposed on the full B ± decay chain [62]. The D and K 0 S candidates are constrained to their known masses [59] and the B ± candidate momentum vector is required to point towards the associated PV. The data sample is split into 8 categories depending on the reconstructed B decay, D decay mode, and K 0 S category, since the latter exhibits slightly different mass resolutions. The fit is performed simultaneously for all categories in order to allow parameters to be shared.
The peaks centered around 5280 MeV/c 2 correspond to the signal B ± → DK ± and B ± → Dπ ± candidates. The parameterisation for the signal invariant-mass shape is determined from simulation; the invariant-mass distribution is modelled with a sum of   LHCb LHCb The top (bottom) plots show data where the K 0 S candidate is long (downstream). A particle within square brackets in the legend denotes the particle that has not been reconstructed.
the probability density function (PDF) for a Gaussian distribution, f G (m|m B , σ), and a modified Gaussian PDF that is used to account for the radiative tail and the wider resolution of signal events that are poorly reconstructed. The modified Gaussian has the form which is Gaussian when ∆m 2 σ 2 /α L/R or ∆m 2 β −1 (with widths of σ and α L/R /β, respectively), with an exponential-like transition that is able to model the effect of the experimental resolution of LHCb. Thus, the signal PDF has the form The values of the tail parameters (α L , α R , β) and k are fixed from simulation and are common for the two D decays (which is possible due to the applied kinematic constraints) but different for each B decay and type of K 0 S candidate. The signal mass, m B , is determined in data and is the same for all categories. The width, σ, of the signal PDF is determined by the data and allowed to be different for each B decay and type of K 0 S candidate. The width is narrower in B ± → DK ± decays compared to B ± → Dπ ± decays due to the smaller free energy in the decay. The width is approximately 3% narrower in decays with long K 0 S candidates. The signal yield is determined in each of the categories where the candidates are reconstructed as B ± → Dπ ± . The signal yield in the corresponding category where the candidates are reconstructed as B ± → DK ± is determined by multiplying the B ± → Dπ ± yield by the parameters B × . The parameter B corresponds to the ratio of the branching fractions for B ± → DK ± and B ± → Dπ ± decays, while the correction factor, , takes into account the ratio of PID and selection efficiencies, and is determined for each pair of B ± → DK ± and B ± → Dπ ± categories. The parameter B is shared across all categories and is found to be consistent with Ref. [59].
To the right of the B ± → DK ± peak there is a visible contribution from B ± → Dπ ± decays that are reconstructed as B ± → DK ± decays. The corresponding contribution in the B ± → Dπ ± category is minimal due to the smaller branching fraction of B ± → DK ± , but is accounted for in the fit. The rates of these cross-feed backgrounds are fixed from PID efficiencies determined in calibration data, which is reweighted to match the momentum and pseudorapidity distributions of the companion track of the signal. A data-driven approach is used to determine the PDF of B ± → Dπ ± decays that are reconstructed as B ± → DK ± candidates, as described in Ref. [10]. The same procedure is implemented to determine the PDF of B ± → DK ± decays reconstructed as B ± → Dπ ± candidates.
The background observed at invariant masses smaller than the signal peak are candidates that originate from other B-meson decays where not all decay products have been reconstructed. Due to the selected invariant-mass range it is only necessary to consider B meson decays where a single photon or pion has not been reconstructed. This background type is split into three sources; the first where the candidate originates from a B ± or B 0 meson, referred to as partially reconstructed background, the second where the candidate originates from a B 0 s meson, and the third where the candidate originates from a B ± or B 0 and furthermore one of the reconstructed tracks is assigned the kaon hypothesis, when the true particle is a pion. The latter type of background appears in the B ± → DK ± candidates and is referred to as misidentified partially reconstructed background. The corresponding type of background is not modelled in the B ± → Dπ ± candidates, since it is suppressed due to the branching fractions involved and the majority is removed by the lower invariant-mass requirement.
There are contributions from B 0 → D * ± h ∓ and B ± → D * 0 h ± decays in all categories, where the pion or the photon originating from the D * meson is not reconstructed. The invariant-mass distributions of these decays depend on the spin and mass of the missing particle as described in Ref. [25]. The parameters of these shapes are determined from simulation, with the exception of a free parameter in the fit to characterise the resolution. The decays B ±,0 → Dπ ± π 0,∓ contribute to the B ± → Dπ ± candidates where one of the pions from the B decay is not reconstructed. The shape of this background is determined from simulated B ± → Dρ ± and B 0 → Dρ 0 decays. The decays B ± → DK ± π 0 and B 0 → DK + π − contribute to the B ± → DK ± candidates where the pion is not reconstructed. The invariant-mass distribution for these events is based on the amplitude model of B 0 → DK + π − decays [63]. The model is used to generate four-vectors of the decay products, which are smeared to account for the LHCb detector resolution.
The invariant mass is then calculated omitting the particle that is not reconstructed, and this distribution is subsequently fit to determine the fixed distribution for the fit. The same shape is used for the B ± → DK ± π 0 decay as the corresponding amplitude model is not available. Finally, the B ± → DK ± candidates also have a contribution from B 0 s → D 0 π + K − decays where the pion is not reconstructed. The shape of this contribution is determined in a similar manner to that of B 0 → DK + π − decays using the B 0 s → D 0 π + K − amplitude model determined in Ref. [64]. The yield of the partially reconstructed background is a floating parameter in each B ± → Dπ ± sample and related to the yield in the corresponding B ± → DK ± sample via the floating parameter B L and correction factors from PID and selection efficiencies. Analogously to the signal-yield parameterisation, B L is a single parameter, common to all categories, but in this case has no direct physical meaning. The relative yield of B ± → D * (→ D[γ])π ± and B 0 → D * (→ D[π ∓ ])π ± decays, where the particle within the square brackets is the one not reconstructed, are fixed from branching fractions [59], and selection efficiencies determined from simulation. The fractional yields of B ± → D * 0 (→ D[γ])π ± , and B ±,0 → D[π 0,π ∓ ]π ± decays are determined in the fit and are constrained to be the same for each B ± → Dπ ± sample. Due to the lower yields in the B ± → DK ± category and presence of additional backgrounds, the relative fractions of the various B ± and B 0 components are all fixed using information from branching fractions [59] and selection efficiencies from simulation. The yield of the B 0 s → D 0 π + K − decays is fixed relative to the yield of B ± → Dπ ± decays in the corresponding category using branching fractions [59], the fragmentation fraction [65], and relative selection efficiencies.
The shapes for the misidentified partially reconstructed backgrounds are determined from simulation, weighted by the PID efficiencies from calibration data. The yield of these backgrounds are determined from the partially reconstructed yields in the B ± → Dπ ± candidates, and the relative selection efficiencies, which include the PID efficiencies from calibration data and the selection efficiency due to requiring the reconstructed invariant mass to be above 5080 MeV/c 2 . The final component of background is combinatorial which is parameterised by an exponential function. The yield and slope of this background in each category are free parameters. The yields of the different signals and background types are integrated in the signal region 5249-5309 MeV/c 2 and reported in Table 1. The B ± → DK ± yields in categories of different D decay and type of K 0 S candidate have uncertainties that are smaller than their Poisson uncertainty since they are determined using the value of B, which is measured from all B ± → DK ± candidates.

CP observables
To determine the CP observables the data are divided into 16 categories (B decay, B charge, D decay, type of K 0 S candidate) and then further split into each Dalitz plot bin. A simultaneous fit to the invariant-mass distribution is performed in all categories and Dalitz plot bins. The mass shape parameters are all fixed from the global mass fit. The lower limit of the invariant mass is increased to 5150 MeV/c 2 to remove a large fraction of the partially reconstructed background. The composition of the remaining background is determined from the global fit described in Sect. 5. The signal yield in each bin is parameterised using Eq. (5) or the analogous set of expressions for B ± → Dπ ± . These equations are normalised such that the parameters h B ± represent the total observed signal Reconstructed as: yield in each category, and these are measured independently.
The parameters x DK ± , y DK ± , x Dπ ξ , and y Dπ ξ are free parameters in the fit and common to the K 0 S and D decay categories. The parameters c i and s i are fixed to those determined from the combination of BESIII and CLEO data in Ref. [30] for the D → K 0 S π + π − decays and in Ref. [31] for the D → K 0 S K + K − decays. The F i parameters for each D decay are determined in the fit; separate sets of F i parameters are determined for the two types of K 0 S candidates because the efficiency profile over the Dalitz plot differs between the K 0 S selections. Since the F i parameters must satisfy the constraints i F i = 1, F i ∈ [0, 1], the fit can suffer from instability if they are included in a naive way due to large correlations. Therefore, the F i parameters are reparameterised as a series of recursive fractions with parameters, R i , determined in the fit. The relation between the F i and R i parameters is given by for a binning scheme with 2 × N bins. The yield of the combinatorial background in each bin is a free parameter. The yield of the partially reconstructed background from B ± or B 0 decays in the B ± → Dπ ± and B ± → DK ± samples is also a free parameter in each bin. The yield of the misidentified partially reconstructed background in the B ± → DK ± samples is determined via the background yield in the corresponding B ± → Dπ ± bin and the relative PID and selection efficiencies. The yield of the B 0 s → D 0 K − π + background is fixed from the global fit and is divided into the Dalitz plot bins according to the F i such that it has the distribution of a D 0 decay in the B + categories and the distribution of a D 0 decay in the B − categories.
There is a small fraction of bins where either the partially reconstructed background or combinatoric background yield is less than one. These bins are identified in a preliminary fit and the background yield is fixed to zero. This procedure is carried out to improve the fit stability.
Pseudoexperiments are performed to investigate any potential biases or remaining instabilities in the fit. The candidate yields and mass distributions in these pseudoexperiments are based on the global fit results. The pull distributions are well described for B − and B + decays. The signature for CP violation is that these vectors must have non-zero length and have a non-zero opening angle between them, since this angle is equal to 2γ, as illustrated on the figure. Therefore, the data exhibit unambiguous features of CP violation as expected. The relation between the hadronic parameters in B ± → Dπ ± and B ± → DK ± decays is also illustrated in Fig 5, where the vector defined by the coordinates (x Dπ ξ ,y Dπ ξ ) is the relative magnitude of r B between the two decay modes. It is consistent with the expectation of 5% [42].
A series of cross checks is carried out by performing separate fits by splitting the data sample into data-taking periods by year, type of K 0 S candidate, D-decay, and magnet polarity. The results are consistent between the datasets. As an additional cross check, the two-stage fit procedure is repeated with a number of different selections applied to the data. Of particular interest are the alternative selections that significantly affect the presence of specific backgrounds: the fits where the value of the BDT threshold value is varied to decrease the level of combinatorial background and those where the choice of PID selection is changed to result in a substantially lower level of misidentified B ± → Dπ ± decays and misidentified partially reconstructed background in the B ± → DK ± candidates. The variations in the central values for the CP observables are consistent within the statistical uncertainty associated with the change in the data sample.
In order to assess the goodness of fit and to demonstrate that the equations involving the CP parameters provide a good description of the signal yields in data, an alternative fit is performed where the signal yield in each B ± → DK ± and B ± → Dπ ± bin is measured independently. These yields are compared with those predicted from the values of   (x DK ± , y DK ± ) in the default fit and a high level of agreement is found. In order to visualise the observed CP violation, the asymmetry, , is computed for effective bin pairs, defined to comprise bin i for a B + decay and bin −i for a B − decay. Figure 6 shows the obtained asymmetries and those predicted by the values of the CP observables obtained in the fit. A further fit that does not allow for CP violation is carried out by imposing the conditions x DK This determines the predicted asymmetry arising from detector and production effects. In the B ± → DK ± sample the CP violation is clearly visible as the data are inconsistent with the CP -conserved hypothesis. The predicted asymmetries in the B ± → Dπ ± decay are an order of magnitude smaller. The data in this analysis cannot distinguish between the CP -violating and CP -conserving predictions for B ± → Dπ ± due to the relatively large statistical uncertainties.

Systematic uncertainties
Systematic uncertainties on the measurements of the CP observables are evaluated and are presented in Table 2. The limited precision on (c i , s i ) coming from the combined BESIII and CLEO [30,31] results induces uncertainties on the CP parameters. These uncertainties are evaluated by fitting the data multiple times, each time with different (c i , s i ) values sampled according to their experimental uncertainties and correlations. 2 The resulting standard deviation of each distribution of the CP observables is assigned as the systematic uncertainty. The size of the systematic uncertainty is notably much smaller than the corresponding uncertainty in Ref. [10] due to the improvement in the knowledge of these strong-phase parameters [30,31].
The non-uniform efficiency profile over the Dalitz plot means that the values of (c i ,s i ) appropriate for this analysis can differ from those measured in Refs. [30,31], which correspond to the case where there is no variation in efficiency over the Dalitz plot. Amplitude models from Refs. [47,66]  The assumption that the relative variation of efficiency over the Dalitz plot is the same in selected B ± → DK ± and B ± → Dπ ± candidates is verified in simulated samples of similar size to the B ± → Dπ ± yields observed in data. No statistically significant difference is observed and no systematic uncertainty is assigned. The uncertainties from the fixed invariant-mass shapes determined in the global fit are propagated to the CP observables through a resampling method [67]. The following procedure, which takes into account the fact that some parameters are determined in simulation and others in data, is carried out a hundred times. First, the simulated decays that were used to determine the nominal mass shape parameters are each resampled with replacement and fit to determine an alternative set of parameters. Then, the final dataset is resampled with replacement and the global fit is repeated using the alternative fixed shape parameters, to determine alternative values for the parameters that are determined from real data. Finally, the CP fit is performed using the alternative invariant-mass parameterisations, without resampling the final dataset. The standard deviation of the CP observables obtained via this procedure is taken as the systematic uncertainty due to the fixed parameterisation.
The PID efficiencies are varied within their uncertainties in the global and CP fit and the standard deviation of the CP parameters is taken as the systematic uncertainty. A similar method is used to determine the uncertainties due to the fixed fractions between different partially reconstructed backgrounds where the uncertainties on the fixed fractions are those from the branching fractions [59] and the selection efficiencies.
The CP fit assumes the same mass shape for each component in each Dalitz plot bin. For the signal and cross-feed backgrounds the shapes are redetermined in each bin using the same procedures described in Sect. 5. The variance is very small due to weak correlations between phase-space coordinates and particle kinematics. The combinatorial slope can also vary from bin to bin, as the relative rate of combinatorial background with and without a real D 0 meson will not be constant. The size of this effect is determined through the study of the high invariant-mass sideband where only combinatorial background contributes. Pseudodata are generated where this variation in mass shape across the Dalitz plot bins is replicated for signal, cross-feed and combinatorial backgrounds, and the generated samples are fit with the default fit assumptions of the same shape in each bin. The mean bias is assigned as the systematic uncertainty.
The partially reconstructed background shape is also expected to vary in each bin, however the leading source of this effect is due to the individual components of this background having a different distribution over the Dalitz plot. Some partially reconstructed backgrounds will be distributed as D 0 (D 0 ) → D → K 0 S h + h − for reconstructed B − (B + ) candidates, while others will be distributed as a D 0 -D 0 admixture depending on the relevant CP -violation parameters. Pseudodata are generated, where the D-decay phase-space distributions for B ± → D * K ± and B ± → DK * + background events are based on the CP parameters reported in Ref. [68]. No CP violation is introduced into the partially reconstructed background in the B ± → Dπ ± samples since it is expected to be small, and the B 0 → Dρ 0 background is treated as an equal mix of D 0 and D 0 since either pion can be reconstructed. The generated pseudodata are fitted with the default fit and the mean bias is assigned as the systematic uncertainty.
Systematic uncertainties are assigned for small residual backgrounds that contaminate the data sample but are not accounted for in the fit. Their impact is assessed by generating pseudoexperiments that contain these backgrounds and are fit with the default model. The mean bias is assigned as the uncertainty. One source of background is from Λ 0 b → Dpπ − decays where the pion is not reconstructed and the proton is misidentified as a kaon. This background is modelled as a D 0 -like contribution in B − decays, and has an expected yield of 0.5% of the B ± → DK ± signal. A further, even smaller, background is Λ 0 b → Λ + c (→ pK 0 S π + π − )π − decays where the π + meson in the Λ + c decay is missed, and the p reconstructed as the π + from the D-decay. The effective distribution of the reconstructed D meson is unknown and is assigned to be D 0 -like in B − decays to be conservative. The mass shapes and rates of these backgrounds are determined from simulation. Another source of background comes from residual B → Dµν decays, where the rate (less than 0.2 % relative to the signal mode, after the applied veto) and shape are determined from simulation with PID efficiencies from calibration data. The residual semileptonic D decay background has a rate of less than 0.1% of signal and the distribution of these events on the Dalitz plot is determined through a simplified simulation [58] taking into account various K * mesons. Finally, a small peaking background from B ± → D(→ K ± π ∓ )K 0 S π ± decays where the kaon is reconstructed as the companion and the other particles are assigned to the D decay is considered. The yield of this background is determined to be 0.5% of the signal yield in B ± → DK ± by a data driven study of the invariant-mass distribution of switched tracks. The distribution on the Dalitz plot is determined through the simplified simulation [58] where different K * ± → K 0 S π ± resonances are generated.
The main effect of migration from one Dalitz plot bin to another is implicitly taken into account by using the data to determine the F i , which thus include the effects of the net bin migration. However, a small effect arises because of the differences in the distributions of the B ± → DK ± and B ± → Dπ ± decays due to the differing hadronic decay parameters. To investigate this, data points are generated according to the amplitude model in Ref. [66] with CP observables consistent with expectation [5,68]. To smear these data points on the Dalitz plot, an event is selected from full LHCb simulation and the difference in m 2 + and m 2 − between its true and reconstructed quantities is applied to the data point in order to determine its reconstructed bin. The difference between true and reconstructed quantities is multiplied by a factor of 1.2 to account for differences in resolution between data and simulation. Pseudoexperiments are generated based on the expected reconstructed yields in each bin and fit with a nominal fit where the c i and s i parameters are determined by the amplitude model [66]. The mean bias in the CP violation parameters is taken as the systematic uncertainty, which is small.
The impact of ignoring the CP violation and matter effects in K 0 S decays is determined through generating pseudoexperiments taking into account all these effects as detailed in Ref. [40], where LHCb simulation is used to obtain the K 0 S lifetime acceptance and momentum distribution. The size of the bias found is consistent with those expected from Ref. [40], where it was also predicted that the relative uncertainties on B ± → Dπ ± observables are be expected to be larger than for B ± → DK ± observables. This is found to be true, but even the most significant uncertainty, on y Dπ ξ , is an order of magnitude smaller than the corresponding statistical uncertainty. The effect of ignoring charm mixing is expected to be minimal, given that the first-order effects are inherently taken into account when the F i parameters are measured as a part of the fit [39]. This is verified by generating pseudoexperiments that include charm mixing and fitting them with the nominal fit.
In previous studies, a bias correction has been necessary when similar measurements have been performed with lower signal yields [10] leading to some fit instabilities. In this case, the higher yields have resulted in a bias that is of negligible size and hence no correction is applied. Nonetheless, the uncertainty on the biases are assigned as the systematic uncertainties.
In general, all the systematic uncertainties are small in comparison to the statistical uncertainties. There is no dominant source of systematic uncertainty for all CP observables, however the description of backgrounds, either those not modelled or the modelling of the partially reconstructed backgrounds are some of the larger sources. The uncertainty attributed to the precision of the strong-phase measurements is of similar size to the total LHCb-related systematic uncertainty.

Interpretation
The CP observables are measured to be where the first uncertainty is statistical, the second arises from systematic effects in the method or detector considerations, and the third from external inputs of strong-phase measurements from the combination of CLEO and BESIII [28,30] results. The correlation matrices for each source of uncertainty are available in the appendices in Tables 3-5.
The CP observables are interpreted in terms of the underlying physics parameters γ, and r B and δ B for each B ± decay mode. The interpretation is done via a maximum likelihood fit using a frequentist treatment as described in Ref. [45]. The solution for the physics parameters has a two-fold ambiguity as the equations are invariant under the simultaneous substitutions γ → γ + 180 • and δ B → δ B + 180 • . The solution that satisfies 0 < γ < 180 • is chosen, and leads to Pseudoexperiments are carried out to confirm that the value of γ is extracted without bias. This is the most precise single measurement of γ to date. The result is consistent with the indirect determination γ = 65.66 +0.90 are consistent with their current world averages [5,6] which include the LHCb results obtained with the 2011-2016 data. The knowledge of r Dπ B and δ Dπ B from other sources is limited, with the combination of many observables presented in Ref. [45] providing two possible solutions. The results here have a single solution, and favour a central value that is consistent with the expectation for r Dπ B , given the value of r DK B and CKM elements [42]. This is likely to remove the two-solution aspect in future combinations of γ and associated hadronic parameters. The low value of r Dπ B means that the direct contribution to γ from B ± → Dπ ± decays in this measurement is minimal. However the ability to use this decay mode to determine the efficiency has approximately halved the total LHCb related experimental systematic uncertainty in comparison to Ref. [10]. The new inputs from the BESIII collaboration have led to the strong-phase related uncertainty on γ to be approximately 1 • , which is a significant reduction compared to the propagated uncertainty when only CLEO measurements were available.

Conclusions
In summary, the decays B ± → DK ± and B ± → Dπ ± with D → K 0 S π + π − or D → K 0 S K + K − obtained from the full LHCb dataset collected to date, corresponding to an integrated luminosity of 9 fb −1 , have been analysed to determine the CKM angle γ. The sensitivity to γ comes almost entirely from B ± → DK ± decays where the signal yields of reconstructed events are approximately 13600 (1900) in the D → K 0 S π + π − (D → K 0 S K + K − ) decay modes. The B ± → Dπ ± data is primarily used to control effects due to selection and reconstruction of the data, which leads to small experimental systematic uncertainties. The analysis is performed in bins of the D-decay Dalitz plot and a combination of measurements performed by the CLEO and BESIII collaborations presented in Refs. [30,31] 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. The analysis also determines the hadronic parameters r B and δ B for each B ± decay mode. Those of the B ± → DK ± decay are consistent with current averages, and those of the B ± → Dπ ± decay are obtained with the best precision to date, and have not previously been measured using these D-decay modes. The CKM angle γ is determined to be γ = (68.7 +5. 2 −5.1 ) • , where the result is limited by statistical uncertainties. This is the most precise measurement of γ from a single analysis, and supersedes the results in Refs. [10,32].

Appendices A Correlation matrices
The correlations matrices for the measured observables are shown in Tables 3-5 for the statistical uncertainties, the experimental systematic uncertainties, and the strong-phaserelated uncertainties, respectively. Correlation matrix x DK