Study of the decay mode D^0 ->K-K-K+pi+

Using data from the FOCUS (E831) experiment at Fermilab, we present a new measurement of the branching ratio for the Cabibbo-favored decay mode $D^0 \to K^-K^-K^+\pi^+$. From a sample of $143 \pm 19$ fully reconstructed $D^0 \to K^-K^-K^+\pi^+$ events, we measure $\Gamma(D^0 \to K^-K^-K^+\pi^+)/\Gamma(D^0 \to K^-\pi^-\pi^+\pi^+) = 0.00257 \pm 0.00034(stat.) \pm 0.00024(syst.)$. A coherent amplitude analysis has been performed to determine the resonant substructure of this decay mode. This analysis reveals a dominant contribution from $\phi$ and $\bar K^{*0}(890)$ states.


Introduction
Hadronic decays of charm mesons have been extensively studied in recent years. Dalitz plot analyses of D meson decays in three-body final states have revealed a rich resonant substructure, showing a dominance of quasi twobody modes. However, much less information is available on the resonant substructure of decays with more than three final state particles. These multibody D decays account for a large fraction of the hadronic decay width. If the amplitude analyses of multi-body decays confirm the picture drawn by the Dalitz plot analyses, then one can make a more complete comparison with theoretical models, which have been developed mainly to describe the two-body and quasi-two-body decay modes [1,2,3,4,5].
We present a new study of the D 0 → K − K − K + π + decay using data from the FOCUS experiment. This is an interesting decay mode: although Cabibbo favored, it is strongly suppressed by phase-space and requires the production of an ss pair, either from the vacuum or via final state interactions (FSI). We measure the branching ratio Γ(D 0 →K − K − K + π + ) Γ(D 0 →K − π − π + π + ) and perform a coherent amplitude analysis to determine its resonant substructure.
FOCUS is a charm photoproduction experiment [6] which collected data during the 1996-97 fixed-target run at Fermilab. Electron and positron beams (with typically 300 GeV endpoint energy) obtained from the 800 GeV/c Tevatron proton beam produce, by means of bremsstrahlung, a photon beam which interacts with a segmented BeO target. The mean photon energy for triggered events is ∼ 180 GeV. A system of three multicell thresholdČerenkov counters performs the charged particle identification, separating kaons from pions with momenta up to 60 GeV/c. Two systems of silicon microvertex detectors are used to track particles: the first system consists of 4 planes of microstrips interleaved with the experimental target [7]; the second system consists of 12 planes of microstrips located downstream of the target. These detectors provide high resolution in the transverse plane (approximately 9 µm), allowing the identification and separation of charm primary (production) and secondary (decay) vertices. The charged particle momenta are determined by measuring their deflections in two magnets of opposite polarity through five stations of multiwire proportional chambers.

Analysis of the decay mode
The final states are selected using a candidate driven vertex algorithm [6]. A secondary vertex is formed from the four candidate tracks. The momentum of the resultant D 0 candidate is used as a seed track to intersect the other reconstructed tracks to search for a primary vertex. The confidence levels of both vertices are required to be greater than 1%. Once the production and decay vertices are determined, the distance ℓ between them and its error σ ℓ are computed. The quantity ℓ/σ ℓ is an unbiased measure of the significance of detachment between the primary and secondary vertices. This is the most important variable for separating charm events from non-charm, prompt backgrounds. Signal quality is further enhanced by cutting on Iso2, which is the confidence level that other tracks in the event might be associated with the secondary vertex. To minimize systematic errors on the measurements of the branching ratio, we use identical vertex cuts on the signal and normalizing mode, namely ℓ/σ ℓ > 6, and Iso2 < 1 %. We also require the D 0 momentum to be in the range 25-250 GeV/c (a very loose cut) and the primary vertex to be formed with at least two reconstructed tracks in addition to the D 0 seed.
The only difference in the selection criteria between the two decay modes lies in the particle identification cuts. TheČerenkov identification cuts used in FOCUS are based on likelihood ratios between the various particle identification hypotheses. These likelihoods are computed for a given track from the observed firing response (on or off) of all cells within the track's (β = 1) Cerenkov cone for each of our threeČerenkov counters. The product of all firing probabilities for all the cells within the threeČerenkov cones produces a χ 2 -like variable W i = −2 ln(Likelihood) where i ranges over the electron, pion, kaon and proton hypotheses [8]. All kaon tracks are required to have ∆ K = W π − W K (kaonicity) greater than 2, whereas the pion tracks are required to have ∆ π = W K − W π (pionicity) exceeding 0.5.
Using the set of selection cuts just described, we obtain the invariant mass distributions for K − K − K + π − and K − π − π + π + shown in Fig. 1. In Fig. 1a the K − K − K + π − mass plot is fit with two Gaussians with the same mean but different sigmas to take into account the change in resolution with momentum of our spectrometer [6] plus a second-order polynomial. A log-likelihood fit returns a signal of 143 ± 19 D 0 → K − K − K + π − events. The large statistics K − π − π + π + mass plot of Fig. 1b is fitted in the same way (two Gaussians plus a second-order polynomial). The fit gives a signal of 64576 ± 360 D 0 → K − π − π + π + events.

