Model-independent determination of the strong phase difference between $D^0$ and $\bar{D}^0 \to\pi^+\pi^-\pi^+\pi^-$ amplitudes

For the first time, the strong phase difference between $D^0$ and $\bar{D}^0\to\pi^+\pi^-\pi^+\pi^-$ amplitudes is determined in bins of the decay phase space. The measurement uses $818\,\mathrm{pb}^{-1}$ of $e^+e^-$ collision data that is taken at the $\psi(3770)$ resonance and collected by the CLEO-c experiment. The measurement is important for the determination of the $C P$-violating phase $\gamma$ in $B^{\pm}\to D K^{\pm}$ (and similar) decays , where the $D$ meson (which represents a superposition of $D^0$ and $\bar{D}^0$) subsequently decays to $\pi^+\pi^-\pi^+\pi^-$. To obtain optimal sensitivity to $\gamma$, the phase space of the $D \to \pi^+\pi^-\pi^+\pi^-$ decay is divided into bins based on a recent amplitude model of the decay. Although an amplitude model is used to define the bins, the measurements obtained are model-independent. The $CP$-even fraction of the $D \to \pi^+\pi^-\pi^+\pi^-$ decay is determined to be $F_{+}^{4\pi} = 0.769 \pm 0.021 \pm 0.010$, where the uncertainties are statistical and systematic, respectively. Using simulated $B^{\pm}\to D K^{\pm}, D \to \pi^+\pi^-\pi^+\pi^-$ decays, it is estimated that by the end of the current LHC run, the LHCb experiment could determine $\gamma$ from this decay mode with an uncertainty of $(\pm10\pm7)^\circ$, where the first uncertainty is statistical based on estimated LHCb event yields, and the second is due to the uncertainties on the parameters determined in this paper.


Introduction
A primary goal in modern flavour physics is to constrain the unitarity triangle (UT); an abstract representation of the famous Cabibbo-Kobayashi-Maskawa matrix that describes transitions between different quark flavours [1,2]. Key to determining the UT is are better experimental constraints on the angle γ (or φ 3 ), which is related to the phase difference between b → u W − and b → c W − quark transitions. Currently, γ is the least-well constrained angle of the UT, and can be determined, for example, using B − → DK − decays 1 , where D represents a superposition of D 0 and D 0 states [3][4][5][6][7][8].
The amplitudes B − → D 0 K − and B − → D 0 K − are overwhelmingly dominated by the tree-level transitions b → u cs and b → c us, respectively, and therefore offer an extremely clean method to measure γ. In order to obtain the necessary interference between B − → D 0 K − and B − → D 0 K − amplitudes, a final state f must be chosen that is accessible from both D 0 and D 0 , such as π + π − π + π − (4π ± ).
To determine γ in B − → DK − decays, one must know the relative magnitude and phase of D 0 → f and D 0 → f amplitudes, collectively known as the D → f hadronic parameters. The relative magnitudes can be determined by measuring D * + → D 0 π + decays that are subsequently followed by a D 0 → f decay; this is possible at a large variety of collider experiments, such as LHCb and the B-factories. Measuring the relative phase, however, is more challenging. One method is to infer the relative phase through use of an amplitude model; in principle this is the best way to exploit the available statistics, but theoretical uncertainties in determining the model can lead to large systematic uncertainties on γ. The relative phase can also be determined modelindependently by using samples of D → f decays, where the D meson is in a known superposition of D 0 and D 0 states. Previously, such data samples have been obtained from two sources: correlated DD pairs from the decay of a ψ(3770) meson [9][10][11][12][13][14][15][16] (the first charmonia resonance above the charm threshold); and the decay D * + → Dπ + , where the superposition of D 0 and D 0 states depends on the D meson decay-time [17][18][19]. In this paper we determine the relative magnitude and phase of D 0 → 4π ± and D 0 → 4π ± amplitudes using ψ(3770) decays collected by the CLEO-c experiment.
In multi-body D decays, such as D → 4π ± , there are infinitely many configurations of the final state momenta, each with a different amplitude. The parameter space that describes these final state configurations is known as the phase space of the decay. For the 4π ± final state, a phase space-integrated measurement was performed in Ref. [20] to determine the CP -content of the inclusive decay, and then applied in a B ± → DK ± study at LHCb [21]. However, to better exploit the information available in multi-body D decays, the phase space can be divided into bins such that regions of constructive and destructive interference do not dilute each other. Such a method has already been applied to the K 0 S π + π − final state [22] which gives the best single measurement of γ to date [23]; here an amplitude-model was used to group regions of phase space that have a similar phase difference between D 0 → K 0 S π + π − and D 0 → K 0 S π + π − amplitudes [24]. Recently an amplitude model for D → 4π ± has become available [25], so in this paper a similar technique is applied to the 4π ± final state. It is important to note that although the binning scheme is defined by an amplitude model, this will not result in any modeldependent bias. If the model is incorrect, this will just result in an increased statistical uncertainty.
This paper is organised as follows: Sec. 2 gives an overview of the formalism for correlated ψ(3770) → DD decays; Sec. 3 introduces the D 0 → 4π ± amplitude model that is used in Sec. 4 to inspire the phase space binning schemes; Sec. 5 discusses the dataset used in the analysis and the selection criteria applied; Sec. 6 describes the fit used to obtain constraints on the D → 4π ± hadronic parameters; Sec. 7 discusses the systematic uncertainties associated to the results in Sec. 8; Sec. 9 uses the measured hadronic parameters to estimate the γ constraints that are possible with current and future LHCb datasets; finally a summary is given in Sec. 10.

Formalism
The mass eigenstates of the D meson, |D 1,2 , can be written in terms of the flavour eigenstates, where the convention CP |D 0 = +|D 0 is followed such that |D 1 and |D 2 are the CP + and CPeigenstates, respectively. Throughout this paper CP violation in the D meson system is neglected, which is a good assumption given current experimental limits [26]. The masses and widths of D 1,2 are given by m 1,2 and Γ 1,2 respectively, which allows the average width, Γ D = 1 2 (Γ 1 + Γ 2 ), and the charm mixing parameters, x D = (m 1 − m 2 )/Γ D and y D = (Γ 1 − Γ 2 )/2Γ D , to be defined. Due to the effects of D-mixing, a D meson produced in a |D 0 eigenstate at t 0 = 0 evolves to an admixture of |D 0 and |D 0 states, denoted |D 0 (t) , after time t. Similarly, the |D 0 eigenstate evolves to |D 0 (t) .
The D 0 and D 0 decay amplitudes for a final state f are defined A f p = f p |H|D 0 andĀ f p = f p |H|D 0 , where H is the relevant Hamiltonian. The parameter p describes a point in the phase space of the D → f decay, and has a dimensionality that depends on the number of final state particles and their spin. For two-, three-and four-body pseudo-scalar final states the phase space dimensionality is 0, 2 and 5, respectively.
In this paper, the measured observables will always be integrated over bins of phase space. For the final states f and g, these regions are labeled by i and j, respectively 2 . The branching fraction for D 0 → f i and D 0 → f i decays are defined, where φ(p) gives the density of states at p. From these follow the quantities T f give the fraction of D 0 → f and D 0 → f decays that populate phase space bin i, respectively 3 . To describe the interference of D 0 → f and D 0 → f amplitudes integrated over the region i, the bin-averaged sine and cosine are defined, 3) where ∆δ f p = arg(A f p ) − arg(Ā f p ). Collectively, the parameters c f i , s f i , K f i andK f i are referred to as the hadronic parameters of the D → f decay.
Using the formalism above, the decay ψ(3770) → DD → f i g j is now considered. The strong decay ψ(3770) → DD results in a correlated DD pair in a C = −1 state. Therefore, |ψ(3770) → |D 0 D 0 − |D 0 D 0 . (2.5) Since the two D mesons evolve coherently, D-mixing has no observable consequences until one meson decays. Therefore, when studying such decays, what is important is the time difference, δt, between the D → f and D → g decays. The decay amplitude for ψ(3770) → DD → f p g q is given by [27], To obtain the decay rate, the magnitude of this amplitude is squared and integrated over the phase space regions i and j, and all decay-times. Expanding to second order in the small parameters x D and y D gives, This single formula is used to describe all decays studied in this paper. Note that Eq. 2.7 can be significantly simplified for some final states; for example CP eigenstates such as K + K − (CP + ) and K 0 S π 0 (CP -) have K g j ≡K g j , s g j ≡ 0 and c g j = η CP , where η CP = ±1 for CP + and CPeigenstates, respectively 4 .
When only one of the D meson final states is reconstructed it is known as a singletag. In this case, the final state g represents all possible D meson final states and K g ≡K g ≡ 1, s g ≡ 0 and c g = y D , leading to, The s g ≡ 0 can be understood by realising that for every final state g, there is a chargeconjugate final state g that has s g = −s g . The c g = y D can be understood by rewriting Eq. 2.3 as, (2.10) Therefore if g represents all final states, c g j = (Γ 1 − Γ 2 )/2Γ D = y D . Although the decay B ± → DK ± , D → f is not measured in this paper, it is important to consider its decay rate so that the D → 4π ± binning schemes defined in Sec. 4 give optimum sensitivity to γ in a future measurement of B ± → DK ± decays. The ratio of , where r B is that ratio of their magnitudes, and δ B is the strong phase difference. The B ± → DK ± , D → f p decay rates are then given by, (2.12) where x ± = r B cos(δ B ± γ) and y ± = r B sin(δ B ± γ). Integrating this expression over a phase space bin i then gives, An amplitude model is used to define how the five-dimensional phase space is divided into bins. Such a model has recently become available [25], which was determined from a fit to flavour tagged D 0 → 4π ± decays collected by the CLEO-c experiment.
To construct the total amplitude, the isobar approach was used, which assumes the decay can be factorised into consecutive two-body decay amplitudes. The dominant contributions to the model are D 0 → a 1 (1260) + π − , D 0 → σf 0 (1370) and D 0 → ρρ. In addition to the main ('nominal') model, Ref. [25] also includes a further 8 alternative models which use a different set of amplitude components -these are used for systematic studies.
Since CP conservation in D → 4π ± decays is assumed, the D 0 → 4π ± model implies Here p is the CP conjugate point of p, which has all charges reversed (C) and three-momenta flipped (P ). The assumption of CP conservation in D → 4π ± decays is explicitly tested in Ref. [25] by determining A f p and A f p independently from samples of D 0 and D 0 tagged decays, respectively. The results are consistent with the CP conservation hypothesis.

