Measurement of CP observables in B ±-- > DCPK ± decays and constraints on the CKM angle gamma

Using the entire sample of 467 x 10(6) Y(4S) -> B (B) over bar decays collected with the BABAR detector at the PEP-II asymmetric-energy B factory at the SLAC National Accelerator Laboratory, we perform an analysis of B-+/- -> DK +/- decays, using decay modes in which the neutral D meson decays to either CP-eigenstates or non-CP-eigenstates. We measure the partial decay rate charge asymmetries for CP-even and CP-odd D final states to be A(CP+) = 0.25 +/- 0.06 +/- 0.02 and A(CP-) = 0.09 +/- 0.07 +/- 0.02, respectively, where the first error is the statistical and the second is the systematic uncertainty. The parameter A(CP+) is different from zero with a significance of 3.6 standard deviations, constituting evidence for direct CP violation. We also measure the ratios of the charged-averaged B partial decay rates in CP and non-CP decays, RCP+ 1.18 +/- 0.09 +/- 0.05 and RCP- = 1.07 +/- 0.08 +/- 0.04. We infer frequentist confidence intervals for the angle gamma of the unitarity triangle, for the strong phase difference delta(B), and for the amplitude ratio r(B), which are related to the B- -> DK- decay amplitude by r(B)e(i(delta B-gamma)) = A(B- -> (D) over bar K-0(-)) = A(B- -> (D) over bar K-0(-))/A(B- -> (DK-)-K-0). Including statistical and systematic uncertainties, we obtain 0: 24 < rB < 0: 45 ( 0: 06 < rB < 0: 51) and, modulo 180 degrees, 11.3 degrees < gamma < 22.7 degrees or 80.8 degrees < gamma < 99.2 degrees or 157.3 degrees < gamma < 168.7 degrees (7.0 degrees < gamma < 173.0 degrees) at the 68% ( 95%) confidence level.


