Measurement of Branching Fractions and CP-Violating Charge Asymmetries for B Meson Decays to D(*)D(*), and Implications for the CKM Angle gamma

We present measurements of the branching fractions and charge asymmetries of B decays to all D(*)D(*) modes. Using 232 million BBbar pairs recorded on the Upsilon(4S) resonance by the BaBar detector at the e+e- asymmetric B factory PEP-II at the Stanford Linear Accelerator Center, we measure the branching fractions BF(B0 ->D*+D*-) = (8.1 +- 0.6 +- 1.0) x 10^-4, BF(B0 ->D*+-D-+) = (5.7 +- 0.7 +- 0.7) x 10^-4, BF(B0 ->D+D-) = (2.8 +- 0.4 +- 0.5) x 10^-4, BF(B+ ->D*+D*0) = (8.1 +- 1.2 +- 1.2) x 10^-4, BF(B+ ->D*+D0) = (3.6 +- 0.5 +- 0.4) x 10^-4, BF(B+ ->D+D*0) = (6.3 +- 1.4 +- 1.0) x 10^-4, and BF(B+ ->D+D0) = (3.8 +- 0.6 +- 0.5) x 10^-4, where in each case the first uncertainty is statistical and the second systematic. We also determine the limits BF(B0 ->D*0D*0)<0.9 x 10^-4, BF(B0 ->D*0D0)<2.9 x 10^-4, and BF(B0 ->D0D0)<0.6 x 10^-4, each at 90% confidence level. All decays above denote either member of a charge conjugate pair. We also determine the CP-violating charge asymmetries A(B0 ->D*+-D-+) = 0.03 +- 0.10 +- 0.02, A(B+ ->D*+D*0) = -0.15 +- 0.11 +- 0.02, A(B+ ->D*+D0) = -0.06 +- 0.13 +- 0.02, A(B+ ->D*+D*0) = 0.13 +- 0.18 +- 0.04, and A(B+ ->D+D0) = -0.13 +- 0.14 +- 0.02. Additionally, when we combine these results with information from time-dependent CP asymmetries in B0 ->D*+D*- decays and world-averaged branching fractions of B decays to Ds(*)D(*) modes, we find the CKM phase gamma is favored to lie in the range [0.07-2.77] radians (with a +0 or +pi radians ambiguity) at 68% confidence level.

Additionally, when we combine these results with information from time-dependent CP asymmetries in B 0 → D ( * )+ D ( * )− decays and world-averaged branching fractions of B decays to D ( * ) s D ( * ) modes, we find the CKM phase γ is favored to lie in the range [0.07 − 2.77] radians (with a +0 or +π radians ambiguity) at 68% confidence level.

I. INTRODUCTION
We report on measurements of branching fractions of neutral and charged B-meson decays to the ten doublecharm final states D ( * ) D ( * ) . For the four charged B decays to D ( * ) D ( * ) and for neutral B decays to D * ± D ∓ , we also measure the direct CP -violating time-integrated charge asymmetry where in the case of the charged B decays, the superscript on Γ corresponds to the sign of the B ± meson, and for D * ± D ∓ , Γ + refers to D * − D + and Γ − to D * + D − . In the neutral B → D ( * )+ D ( * )− decays, the interference of the dominant tree diagram (see Fig. 1a) with the neutral B mixing diagram is sensitive to the Cabibbo-Kobayashi-Maskawa (CKM) phase β ≡ arg [ −V cd V * cb /V td V * tb ], where V is the CKM quark mixing matrix [1]. However, the theoretically uncertain contributions of penguin diagrams (Fig. 1b) with different weak phases are potentially significant and may shift both the * Also at Laboratoire de Physique Corpusculaire, Clermont-Ferrand, France † Also with Università di Perugia, Dipartimento di Fisica, Perugia, Italy ‡ Also with Università della Basilicata, Potenza, Italy observed CP asymmetries and the branching fractions by amounts that depend on the ratios of the penguin to tree contributions and their relative phases. A number of theoretical estimates exist for the resulting values of the branching fractions and CP asymmetries [2][3][4][5][6].
The penguin-tree interference in neutral and charged B → D ( * ) D ( * ) decays can provide sensitivity to the angle γ = arg [ −V ud V * ub /V cd V * cb ] [7,8]. With additional information on the branching fractions of B → D ( * ) s D ( * ) decays, the weak phase may be extracted, assuming SU(3) flavor symmetry between B → D ( * ) D ( * ) and B → D ( * ) s D ( * ) . For this analysis, we assume that the breaking of SU(3) can be parametrized via the ratios of decay constants f D ( * ) s /f D ( * ) , which are quantities that can be determined either with lattice QCD or from experimental measurements [9].
In addition to presenting measurements of the B 0 → D ( * )+ D ( * )− and B + → D ( * )+ D ( * )0 branching fractions, and the CP -violating charge asymmetries for the latter modes and for B 0 → D * ± D ∓ , we search for the colorsuppressed decay modes B 0 → D ( * )0 D ( * )0 , which have not been previously measured, and determine limits on those branching fractions [10]. If observed, the decays B 0 → D ( * )0 D ( * )0 would provide evidence of W -exchange or annihilation contributions (see Fig. 1c,1d). In principle, these decays could also provide sensitivity to the CKM phase β if sufficient data were available. By combining all of these results with information from time-dependent CP asymmetries in s D ( * ) modes, we determine the implications for γ using the method of Refs. [7,8].