Binning
The definition of the 4π ± phase space bins strongly influences sensitivity to γ in B ± → DK ± , D → 4π ± decays. To best exploit the symmetries of the self-conjugate 4π ± final state, phase space bins are defined in pairs that map to each other under the CP operation. The bins are labeled such that bin +i is paired with bin −i, therefore, for any point p in +i, the CP conjugate point p will fall into bin −i. This choice of binning means that the following relations exist between the hadronic parameters of +i and −i bins: Since the relative magnitude and phase of A 4π p andĀ 4π p varies over the D → 4π ± phase space, so will the relative size of the interference term in B ± → DK ± , D → 4π ± decays. If a single bin contains regions of phase space with differing levels of interference (for example, constructive and destructive interference) the overall interference is diluted, and the sensitivity to γ is reduced. It is therefore preferable for both r 4π to be approximately constant within each bin. This is possible by using an amplitude model to assign each point in phase space a value of r 4π p and ∆δ 4π p , which are used to determine the bin number. Although a model is used to determine the bin number, this will not introduce any model-dependent systematic uncertainties, since the hadronic parameters will still be determined modelindependently. An incorrect model will only lead to a non-optimal binning, and an increased statistical uncertainty.
Before discussing the D → 4π ± binning scheme used in this paper, it is informative to review previous work on the final state K 0 S π + π − in Ref. [22]. This decay has a twodimensional phase space (the Dalitz plot) which can be parameterised by the variables . The region m 2 + > m 2 − is divided into N bins, labelled +1 to +N , which are reflected over the line m 2 + = m 2 − to obtain the −N to −1 bins (a reflection over this line is equivalent to CP ). Using the line m 2 + = m 2 − to divide the Dalitz plot is a good choice since most Cabibbo favoured (CF) amplitudes, such as D 0 → K * − π + , fall in the region m 2 + > m 2 − , whereas most Doubly Cabibbo suppressed (DCS) amplitudes fall into the region m 2 − > m 2 + . This is beneficial since it makes r K 0 S ππ p consistently large (small) over the +i (−i) bins. To determine the absolute bin numbers, the model prediction for ∆δ K 0 S ππ p is divided into 8 equal regions. The K 0 S π + π − binning scheme for N = 8 is shown in Fig. 1. The authors of Ref. [22] also provide a fine granularity lookup table that describes the binning shown in Fig. 1; this is very useful because the amplitude model is not necessary to reproduce the binning scheme. A similar idea will be used for the 4π ± binning schemes.
In order to remove the majority of this background, a K 0 S -veto bin is included in all D → 4π ± binning schemes that are later described in Sections 4.2 -4.5. The region of phase space that contains any π + π − pair satisfying 480 MeV < m(π + π − ) < 505 MeV is designated as the K 0 S -veto bin. Using the nominal D 0 → 4π ± amplitude model, the K 0 S -veto was found to remove approximately 10% of signal.

Equal / variable ∆δ 4π
p binning When comparing the K 0 S π + π − to the 4π ± final state, one clear difference is the decay amplitudes that contribute. As discussed, K 0 S π + π − has contributions from both CF and DCS amplitudes, whereas 4π ± only has contributions from singly Cabibbo suppressed (SCS) amplitudes. This means that there is no clear way to divide the phase space, like the line m 2 + = m 2 − in the K 0 S π + π − Dalitz plot. A different approach is therefore followed. The baseline amplitude model from Ref. [25] is used to assign each point p a value of ∆δ 4π p , then a bin number is assigned using, where δ 0 ≡ 0, δ N ≡ π and δ i < δ i+1 . This automatically fulfils the requirement that bin +i maps to bin −i under CP , since ∆δ 4π p ≡ −∆δ 4π p . The values of δ i are chosen using two methods: the equal ∆δ 4π p binning, for which δ i = iπ/N ; and the variable ∆δ 4π p binning, for which the values of δ i are chosen such that K 4π i +K 4π i is approximately the same in each bin.
Since amplitude models are difficult to reproduce, it is desirable to have a modelimplementation-independent binning scheme. This is possible by splitting the five dimensional phase space into many small hypervolumes, each of which is assigned a bin number. The overall bin is then formed from the combination of all hyper-volumes with that bin number. To create a model-implementation-independent binning scheme, referred to as a hyper-binning, a set of variables must be defined that parameterises the five-dimensional phase space of D → 4π ± decays. The variables {m + , m − , cos θ + , cos θ − , φ} are chosen, where m + (m − ) is the invariant mass of the π + π + (π − π − ) pair; θ + (θ − ) is the helicity angle of the π + π + (π − π − ) pair; and φ is the angle between the π + π + and π − π − decay planes (a full definition of these variables can be found in Appendix A). Since the hyper-binning is most easily implemented with square phase space boundaries, the following transformation is made, where m min is the minimum value kinematically possible for m + (or m − ). When using the variables {m + , m − , cos θ + , cos θ − , φ}, the kinematically allowed region of phase space is a hypervolume defined by the corners {m min , m min , −1, −1, −π} and {m max , m max , 1, 1, π}. This set of variables has been chosen to exploit the symmetries of the system, these being CP -conjugation and identical particle interchange: The symmetries for identical particle exchange allow the phase space to be 'folded' twice along the lines cos θ + = 0 and cos θ − = 0, reducing the phase space volume by a factor of four. A further folding is also possible by considering the CP operation; for a point p with bin number i, it follows that point p has bin number −i. An adaptive binning algorithm is used to create a hyper-binning scheme. At the beginning of the algorithm one hypervolume is defined with corners {m min , m min , 0, 0, 0} and {m max , m max , 1, 1, π}. At each iteration of the algorithm, the hypervolumes from the previous iteration are split in two, choosing to split in the dimension that has the fastest varying ∆δ 4π p , and picking a split point that is as close as possible to one of the bin boundaries defined in Eq. 4.1. The algorithm runs until either: splitting a hypervolume will always result in two hypervolumes with the same bin number; splitting a hypervolume will always result in a hypervolume that has an edge length narrower than the minimum allowed. Several minimum edge lengths were tested and the values {39 MeV, 39 MeV, 0.06, 0.06, 0.19 rad} were chosen since this results in a reasonable number of volumes (∼ 250, 000) while reproducing the parameters c i and s i to within 2% compared to a binning scheme that uses the model directly. It is possible to visualise the hyper-binning by taking two-dimensional slices of the five-dimensional phase space. Some examples are shown for the equal ∆δ 4π p binning with N = 5 in

