Determining the unitarity triangle angle $\gamma$ with a four-body amplitude analysis of $\rm B^\pm \to (K^+K^-\pi^+\pi^-)_D K^\pm$ decays

We explain how a four-body amplitude analysis of the $\rm D$ decay products in the mode $\rm B^\pm \to (K^+K^-\pi^+\pi^-)_D K^\pm$ is sensitive to the unitarity triangle angle $\gamma$. We present results from simulation studies which show that a precision on $\gamma$ of $15^\circ$ is achievable with 1000 events and assuming a value of 0.10 for the parameter $r_B$.


Introduction
A precise measurement of the unitarity triangle angle γ is one of the most important goals of CP violation experiments. γ is defined as arg(−V * ub V ud /V * cb V cd ), where V ij are the elements of the Cabibbo-Kobayashi-Maskawa (CKM) mixing matrix. In the Wolfenstein convention [1] γ = arg(V * ub ). A class of promising methods to measure γ exists which exploits the interference between the amplitudes leading to the decays B − → D 0 K − and B − →D 0 K − (Figure 1), where the D 0 andD 0 are reconstructed in a common final state. This final state may be, for example, a CP eigenstate such as K + K − ('GLW method') [2], or a non-CP eigenstate such as K + π − , which can be reached both through a doubly Cabibbosuppressed D 0 decay and a Cabibbo-favouredD 0 decay ('ADS method') [3]. Recent attention has focused on self-conjugate three-body final states, in particular D → K S π + π − 1 . Here a Dalitz analysis of the resonant substructure in the K S π + π − system allows γ to be extracted [4]. The B-factory experiments have used this method to obtain the first interesting direct constraints on γ [6,7].
Here we explore the potential of determining γ through a four-body amplitude analysis of the D decay products in the mode B ± → (K + K − π + π − ) D K ± . CP studies involving amplitude analyses of four-body systems have been 1 Here and subsequently D signifies either a D 0 or aD 0 .
proposed elsewhere [4], and strategies already exist for B ± → DK ± approaches exploiting singly-Cabibbo suppressed decays [5]. Our method benefits from a final state that involves only charged particles, which makes it particularly suitable for experiments at hadron colliders, most notably LHCb.
This paper is organised as follows. In Section 2 we summarise the essential features of B ± → DK ± decays, and state the present knowledge of the parameters involved, and of the decay D → K + K − π + π − . In Section 3 a full model of B ± → (K + K − π + π − ) D K ± decays is developed, which is then used within a simulation study to estimate the precision on γ which may be obtained through a four-body amplitude analysis. We conclude in Section 4.
2. B ± → DK ± decays and D → K + K − π + π − Let us define the amplitudes of the two diagrams illustrated in Figure 1 as follows Here the strong phase of A B is set to zero by convention, and δ B is the difference of strong phases between the two amplitudes. γ represents the weak phase difference between the amplitudes, where contributions to the CKM elements of order λ 4 and higher (with λ being the sine of the Cabibbo angle) have been neglected. In the CP- conjugate transitions γ → −γ, whereas δ B remains unchanged. r B is the relative magnitude of the colour-suppressed B − →D 0 K − process to the colour-favoured B − → D 0 K − transition. Preliminary indications as to the values of γ, δ B and r B come from the B ± → DK ± , D → K S π + π − analyses performed at the B-factories [6,7]. Fits to the ensemble of hadronic flavour data also provide indirect constraints on the value of γ [8,9]. These results lead us to assume values of γ = 60 • and δ B = 130 • for the illustrative sensitivity studies presented in Section 3. We set r B to 0.10, which is the approximate average of the Dalitz results and the lower values favoured by the ADS and GLW analyses [10,11].
Results have recently been reported from an amplitude analysis of the decay D → K + K − π + π − [12], which shows that the dominant contributions come from D → AP and D → VV modes. Earlier measurements of D → K + K − π + π − were published in [13]. Our sensitivity studies for the γ extraction, presented in Section 3, are based on the results found in [12].
The branching ratio of the mode B ± → (K + K − π + π − ) D K ± can be estimated as the product of the two meson decays, and found to be 9.2 × 10 −7 [14]. This channel is particularly well matched to the LHCb experiment, on account of the kaon-pion discrimination provided by the RICH system, and the absence of any neutrals in the final state, which allows for good reconstruction efficiency and powerful vertex constraints. Consideration of the trigger and reconstruction efficiencies of similar topology decays reported in [15] leads to the expectation of sample sizes of more than 1000 events per year of operation.