Relative Branching Ratio
The evaluation of relative branching ratios requires yields from the fits to be corrected for detection efficiencies, which differ among the various decay modes due to differences in both spectrometer acceptance (due to different Q values and resonant substructure of the two decay modes) andČerenkov identification efficiency.
From Monte Carlo simulations, we compute the relative efficiency to be: ǫ(D 0 →K − π − π + π + ) = 0.862 ± 0.010. Combined with the yields measured above we find a branching ratio of Our final measurement has been tested by modifying each of the vertex anď Cerenkov cuts individually. The branching ratio is stable versus several sets of cuts as shown in Fig. 2 (we vary the confidence level of the secondary vertex from 1% to 10%, Iso2 from 10 −6 to 1, ℓ/σ ℓ from 5 to 15, ∆ K from 0.5 to 4.5 and finally ∆ π from 0.5 to 4.5).
Systematic uncertainties on branching ratio measurements can come from different sources. We determine four independent contributions to the systematic uncertainty: the split sample component, the fit variant component, the component due to the particular choice of the vertex andČerenkov cuts (discussed previously), and the limited statistics of the Monte Carlo.
The split sample component takes into account the systematics introduced by a residual difference between data and Monte Carlo, due to either a possible mismatch in the reproduction of the D 0 momentum or the change in the experimental conditions of the spectrometer during data collection. This component has been determined by splitting data into four independent subsamples, according to the D 0 momentum range (high and low momentum) and the configuration of the vertex detector, that is, before and after the insertion of an upstream silicon system [7]. A technique employed in FOCUS and in the predecessor experiment E687, modeled after the S-factor method from the Particle Data Group [9], is used to try to separate true systematic variations from statistical fluctuations. The branching ratio is evaluated for each of the 4 statistically independent sub-samples and a scaled variance is calculatedσ (the errors are boosted when χ 2 /(N − 1) > 1). The split sample variance σ split is defined as the difference between the reported statistical variance and the scaled variance, if the scaled variance exceeds the statistical variance [10].
Another possible source of systematic uncertainty is the fit variant. This component is computed by varying, in a reasonable manner, the fitting conditions on the whole data set. In our study we fixed the widths of the Gaussians to the values obtained by the Monte Carlo simulation, changed the background parametrization (varying the degree of the polynomial), and used one Gaussian instead of two. In addition we considered the variation of the computed efficiency, both for D 0 → K − K − K + π − and the normalizing decay mode, due to the different resonant substructure simulated in the Monte Carlo. The BR values obtained by these variants are all a priori likely; therefore this uncertainty can be estimated by the rms of the measurements [10].
Analogous to the fit variant, the cut component is estimated using the standard deviation of the values obtained from the many sets of cuts shown in Fig. 2. Actually this is an overestimate of this component because the event samples of the various cut sets are different.
Finally, there is a contribution due to the limited statistics of the Monte Carlo simulation used to determine the efficiencies. The resulting systematic errors are summarized in Table 1. Adding in quadrature the four components, we obtain the total systematic error also shown in Table 1.
The final branching ratio result is shown in Table 2 along with a comparison to previous measurements.
4. Amplitude analysis of D 0 → K − K − K + π + A fully coherent amplitude analysis was performed to determine the reso- Events E687 [11] 0.0028 ± 0.0007 ± 0.0001 20 ± 5 E791 [12] 0.0054 ± 0.0016 ± 0.0008 18.4 ± 5.3 FOCUS (this result) 0.00257 ± 0.00034 ± 0.00024 143 ± 19 Table 2 Branching ratio measurement and comparison with other experiments. nant substructure of the D 0 → K − K − K + π + decay. While a large number of intermediate states could lead to the K − K − K + π + final state, phase space limitations restrict the possible contributions. A plot of the K − K + invariant mass shows a clear φ contribution (Fig. 3, top left plot). It is also possible, in principle, to have K − K + contribution via f 0 (980) and a 0 (980). There are large uncertainties in the line shape of these resonances and in their coupling to K − K + . They do not appear to be required by the fit and are therefore not included.
Contributions from κ(800)K − K + and K * 0 (1430)K − K + might also be present. However, these resonances are broad scalar states, with no characteristic angular distribution that could distinguish them from the non-resonant mode, with the present level of statistics.
Modes containing a K * 0 (892) can also contribute, even though the nominal K * 0 (892) mass is just above the kinematical limit of the K − π + spectrum. Since K * 0 (892) is a narrow vector meson, its contribution can be distinguished by the angular distribution of the decay products, even without a clear mass peak.
The formalism used in this amplitude analysis is a straightforward extension to four-body decays of the usual Dalitz plot fit technique. The D 0 is a spin zero particle, so the four-body decay kinematics are defined by five degrees of freedom.
Individual amplitudes, A k , for each resonant mode are constructed as a product of form factors, relativistic Breit-Wigner functions, and spin amplitudes which account for angular momentum conservation. We use the Blatt-Weisskopf damping factors [14], F l , as form factors (l is the orbital angular momentum of the decay vertex). For the spin amplitudes we use the Lorentz invariant amplitudes [13], which depend both on the spin of the resonance(s) and the orbital angular momentum. The relativistic Breit-Wigner is In the above equations m is the two-body invariant mass, m 0 and s are the resonance nominal mass and spin, and p * = p * (m) is the breakup momentum at resonance mass m.
Since there are two identical kaons, each amplitude A k is Bose-symmetrized. The overall signal amplitude is a coherent sum of the individual amplitudes, A = k c k A k , assuming a constant complex amplitude for the non-resonant mode. The coefficients c k are complex numbers to be determined by the fit. The overall signal amplitude is corrected on an event-by-event basis for the acceptance, which is nearly constant across the phase space.
In the amplitude analysis we have taken events having a K − K − K + π + invariant mass within the interval M D ± 10 MeV/c 2 . In this interval there are 139 signal and 65 background events. The finite detector resolution causes a smearing of the edges of the five-dimensional phase space. This effect is accounted for by multiplying the overall signal distribution by a Gaussian factor, g(M), where M is the K − K + K − π + mass. The normalized signal probability distribution is where φ represents the coordinates of an event in the five dimensional phase space, ε(φ) is the acceptance function, and ρ(φ) is the phase space density.
Two types of background events were considered: random φ's combined with a K − π + pair, and random combinations of K − K − K + π + . Inspection of the side bands of the K − K − K + π + mass spectrum indicate that nearly 30% of the background events are of the former type. We assume the random K − K − K + π + combinations to be uniformly distributed in phase space, while for the φ background we assume an incoherent sum of Breit-Wigners with no form factors and no angular distribution. The overall background distribution, kept fixed in the fit, is a weighted, incoherent sum of these two components. The overall background distribution is also corrected for the acceptance (assumed to be the same as for the signal events) on an event-by-event basis, and multiplied by an exponential function b(M), to account for the detector resolution. The normalized background probability distribution is An unbinned maximum likelihood fit was performed, minimizing the quantity w ≡ −2 ln L. The likelihood function, L, is Neither the acceptance function nor the phase space density depend on the fit parameters c k , so the term −2ln [ε(φ)ρ(φ)] is irrelevant to the minimization. The acceptance correction is important only for the normalization integrals N S and N B .
Decay fractions are obtained from the coefficients c k determined by the fit, after integrating the overall signal amplitude over the phase space: We fit the data to both models A and B. We find that model B, with only the non-resonant and φK − π + contributions, does not provide a good description of the data. The inclusion of the two K * 0 (892) amplitudes results in a much better fit, with an improvement in ∆w of 138.
Results from the fit with model A are shown in Table 3. The dominant contribution comes from D 0 → φK * 0 (892). Adding the contributions from the two K * 0 (892) amplitudes, we see that they account for nearly 70% of the total decay width. This is somewhat surprising, given the very small phase space. In Fig. 3 the K + K − and K − π + projections of events used in the amplitude analysis are superimposed on the fit result (top two plots). The projections from the background model are shown in the shaded histograms. The remaining plots are two-dimensional projections of the events in the signal region.
We have performed a log-likelihood test to check our ability to distinguish between models A and B. Two ensembles of 10,000 mini-MC samples were generated, one simulated according to model A and another according to model B. For each sample in each ensemble we compute the quantity ∆w = 2 ln L A − 2 ln L B , where L B and L A are the likelihoods calculated with models B and A, respectively. The resulting ∆w distributions are shown in Fig. 4. On the right we see the ∆w distribution computed from the model A ensemble, and, on the left, ∆w with the model B ensemble. The two distributions are well separated showing that we can easily distinguish between these two models. Moreover, the value of ∆w obtained from the real data (138) is consistent with the distribution from model A, showing that this is indeed a better description of the data.
The goodness-of-fit was assessed in two ways. We have estimated the confidence level of the fit using a χ 2 test. The five invariants used to define the kinematics of this decay are the four K − K + and K − π + masses squared, plus either the K − K − or K + π + mass squared. Due to the limited statistics we have integrated over the latter invariant and divided the other four into two bins, yielding a total of sixteen cells. A χ 2 was computed and the estimated confidence level was 35%. The confidence level obtained with model B was 6×10 −11 .
Given the limited statistics we have also estimated the confidence level using a method which is often less stringent than the χ 2 . In this second method the confidence level is estimated using the distribution of w = −2 ln L from an ensemble of mini-MC samples generated with the parameters of model A. This distribution is approximately Gaussian. A confidence level can be estimated by the fraction of samples in which the value of w exceeds that of the data. With this technique we estimate a CL of 86% for the fit with model A.
As in the branching ratio measurement, we consider split sample and fit variant systematic uncertainties, using for the former the same sub-samples described previously. The dominant contributions from fit variant systematic errors come from variations on the signal amplitudes (removing the Blatt-Weisskopf form factors and replacing the non-resonant amplitude by K * 0 (1430)K + K − ), variations in the background relative fractions, and variations of the analysis cuts. Split sample and fit variant errors are added in quadrature. Table 3 shows the Model A fit results. The systematic errors are the second errors in the fractions and phases. φK − π + 0.60 ± 0.12 194 ± 24 ± 8 18 ± 6 ± 4 K * 0 (892)K + K − 0.65 ± 0.13 255 ± 15 ± 4 20 ± 7 ± 2 non-resonant 0.55 ± 0.14 278 ± 16 ± 42 15 ± 6 ± 2 Table 3 Results from the best fit (Model A). The second error on the fractions and phases is systematic.

Conclusions
Using data from the FOCUS (E831) experiment at Fermilab, we have studied the Cabibbo-favored decay mode D 0 → K − K − K + π + .
A comparison with the two previous determinations of the relative branching ratio Γ(D 0 →K − K − K + π + ) Γ(D 0 →K − π − π + π + ) shows an impressive improvement in the accuracy of this measurement.