Model predictions of the hadronic parameters
Using the integral expressions in Eqs. 2.2 -2.4 it is possible to calculate the hadronic parameters for a given amplitude model and binning scheme. This is done using the baseline and alternative amplitude models given in Ref. [25]. Since the baseline-model is used to determine the D → 4π ± binning schemes, using the hadronic parameters predicted with this model could result in a bias. Therefore, the arithmetic-mean of the hadronic parameters from all alternative models is used as the model prediction, and the covariance of the results is used to determine a model-uncertainty. To determine the statistical and systematic uncertainties, the hadronic parameters are calculated many times using the baseline model, each time varying the model parameters within their statistical and systematic uncertainties. The covariance of the results is used to determine a combined statistical and systematic uncertainty, which is added to the model-uncertainty in quadrature to obtain the total uncertainty. The model predictions for the equal / variable ∆δ 4π p binning are shown in Fig. 3.

Alternate binning
One drawback of the ∆δ 4π p binning schemes is that the variation of r 4π p across each bin is not considered, leading to K f i ∼K f i , as seen in Fig. 3. This means that the interference term in the B − → DK − decay rate, given in Eq. 2.13, is relatively small in all phase space bins. Ideally, one would choose to have r 4π p 1 in half of the phase space bins, enhancing the interference in these regions (and therefore the sensitivity to  γ). The r 4π p 1 condition is satisfied in the K 0 S π + π − final state, where many bins are dominated by DCS amplitudes. Although the SCS 4π ± final state has no clear line of symmetry that divides favoured from suppressed phase space regions, the amplitude model can be used to define such a split. Any point p that that satisfies r 4π p < 1 is assigned a bin number +i, whereas any satisfying r 4π p > 1 is assigned a bin number −i. The +i bin numbers are assigned using, which also uniquely defines the −i bin numbers i.e.
The same hypervolumes from the equal ∆δ 4π p binning schemes are used for the alternative binning schemes, but the bin number associated to each hypervolume is reassigned using Eq. 4.6. The model predictions for the alternative binning with N = 5 is shown in Fig. 3.

Optimal binning
To determine how sensitive a binning scheme is to a measurement of γ, the Q ± values are defined [24], 8) where N i B ± is the number of B ± → DK ± , D → f decays expected in bin i (Eq. 2.11 and Eq. 2.12), and Γ B ± (p) gives the differential decay rate (Eq. 2.13 and Eq. 2.14). The value of Q ± gives the statistical sensitivity on the parameters x ± and y ± from a binned analysis of B ± → DK ± , D → f decays, divided by the statistical sensitivity from an analysis with infinitely many bins. Substituting Eqs. 2.11 -2.14 into Eq. 4.8 gives, The Q value, Q 2 = 1 2 (Q 2 + + Q 2 − ), is then used to rank the sensitivity of different binning schemes to γ. The values δ B = 140 • , γ = 70 • and r B = 0.1 are used to determine Q. For the optimisation of the K 0 S π + π − binning schemes in Ref. [22], a simplified Q value was used where it was assumed r B = 0. Since the relative size of K f i andK f i does not need to be optimised for K 0 S π + π − (due to the division at m 2 + = m 2 − ), this assumption works well. For 4π ± decays, the simplified expression gives solutions where K f i ∼K f i , so the full expression is used instead.
An iterative algorithm is used to take any hyper-binning scheme (i.e. a collection of hypervolumes, each with a bin number, that span the D → 4π ± phase space) and reassign the bin numbers in order to maximise the model-prediction of Q. Each iteration of the algorithm involves looping over every hypervolume in the hyper-binning. For each hypervolume, every possible bin number (−N , ..., −1, +1, +N ) is assigned, and Q is recalculated; the bin number that gave the largest Q is then kept. The algorithm keeps running until no hypervolumes change their bin number, typically taking around 20 − 50 iterations.
Since the number of free parameters being optimised is so large, it is unavoidable that the optimisation procedure will fall into a local maximum. The outcome is therefore dependent on the starting values (i.e. the bin numbers assigned to each hypervolume). The starting bin numbers are therefore assigned using two methods: the equal ∆δ 4π p binning scheme (Eq. 4.1); and the alternate binning scheme (Eq. 4.6). The two sets of starting values give the 'optimal binning' and 'optimal-alternative binning', respectively. The set of hypervolumes used for the optimisation must have sufficient flexibility to describe the optimal binning. For all optimal binning schemes, the hypervolumes are first taken from the equal ∆δ 4π p binning scheme with N = 8, then further divided so that, for the sample sizes used in this paper, the probability of any single hypervolume being populated is less than 1/50.
After running the Q optimisation procedure it was found that occasionally the results had very small values of K f i +K f i for one or more bin pairs. For this reason a small change was made to the optimisation metric, is the lower threshold at which a constraint is applied to K f i +K f i . The Q value for the optimal and optimal-alternative binning schemes is shown in Fig. 4 for N = 1 − 8. Also shown are the Q values for the other binning schemes discussed in this paper.