II. DETECTOR AND DATA
The results presented in this paper are based on data collected with the BABAR detector [11] at the PEP-II asymmetric-energy e + e − collider [12] located at the Stanford Linear Accelerator Center. The integrated luminosity is 210.5 fb −1 , corresponding to 231.7 million BB pairs, recorded at the Υ (4S) resonance ("on-peak", at a center-of-mass (c.m.) energy √ s = 10.58 GeV). The asymmetric beam configuration in the laboratory frame provides a boost of βγ = 0.56 to the Υ (4S). Charged particles are detected and their momenta measured by the combination of a silicon vertex tracker (SVT), consisting of five layers of double-sided detectors, and a 40-layer central drift chamber (DCH), both operating in the 1.5-T magnetic field of a solenoid. For tracks with transverse momentum greater than 120 MeV/c, the DCH provides the primary charged track finding capability. The SVT provides complementary standalone track finding for tracks of lower momentum, allowing for reconstruction of charged tracks with transverse momentum p T as low as 60 MeV/c, with efficiencies in excess of 85%. This ability to reconstruct tracks with low p T efficiently is necessary for reconstruction of the slow charged pions from D * + → D 0 π + decays in B → D ( * ) D ( * ) signal events. The transverse momentum resolution for the combined tracking system is σ pT /p T = 0.0013p T +0.0045, where p T is measured in GeV/c. Photons are detected and their energies measured by a CsI(Tl) electromagnetic calorimeter (EMC). The photon energy resolution is σ E /E = 2.3/E( GeV) 1/4 ⊕ 1.4 %, and their angular resolution with respect to the interaction point is σ θ = (4.2 mrad)/ E( GeV). The measured π 0 mass resolution for π 0 's with laboratory momentum in excess of 1 GeV/c is approximately 6 MeV/c 2 .
Charged-particle identification (PID) is provided by an internally reflecting ring-imaging Cherenkov light detector (DIRC) covering the central region, and the most probable energy loss (dE/dx) in the tracking devices. The Cherenkov angle resolution of the DIRC is measured to be 2.4 mrad, which provides over 5σ separation between charged kaons and pions at momenta of less than 2 GeV/c. The dE/dx resolution from the drift chamber is typically about 7.5% for pions. Additional information to identify and reject electrons and muons is provided by the EMC and detectors embedded between the steel plates of the magnetic flux return (IFR).

III. CANDIDATE RECONSTRUCTION AND B MESON SELECTION
Given the high multiplicity of the final states studied, very high combinatorial background levels are expected. Selection criteria (described in Sec. III A-E) are designed to minimize the expected statistical error on the B branching fractions (as described in Sec. III F). A GEANT4-based [13] Monte Carlo (MC) simulation of the material composition and the instrumentation response of the BABAR detector is used to optimize signal selection criteria and evaluate signal detection efficiency. We retain sufficient sidebands in the discriminating variables to characterize the background in subsequent fits.
A. Charged track and K 0 S selection Charged particle tracks are selected via pattern recognition algorithms using measurements from the SVT and DCH detectors. We additionally require all chargedparticle tracks (except for those from K 0 S → π + π − decays) to originate within 10 cm along the beam axis and 1.5 cm in the plane perpendicular to the beam axis of the center of the beam crossing region. To ensure a well-measured momentum, all charged-particle tracks except those from K 0 S → π + π − decays and π + from D * + → D 0 π + decays must also be reconstructed from at least 12 measurements in the DCH. All tracks that meet these criteria are considered as charged pion candidates.
Tracks may be identified as kaons based on a likelihood selection developed from Cherenkov angle and dE/dx information from the DIRC and tracking detectors respectively. For the typical laboratory momentum spectrum of the signal kaons, this selection has an efficiency of about 85% and a purity of greater than 98%, as determined from control samples of D * + → D 0 π + , D 0 → K − π + decays.
We require K 0 S → π + π − candidates to have an invariant mass within 15 MeV/c 2 of the nominal K 0 S mass [14]. The probability that the two daughter tracks originate from the same point in space must be greater than 0.1%. The transverse flight distance of the K 0 S from the primary event vertex must be both greater than 3σ from zero (where σ is the measured uncertainty on the transverse flight length) and also greater than 2 mm.

B. Photon and π 0 selection
Photons are reconstructed from energy deposits in the electromagnetic calorimeter which are not associated with a charged track. To reject backgrounds from electronics noise, machine background, and hadronic interactions in the EMC, we require that all photon candidates have an energy greater than 30 MeV in the laboratory frame and to have a lateral shower shape consistent with that of a photon. Neutral pions are reconstructed from pairs of photon candidates whose energies in the laboratory frame sum to more than 200 MeV. The π 0 candidates must have an invariant mass between 115 and 150 MeV/c 2 . The π 0 candidates that meet these criteria, when combined with other tracks or neutrals to form B candidates, are then constrained to originate from their expected decay points, and their masses are constrained to the nominal value [14]. This procedure improves the mass and energy resolution of the parent particles.

C. Event selection
We select BB events by applying criteria on the track multiplicity and event topology. At least three reconstructed tracks, each with transverse momentum greater than 100 MeV/c, are required in the laboratory polar angle region 0.41 < θ lab < 2.54. The event must have a total measured energy in the laboratory frame greater than 4.5 GeV to reject beam-related background. The ratio of Fox-Wolfram moments H 2 /H 0 [15] is a parameter between 0 (for "perfectly spherical" events) and 1 (for "perfectly jet-like" events), and we require this ratio to be less than 0.6 for each event, in order to help reject non-BB background. This criterion rejects between 30 and 50 percent of non-BB background (depending on the decay mode), while keeping almost all of the signal decays.

