Improved measurement of the CKM angle gamma in B^-+ ->D^(*) K^(*)-+ decays with a Dalitz plot analysis of D decays to K_S^0 pi^+ pi^- and K_S^0 K^+ K^-

We report on an improved measurement of the Cabibbo-Kobayashi-Maskawa {\it CP}-violating phase $\gamma$ through a Dalitz plot analysis of neutral $D$ meson decays to $K_S^0 \pi^+ \pi^-$ and $K_S^0 K^+ K^-$ in the processes $B^\mp \to D K^\mp$, $B^\mp \to D^* K^\mp$ with $D^* \to D\pi^0,D\gamma$, and $B^\mp \to D K^{*\mp}$ with $K^{*\mp} \to K_S^0 \pi^\mp$. Using a sample of 383 million $B\bar{B}$ pairs collected by the BABAR detector, we measure $\gamma=(76 \pm 22 \pm 5 \pm 5)^\circ$ (mod $180^\circ$), where the first error is statistical, the second is the experimental systematic uncertainty and the third reflects the uncertainty on the description of the Dalitz plot distributions. The corresponding two standard deviation region is $29^\circ<\gamma<122^\circ$. This result has a significance of direct {\it CP} violation ($\gamma \ne 0$) of 3.0 standard deviations.


I. INTRODUCTION AND OVERVIEW
In the Standard Model (SM) the phase in the Cabibbo-Kobayashi-Maskawa (CKM) quark-mixing matrix [1] is the sole source of CP violation in the quark sector of the electroweak interactions. This phase can be directly determined using a variety of methods, involving either the interference between decays with and without mixing in time-dependent CP asymmetries in neutral B meson decays, or interference between neutral B (self-tagged) or charged B decays yielding the same final state (direct CP violation). These multiple determinations in CP -violating tree-level processes as well as in decays involving penguin diagrams test the CKM mechanism, thus probing the presence of physics beyond the SM [2]. Among these determinations, the measurement of the angle γ, defined as arg [−V ud V * ub /V cd V * cb ], where V ij are the elements of the CKM matrix, is one of the most difficult to achieve and constitutes an important goal of present and future B physics experiments. Several methods have been proposed to extract γ. However, those using B ∓ →D ( * )0 K ( * )∓ decays [3,4] (the symbolD ( * )0 indicates either a D ( * )0 or a D ( * )0 meson) are theoretically clean and are unlikely to be affected by new physics because the main contributions to the amplitudes come from tree-level diagrams, as shown in Fig. 1. This is an important distinction from most of other direct measurements of phases of CKM elements. The decay amplitudes for the color allowed B − → D ( * )0 K ( * )− (b → cus) and the color suppressed B − → D ( * )0 K ( * )− (b → ucs) transitions [5] differ by a factor r ( * ) B e i(δ ( * ) B ∓γ) . Here, r ( * ) B is the magnitude of the ratio of the amplitudes A(B − → D ( * )0 K − ) and A(B − → D ( * )0 K − ) and δ ( * ) B is their relative strong phase. The weak phase γ leads to different B − and B + decay rates (direct CP violation) and, when the D ( * )0 and D ( * )0 decay to a common final state [6][7][8][9], the phases become observable. The uncertainty in γ scales roughly as 1/r ( * ) B . From the ratio of CKM matrix elements we expect r ( * ) B ≈ c F | V cs V * ub | / | V us V * cb | to be approximately in the range 0.1−0.2, where c F ∼ 0.2−0.4 is the color suppression factor [10,11].
Here, A c (p) and A u (p) are the magnitudes of the b → c and b → u amplitudes as a function of the B ∓ → D 0 K 0 S π ∓ phase space position p, and δ(p) is the relative strong phase. The parameter κ accounts for the interference between B ∓ →D 0 K * ∓ and other B ∓ →D 0 K 0 S π ∓ amplitudes with 0 < κ < 1 in the most general case. This effective parameterization also accounts for efficiency variations as a function of the kinematics of the B decay.
In this paper we present an improved measurement of γ based on the analysis of the Dalitz plot distribution of D 0 → K 0 S π + π − and, for the first time,D 0 → K 0 S K + K − , using a sample of 351 fb −1 of integrated luminosity recorded at the Υ (4S) resonance, corresponding to 383 million BB pairs. We analyze a total of seven signal samples (also referred to as CP samples), B − → D ( * )0 K − and B − →D 0 K * − , withD * 0 →D 0 π 0 ,D 0 γ, [5]. Due to lack of statistics, the decay B − →D 0 K * − with D 0 → K 0 S K + K − has been excluded from the analysis. We also reconstruct high statistics control samples, one for each signal B decay channel: B − → D ( * )0 π − and (for [5,19]. We exploit the same dataset to determine A D∓ for D 0 → K 0 S π + π − and D 0 → K 0 S K + K − decays from analyses of the respective Dalitz plots for high-statistics samples of flavor-tagged D 0 mesons from D * + → D 0 π + decays [5], produced in e + e − → cc events. Additional improvements compared to our previous publication [18] include a higher reconstruction efficiency, an optimized treatment of the e + e − → qq, q = u, d, s, c background and an improved D 0 → K 0 S π + π − description of the Dalitz plot distribution (referred to hereafter as Dalitz model), resulting in significant decrease of statistical, systematic and model uncertainties. This measurement supersedes our previous result based on 227 million BB pairs [18].
The paper is organized as follows. In Sec. II we describe the reconstruction and selection of the signal and control samples. Sec. III is devoted to the determination of A D∓ for D 0 → K 0 S π + π − and D 0 → K 0 S K + K − decays. In Sec. IV we describe the simultaneous maximum likelihood fit to the distributions Γ ∓ , x s∓ , and y s∓ . In that section we also present the experimental results, including systematic uncertainties. We extract these CP parameters since they have a good Gaussian behavior for small values of r B , κr s and δ s , using a statistical (frequentist) analysis.

II.
EVENT SELECTION

A. BABAR detector
This analysis is based on a data sample collected by the BABAR detector at the Stanford Linear Accelerator Center PEP-II e + e − asymmetric-energy storage ring. The BABAR detector is described in detail elsewhere [20]. We summarize briefly the components that are crucial to this analysis. Charged-particle tracking is provided by a fivelayer silicon vertex tracker (SVT) and a 40-layer drift chamber (DCH). In addition to providing precise space coordinates for tracking, the SVT and DCH also measure the specific ionization (dE/dx), which is used for particle identification of low-momentum charged particles. At higher momenta (p > 0.7 GeV/c) pions and kaons are identified by Cherenkov radiation detected in a ring-imaging device (DIRC). The typical separation between pions and kaons varies from 8σ at 2 GeV/c to 2.5σ at 4 GeV/c, where σ denotes here the standard deviation. The position and energy of photons are measured with an electromagnetic calorimeter (EMC) consisting of 6580 thallium-doped CsI crystals. These systems are mounted inside a 1.5-T solenoidal super-conducting magnet. We use a GEANT4-based Monte Carlo (MC) simulation to model the response of the detector, taking into account the varying accelerator and detector conditions, and to generate large samples of signal and background for the CP and control modes considered in the analysis.

B. Event reconstruction and selection
The B − candidates are formed by combining aD ( * )0 candidate with a track identified as a kaon [20] or with a K * − candidate formed as a combination of a K 0 S and a negatively charged pion, with an invariant mass within 55 MeV/c 2 of the nominal K * − mass. Here and in the following, nominal mass values are taken from [21]. Thẽ D 0 candidates are selected by requiring the K 0 S π + π − or K 0 S K + K − invariant mass to be within 12 MeV/c 2 of the nominal D 0 mass, and the momentum in the center-ofmass (CM) frame to be greater than 1.3 GeV/c. The K ∓ tracks inD 0 → K 0 S K + K − are required to be positively identified as kaons in the DCH and DIRC. The π 0 candidates fromD * 0 →D 0 π 0 decays are formed from pairs of photons with invariant mass in the range [115,150] MeV/c 2 , and with photon energy greater than 30 MeV. Photon candidates fromD * 0 →D 0 γ decays are selected if their energy is greater than 100 MeV. TheD 0 candidates are combined with a π 0 (γ) to form theD * 0 candidate, and are required to have aD * 0 -D 0 mass difference within 2.5 (10) MeV/c 2 of its nominal value. The K 0 S candidates are reconstructed from pairs of oppositelycharged pions constrained to originate from the same point and with an invariant mass within 9 MeV/c 2 of the K 0 S nominal mass. The cosine of the collinearity an-gle between the K 0 S momentum and the line connecting its parent particle (theD 0 or the K * − ) and the K 0 S decay points in the plane transverse to the beam, is required to be larger than 0.990 (0.997 for K 0 S from K * − decays). This cut helps to significantly reduce background contributions fromD 0 → ππππ decays and from a − 1 misreconstructed as K * − . The kinematic variables of thẽ D ( * )0 , K 0 S , and π 0 when forming the B − ,D 0 , andD * 0 , respectively, are fitted with masses constrained to nominal values. For B − →D 0 K * − decays we also require | cos θ H | ≥ 0.35, where θ H is the angle between the momentum of the K * − daughter pion and the parent B − in the K * − rest frame. The distribution of cos θ H is proportional to cos 2 θ H for B − →D 0 K * − while it is approximately flat for e + e − → qq, q = u, d, s, c (continuum) background.
We characterize B mesons using two almost independent variables, the beam-energy substituted mass, where the subscripts 0 and B refer to the initial e + e − system and the B candidate, respectively, and the asterisk denotes the CM frame. The signal events peak at the B mass in m ES and at zero in ∆E. The m ES resolution is about 2.6 MeV/c 2 and does not depend on the decay mode or on the nature of the prompt particle (K − or K * − candidate). In contrast, the ∆E resolution depends on the momentum resolution of the D ( * )0 meson and the prompt particle, and ranges between 15 MeV and 18 MeV, depending on the decay mode. We select events with m ES > 5.2 GeV/c 2 and −80 MeV < ∆E < 120 MeV. We discriminate against the main background contribution coming from continuum events through the fit to the data, as described in Sec. II C.
For events in which multiple B candidates satisfy the selection criteria, the one whose measuredD 0 mass differs from the nominal value by the least number of standard deviations, is accepted as signal candidate. For B − →D 0 K * − decays we select the candidate with the smallest value for the sum of the squares of the differences from nominal values, in standard deviations, of both K * − andD 0 masses. The fraction of events in which we reconstruct more than one candidate is less than 1% for B − →D ( * )0 K − samples and about 6% for candidates are selected in the same event, only the B − →D * 0 [D 0 π 0 ]K − is kept. This contamination has a negligible effect on the measurement of the CP parameters. Figure 2 shows the m ES distributions in the ∆E signal region defined through the requirement |∆E| < 30 MeV, after all selection criteria are applied. The reconstruction efficiencies are 20%, 9%, 12%, and 12%, for and B − →D 0 K * − decay modes, respectively, with D 0 → K 0 S π + π − . Similarly, forD 0 → K 0 S K + K − channels we obtain 19%, 8%, and 11% (B − →D 0 K * − is not reconstructed).
The same selection criteria are applied to select B − → D ( * )0 π − control samples, apart from the particle identification requirement on the prompt track, which is replaced by a kaon identification veto. Since we evaluate ∆E with the kaon mass hypothesis, the ∆E distributions are shifted by approximately +50 MeV, as given by Eq. (4). B − → D 0 a − 1 candidates are reconstructed similarly to B − →D 0 K * − , with a − 1 → ρ(770) 0 π − candidates made using combinations of three charged tracks with the requirements that the a − 1 invariant mass must be in the range [1.0, 1.6] GeV/c 2 , and that of the ρ(770) 0 within 150 MeV/c 2 of its nominal mass. As for signal samples, we select events with −80 MeV < ∆E < 120 MeV. Figure 3 shows the m ES distributions for all control samples in the ∆E signal region defined through the requirement 20 MeV < ∆E < 80 MeV, after all selection criteria. The corresponding reconstruction efficiencies are similar to those estimated for the CP samples.

C. Background composition and signal yields
The largest background contribution is from continuum events, where a fake or trueD ( * )0 is combined with a random track (B − →D ( * )0 K − samples), or a fake or trueD 0 is combined with a random or fake K * (B − →D 0 K * − sample). To separate continuum from BB events in a likelihood fit (discussed below and in Sec. IV), variables that characterize the event shape are used. We construct a Fisher discriminant F [22] from a linear combination of four topological variables: the monomials L 0 = i p * i and L 2 = i p * i | cos θ * i | 2 , | cos θ * T |, and | cos θ * B |. Here, p * i and θ * i are the CM momentum and the angle of the remaining tracks and clusters in the event, with respect to the B candidate thrust axis [23]. θ * T is the angle between the thrust axis of the B candidate and that of the rest of the event, and θ * B is the polar angle of the B candidate momentum, in the CM frame. The first three variables account for the jet-like shape of continuum events, in comparison to the spherical topology of BB events. In particular, the variable | cos θ * T | peaks close to one for continuum while for BB it is essentially uniformly distributed. The angular distribution of the variable | cos θ * B | follows 1 − cos 2 θ * B for BB events and 1+cos 2 θ * B for e + e − → qq [24]. The strategy of using the Fisher discriminant in a likelihood fit enhances significantly, typically about 25%, the signal reconstruction efficiency compared to our previous analysis [18], where we required | cos θ * T | < 0.8. At the same time it provides a larger sample of continuum events in the m ES sidebands, thus allowing the determination of the background properties directly from data (see Sec. IV).
Another source of background is related to BB decays where a fake or trueD ( * )0 is combined with a random or The curves superimposed represent the projections of the fit described in Sec. II C: signal plus background (solid black lines), the continuum plus BB background contributions (dotted red lines), and the sum of the continuum, BB, and K/π misidentification background components (dashed blue lines). misidentified track or K * . The main single contribution for B − →D ( * )0 K − signal comes from B − → D ( * )0 π − decays when the pion is misidentified as a kaon. This source is accounted for separately from other BB backgrounds. This contribution can be discriminated from the signal due to the shift in the ∆E distribution relative to that of signal events. Since ∆E is computed by The curves superimposed represent the projections of the fit described in Sec. II C: signal plus background (solid black lines), the continuum plus BB background contributions (dotted red lines), and the sum of the continuum, BB, and K/π misidentification background components (dashed blue lines).
assigning the kaon mass hypothesis to the prompt track, it is shifted by a quantity which depends on the momentum p of the prompt track in the laboratory frame and the Lorentz parameter γ PEP−II characterizing the boost of the CM relative to the laboratory frame, estimated from the PEP-II beam energies. For B − →D 0 K * − signal the main BB background source comes from B − → D 0 a − 1 decays. Since this contribution is highly suppressed by the cut on the cosine of the collinearity angle at 0.997, it is not treated separately from other BB backgrounds. Non-K * decays contributing to the B − →D 0 K * − sample are considered as signal, and their effect is accounted for by the factor κ defined in Eq. (3).
We fit the seven signal samples B − →D ( * )0 K − and B − →D 0 K * − , and their control samples B − → D ( * )0 π − and B − → D 0 a − 1 , using an unbinned extended maximum likelihood method to extract signal and background yields, and probability density functions (PDFs) for the variables m ES , ∆E, and the F discriminant, in the ∆E selection region. Three different background components are considered: continuum events, K/π misidentification (for B − →D ( * )0 K − and B − → D ( * )0 π − samples only), and other Υ (4S) → BB decays. The log-likelihood for each of the CP and control samples is where The signal m ES distributions for each CP sample and its corresponding control sample are parameterized using a common single Gaussian. Similarly, the F PDF makes use of a double Gaussian with different widths for the left and right parts of the curve (bifurcated Gaussian), and is assumed common for all CP and control samples. The signal ∆E distribution for B − →D ( * )0 K − events is parameterized with a double Gaussian function, while for B − → D ( * )0 π − events we use the same function, shifted event-by-event using Eq. (4). For B − →D 0 K * − and B − →D 0 a − 1 signal events a common double Gaussian is used instead.
The continuum background in the m ES distribution is described by a threshold function [25] while the continuum ∆E distribution is described using a first order polynomial parameterization. The free parameters are different for each CP sample but common to the corresponding control sample. The F distribution for continuum background is parameterized with the sum of two Gaussian functions and assumed common for all samples.
The shape of the m ES distribution for Υ (4S) → BB background (excluding the K/π misidentification contribution, as indicated previously) is taken from generic BB simulated events for each CP and control sample independently, and uses a threshold function [25] to describe the combinatorial component plus a bifurcated Gaussian to parameterize the contribution peaking at the B mass. The fraction of the peaking contribution is ex-tracted directly from the fit to the data, except for thẽ D 0 → K 0 S K + K − CP samples, where it is taken from the generic BB MC due to lack of statistics. The ∆E distribution for BB background is taken similarly from simulation and is parameterized with the sum of a second order polynomial and an exponential function that takes into account the increase of combinatorial background at negative ∆E values. A Gaussian function is also included to account for potential ∆E peaking background, although we find no significant peaking structure in any of our samples. The F distributions for signal and control samples, and for the generic BB background, are assumed to be the same as found in the simulation.
decay amplitudes are determined from Dalitz plot analyses of D 0 mesons from D * + → D 0 π + decays produced in e + e − → cc events. The charge of the low momentum π + from the D * + decay identifies ("tags") the flavor of the D 0 . Reconstruction and selection of D 0 → K 0 S h + h − candidates from D * + → D 0 π + decays are similar to those from B − →D ( * )0 K − decays, the only exception being the kaon identification of only one charged kaon for The D * + candidates are formed by combining the D 0 with the low momentum charged track. The two D * + decaying daughters are constrained to originate from the same point inside the PEP-II luminous region. To reduce combinatorial background and contamination from BB decays, the D 0 candidates are required to have a CM momentum greater than 2.2 GeV/c.
Each D 0 sample is characterized by the distributions of two variables, the invariant mass of the D 0 candidate m D and the ∆m = D * + − D 0 mass difference. We select D 0 candidates within ±0.64 MeV/c 2 and ±0.61 MeV/c 2 , corresponding to ±2 standard deviations, around the nominal ∆m [21], for D 0 → K 0 S π + π − and D 0 → K 0 S K + K − , respectively. Figure 4 shows the resulting D 0 → K 0 S h + h − mass distributions. The m D lineshape is described using a two Gaussian function for the signal and a linear background, as also shown in Fig. 4. The m D resolutions are 6.7 MeV/c 2 and 3.9 MeV/c 2 , for D 0 → K 0 S π + π − and D 0 → K 0 S K + K − . The mass resolution for the latter is better than that of the former because of the much smaller Q-value involved. The signal purity in the signal box (±2σ cutoff on m D , where σ stands for the m D resolution) is 97.7% and 99.3%, with about 487000 and 69000 candidates, The curves superimposed represent the result from the mD fit (solid blue lines) and the linear background contribution (dotted red lines).
events after all selection criteria, in the D 0 mass signal signal region. The contours (solid red lines) represent the kinematical limits of the

B. Dalitz plot analysis
Three-body charm decays are expected to proceed through intermediate quasi-two body modes [26] and this is the observed pattern. We therefore use, as a baseline model to describe A D (m 2 ∓ , m 2 ± ), an isobar approach consisting of a coherent sum of two-body amplitudes (sub-script r) and a "non-resonant" (subscript NR) contribution [27], where we have introduced the notation m ≡ (m 2 − , m 2 + ). The parameters a r (a NR ) and φ r (φ NR ) are the magnitude and phase of the amplitude for component r (NR).
as a function of position in the Dalitz plane. Here, F D (F r ) is the Blatt-Weisskopf centrifugal barrier factor for the D (resonance) decay vertex [28] with radius R = 1.5 GeV −1h c ≡ 0.3 fm, T r is the resonance propagator, and W r describes the angular distribution in the decay. For T r we use a relativistic Breit-Wigner (BW) parameterization with mass-dependent width [27], except for r = ρ(770) 0 and ρ(1450) 0 resonances where we use the Gounaris-Sakurai functional form [29]. The angular dependence W r is described using either Zemach tensors [30,31] where transversality is enforced or the helicity formalism [32][33][34] when we allow for a longitudinal component in the resonance propagator (see Ref. [27] for a comprehensive summary). Mass and width values are taken from [21], unless otherwise specified.
The complex ππ S-wave dynamics in the D 0 → K 0 S π + π − reaction [35], with the presence of several broad and overlapping scalar resonances, is more adequately described through the use of a K-matrix formalism [36] with the P-vector approximation [37]. This approach offers a direct way of imposing the unitarity constraint of the scattering matrix, not guaranteed in the case of the isobar model. The Dalitz plot amplitude A D (m) given by Eq. (6) is then modified as a r e iφr A r (m) + a NR e iφNR , (7) where F 1 (s) is the contribution of ππ S-wave states written in terms of the K-matrix formalism, Here, s = m 2 0 is the squared invariant mass of the π + π − system, I is the identity matrix, K is the matrix describing the S-wave scattering process, ρ is the phase-space matrix, and P is the initial production vector (P-vector). The index u (and similarly v) represents the u th channel (1 = ππ, 2 = KK, 3 = ππππ, 4 = ηη, 5 = ηη ′ ). In this framework, the production process can be viewed as the initial preparation of several states, which are then propagated by the [I − iK(s)ρ(s)] −1 term into the final one. The propagator can be described using scattering data, provided that the two-body system in the final state is isolated and does not interact with the rest of the final state in the production process. The P-vector has to be determined from the data themselves since it depends on the production mechanism. Only the F 1 amplitude appears in Eq. (7) since we are describing the ππ channel. See Sec. III C for more details.
The decay amplitude A D (m) is then determined from a maximum likelihood fit to the (9) where f sig represents the fraction of signal obtained from the fit to the mass spectrum, and D sig(bkg),∓ (m j ) is the signal (background) Dalitz plot PDF for event j, satisfying the condition D sig(bkg),∓ (m)dm = 1.
Here ǫ(m) represents the efficiency variations on the Dalitz plot, evaluated using high statistics signal MC samples. These are generated according to a uniform distribution and parameterized using third-order polynomial functions in two dimensions, symmetric for D 0 → K 0 S π + π − and asymmetric for D 0 → K 0 S K + K − to account for possible charge asymmetries in the K − and K + detection efficiencies. The Dalitz plot distributions for the background, D bkg,∓ (m), are determined using D 0 mass sideband data.
For each contribution r we evaluate the fit fraction as the normalized integral of a 2 r |A r (m)| 2 over the Dalitz plane [27], The sum of fit fractions does not necessarily add up to unity because of interference effects among the amplitudes.
The P-and D-waves of the D 0 → K 0 S π + π − decay amplitude are described using a total of 6 resonances leading to 8 two-body decay amplitudes: the Cabibbo allowed 1430) + , and the CP eigenstates ρ(770) 0 , ω(782), and f 2 (1270). Since the Kπ P-wave is largely dominated by the K * (892) ∓ , the mass and width of this resonance are simultaneously determined from our fit to the tagged D 0 sample, M K * (892) ∓ = 893.61 ± 0.08 MeV/c 2 and Γ K * (892) ∓ = 46.34 ± 0.16 MeV/c 2 (errors are statistical only). The mass and width values of the K * (1680) − are taken from [38], where the interference between the Kπ S-and P-waves is properly accounted for.
We adopt the same parameterizations for K, ρ, and P in Eq. (8) as in Refs. [27,39,40]. For the K matrix we have (11) where g α u is the coupling constant of the K-matrix pole m α to the u th channel. The parameters f scatt uv and s scatt 0 describe the slowly-varying part of the K-matrix. The factor suppresses the false kinematical singularity at s = 0 in the physical region near the ππ threshold (the Adler zero [41]). The parameter values used in this analysis are listed in Table I, and are obtained from a global analysis of the available ππ scattering data from threshold up to 1900 MeV/c 2 [39]. The parameters f scatt uv , for u = 1, are all set to zero since they are not related to the ππ scattering process. Similarly, for the P vector we have Note that the P-vector has the same poles as the Kmatrix, otherwise the F 1 vector would vanish (diverge) at the K-matrix (P-vector) poles. The parameters β α , f prod 1v and s prod 0 of the initial P-vector are obtained from our fit to the tagged D 0 → K 0 S π + π − data sample.  For the Kπ S-wave contribution to Eq. (7) we use a parameterization extracted from scattering data [38] which consists of a K * 0 (1430) − or K * 0 (1430) + BW (for CA or DCS contribution, respectively) together with an effective range non-resonant component with a phase shift, with The parameters a and r play the role of a scattering length and effective interaction length, respectively, F (φ F ) and R (φ R ) are the amplitudes (phases) for the nonresonant and resonant terms, and q is the momentum of the spectator particle in the Kπ system rest frame. Note that the phases δ F and δ R depend on m 2 Kπ . M and Γ(m 2 Kπ ) are the mass and running width of the resonant term. This parameterization corresponds to a Kmatrix approach describing a rapid phase shift coming from the resonant term and a slow rising phase shift governed by the non-resonant term, with relative strengths R and F [42]. The parameters M , Γ, F , φ F , R, φ R , a and r are determined from our fit to the tagged D 0 sample, along with the other parameters of the model. Other recent experimental efforts to improve the description of the Kπ S-wave using K-matrix and model independent parameterizations from high-statistics samples of D + → K − π + π + decays are described in Ref. [43]. Table II summarizes the values obtained for all free parameters of the D 0 → K 0 S π + π − Dalitz model: CA, DCS, and CP eigenstates complex amplitudes a r e iφr , π + π − Swave P-vector parameters, and Kπ S-wave parameters, along with the fit fractions. The non-resonant term of Eq. (7) has not been included since the ππ and Kπ Swave parameterizations naturally account for their respective non-resonant contributions. The fifth P-vector channel and pole have also been excluded since the ηη ′ threshold and the pole mass m 5 are both far beyond our ππ kinematic range, and thus there is little sensitivity to the associated parameters, f prod 15 and β 5 , respectively. The amplitudes are measured with respect to D 0 → K 0 S ρ(770) 0 which gives the second largest contribution. We report statistical errors only for the amplitudes, but for the fit fractions we also include systematic uncertainties (see Sec. III E), which largely dominate. The Kπ and ππ P-waves dominate the decay, but significant contributions from the corresponding Swaves are also observed (above 6 and 4 standard deviations, respectively). We obtain a sum of fit fractions of (103.6 ± 5.2)%, and the goodness of fit is estimated through a two-dimensional χ 2 test performed binning the Dalitz plot into square regions of size 0.015 GeV 2 /c 4 , yielding a reduced χ 2 of 1.11 (including statistical errors only) for 19274 degrees of freedom. The variation of the contribution to the χ 2 as a function of the Dalitz plot position is approximately uniform. Figure 6(a,b,c) shows the Dalitz fit projections overlaid with the data distributions. The Dalitz plot distributions are well reproduced, with some small discrepancies in low and high mass regions of the m 2 0 projection, and in the ρ(770) 0 − ω(782) interference region. As a cross-check, we alternatively parameterize the ππ and Kπ S-waves using the isobar approximation with the following BW amplitudes (plus the non-resonant contribution): the CA K * 0 (1430) − , the DCS K * 0 (1430) + , and the CP eigenstates f 0 (980), f 0 (1370), σ and an ad hoc σ ′ . This model is very similar to that used in our previous measurement of γ [18], except that the K * (1410) − and ρ(1450) 0 resonances have been removed because of their negligible fit fractions. Masses and widths of the σ and σ ′ scalars are obtained from the fit, M σ = 528 ± 5, Γ σ = 512 ± 9, M σ ′ = 1033 ± 4, and Γ σ ′ = 99 ± 6, given in MeV/c 2 . Mass and width values for the K * 0 (1430) ∓ , f 0 (980), and f 0 (1370) are taken from [44,45]. We obtain a sum of fit fractions of 122.5%, and a reduced χ 2 of 1.20 (with statistical errors only) for 19274 degrees of freedom, which strongly disfavors the isobar approach in comparison to the K-matrix formalism.
The curves are the reference model fit projections.
The description of the D 0 → K 0 S K + K − decay amplitude uses Eq. (6) and consists of five distinct resonances leading to 8 two-body decays: K 0 S a 0 (980) 0 , K 0 S φ(1020), K − a 0 (980) + , K 0 S f 0 (1370), K + a 0 (980) − , K 0 S f 2 (1270) 0 , K 0 S a 0 (1450) 0 , and K − a 0 (1450) + . This isobar model is essentially identical to that used in our previous analysis of the same reaction [46], but for the addition of the a 0 (1450) scalar, whose contribution is strongly supported by the much larger data sample, as well as of a D-wave contribution parameterized with the f 2 (1270) tensor. Attempts to improve the model quality by adding other contributions (including the non-resonant term) did not give better results.
The φ(1020) resonance is described using a relativistic BW, with mass and width left free in our fit to the D 0 tagged sample in order to account for mass resolution effects. The a 0 (980) resonance has a mass very close to the KK threshold and decays mostly to ηπ. Therefore it is described using a coupled channel BW [27,46], where the mass pole and coupling constant to ηπ are taken from [47], while the coupling constant to KK, g KK , is determined from our fit. Table III summarizes the values obtained for all free parameters of the D 0 → K 0 S K + K − Dalitz model, the complex amplitudes a r e iφr , the mass and width of the φ(1020) and the coupling constant g KK , together with the fit fractions. The value of g KK is consistent with our previous result [46], and differs significantly from the measurement reported in [47]. All amplitudes are measured with respect to D 0 → K 0 S a 0 (980) 0 , which gives the largest contribution. The sum of fit fractions is 152.3%, and the reduced χ 2 is 1.09 (with statistical errors only) for 6856 degrees of freedom, estimated from a binning of the Dalitz plot into square regions of size 0.045 GeV 2 /c 4 . The variation of the contribution to the χ 2 as a function of the Dalitz plot position is approximately uniform in the regions where most of the decay dynamics occurs. Figure 6(d,e,f) shows the fit projections overlaid with the data distributions. The Dalitz plot distributions are well reproduced, with some small discrepancies at the peaks of the m 2 − and m 2 + projections.