Estimating the γ sensitivity in
In this section we formulate a model to describe B ± → (K + K − π + π − ) D K ± decays. This model neglects D 0 −D 0 oscillations and CP violation in the D system, which is a good approximation in the Standard Model. The model is then used in a simulation study to estimate the sensitivity with which γ can be determined from an analysis of B ± → (K + K − π + π − ) D K ± events.

Decay Model
In the same way as the kinematics of a threebody decay can be fully described by two variables (Dalitz Plot), typically s 12 = (p 1 + p 2 ) 2 , s 23 = (p 2 + p 3 ) 2 , where p 1 , p 2 , p 3 are the 4momenta of the final state particles, so can a four-body decay be described by five variables. In this paper we use the following convention for labelling the particles involved in the D decay and their 4-momenta: We also define: We then choose a set of five variables to describe the decay kinematics: t 01 = s 234 , s 12 , s 23 , s 34 and t 40 = s 123 . From these variables all other invariant masses s ij , s ijk , and, for a given frame of reference, all momenta p i can be calculated.
In contrast to the phase space density for threebody decays, which is uniform in terms of the usual parameters s 12 , s 23 , four-body phase space density, dφ/dt 01 ds 12 ds 23 ds 34 dt 40 , is not flat in 5 dimensions, but proportional to the square-root of the inverse of the 4-dimensional Grahm determinant [17]: The total decay amplitude for the D 0 decay to the K + K − π + π − final state is the sum over all individual amplitudes A k to each set of intermediate states k, weighted by a complex factor |c k |e iφ k An analysis of the D 0 → K + K − π + π − decay amplitude is reported in [12], which fits 10 separate contributions. In this analysis, however, no distinction is made between the modes D 0 → K 1 (1270) + K − , K 1 (1400) − K + and K * (892) 0 K − π + , and those decays to the CPconjugate final states. In our study we base |c k | and φ k on the values found in [12], but consider different scenarios for the relative contributions of the above modes. In order to label these scenarios we make the definitions and We define similar variables for the D 0 → Our default scenario assumes the arbitrary values The amplitudes A k are constructed as a product of form factors (F l ), relativistic Breit-Wigner functions (BW ), and spin amplitudes (s l ) which account for angular momentum conservation, where l is the angular momentum of the decay vertex. Therefore the decay amplitude with a single resonance is given by (where the subscript k has now been omitted), and a decay amplitude with two resonances α and β is written For F l we use Blatt-Weisskopf damping factors [16] and for s l we use the Lorentz invariant amplitudes [18], which depend both on the spin of the resonance(s) and the orbital angular momentum.
With these definitions, the total decay amplitude for B − → DK − , D → K + K − π + π − is given by The corresponding expression for the CP conjugate decay is s 12 , s 23 , s 34 , t 40 ) .
The total probability density function for a B − → (K + K − π + π − ) D K − event is then given by (with an equivalent expression for B + decays) where N is an appropriate normalisation factor which may be obtained through numerical integration.

