Measurement of time dependent CP asymmetry parameters in B0 meson decays to ω KS0, ν K0, and n0 KS0

We present measurements of the time-dependent CP -violation parameters S and C in the decays B 0 ! !K 0 S , B 0 ! (cid:1) 0 K 0 , reconstructed as (cid:1) 0 K 0 S and (cid:1) 0 K 0 L , and B 0 ! (cid:2) 0 K 0 S . The data sample corresponds to the full BABAR dataset of 467 (cid:1) 10 6 B (cid:1) B pairs produced at the PEP-II asymmetric-energy e þ e (cid:2) collider at the Stanford Linear Accelerator Center. The results are S !K 0 S ¼ 0 : 55 þ 0 : 26 (cid:2) 0 : 29 (cid:3) 0 : 02 , C !K 0 S ¼ (cid:2) 0 : 52 þ 0 : 22 (cid:2) 0 : 20 (cid:3) 0 : 03 , S (cid:1) 0 K 0 ¼ 0 : 57 (cid:3) 0 : 08 (cid:3) 0 : 02 , C (cid:1) 0 K 0 ¼ (cid:2) 0 : 08 (cid:3) 0 : 06 (cid:3) 0 : 02 , S (cid:2) 0 K 0 S ¼ 0 : 55 (cid:3) 0 : 20 (cid:3) 0 : 03 , and C (cid:2) 0 K 0 S ¼ 0 : 13 (cid:3) 0 : 13 (cid:3) 0 : 03 , where the ﬁrst errors are statistical and the second systematic. These results are consistent with our previous measurements and the world average of sin2 (cid:3) measured in B 0 ! J= c K 0 S .