E. Systematic uncertainties
Systematic uncertainties on A D (m) are evaluated by repeating the fit to the tagged D 0 samples with alternative assumptions to those adopted in the reference D 0 → K 0 S π + π − and D 0 → K 0 S K + K − amplitude analyses. These uncertainties can then be directly propagated to the measurement of the CP parameters, as discussed in Sec. IV B, and the total systematic error can be obtained III: CP eigenstates, CA, and DCS complex amplitudes are iφr and fit fractions, obtained from the fit of the D 0 → K 0 S K + K − Dalitz plot distribution from D * + → D 0 π + . The mass and width of the φ(1020), and the g KK coupling constant are simultaneously determined in the fit, yielding M φ(1020) = 1.01943 ± 0.00002 GeV/c 2 , Γ φ(1020) = 4.59319 ± 0.00004 MeV/c 2 , and g KK = 0.550±0.010 GeV/c 2 . Errors for amplitudes are statistical only. Uncertainties (largely dominated by systematic contributions) are not estimated for the fit fractions. from the sum square of the individual contributions. In this paper we have also propagated these systematic uncertainties to the measurement of the D 0 → K 0 S π + π − fit fractions, as reported in Table II. In general, each of the considered alternative models has a reduced χ 2 poorer than that of the reference. Therefore, our systematic uncertainties do not include potential contributions due to the residual poor quality of the reference model fit, as reported in Secs. III C and III D.