D. D and D * meson selection
We reconstruct D 0 mesons in the four decay modes D 0 → K − π + , D 0 → K − π + π 0 , D 0 → K − π + π − π + , and D 0 → K 0 S π + π − , and D + mesons in the two decay modes D + → K − π + π + and D + → K 0 S π + . We require D 0 and D + candidates to have reconstructed masses within ±20 MeV/c 2 of their nominal masses [14], except for D 0 → K − π + π 0 , for which we require ±40 MeV/c 2 due to the poorer resolution for modes containing π 0 's. These criteria correspond to approximately 2.5σ of the respective mass resolutions. The D 0 → K − π + π 0 decays must also satisfy a criterion on the reconstructed invariant masses of the K − π + and K − π 0 pairs: the combination of reconstructed invariant masses must lie at a point in the K − π + π 0 Dalitz plot [16] for which the expected density normalized to the maximum density ("Dalitz weight") is at least 6%. Additionally, the daughters of D 0 and D + candidates must have a probability of originating from a common point in space greater than 0.1%, and are then constrained both to originate from that common spatial point and to have their respective nominal invariant masses.
Candidate D * + and D * 0 mesons are reconstructed in the decay modes D * + → D 0 π + , D * + → D + π 0 , D * 0 → D 0 π 0 , and D * 0 → D 0 γ, using pairs of selected D 0 , D + , π 0 , π + , and γ candidates. The π + from D * + → D 0 π + decays is additionally required to have a c.m. momentum of less than 450 MeV/c. Candidate π 0 mesons from D * + → D + π 0 and D * 0 → D 0 π 0 are required to have c.m. momenta p * in the range 70 < p * < 450 MeV/c. Photons from D * 0 → D 0 γ decays are required to have energies in the laboratory frame greater than 100 MeV and c.m. energies less than 450 MeV. The D * daughter particles are constrained to originate from a common point in space. After this constraint is applied, the mass differences ∆m of the reconstructed masses of the D * and D candidates are required to be within the ranges shown in Table I. As shown in Fig. 2, the excellent res-  olution in ∆m for signal candidates makes the ∆m requirement a very powerful criterion to reject background (see next section), especially for decay modes containing a D * + → D 0 π + .

E. Variables used for B meson selection
A B-meson candidate is constructed by combining two D ( * ) candidates that have both passed the selection criteria described previously. The pairs of D ( * ) candidates are constrained to originate from the same point in space.
We form a likelihood variable, L Mass , that is defined by a product of Gaussian distributions for each D mass and D * − D mass difference.
For example, in the decay B 0 → D * + D * − , L Mass is the product of four terms: Gaussian distributions for each D mass and double Gaussian (i.e. the sum of two Gaussian distributions) terms for each ∆m term (the D * − D mass difference). Defining G(x; µ, σ) as a normalized Gaussian distribution where x is the independent variable, µ is the mean, and σ is the resolution, L Mass for B 0 → D * + D * − decays is defined as: where the subscript "PDG" refers to the nominal value [14], and all reconstructed masses and uncertainties are determined before mass constraints are applied. For σ mD , we use errors calculated candidate-by-candidate. The parameter f core is the ratio of the area of the core Gaussian to the total area of the double Gaussian distribution. This, along with σ ∆mcore and σ ∆m tail , is determined separately for each of the four D * decay modes given above, using MC simulation of signal events that is calibrated to inclusive samples of the D * decay modes in data. For each of the B decay modes, a higher value of L Mass tends to indicate a greater signal likelihood. The distributions of − ln(L Mass ) for the representative signal mode B 0 → D 0 D 0 and for the corresponding combinatorial background from generic B 0 B 0 , B + B − , cc, and (uu +dd +ss) decays, are shown in Fig. 3a. We use L Mass in selecting signal candidates, as will be described in the upcoming section. We also use the two variables for fully-reconstructed B meson selection at the Υ (4S) energy: the beam-energy- , where the initial total e + e − four-momentum (E i , p i ) and the B momentum p B are defined in the laboratory frame; and ∆E ≡ E cm B − √ s/2 is the difference between the reconstructed B energy in the c.m. frame and its known value. The normalized distribution of ∆E for the representative signal mode B 0 → D 0 D 0 , and for the corresponding combinatorial background components, is shown in Fig. 3b.
In addition to L Mass , m ES , and ∆E, a Fisher discriminant F [17] and a D-meson flight length variable L are used to help separate signal from background. The  D ( * ) decay mode, which are used for the purpose of determining selection criteria that minimize the expected uncertainty on the measured branching fraction for each mode; also, optimized F and L selection criteria for each mode. An "-" indicates no cut is made in F or L for that decay mode.
Fisher discriminant assists in the suppression of background from continuum events by incorporating information from the topology of the event. The discriminant is formed from the momentum flow into nine polar angular intervals of 10 • centered on the thrust axis of the B candidate, the angle of the event thrust axis with respect to the beam axis (θ T ), and the angle of the B candidate momentum with respect to the beam axis (θ B ): The values x i (i = 1, ..., 9) are the scalar sums of the momenta of all charged tracks and neutral showers in the polar angle interval i, x 10 is |cosθ T |, and x 11 is |cosθ B |.
The coefficients α i are determined from MC simulation to maximize the separation between signal and background [17]. The normalized distribution of F for the representative signal mode B 0 → D 0 D 0 , and for the corresponding background components, is shown in Fig. 3c. The flight length variable L that we consider is defined as (ℓ 1 + ℓ 2 )/ σ 2 1 + σ 2 2 , with the decay lengths ℓ i of the two D mesons defined as where x D and x B are the measured decay vertices of the D and B, respectively, and p D is the momentum of a D. The σ i are the measured uncertainties on ℓ i . This observable exploits the ability to distinguish the long D lifetime. Thus, background events have an L distribution centered around zero, while events with real D mesons have a distribution favoring positive values. The normalized distribution of L for the representative signal mode B 0 → D 0 D 0 , and for the corresponding background components, is shown in Fig. 3d.