I. INTRODUCTION
Measurements of time-dependent CP asymmetries in B 0 meson decays through b ! c " cs amplitudes have provided crucial tests of the mechanism of CP violation in the standard model (SM) [1]. These amplitudes contain the leading b-quark couplings, given by the Cabibbo-Kobayashi-Maskawa [2] (CKM) flavor mixing matrix, for kinematically allowed transitions. Decays to charmless final states such as 0K 0 , % 0 K 0 , 0 K 0 , !K 0 , K þ K À K 0 , f 0 ð980ÞK 0 are CKM-suppressed b ! q " qs (q ¼ u, d, s) processes dominated by a single loop (penguin) amplitude. This amplitude has the same weak phase ¼ argðÀV cd V Ã cb =V td V Ã tb Þ of the CKM-mixing matrix as that measured in the b ! c " cs transition, but is sensitive to the possible presence of new heavy particles in the loop [3]. Because of the different nonperturbative strong-interaction properties of the various penguin decays, the effect of new physics is expected to be channel dependent.
The CKM phase is accessible experimentally through interference between the direct decay of the B meson to a CP eigenstate and B 0 " B 0 mixing followed by decay to the same final state. This interference is observable through the time evolution of the decay. In the present study, we reconstruct one B 0 from Çð4SÞ ! B 0 " B 0 , which decays to the CP eigenstate !K 0 S , 0 K 0 S , 0 K 0 L , or % 0 K 0 S (B CP ). From the remaining particles in the event we also reconstruct the decay vertex of the other B meson (B tag ) and identify its flavor. The distribution of the difference Át t CP À t tag of the proper decay times t CP and t tag of these mesons is given by fðÁtÞ ¼ e ÀjÁtj=( 4( f1 AE ½À f S f sinðÁm d ÁtÞ where f is the CP eigenvalue of final state f ( À 1 for !K 0 S , 0 K 0 S , and % 0 K 0 S ; þ1 for 0 K 0 L ). The upper (lower) sign denotes a decay accompanied by a B 0 ð " B 0 Þ tag, ( is the mean B 0 lifetime, and Ám d is the mixing frequency. A nonzero value of the parameter C f would indicate direct CP violation. In these modes we expect C f ¼ 0 and À f S f ¼ sin2, assuming penguin dominance of the b ! s transition and neglecting other CKM-suppressed amplitudes with a different weak phase. However, these CKMsuppressed amplitudes and the color-suppressed tree diagram introduce additional weak phases whose contributions may not be negligible [4][5][6][7]. As a consequence, the measured S f may differ from sin2 even within the SM. This deviation ÁS f ¼ S f À sin2 is estimated in several theoretical approaches: QCD factorization (QCDF) [4,8], QCDF with modeled rescattering [9], soft collinear effective theory [10], and SU(3) symmetry [5,7,11]. The estimates are channel dependent. Estimates of ÁS from QCDF are in the ranges (0.0, 0.2), ðÀ0:03; 0:03Þ, and (0.01, 0.12) for !K 0 S , 0 K 0 , and % 0 K 0 S , respectively [8,10,12]; SU(3) symmetry provides bounds of ðÀ0:05; 0:09Þ for 0 K 0 and ðÀ0:06; 0:12Þ for % 0 K 0 S [11]. Predictions that use isospin symmetry to relate several amplitudes, including the I ¼ 3 2 B ! K% amplitude, give an expected value for S % 0 K 0 S near 1.0 instead of sin2 [13].
We present updated measurements of mixing-induced CP violation in the B 0 decay modes !K 0 S , 0 K 0 , and % 0 K 0 S , which supersede our previous measurements [14][15][16]. Significant changes to previous analyses include twice as much data for !K 0 S , 20% more data for 0 K 0 and % 0 K 0 S , improved track reconstruction, and an additional decay channel in 0 K 0 L . Despite the modest increase in data, the uncertainties on S 0 K 0 and C 0 K 0 decrease by 20% and 25%, respectively. Measurements in these modes have also been made by the Belle Collaboration [17,18].

II. THE BABAR DETECTOR AND DATASET
The results presented in this paper are based on data collected with the BABAR detector at the PEP-II asymmetric-energy e þ e À storage ring, operating at the Stanford Linear Accelerator Center. At PEP-II, 9.0 GeV electrons collide with 3.1 GeV positrons to yield a centerof-mass energy of ffiffi ffi s p ¼ 10:58 GeV, which corresponds to the mass of the Çð4SÞ resonance. The asymmetric energies result in a boost from the laboratory to the e þ e À center-ofmass frame of % 0:56. We analyze the entire BABAR dataset collected at the Çð4SÞ resonance, corresponding to an integrated luminosity of 426 fb À1 and ð467 AE 5Þ Â 10 6 B " B pairs. We use an additional 44 fb À1 of data recorded about 40 MeV below this energy (off-peak) for the study of the non-B 0 " B 0 background. A detailed description of the BABAR detector can be found elsewhere [19]. Surrounding the interaction point is a five-layer double-sided silicon vertex tracker (SVT) that provides precision measurements near the collision point of charged particle tracks in the planes transverse to and along the beam direction. A 40-layer drift chamber surrounds the SVT. Both of these tracking devices operate in the 1.5 T magnetic field of a superconducting solenoid to provide measurements of the momenta of charged particles. Charged hadron identification is achieved through measurements of particle energy loss in the tracking system and the Cherenkov angle obtained from a detector of internally reflected Cherenkov light. A CsI(Tl) electromagnetic calorimeter (EMC) provides photon detection, electron identification, and % 0 , , and K 0 L reconstruction. Finally, the instrumented flux return (IFR) of the magnet allows discrimination of muons from pions and detection of K 0 L mesons. For the first 214 fb À1 of data, the IFR was composed of a resistive plate chamber system. For the most recent 212 fb À1 of data, a portion of the resistive plate chamber system has been replaced by limited streamer tubes [20].

III. VERTEX RECONSTRUCTION
In the reconstruction of the B CP vertex, we use all charged daughter tracks. Daughter tracks that form a K 0 S are fit to a separate vertex, with the resulting parent momentum and position used in the fit to the B CP vertex. The vertex for the B tag decay is constructed from all tracks in the event except the daughters of B CP . An additional constraint is provided by the calculated B tag production point and three-momentum, with its associated error matrix. This is determined from the knowledge of the three-momentum of the fully reconstructed B CP candidate, its decay vertex and error matrix, and from the knowledge of the average position of the e þ e À interaction point and Çð4SÞ average boost. In order to reduce bias and tails due to long-lived particles, K 0 S and Ã 0 candidates are used as input to the fit in place of their daughters. In addition, tracks consistent with photon conversions ( ! e þ e À ) are excluded from the fit. To reduce contributions from charm decay products that bias the determination of the vertex position the tracks with a vertex 1 2 contribution greater than 6 are removed and the fit is repeated until no track fails the 1 2 requirement. We obtain Át from the measured distance Áz between the B CP and B tag vertex with the relation Áz ' cÁt.
Because there are no charged particles present at the B 0 ! % 0 K 0 S decay vertex, the % 0 K 0 S vertex reconstruction differs significantly from that of the !K 0 S and 0 K 0 analyses. In % 0 K 0 S we identify the vertex of the B CP using the single K 0 S trajectory from the % þ % À momenta and the knowledge of the average interaction point (IP) [21], which is determined several times per hour from the spatial distribution of vertices from two track events. The average transverse size of the IP is 180 "m Â 4 "m. We compute Át and its uncertainty with a geometric fit to the Çð4SÞ ! B 0 " B 0 system that takes this IP constraint into account. We further improve the accuracy of the Át measurement by constraining the sum of the two B decay times (t CP þ t tag ) to be equal to 2( (( is the mean B 0 lifetime) with an uncertainty ffiffiffi 2 p (, which effectively improves the determination of the decay position of the Çð4SÞ. We have verified in a full detector simulation that this procedure provides an unbiased estimate of Át. The estimate of the uncertainty on Át for each % 0 K 0 S event reflects the strong dependence of the Át resolution on the K 0 S flight direction and on the number of SVT layers traversed by the K 0 S decay daughters. When both pion tracks are reconstructed with information from at least the first three layers of the SVT in the coordinate along the collision axis (axial) as well as on the transverse plane (azimuthal), we obtain Át with resolution comparable to that of the !K 0 S and 0 K 0 analyses. The average Át resolution in these modes is about 1.0 ps. Events for which there is axial and azimuthal information from the first three layers of the SVT and for which Át and the error on Át satisfy jÁtj < 20 ps and ' Át < 2:5 ps are classified as ''good'' (class g), and their Át information is used in the time-dependent part of the likelihood function Eq. (10). About 60% of the events fall in this class. Otherwise events are classified as ''bad'' (class b). Since C f can also be extracted from flavor tagging information alone, events of class b contribute to the measurement of C f Eq. (11) and to the signal yield in the % 0 K 0 S analysis. In !K 0 S and 0 K 0 decays, the determination of the B decay vertex is dominated by the charged daughters of the ! and 0 , so we do not require information in the first three SVT layers from K 0 S daughter pions for events in class g. Also, since about 95% of events in these modes are of class g, the precision of the measurement of C f is not improved by including class b events. We maintain simplicity of these analyses by simply rejecting class b events.

IV. FLAVOR TAGGING AND Át RESOLUTION
In the measurement of time-dependent CP asymmetries, it is important to determine whether at the time of decay of the B tag the B CP was a B 0 or a " B 0 . This ''flavor tagging'' is achieved with the analysis of the decay products of the recoiling B tag meson. Most B mesons decay to a final state that is flavor specific; i.e., only accessible from either a B 0 or a " B 0 , but not from both. The purpose of the flavor tagging algorithm is to determine the flavor of B tag with the highest possible efficiency and lowest possible probability w of assigning a wrong flavor to B tag . The figure of merit for the performance of the tagging algorithm is the effective tagging efficiency which is approximately related to the statistical uncertainty ' in the coefficients S and C through It is not necessary to reconstruct B tag fully to determine its flavor. We use a neural network based technique [22] to exploit signatures of B decays that determine the flavor at decay of the B tag . Primary leptons from semileptonic B decays are selected from identified electrons and muons as well as isolated energetic tracks. The charges of identified kaon candidates define a kaon tag. Soft pions from D Ãþ decays are selected on the basis of their momentum and direction with respect to the thrust axis of B tag . Based on the output of this algorithm, candidates are divided into seven mutually exclusive categories. These are (in order of decreasing signal purity) Lepton, Kaon I, Kaon II, Kaon-Pion, Pion, Other, and Untagged. We apply this algorithm to a sample of fully reconstructed, self-tagging, neutral B decays (B flav sample). We use B decays to D ðÃÞÀ ð% þ ; & þ ; a þ 1 Þ to measure the tagging efficiency , mistag rate w, and the difference in mistag rates for B 0 and " B 0 tagside decays Áw wðB 0 Þ À wð " B 0 Þ. The results are shown in Table I. The Untagged category of events contains no flavor information and therefore carries no weight in the time-dependent analysis. The total effective tagging efficiency Q for this algorithm is measured to be ð31:2 AE 0:3Þ%.
Including the effects of the mistag rate, Eq. (1) becomes Finally, to account for experimental Át resolution, we convolve Eq. (4) with a resolution function, the parameters of which we obtain from fits to the B flav sample. The Át resolution function is represented as a sum of three Gaussian distributions with different widths. For the core and tail Gaussians, the widths are scaled by ' Át . In addition we allow an offset for the core distribution in the hadronic tagging categories (Sec. IV) separate from that of the Lepton category, to allow for a small bias of Át from secondary D-meson decays; a common offset is used for the tail component. The third Gaussian (of fixed 8 ps width) accounts for the few events with incorrectly reconstructed vertices. Identical resolution function parameters are used for all B CP modes, since the B tag vertex precision dominates the Át resolution. Events without reliable Át information (class b) are sensitive to the parameter C f and are used to constrain this parameter in the % 0 K 0 S analysis. Integrating Eq. (4) over Át we get We also account for the asymmetry in tagging efficiency for B 0 and " B 0 decays, but, for simplicity, we assume the asymmetry is zero in the above equations.

V. EVENT RECONSTRUCTION AND SELECTION
We choose event selection criteria with the aid of a detailed Monte Carlo (MC) simulation of the B production and decay sequences, and of the detector response [23]. These criteria are designed to retain signal events with high efficiency while removing most of the background.
We reconstruct the B CP candidate by combining the four-momenta of the two daughter mesons, with a vertex constraint. The B-daughter candidates are reconstructed with the following decays: In the 0 K 0 S analysis we also reconstruct K 0 S via its decay to two neutral pions (K 0 % 0 % 0 ). The requirements on the invariant masses of these particle combinations are given in Table II. We consider as photons energy depositions in the EMC that are isolated from any charged tracks, carry a minimum TABLE I. Efficiencies ", average mistag fractions w, mistag fraction differences Áw wðB 0 Þ À wð " B 0 Þ, and effective tagging efficiency Q ð1 À 2wÞ 2 for each tagging category from the B flav data. Category L , and !K 0 % 0 % 0 preclude these modes from improving the precision of the measurement of CP parameters in 0 K 0 and !K 0 ; % 0 K 0 L and % 0 K 0 % 0 % 0 events lack the minimum information for reconstruction of the decay vertex.
For decays with a K 0 S ! % þ % À candidate we perform a fit of the entire decay tree which constrains the K 0 S flight direction to the pion pair momentum direction and the K 0 S production point to the B CP vertex determined as described in Sec. III. In this vertex fit we also constrain the , 0 , and % 0 candidate masses to world-average values [24], since these resonances have natural widths that are negligible compared to the resolution. Given that the natural widths of the ! and & mesons are comparable to or greater than the detector resolution, we do not impose any constraint on the masses of these candidates; constraining the mass of the K 0 S does not improve determination of the vertex. In the !K 0 S and 0 K 0 S analyses, we require the 1 2 probability of this fit to be greater than 0.001. We also require that the K 0 S flight length divided by its uncertainty be greater than 3.0 (5.0 for L candidates are reconstructed from clusters of energy deposited in the EMC or from hits in the IFR not associated with any charged track in the event [25]. Taking the flight direction of the K 0 L to be the direction from the B 0 decay vertex to the cluster centroid, we determine the K 0 L momentum p K 0 L from a fit with the B 0 and K 0 L masses constrained to world-average values [24].
In this experiment the energy of the initial e þ e À state is known within an uncertainty of a few MeV. For a final state with two particles we can determine four kinematic variables from conservation of energy and momentum. These may be taken as polar and azimuthal angles of the line of flight of the two particles, and two energy, momentum, or mass variables, such as the masses of the two particles. In practice, since we fully reconstruct one B meson candidate, we make the assumption that it is one of two final-state particles of equal mass. We compute two largely uncorrelated variables that test consistency with this assumption, and with the known value [24] of the B-meson mass. The choice of these variables depends on the decay process, as we discuss below.
In the reconstruction of B 0 ! !K 0 S and B 0 ! 0 K 0 S the kinematic variables are the energy-substituted mass and the energy difference where ðE 0 ; p 0 Þ and ðE B ; p B Þ are the laboratory fourmomenta of the Çð4SÞ and the B CP candidate, respectively, and the asterisk denotes the Çð4SÞ rest frame. The resolution is 3 MeV in m ES and 20-50 MeV in ÁE, depending on the decay mode. We require 5:25 < m ES < 5:29 GeV and jÁEj < 0:2 GeV, as distributions of these quantities for signal events peak at the B-meson mass in m ES and zero in ÁE.
For the B 0 ! 0 K 0 L channel only the direction of the K 0 L momentum is measured. For these candidates m ES is not determined; instead we obtain ÁE from a calculation with the B 0 and K 0 L masses constrained to world-average values. Because of the mass constraint on the B 0 , the ÁE distribution for K 0 L events, which peaks at zero for signal, is asymmetric and narrower than that of K 0 S events; we require À0:01 < ÁE < 0:08 GeV for 0 K 0 L . For the % 0 K 0 S analysis we use the kinematic variables m B and m miss . The variable m B is the invariant mass of the reconstructed B CP . The variable m miss is the invariant mass of the B tag , computed from the known beam energy and the measured B CP momentum with mðB CP Þ constrained to the nominal B-meson mass m PDG B [26]. For signal decays, m B and m miss peak at the B 0 mass and have resolutions of $47 MeV and $5:4 MeV, respectively; the distribution of m B exhibits a lowside tail due to leakage of energy out of the EMC. To compare the m miss resolution with the m ES resolution a factor of 2 from the approximate relation m ES $ ðm miss þ m PDG B Þ=2 should be taken into account. The beam-energy constraint in m miss helps to eliminate the correlation with m B . We select candidates within the window 5:11 < m miss < 5:31 GeV and 5:13 < m B < 5:43 GeV, which includes the signal peak and a sideband region for background characterization.

B. Background reduction
Background events arise primarily from random combinations of particles in continuum e þ e À ! q " q events (q ¼ u, d, s, c). For some of the decay chains we must also consider cross feed from B-meson decays by modes other than the signal; we discuss these in Secs. VI A and VI B below. To reduce the q " q backgrounds we make use of additional properties of the event that are consequences of the decay.
For the B 0 ! !K 0 S and B 0 ! 0 K 0 channels we define the angle T between the thrust axis [27] of the B CP candidate in the Çð4SÞ frame and that of the charged tracks and neutral calorimeter clusters in the rest of the event. The event is required to contain at least one charged track not associated with the B CP candidate. The distribution of j cos T j is sharply peaked near 1 for q " q jet pairs, and nearly uniform for B-meson decays. The requirement is j cos T j < 0:9 for all modes.
For the 0 & decays we also define the angle & dec between the momenta of the & 0 daughter % À and of the 0 , measured in the & 0 rest frame. We require j cos & dec j < 0:9 to reduce the combinatorial background of & 0 candidates incorporating a soft pion that are reconstructed as decays with j cos & dec j ' 1. For 0 K 0 L candidates we require that the cosine of the polar angle of the total missing momentum in the laboratory system be less than 0.96 to reject very forward q " q jets. We construct the missing momentum p miss as the difference of p 0 and the momenta of all charged tracks and neutral clusters not associated with the K 0 L candidate. We project p miss onto p K 0 L , and require the component perpendicular to the beam line, p proj miss? , to satisfy p proj miss? À p K 0 L ? > À0:8 GeV. These values are chosen to minimize the uncertainty on S and C in the presence of background.
The purity of the K 0 L candidates reconstructed in the EMC is further improved by a requirement on the output of a neural network (NN) that takes cluster-shape variables as inputs. For the NN, we use the following eight variables: the number of crystals in the EMC cluster; the total energy deposited in the EMC cluster; the second moment of the cluster energy where E i is the energy deposited in the i th crystal and r i its distance from the cluster centroid; the lateral moment where E 0 refers to the most energetic crystal and E n to the least energetic one; the ratio S1=S9 of the energy of the most energetic crystal (S1) to the sum of energy of the 3 Â 3 crystal block with S1 in its center (S9); the ratio S9=S25, where S25 is the sum of energy of the 5 Â 5 crystal block with S1 in its center; and the absolute value of the expansion coefficients jZ 20 j and jZ 42 j of the spatial energy distribution of the EMC cluster expressed as a series of Zernike polynomials ðÞ : Eðx; yÞ ¼ P n;m Z n;m Á n;m ðr; 0Þ, where ðx; yÞ are the Cartesian coordinates in the plane of the calorimeter, ðr; 0Þ are the polar coordinates of the Zernike polynomials (0 r 1) and n, m are nonnegative integers. The NN is trained on MC signal events and off-peak data reconstructed for the 0 ðÞ %% K 0 L decay mode to return þ1 if the event is signal-like and À1 if it is background-like. We check the performance of the NN on an independent sample of MC signal events and off-peak data reconstructed for the 0 ðÞ %% K 0 L decay mode. Using ensembles of simulated experiments, as discussed in Sec. VI D, we find that requiring the output of the NN to be greater than À0:2 minimizes the average statistical uncertainty on S and C.
For the % 0 K 0 S channel we require the 1 2 probability of the kinematic fit to be greater than 0.001. We exclude events in which the absolute value of cosine of the angle between the beam axis and the B CP momentum in the Çð4SÞ frame ( cos B ) is greater than 0.9. Finally, we apply a cut on the event shape, selecting events with L 2 =L 0 < 0:55; L i is the i th angular moment defined as L i ¼ P j p j Â j cos j j i , where j is the angle with respect to the B thrust axis of track or neutral cluster j, p j is its momentum, and the sum excludes the daughters of the B candidate.
The average number of candidates found per selected event in 0 K 0 and !K 0 S is between 1.08 and 1.32, depending on the final state. For events with multiple candidates we choose the one with the largest decay vertex probability for the B meson. Furthermore, in the B 0 ! 0 K 0 L sample, if several B candidates have the same vertex probability, we choose the candidate with the K 0 L information taken from, in order, EMC and IFR, EMC only, or IFR only. From the simulation we find that this algorithm selects the correctcombination candidate in about two-thirds of the events containing multiple candidates.
For the % 0 K 0 S channel the average number of candidates found per selected event is 1.03. In events with multiple candidates we choose the one with the smallest value of 1 2 obtained from the reconstructed mass of the K 0 S and the % 0 candidates and their respective errors.
In the 0 K 0 analysis, we estimate from MC the fraction of events in which we misreconstruct charged daughters of the 0 (self-crossfeed events), which dominate the determination of the B CP decay vertex. We find that (1-4)% of events are self-crossfeed, depending on the 0 and K 0 decay channel.

VI. MAXIMUM LIKELIHOOD FIT
The selected sample sizes are given in the table of results in Sec. VII. We perform an unbinned maximum likelihood (ML) fit to these data to obtain the common CP-violation parameters and signal yields for each channel. For each signal or background component j, tagging category c, and event class t ¼ g (Sec. III), we define a total probability density function (PDF) for event i as where T is the function FðÁtÞ defined in Eq. (4) convolved with the Át resolution function, and ' ¼ AE1 is the B flavor defined by upper and lower signs in Eq. (1). For event class t ¼ b (for the % 0 K 0 S analysis), we define a total PDF for event i as where F C is the function defined in Eq. (5). The set of variables x i k , which serve to discriminate signal from background, are described along with their PDFs Q k;j ðx k Þ in Sec. VI B below. The factored form of the PDF is a good approximation since linear correlations are smaller than 5%, 7%, and 6% in the !K 0 S , 0 K 0 , and % 0 K 0 S analyses, respectively. The effects of these correlations are estimated as described in Sec. VI D.
We write the extended likelihood function for all candidates for decay mode d as where n j is the yield of events of component j, f j;t is the fraction of events of component j for each event class t, j;c is the efficiency of component j for each tagging category c, and N c;t is the number of events of category c with event class t in the sample. When combining decay modes we form the grand likelihood L ¼ Q L d .
We fix j;c for all components except the q " q background to B flav ;c , which are listed in Table I. For the % 0 K 0 S channel we assume the same j;c for class g and class b events. For !K 0 S and 0 K 0 we fix f j;g ¼ 1 because we accept only events of class g.

A. Model components
For all of the decay chains we include in the ML fit a component for q " q combinatorial background (j ¼ q " q), in addition to the one for the signal. The functional forms of the PDFs that describe this background are determined from fits of one observable at a time to sidebands of the data in the kinematic variables that exclude signal events. Some of the parameters of this PDF are free in the final fit. Thus, the combinatorial component receives contributions from all nonsignal events in the data.
We estimate from the simulation that charmless B decay modes contribute less than 2% of background to the input sample. These events have final states different from the signal, but similar kinematics, and exhibit broad peaks in the signal regions of some observables. We find that the charmless B " B background component (j ¼ chls) is needed for the final states !K 0 S , 0 & K 0 % þ % À , and 0 & K 0 % 0 % 0 . We account for these with a separate component in the PDF. Unlike the other fit components, we fix the charmless B " B yields using measured branching fractions, where available, and detection efficiencies determined from MC. For unmeasured background modes, we use theoretical estimates of branching fractions.
We also consider the presence of B decays to charmed particles in the input sample. The charmed hadrons in these final states tend to be too heavy to be misreconstructed as the two light bodies contained in our signals, and their distributions in the B kinematic variables are similar to those for q " q. However, in the event shape variables and Át they are signal-like. We have found that biases in the fit results are minimized for the modes with 0 & by including a component specifically for the B decays to charm states (j ¼ chrm). Finally, for 0 ð3%Þ%% K 0 L we divide the signal component into two categories for correctly reconstructed and self-crossfeed events; we fix the fraction of the selfcrossfeed category to the value obtained from MC.

B. Probability density functions
The set of variables x k of Eqs. (10) and (11) is defined for each family of decays as (i) B 0 ! !K 0 S : fm ES ;ÁE;F ;mð% þ % À % 0 Þ m ! ;H g, (ii) B 0 ! 0 K 0 S : fm ES ; ÁE; F g, (iii) B 0 ! 0 K 0 L : fÁE; F g, (iv) B 0 ! % 0 K 0 S : fm B ; m miss ; L 2 =L 0 ; cos B g. Here, F is a Fisher discriminant described below, and H is the cosine of the polar angle of the normal to the ! decay plane in the ! helicity frame, which is defined as the ! rest frame with polar axis opposite to the direction of the B. From Monte Carlo studies we find that including the 0 mass, & mass, & helicity, or ! Dalitz plot coordinates does not improve the precision of the measurements of S and C.
In Fig. 1 we show PDFs for the signal and q " q components for the !K 0 S analysis, which are similar to those for the 0 K 0 S analysis. We parameterize the PDFs for the signal component using simulated events, while the background distributions are taken from sidebands of the data in the kinematic variables that exclude signal events. The parameters used in the PDFs are different for each mode.
For the background PDF shapes we use the following: the sum of two Gaussians for Q sig ðm ES Þ and Q sig ðÁEÞ; a quadratic dependence for Q q " q ðÁEÞ, Q chrm ðÁEÞ, and Q q " q ðm B Þ; and the sum of two Gaussians for Q chls ðÁEÞ. For Q q " q ðm ES Þ and Q q " q ðm miss Þ we use the function with x 2m ES = ffiffi ffi s p ð2m miss = ffiffi ffi s p for % 0 K 0 S Þ and $ a free parameter [28], and the same function plus a Gaussian for Q chrm ðm ES Þ and Q chls ðm ES Þ. For Q sig ðm B Þ and Q sig ðm miss Þ we use the function where " is the peak position of the distribution, ' L;R are the left and right widths, and L;R are the left and right tail parameters. For Q q " q ðÁEÞ in the 0 K 0 L analysis, we use the function where x ÁE À ðÁEÞ min , with ðÁEÞ min fixed to À0:01, and $ 0 is a free parameter. To reduce q " q background beyond that obtained with the cos T requirement described above for !K 0 S and 0 K 0 (and the B and L 2 =L 0 requirements for % 0 K 0 S ), we use additional event topology information in the ML fit. The variables used include B , L 0 , L 2 , and the angle with respect to the beam axis in the Çð4SÞ frame of the signal B thrust axis ð S Þ. For the % 0 K 0 S analysis, we use cos B and the ratio L 2 =L 0 directly in the fit, parameterized by a second-order polynomial and a seven-bin histogram, respectively. The parameters of the L 2 =L 0 PDF depend on the tagging category c in the signal component. In Fig. 2   distribution for data where background (signal) events are subtracted using an event weighting technique [29]. The bin widths of the L 2 =L 0 histogram have been adjusted to be coarser where the background is small to reduce the number of free parameters of the PDF.
For the other decay modes we construct a Fisher discriminant F , which is an optimized linear combination of L 0 , L 2 , j cos B j, and j cos S j. For the K 0 L modes we also use the continuous output of the flavor tagging algorithm as a variable entering the Fisher discriminant. The coefficients used to combine these variables are chosen to maximize the separation (difference of means divided by quadrature sum of errors) between the signal and continuum background distributions of F , and are determined from studies of signal MC and off-peak data. We have studied the optimization of F for a variety of signal modes, and find that a single set of coefficients is nearly optimal for all.
The PDF shape for F is an asymmetric Gaussian with different widths below and above the peak for signal, plus a broad Gaussian for q " q background to account for a small tail in the signal F region. The background peak parameter is adjusted to be the same for all tagging categories c. Because F describes the overall shape of the event, the distribution for B " B background is similar to the signal distribution.
For Q sig ðm ! Þ we use the sum of two Gaussians; for Q q " q ðm ! Þ and Q chls ðm ! Þ the sum of a Gaussian and a quadratic. For Q sig ðH Þ and Q q " q ðH Þ we use a quadratic dependence, and for Q chls ðH Þ a fourth-order polynomial. As described in Sec. IV, the resolution function in T j ðÁtÞ is a sum of three Gaussians for all fit components j. For q " q background we use the same functional form T j ðÁtÞ as for signal, but fix the B lifetime ( to zero so that T j ðÁtÞ is effectively just the resolution model. For the signal and B " B background components we determine the parameters of Q k;j ðx i k Þ from simulation, and the q " q background parameters are free in the final fit. For the signal resolution function we fix all parameters to values obtained from the B flav sample; we obtain parameter values from MC for the charm and charmless B " B resolution models; we leave parameters of the Át resolution model for q " q free in the final fit. For the !K 0 S and 0 K 0 analyses, we use large control samples to determine any needed adjustments to the signal PDF shapes that have initially been determined from Monte Carlo. For m ES and ÁE in !K 0 S and 0 K 0 S , we use the decay B À ! % À D 0 with D 0 ! K À % þ % 0 , which has similar topology to the modes under study here. We select this sample by making loose requirements on m ES and ÁE, and requiring for the D 0 candidate mass 1845 < m D < 1885 MeV. We also place kinematic requirements on the D and B daughters to force the charmed decay to look as much like that of a charmless decay as possible. These selection criteria are applied both to the data and to MC. For F , we use a sample of B þ ! 0 & K þ decays selected with requirements very similar to those of our signal modes. From these control samples, we determine small adjustments (of the order of few MeV) to the mean value of the signal ÁE distribution. The means and widths of the other distributions do not need adjustment. For the ! mass line shape, we use ! production in the data sidebands. The means and resolutions of the invariant mass distributions are compared between data and MC, and small adjustments are made to the PDF parameterizations. These studies also provide uncertainties in the agreement between data and MC that are used for evaluation of systematic errors. For the % 0 K 0 S analysis, we apply no correction to the signal PDF shapes, but we evaluate the related systematic error as described in Sec. VIII.

C. Fit variables
For the !K 0 S analysis we perform a fit with 25 free parameters: S, C, signal yield, continuum background yield and fractions (6), and the background PDF parameters for Át, m ES , ÁE, F , m ! , and H (15). For the five 0 K 0 S channels we perform a single fit with 98 free parameters: S, C, signal yields (5), 0 & K 0 S charm B " B background yields (2), continuum background yields (5) and fractions (30), and the background PDF parameters for Át, m ES , ÁE, and F (54). Similarly, the two 0 K 0 L channels are fit jointly, with 34 parameters: S, C, signal yields (2), background yields (2), fractions (12), and PDF parameters (16). For the % 0 K 0 S analysis we perform a fit with 36 free parameters: S, C, signal yield, the means of m miss and m B signal PDFs, background yield, background PDF parameters for Át, m B , m miss , L 2 =L 0 , cos B (16), background tagging efficiencies (12), and the fraction of good events (2). For the signal, charm B " B, and charmless B " B components the parameters ( and Ám d are fixed to world-average values [24]; for the B " B components S and C are fixed to zero and then varied to obtain the related systematic uncertainty as described below; for the q " q model ( is fixed to zero.

D. Fit validation
We test the fitting procedure by applying it to ensembles of simulated experiments with q " q and B " B charmed events drawn from the PDF into which we have embedded the expected number of signal and B " B charmless background events (with the expected values of S and C) randomly extracted from the fully simulated MC samples. We find biases (measured À expected) for S !K 0 S , S 0 K 0 S , C 0 K 0 S , S 0 K 0 L , and C 0 K 0 L of 0:034 AE 0:010, 0:006 AE 0:006, À0:008 AE 0:005, À0:022 AE 0:014, and À0:013 AE 0:009, respectively. These small biases are due to neglected correlations among the observables, contamination of the signal by self-crossfeed, and the small signal event yield in !K 0 S . We apply additive corrections to the final results for these biases. For C !K 0 S , S % 0 K 0 S , and C % 0 K 0 S we make no correction but assign a systematic uncertainty as described in Sec. VIII.

VII. FIT RESULTS
Results from the fits for the signal yields and the CP parameters S f and C f are presented in Table III. In Figs. 3-9, we show projections onto the kinematic variables and Át for subsets of the data for which the ratio of the likelihood to be signal and the sum of likelihoods to be signal and background (computed without the variable plotted) exceeds a mode-dependent threshold that optimizes the statistical significance of the plotted signal. In !K 0 S the fraction of signal events with respect to the total after this requirement has been applied is $70%, while in 0 K 0 S and 0 K 0 L , the fraction of signal events is in the (42-85)% and (22-55)% range, respectively, depending on the decay mode. In Fig. 3     background-subtracted distributions for m B and m miss in Fig. 2.
In Figs. 6-9, we show the Át projections and the asymmetry ðN B 0 À N " B 0 Þ=ðN B 0 þ N " B 0 Þ for each final state. In the !K 0 S , 0 K 0 S , 0 K 0 L , and % 0 K 0 S analyses, we measure the correlation between S f and C f in the fit to be 2.9%, 3.0%, 1.0%, and À6:2%, respectively.

Crosschecks
We perform several additional crosschecks of our analysis technique including time-dependent fits for B AE decays to the final states 0 ðÞ%% K AE , 0 & K AE , and 0 ð3%Þ%% K AE in which measurements of S and C are consistent with zero. There are only small changes in the results when we do any of the following: fix C ¼ 0 or S ¼ 0, allow S and C to be different for each tagging category, remove each of the discriminating variables one by one, and allow the signal resolution model parameters to vary in the fit.
To validate the IP-constrained vertexing technique in % 0 K 0 S , we examine B 0 ! J=c K 0 S decays in data where J=c ! " þ " À or J=c ! e þ e À . In these events we determine Át in two ways: by fully reconstructing the B 0 decay vertex using the trajectories of charged daughters of the J=c and the K 0 S mesons (standard method), or by neglecting the J=c contribution to the decay vertex and using the IP constraint and the K 0 S trajectory only. This study shows that within statistical uncertainties of order 2% of the error on Át, the IP-constrained Át measurement is unbiased with respect to the standard technique and that the fit values of S J=c K 0 S and C J=c K 0 S are consistent between the two methods.

VIII. SYSTEMATIC UNCERTAINTIES
A number of sources of systematic uncertainties affecting the fit values of S f and C f have been considered.
We vary the parameters of the signal PDFs that are kept fixed in the fit within their uncertainty and take as systematic error the resulting changes of S f and C f . These parameters include ( and Ám d , the mistag parameters w and Áw, the efficiencies of each tagging category, the parameters of the resolution model, and the shift and scale factors applied to the variables related to the B kinematics and event shape variables that serve to distinguish signal from background. The deviations of S f and C f for 0 K 0 S and 0 K 0 L for variations of ( and Ám d are less than 0.002. For the % 0 K 0 S channel as an additional systematic error associated with the shape of the PDFs we also use the largest deviation observed when the parameters of the individual PDFs are free in the fit. As a systematic uncertainty related to the fit bias on S f and C f we assign the statistical uncertainty on the bias obtained from simulated experiments during the fit validation. As explained in Sec. IV, we obtain parameters of the signal resolution model from a fit to the B flav sample instead of from a fit to signal MC. We evaluate the systematic uncertainty of this approach with two sets of simulated experiments that differ only in the values of resolution model parameters (one set with parameters from the B flav sample and one set with parameters from MC). We take the difference in the average S f and C f from these two sets of experiments as the related systematic error.
We evaluate the impact of potential biases arising from the interference of doubly Cabibbo-suppressed decays with the Cabibbo-favored decays on the tagside of the event [30] by taking into account realistic values of the ratio between the two amplitudes and the relative phases. For !K 0 S and 0 K 0 we estimate using MC, published measurements, and theoretical predictions that conservative ranges of the net values for CP parameters in the B " B background are S ¼ ½0; 0:2 and C ¼ AE0:1 for the charmless background and S ¼ AE0:1 and C ¼ AE0:1 for the charm background. We perform a fit in which we fix the parameters to these values and take the difference in signal CP parameters between this fit and the nominal fit as the systematic error.
For the !K 0 S and 0 K 0 channels we also vary the amount of the charmless B " B background by AE20%. For % 0 K 0 S we do not include a B " B background component in the fit but we embed B " B background events in the data sample and extract the peaking background from the observed change B 0 tags. Points with error bars represent the signal where background is subtracted using an event weighting technique [29]; the solid line displays the signal fit function. In (c) we show the raw asymmetry, in the yield. We use this yield to estimate the change in S and C due to the CP asymmetry of the peaking background. We also measure the systematic error associated with the vertex reconstruction by varying within uncertainties the parameters of the alignment of the SVT and the position and size of the beam-spot. We quantify the effects of self-crossfeed events in the 0 K 0 analysis. For 0 K 0 S we perform sets of simulated experiments in which we embed only correctly reconstructed signal events and compare the results to the nominal simulated experiments (Sec. VI D) in which we embed both correctly and incorrectly reconstructed signal events. We take the difference as the systematic uncertainty related to self-crossfeed. For the 0 ð3%Þ%% K 0 L analysis, in which we include a self-crossfeed component in the fit, we perform a fit in which we take parameter values for the self-crossfeed resolution model from self-crossfeed MC events instead of the nominal B flav sample. We take the difference of the results from this fit and the nominal fit as the self-crossfeed systematic for 0 K 0 L . The effects of self-crossfeed are negligible for !K 0 S and % 0 K 0 S . Finally, for the % 0 K 0 S analysis we examine large samples of simulated B 0 ! K 0 S % 0 and B 0 ! J=c K 0 S decays to quantify the differences between resolution function parameter values obtained from the B flav sample and those of the signal channel; we use these differences to evaluate uncertainties due to the use of the resolution function extracted from the B flav sample. We also use the differences between resolution function parameters extracted from data and MC in the B 0 ! J=c K 0 S decays to quantify possible problems in the reconstruction of the K 0 S vertex. We take the sum in quadrature of these errors as the systematic error related to the vertexing method.
The contributions of the above sources of systematic uncertainties to S f and C f are summarized in Table IV. IX. S AND C PARAMETERS FOR B 0 ! 0 K 0 As noted in Sec. I, the final states 0 K 0 S and 0 K 0 L have opposite CP eigenvalues, and in the SM, if ÁS f ¼ 0, then À f S f ¼ sin2. We therefore compute the values of S 0 K 0 and C 0 K 0 from our separate measurements with B 0 ! 0 K 0 S and B 0 ! 0 K 0 L , taking ÀS 0 K 0 L in combination with S 0 K 0 S , and C 0 K 0 L with C 0 K 0 S . To represent the results of the individual fits, we project the likelihood by maximizing L (Sec. VI) at a succession of fixed values of S f to obtain LðÀ f S f Þ. We then convolve this likelihood with a Gaussian function representing the independent systematic errors for each mode. The product of these convolved one-dimensional likelihood functions for the two modes, shifted in À f S f by their respective corrections (Sec. VI D), gives the joint likelihood for S 0 K 0 . The likelihood for C 0 K 0 is computed similarly. Since the measured correlation between S f and C f is small in our fits (Sec. VII), we extract the central values and total uncertainties of these quantities from these one-dimensional likelihood functions. Applying the same procedure without the convolution over systematic errors yields the statistical component of the error. The systematic component is then extracted by subtraction in quadrature from the total error.

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 U.S. 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.