Model contributions
Dalitz model systematic uncertainties on A D (m) are related to the model dependence of the strong charm decay phase as a function of the Dalitz plot position when it is determined from the Dalitz plot density, which only depends on decay rates.
We use alternative models where the BW parameters are varied according to their uncertainties or changed by values measured by other experiments. This is the case of the f 0 (1370), where the reference values [44] are replaced by alternative measurements [48], and the K * (1680) − , where the reference parameters [38] are replaced by those from [21]. We also build models using alternative parameterizations, as in the case of the ρ(770) 0 where the reference Gounaris-Sakurai form is replaced by the standard relativistic BW.
To estimate the ππ S-wave systematic error we replace the reference K-matrix solution (Table I) by all alternative solutions analyzed in Ref. [39]. Analogously, the uncertainty on the parameterization of the Kπ S-wave is estimated using a standard relativistic BW describing the K * 0 (1430) ∓ with parameters taken either from [44] or simultaneously determined from our fit to the D 0 sample. Additionally, the isobar model is used as a cross-check of the combined ππ and Kπ S-wave effect.
Uncertainties due to our choice of the angular dependence are estimated by replacing the reference Zemach tensors by the helicity formalism. The effect is negligible for S-waves, very small for P-waves, but larger for D-waves [31]. Other alternative models are built by changing the Blatt-Weisskopf radius between 0 and 3 GeV −1h c, and removing and adding resonances with small or negligible fit fractions. For D 0 → K 0 S π + π − , we added ρ(1450) 0 and CA K * (1410) − . Similarly, for D 0 → K 0 S K + K − , we removed all the a 0 (1450) charged states and the f 2 (1270), and added the f 0 (980) and the charged DCS a 0 (1450). The f 0 (980) resonance is described using a coupled channel BW with parameters taken from a variety of experiments [45,48,49].