F. Analysis optimization and signal selection
We combine information from the L Mass , ∆E, F , and L variables to select signal candidates in each decay mode.
The fractional statistical uncertainty on a measured branching fraction is proportional to (N s + N b )/N s , where N s is the number of reconstructed signal events and N b is the number of background events within the selected signal region for a mode. The values N s and N b are calculated, using detailed MC simulation of the signal decay modes as well as of BB and continuum background decays, by observing the number of simulated B decay candidates that satisfy the selection criteria for − ln(L Mass ), |∆E|, F , and L. We choose criteria which minimize the expected (N s + N b )/N s for each mode. Note that to calculate the expected number of signal events N s , one must assume an expected branching fraction, as well as the ratios of BB and continuum events using their relative crosssections. These are given, along with the requirements on F and L, in Table II.
For each possible combination of D * + , D * 0 , D + , and D 0 decay modes, we determine the combination of selection criteria on − ln(L Mass ) and |∆E| that minimizes the overall expected (N s + N b )/N s for each B decay mode (see Tables III, IV, and V). The selection criteria for F and L are chosen, however, only for each B decay mode and not separately for each D ( * ) mode combination. The restrictiveness of the kaon identification selection is also optimized separately for each charged and neutral D ( * ) mode.
Between 1% and 34% of selected B → D ( * ) D ( * ) events have more than one reconstructed B candidate that passes all selection criteria in L Mass , ∆E, F , and L, with the largest percentages occurring in the decay modes B 0 → D * 0 D * 0 and B 0 → D 0 D 0 , and the smallest occurring in B 0 → D * ± D ∓ and B 0 → D + D − . In such events, we choose the reconstructed B with the largest value of L Mass as the signal B candidate.

IV. EFFICIENCY AND CROSSFEED DETERMINATION
The efficiencies are determined using fits to m ES distributions of signal MC events that pass all selection criteria in L Mass , |∆E|, F , and L. There is a small, but non-negligible probability that a signal B decay of mode i is reconstructed as a different signal decay mode j. We refer to this as crossfeed. Thus, efficiencies can be represented as a matrix ǫ ij . where each contributing generated event is weighted by the D and D decay mode branching fractions. To determine the elements of ǫ ij , we fit the m ES distributions of signal MC events generated as B decay mode i and reconstructed as B decay mode j. The distributions are modeled as the sum of signal and background probability distribution functions (PDFs), where the PDF for the signal is a Gaussian distribution centered around the B mass, and the PDF for background is an empirical function [18] of the form where we define x ≡ 2m ES / √ s, and κ is a parameter determined by the fit. In BB MC samples containing signal and background decays, we find that the m ES distribution is well-described by adding a simple Gaussian function to the empirical shape in Eq. 5. We fit the m ES distributions of signal MC events generated as mode i and passing selection criteria in mode j to the above distribution by minimizing the χ 2 ij of each fit with respect to κ ij (the κ parameter for each mode (i, j)), the number of signal events N s ij , and the number of background events N b ij . We determine the efficiencies ǫ ij as N s ij /N g i , where N g i is the total number of signal MC events that were generated in mode i. The diagonal elements of the ǫ ij matrix (i.e. the numbers typically denoted as "efficiencies") are in the range (0.2 -1.5)×10 −3 . The main crossfeed source is misidentification between D * 0 and D * ± candidates. The matrix ǫ ij and the uncertainties on the elements of this matrix are given in Table VI. Crossfeed between different D submodes (i.e. mode numbers [12][13][14][15] in Table III) is negligible.    Table III above. Elements with "-" above and on the diagonal are modes that are unused since, due to high backgrounds, they do not help to increase signal sensitivity.  Table III above. Elements with "-" above and on the diagonal are modes that are unused since, due to high backgrounds, they do not help to increase signal sensitivity.  Mode (3)

V. BRANCHING FRACTION RESULTS
In order to determine the number of signal events in each mode, one must not only account for background which is distributed according to combinatorial phase space, but also for background which can have a different distribution in m ES . It is possible for a component of the background to have an m ES distribution with a PDF that is more similar to signal (i.e. a Gaussian distribution centered around the B mass) than to a phasespace distribution. Such a component is known as "peaking" background and typically derives from background events that have the same or similar final state particles as the signal decay mode. For example, in B 0 → D + D − , peaking background primarily comes from the decays B 0 → DKX or B 0 → DπX, where D → Kππ and X is K 0 , ρ, a 1 or ω, and the light mesons (KX) or (πX) fake a D → Kππ decay. The optimization procedure that was detailed in Sec. III F eliminates decay submodes that have a large enough amount of peaking (in addition to combinatorial) background to decrease, rather than increase, the sensitivity for a particular decay; the final selection was detailed in Tables II, IV, and V. We determine the amount of peaking background P i in each B decay mode i via fitting the m ES distributions of BB MC simulated events. We minimize the χ 2 i of each fit, allowing the variables κ P i (representing the "ARGUS parameter" described earlier), the number of expected peaking background events in data P i , and the number of phasespace background events N MCbkg i , to float. The fitted number of peaking background events P i is compatible with zero, within two standard deviations, for all modes i.
We then fit the actual data to determine the number of reconstructed signal events in each mode. We fit the m ES distributions of reconstructed B decays that pass all selection criteria in each mode i to a sum of a Gaussian distribution and a phase space distribution (Eq. 5), similar to the PDFs used for efficiency and peaking background fits described above. We minimize the χ 2 i of each data fit, allowing the parameter κ i , the number of signal events in data N sig i , and the number of background events in data N bkg i , each to float. The m ES distributions and the results of the fits are shown in Fig. 4. The branching fractions B i are then determined via the equation where N B = N BB = (231.7 ± 2.6) × 10 6 is the total number of charged or neutral B decays in the data sample, assuming equal production rates of charged and neutral B pairs. We determine the branching fractions as (where ǫ −1 ij is the inverse of matrix ǫ ij ) yields the branching fractions given in Table VII. Note that the measured branching fractions for the three modes B 0 → D ( * )0 D ( * )0 are not significantly greater than zero. Thus, we have determined upper limits on the branching fractions for these modes. The 90% confidence level (C.L.) upper limits quoted in Table VII are determined using the Feldman-Cousins method [19] and include all systematic uncertainties detailed below. Since the branching fractions can be correlated through the use of Eq. 6, we also provide the covariance matrix, with all systematic uncertainties included, in Table VIII. The covariance matrix is obtained via the approximation given in [20].