Simulation Study
To estimate the statistical precision achievable with this method, we generated several simulation samples which we then fitted to determine the parameters of interest, most notably γ. The samples were generated neglecting background and detector effects. Figure 2 shows the projections of the chosen kinematical variables for 200k events, separately for B + and B − decays, and the CP-asymmetry, defined as the number of B + events minus the number of B − events, normalised by the sum. In these projections the observable CP violation is small, typically being at the few percent level only. Full sensitivity to γ is obtained through a likelihood fit to all five variables.
A log-likelihood function is defined as: where the probability density functions are defined as in expression 12, and the sums run over all B candidates in the sample. The function was maximised for each sample using the MINUIT package [19], with r B , δ B and γ as the free parameters. (It is assumed that all parameters associated with the D decay model are known.) Each sample contained 1000 events. A scan of the negative log-likelihood, plotted for γ against δ B is shown for a typical sample in Figure 3. The function is well behaved with a minimum close to the input value and a second solution at γ − 180 • and δ B − 180 • . Also shown is a scan for γ against r B . For each sample the fitted parameters and the assigned errors were recorded. The reliability of the fit result was studied for each variable by constructing the 'pull distribution', which is the   Table 2 Dependence of γ fit results on the value of r B , showing the average assigned error and the means and widths of the pull distributions. difference between the fitted and input parameter, divided by the assigned error. The means and RMS widths of the pull distributions are displayed in Table 1 and are seen to be compatible with 0 and 1 respectively. This indicates that the log-likelihood fit is unbiased and the returned errors are reliable. The fit errors are also included in Table 1, averaged over all fits. γ is extracted with a precision of 14 • . There is very little correlation between the three fit parameters, as is clear from the contours in Figure 3. The size of the interference effects in B ± → DK ± decays, and hence the sensitivity of the fit to γ, depends on the value of r B . To investigate this dependence several 1000 event samples were generated with different values of r B between 0.05 and 0.20. These samples were then fitted as previously. The fit result on γ and associated uncertainty for each r B value are shown in Table 2. It can be seen that the γ error varies approximately linearly with the inverse of r B .
As explained in Section 3.1 the fitted model reported in [12] does not distinguish between the Table 3 Statistical uncertainty on γ for various values of R and ∆φ. These parameters are defined in expressions 6 and 7 with the same values being used for K 1 (1270)K, K 1 (1270)K and K ⋆ (892) 0 Kπ. relative contribution of certain D decay amplitudes and their CP-conjugate final states. The importance of this unknown information on the fit sensitivity was assessed by generating and fitting 1000 event simulated datasets with different values of the R and ∆φ parameters defined in expressions 6 and 7. In varying these parameters the overall contribution of each mode and its CPconjugate state, eg. |A(D 0 → K 1 (1270) The results for the uncertainty on γ are shown in Table 3 in the case where a common value of R and ∆φ is taken for the three final states under consideration. In a further study the phase shift ∆φ was allowed to take different values between the three modes. Four scenarios were considered with the following arbitrary (randomly chosen) sets of values for ∆φ K1(1270)K , ∆φ K1(1400)K and ∆φ K ⋆ (892) 0 Kπ respectively: The statistical uncertainties found on γ for these scenarios are given in Table 4. For both Table 3 and Table 4 only a single experiment was performed at each point in parameter space, hence the stated error carries an uncertainty of a few degrees. However any minor variation in result arising from the exact value of the fitted r B parameter, experiment-to-experiment, has been corrected for by using the dependence observed in the study reported in Table 2.
It can be seen that the precision of the fit is fairly uniform over parameter space, with a typical value of 15 • . In certain cases however the precision is worse, particularly when R = 1 and/or ∆φ = 0. More detailed studies of D → K + K − π + π − decays are therefore needed to reliably estimate the instrinsic sensitivity of B ± → (K + K − π + π − ) D K ± for a γ measurement. However, variations in other aspects of the D → K + K − π + π − decay structure were found to have limited consequences for the fit precision.
Finally it was investigated what biases would be introduced in the γ extraction through incorrect knowledge of the decay model. Experiments were performed in which the datasets were generated with the full model in the default scenario, but fitted with a model which omitted all the decay amplitudes with a contribution less than 3% to the overall rate. Shifts of up to 8 • were observed in the measured value of γ. This value can be considered as an upper bound to any final systematic uncertainty, as it will be possible to accumulate very large samples of D → K + K − π + π − events at the LHC, which will allow the decay model to be refined and improved with respect to the one assumed here. Additional information will also become available from CP-tagged D decays at facilities operating at the ψ(3770) resonance [20].

Conclusions
We have shown that the decay B ± → (K + K − π + π − ) D K ± can be used to provide an interesting measurement of the unitarity triangle angle γ. With 1000 events and assuming a value of r B = 0.10 it is possible to measure γ with a precision of around 15 • . The exact sensitivity achievable depends on the relative contributions of certain unmeasured modes in the D decay model. The final state, involving only charged particles, and kaons in particular, is well suited to LHCb. A full reconstruction study is necessary to estimate reliably the expected event yields and the level of background.
Finally we remark that the same technique of a four-body amplitude analysis in B ± → DK ± decays can be applied to other modes, most notably the 'ADS' channel D → K ± π ∓ π + π − .