Experimental contributions
Experimental systematic errors come from uncertainties in the knowledge of variations of the reconstruction efficiency on the Dalitz plot, background Dalitz plot shapes, mass resolution, mistag rate, and binning.
The uncertainty from the efficiency variations on the Dalitz plot ǫ(m) has been evaluated assuming the efficiency to be flat. Tracking efficiency studies in data and MC show that this method gives a conservative estimate of the imperfections of the detector simulation, which appear mainly at the boundaries of the phase space because of the presence of very low momentum tracks.
Systematic errors related to the background Dalitz plot profile D bkg,∓ (m) are determined assuming a flat shape, which gives the largest effect among other alternative profiles obtained using either the m D sideband from continuum MC, or the m D signal region from continuum MC after removal of true D 0 mesons.
All the resonances, except for the ω(782) and φ(1020), have intrinsic width significantly larger than possible bias on invariant mass measurement and resolution. We estimate the systematic uncertainty associated with ω(782) by repeating the model fit using an overall width resulting from adding in quadrature its natural width and the mass resolution in the 782 MeV/c 2 ππ mass region. No systematic error is assigned for the φ(1020), since the reference model has been extracted with its mass and width as free parameters.
The uncertainty due to a wrong identification of the flavor of the D 0 (D 0 ) meson from the D * + → D 0 π + (D * − → D 0 π − ) decay, due to the association of the D meson with a random soft pion of incorrect charge has been evaluated taking into account explicitly the rate of mistags observed in the MC, at 0.7% level.
Effects from limited numerical precision in the computation of normalization integrals and binning in the Dalitz plane have been evaluated using coarser and thinner bins.