VI. BRANCHING FRACTION SYSTEMATIC UNCERTAINTIES
4.83 ± 0.78 ± 0.58 [28] < 67 [29] −0.13 ± 0.14 ± 0.02  b. Charged track finding efficiency From studies of absolute tracking efficiency, we assign a systematic uncertainty of 0.8% per charged track on the efficiency of finding tracks other than slow pions from charged D * decays and daughters of K 0 S decays. For the slow pions, we as-sign a systematic uncertainty of 2.2% each, as determined from a separate efficiency study (using extrapolation of slow tracks found in the SVT into the DCH tracking detector and vice-versa). Track finding efficiency uncertainties are treated as 100% correlated among the tracks in a candidate. These uncertainties are weighted by the D and D * branching fractions.
c. K 0 S reconstruction efficiency From a study of the K 0 S reconstruction efficiency (using an inclusive data sample of events containing one or more K 0 S , as well as corresponding MC samples), we assign a 2.5% systematic IX: Estimates of branching fraction systematic uncertainties (as percentages of the absolute values of the branching fraction central values) for all B modes, after propagating the errors through Eq. 6. The totals are the sums in quadrature of the uncertainties in each column. Note that the term "Dalitz weight" refers to the selection on the reconstructed invariant masses of the K − π + and K − π 0 pairs for D 0 → K − π + π 0 decays that was described in Sec. III D.
Mode uncertainty for all modes containing a K 0 S . The value 2.5% comes from the statistical uncertainty in the ratio of data to MC yields and the variation of this ratio over different selection criteria. The uncertainty is weighted by the D and D * branching fractions.
d. π 0 and γ finding efficiency From studies of the neutral particle finding efficiency through the ratios of τ + → ρ + (π + π 0 )ν to τ + → π + ν between data and MC, we assign a 3% systematic uncertainty per π 0 , including the slow π 0 from D * and D * 0 decays. For isolated photons from D * 0 decays, we assign a 1.8% systematic uncertainty, 100% correlated with the π 0 efficiency uncertainty. These uncertainties are weighted by the D and D * branching fractions.
e. Charged kaon identification We assign a systematic uncertainty of 2.5% per charged kaon, according to a study of kaon particle identification efficiency (using kinematically-reconstructed D 0 → K − π + candidates). The uncertainty is weighted by the D and D * branching fractions.
f. Other selection differences between data and MC Differences in momentum measurement, decay vertex finding efficiency, etc., can result in additional differences between efficiencies in data and in MC. We use a sample of the more abundant B 0 → D * + s D * − events in data, selected in a similar manner as the B → D ( * ) D ( * ) modes, to determine these uncertainties. To estimate the sys-tematic error arising from differences between the data and MC D and D * mass resolutions, we calculate the number of D * s D * events seen in the data and MC as a function of the L Mass cut, while fixing the other selection criteria to their nominal values. The number of observed events is extracted from a fit to the m ES distribution. We then plot the ratio of the data yield (N data ) to the MC yield (N MC ) as a function of the L Mass cut over a range of values that gives the same efficiencies as in the D ( * ) D ( * ) analyses. We find the rms of the N data /N MC ratio and assign this as a systematic uncertainty for applying this cut. The same technique is used to determine the systematic uncertainties from all other selection criteria in Table IX: the selections on F , L, ∆E, the reconstructed invariant masses for D 0 → K − π + π 0 ("Dalitz weight"), and vertex P(χ 2 ). g. Fit model The data yield is obtained from an m ES fit where the mean (µ) and width (σ) of the B mass and the end-point ( √ s/2) of the phase-space distribution (Eq. 5) are fixed. These parameters are estimated and have associated uncertainties. The nominal value of σ is determined from signal MC for each B decay mode. To estimate the systematic uncertainty due to possible differences between the m ES resolutions in data and signal MC, we first look at this difference (∆σ = σ data − σ MC ) for those modes with high purity, including our control sample. These differences are consistent with zero, jus-tifying our use of σ MC in obtaining the data yield. We then find the weighted average of ∆σ, which is given by (0.11 ± 0.08) MeV/c 2 . As a conservative estimate, we repeat the data yield determinations by moving σ up and down by 0.2 MeV/c 2 , and take the average of the absolute values of the changes in each data yield as the systematic uncertainty of fixing σ to the MC value for that B mode. A combined fit of common modes in data is used to determine the nominal values for µ and for the endpoint of the m ES distribution √ s/2. Hence, we move the parameters up and down by their fitted errors (0.2 MeV/c 2 for µ and 0.1 MeV/c 2 for √ s/2) to obtain their corresponding systematic uncertainties. The quadratic sum of the three uncertainties from µ, σ and √ s/2 gives the systematic uncertainty of the fit model for each B mode.
h. Spin-alignment dependence The B 0 → D * + D * − , B 0 → D * 0 D * 0 , and B + → D * 0 D * + decays are pseudoscalar → vector vector (VV) transitions described by three independent helicity amplitudes A 0 , A , and A ⊥ [31]. The lack of knowledge of the true helicity amplitudes in the B → V V final states contributes a systematic uncertainty to the efficiency. The dominant source of this effect originates from the p T -dependent inefficiency in reconstructing the low-momentum "soft" pions in the D * + and D * 0 decays, and the fact that the three helicity amplitudes contribute very differently to the slow pion p T distributions. To estimate the size of this effect, MC samples are produced with a phase-space angular distribution model for the decay products. Each event is then weighted by the angular distribution for given input values of the helicity amplitudes and phase differences. The efficiency is then determined for a large number of amplitude sets and the observed distributions in efficiencies are used to estimate a systematic uncertainty. For a given iteration, a random number, based on a uniform PDF, is generated for each of the three parameters: R ⊥ , α, and η, where and η is the strong phase difference between A 0 and A || . Since R ⊥ for B 0 → D * + D * − has already been measured [36], a Gaussian PDF with mean and width fixed to the measured values is used instead for that mode.
The events of the MC sample are weighted by the corresponding angular distribution and the efficiency is determined (after applying all selection cuts) by fitting the m ES distribution and dividing by the number of generated events. The procedure is repeated 1000 times for each B → V V sample. The relative spread in efficiencies (rms divided by the mean) is used to estimate the systematic uncertainty due to a lack of knowledge of the true amplitudes. i. Peaking background and crossfeed The uncertainties on the peaking background vector P i and on the efficiency matrix ǫ ij are dominated by the available MC statistics. The resulting uncertainties on each element of the vector and matrix are propagated through to the branching fraction results via the formalism of Eq. 6.
j. Number of BB The number of BB events in the full data sample, and the uncertainty on this number, are determined via a dedicated analysis of charged track multiplicity and event shape [15]. The uncertainty introduces a systematic uncertainty of 1.1% on each of the branching fractions.

