Unbinned model-independent measurements with coherent admixtures of multibody neutral D meson decays

Various studies of Standard Model parameters involve measuring the properties of a coherent admixture of \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${D} ^0$$\end{document}D0 and \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${\overline{D}{}} {}^0$$\end{document}D¯0 states. A typical example is the determination of the Unitarity Triangle angle \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\gamma $$\end{document}γ in the decays \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$B\rightarrow DK$$\end{document}B→DK, \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$D\rightarrow {{K} ^0_\mathrm{\scriptscriptstyle S}} {{\pi } ^+} {{\pi } ^-} $$\end{document}D→KS0π+π-. A model-independent approach to perform this measurement is proposed that has superior statistical sensitivity than the well-established method involving binning of the \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$D\rightarrow {{K} ^0_\mathrm{\scriptscriptstyle S}} {{\pi } ^+} {{\pi } ^-} $$\end{document}D→KS0π+π- decay phase space. The technique employs Fourier analysis of the complex phase difference between \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${D} ^0$$\end{document}D0 and \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${\overline{D}{}} {}^0$$\end{document}D¯0 decay amplitudes and can easily be generalised to other similar measurements, such as studies of charm mixing or determination of the angle \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\beta $$\end{document}β from \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${{B} ^0} \rightarrow D h^0$$\end{document}B0→Dh0 decays.