I. INTRODUCTION
In the standard model (SM) of fundamental particles, CP violation in weak interactions is allowed by a single, irreducible phase in the 3 Â 3 Cabibbo-Kobayashi-Maskawa (CKM) quark flavor-mixing matrix [1,2]. The unitarity of the CKM matrix, V, implies a set of relations among its elements V ij , in particular, the condition V ud V Ã ub þ V cd V Ã cb þ V td V Ã tb ¼ 0, which can be depicted in the complex plane as a ''unitarity'' triangle, whose sides and angles are related to the magnitudes and phases of the six elements V id and V ib , where i ¼ u, c, t. Overconstraining the unitarity triangle by means of precise measurements of all its sides and angles allows tests of whether the CKM mechanism is the correct description of CP violation. Any inconsistencies among the various experimental constraints would reveal effects of physics beyond the standard model.
After a decade of successful operation and a total of about 1:3 Â 10 9 B " B pairs collected by the BABAR and Belle experiments, the three CKM angles have been measured with varied precision. The angle has been measured with the highest precision, to around 1 , using B 0 ! ðc " cÞK ðÃÞ0 decays. Using a variety of two-body B decays (B ! , , and a 1 ð1260Þ) the angle has been measured to a precision of around 4 . The angle has a relatively large uncertainty, around 14 , compared with and . The lack of precision in our knowledge of reflects the difficulty in measuring this angle. The uncertainties of the CKM angles quoted in this paragraph are taken from [3].
Several techniques for measuring in a theoretically clean way are based on B meson decays to open-charm final states, D ðÃÞ0 X s and " D ðÃÞ0 X s (X s ¼ K ðÃÞAE , K ðÃÞ0 ). In these decays, the interference between the b ! c " us and b ! u " cs tree amplitudes, when the D 0 and " D 0 decay to a common final state, leads to observables that depend on the relative weak phase . The size of the interference also depends on the magnitude of the ratio r B and the relative strong phase B of the two amplitudes, which cannot be precisely calculated from theory. They can be extracted directly from data by simultaneously reconstructing several related B ! DK decays. Many methods have been proposed to extract from B decays using D ðÃÞ K ðÃÞAE and D ðÃÞ K ðÃÞ0 final states (here and in the following D refers to any admixture of the neutral D 0 meson and its CP-conjugate " D 0 ). The three methods that have been used most productively to date are the ''GLW'' method [4,5], based on Cabibbo-suppressed D decays to CP-eigenstates, such as K þ K À or K 0 S 0 ; the ''ADS'' method [6,7], where the D is reconstructed in Cabibbofavored and doubly-Cabibbo-suppressed final states such as K AE Ç ; and the ''GGSZ'' method [8], which studies the Dalitz-plot distribution of the products of D decays to multibody self-conjugate final states, such as K 0 S þ À . A common problem with these methods is the small overall branching fraction of these decays ranging from 5 Â 10 À6 to 5 Â 10 À9 . Therefore a precise determination of requires a very large data sample. BABAR has published several related measurements: GLW analyses of B AE ! DK AE [9], D Ã K AE [10] and DK ÃAE [11] decays; ADS analyses of B AE ! D ðÃÞ K AE [12,13], DK ÃAE [11] and B 0 ! DK Ã0 [14]; and GGSZ analyses of B AE ! D ðÃÞ K ðÃÞAE [15,16] and B 0 ! DK Ã0 decays [17]. To date, the single most precise experimental determination of from BABAR is ¼ ð68 AE 14 AE 4 AE 3Þ and 39 < < 98 , obtained from the GGSZ analysis of B AE ! D ðÃÞ K ðÃÞAE decays [16]. In this measurement, the first error represents the statistical uncertainty, the second is the experimental systematic uncertainty, and the third reflects the uncertainty on the description of the D Dalitz-plot distributions.

II. GLW ANALYSIS OF B ! DK DECAYS
In this paper we present the update of the GLW analysis of B AE ! DK AE decays based on the full BABAR data set collected near the Çð4SÞ resonance. In addition to a 22% increase in statistics of the data sample, this study benefits from other significant improvements compared to our previous result [9]: (i) More refined charged track reconstruction and particle identification algorithms, with higher purity and efficiency, have been employed; (ii) The event shape variable F , used to discriminate the signal from the continuum e þ e À ! q " q background (described in detail in Sec. IV) has been removed from the selection criteria and has instead been included in the final fit to the selected B candidates. This allows us to increase the signal efficiency by 40% to 60%. At the same time it provides a larger sample of continuum background events, thus allowing for the determination of the background properties directly from data (see Sec. V); (iii) Better kaon/pion separation, which is needed to distinguish B AE ! DK AE candidates from the 12 times more abundant B AE ! D AE decays, is achieved through the use of a global likelihood based not only on the Cherenkov angle C reconstructed by the Cherenkov detector, but also on the specific energy loss dE=dx measured by the tracking devices. The inclusion of dE=dx in the likelihood increases the kaon identification efficiency and decreases the pion misidentification both at low momentum and outside of the geometrical acceptance of the Cherenkov detector (which is 10% lower than the acceptance of the tracking devices). In order to determine from B AE ! DK AE decays with the GLW method, we measure the two direct-CP-violating partial decay rate asymmetries,

A CPAE
ÀðB À ! D CPAE K À Þ À ÀðB þ ! D CPAE K þ Þ ÀðB À ! D CPAE K À Þ þ ÀðB þ ! D CPAE K þ Þ ; (1) and the two ratios of charge averaged partial rates using D decays to CP and flavor eigenstates, where D CPAE refer to the CP eigenstates of the D meson system. We then extract , together with the other two unknowns r B and B , by means of a frequentist procedure, which exploits the following relations [4,5], neglecting D 0 -" D 0 mixing [18]: Here, r B jAðB À ! " D 0 K À Þ=AðB À ! D 0 K À Þj is the magnitude of the ratio of the amplitudes for B À ! " D 0 K À and B À ! D 0 K À and B the difference of their strong phases. Taking into account the CKM factor (jV ub V cs = V cb V us j % 0:4) and color-suppression of the B À ! " D 0 K À amplitude, r B is expected to be around 0.1. The current world averages for the B AE ! DK AE GLW observables from the measurements in [9,19,20] are summarized in Table I. The world averages for the parameters r B and B are r B ¼ 0:104 þ0:015 À0:025 and B ¼ ð117 þ17 À24 Þ at 68% confidence level (CL) [3].
To reduce the systematic uncertainties from branching fractions and reconstruction efficiencies of different D channels appearing in the numerator and denominator of Eq. (2), we approximate R CPAE with the double ratios where and Equation (5) would be exact in the limit in which the Cabibbo-suppressed contributions to the B AE ! D AE amplitudes vanish, as well as terms proportional to r B r D % 5 Â 10 À3 , as we will discuss in Sec. VII. This approximation results in a systematic uncertainty on the final values of R CPAE . The paper is organized as follows. In Sec. III we describe the data sample used for these measurements and the main features of the BABAR detector and of the PEP-II storage rings. In Sec. IV we summarize the procedure adopted to select B AE ! Dh AE candidates and suppress the main backgrounds. In Sec. V we introduce the simultaneous extended maximum likelihood fit used to extract the observables R CPAE and A CPAE . In Sec. VI we explain how, by applying the same fit procedure to selected control samples, we estimate the irreducible background present in the final samples. A discussion of the sources of systematic uncertainties and the evaluation of the uncertainties is presented in Sec. VII. Section VIII lists the final results on the GLW observables R CPAE and A CPAE , including statistical and systematic uncertainties. It also contains a description of the statistical method used to construct frequentist confidence intervals for the parameters , B , and r B . Section IX gives a summary of our results.

III. DATA SAMPLE AND DETECTOR
The measurements presented in this paper use the entire B " B data sample collected with the BABAR detector at the PEP-II asymmetric-energy B factory at the SLAC National Accelerator Laboratory. The B " B pairs are produced from the decays of Çð4SÞ mesons that originate in collisions of 9.0 GeV electrons and 3.1 GeV positrons ( ffiffi ffi s p ¼ 10:58 GeV ¼ M Çð4SÞ c 2 Þ. In total, ð467 AE 5Þ Â 10 6 B " B pairs, approximately equally divided into B 0 " B 0 and B þ B À , have been collected in the years from 1999 until early 2008. The B meson pairs are produced almost at rest in the Çð4SÞ center-of-mass (CM) frame, but the asymmetric beam energies boost them in the laboratory frame by ðÞ CM % 0:56. The BABAR detector is described in detail elsewhere [22]. Primary and secondary vertex reconstruction and charged-particle tracking are provided by a five-layer double-sided silicon vertex tracker and a 40-layer drift chamber. Charged particle identification (PID) is provided by measurement of specific ionization energy loss in the tracking devices and of the Cherenkov radiation cone in a ring-imaging detector. Photons and electrons are identified by combining the information from the tracking devices and the energy deposits in the electromagnetic calorimeter, which consists of 6580 thallium-doped CsI crystals. These systems are located inside a 1.5 T solenoidal superconducting magnet. Finally, the flux return of the magnet is instrumented with resistive plate chambers and limited streamer tubes in order to discriminate muons from pions. We use the GEANT4 [23] software toolkit to simulate interactions of particles in the detector, taking into account the varying accelerator and detector conditions.

IV. EVENT SELECTION
We reconstruct B AE ! Dh AE decays, where the charged track h is either a kaon or a pion. Neutral D mesons are reconstructed in the CP-even eigenstates À þ and is assumed to be a pure CP ¼ þ1 eigenstate. The D CP daughters are reconstructed in the decay modes K 0 S ! þ À , ! K þ K À and ! ! À þ 0 . We optimize all our event selection requirements by maximizing the significance of the expected B AE ! DK AE signal yield, defined as N sig = ffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi N sig þ N bkg p , where N sig (N bkg ) is the expected signal (background) yield. The optimization is done for each D decay channel using simulated signal and background events, which are generated with the EVTGEN software package [24].
Neutral pions are reconstructed by combining pairs of photon candidates with energy deposits larger than 30 MeV that are not matched to charged tracks and whose energy deposition profile is consistent with that expected from a photon. The photon pair invariant mass is required to differ from the nominal 0 mass [25] by less than 2.5 times its resolution ( % 6 MeV=c 2 ) and the total 0 energy in the laboratory frame must be greater than 240 MeV for D ! K 0 S 0 and 210 MeV for ! ! þ À 0 . Neutral kaons are reconstructed from pairs of oppositely charged tracks with invariant mass within 2:5 ( % 2:1 MeV=c 2 ) of the nominal K 0 S mass [25]. The ratio between the K 0 S signed 3-dimensional flight length and its uncertainty, determined from the position of the K 0 S and the D decay vertices and the K 0 S momentum direction, must be greater than 1.9, 2.0, and 2.2 for D ! K 0 S 0 , D ! K 0 S , and D ! K 0 S !, respectively. The candidates are reconstructed from pairs of oppositely charged tracks passing kaon identification criteria with typical kaon selection efficiency of % 98% and pion misidentification of % 15%. The two tracks are assigned the kaon mass hypothesis and their invariant mass is required to be within 6:5 MeV=c 2 of the nominal mass [25] (the resolution is ¼ 1:0 MeV=c 2 and the natural width is À ¼ 4:3 MeV). We also require that the helicity angle H between the flight direction of one of the two kaons and the D flight direction, in the rest frame, satisfies the condition j cos H j > 0:4. This requirement exploits the fact that in D ! K 0 S decays the is produced in a longitudinally polarized state, thus cos H follows a cos 2 H distribution, while in candidates not from D ! K 0 S decays, cos H is approximately uniformly distributed.
The ! candidates are reconstructed from þ À 0 combinations with invariant mass within 17 MeV=c 2 (2À ! ) of the nominal ! mass [25] (the resolution is ¼ 6:9 MeV=c 2 ). The charged pion candidates are required to pass pion identification criteria with pion selection efficiency around 98% and kaon misidentification rate around 12%. To improve the ! momentum resolution, the invariant mass of the two photons forming the 0 candidate is constrained to the nominal 0 mass. We define N as the angle between the normal to the ! decay plane and the D momentum in the ! rest frame, and as the angle between the flight direction of one of the three pions in the ! rest frame and the flight direction of one of the other two pions in the two-pion rest frame. The quantities cos N and cos follow cos 2 N and ð1 À cos 2 Þ distributions, respectively, for the signal, and are almost uniformly distributed for wrongly reconstructed ! candidates. We require the product cos 2 N sin 2 > 0:046.
Neutral D candidates are formed from two-body combinations of K AE , AE , K 0 S , 0 , and ! candidates consistent with one of the six D decay channels under study. To improve the D CPÀ momentum resolution, the invariant masses of the 0 and K 0 S daughters are constrained to the nominal 0 and K 0 S masses. To suppress poorly reconstructed D candidates and candidates from random combinations, we perform a geometric fit of the D daughters to a common origin, and reject D candidates for which the 2 probability of the vertex fit is lower than 0.01%. The invariant mass of a D candidate M D must be within a range that corresponds to slightly more than twice the M D resolution: the range varies from about AE6 MeV=c 2 for the K 0 S channel to about AE44 MeV=c 2 for the K 0 S 0 channel. We apply the following particle identification criteria to the charged daughters of the D meson: in D ! þ À , the two-pion candidates must pass the same pion identification criteria adopted in the reconstruction of ! ! À þ 0 ; in D ! K þ K À , the two kaon candidates are required to pass tighter kaon identification criteria than those applied to the daughters (typical kaon selection efficiency around 94%, and pion misidentification rate around 6%); in D ! K À þ , the kaon candidate must pass the same kaon identification criteria required for the daughters. In order to reduce the large combinatorial background from random combinations of tracks and photons in e þ e À ! q " q events (q ¼ u, d, s, c), we put requirements on the cosine of the D decay angle, j cos D j. We define D as the angle between one of the D daughters in the D rest frame, and the direction of the D meson in the B rest frame. Because of angular momentum conservation we expect the distribution of cos D to be uniform for B AE ! Dh AE , D ! þ À and D ! K 0 S 0 signal events, while for q " q events the distribution is strongly peaked at AE1. We require j cos D j < 0:74 (0.99) for the B AE ! Dh AE , D ! þ À (D ! K 0 S 0 ) channel. The invariant mass distributions of the reconstructed D candidates, after all the other selection criteria described in this section have been applied, are shown in Fig. 1.
We reconstruct B AE meson candidates by combining a neutral D candidate with a track h AE . For the D ! K mode, the charge of the track h must match that of the kaon from the D meson decay. This selects b ! c mediated B decays B À ! D 0 h À and B þ ! " D 0 h þ . The contamination from b ! u mediated B decays followed by doubly-Cabibbo-suppressed D decay, i.e. B À ! " D 0 K À , " D 0 ! K À þ , and from D 0 -" D 0 mixing is negligible. In the B AE ! Dh AE , D ! þ À channel we require that the invariant mass of the ðh AE Ç Þ system is greater than 1:9 GeV=c 2 to reject background from B À ! D 0 À , D 0 ! K À þ and B À ! K Ã0 À , K Ã0 ! K À þ decays and their CP conjugates. Here is the pion from the D and h is the track from the B candidate taken with the kaon mass hypothesis. To improve the B momentum resolution, the neutral D invariant mass is constrained to the nominal D 0 mass [25] for all D decay channels.
We identify signal B ! DK and B ! D candidates using two kinematic variables: the difference between the CM energy of the B meson (E Ã B ) and the beam energy, and the beam-energy-substituted mass, where ðE B ; p B Þ and ðE ee ; p ee Þ are the four-momenta of the B meson and of the initial e þ e À system, respectively, measured in the laboratory frame. The m ES distributions for B AE ! Dh AE signals are centered at the B mass [25], have a root-mean-square of approximately 2:6 MeV=c 2 , and do not depend strongly on either the D decay mode or the nature of the track h. In contrast, the ÁE distributions depend on the mass assigned to the track h. We evaluate ÁE with the kaon mass hypothesis so that the peaks of the distributions are centered near zero for B AE ! DK AE events and are shifted by approximately þ50 MeV for B AE ! D AE events. The ÁE resolution depends on the kinematics of the decay, and is typically 16 MeV for all D decay modes under study after the D invariant mass is constrained to its nominal value. We retain B candidates with m ES and ÁE within the intervals 5:20 < m ES < 5:29 GeV=c 2 and À80 < ÁE < 120 MeV, which define the region for the fit described later.
In order to discriminate the signal from e þ e À ! q " q background events, denoted q " q in the following, we construct a Fisher discriminant F based on the four eventshape quantities L ROE 20 , j cos Ã T j, j cos Ã B j and H ROE 20 . These quantities, evaluated in the CM frame, are defined as (i) L ROE 20 ¼ L 2 =L 0 is the ratio of the second and zeroth event shape moments of the energy flow in the rest of event (ROE), i.e. considering all the charged tracks and neutral clusters in the event that are not used to reconstruct the B candidate. They are defined as where p i are the momenta and i the angles of the charged and neutral particles in the ROE, with respect to the thrust axis of the B candidate's decay products. The thrust axis is defined as the direction that maximizes the sum of the longitudinal momenta of the particles used to define it; (ii) Ã T is the angle between the thrust axis of the B candidate's decay products and the beam axis; (iii) Ã B is the angle between the B candidate momentum and the beam axis; (iv) H ROE 20 ¼ H 2 =H 0 is the ratio of the second and zeroth Fox-Wolfram moments H 2 and H 0 [26], computed using charged tracks and photons in the ROE.
The quantity F is a linear combination of the four aforementioned event-shape variables: The values of the coefficients c i are the ones which maximize the separation between simulated signal events and a continuum background sample provided by off-resonance data, taken % 40 MeV below the Çð4SÞ resonance. The maximum likelihood fit described in Sec. V is restricted to events with F within the interval À1:5 < F < 1:5, to remove poorly reconstructed candidates. For events with multiple B AE ! Dh AE candidates (about 16% of the selected events), we choose the B candidate with the smallest 2 c Þ formed from the measured and true masses, M c and M PDG c , of all the unstable particles c produced in the B decay tree (D, 0 , K 0 S , , !), scaled by the sum in quadrature of the resolution M c of the reconstructed mass and the intrinsic width À c . From simulated signal events, we find that this algorithm has a probability to select the correct candidate between 98.2% and 99.9% depending on the D decay mode. We also find that the algorithm has negligible effect on the M D distributions.
We compare the distribution of each selection variable in data and simulated events after the requirements on all other variables have been applied. In order not to introduce biases that may artificially enhance the signal yield, we perform a blind study by explicitly removing, in this comparison, events consistent with the B AE ! DK AE signal, i.e. those with jm ES À m B j < 10 MeV=c 2 , jÁEj < 40 MeV, F > À0:8 and track h passing kaon identification criteria. We find excellent agreement between data and simulated events, both for events consistent with the B AE ! D AE signal (jm ES À m B j < 10 MeV=c 2 , jÁE À 50 MeVj < 40 MeV, F > À0:8 and track h failing the kaon identification criteria) and for backgroundlike events.
We correct for small differences in the means and widths of the distributions of the invariant masses of the unstable particles and of m ES and ÁE both when applying to data the selection criteria obtained from simulated events and in the final fit described in the next section.
The total reconstruction efficiencies, based on simulated B AE ! DK AE events, are summarized in the second column of Table II. For the reasons explained in Sec. II, the efficiencies are 40% to 60% higher than in our previous study of the same decay channels [9]. The efficiencies obtained for B AE ! D AE events from the simulation are TABLE II. Reconstruction efficiency for B ! DK from simulated events. We also quote the efficiency and purity in a signalenriched subsample (see text for details).
statistically consistent with those for B AE ! DK AE , where the D meson is reconstructed in the same final state. For illustration purposes we define a signal-enriched sample for each D decay mode, containing all B AE ! Dh AE candidates satisfying the criteria À40 < ÁE < 100 MeV, 0:2 < F < 1:5, 5:275 < m ES < 5:285 GeV=c 2 , and whose daughter track h passes charged kaon identification criteria. The typical kaon efficiency is % 77% and the pion misidentification rate is % 2%. The reconstruction efficiencies and the expected purities for the signal-enriched subsamples, determined on simulated data, are listed in Table II.

V. MAXIMUM LIKELIHOOD FIT
We measure R ðAEÞ K= and A CPAE using simultaneous extended and unbinned maximum likelihood fits to the distributions of the three variables ÁE, m ES , and F of B candidates selected in data. The data set is split into 24 subgroups by means of three discrete variables: the charge ¼ AE1 of the reconstructed B meson ( Â 2 subgroups); the two-body D decay final state X ( Â 6), allowing for a more accurate description of the corresponding probability density functions compared to the larger CPAE subgroups; and a PID variable denoting whether or not the track h from the B passes (p) or fails (f) charged kaon identification criteria ( Â 2). The pion misidentification rate of these criteria is determined directly from data as described later, and is expected from simulation to be around 2%. The corresponding kaon identification efficiency is ð77 AE 1Þ%, as determined from the signal MC samples after weighting the bidimensional distribution of the momentum and polar angle of the track h by the ratio of the analogous distributions observed in MC and data kaon control samples. The uncertainty on the kaon identification efficiency is dominated by the systematic contribution from the uncertainties on the weights. We perform in total three simultaneous fits to these 24 subgroups: one fit for the two CP-even D final states (8 subgroups), one for the three CP-odd D final states (12 subgroups), and one for the D ! K decay (4 subgroups).
The likelihood function L for each of these simultaneous fits has the form where s ranges over the subgroups under consideration, N s is the number of events in subgroup s, n is the total number of events in the fit n ¼ P s N s , and N is the expected number of events. We minimize À lnL with respect to the set of fit parameters specified later. The probability P s;i P s ðm ESi ; ÁE i ; F i Þ for an event i is the sum of six signal and background components: B AE ! DK AE signal, B AE ! D AE signal, background candidates from e þ e À ! q " q events, irreducible background arising from charmless B AE ! XK AE and B AE ! X AE decays, and background candidates from other B " B events (reducible B " B background): where the N j s are the expected yields in each component j. In case of negligible correlations among the fit variables, each probability density function (PDF) P factorizes as The irreducible B " B background originates from events where a B meson decays to the same final state Xh as the signal, but without the production of an intermediate charmed meson in the decay chain. When exploiting the ÁE, m ES , and F variables, this background is therefore indistinguishable from the signal. As an example, the decay As described later in Sec. VI, the irreducible background yield can be estimated by studying sideband regions of the D candidate invariant mass distribution, and can then be fixed in the final fit, where we assume P Dh i ¼ P Xh i . We express the signal yield parameters N DK s and N D s through the CP asymmetries A X DK and A X D of B AE ! DK AE , D ! X and B AE ! D AE , D ! X, their branching fraction ratios, R X K= , the total number N D tot;X of B AE ! D AE , D ! X signal events, the true kaon identification efficiency " of the PID selector, and the pion misidentification rate m of the PID selector: N DK ;f;X ¼ Because the ratios R X K= are small, the fit is not able to determine the value of ". Therefore we fix it to the aforementioned value of " ¼ ð77 AE 1Þ%. The reconstruction and selection efficiencies for true B AE ! DK AE and B AE ! D AE candidates, where the D meson decays to the same final state, are assumed to be identical. A systematic uncertainty is assigned due to this assumption (see Sec.VII). The simultaneous fit to the two CP-even modes constrains while the simultaneous fit to the three CP-odd modes constrains The m ES distributions of the signal components are parametrized using an asymmetric Gaussian shape, i.e. a Gaussian with different widths on both sides of the peak. We use the same shape for B AE ! DK AE and B AE ! D AE , so the m ES B AE ! DK AE signal shape (whose parameters are floating in the fit) will mostly be determined by the much more abundant B AE ! D AE control sample. Since the selection efficiencies for the two channels are the same, we expect the number of reconstructed candidates from B AE ! D AE to be about 12 times higher than for B AE ! DK AE . We have checked that the m ES shapes for B AE ! DK AE and B AE ! D AE are consistent, and that the assumption that they are identical does not bias the parameters of interest.
The ÁE distribution of the B AE ! DK AE signal component is parametrized with a double Gaussian shape. The core Gaussian has a mean close to zero, a width around 16 MeV and, according to the simulation, accounts for about 90% of the true B AE ! DK AE candidates. The second Gaussian accounts for the remaining 10% of candidates whose energy has been poorly measured. The mean and width of the core Gaussian are directly determined from data, while the remaining three parameters (the difference between the two means, the ratio between the two widths and the ratio of the integrals of the two Gaussian functions) are fixed from the simulation. In contrast to the m ES case, the B AE ! D AE ÁE shape is not the same as for B AE ! DK AE . This is due to the fact that we always assign the kaon mass hypothesis to the track h: the wrong mass assignment, in the case of B AE ! D AE , introduces a shift to the reconstructed energy of the pion and thus to ÁE, since ÁE ¼ The shift depends on the magnitude of the momentum p of the track h in the laboratory frame, We parametrize the m ES -ÁE distribution of the B " B background by means of two factorizing components: The m ES component of the peaking part, g peak ðm ES Þ, is parametrized with a Gaussian function for X ¼ þ À , K 0 S !, K 0 S . For X ¼ K þ K À , K 0 S 0 we use the ''Crystal Ball'' lineshape [27], an empirical smooth function that better describes the non-Gaussian tail on the negative side of the distribution, n n jj n e Àðjj 2 =2Þ n jj À jj À " x Àn " x < Àjj; with " x ¼ ðx À Þ= and " x ! À" x for < 0. For X ¼ K we use an empirical function of the form: Here is the position of the peak, while and are parameters related to the width of the distribution on the two sides of the peak. The ÁE component of the peaking part h peak ðÁEÞ is described with a simple exponential function for the five CP self-conjugate D final states, and with a Landau function for the non-CP-eigenstate final state. The B " B purely combinatorial background component is described by the 2dimensional product of a linear background, h cont ðÁEÞ, and an empirical function introduced by the ARGUS Collaboration [28], g cont ðm ES Þ ¼ Aðm ES =m 0 Þ: where m 0 ¼ ffiffi ffi s p =ð2c 2 Þ ¼ 5:29 GeV=c 2 is the kinematic endpoint of the m ES distribution. All the parameters of the B " B background m ES -ÁE distribution are fixed from simulated B " B events. The only exception is the width of the Landau function used for h peak ðÁEÞ in X ¼ K À þ . This parameter controls the behavior at low ÁE values, ÁE % À80 MeV, where we find the simulation not to be sufficiently precise given the high statistics of this channel. We note that the shape parameters differ across the six final states, but are similar across the charge and PID selector subgroups belonging to one final state.
In q " q events, B candidates arise from random combinations of charged tracks and neutral particles produced in the hadronization of the light quark-antiquark pairs produced in e þ e À collisions. Similarly to the combinatorial component of the B " B background, the q " q background distribution in the m ES -ÁE plane is parameterized by the product of an ARGUS function in m ES and a linear background in ÁE. We float the slope of the linear components, while the parameters of the ARGUS function are fixed, in each D final state, from simulated q " q events. They are in good agreement across the final states and other subgroups.
The F distributions are parametrized in a similar way for all fit components. We find that the distributions of B AE ! DK AE and B AE ! D AE signal events are consistent with each other, as expected since their kinematics are very similar, and choose to parametrize them with the same shape. For this we use the sum of two asymmetric Gaussian functions. Some channels with lower statistics do not require the full complexity of this parametrization: in those cases we use a single asymmetric Gaussian, a double Gaussian, or a single Gaussian. In particular we use: for the signal components a double asymmetric Gaussian, except for X ¼ K 0 S , where a double Gaussian function is adopted; for the B " B background components a double asymmetric Gaussian in case of X ¼ K 0 S 0 , K À þ , a double Gaussian in case of X ¼ K 0 S !, and a single Gaussian otherwise; for the q " q background components a double asymmetric Gaussian, except for X ¼ K 0 S , where we use a single Gaussian.
In summary, the floating parameters of the fits are all parameters related to the signal yields, and therefore to the GLW parameters, as given in Eqs. (14)- (17), except "; all background yields and CP-asymmetries except the irreducible background yields and asymmetries, the B " B asymmetries for CPÀ modes and for the ðCPþ; pÞ subgroups (B ! D CPþ h candidates where the track h passes the kaon identification criteria), and the B " B yield in the ðK 0 S ; pÞ subgroup; selected shape parameters, namely, the overall width and mean of the ÁE signal, the m ES signal shape, and the ÁE and F shape for q " q background. The nonfloating parameters are fixed to their expectations obtained from simulation or, in case of the irreducible background yields, to values obtained from data control samples (see next section). Nonfloating CP asymmetries are fixed to zero. We assign systematic uncertainties due to the fixed parameters.
We check that the fitter is correctly implemented by generating and fitting a large number of test data sets using the final PDFs. In this study, we include an analytic description for the conditional variable ÁE shift . The residuals for a given parameter, divided by the measured parameter error, should follow a Gaussian distribution with zero mean () and unitary width (). We observe no significant deviations from the expected distribution. In particular, R þ K= shows the largest shift from zero mean ( ¼ À0:06 AE 0:07) and A CPþ shows the largest deviation from unity width ( ¼ 1:13 AE 0:06) among the parameters of interest.
We investigate fit biases, arising from possible discrepancies between the true signal distribution and the chosen fit model, by fitting a large number of test data sets, in which the B AE ! DK AE and B AE ! D AE signal components are taken from simulated samples of sufficient statistics, while the background components are randomly generated according to their PDFs. Of all floating parameters, only R K= acquires a significant bias, resulting in corrections of 0.5 and 1.0 times the expected statistical uncertainties on these parameters in the CP and flavor modes, respectively. This bias is caused by small differences in the ÁE distributions of the signal components across the kaon PID subgroups (p and f), which the final PDF does not account for. A second, smaller contribution to this bias is a small discrepancy between the ÁEðÞ signal shape of B AE ! D AE events and the ÁEðKÞ shape of B AE ! DK AE events. The biases in the R K= parameters are correlated, and partly cancel in the ratio, resulting in a smaller bias on the GLW parameters R CPAE . The largest (smallest) remaining bias is 0.12 (0.05) times the expected statistical uncertainty for R CPþ (A CPÀ ). We correct the final values of the parameters A CP and R K= for the observed biases, and assign systematic uncertainties to these corrections.

VI. IRREDUCIBLE BACKGROUND DETERMINATION
As discussed in the previous section, the irreducible background arises from charmless B AE ! Xh AE decays, which have the same final states as the B AE ! Dð! XÞh AE signal and therefore the same distribution of the three fit variables ÁE, m ES , and F .
In the D 0 ! K À þ flavor mode, the irreducible background-taking into account the measured branching fractions for B AE ! K AE Ç K AE and B AE ! K AE Ç AE [25] and a selection efficiency of % 1%, estimated from simulated events-is negligible compared to the expected signal yields (about 3400 B AE ! DK AE and 45000 B AE ! D AE expected signal events). On the other hand, in the CP modes, where the signal yields are expected to be an order of magnitude lower than in K À þ , and the upper limits for the branching ratios of B AE ! Xh AE decays are at the 10 À5 level, we cannot a priori exclude a relevant irreducible background contribution.
We estimate the irreducible background yields in our sample by exploiting the fact that the D invariant mass distribution for this background is approximately uniform, while for the signal it is peaked around the nominal D mass. Therefore we can select a control sample containing irreducible background candidates, but with the signal strongly suppressed, by applying the same selection as for the signal, with the only difference that the D invariant mass is required to lie in a region (D invariant mass sidebands) which is separated by at least a few M D from the nominal D mass (see Table III). We then perform an extended maximum likelihood fit to the m ES , ÁE, and F distributions of the control sample in order to measure the irreducible background yields in the D invariant mass sidebands. The fit is similar to the nominal one described in the previous section. However, due to the limited statistics available in the sidebands, we are forced to fix more parameters compared to the nominal fit; in particular, we fix any possible charge asymmetry of the B AE ! Xh AE decays to zero (a systematic uncertainty is assigned to this assumption). Finally, since the D candidate invariant mass distribution of the irreducible background is approximately uniform, we scale the obtained yields by the ratio of the widths of the D signal and control sideband mass regions to obtain the irreducible background yield N Xh (scale factor in Table III). Table IV shows the scaled irreducible background yields that enter the final fit.

VII. SYSTEMATIC UNCERTAINTIES
We consider nine sources of systematic uncertainty that may affect the GLW parameters A CPAE and R CPAE . Their contributions are summarized in Table V. First, we estimate the influence of fixed parameters of the nominal PDF. We perform a large number of test fits to the data, similar to the nominal fit. In each of these test fits the fixed parameters are varied according to their covariance matrices. From the resulting distributions we calculate the systematic covariances of the fit parameters A CPAE and R K= . The parameters responsible for the largest uncertainty are the m ES endpoint m 0 , and parameters related to the measured yields, e.g. BB background asymmetries and the efficiency of the kaon selector.
The uncertainties in the irreducible background event yields introduce a systematic uncertainty in the B AE ! D CP h AE yields and therefore in R CPAE . Likewise, any charge asymmetry in this background would affect the measured values of A CPAE . We again perform a series of test fits to onpeak data, where we vary the B AE ! Xh AE yields and asymmetries by their uncertainties. For the latter, we take the uncertainties to be AE10% for X ¼ K þ K À and AE20% for the other modes, which are conservative estimates consistent with the existing upper limits on the CP asymmetries in those decays [21].
As explained in Sec. V, we correct the fit results for biases observed in Monte Carlo studies. We take the  associated systematic uncertainties to be half the size of the bias corrections, summed in quadrature with the statistical uncertainties on the biases. The latter are due to the limited number of test fits used to estimate the corrections. We investigate a potential charge asymmetry of the BABAR detector, due to a possible charge bias in tracking efficiency (e.g. K þ vs K À ) and/or particle identification. Our analysis includes a number of control samples, in which the CP asymmetry is expected to be negligible: the six B AE ! D AE samples and the B AE ! DK AE flavor mode (D ! K). The weighted average of the charge asymmetry in the control samples is ðÀ0:95 AE 0:44Þ%, from which we assign uncertainties of 1.4% to both A CPþ and A CPÀ . We consider these uncertainties to be 100% correlated.
The measured CP asymmetry in B AE ! DK AE , D ! K 0 S , can be diluted by the presence of B AE ! DK AE decays followed by D decays to the same final state K 0 S K þ K À as the signal but with opposite CP content, such as D ! K 0 S a 0 , a 0 ! K þ K À . The same can happen in the This background can also affect the ratios R CPÀ . It is possible to obtain correction factors to both A CPÀ and R CPÀ from a fit to the distributions of the relevant helicity angles, cos N and cos H for K 0 S ! and K 0 S , respectively. The fit is performed on dedicated B AE ! D AE samples, in which the selection requirements on the helicity angles have not been applied. It can be shown [29] that for these two final states the observed charge asymmetries and ratios should be corrected by a factor Here, R 0 is the ratio of the R AE K= values, where R À K= is taken from a single fit to the D 0 ! K 0 S 0 final state only (as opposed to using all three CPÀ final states under The curves are the full PDF (solid, blue), and B ! D (dash-dotted, green) stacked on the remaining backgrounds (dotted, purple). The region between the solid and the dash-dotted lines represents the B ! DK contribution. We show the subsets of the data sample in which the track h from the B decay is identified as a kaon. We require candidates to lie inside the signal-enriched region defined in Sec. IV, except for the plotted variable. angles: f ;K 0 S ! ¼ 0:71 and f ;K 0 S ¼ 0:64. To apply these corrections, we first perform a fit of the K 0 S 0 final state alone to obtain R K 0 S 0 K= . We then perform the simultaneous fit of the CPþ final states, from which we take the value of R þ K= . Finally, we include the correction factors into the CPÀ final PDF, which will allow the likelihood fitter to correctly estimate their influence. The parameter jzj 2 in Eqs. (27) and (28) is extracted from fits of the helicity angle distributions in the D 0 ! K 0 S ! and D 0 ! K 0 S subsamples to the function jzj 2 þ 3cos 2 [29]. We subtract the background expected from the Monte Carlo simulation, which has been rescaled to match the data. We find jzj 2 ¼ 0:065 AE 0:033 in the case of K 0 S !, and jzj 2 ¼ 0:217 AE 0:063 in the case of K 0 S . The uncertainties contain propagated uncertainties due to the background subtraction. The resulting corrections are In order to assign systematic uncertainties, we propagate the uncertainties on the correction factors into the final result.
When calculating R CP through Eq. (5) one has to take into account that this equation is an approximation. We define the double ratios used to approximate R CPAE as R AE . They are given by where D f denotes the K À þ final state. These can be written as  3 (color online). m ES projections of the fits to the data, split into subsets of definite CP of the D candidate and charge of the B candidate: (a) B À ! D CPþ K À , (b) B þ ! D CPþ K þ , (c) B À ! D CPÀ K À , (d) B þ ! D CPÀ K þ . We show the subsets of the data sample in which the track h from the B decay is identified as a kaon. See caption of Fig. 2 for line definitions. Only a subrange of the whole fit range is shown in order to provide a closer view of the signal peak.
where r B and B are defined, in analogy to r B and B , as r B e ið B ÀÞ ¼ AðB À ! " D 0 À Þ=AðB À ! D 0 À Þ, while r D and D are defined as r D e i D ¼ Að " D 0 ! K À þ Þ=AðD 0 ! K À þ Þ. We write Eq. (34) in the form R AE ¼ R CPAE Â ð1 þ R c Þ, and we assign a relative systematic uncertainty based on the value of the correction R c . Taking sin C ¼ 0:2257 AE 0:0010 (where C is the Cabibbo angle) and r B ¼ 0:104 þ0:015 À0:025 from [3], and expressing r D ¼ jV cd V us j=jV ud V cs j ¼ tan 2 C , and r B ¼ r B tan 2 C , we find R c % 4r B tan 2 C % 2:2%. Here, we have conservatively assumed values for the cosine terms which maximize R c . We thus assign a relative uncertainty of 2.2% to the values of R CP , fully correlated between R CPþ and R CPÀ .
We also consider the influence on the measured value of A CP of misreconstructed signal B candidates, i.e. candidates reconstructed, in events containing a true B ! DK decay with D decaying to the same final state X as the reconstructed candidate, from random combinations of particles produced in the true B ! DK decay and the particles of the ROE. The fraction of these candidates ranges from 0.3% to 12% in simulated B AE ! D CP K AE events, depending on the channel. Since we treat this component like the signal, we implicitly assume that its charge asymmetry is equal to the asymmetry in the signal component. We use simulated signal events to estimate the ratio between misreconstructed and true B þ ! DK þ candidates and the ratio between misreconstructed and true B À ! DK À candidates, and find these two quantities to differ by less than 0.1%, from which we derive an upper limit on the difference between the observed and the true value of A CP .
The yield double ratios R CPAE should be corrected by the corresponding double ratio of selection efficiencies. We find from simulated events that the efficiency double ratios are compatible with each other, and their average value is very close to unity, ð99:46 AE 0:23Þ%. Thus we do not correct the central values but conservatively assign a relative uncertainty equal to 1 À ð0:9946 À 0:0023Þ ¼ 0:0077.
The final PDF does not contain an explicit description of the conditional parameter ÁE shift , assuming implicitly that the distribution of ÁE shift observed in data is the same for all the components of the fit. However, the distributions are found to be slightly different across the components, thus introducing a possible bias in the fit results. To estimate the size of this bias, we use simulated events to obtain parametrizations of the ÁE shift distributions of all the fit components and repeat the fits to data. We assign the differences compared to the results of the nomi- nal fits as the systematic uncertainty. We expect this effect to be highly correlated between A CP parameters, because the PDFs are similar in each D decay channel. Thus they are affected by nonuniform ÁE shift distributions in a similar way. The same argument holds for the R K= parameters. We studied the effect of assigning a 0%, 50%, and 100% correlation. The uncorrelated case gave the largest deviations from the nominal results, the fully correlated case gave the smallest. However, the variation was found to be at the 10% level. We assign the systematic uncertainty corresponding to a correlation of 50%. Table V lists the contributions of the effects discussed above. Compared to our previous analysis [9], the systematic uncertainty on A CPþ is reduced due to better understanding of the detector intrinsic charge asymmetry (the determination of which benefits from the larger data set) and due to improved evaluation of the correlations among the different sources of systematic uncertainties. The uncertainty on A CPÀ is only slightly reduced. By contrast, the systematic uncertainties on R CPAE are increased due to two additional sources of uncertainty that were not considered previously: the bias correction and the differences of the ÁE shift distributions among the fit components. The systematic correlations between the GLW parametersỹ ¼ ðA CPþ ; A CPÀ ; R CPþ ; R CPÀ Þ T are

VIII. RESULTS
The signal yields returned from the fit for each of the D decay modes under study are listed in Table VI. We reconstruct almost 1000 B AE ! D CP K AE decays and about 4 times more B AE ! DK AE , D ! K decays.
The final values of the GLW parameters that we measure are The statistical correlations among these four quantities are The results are in good agreement with those from our previous analysis [9] and the current world averages [21]. The statistical significance of a nonzero A CPþ value is determined from the maximum value of the likelihood function of the nominal fit and that of a dedicated nullhypothesis fit, where A CPþ was fixed to zero, Taking into account systematic uncertainties, the statistical significance of A CPþ is slightly decreased to E (GeV) ∆ This constitutes evidence for direct CP violation in charged B decays and the first evidence of direct CP violation in B ! DK.
We constrain the CKM angle , the strong phase B , and the amplitude ratio r B from the present measurement by adopting the frequentist procedure also exploited in [15]. We define a multivariate Gaussian likelihood function relating the experimentally measured observablesỹ and their statistical and systematic covariance matrices V cov ¼ V stat þ V syst with the corresponding truth parametersỹ t 1 y t ð; B ; r B Þ calculated using Eqs. (3) and (4). The matrices V stat and V syst are constructed from Eqs. (35)-(40). The normalization is N ¼ ð2Þ 2 ffiffiffiffiffiffiffiffiffiffiffiffi jV cov j p . We then define a 2 -function as 2 ð; B ; r B Þ ¼ À2 lnLð; B ; r B Þ: Because of the inherent eight-fold ambiguity of the GLW method there are eight equivalent minima of the 2 -function, 2 min , which correspond to the same value of r B and to eight alternative solutions for ð; B Þ. To evaluate the confidence level of a certain truth parameter (for example ) at a certain value ( 0 ) we consider the value of the 2 -function at the new minimum, 2 min ð 0 ; 0 B ; r 0 B Þ, satisfying Á 2 ¼ 2 min ð 0 ; 0 B ; r 0 B Þ À 2 min ! 0. In a purely Gaussian situation for the truth parameters the CL is given by the probability that Á 2 is exceeded for a 2 -distribution with one degree of freedom: A more accurate approach is to take into account the nonlinearity of the GLW relations, Eqs. (3) and (4)   In order to facilitate the future combination of these measurements with the results of the Dalitz plot analysis of B AE ! DK AE , D ! K 0 S h þ h À decays (h ¼ , K) [16], we recompute the GLW parameters after excluding from the nominal fit the D CPÀ ! K 0 S ( ! K þ K À ) subsample. The sample obtained in this way is statistically independent of that selected in [16].
The statistical correlations among these four quantities are   n/a n/a n/a n/a n/a 0:01048 AE 0:00057 a q " To compare the results obtained after removing the D CPÀ ! K 0 S subsample with those from the B AE ! DK AE , D ! K 0 S h þ h À analyses, which are expressed in terms of the variables x AE ¼ r B cosð B AE Þ and y AE ¼ r B sinð B AE Þ, we use the GLW parameters measured in this way to determine the quantities x AE through the relations: We obtain x þ ¼ À0:057 AE 0:039ðstatÞ AE 0:015ðsystÞ; (54) x À ¼ 0:132 AE 0:042ðstatÞ AE 0:018ðsystÞ: These results are in good agreement with the current world averages [21] and have precision close to the single most precise measurements [16]. We also measure r 2 B , which provides a constraint on x AE and y AE via r 2 B ¼ x 2 AE þ y 2 AE , from We determine: The constraints that could be placed on the quantities y AE from these measurements, by exploiting the relation r 2 B ¼ x 2 AE AE y 2 AE , are much weaker than those provided by the B AE ! DK AE , D ! K 0 S h þ h À analysis. As a final check of consistency we consider the quantity a, a ¼ A CPþ R CPþ þ A CPÀ R CPÀ : From Eqs. (3) and (4) one expects a to satisfy a ¼ 0. We measure a ¼ 0:19 AE 0:11ðstat þ systÞ, which is compatible with 0.

IX. SUMMARY
Using the entire data set collected by BABAR at the e þ e À center-of-mass energy close to the Çð4SÞ mass, we have reconstructed B AE ! DK AE decays, with D mesons decaying to non-CP (K), CP-even (K þ K À , þ À ) and CP-odd (K 0 S 0 , K 0 S , K 0 S !) eigenstates. Through an improved analysis method compared to the previous BABAR measurement [9] and through an enlarged data set, corresponding to an increase in integrated luminosity at the Çð4SÞ peak from 348 fb À1 to 426 fb À1 , we obtain the most precise measurements of the GLW parameters A CPAE and R CPAE to date: A CPþ ¼ 0:25 AE 0:06ðstatÞ AE 0:02ðsystÞ; A CPÀ ¼ À0:09 AE 0:07ðstatÞ AE 0:02ðsystÞ; R CPþ ¼ 1:18 AE 0:09ðstatÞ AE 0:05ðsystÞ; R CPÀ ¼ 1:07 AE 0:08ðstatÞ AE 0:04ðsystÞ: We measure a value of A CPþ which is 3.6 standard deviations from zero, which constitutes the first evidence for direct CP violation in B ! DK decays.
From the measured values of the GLW parameters, we extract confidence intervals for the CKM angle , the strong phase B , and the amplitude ratio r B , using a frequentist approach, taking into account both statistical and systematic uncertainties. At the 68% CL we find that both and B (modulo 180 ) belong to one of the three intervals [11.3 , 22.7 ] [25].
To facilitate the combination of these measurements with the results of our Dalitz plot analysis of B AE ! DK AE , D ! K 0 S h þ h À ðh ¼ K; Þ [16], we exclude the D ! K 0 S , ! K þ K À channel from this analysis-thus removing events selected also in [16]-and then determine A CPÀ ¼ À0:08 AE 0:07ðstatÞ AE 0:02ðsystÞ; R CPÀ ¼ 1:03 AE 0:09ðstatÞ AE 0:04ðsystÞ: For comparison with the results of the B AE ! DK AE , D ! K 0 S h þ h À analyses, which are expressed in terms of the variables x AE ¼ r B cosð B AE Þ and y AE ¼ r B sinð B AE Þ, we express our results for the GLW observables in terms of x þ and x À . We measure x þ ¼ À0:057 AE 0:039ðstatÞ AE 0:015ðsystÞ; x À ¼ 0:132 AE 0:042ðstatÞ AE 0:018ðsystÞ; at 68% CL. These results are in good agreement with the current world averages [21] and have precision comparable to the single most precise measurements [16]. We also evaluate r B after the exclusion of the D ! K 0 S channel, and obtain a weak constraint on r 2 B ¼ x 2 AE AE y 2 AE : r 2 B ¼ 0:105 AE 0:067ðstatÞ AE 0:035ðsystÞ at 68% CL.