VII. MEASUREMENT OF CP -VIOLATING CHARGE ASYMMETRIES
To obtain the charge asymmetries A CP (defined in Eq. 1), we perform unbinned extended maximum likelihood fits to the m ES distributions of the selected events in each of the four charged-B decay modes D * + D * 0 , D * + D 0 , D + D * 0 , D + D 0 , and their respective charge conjugates, and in the neutral-B decay mode D * ± D ∓ , using Eq. 5 as the PDF for the combinatorial background for both charges in each pair. The free parameters of each of the five fits individually are: 1) the combinatorial background shape parameter κ, 2) the total number of signal events, 3) the total number of background events, and 4) the "raw" charge asymmetry A. Parameters 1 and 3 are considered (and thus constrained to be) the same for both charge states in each mode; this assumption is validated in MC simulation of the background as well as in control samples of B 0 → D * − ρ + and B 0 → D * − a + 1 decays in data. The results of the fits are shown in Fig. 5. Two potentially biasing effects must be considered: there can be a asymmetry in the efficiencies for reconstructing positively-and negatively-charged tracks, and peaking background and crossfeed between the modes can cause a small difference between the measured ("raw") asymmetry and the true asymmetry. The former of those two effects is discussed in Sec. VIII below. Regarding the latter, to obtain the charge asymmetries A CP from the "raw" asymmetries A, very small corrections for peaking background and crossfeed between modes must be made. Using the terminology of Eq. 6, and considering the branching fractions B i to be sums of a "+" mode (with a B 0 or B + , containing ab quark, as the initial state) and a "−" mode (with a B 0 or B − , which contain a b quark, as the initial state): for the "+" and "−" states respectively, which imply As we have . (12) Since N sig j ≡ N sig− j + N sig+ j and the "raw" asymmetry where are the charge asymmetries of the peaking backgrounds. The total yields N sig j , peaking backgrounds P j , and efficiency matrix ǫ ij are identical to those used for the branching fraction measurements and are given in Tables VII and VI. The values A P j are nominally set to 0 and are varied to obtain systematic uncertainties due to the uncertainty on the charge asymmetry of the peaking background (see Sec. VIII). Thus, Eq. 13 is used to determine the final A CP values from the measured asymmetries, in order to account for the small effects due to peaking background and crossfeed between modes. The measured A CP values are given in Table VII. They are all consistent with zero, and their errors are dominated by statistical uncertainty. Table X shows the results of our evaluation of the various sources of systematic uncertainty that are important for the A CP measurements.