IV. DALITZ PLOT ANALYSIS OF
Once the decay amplitudes A D (m) for D 0 → K 0 S π + π − and D 0 → K 0 S K + K − are known, they are fed into Γ Here, D c,∓ (m j ) is the Dalitz plot PDF for event j satisfying the normalization condition D c,∓ (m)dm = 1, and A c accounts for any asymmetry in the absolute number of B − and B + candidates (charge asymmetry) for component c.
For B ∓ →D 0 K ∓ signal, D sig,∓ (m) = Γ ∓ (m)ǫ(m), where the efficiency map in the Dalitz plot ǫ(m) is determined as for D * + → D 0 π + events (Sec. III B). We replace r 2 B in Eq. (2) by r 2 B ∓ = x 2 ∓ + y 2 ∓ . The physical condition r B − = r B + is recovered in the statistical procedure to extract γ from x ∓ , y ∓ , as discussed in Sec. V. The same procedure is applied analogously to the other signal samples.
We consider the same background components as in the selection fit, with some important modifications. First, events falling into the continuum and BB background components are divided into events with a real or a fake (combinatorial)D 0 meson. Dalitz plot shapes for fakeD 0 mesons from continuum are extracted as described in Sec. III B, using events in the continuum enriched region (m ES and m D sideband regions), while those from BB are determined from MC events. Events containing a realD 0 are further divided into "right-sign" and "wrong-sign" flavor categories depending on whether they are combined with a negative or positive kaon (or K * ). We pay special attention to this charge-flavor correlation in the background since it can mimic either the b → c or the b → u signal component. Second, we have included a background contribution due to signal events where the kaon (or K * ) comes from the other B decay; this amounts to 9% of the B − →D 0 K * − signal, but is negligible for B − →D ( * )0 K − .
and (h) B + →D 0 K * + , for the ∆E signal region. The requirements mES > 5.272 GeV/c 2 and F > −0.1 have been applied to reduce the background contamination, mainly from continuum events. The contours (solid red lines) represent the kinematical limits of theD 0 → K 0 S π + π − decay.
The m ES , ∆E, and F PDF parameters in the CP fit are the same as those used in or obtained from the selection fit, except for the m ES peaking fractions for D 0 → K 0 S π + π − channels, which are allowed to vary since their values depend on the ∆E region used for the fit. Other parameters simultaneously determined ∓ , x s∓ , and y s∓ , are: signal and background yields, signal charge asymmetries, and fractions of true D 0 mesons for all decay modes and right-sign fractions forD 0 → K 0 S π + π − channels in continuum background. Right-sign fractions for the modes withD 0 → K 0 S K + K − are fixed from MC simulation due to lack of statistics and the limited discriminating power between D 0 and D 0 Dalitz plot distributions. Similarly, fractions of truẽ D 0 mesons and charge-flavor correlation for the BB component are determined using MC events, because of the lack of BB background statistics.
A variety of studies using data, parameterized fast Monte Carlo, and full GEANT4-simulated samples have been performed to test the consistency of the results and to verify the analysis chain and fitting procedure, as described below.
The CP fit to the B − →D ( * )0 K − samples has been performed separately forD 0 → K 0 S π + π − andD 0 → K 0 S K + K − samples. Figure 10 shows the resulting oneand two-standard deviation regions in the (x ( * ) ∓ , y ( * ) ∓ ) planes. We find statistically consistent results between the different subsets. The same fitting procedure has been applied to the B − → D ( * )0 π − control samples. In this case we expect r ∓,π ) contours for B − and B + decays should be close to the origin up to ∼ 0.01. Deviations from this pattern could be an indication that the Dalitz plot distributions are not well described by the models. Figure 11 shows the resulting one-and two-standard deviation regions for (x ( * ) ∓,π , y ( * ) ∓π ), consistent with the expected values. Moreover, we find statistically consistent results between the D 0 → K 0 S π + π − and D 0 → K 0 S K + K − samples. An additional test of the fitting procedure is performed with parameterized MC simulations consisting of about 500 experiments generated with a sample size and composition corresponding to that of the data. The CP parameters are generated with values close to those found in the data and the reference CP fit is performed on each of these experiments. The r.m.s. of the residual distributions for all the CP parameters (where the residual is ∓ , xs∓, and ys∓, as obtained from the CP fit. The first error is statistical, the second is experimental systematic uncertainty and the third is the systematic uncertainty associated with the Dalitz models. 090 ± 0.043 ± 0.015 ± 0.011 −0.111 ± 0.069 ± 0.014 ± 0.004 0.115 ± 0.138 ± 0.039 ± 0.014 y− , y * − , ys− 0.053 ± 0.056 ± 0.007 ± 0.015 −0.051 ± 0.080 ± 0.009 ± 0.010 0.226 ± 0.142 ± 0.058 ± 0.011 x+ , x * + , xs+ −0.067 ± 0.043 ± 0.014 ± 0.011 0.137 ± 0.068 ± 0.014 ± 0.005 −0.113 ± 0.107 ± 0.028 ± 0.018 y+ , y * + , ys+ −0.015 ± 0.055 ± 0.006 ± 0.008 0.080 ± 0.102 ± 0.010 ± 0.012 0.125 ± 0.139 ± 0.051 ± 0.010 defined as the difference between the fitted and generated values) is found to be consistent with the mean (Gaussian) statistical errors reported by the fits. The mean values of the residual distributions are consistent with zero.
Only for x s∓ and y s∓ we observe small biases (at 10% level of the statistical uncertainty), as a consequence of the non-Gaussian behavior of samples with small statistics. This small deviation from Gaussian behavior is also observed in the data, as shown in Fig. 9(c). The statistical errors on the CP parameters and the calculated correlation coefficients among them extracted from the fit are consistent with the range of values obtained from these experiments. We also observe that the fit errors are independent of the truth values.
Finally, samples of signal and background GEANT4simulated MC events with a full detector simulation are used to validate the measurement. We performed fits to signal samples, using the true and reconstructed B meson charge andD 0 Dalitz plot distributions, obtaining in all cases results consistent with those generated.