Introduction
Precise measurements of CP violation in decays of beauty hadrons is one of the key methods to search for effects of physics beyond the Standard Model. The phenomenon of CP violation is described in the Standard Model (SM) by the Cabibbo-Kobayashi-Maskawa (CKM) mechanism [1,2], where CP violation enters as a complex phase in the unitary 3 × 3 matrix (CKM matrix) describing transitions between quarks of the three generations due to charged-current weak interactions. A common representation of the CKM matrix is the Unitarity Triangle (UT), the sides and angles of which are experimentally observable parameters. The fundamental CP-violating phase, the angle γ of the UT (also known in the literature as φ 3 ), can be obtained with extremely low theoretical uncertainty [3] from tree-dominated b hadron decays and thus serves as a "standard candle" for searches of effects beyond the Standard Model in other heavy flavour processes.
Various techniques have been proposed to measure γ experimentally in the decays of B mesons into final states with neutral D mesons [4][5][6][7]. The CP violation in these a e-mail: anton. poluektov@cern.ch decays is generated by interference of b → c and b → u quark level transitions once the neutral D meson is reconstructed in a final state accessible to both D 0 and D 0 decays. The neutral D meson in this case forms a coherent admixture of D 0 and D 0 states which is denoted here as D. One of the most sensitive techniques involves analysis of the Dalitz plot density of multibody D decays such as D → K 0 S π + π − [8,9]. Two different techniques have been developed and implemented experimentally to extract γ from B → DK decays using multibody D meson final states. One is modeldependent, with the complex amplitude of the D decay obtained by fitting the flavour-specific D 0 decay density to a model [10][11][12][13][14][15][16]. This technique offers optimal statistical precision since the fit can be performed in an unbinned fashion, however, it suffers from uncertainty, which is difficult to quantify, due to modelling of the D 0 amplitude. Another method is a binned model-independent approach, where information as regards the behaviour of the strong phase across the phase space of the D 0 decay is obtained from samples of quantum-correlated D 0 D 0 decays produced near kinematic threshold [8,[17][18][19][20].
In the model-independent technique, one needs to determine the relation between the decay densities of quantumcorrelated D 0 D 0 decays and D decays from B → DK . This necessarily requires estimation of the decay density from scattered data, which is achieved by binning both decay densities. Each bin is assigned a number of parameters that characterise the averaged behaviour of the amplitude (its magnitude and phase) over the bin; these parameters are obtained by solving a system of equations that also includes the value of γ . In general, the binned approach reduces statistical sensitivity compared to the unbinned model-dependent technique, but the procedure is developed in such a way that it produces an unbiased measurement even in the case of a very rough binning.
In this paper, a method to extract γ is proposed which does not involve binning and aims to combine the advantages of the model-dependent and model-independent approaches. Like the binned approach with optimal binning, it uses a construction inspired by a D 0 amplitude model, but provides an unbiased measurement even if the wrong model is used. It is shown to offer better statistical sensitivity than the binned approach. The method employs Fourier analysis of the distribution of the complex phase difference between the D 0 and D 0 amplitudes. The method is illustrated using the "golden" channel B → DK with subsequent D → K 0 S π + π − decay, but can easily be generalised to other cases of γ determination where the binned model-independent technique is applicable: analyses using other three-or four-body D 0 decays [21][22][23][24][25][26], multibody B decays [27,28] or analyses using correlated Dalitz plots of multibody B-and D-meson decays [29,30].
Apart from measurements of γ , similar model-independent techniques, which employ interference between D 0 and D 0 amplitudes, have been developed for other kinds of measurements: studies of CP violation and mixing parameters of D 0 mesons [31][32][33], measurements of the UT angle β in B 0 → Dh 0 (where h 0 is a neutral light meson) and B 0 → Dπ + π − decays [34,35]. In all these cases, the technique proposed can be applied instead of the binned methods.

Model-independent formalism with weight functions
In this section, the formalism for γ measurement is recalled to introduce the notation, and the established model-independent technique is reformulated in slightly different terms. This allows a demonstration that the binned approach is not the only possible method to perform such a measurement.
Measurements of γ based on B → DK processes use the fact that the decay involves the interference of tree-dominated b → c and b → u diagrams, which produce neutral D mesons with opposite flavours. In the case of B + → DK + decays followed by D → K 0 S π + π − , the amplitude as a function of two variables of the D decay Dalitz plot, the squared invariant masses m 2 where the first term is due tob →c and the second due tob →ū transition. Here is that, for the D 0 → K 0 S π + π − decay, r B is the relative magnitude of the two contributions and δ B is the CP-conserving strong phase between them. The amplitude A B for the CP-conjugated decay B − → DK − can be obtained by replacing γ → −γ and swapping the D decay amplitudes: A D ↔ A D . A simultaneous analysis of the two amplitudes A B and A B provides information on the unknown parameters γ , r B , and δ B .
Experimentally, one deals with probability densities rather than amplitudes. The decay density for B + → DK + decays as a function of z ≡ (m 2 where the Cartesian CP-violating observables are introduced: x + = r B cos(δ B + γ ) and y + = r B sin(δ B + γ ).
The decay density p B (z) for B − → DK − decay involves the corresponding parameters x − = r B cos(δ B − γ ) and The expressions for the decay densities can be rewritten as where h B andh B are the normalisation factors, and p D (z) andp D (z) are the Dalitz plot densities of flavour-tagged i.e. the Dalitz plot distributions for D 0 and D 0 decays are symmetric under the exchange m 2 + ↔ m 2 − assuming CP conservation in D 0 decays. 1 The functions C(z) and S(z) contain information as regards the motion of the complex strong phase over the Dalitz plot which cannot be obtained from flavour-specific D meson decays: One needs to know them to obtain the values of CP violating parameters x ± and y ± fromp B (z) and p B (z).
In the model-dependent approach to measure γ , the strong phase motion is fixed by an amplitude model. The modelindependent technique employs pairs of neutral D mesons produced at the kinematic threshold in the e + e − → D 0 D 0 process to obtain this information. In this case, the two D mesons are produced in a P-wave such that their wave function is antisymmetric. As a result, if both D mesons are reconstructed in the K 0 S π + π − final state, the densities of two Dalitz plots will be correlated: Here the indices "1" and "2" correspond to the two decaying D mesons and h D D is a normalisation factor. The necessary information as regards C(z) and S(z) is present in Eq. (8), but it is not straightforward to obtain the explicit expressions for the functions C(z) and S(z) from the observable distributions p D (z),p D (z) and p D D (z 1 , z 2 ). Equation (8) contains an ambiguity: it is invariant under rotation of the pair C(z), S(z) by an arbitrary phase : This does not constitute a significant problem since it effectively results in the redefinition of the strong phase δ B , leaving the CP-violating phase γ unaffected. The other abiguity is the change of sign of C(z) or S(z), which results in the change of sign for γ . Other decays of D mesons from correlated D 0 D 0 pairs can offer additional information to resolve these ambiguities. For instance, decays where one of the D mesons is reconstructed in a CP eigenstate and the other is reconstructed as K 0 S π + π − constrain C(z), and they resolve the ambiguity (9), as well as fix the sign for C(z). The remaining ambiguity, the sign of S(z), can be resolved by a weak model assumption using isobar parametrisation of the D decay amplitude [18]. In practice, several D decay modes are combined to measure the same strong phase parameters [36], but the description below will concentrate only on D 0 D 0 pairs where both D mesons are decaying to K 0 S π + π − . The model-independent technique can be built based on the observation that explicit expressions for the functions C(z) and S(z) are not needed to obtain x ± , y ± . One can derive a number of independent equations from the expressions (8) and (4) (10) while Eq. (4) becomē where and The integration in Eqs. (11)-(13) is performed over the entire Dalitz plot D of the D decay, while for Eq. (10) double integral is performed over the Dalitz plots D 1 and D 2 of two decaying D mesons. Unlike in the binned formalism described in Refs. [17,18], here the terms proportional to |A D (z)| · |A D (z)| are not factored out, thus capital letters are used to distinguish the expressions of Eq. (13) from c i and s i coefficients commonly used in the binned formalism. The values of weighted integrals for the flavour-specific D sample ( p n andp n ), B sample (p B,n and p B,n ) and correlated D 0 D 0 sample ( p D D,mn ) can be obtained directly from each of the corresponding scattered data samples by replacing the integrals with sums over individual observed events. The values of the weighted integrals for the phase terms C n and S n are considered as free parameters constrained by Eq. (10). This allows the values of x ± and y ± to be obtained by solving the system of Eq. (10) and (11).
The family of weight functions w n can be chosen arbitrarily, but the performance of the method with a limited data sample will depend on this choice. The binned modelindependent approach is a particular case of the considered formalism where the weight functions are of the form Here D n are non-overlapping regions of the Dalitz plot which define the bins.
To reach optimal statistical sensitivity, the binning has to be chosen in such a way as to maximise the interference term in Eq. (11). A good approximation to the optimum is known to be the binning based on the strong phase difference between the favoured and suppressed D decay amplitudes [18]. Specifically, if one defines the phase difference then the bin D n (1 ≤ n ≤ M) is the region of the phase space which satisfies The bins in the region with m 2 + > m 2 − are defined symmetrically with respect to exchange m 2 is an amplitude model that ideally should approach the true amplitude A D (z) to reach optimal statistical precision, but does not need to match it exactly to provide an unbiased measurement.
The following section shows how to construct an unbinned model-independent formalism using a model-based phasedifference function (m 2 + , m 2 − ) which will be a generalisation of the technique with phase-difference binning. For reasons which will become obvious, this approach will not be optimal from the point of view of statistical uncertainty, and it will serve solely as a demonstration. Subsequently, a more optimal approach based on a similar construction will be presented.

Unbinned technique using Fourier series expansion of phase difference
Let (z) ≡ (m 2 + , m 2 − ) be the function defined by Eq. (15) that maps two-dimensional Dalitz plot coordinates z to the one-dimensional space represented by a phase difference φ between the D 0 and D 0 amplitudes at the same Dalitz plot point. One can now define probability densities as functions of φ = (z). The density of the flavour-specific D decay becomes From the experimentalist's point of view, this function is the probability density (PDF) of the (z) value for a sample of flavour-specific D → K 0 S π + π − decays, and is a continuous generalisation of the number of events K n that enter the n th bin in the approach with binning based on equal phase difference [18]. Following Eqs. (6) and (15), the density for the CP-conjugate decay is After a similar mapping is applied to the correlated densities of the two D → K 0 S π + π − Dalitz plots of the D 0 D 0 sample (8), the following PDF of the variables φ 1 = (z 1 ) and φ 2 = (z 2 ) is obtained: where From the definitions (20) and (15) it follows that C(φ) is an even function, while S(φ) is odd: Switching to the phase-difference representation for the B ± → DK ± densities (4), one obtains The next step is to choose the family of weight functions to construct a system of equations which allow the determination of x ± and y ± from Eqs. (19) and (22). Since the densities as functions of φ are periodic by construction, it appears that the natural choice is to use Fourier expansion of the functions of the phase difference, i.e. use weight functions of the form cos(nφ) and sin(nφ), where n is an integer number. The unknowns x ± and y ± will then enter the system of equations which relates the coefficients of the Fourier expansions of the p D , p D D ,p B and p B densities.
Specifically, the functions p D (φ), C(φ) and S(φ) can be represented as [a D D mn cos(mφ 1 ) cos(nφ 2 ) Strictly speaking, the equations above are exact only in the limit M → ∞; however, in practice one has to truncate the Fourier series at a certain finite M.
For p D (φ), the values of the Fourier coefficients can be calculated directly from scattered data φ (i) where N D is the number of events in the data sample and φ (i) = (z (i) ) are the calculated phase-difference values for the data sample entries z (i) . Similarly, the coefficients of the Fourier expansion for the correlated D 0 D 0 sample can be calculated from the 2D scattered data φ On the other hand, from Eq. (19) one can obtain a set of relations between the Fourier coefficients for flavour-specific and D 0 D 0 densities: Equations (30) can be used to obtain the unknown coefficients a C n and b S n from the known values of (a, b) D n and (a, b, c, d) D D mn . The system of Eq. (30) is solvable for any M ≥ 1 (there are 2M 2 + M + 1 independent equations and 2M + 2 unknown parameters). In practice, since the system of equations is overconstrained for M > 1, it should be solved using a maximum likelihood fit, which will also provide estimate of the covariance matrix.
A maximum likelihood fit needs uncertainties for the coefficients that enter the equations. These can be calculated analytically by applying a Poisson bootstrapping technique [37]. Each term entering the sum in Eq. (28) or (29) is multiplied by a random number which follows the Poisson distribution with unit mean value. The variances for the sums can then be obtained assuming they have a Gaussian distribution (which is a valid assumption for large N D ): and In addition, unlike in the binned case where the yields in each of the bins are statistically independent, the coefficients of the Fourier series are in general correlated. The covariance matrix can be calculated similarly using Poisson bootstrapping, e.g. the covariance between the a n and b m coefficients can be calculated as Similarly, the expressions for covariances between a n and a m , b n and b m , or between the coefficients (a, b, c, d)   and plugging them into Eq. (22), one obtains the following system of equations: which can be solved again using a maximum likelihood fit for any M ≥ 1, after the extraction of the coefficients (ā,b) B n and (a, b) B n and their uncertainties and correlations from the B → DK sample in a similar way. Alternatively, both sets of equations (30) and (35) can be solved simultaneously using a single combined likelihood.
As an illustration, the functions p D (φ), p D D (φ 1 , φ 2 ), C(φ) and S(φ) obtained using the D → K 0 S π + π − amplitude model A D (m 2 + , m 2 − ) from the Belle measurement [14] and their respective coefficients of the Fourier expansion are shown in Figs. 1, 2 and 3. The function p D (φ) shown in Fig. 1a is obtained by plotting the distribution of the function φ = (z) (black points) for events generated according to PDF p(z). Its Fourier coefficients a D n and b D n up to n = 19 are shown in Fig 1b, c, respectively. Since the normalisation is arbitrary, the coefficients are normalised such that a D 0 = 1.
The solid red line in Fig. 1a shows the result of Fourier expansion up to n = 19, and the dashed blue line shows the first harmonic (expansion up to n = 1). In the p D D (φ) function that is obtained similarly by plotting the two-dimensional distribution of φ 1 = (z 1 ), φ 2 = (z 2 ) for the correlated Dalitz plot points generated according to the p D D (z 1 , z 2 ) density ( Fig. 2), only the a D D mn and d D D mn coefficients are nonzero, while b D D mn and c D D mn are consistent with zero as expected from Eq. (30). The normalisation a D D 00 = 1 is used. The true functions C(φ) and S(φ) can be obtained, on one hand, from the known amplitude A D (z) by plotting the distribution of the function φ = (z) for events generated uniformly across the Dalitz plot with event-by-event weights √ p D (z)p D (z) cos (z) and √ p D (z)p D (z) sin (z), respectively. These functions are shown in Fig. 3a, b as black points. On the other hand, the functions can be reconstructed from the spectral coefficients a C n and b S n obtained from Eq. (30). The fitted coefficients a C n and b S n are plotted in Fig. 3c, d, while the functions C(φ) and S(φ) reconstructed from them are shown in Fig. 3a, b as solid red lines (from the coefficients up to n = 19) and dashed blue line (only one harmonic, n = 1). It can be seen from Fig. 3c, d that the highest "power" of the C(φ) and S(φ) spectrum is contained in the first harmonic, n = 1. As a result, as will be seen from further studies with pseudoexperiments, limiting M = 1 is sufficient to reach good sensitivity to γ .

(d) ( e ) 4 Strategy with Fourier expansion on split Dalitz plot
The strategy outlined above is the simplest example of the approach using Fourier expansion of the phase-difference distribution to measure γ . However, it is clear that this approach is not optimal from the point of view of statistical precision. The reason is that one integrates over all points of the phase space with the same expected phase difference, regardless of the magnitudes of the interfering D 0 → K 0 S π + π − and D 0 → K 0 S π + π − amplitudes. This effectively reduces the interference term and, as a consequence, the sensitivity to the relative phase between the two amplitudes. For similar reasons, the "optimal" binning scheme was introduced for the binned model-independent approach to improve the precision with the equal phasedifference binning [18].
The simplest way to improve the situation with the technique described here is to split the Dalitz plot into regions with comparable ratios between the absolute values of the interfering amplitudes, and to perform Fourier expansion in those regions separately. This approach is illustrated below in an example with the Dalitz plot split into two regions.
Two regions of the Dalitz plot are considered: one with p D (z) > p D (z) (denoted region D + ) and the other with p D (z) < p D (z) (region D − ). These are shown in Fig. 4 for the same D → K 0 S π + π − amplitude model as used in the previous section. Clearly, the exchange m 2 + ↔ m 2 − transforms D + into D − . Now one has to deal with two independent dis- tributions for the D → K 0 S π + π − density, p + D and p − D , as functions of the phase difference φ, defined as The corresponding distributions for the CP-conjugated decays arē It should be stressed that the superscripts "+" and "−" denote two Dalitz plot regions rather than B meson flavours. Throughout this paper, the flavour (b orb) is consistently denoted by the absence or presence of a "bar" in the corresponding quantities, for examplep B and p B , except for the subscript for CP-violating parameters x ± , y ± , which is a commonly used notation.
With the Dalitz plot split in this way, one needs to define two sets of functions C(φ) and S(φ) in the two Dalitz plot regions: These functions will not be even and odd, as in the previous example, but instead they will satisfy the following properties: The two-dimensional density of the D 0 D 0 sample will be described by a set of four functions p ++ D D , p +− D D , p −+ D D and p −− D D defined as The Fourier expansion coefficients (a, b) C± n and (a, b) S± n for the C ± (φ) and S ± (φ) functions, respectively, are defined as and they are related as The relations between the coefficients of the Fourier expansion of the D 0 D 0 and flavour D decay densities in that case take the following form: where the coefficients (a, b, c, d)  Finally, the equations for the densities of the D decay from B ± → DK ± take the following form in the split Dalitz plot case: The number of unknown D phase parameters in the equations has now increased: there are 4M +2 independent coefficients (a, b) C,S+ n (0 ≤ n ≤ M for a and 1 ≤ n ≤ M for b) plus a common normalisation factor h D D in the system of equations (46). Nevertheless, the statistical precision in this approach appears to be better as will be seen in the feasibility study.
In principle, one can even consider splitting the Dalitz plot into more regions, but certainly the increase in the number of free parameters can diminish the possible gain in statistical precision. Any strategy involving splitting the Dalitz plot should be optimised taking into account the size of experimentally available samples of correlated D 0 D 0 and B → DK decays.

Simulation results
To test the feasibility of the proposed method, simulation studies using pseudoexperiments are performed. Samples of flavour-specific D → K 0 S π + π − decays, correlated D 0 D 0 pairs decaying to K 0 S π + π − , and D → K 0 S π + π − decays from B → DK are generated using the D decay amplitude measured by Belle collaboration [14]. Samples are simulated with r B = 0.1, γ = 60 • and δ B = 130 • , which is close to the results of the recent model-independent measurement of the B → DK , D → K 0 S π + π − channel by the LHCb collaboration [19]. For each of those event samples, the Fourier series coefficients and their covariance matrices are calculated as described in Sect. 3. Systems of equations which contain relations between Fourier spectrum coefficients of flavour-specific D, D 0 D 0 and B → DK densities are then solved by maximising the combined likelihood to obtain the value of γ .
The formalism in Sects. 3 and 4 involved Cartesian CPviolating parameters x ± and y ± . This approach is likely more suitable when dealing with real data when one has to combine the results of different γ -sensitive analyses. In the simulation study presented here, the free parameters are chosen to be γ , r B and δ B .
For the flavour-specific D → K 0 S π + π − mode, a large sample of 10 7 generated events is used. This sample is not expected to contribute significantly to the uncertainty on γ since high-statistics data sets are available at both the B factories and the LHCb. The size of the B → DK , D → K 0 S π + π − sample generated is 10 4 events for each B meson flavour, which corresponds roughly to 10 times the data sample from LHCb Run 1 [19]. Three scenarios with different correlated D 0 D 0 sample sizes are considered, 10 5 , 10 4 and 10 3 events. For comparison, the e + e − → D 0 D 0 data sample collected by CLEO experiment where both D mesons decay into K 0 S π + π − contains 473 events, however, many other D decay modes are used in the combined fit to obtain the phase coefficients (notably, the modes where one of the D mesons is reconstructed in a CP eigenstate or as K 0 L π + π − ) [36]. It is expected that the statistical uncertainty of the D 0 D 0 sample of 10 5 events will contribute negligibly to the uncertainty on γ , thus pseudoexperiments with this sample size probe uniquely how the approximation of the amplitude with a finite number of parameters (i.e. truncated Fourier series or limited number of bins) affects γ sensitivity. The low-statistics sample of 10 3 events, on the other hand, will demonstrate the contribution of a limited D 0 D 0 sample to the sensitivity.
Each ensemble of pseudoexperiments is fitted with the binned model-independent procedure with 3, 5, 8, 12, and 20 bins using both the phase-difference and the "optimal" binning schemes [18], and with the two Fourier analysis techniques outlined above, using the entire Dalitz plot or the Dalitz plot split in two regions, respectively. In the approaches with Fourier expansion, the limit M on the number of harmonics is set to M = 1, 2, 4, 7, 11, or 19. In addition, an unbinned model-dependent fit is performed to serve as a reference for the best possible statistical γ precision that can be reached.
The Fourier expansion approach is verified to produce unbiased results if different D → K 0 S π + π − amplitudes are used for event generation and calculation of the phase difference (z). This is certainly a requirement for a technique to be model-independent. This check is performed by using a reduced D → K 0 S π + π − model where a subset of two-body amplitudes is present (ρ(770) 0 , ω(782) and f 0 (980) in the π + π − amplitude, and K * (892) ± in the K 0 S π amplitudes) plus a flat non-resonant term. The results in Fig. 5 are shown for M = 1 and 10 5 D 0 D 0 events, but a similar check is performed for each value of M. Figure 6 shows the γ and r B resolutions as functions of the number of bins (for the binned scenarios) and the number of Fourier expansion terms (for the unbinned scenarios) with the four fit strategies described above and for the three different D 0 D 0 sample sizes. For comparison, the uncertainty of the unbinned model-dependent fit is also shown. While the precision of the binned approaches depends on the number of bins, the uncertainty of the Fourier expansion techniques practically does not depend on the number of harmonics M for relatively large D 0 D 0 samples sizes, while for a small D 0 D 0 sample size of 10 3 the optimum is reached for M = 1 (i.e. for the smallest possible number of free parameters, which is three for non-split and six for split Dalitz plot). It is possible that other multibody D decays may require higher harmonics to reach optimal sensitivity. Another case when Fourier terms with n > 1 might be required is if the amplitude model A The γ uncertainties for the optimal scenarios with the binned and unbinned techniques are compared in Table 1.  The numbers correspond to the best γ resolution obtained in a range of M (see Fig. 6). For comparison, the γ uncertainty for unbinned model-dependent fit is σ (γ ) = 2.91 ± 0.07 • The uncertainty of the approach with split Dalitz plot is significantly better than when the Dalitz plot is taken as a whole.
It is also clear that the Fourier expansion technique with split Dalitz plot shows better sensitivity than the binned method using "optimal" binning, with the gain being the most significant for smaller D 0 D 0 sample size. The technique, however, is still about 10% less sensitive than the unbinned modeldependent approach. The possibilities to further improve the sensitivity of the unbinned model-independent method are discussed in Sect. 7.

Practical considerations
To be applicable to real data, the technique should be able to deal with experimental effects such as backgrounds and non-uniform detection efficiency across the Dalitz plot. Since the background enters the decay density additively, it can be treated at the level of Fourier-transformed variables, by calculating the Fourier expansion of the background density and subtracting it from the coefficients calculated on data.
On the other hand, the efficiency enters the density in a multiplicative way, thus Fourier expansion need to be applied to efficiency-corrected data. The correction can be applied on an event-by-event basis, by assigning each event a weight proportional to the inverse of efficiency while calculating the Fourier coefficients. The studies presented above have been performed using a combined likelihood fit to both B → DK and correlated D 0 D 0 samples. It is also possible to perform the analysis in two stages, by first calculating the coefficients of Fourier transformation of the functions C(φ) and S(φ) from the D 0 D 0 data, followed by a fit to B → DK sample using the coefficients, their correlations and uncertainties from the first stage. This is likely to be more convenient in practice, since the data samples come from different experiments.

Further directions of development
Using notation of the generalised model-independent formalism presented in Sect. 2, the Fourier analysis technique proposed above uses a family of 2M + 1 weight functions where 1 ≤ n ≤ M. The use of the function (z) ensures that different points in the phase space do not cancel each other out while calculating the integral, and thus the interference term that provides sensitivity to CP-violating observables is large (assuming, of course, that the amplitude model that provides (z) is close to the true amplitude). However, information about the absolute value of the amplitude is ignored in the formalism presented in Sect. 3 and is taken into account only rather roughly in Sect. 4. Alternatively, one could consider a weight function that in addition takes into account the magnitudes of the favoured and suppressed amplitudes from the model |A (z)|, and thus adds more information to maximise the interference term. In the presence of background, the family of weight functions should also take into account the distribution of background events over the phase space. Further optimisation of the family of weight functions needs additional study.
The proposed technique could be especially useful in the cases where a binned approach will limit precision due to small sample sizes of decays which determine the phase information. Examples are the D 0 → K 0 S K + K − mode, where the sample of quantum-correlated decays is small and currently only two bins are used in the γ measurement [19]. Another example is B → DK π decays, where the phase coefficients corresponding to the three-body B decay are free parameters together with γ [29,30]. Having an amplitude model which describes the strong phase variation across the B decay Dalitz plot with a small number of parameters should improve the statistical sensitivity.
Other analyses, where the coherent D 0 -D 0 admixtures are involved, are measurements of charm mixing and CP violation in mixing and measurement of the UT angle β in B → Dh 0 decays. These classes of measurements utilise oscillations of D 0 and B 0 mesons, respectively, and thus the parameters of the D 0 -D 0 admixture are functions of decay time. In the proposed formalism, the coefficients of the Fourier series will be functions of decay time as well. While such analyses will certainly be more complicated than the case with constant coefficients, they are conceptually similar to the measurements using the binned technique which have already been carried out [33,35].

Conclusion
A technique to perform unbinned model-independent analysis of a coherent admixture of D 0 and D 0 states decaying to a multibody final state is proposed. It is illustrated in detail using the measurement of the UT angle γ from B → DK decays. Unlike the well-known technique with Dalitz plot binning, the proposed method employs Fourier analysis of the spectrum of the strong phase difference between the D 0 and D 0 amplitudes. While the method relies on an amplitude model to reach optimal statistical precision, it is unbiased by construction even if the wrong model is used.
A study of the feasibility of the proposed method has been performed with simulated pseudoexperiments. The precision of the method does not depend strongly on the number of Fourier expansion terms used, and even with only the single leading term yields sensitivity comparable to that of the binned model-independent approach. A modification of the procedure, where Fourier expansion is performed in two regions of the Dalitz plot separated according to the ratio of the suppressed and favoured amplitudes, provides γ sensitivity better than the most optimal binned strategy. The gain compared to the binned approach is especially significant if the size of the correlated D 0 D 0 sample, which determines the strong phase in D meson decay, is small. Possible ways of improving the sensitivity of the proposed technique even further are identified and need further study.
The method is not limited to γ measurements with threebody D decays and can be generalised to any analysis where the parameters of a coherent admixture of D 0 and D 0 in a multibody final state need to be determined, such as measurements of charm mixing and CP violation, and measurements of the UT angle β in B → Dh 0 decays. The technique could also be useful in γ measurements with a double Dalitz plot analysis of B → DK π , D → K 0 S π + π − decay; in that case the Fourier expansion can be applied to both the B and the D Dalitz plots.