VIII. SYSTEMATIC UNCERTAINTIES ON CHARGE ASYMMETRY MEASUREMENTS
a. Slow π ± charge asymmetry A charge asymmetry in the reconstruction efficiency of the low-transversemomentum charged pions from D * ± → D 0 π ± decays can cause a shift in A CP by biasing the rates of positively charged vs. negatively charged decays for each mode. We estimate this systematic uncertainty by using data control samples of B 0 → D * − X + and B 0 → D * + X − decays, where X is either π, ρ, or a 1 , and determining if there is an asymmetry in the number of D * + vs. D * − reconstructed. There are two potential biases of this technique: 1) a charge asymmetry in tracks other than the slow charged pions, and 2) the presence of doubly-Cabibbo-suppressed B 0 (B 0 ) → D * ∓ X ± decays which could potentially introduce a direct-CP -violating asymmetry between the two states in the control sample. Discussion of 1) is detailed in the paragraph below, and the rate of 2) has been determined in analyses such as Refs. [32] and [33] to be of order 0.1%, well below the sensitivity for this measurement. We combine the information from the control sample modes and determine an uncertainty of 0.5% for each A CP measurement for modes with a charged slow pion.
b. Charge asymmetry from tracks other than slow π ± Auxilliary track reconstruction studies place a stringent bound on detector charge asymmetry effects at transverse momenta above 200 MeV/c. Such tracking and PID systematic effects were studied in detail in the analysis of B → φK * [34]. We assign a 0.2% systematic per charged track, thus an overall systematic of 0.4% per mode (as the positively charged and negatively charged decays for each mode have, on balance, one positive vs. one negative track respectively). This systematic uncertainty is added linearly to the slow pion charge asymmetry systematic due to potential correlation.
c. Amount of peaking background Peaking background can potentially bias A CP measurements in two ways: 1) a difference in the total amount of peaking background from the expected total amount can, to second order, alter the measured asymmetry between the positively charged and negatively charged decays, 2) a more direct way for peaking background to alter the measured A CP would be if the peaking background itself were to have an asymmetry between the amount that is reconstructed as positively charged and the amount reconstructed as negative. 1) is discussed here; 2) is discussed in the paragraph below. The systematic uncertainty due to the uncertainty on the total amount of peaking background in the five decays is determined via the formalism of Eq. 13. Namely, the uncertainty is given by where (δP ) j are the uncertainties on the amount of peaking background (which are given, along with the other parameters in the equation, in Table VII).

d.
A CP of peaking background The systematic uncertainty due to the A CP of the peaking background is also determined using the formalism of Eq. 13. Namely, the uncertainty is given by .
Investigation of the sources of the peaking background in these modes motivates a conservative choice of 0.68 for the (δA P ) j values.
e. Amount of crossfeed The systematic error due to uncertainties in the amount of crossfeed between the modes is also determined via the formalism of Eq. 13. Namely, the uncertainty is given by The covariance between the elements of the inverse efficiency matrix is obtained using the method of Ref. [20]. The very small systematic uncertainty due to crossfeed is thus obtained using Eq. 16 and the amounts of crossfeed and their uncertainties that are given in Table VI. f. Uncertainty in m ES resolution, B mass, and √ s The uncertainties in m ES resolution and the beam energy √ s are determined by varying these parameters within their fitted ±1σ ranges and observing the resulting changes in A CP . The uncertainty in the reconstructed B mass can also have an impact on the fitted m ES distributions and thus on the fitted A CP values. Varying the B mass between the fitted value and the ±1σ range of the nominal B 0 or B + invariant mass allows the determination of the resulting effect on the A CP values.
g. Potential fit bias Uncertainties in the potential biases of the A CP fits are determined by performing the fits on large samples of MC simulation of the signal decay modes and of BB and continuum background decays. All results are consistent with zero bias, and the uncertainties of the fitted asymmetries on the simulated data samples are conservatively assigned as systematic uncertainties from biases of the fits. s D ( * ) decays [7,8]. For this analysis, we assume that the breaking of SU(3) can be parametrized via the ratios of decay constants f D ( * ) s /f D ( * ) , which are quantities that can be determined either with lattice QCD or from experimental measurements [9].
In this model, one obtains the relation (for B 0 → D + D − and individual helicity states of B 0 → D * + D * − ): where and A D andĀ D represent amplitudes of a given B 0 and B 0 → D ( * )+ D ( * )− decay respectively, B represents the corresponding average branching fraction, and a dir and a indir represent the corresponding direct and indirect CP asymmetries respectively. The phases β and γ are the CKM phases and δ is a strong phase difference.
are the magnitudes of the combined B → D ( * ) D ( * ) decay amplitudes containing V * cb V cd and V * ub V ud terms respectively, and the T , P , and E terms are the tree, penguin, and the sum of exchange and annihilation amplitudes respectively [7]. One can directly measure the parameters B, a dir , and a indir using information from B → D ( * ) D ( * ) decays; the parameter A ct using information from B → D ( * ) s D ( * ) decays; and the weak phase β can be obtained from the measurements of sin 2β based on B 0 → ccK 0 S decays [35] thus allowing for solution of γ (up to two discrete ambiguities) via Eq. 17. As the vector-pseudoscalar modes B 0 → D * ± D ∓ are not CP eigenstates, a slightly more complicated analogue to Eq. 17 is needed for these modes [8]. Measurement of A CP for D * ± D ∓ is also necessary to obtain information on γ from the vector-pseudoscalar modes.
Using these relations, there are four variables besides β for each B → D ( * ) D ( * ) decay for which to solve: A ct , A ut , δ, and γ. The branching fraction and the direct and indirect CP asymmetries of the B → D ( * ) D ( * ) decay provide three measured quantities. The other measurement that can be used is the branching fraction of the corresponding B → D where θ c is the Cabibbo angle [14] and the parentheses around the asterisks correspond to the B → D ( * ) D ( * ) and B → D ( * ) s D ( * ) decays that are used. The theoretical uncertainty of this relation is determined to be 10% [7].
To use the VV decays, we must make the assump-tion that the strong phases for the 0 and helicity amplitudes are equal. The constraints from the PP decays require no such assumption. The assumption of equal 0 and helicity amplitudes is theoretically supported by a QCD factorization argument described in [8].
We use a fast parametrized MC method, described in Ref. [8], to determine the confidence intervals for γ. We consider 500 values for γ, evenly spaced between 0 and 2π. For each value of γ considered, we generate 25000 MC experiments, with inputs that are generated according to Gaussian distributions with widths equal to the experimental errors of each quantity. For each experiment, we generate random values of each of the experimental inputs according to Gaussian distributions, with means and sigmas according to the measured central value and total errors on each experimental quantity. We make the assumption that the ratio f D * s /f D * is equal to f Ds /f D = 1.20 ± 0.06 ± 0.06 [9], allowing for the additional 10% theoretical uncertainty [7]. We then calculate the resulting values of A ct , a dir , a indir , and B, given the generated random values (based on the experimental values). When the quantities a dir , a indir , and B, along with β and the value of γ that is being considered, are input into Eq. (17), we obtain a residual value for each experiment, equal to the difference of the left-and right-hand sides of the equation. Thus, using Eq. 17, the 25000 trials per value of γ provide an ensemble of residual values that are used to create a likelihood for γ to be at that value, given the experimental inputs. The likelihood, as a function of γ, can be obtained from χ 2 (γ), where χ 2 ≡ (µ/σ) 2 , µ is the mean of the above ensemble of residual values, and σ is the usual square root of the variance. The value of χ 2 (γ) is then considered to represent a likelihood which is equal to that of a value χ standard devations of a Gaussian distribution from the most likely value(s) of γ. We define the "exclusion level," as a function of the value of γ, as follows: the value of γ is excluded from a range at a given C.L. if the exclusion level in that range of γ values is greater than the given C.L.
We now turn to the VP decays. The method using VP decays shares the advantage with PP decays that no assumptions on strong phases are required. The disadvantage is that, as we will see, the constraints from the VP modes are weak.
We combine the information given above on the B 0 → D * ± D ∓ , B + → D * + D 0 , and B + → D + D * 0 branching fractions and A CP information with measurements of the B 0 → D * − s D + , B 0 → D − s D * + , B + → D * + s D 0 , and B + → D + s D * 0 branching fractions [14], measurements of the B 0 → D * ± D ∓ time-dependent CP asymmetries [24,37], and the world-average values of sin 2β [35] and sin 2 θ c [14]. Similar to the MC γ determination for the VV and PP modes, we generate random values of each of the experimental inputs according to Gaussian distributions, with means and sigmas according to the measured central value and total errors on each experimental quantity. We again obtain a confidence level distribution as a function of γ.
Finally, we can combine information from the VV, PP, and VP modes. The resulting measured exclusion level as a function of γ from each of the three sets of modes, as well as from their combination, is shown in Fig. 6. From the combined fit, we see that γ is favored to lie in the range [0.07 − 2.77] radians (with a +0 or +π radians ambiguity) at 68% confidence level. This corresponds to [4.1 • − 158.6 • ](+0 • or 180 • ).
These constraints are generally weaker than those found in Ref. [8] due to the fact that the measured CP asymmetry in B 0 → D * + D * − has moved closer to the world-average sin2β, with the newer B 0 → D * + D * − measurements in Ref. [38]. The closer this CP asymmetry is to sin2β, the weaker the resulting constraints are on γ, due to the fact that the closeness of the CP asymmetry to sin2β favors the dominance of the tree amplitude, rather than the penguin amplitude whose phase provides the sensitivity to γ. Although the constraints are not strong, they contribute to the growing amount of information available on γ from various sources.