Dalitz model contributions
Dalitz model uncertainties are evaluated by repeating the fit to the tagged D 0 samples with alternative assumptions to those adopted in the reference D 0 → K 0 S π + π − and D 0 → K 0 S K + K − amplitude analyses (Sec. III E), and then are propagated to the CP parameters. To propagate each systematic uncertainty on A D (m) to the CP parameters we have generated samples of B − → D ( * )0 K − and B − →D 0 K * − signal events that are one hundred times larger than each measured signal yield in data. These virtually infinite samples reduce to a negligible level statistical differences between the models. The D 0 Dalitz plot distributions are generated according to the reference models and to CP parameters consistent with the values found in data. The CP parameters are then extracted by fitting the generated Dalitz plot distributions using the reference or one of the alternative models. The difference is taken as the systematic uncertainty associated with each alternative model, and the sign of the variation is used to estimate whether the different contributions are positively or negatively correlated (Appendix A). When two alternative models are built from an up and down variation of the same param- eter, we take the maximum variation as the systematic error. Assuming the contributions are uncorrelated, we sum in quadrature to obtain the total systematic uncertainty.
The statistical errors in the Dalitz model parameters obtained from the tagged D 0 samples have been propagated to the CP parameters by repeating the CP fit with those parameters randomized according to their covariance matrix. Table V summarizes the main contributions from all the alternative models considered and discussed in Sec. III E. Contributions from other models are found to be negligible.
We have also evaluated the effect on the measured CP parameters when we parameterize the ππ and Kπ Swaves in D 0 → K 0 S π + π − using the isobar model instead of the K-matrix model (plus the non-resonant contribution), as described in Sec. III C. The variations are found to be smaller than the sum of the ππ and Kπ S-wave systematic uncertainties, and are used as a cross-check of the procedure adopted for assigning this contribution to the total Dalitz model error.