Event Selection
The data set analysed consists of e + e − collisions produced by the Cornell Electron Storage Ring (CESR) at √ s = 3.77 GeV corresponding to an integrated luminosity of 818 pb −1 and collected with the CLEO-c detector. The CLEO-c detector is described in detail elsewhere [28][29][30][31]. Monte Carlo (MC) simulated samples of signal decays are used to estimate selection efficiencies. Possible background contributions are determined from a generic D 0 D 0 simulated sample corresponding to approximately fifteen times the integrated luminosity of the data set. The EVTGEN generator [32] is used to simulate the decays. The detector response is modelled using the GEANT software package [33]. Table 1 lists all D decay final states that are reconstructed in conjunction with a D → 4π ± decay, referred to as double-tagged decays. Underlined in Tab. 1 are the D decay final states that are also reconstructed alone, referred to as single-tagged decays. Unstable final state particles are reconstructed in the following decay modes: π 0 → γγ; Table 1.
List of all D decay final states that are reconstructed in conjunction with a D → 4π ± decay (double-tag modes). The underlined final states are also reconstructed alone (single-tag modes).
The selection procedure used for this paper is intended to be almost identical to that in Ref. [20]. The only change is to the selection criteria used to reject peaking background from D → K 0 S π + π − decays that are reconstructed as D → 4π ± ; henceforth referred to as K 0 S π + π − background. In Ref. [20] any π + π − pair with an invariant mass in the range [0.470, 0.530] GeV is required to have a reconstructed vertex that is compatible with the e + e − collision point. In this paper, any π + π − pair with an invariant mass in the range [0.480, 0.505] GeV is rejected, regardless of its compatibility with the e + e − collision point. The 4π ± phase space bins defined in Sec. 4 have the same region of phase space removed, so no corrections to the measured hadronic parameters are needed. In addition to the tags in Ref. [20], this analysis also uses the flavour-tags K ∓ e ± ν, and the quasi-flavour-tags K ∓ π ± , K ∓ π ± π 0 and K ∓ π ± π ∓ π ± . These decays are selected following the same criteria as Ref. [22].
The final states that do not include a neutrino or a K 0 L are fully reconstructed using the beam-constrained candidate mass, m bc ≡ s/(4c 4 Requirements are first placed on the value of ∆E, then m bc is used as the discriminating variable to distinguish signal from non-peaking backgrounds. For doubletags that are dominated by background from continuum production of light quarkantiquark pairs (π + π − , K + K − , π + π − π 0 and 4π ± ), the signal yield is determined using an unbinned maximum likelihood fit to the average m bc of the two D decays, m ave bc ≡ 1 2 (D 1 m bc + D 2 m bc ). The signal probability density function (PDF) is parameterised using the sum of a bifurcated Gaussian and a Gaussian, which have shape parameters fixed from a fit to samples of simulated signal decays 5 . The background PDF is parameterised using an Argus function [34]. Figure 5 shows an example of this fit for double-tagged π + π − π 0 candidates -the signal yield is determined in the m ave bc window [1.86, 1.87] GeV. For fully-reconstructed decays that are not continuum dominated, the double-tag yield is determined by counting events in signal and sideband regions of the two dimensional D 1 m bc vs. D 2 m bc plane, as indicated in Figure 5 for double-tagged K ∓ π ± π 0 candidates.  Figure 5.
(left) The m ave bc distribution of selected double-tagged π + π − π 0 candidates. Superimposed with a red line is the result of a unbinned maximum likelihood fit described in the text. The shaded purple area shows the background PDF, and the blue line shows the signal PDF. The vertical lines show the signal region. (right) Two dimensional D 1 m bc vs. D 2 m bc distribution of selected double tagged K ∓ π ± π 0 candidates. The square box covering the range 1.86 − 1.87 GeV shows the signal region, and the remaining boxes show the various sideband regions that are used to determine the combinatorial background contribution.
The final states containing a neutrino or a K 0 L cannot be fully-reconstructed; the energy and momentum, p miss and E miss , of the missing particle is inferred by using knowledge of the initial e + e − state and conservation of energy and momentum. The missing-mass squared, m 2 miss ≡ E 2 miss /c 4 − p 2 miss /c 2 , and the quantity U miss ≡ E miss − |p miss |c, are used to discriminate signal from background for decays involving a K 0 L or a neutrino, respectively. The double-tag yields are determined using an unbinned maximum likelihood fit to the discriminating variable, where the signal and background PDFs are taken from histograms of simulated data samples. Figure 6 shows an example of this fit for double-tagged K ∓ e ± ν and K 0 L π + π − candidates -the signal yields are determined within the signal windows indicated.  Figure 6. Distribution of U miss (m 2 miss ) for selected double tags containing a neutrino (K 0 L ). Superimposed is the result of an unbinned maximum likelihood fit that is described in the text. The blue/green/red shaded area shows the distribution of combinatoric/continuum/peaking background respectively. The red line shows the total signal + background PDF. The black vertical lines indicate the events that fall within the signal region that are used for further analysis.
The dominant peaking background contribution to all double-tags is from K 0 S π + π − background, which is estimated from the generic MC sample of DD events, and typically constitutes about 5 − 10% of the selected events. A data-driven estimate of this background is also calculated using the events that are rejected by the π + π − mass cut -this shows good agreement with the estimates from generic MC. All decays involving a K 0 L decay have a peaking background from the equivalent decay with a K 0 S instead of a K 0 L -these are referred to as cross-feed backgrounds. Using the simulated samples of D → K 0 S X decays it is possible to find the ratio of D → K 0 S X decays that are incorrectly reconstructed as D → K 0 L X to those correctly reconstructed as D → K 0 S X. Since for every D → K 0 L X decay considered in this paper, the equivalent D → K 0 S X decay is also considered, this allows the background to be estimated using the measured D → K 0 S X yields. The decay π + π − π 0 has a peaking background from K 0 S π 0 that is largely suppressed by requiring the π + π − vertex to be consistent with with e + e − collision point. Since the decay K 0 S π 0 is also considered in this paper, the K 0 S π 0 signal yield can be used (in the same manner as for the cross-feed backgrounds) to estimate the background contribution. All remaining peaking backgrounds are either negligible, or considered in the systematics uncertainties in Sec. 7.
Single-tagged candidates are selected using identical criteria to the corresponding double tags, with the exception of π + π − , K + K − and K ∓ π ± decays that have additional cuts to veto cosmic ray muon and radiative Bhabha events [35]. The number of singletags is estimated from a fit to the m bc distribution. The signal and background PDFs are the same as those used in the fit to the m ave bc distribution of continuum dominated double-tags. The signal shape parameters are fixed from a sample of simulated signal decays. Figure 7 shows an example of this fit for single-tagged π + π − π 0 candidates -the signal yield is determined in the signal region indicated. Following Ref. [35], a further uncertainty is assigned to each of the single-tag yields to account for any mismodelling of the signal PDF. For final states with no electromagnetically neutral final state particles (K + K − , π + π − , K ∓ π ± ) the uncertainty assigned is 2.0% of the measured signal yield. For final states where the neutrals are relatively hard (K 0 S π 0 , K 0 S η(γγ)) or soft (all other modes), uncertainties of 2.5% and 5.0% are assigned, respectively. In events where more than one single-or double-tagged candidate is reconstructed, an algorithm is used to select a single candidate based on information provided by the ∆E and m bc variables. The particular choice of metric varies depending on the category of double-tag, and is optimised through simulation studies.
For double-tagged decays, the signal yields are evaluated in bins of 4π ± , K 0 S π + π − and K 0 L π + π − phase space. For these final states, the four-momenta of the D daughters Decay Mode are determined with a constraint on the previously measured D 0 mass [36], ensuring that all signal candidates fall within the kinematically allowed region of phase space. The 4π ± final state is binned using the schemes in Sec. 4. The K 0 S π + π − and K 0 L π + π − final states are binned according to the 'Equal ∆δ D BABAR 2008' scheme from Ref. [22], which is shown in Fig. 1. For non-continuum dominated decays, the binned yields are determined by counting the number of candidates in the signal region of the D 1 m bc vs. D 2 m bc plane -the background estimates are discussed in Sec. 6. For continuumdominated final states a fit the m ave bc distribution is performed in each phase space bin. In the case where 2 or more phase space bins have an identical decay rate (e.g. 4π ± vs. π + π − π 0 has the same decay rate in bin +i and −i) they are merged before determining the signal yield. The samples of flavour and quasi-flavour double-tags are split using the charge of the kaon before the binned yields are determined. The phase space-integrated background subtracted event yields for all single-and double-tagged decays are given in Tab. 2.

Fit for 4π ± hadronic parameters
This section describes the fitting algorithm used to determine constraints on the D → 4π ± hadronic parameters. Following from Eq. 2.7, the expected number of ψ(3770) → DD → f i g j signal decays is given by, where N DD is the total number of ψ(3770) → DD decays in the data sample. In the literature, different parameterisations of the hadronic parameters are used for different categories of final state, which sometimes differ from K f i ,K f i , c f i and s f i parameterisation used to derive the formalism in this paper. The different parameterisations used are summarised in Tab. 3, which are used as free parameters in the fit for the relevant final states. The new parameters introduced are: the CP -even fraction, F f + ; the coherence factor, R f D ; the average strong phase difference, δ f D ; and the ratio of The relationship between these and the K Substituting the various parameterisations in Tab. 3 into Eq. 6.1, it is clear that different categories of tag provide sensitivity to different hadronic parameters. The flavour and quasi-flavour tags give sensitivity to K f i andK f i ; the CP tags and π + π − π 0 tags give sensitivity to K f i ,K f i , and c f i ; and the K 0 S π + π − and K 0 L π + π − tags give sensitivity to all hadronic parameters.
The expected efficiency and background corrected yield is given by, where f i g j is the reconstruction and selection efficiency for the decay in question, and N bkg f i g j is the expected number of background. The quantity f i g j is determined from large samples of simulated signal decays, correcting for known discrepancies between data and simulation. Before efficiencies are calculated, the simulated samples containing K 0 S π + π − and 4π ± decays are reweighted to their model expectations (using the D 0 → K 0 S π + π − BABAR model [37] and the nominal D 0 → 4π ± model [25]) including the effects of quantum correlations. The simulated sample of K 0 L π + π − decays is also reweighted to ; this approximation holds in the scenario that only CF and DCS amplitudes contribute, and the two do not overlap in the Dalitz plot. A systematic uncertainty is later assigned to account for any model dependence in the efficiency determination.
The total background estimate is broken down into the following expression, f g are the total number of K 0 S π + π − and combinatoric background in the DD → f g decay, respectively. The quantities κ and κ flat f i g j give the fraction of background that falls into the phase space bins i and j. The final term, N sig gives the number of cross-feed background from the decay DD → f i h j . The quantity f h g gives the fraction of DD → f i h j decays that are incorrectly reconstructed as DD → f i g j , to those correctly reconstructed. The value of N K 0 S ππ f g is taken from generic MC, as was used for the determination of the of background subtracted yields in Tab. 2. The value of κ K 0 S ππ f i g j is found using a large sample of simulated D → K 0 S π + π − decays that are reconstructed as D → 4π ± . Before calculating κ K 0 S ππ f i g j , the simulated sample is first reweighted to the model expectation, based on the phase space location of the generated K 0 S π + π − decay, and including quantum correlations. Since the K 0 S π + π − model has been shown to give good agreement with model-independent measurements [22], any model dependent bias should be small, but this is considered as a systematic uncertainty later. The value of N flat f g is determined from the sideband regions, as described in Sec. 5. For continuum-dominated and single-tagged decays N flat f g = 0, since the signal yields are determined from a fit to m ave bc , so already have the combinatoric background component subtracted. The value of κ flat f i g j is determined using simulated signal decays distributed according to the density of states (phase space). Where possible, this assumption is checked using the sideband regions, which shows good agreement. Systematic uncertainties are assigned to cover any bias from this assumption.
The values of the 4π ± hadronic parameters are obtained by maximising the loglikelihood, log L. The Poisson distribution, P (k; λ) ≡ λ k e −λ /k!, gives the probability of observing k events when λ are expected. For double-tagged decays that are not continuum-dominated the log-likelihood receives a term, where M f i g j is the number of events counted in the signal region of the decay DD → f i g j . For continuum-dominated double-tags and single-tags, the signal yield is obtained from a fit, which has an associated uncertainty σ f i g j . Therefore, the log-likelihood receives a term, where G(k; µ, σ) is a Gaussian distribution with mean µ and width σ. External inputs are needed to constrain various parameters in the fit. For the partially-reconstructed CP final states K 0 L π 0 and K 0 L ω it is not possible to obtain a single-tagged sample, which would provide the fitter with constraints on the product N DD ×BF(D 0 → f ). This constraint is important for normalising the respective doubletag yield, so an alternative method is followed for the K 0 L π 0 and K 0 L ω final states. In order to constrain N DD , the single-tagged K ∓ π ± yield is measured in conjunction with an external constraint on BF( [36]. For the quasi-flavour tags, external constraints are provided for the hadronic parameters r f D , R f D and δ f D , which are taken from Ref. [26] and Ref. [38]. The selfconjugate final state π + π − π 0 is not a CP eigenstate, so its CP -even fraction, F πππ 0 + , is constrained to its previously measured value [20]. The charm-mixing parameters are constrained to their world-average values [26]. The central values and uncertainties of the constraints are listed in Tab. 4. All constraints are applied by including a Gaussian constraint, similar to Eq. 6.5, in the log L; where available, correlations between the parameters are also included.
The hadronic parameters of the K 0 S π + π − and K 0 L π + π − final states are also constrained. The parameters c are constrained using the covariance matrix for the BABAR equal ∆δ D binning given in Ref. [22]. An adjustment must be made to the constraints on c L ππ i are taken directly from Ref. [39]; since it is the parameters K

Constraint
Source Ref. [20] * The constraint on BF(D 0 → K 0 L ω(3π)) is taken from BF(D 0 → K 0 S ω(3π)) with a systematic uncertainty of 20%. † Ref. [26] uses the convention CP |D 0 = −|D 0 , so the transformation δ K − π + → δ K − π + + π is applied. Table 4. List constraints used in the analysis. The right hand column gives the source of the constraint, along with any conventional adjustments that have to be made to use them in this analysis.
predictions [40,41], as determined in Ref. [42]. Since the amplitude models are fit to D decay-time integrated samples of D * + → D 0 π + decays, small corrections must be made for D-mixing using the expression [19], where B K 0 S ππ i is the fraction of D * + → D 0 π + , D 0 → K 0 S π + π − decays in phase space bin i. Using external inputs from Refs. [22,26], the system of equations is solved to find F In principle, the normalisation parameter N DD can be shared for every decay mode considered in the analysis, since the same e + e − collision data are used. In reality, however, this is not always desirable since the estimation of N DD relies on the absolute efficiencies (rather than the relative, bin-to-bin, efficiencies) determined from simulated samples. For the double-tagged samples of K 0 S π + π − , K 0 L π + π − , K ∓ π ± π 0 , K ∓ π ± π ∓ π ± , and K ∓ e ± ν decays, almost all information comes from the relative bin-to-bin yields, so sharing a normalisation parameter provides little benefit while introducing a potential source of systematic uncertainty. Therefore, these final states each have their own normalisation parameter, N f DD , in the fit. On the other hand, the double-and single-tagged K ∓ π ± samples share a normalisation constant, which allows the fitter to constrain BF(D 0 → 4π ± ), since BF(D 0 → K − π + ) has an external constraint (Tab. 4). This normalisation constant is also shared with all single-and double-tagged CP and π + π − π 0 final states.
The log L expression is maximised numerically using the MINUIT software [43]. The maximisation procedure is repeated 5 times with different starting values to ensure the global maximum of log L has been found (as opposed to a local maximum). Statistical uncertainties and correlations between fit parameters are provided by Minuit from evaluating the second derivatives of log L with respect to the fit parameters.
The fitting procedure is tested using 400 simulated experiments that use the background and efficiency estimates from the fit to data. The D → 4π ± hadronic parameters used to generated the pseudo-experiments are taken from model predictions. The hadronic parameters of other final states are taken from their previously measured values, and randomly sampled from their associated uncertainties. No statistically significant bias was found in the fit procedure.
The central values and statistical uncertainties of the D → 4π ± hadronic parameters from the fit to data are given in Tab. 7, and the statistical correlations in Appendix B. In this paper only results using 4π ± binning schemes with N = 5 are presented, although the results for N = 1 − 5 can be found in the supplementary material.

Systematics
The systematic uncertainties on the 4π ± hadronic parameters are broken down into several components, as listed in the systematic uncertainty breakdown in Tab. 5. Each of these components will be discussed in the following.
Bin migration Due to the finite detector resolution, it is possible for an event occurring in one phase space bin to be reconstructed in another; this bin-migration is relevant to the 4π ± , K 0 S π + π − and K 0 L π + π − final states. Since decays to these final states do not proceed by any narrow resonances, bin migration is not expected to significantly bias the result. Using samples of simulated signal events (that are reweighted to their model expectations), a migration matrix is calculated, whose elements M ik give the probably of an event generated in bin k to be reconstructed in bin i. For the fullyreconstructed final states 4π ± and K 0 S π + π − , the diagonal elements of M ik are typically ∼ 95%, whereas for the partially reconstructed K 0 L π + π − final state they are ∼ 85%. The migration matrices are used in the calculation of the expected yields (Eq. 6.2) and the fit is rerun. The absolute difference between this result and the nominal result, which is obtained without correcting for bin-migration, is taken as a systematic uncertainty.
Multiple candidate selection To check that the multiple candidate selection (MCS) procedure does not bias the result, an alternative MCS procedure is followed where one candidate is chosen at random (rather than based on a metric). The difference between the hadronic parameters determined using this selection and the nominal selection is taken as a systematic uncertainty.
Relative efficiencies In the nominal fit, the relative efficiency between phase space bins is determined using simulated signal samples that are reweighted to their model expectation. To estimate an upper limit on the systematic uncertainty introduced by the model uncertainty, the efficiency estimates are redetermined with the simulated samples reweighted to phase space. The absolute difference between the result using alternative efficiency estimates and the nominal result is taken as a systematic uncertainty.
Relative K 0 S π + π − background distribution To determine the relative distribution of K 0 S π + π − background, a sample of simulated D → K 0 S π + π − decays, reconstructed as D → 4π ± , is reweighted to its model prediction, including quantum correlations. In order to determine a systematic uncertainty, the quantum correlations are neglected (equivalent to setting c f i = s f i = 0 in Eq. 6.1) and the K 0 S π + π − background distribution is recalculated. The absolute difference between this result and the nominal result is taken as a systematic uncertainty.
Absolute K 0 S π + π − background yields In the nominal fit, the total number of K 0 S π + π − background events are estimated using the generic sample of simulated data. Alternatively, this is determined using a data-driven technique. The relative event numbers in the K 0 S veto region and the signal region is determined from simulation for both K 0 S π + π − background and 4π ± signal. These numbers are used to estimate the K 0 S π + π − background contamination in the signal region based on the observed number of events in the K 0 S veto region. The fit is rerun with the alternative K 0 S π + π − background yields, and the absolute difference between this result and the nominal result is taken as a systematic uncertainty.
Relative flat background distribution The relative number of combinatorial background events across phase space bins is assumed to be distributed according to phase space. As an alternative method, the relative numbers are taken from the combinatorial background events in the generic MC sample. The absolute difference between this result and the nominal result is taken as a systematic uncertainty.
Absolute flat background yields For fully-reconstructed non-continuum dominated double-tagged decays, the total number of combinatorial background is estimated from the number of events in five sideband regions of the two dimensional D 1 m bc vs. D 2 m bc plane (see Fig. 5). Each sideband region is associated with a particular background type, which is assumed to have the same density in the sideband and signal regions. Alternatively, the relative density of background between the sideband and signal regions is taken from generic MC. The alternative background estimates are used in the fit, and the difference between this result and the nominal result is taken as a systematic uncertainty.
For partially-reconstructed double-tagged decays, the total number of combinatorial background is determined from a fit to the m 2 miss or U miss distribution (see Fig. 6). Alternatively, the combinatorial background yield is determined using a simpler sideband-subtraction approach. The alternative background estimates are used in the fit, and the difference between this result and the nominal result is taken as a systematic uncertainty.
Continuum dominated signal yields For continuum-dominated double-tagged decays, the signal yield in each phase space bin is determined from a fit to the m ave bc distribution. The fits are repeated with an alternative signal (sum of a Johnson function [44] and a Gaussian) and background (second order polynomial) parameterisation, in a reduced m ave bc range. The alternative signal yields and uncertainties are used to determine the hadronic parameters, and the difference between this result and the nominal result is taken as a systematic uncertainty.
Non-resonant dilution The final states K 0 S ω and K 0 S η(π + π − π 0 ) have small contributions from non-resonant K 0 S π + π − π 0 decays, which are estimated from generic MC. Since this background contributes to both the single-tagged and double-tagged modes, it can be accounted for by making a small adjustment to the CP -even fraction of each final state, which would be identically zero (CP -odd) in the case of no background. Since the CP -content of this background is not known, it is conservatively assumed to be CP -even. The fit is rerun with the updated CP -even fractions, and the difference between this result and the nominal result is taken as a systematic uncertainty.
Simulated sample statistics In the nominal fit, the background and efficiency estimates all have an uncertainty due to limited statistics in simulated data samples. The fit is rerun twenty times, each time randomly varying the efficiency and background estimates within their uncertainties. The covariance of the results obtained is used to determine a systematic uncertainty.
A breakdown of the systematic uncertainties for the optimal alternative binning with N = 5 is given in Tab. 5. The largest systematic uncertainty comes from imperfect knowledge of the combinatorial background. The total systematic uncertainties for all binning schemes with N = 5 are given in Tab. 7, and the systematic correlations in Appendix B. The equivalent information for the other binning schemes considered is provided in the supplementary material. For all parameters the total uncertainty is statistically dominated.  Table 5. A breakdown of the systematic uncertainties for the optimal alternative binning scheme with N = 5.
Binning scheme  Table 6. The compatibility of the measured 4π ± hadronic parameters with the model predictions for all binning schemes with N = 5.

Results and consistency checks
The measurement of the 4π ± hadronic parameters with statistical and systematic uncertainties is given in Tab. 7, with correlations in Appendix B. The results are compared to the model predictions in Fig. 8 and Fig. 9. The compatibility between the results and the model predictions is quantified by calculating the χ 2 between to two, where all correlations are included. This is done independently for the c 4π i /s 4π i , and T 4π i /T 4π i parameters, and for the combination, with the results in Tab. 6. The parameters T 4π i andT 4π i show good agreement with the model predictions, which is expected since the model was determined from a fit to D 0 and D 0 tagged data. The parameters c f i and s f i are in slight tension with the model predictions, with p-values ranging from 0.03 to 0.18, but they clearly follow the same general trend in the c f i -s f i plane. It is worth repeating here that any incompatibility with the model will not introduce additional systematic uncertainties to a measurement of γ, but will only increase the statistical uncertainty.
Using the measured 4π ± hadronic parameters, the CP -even fraction of all phase space bins,F 4π + , is calculated using the formula, where the tilde indicates that a π + π − mass window is excluded from the D → 4π ± phase space i.e. F 4π + represents the CP -even fraction for the entire phase space. The values ofF 4π + are presented in Tab. 7, and are consistent among binning schemes. The nominal model is used to determine F 4π + −F 4π + = −0.002 ± 0.002 which can be used as a correction factor to determine F 4π + from the values ofF 4π + in Tab. 7. The Q value of each binning scheme is determined using Eq. 4.8, and presented in Tab. 7; as expected, the optimal binning schemes give the largest Q values. The Q value for a single phase space bin is calculated, usingF 4π + , to be 0.505. Therefore, based on the relative Q values, and using the optimal-alternative binning scheme with N = 5, the increase in statistical power for a measurement of γ is increased by ∼ 2.2 times with respect to the phase space integrated case. 6 The consistency of the c 4π i and s 4π i constraints obtained using different categories of final state is shown in Fig. 10 and Fig. 11 for the optimal alternative binning scheme with N = 5. For Fig. 10, each fit to one of the five categories (CP + , CP -, π + π − π 0 , K 0 S π + π − and K 0 L π + π − ) uses all flavour and quasi-flavour tags. The constraints obtained are consistent between all categories of final state.
The fit is also run using a single 4π ± phase space bin, which givesF 4π + = 0.760 ± 0.021 ± 0.021. The consistency of this result is checked between all final states in Fig. 12, following a similar method to the one used to obtain Fig. 10. Good consistency is observed.
As a 'default' binning scheme, we take the optimal-alternative binning with N = 5, as this has highest predicted Q value. The default binning scheme also has the largest measured Q value, although this information was not used to pick the default binning since it could bias the results. The value ofF 4π + determined using the default model is 0.771 ± 0.021 ± 0.010, which leads to F 4π + = 0.769 ± 0.021 ± 0.010 ± 0.002, where the final uncertainty is due to the K 0 S veto. The value of F 4π + is an important input for determining the total CP content of the neutral D meson, which is related to the charm mixing parameter y D through Eq. 2.10 [45].  Table 7. The hadronic parameters measured for each of the 4π ± binning schemes discussed in Sec. 4 where N = 5. The first uncertainty given is statistical, and the second systematic. Also given is the CP-even fraction,F 4π + , and the Q value, defined in Sec. 4; the uncertainties on these parameters are propagated from the statistical and systematic uncertainties on the hadronic parameters.   Figure 10. Constraints on the c 4π i and s 4π i parameters using the optimal alternative binning scheme with N = 5, determined using different subsets of tags. The grey bands show the model predictions and uncertainties. The red lines show the measured values and uncertainties when using a single subset of tags -the inner error bar shows the statistical uncertainty, and the outer error bar shows the combined statistical and systematic uncertainty. The yellow band shows the combined result using all subsets of tags -the lighter shade of yellow represents the statistical uncertainty, and the darker shade of yellow shows the combined statistical and systematic uncertainties.  Figure 12. The CP -even fraction over all phase space bins,F 4π + , determined using different subsets of tags. The grey bands show the model predictions and uncertainties. The red/grey lines show the measured values and uncertainties when using a single subset of tags -the inner error bar shows the statistical uncertainty, and the outer error bar shows the combined statistical and systematic uncertainty. The yellow band shows the combined result using all subsets of tags -the lighter shade of yellow represents the statistical uncertainty, and the darker shade of yellow shows the combined statistical and systematic uncertainties.

Sensitivity studies
In this section the measured 4π ± hadronic parameters from Sec. 8 are used to simulate B ± → DK ± , D → 4π ± datasets, which in turn are used to estimate the sensitivity to γ. Three scenarios with different event yields are studied, based on measured and extrapolated B ± → DK ± , D → 4π ± event yields from LHCb: "LHCb run I", with event yields of ∼ 1, 500 already recorded by LHCb with 3 fb −1 [21] of data; "LHCb run II", with plausible event yields of ∼ 6, 500 at the end of the next LHC data taking period with approximately twice the collision energy and an estimated 8 fb −1 of data; and "LHCb phase 1 upgrade", with plausible event yields of ∼ 100, 000 after phase 1 of the LHCb upgrade. The increase in the heavy flavour cross section at higher collision energies is accounted for, along with the expected improvement in trigger efficiency at the LHCb phase 1 upgrade [46]. The extrapolations have of course large uncertainties. The presence of background and systematic effects has been neglected in these studies, which is a reasonable assumption given previous measurements [21].
Toy datasets of B ± → DK ± , D → 4π ± decays are generated using Eq. 2.13 and Eq. 2.14 with δ B = 140 • , γ = 70 • and r B = 0.1. For each toy dataset, the central values of the 4π ± hadronic parameters are randomly sampled from the measured values and uncertainties. When fitting the toy datasets, the parameters δ B , γ, r B and an overall normalisation parameter are allowed to float, whereas the 4π ± hadronic parameters are fixed to their measured values. Therefore, the uncertainties obtained from the fit only account for the finite B ± → DK ± , D → 4π ± statistics, σ stat . The uncertainties on the parameters c 4π i and s 4π i are propagated to δ B , γ and r B by repeating the fit 200 times, where for each fit c 4π i and s 4π i are randomly sampled from their associated covariance matrix. The covariance of the values obtained is used to assign an uncertainty, σ had . The parameters K 4π i andK 4π i can be determined to an arbitrarily high precision at LHCb using D * + → D 0 π + decays, so the uncertainties on these parameters are neglected. As an alternative approach, the c 4π i and s 4π i parameters are Gaussian constrained in the fit, but this method was found to give a heavily biased estimate of γ, up to 70% of the statistical uncertainty. The nominal fit method gives good coverage and small biases of less than 10%.
The expected γ uncertainties are presented in Tab. 8 for several binning schemes. For each case the expected γ uncertainty is median uncertainty determined from fits to 100 simulated datasets. For each binning scheme type the uncertainty on γ generally decreases with increasing numbers of bins -for illustration, the uncertainty on γ is shown for the optimal-alternative binning scheme for N = 2 − 5 in Tab. 8. The γ uncertainties are also compared between different binning scheme types with N = 5; all result in similar values of σ stat (γ), although the values of σ had (γ) are notably larger for the 'variable ∆δ 4π p ' and the 'alternative' binning schemes. This is likely due to the measured central values of the s 4π i parameters being consistent with zero for these schemes. For the default binning (optimal alternative with N = 5) the expected uncertainties are (18 ⊕ 13) • , (10 ⊕ 7) • and (2.5 ⊕ 4.4) • for the LHCb "Run I", "Run II" and "Phase 1 Upgrade" scenarios, respectively, where the uncertainties are given in the form σ stat (γ) ⊕ σ had (γ).
Since σ had (γ) ≈ σ stat (γ) for the LHCb "Run I" and "Run II" scenarios, and σ had (γ) > σ stat (γ) for the "Phase 1 Upgrade" scenario, it is interesting to consider the impact that BESIII could have on reducing σ had (γ). Currently BESIII have collected 2.9 fb −1 of e + e − collisions at the ψ(3770) resonance, and a further ∼ 7 fb −1 is planned for the future. These datasets correspond to approximately 3.5 and 12 times the amount collected by CLEO-c, respectively. It is assumed that the uncertainties on the 4π ± hadronic parameters would be reduced by 1/ √ 3.5 and 1/ √ 12, respectively, compared to the constraints obtained in Sec. 8. The central values of the estimated BESIII measurements are different for each simulated dataset, and are randomly sampled from the constraints obtained in Sec. 8. Fig. 13 shows, for the default binning scheme, the expected values of σ had (γ) for different numbers of B ± → DK ± , D → 4π ± decays. This is shown for the hadronic parameter constraints measured in this paper, and the expected constraints for the two BESIII data taking periods. With 10.0 fb −1 of BESIII data, the expected γ uncertainties become (18 ⊕ 3) • , (10 ⊕ 1.7) • and (2.5 ⊕ 1.2) • for the LHCb "Run I", "Run II" and "Phase 1 Upgrade" scenarios, respectively. It is also possible that BESIII could make further gains in sensitivity by using additional numbers of phase space bins. Improved constraints on the 4π ± hadronic parameters could be obtained using D-mixing, as has been done for the K ± π ∓ π ± π ∓ final state in Ref. [19]; this would require further investigation. 10.0 ⊕ 7.9 2.6 ⊕ 5.0 Table 8. Expected γ sensitivity determined from simulated samples of B ± → DK ± , D → 4π ± decays for a variety of D → 4π ± binning schemes. Details of the simulation and fitting procedure can be found in the text. The uncertainties are given for three different data taking periods of the LHCb experiment, where the number of signal decays in each case it taken/extrapolated from existing measurements. The uncertainty on γ comes from two sources: the uncertainty due to to limited B ± → DK ± , D → 4π ± statistics, σ stat (γ); and the uncertainty due to limited knowledge of the 4π ± hadronic parameters that are measured in this paper, σ had (γ). Both uncertainties are shown in the table, and are given in the format (σ stat (γ) ⊕ σ had (γ)). All expected uncertainties are the median uncertainty from 100 simulated experiments.  Figure 13. Expected γ uncertainties obtained using different numbers of B ± → DK ± , D → 4π ± decays and the default binning scheme. The black line shows the estimated uncertainty due to limited B ± → DK ± , D → 4π ± statistics. The red, green and blue lines shows the estimated uncertainty due to the measured/predicted constraints on the 4π ± hadronic parameters from CLEO-c data with 0.818 fb −1 , BES III with 2.9 fb −1 , and BES III with 10.0 fb −1 , respectively. The grey bands highlight the event numbers that correspond to different LHCb data taking periods.
Using 818 pb −1 of e + e − collision data collected by the CLEO-c detector, the hadronic parameters of the D → 4π ± decay are measured in bins of phase space for the first time. This allows the UT angle γ to be determined using only B ± → DK ± decays where D decays to the 4π ± final state; previously only phase space integrated measurements have been possible [20,21], which need to be combined with other final states to obtain constraints on γ [21,47].
The phase space of the D → 4π ± decay is divided into bins based on the nominal amplitude model from Ref. [25]. The equal and variable ∆δ 4π p binning schemes are based on an equal/variable division of ∆δ 4π p , whereas the alternate binning scheme also uses the relative magnitude of D 0 → 4π ± to D 0 → 4π ± amplitudes. The optimal and optimal alternative binning schemes are defined to optimise the expected sensitivity to γ in B ± → DK ± decays. Although an amplitude model is used to inspire the binning schemes, the results are model-unbiased; any modelling deficiencies will only result in an increased statistical uncertainty on γ.
Since amplitude models can be notoriously difficult to reproduce, it is useful to have a model-implementation independent method to represent a binning scheme. The phase space of the D → 4π ± decay is five-dimensional, so using traditional techniques to divide the phase space into N 5 equally sized hypervolumes, where each is assigned a bin-number, would result in an unmanageable number of hypervolumes. An adaptive binning scheme is developed that uses an array of differently sized hypervolumes to drastically reduce the total number of hypervolumes needed, typically around 250, 000.
The measured values of the hadronic parameters are compared to the modelpredictions, which show good agreement for the parameters T 4π i andT 4π i , but a slight tension for c 4π i and s 4π i . This could either be due to statistical fluctuations, which could be tested with larger datasets at BESIII, or a possible residual mismodelling of the phase motion across the D → 4π ± phase space in Ref. [25].
The consistency of the results is checked using different subsets of final states, which give statistically compatible results. The CP even fraction over all phase space bins,F 4π + , is observed to be consistent between all binning schemes. Using the 'default' binning scheme, F 4π + is determined as 0.769 ± 0.021 ± 0.010 ± 0.002 where the uncertainties are statistical, systematic, and from the K 0 S veto, respectively. This is the most precise determination F 4π + to date. Using the 4π ± hadronic parameters measured in this paper, samples of B ± → DK ± , D → 4π ± decays are simulated, then used to estimate the potential sensitivity to γ. It is shown that, using estimated sample sizes from LHCb at the end of its current running period ("Run II") and the hadronic parameter constraints from this paper, constraints of σ(γ) = (10 ⊕ 7) • could be obtained, potentially making 4π ± one of the most sensitive final states for a measurement of γ. The first uncertainty is due to limited B ± → DK ± statistics, and the second is due to uncertainties on the 4π ± hadronic parameters. It is shown that the latter uncertainty could be reduced to around 1.7 • by using current and future BESIII datasets.

Acknowledgements
This analysis was performed using CLEO-c data. The authors, some of whom were members of CLEO, are grateful to the collaboration for the privilege of using these data. We also thank Guy Wilkinson, Sneha Malde, Jim Libby, Tim Gershon, Philippe D'Argent, Kostas Petridis and Vincenzo Vagnoni for their valuable feedback on the paper draft. We also acknowledge the support of the UK Science and Technology Facilities Council (STFC), and the European Research Council's support under Framework 7 / ERC Grant Agreement number 307737.

B Statistical and systematic correlations
Equal ∆δ 4π p binning statistical correlations Equal ∆δ 4π p binning systematic correlations  Table 9. The statistical and systematic correlations between the 4π ± hadronic parameters using the equal ∆δ 4π p binning scheme with N = 5.
Variable ∆δ 4π p binning statistical correlations Variable ∆δ 4π p binning systematic correlations  Table 10. The statistical and systematic correlations between the 4π ± hadronic parameters using the Variable ∆δ 4π p binning scheme with N = 5.
Alternative binning statistical correlations Alternative binning systematic correlations  Table 11. The statistical and systematic correlations between the 4π ± hadronic parameters using the Alternative binning scheme with N = 5. Optimal binning systematic correlations  Table 12. The statistical and systematic correlations between the 4π ± hadronic parameters using the Optimal binning scheme with N = 5.

C.1 List of files
The supplementary material can be found at Ref. [48]. The directory structure is organised so that each phase space binning scheme has its own directory. Each of these directories has the same file structure inside, which is described in Tab. 14. Additionally there is a Root macro loadresults.C, and a collection of C ++ functions in usehypbinning.cpp that can be used to load the supplementary material files that are in Root format. All results are additionally given in text format for greater flexibility.

C.2 Hyper-binning
For flexibility, the hyper-binning schemes are given in three different formats in the supplementary material, which will be discussed in this section. All binning schemes have been produced with a D 0 mass of 1864.84 MeV, and a π ± mass of 139.57 MeV; this defines the boundaries of the m + and m − variables.
It is recommended to use the Root format (hypbinning.root), which can be loaded using the HyperPlot C ++ package located at, http://samharnew.github.io/HyperPlot/index.html, using the HyperHistogram class. An example C ++ function is given in usehypbinning.cpp that can be compiled with the HyperPlot package to load any of the hyper-binning schemes.
The compressed directory hypbinning.zip contains two text files; hypbinning.txt and hypbinningwlinks.txt. Implementing the hyper-binning using the information in hypbinning.txt is significantly easier than hypbinningwlinks.txt, but the resulting code will be up to 10, 000 times slower (although this may still be fast enough for small event numbers). Using the previously discussed Root format will automatically include this speed benefit.
The hypbinning.txt file lists the low and high corner of each hypervolume in the binning scheme with its associated bin content. The bin content gives the phase space bin number ∈ {−N , ..., −1, +1, ..., +N }. The coordinates are given in the order {m + , m − , cos θ + , cos θ − , φ}; where invariant masses are given in units of MeV, and φ is given in radians.
To describe the format of the hypbinningwlinks.txt file, it is useful to revisit how the binning algorithm works. At iteration 0, there is one hypervolume; at iteration 1,

results.txt
The central values, statistical uncertainties, and systematic uncertainties for the measured hadronic parameters.

statcor.txt
The statistical correlations between the measured hadronic parameters. systcor.txt The systematic correlations between the measured hadronic parameters. stat.root The central values, statistical uncertainties, and statistical correlations of the measured hadronic parameters in Root format. This can be loaded with the Root macro loadresults.C. syst.root The central values, systematic uncertainties, and systematic correlations of the measured hadronic parameters in Root format. This can be loaded with the Root macro loadresults.C. statsyst.root The central values, combined statistical and systematic uncertainties, and combined statistical and systematic correlations of the measured hadronic parameters in Root format. This can be loaded with the Root macro loadresults.C. hypbinning.root The hyper-binning scheme in Root format. Further description of how to use this file is described in Appendix C.2.

hypbinning.zip
A compressed directory containing the the hyper-binning scheme in a text file. Further description of how to use this file is described in Appendix C.2.

benchmark.txt
The four-vectors associated to 100 phase space points, and their associated bin numbers. This can be used to check that the phase space binning has been correctly implemented.

modpred.txt
The central values and uncertainties of the hadronic parameter model predictions.

modpredcor.txt
The correlations between the uncertainties of the hadronic parameter model predictions.

modpred.root
The central values, uncertainties, and correlations of the hadronic parameter model predictions in Root format. This can be loaded with the Root macro loadresults.C. modcompat.txt The compatibility between the measured hadronic parameters and the model predictions. this gets split to give two hypervolumes; at iteration 2 each of these gets split to give 4 hypervolumes etc. Rather than discard the hypervolumes from iteration 0 and iteration 1, these can be kept to speed up the binning process later. The final set of hypervolumes that come out of the binning algorithm are known as 'bins' (B). Other hypervolumes that were used during the binning algorithm (but were then further divided) are known simply as 'volumes' (V). The first volume from iteration 0 is known as the primary volume (PV). A simple example of a 2 dimensional binning scheme, iteration-by-iteration, is given in Fig. 14, with bins and volumes labelled. Each volume and bin has a unique identifier called a 'volume number'. Every volume has links to two volume numbers, whereas each bin has a bin content (which gives the phase space bin number). The simple binning scheme in Fig. 14 is described by the information in Fig. 15, which has the same format as hypbinningwlinks.txt. For comparison, the same binning scheme is described in the same format as hypbinning.txt in Fig. 16. The general use case of a binning scheme is to find the bin (and its associated bin content), that an arbitrary phase space point, p, falls into. Using the information in hypbinning.txt requires looping over every bin, and seeing which one contains p; on average this will take ∼ N/2 operations, where N is the number of bins. To use the information in hypbinningwlinks.txt, one would first check if p is within the PV; if it is, then one would see which of the linked volumes p falls into etc. On average this will take ∼ log 2 N operations.  Figure 14. Simple example that demonstrates how the hyper-binning algorithm works. At iteration 0 there is a single primary volume (PV) with volume number 0. At iteration 1, the primary volume is split into two volumes with volume numbers 1 and 2. Volume number 1 is not split any further, so it is labelled as a 'bin' (B) rather than a 'volume' (V) -the content of this bin is -3. In iteration 2, volume number 2 is further divided into volume numbers 3 and 4; since this is the final iteration, these volumes are labelled as bins, which have bin contents of 5 and -1 respectively.