X. CONCLUSIONS
In summary, we have measured branching fractions, upper limits, and charge asymmetries for all B meson decays to D ( * ) D ( * ) . The results are shown in Table VII. This includes observation of the decay modes B 0 → D + D − and B + → D * + D * 0 , evidence for the decay modes B + → D + D * 0 and B + → D + D 0 at 3.8σ and 4.9σ levels respectively, constraints on CP -violating charge asymmetries in the four decay modes B + → D ( * )+ D ( * )0 , measurements of (and upper limits for) the decay modes B 0 → D * 0 D 0 and B 0 → D 0 D 0 , and improved branch-ing fractions, upper limits, and charge asymmetries in all other B → D ( * ) D ( * ) modes. The results are consistent with theoretical expectation and (when available) previous measurements. When we combine information from time-dependent CP asymmetries in B 0 → D ( * )+ D ( * )− decays [38,39] and world-averaged branching fractions of B decays to D ( * ) s D ( * ) modes [14] using the technique proposed in Ref. [7] and implemented in Ref. [8], we find the CKM phase γ is favored to lie in the range [0.07−2.77] radians (with a +0 or +π radians ambiguity) at 68% confidence level.

XI. ACKNOWLEDGMENTS
We are grateful for the extraordinary contributions of our PEP-II colleagues in achieving the excellent luminosity and machine conditions that have made this work possible. The success of this project also relies critically on the expertise and dedication of the computing organizations that support BABAR. The collaborating institutions wish to thank SLAC for its support and the kind hospitality extended to them. This work is supported by the US Department of Energy and National Science Foundation, the Natural Sciences and Engineering Research Council (Canada), Institute of High Energy Physics (China), the Commissariatà l'Energie Atomique and Institut National de Physique Nucléaire et de Physique des Particules (France), the Bundesministerium für Bildung und The measured exclusion level, as a function of γ, from the combined information from vector-vector, vector-pseudoscalar, and pseudoscalar-pseudoscalar modes. The combined information implies that γ is favored to lie in the range [0.07 − 2.77] radians (with a +0 or +π radians ambiguity) radians at 68% confidence level.