Experimental contributions
Experimental systematic uncertainties arise from several sources and their main contributions are summarized in Table VI. They are small compared to the statistical precision, and their sum is similar to the Dalitz model uncertainty. Other sources of experimental systematic uncertainty, e.g. the assumption of perfect mass resolution for the Dalitz plot variables m, are found to be negligible.
Statistical uncertainties due to the m ES , ∆E, and F PDF parameters for signal and background extracted from the selection fit (fixed in the reference CP fit) are estimated by repeating the CP fit with PDF parameters randomized according to their covariance matrix. Possible bias due to differences in the m ES and F shapes for continuum and BB background events between the ∆E selection and signal regions are evaluated applying the selection fit in the ∆E signal region. Other PDF parameters, such as the m ES end-point, BB ∆E peaking fractions, m ES BB peaking fractions forD 0 → K 0 S K + K − channels, and PEP-II boost are varied by one standard deviation. We account for m ES and ∆E differences in BB background for true and fake D mesons, while for continuum events we do not observe differences. We also find the effect of the small correlation between m ES , ∆E, and F variables negligible.
The uncertainties related to the knowledge of theD 0 fractions for the small BB background are estimated from the maximum variations of the CP parameters when the fractions are varied one σ up and down from their MC estimates, or replaced by the values found for the continuum background, or assumed to be zero. Similarly, the uncertainties due to our knowledge of the right-sign fractions forD 0 → K 0 S K + K − continuum events and BB events are evaluated from the maximum variations of the CP parameters after varying these fractions according to their MC values or assuming that theD 0 is randomly  The effect due to reconstruction efficiency variations of the signal across the Dalitz plane, ǫ(m), has been evaluated by varying randomly the coefficients of the polynomial parameterization according to their covariance matrix, including the statistical errors due to the limited MC statistics as well as systematic uncertainties arising from the imperfections of the detector simulation, as discussed in Sec. III E 2.
The uncertainty associated with the knowledge of the Dalitz plot distributions of continuum background events is taken to be the difference in the CP parameters using background Dalitz plot shapes from sideband data instead of signal region backgrounds from MC. We also account for statistical uncertainties adding in quadrature the r.m.s. of the distributions of CP parameters when the two sets of profile distributions are randomized. Uncertainties due to the Dalitz plot shapes of combinatorial D mesons in BB background are conservatively estimated from the variation of CP parameters when the reference shapes are replaced by a flat profile.
The The B − →D 0 K * − sample has two additional sources of uncertainty. The first one comes from signal events where the prompt K * − is replaced by a combinatorial K * − (about 9% of the signal), with either the same or opposite charge. This systematic uncertainty, evaluated by changing by ±10% the fraction of these events and neglecting the charge-flavor correlation, has been found to be negligible.
The second additional uncertainty is due to our knowledge of the parameter κ, as defined in Eq. (3), which accounts for the interference between B − →D 0 K * − and other B − →D 0 K 0 S π − (higher K * resonances plus nonresonant) decays. Since this parameter cannot be extracted from the CP fit and no experimental data analysis is available on the B − →D 0 K 0 S π − decay, we study a B − Dalitz (isobar) model including K * (892) − , K * 0 (1410) − , K * 2 (1430) − , D * (2010) − , D * 2 (2460) − and non-resonant terms, and randomly varying phases in the range [0, 2π] and magnitudes [50]. The magnitude of the contribution from b → c transitions relative to b → u was fixed to be around 3, while the magnitude of the non-resonant contribution was varied between 0 and 1. Since our model has a large uncertainty we made several alternative models adding/removing resonances and changing ranges for b → u amplitudes, keeping the K * pollution (defined as the non-K * fit fraction) below 5-10%, since from earlier studies with very similar selection criteria we estimate that, neglecting higher resonances, the nonresonant K * decays contribute about 5% of the signal events [51]. Evaluating κ from Eq. (3) for the region within 55 MeV/c 2 of the K * mass and | cos θ H | ≥ 0.35 we find quite narrow distributions, centered around 0.9 and with r.m.s. not larger than 0.1, in agreement with previous studies [52]. For this reason we have fixed the value of κ to 0.9 in the reference CP fit, and varied it between 0.8 and 1.

V. INTERPRETATION OF RESULTS
A frequentist procedure [21] has been adopted to transform the measurement of the CP parameters z ∓ ≡ (x ∓ , y ∓ , x * ∓ , y * ∓ , x s∓ , y s∓ ) into the measurement of the physically relevant quantities p ≡ (γ, r B , r * B , κr s , δ B , δ * B , δ s ). Using a large number of pseudo-experiments with probability density functions and parameters as obtained from the fit to the data but with different values of the CP parameters, we construct a multivariate Gaussian likelihood function L(p|z, C) relating the experimentally measured observables z ≡ {z + , z − } (reported in Table IV) and their 12 × 12 statistical and systematic covariance matrices C with the corresponding true values calculated using their definition in terms of the quantities p. The matrices C are constructed from the uncertainties summarized in Table IV and the statistical and systematic correlation coefficients given in Sec. IV A and Appendix A, respectively. For a single B decay channel the procedure is identical to that outlined here but with a reduced space of measured and truth parameters. For example, for B − →D 0 K − , z ± ≡ (x ± , y ± ), C is the corresponding 4 × 4 covariance matrix, and p ≡ (γ, r B , δ B ).
We evaluate the confidence level (CL) as a function of the true value for a given parameter µ from p ≡ {µ, q}, minimizing the function χ 2 (p|z, C) ≡ −2 ln L(p|z, C) with respect to the parameters q. For each given value µ 0 of µ, between its minimum and maximum value, the fit provides a minimum chi-square χ 2 (µ 0 , q 0 ), where q 0 are the best parameters for the given µ 0 and the actual z measurements with covariance matrix C. Then we take the values p best ≡ {µ best , q best } for which χ 2 (µ best , q best ) is minimum and compute the χ 2 -difference ∆χ 2 (µ 0 ) = χ 2 (µ 0 , q 0 ) − χ 2 (µ best , q best ) In a purely Gaussian situation for the truth parameters p, the CL can be obtained by computing the probability that this value is exceeded for a χ 2 -distribution with one degree of freedom, CL = 1 − α = F (∆χ 2 (µ 0 ); ν = 1), where F (∆χ 2 (µ 0 ); ν = 1) is the corresponding cumulative distribution function. In a non-Gaussian situation one has to consider ∆χ 2 (µ 0 ) as a test statistic, and has to rely on a Monte Carlo simulation to obtain its expected distribution. This Monte Carlo simulation is built by generating a large number of samples with truth values p 0 ≡ {µ 0 , q 0 } as determined from the actual data analysis, and then counting the number of experiments for which ∆χ ′2 (µ 0 ) < ∆χ 2 (µ 0 ), where ∆χ ′2 (µ 0 ) = χ ′2 (µ 0 , q ′ 0 ) − χ ′2 (µ ′ best , q ′ best ) is determined by letting the q parameters free to vary for each of the generated (primed) samples. The one-(two-) standard deviation region of the CP parameters is defined as the set of µ 0 values for which α is greater than 31.7% (4.6%).
This technique to obtain the physical parameters takes into account unphysical regions of the parameter space [53], which may arise since in the z measurements we allow B − and B + events to have different r B − and r B + values, while the space of true values is built using a common r B parameter. Moreover, this approach provides 1-dimensional intervals that include the true value as implied by the confidence level, while in previous measurements [18,54] the 1-dimensional intervals were determined from projections of the multidimensional confidence regions onto each of the parameters. Figure 12 shows α = 1 − CL as a function of the parameter γ, for each of the three B decay channels separately and their combination. As expected from Eq. (2), the method has a two-fold ambiguity in the weak and strong phases, (γ; δ cays, respectively. For the combined analysis of the three charged B → DK decay modes we obtain CL = 0.997, corresponding to 3.0 standard deviations.

VI. CONCLUSION
In summary, using 383 million BB decays recorded by the BABAR detector, we have performed a new measurement of the direct CP -violating parameters (x ( * ) ∓ , y ( * ) ∓ ) and (x s∓ , y s∓ ) in B − →D ( * )0 K − and B − →D 0 K * − decays, respectively, using a Dalitz plot analysis ofD 0 → K 0 S π + π − andD 0 → K 0 S K + K − . Compared to our previous analysis based on 227 million BB decays [18], this measurement takes advantage of significant improvements in reconstruction efficiencies, treatment of e + e − → qq, q = u, d, s, c background, and Dalitz models, along with the use, for the first time, ofD 0 → K 0 S K + K − decays. These upgrades result in reduced experimental and Dalitz model systematic uncertainties, and statistical uncertainties improved beyond the increase in data sample size. The results, summarized in Table IV are consistent with, and improve significantly, the previous measurements from BABAR and Belle [18,54].
A significant reduction in Dalitz model systematic uncertainties has been achieved through the detailed study of high-statistics samples of e + e − → cc → D * + → D 0 π + decays, where the D 0 is reconstructed in the K 0 S π + π − and K 0 S K + K − final states. We have adopted a K-matrix formalism to describe the complex ππ and Kπ S-wave dynamics in D 0 → K 0 S π + π − . For this decay, the fit fractions measured and reported in Table II, show that the Kπ and ππ P-waves dominate, but for the first time significant contributions from the corresponding S-waves are observed (above 6 and 4 standard deviations, respectively).
Using a frequentist analysis we interpret the (x ( * ) ∓ , y ∓ ) and (x s∓ , y s∓ ) experimental results in terms of the weak phase γ, the amplitude ratios r B , r * B , and r s , and the strong phases δ B , δ * B , and δ s . We obtain γ = (76 ± 22 ± 5±5) • (mod 180 • ), where the first error is statistical, the second is the experimental systematic uncertainty and the third reflects the uncertainty on the D decay Dalitz models (parabolic errors). The corresponding two standard deviation region is 29 • < γ < 122 • . The combined significance of direct CP violation (i.e. γ = 0) is 99.7%, corresponding to 3.0 standard deviations. This direct determination of γ supersedes and significantly improves our previous constraint [18], and is consistent with that reported by the Belle Collaboration [54]. The latter has a slightly better precision in spite of a larger uncertainty on the measured CP parameters because the error on γ scales roughly as 1/r ( * ) B (1/r s ) and our tighter r ( * ) B , κr s constraints favor smaller values.

VII. ACKNOWLEDGMENTS
We are grateful for the extraordinary contributions of our PEP-II colleagues in achieving the excellent luminosity and machine conditions that have made this work possible. The success of this project also relies critically on the expertise and dedication of the computing organizations that support BABAR. The collaborating institutions wish to thank SLAC for its support and the kind hospitality extended to them. This work is supported by the US Department of Energy and National Science Foundation, the Natural Sciences and Engineering Research Council (Canada), the Commissariatà l'Energie Atomique and Institut National de Physique Nucléaire et de Physique des Particules (France), the Bundesministerium für Bildung und Forschung and Deutsche Forschungsgemeinschaft (Germany), the Istituto Nazionale di Fisica Nucleare (Italy), the Foundation for Fundamental Research on Matter (The Netherlands), the Research Council of Norway, the Ministry of Education and Science of the Russian Federation, Ministerio de Educación y Ciencia (Spain), and the Science and Technology Facilities Council (United Kingdom). Individuals have received support from the Marie-Curie IEF program (European Union) and the A. P. Sloan Foundation.