Empirical parameterization of the K+- ->pi+- pi0 pi0 decay Dalitz plot

As first observed by the NA48/2 experiment at the CERN SPS, the $\p0p0$ invariant mass (M00) distribution from $\kcnn$ decay shows a cusp-like anomaly at M00=2m+, where m+ is the charged pion mass. An analysis to extract the pi pi scattering lengths in the isospin I=0 and I=2 states, a0 and a2, respectively, has been recently reported. In the present work the Dalitz plot of this decay is fitted to a new empirical parameterization suitable for practical purposes, such as Monte Carlo simulations of K+- ->pi+- pi0 pi0 decays.

However, in 2005 the NA48/2 experiment at the CERN SPS first observed a cusp-like anomaly in the π 0 π 0 invariant mass (M 00 ) distribution of this decay in the region around M 00 = 2m + , where m + is the charged pion mass [5]. This anomaly had been predicted in 1961 [6] as an effect due mainly to the destructive interference between the direct amplitude of K ± → π ± π 0 π 0 decay and the final state charge exchange scattering process π + π − → π 0 π 0 in K ± → π ± π + π − decay (see also [7]).
Best fits using two theoretical formulations of rescattering effects [8] and [9,10] have provided a precise determination of a 0 −a 2 , the difference between the S-wave ππ scattering lengths in the isospin I = 0 and I = 2 states, and an independent, though less precise, determination of a 2 [11]. Such an analysis leads to a successful fit of the Dalitz plot using a rather long expression which contains physically meaningful constants (that could be measured better in future) and is affected by theoretical uncertainties. It is not practical to implement these formulae if one just needs to describe the Dalitz plot shape, say, in a Monte Carlo generator. On the other hand, it is known now [5] that the Dalitz plot region near s 3 = (2m + ) 2 cannot be described by Eq. (1), so a model-independent, empirical description of K ± → π ± π 0 π 0 decay is certainly useful to replace the old parameterization.
The main purpose of the NA48/2 experiment at the CERN SPS was to search for direct CP violation in K ± decay to three pions [12,13,14]. The experiment used simultaneous K + and K − beams with a momentum of 60 GeV/c propagating along the same beam line. Data were collected in 2003-04, providing large samples of fully reconstructed K ± → π ± π + π − and K ± → π ± π 0 π 0 decays. Here we report the results from a study of a partial sample of ∼ 30.4 × 10 6 K ± → π ± π 0 π 0 decays recorded in the second half of the 2004 run with the purpose of providing a new empirical, model-independent parameterization of the K ± → π ± π 0 π 0 Dalitz plot. This parameterization describes the K ± → π ± π 0 π 0 experimental data with no distortions from instrumental effects, such as resolution, geometrical acceptance and detection efficiency, as they would be measured by a detector with full acceptance and ideal performance. It could also be useful, therefore, in the development of new theoretical formulations of rescattering effects in K ± → π ± π 0 π 0 decay, or in the refinement of existing ones.

Beam and detectors
The two simultaneous beams are produced by 400 GeV protons impinging on a 40 cm long Be target. Particles of opposite charge with a central momentum of 60 GeV/c and a momentum band of ±3.8% produced at zero angle are selected by a system of dipole magnets forming an "achromat" with null total deflection, focusing quadrupoles, muon sweepers and collimators. With 7 × 10 11 protons per burst of ∼ 4.5 s duration incident on the target the positive (negative) beam flux at the entrance of the decay volume is 3.8 × 10 7 (2.6 × 10 7 ) particles per pulse, of which ∼ 5.7% (∼ 4.9%) are K + (K − ). The decay volume is a 114 m long vacuum tank with a diameter of 1.92 m for the first 66 m, and 2.4 m for the rest.
Charged particles from K ± decays are measured by a magnetic spectrometer consisting of four drift chambers (DCH) and a large-aperture dipole magnet located between the second and third chamber [16]. Each chamber has eight planes of sense wires, two horizontal, two vertical and two along each of two orthogonal 45 • directions. The spectrometer is located in a tank filled with helium at atmospheric pressure and separated from the decay volume by a thin (0.0031 radiation lengths, X 0 ) Kevlar window. A 16 cm diameter aluminium vacuum tube centered on the beam axis runs the length of the spectrometer through central holes in the Kevlar window, drift chambers and calorimeters. Charged particles are magnetically deflected in the horizontal plane by an angle corresponding to a transverse momentum kick of 120 MeV/c. The momentum resolution of the spectrometer is σ(p)/p = 1.02% ⊕ 0.044%p (p in GeV/c), as derived form the known properties of the spectrometer and checked with the measured invariant mass resolution of K ± → π ± π + π − decays. The magnetic spectrometer is followed by a scintillator hodoscope consisting of two planes segmented into horizontal and vertical strips and arranged in four quadrants.
A liquid Krypton calorimeter (LKr) [17] is used to reconstruct π 0 → γγ decays. It is an almost homogeneous ionization chamber with an active volume of ∼ 10 m 3 of liquid krypton, segmented transversally into 13248 2 cm × 2 cm projective cells by a system of Cu-Be ribbon electrodes, and with no longitudinal segmentation. The calorimeter is 27 X 0 thick and has an energy resolution σ(E)/E = 0.032/ √ E ⊕ 0.09/E ⊕ 0.0042 (E in GeV). The space resolution for a single electromagnetic shower can be parameterized as σ x = σ y = 0.42/ √ E ⊕ 0.06 cm for each transverse coordinate x, y. A neutral hodoscope consisting of a plane of scintillating fibers is installed in the LKr calorimeter at a depth of ∼ 9.5 X 0 . It is divided into four quadrants, each consisting of eight bundles of vertical fibers optically connected to photomultiplier tubes.

Event selection
A specific subset (about 50%) of the full data sample (collected in 2003 and 2004) was used, recorded with optimised trigger conditions allowing precise control of the trigger efficiency.
K ± → π ± π 0 π 0 events were recorded by a first level trigger using signals from the scintillator hodoscope (Q1) and LKr (NUT), followed by a second level trigger using drift chamber information (MBX). Events were also recorded using other triggers with different downscaling factors for different periods: a minimum bias NUT trigger (ignoring both Q1 and MBX); and a minimum bias Q1*MBX trigger (ignoring LKr information). Using the event samples recorded with these downscaled triggers, and selecting K ± → π ± π 0 π 0 decays, it was possible to measure separately the efficiency of the minimum bias Q1*MBX trigger using the event sample recorded by the minimum bias NUT trigger and the efficiency of the minimum bias NUT trigger using the events recorded by the minimum bias Q1*MBX trigger. These two efficiencies were multiplied together to obtain the full trigger efficiency, which was always above 94% for the data sample used in this analysis. Details of the trigger efficiency for K ± → π ± π 0 π 0 decay events are given in [12,14].
Events with at least one charged particle track having a momentum above 5 GeV/c, measured with a maximum error of 6%, and at least four energy clusters in the LKr, each consistent with a photon and above an energy threshold of 3 GeV, were selected for further analysis. In addition, the distance between any two photons in the LKr was required to be larger than 10 cm, and the distance between each photon and the impact point of any track on the LKr front face had to exceed 15 cm. Fiducial cuts on the distance of each photon from the LKr edges and centre were also applied in order to ensure full containment of the electromagnetic showers.
Every combination of four clusters and one track was considered as a K ± → π ± π 0 π 0 decay candidate if clusters were in time within 5 ns, and if the track was in time with the cluster average time within 10 ns. The distribution of the difference between the time of each cluster and their average value has an approximately Gaussian shape with σ ≈ 0.73 ns, while the distribution of the difference between the track time and the cluster average time has σ ≈ 1.5 ns, so these cuts accept almost all the time-correlated combinations. At this stage of event selection there is a ∼ 1.5% background associated with accidental LKr clusters. However, after the π 0 π 0 pair selection (see below) the level of residual accidental background, estimated from the distribution of the difference between the track time and the average time of the four clusters, is less than 0.02% and can be safely neglected.
Other rate effects, such as losses caused by mismeasurement of cluster and track parameters due to accidental activity in the detectors, were considered as part of the detector performance. The simulation of relevant resolutions and tails has been tuned to the experimental data, hence our Monte Carlo model includes also these rate effects. Residual discrepancies between experimental and simulated samples were taken into account in the study of systematic uncertainties (see section 7).
Each possible combination of two photon pairs in the event was assumed to originate from the two-photon decays of a pair of neutral pions, and for every π 0 candidate the position of the decay vertex along the beamline was calculated as , where Z LKr is the LKr longitudinal position, and E 1 , E 2 , x 1 , x 2 , y 1 , y 2 are the measured energies and transverse coordinates of the two photons, as measured in the LKr. The K ± → π ± π 0 π 0 decay vertex position Z was taken as the arithmetic average of the two Z π 0 values. The reconstructed decay vertex position Z was further required to be at least 2 m after the downstream end of the final beam collimator. In addition, the reconstructed kaon momentum was required to be between 54 and 66 GeV/c.
For each DCH plane the event energy-weighted center-of-gravity (COG) coordinates were calculated using the photon coordinates and energies, as measured by the LKr, and the track parameters before deflection, so COG represents the intersection of the initial kaon flight line with the DCH plane. Inner acceptance cuts were applied at each DCH plane to reject events with COG radius larger than R COG max (typically between 2 and 3 cm) 1 , and with the charged track closer than R COG−track min (typically between 15.5 and 19 cm) to the event COG. The exact cut values for every DCH plane have been chosen depending on the COG and track impact point distributions on that plane.
In order to reject events with photons emitted at very small angles to the beam and traversing the beam pipe in the spectrometer or the DCH1 central flange, and converting to e + e − before reaching the LKr, for each photon detected in LKr its distance from the nominal beam axis at the DCH1 plane was required to be > 11 cm, assuming an origin on axis at Z + 400 cm.
For every K ± → π ± π 0 π 0 decay candidate in the event, both the reconstructed π ± π 0 π 0 invariant mass (M) and the difference between the two Z π 0 coordinates (δZ) were used. For each K ± → π ± π 0 π 0 decay candidate an estimator χ 2 was defined as where the resolutions RMS z and RMS m have been parameterized from the experimental data as a functions of Z. The combination with the minimum χ 2 was chosen as the reconstructed K ± → π ± π 0 π 0 decay after applying the final loose cut χ 2 < 30.
The π ± π 0 π 0 invariant mass distribution is shown in Fig. 1. Non-gaussian tails, mainly associated with π → µν decays in K ± → π ± π 0 π 0 events, are suppressed by the χ 2 cut. There are also small contributions from wrong photon pairings in the decay of the two π 0 , and from non-gaussian tails of the LKr response due to photonuclear reactions. All these effects are included in the Monte Carlo simulation and are taken into account in the evaluation of the systematic uncertainties (Section 7). Radiative photons from K ± → π ± π 0 π 0 decays produce a slight shift of the measured 1 The beams were focused at DCH1, where the RMS values of their radial distributions were ∼ 0.45cm kaon mass, and thus also contribute to the tails of the χ 2 distribution. Our simulation does not take into account radiative photons, and we assume that the emission of soft real γ leaves the decay kinematics essentially unchanged. There is no limit to the presence of additional clusters in our event selection from the data. We have checked that the replacement of the χ 2 cut with the cut δZ < 500 cm (with no cuts on the measured π ± π 0 π 0 invariant mass) leads to a negligible change of the s 3 spectrum and of the fit results. So, within the present statistical uncertainty our analysis includes all the radiative K ± → π ± γπ 0 π 0 decays.
There are no important physical background sources for the K ± → π ± π 0 π 0 decay mode. Accidental overlaps of two events could produce some background, which, however, is expected to have a flat distribution in the δZ, M plane, hence a flat χ 2 distribution. If one interprets the small differences observed in the tails of the χ 2 distributions of data and MC events as totally due to this background rather than to the quality of the simulation, the accidental background can be conservatively estimated to be < 0.2%.
A total of 30.4 × 10 6 K ± → π ± π 0 π 0 decay candidates have been selected for the present analysis. Fig. 2 a) shows the distribution of the square of the π 0 π 0 invariant mass, M 2 00 , for the final event sample. This distribution is displayed with a bin width of 0.00015 (GeV/c 2 ) 2 , with the 51 st bin centered at M 2 00 = (2m + ) 2 , where m + is the charged pion mass (the M 2 00 resolution is 0.00031 (GeV/c 2 ) 2 at M 2 00 = (2m + ) 2 ). For our fits we use the bin interval 21 − 311 which contains the major part (> 98%) of selected events. The sudden change of slope near M 2 00 = (2m + ) 2 = 0.07792 (GeV/c 2 ) 2 , first observed in this experiment [5] is clearly visible.

Monte Carlo simulation
Samples of simulated K ± → π ± π 0 π 0 events ∼ 10 times larger than the data have been generated using a full detector simulation based on the GEANT-3 package [18]. This Monte Carlo (MC) program takes into account all known detector effects, including the time-dependent efficiencies and resolutions of the detector components.
The MC program also includes the simulation of the beam line. The beam average position and momentum are tuned for each period of few hours using fully reconstructed K ± → π ± π + π − events, which provide precise information on the average beam angles and positions. Furthermore, the requirement that the average reconstructed π ± π + π − invariant mass be equal to the nominal K ± mass for both K + and K − fixes the absolute momentum scale of the magnetic spectrometer.
The Monte Carlo simulation does not include the overlay of two independent K ± → π ± π 0 π 0 events or of a simulated K ± → π ± π 0 π 0 event with a randomly triggered one, so the timing cuts described in section 3 were not applied in the analysis of the simulated event sample. It should be noted that rate effects depend on the time structure of the SPS beam spills, which may vary from spill to spill during data taking and cannot be easily included in the Monte Carlo simulation.
The Dalitz plot distribution of K ± → π ± π 0 π 0 decays has been generated according to Eq.(1). For any given value of the generated π 0 π 0 invariant mass the simulation provides the detection probability and the distribution function for the reconstructed value of M 2 00 . This allows the transformation of any theoretical distribution into an expected distribution which can be compared directly with the measured one. , from K ± → π ± π 0 π 0 decay in the fit region. b) -relative deviation of the experimental spectrum from the best fit result (Data − F it)/F it.

Parameterization
In order to describe the cusp observed in the π 0 π 0 invariant mass distribution, we propose the following empirical parameterization for the square of the K ± → π ± π 0 π 0 decay matrix element: where and Here H is the Heaviside step function (H(x < 0) = 0, H(x ≥ 0) = 1) and δ is the Dirac delta function (in particular, +w/2 −w/2 δ(x)dx = 1). The constant U t = −1.1272 is the U value at the threshold of charged pion pair production, which corresponds to s 3 = 4m 2 π + . The factor f (U) takes into account the additional contribution from π + π − bound states and other narrow peaks from electromagnetic effects, all decaying to π 0 π 0 [19]. All these contributions have widths that are much narrower than our experimental M 2 00 mass resolution.
The s 3 bin width used to store the measured spectrum is denoted as w. With this definition the parameter p is dimensionless, and represents the relative increase of the content of the bin containing the value U = U t with respect to the value calculated with f (U) = 1. In our analysis we use w = 0.00015 (GeV/c 2 ) 2 , and all the p values listed below are written for this value of bin width.
The exponent q could be different, in principle, above and below the cusp point, but our fits show that there is no need for such an additional degree of freedom, because the M 2 00 shape in these two regions is successfully described by the two independent constants a and b.
The parameters describing the K ± → π ± π 0 π 0 Dalitz plot are g, h, k, a, b, p, q. The parameters g, h, k are not equivalent to the corresponding constants G, H, K of the old PDG parameterization (1) [4] and of the physical parameterizations [11], but have a similar meaning. The expression (3) is inspired by the Cabibbo-Isidori physical parameterization of the K ± → π ± π 0 π 0 matrix element at tree level [7,8]. The last two terms of (3) correspond to an empirical description of the ππ rescattering effects [8,9,10].

Fitting the data
The V -dependence of expression (2) is described by the kV 2 2 term which is known from earlier measurements to be rather small, k ≈ 0.01 [2,3,11]. So, ignoring the term ∝ k 2 , the U-dependence of the K ± → π ± π 0 π 0 decay width can be expressed as where V max (U) is the maximum kinematically allowed V for a given U. If k is known, formula (5) can be used to fit the K ± → π ± π 0 π 0 decays U-distribution provided the sensitivity of the acceptance to the small kV 2 2 term is taken into account as a contribution to the systematic uncertainty of the results.
The main U-dependence of the parameterization is described by the parameters g, h, a, b, p, q, which are related to the measurement of s 3 , which is equal to the the square of the π 0 π 0 invariant mass, M 2 00 . So the systematic uncertainties of these parameters depend mainly on the performance of the LKr calorimeter. The measurement of k relies also on the measurement of the π ± track in the DCH, but due to the smallness of the k value its uncertainty affects only weakly the determination of the other parameters. Furthermore, a, b, p, q describe the fine features of the Dalitz plot in the cusp region that require narrow s 3 bins, while the k term of formula (5) is smooth over the Dalitz plot and does not require such narrow bins. So we have decided to measure the V -dependence of the Dalitz plot separately by an iterative procedure.
Assuming an initial value k = 0.01, a first fit to the one-dimensional s 3 distribution has been performed using the MINUIT package. The χ 2 was calculated from the difference between the number of observed events in each bin and the number predicted from the parameterization (5) with the current values of the fit parameters. The predicted number of events was calculated by convoluting the parameterization (5) with the MC distributions of the measured s 3 for each generated ('true') s 3 value. In such a way both acceptance and resolution effects were taken into account.
The parameters a, b, p, q, describing the cusp shape, were then fixed to the values obtained from the first fit and used for the two-dimensional fit to determine k. This fit was performed by implementing the event-weighting technique in the MINUIT package. As a first step, the number of events in each bin of the experimental Dalitz plot was corrected for the trigger inefficiency. Then, at each step of the χ 2 minimization the full MC sample corresponding to the experimental data used in the fit (≈ 280 × 10 6 events) was used to build a simulated Dalitz plot by giving each event a weight equal to the ratio between the parameterization (2) with the current values of the fit parameters, and (1), which was used for the simulation of MC events. In the calculation of the weights the 'true' U, V values were used, while the MC events were binned using the reconstructed U, V values (here V means |V |). The MC Dalitz plot was normalized to the total number of data events. The χ 2 was then calculated from the difference between the MC and data Dalitz plots.
For the two-dimensional U, V histograms we used 50 ×50 bins in the intervals −1.45 < U < 1.35 , 0 < V < 2.8. The χ 2 contribution was calculated for the center of each bin over the U range corresponding to the one-dimensional fit limits, and with the V upper limit set to 0.9V max (U) to avoid tails effect.
This fit was performed with a, b, p, q fixed to the values obtained from the onedimensional fit made with formula (5) under the initial assumption k = 0.01. The result of the two-dimensional fit was k = 0.0081 (2). When the procedure was repeated with k = 0.0081 as the initial assumption, it reproduced the measurement k = 0.0081(2) with χ 2 = 1163.5 for 1249 degrees of freedom (probability 0.96), so no further iteration was needed. The fit without the trigger correction gives k = 0.0086 (2), providing an estimate of the trigger inefficiency effect, which is conservatively taken as the contribution to the systematic error on k. Thus, our result for the k parameter of the Dalitz plot is k = 0.0081 ± 0.0005 T rigger ± 0.0002 Stat = 0.0081(5).
(6) Fig. 3 shows the comparison of the experimental and simulated |V | distributions obtained by projection of the the two-dimensional distribution used in the fit to extract the k value (6). Using the fixed k value given in (6), the values of all other parameters in (2) as well as their systematic uncertainties were obtained from the fit to the one-dimensional s 3 distribution, after correcting the content of each bin for the trigger efficiency. The fit gives χ 2 = 265.1 for 284 degrees of freedom (probability 0.78). The best fit values of the parameters are listed in Table 1. The uncertainty affecting the k value is taken into account as one of the sources of systematics errors for the other parameters, and is denoted as k error in Table 1. The effect of the trigger efficiency is also conservatively taken as the contribution to the systematic error for every parameter, and is denoted as Trigger in Table 1.

Systematic uncertainties
All sources of systematic uncertainties are described in detail in ref. [11]. The detector acceptance to K ± → π ± π 0 π 0 decays depends strongly on the position of the K ± decay vertex along the nominal beam axis, Z. A small difference between the shapes of the experimental and simulated distributions is present in the high Z region (close to the spectrometer) where the acceptance drops because of the increasing probability for the charged pion track to cross the spectrometer too close to the event COG. The effect of this difference has been checked by introducing a small mismatch in the track radius cuts between real and simulated samples, and also by applying a small change to the LKr energy scale (that leads to a shift of the measured Z position). The corresponding small changes of the fit results are considered as the acceptance related contribution to the systematic errors (denoted as Acceptance(Z) in Table 1). The simulated sample from which the acceptance and resolution effects used in the fits are derived, is generated under the assumption that the K ± → π ± π 0 π 0 matrix element does not depend on V . We have studied the sensitivity of the fit results to the presence of the V -dependent term compatible with our data in the simulated sample. The largest variations of the fit results are shown in Table 1 as the contributions to the systematic errors arising from the simplified matrix element used in the MC (they are denoted as Acceptance(V)).
The π 0 π 0 invariant mass, M 00 , is determined using only information from the LKr calorimeter. We find that a convenient variable which is sensitive to all random fluctuations of the LKr response, and hence to its energy resolution, is the ratio m π 0 1 /m π 0 2 , where m π 0 1 and m π 0 2 are the measured two-photon invariant masses for the more and less energetic π 0 , respectively, in the the same event. The width of the distribution for simulated events is slightly larger than that of the data: the RMS value of the simulated distribution is 0.0216, while it is 0.0211 for the data.
In order to check the sensitivity of the fit results to a resolution mismatch of this size, we have smeared the measured photon energies in the data by adding a random energy with a gaussian distribution centered at zero and with σ = 0.06 GeV. Such a change increases the RMS value of the m π 0 1 /m π 0 2 distribution from 0.0211 to 0.0224. A fit is then performed for the data sample so modified, and the values of the fit parameters are compared with those obtained using no energy smearing.
The artificial smearing of the photon energies described above introduces random shifts of the fit parameters within their statistical errors. In order to determine these shifts more precisely than allowed by the statistics of a single fits, we have repeated the   Table 1.
In order to study possible non-linearity effects of the LKr calorimeter response to low energy photons, we select π 0 pairs from K ± → π ± π 0 π 0 events with symmetric π 0 → γγ decays (0.45 < E γ /E π 0 < 0.55), and with the more energetic π 0 (denotes as π 0 1 ) in the energy range 22 GeV < E π 0 1 < 26 GeV. For the π 0 pairs so selected we define the ratio of the two-photon invariant masses, r = M π 0 2 /M π 0 1 , where π 0 2 is the lower energy π 0 . Because of the resolution effects discussed above its average value r depends on the lower pion energy even in the case of perfect LKr linearity. However, for E π 0 2 /2 < 9 GeV the values of r for simulated events are systematically above those of the data, providing evidence for the presence of non-linearity effects of the LKr response at low energies.
To study the importance of these effects, we modify all simulated events to account for the observed non-linearity multiplying each photon energy by the ratio r Data / r M C , where r Data and r M C are the average ratios for data and simulated events, respectively. The values of r for the sample of simulated events so modified are very close to those of the data. The small shifts of the best fit parameters obtained using these non-linearity corrections are taken as contributions to the systematic errors in Table 1, where they are denoted as "LKr non-linearity".
The π ± interaction in LKr may produce multiple energy clusters which are located, in general, near the impact point of the π ± track and in some cases may be identified as photons. To reject such "fake" photons a special cut on the distance d between each photon and the impact point of any charged particle track at the LKr is implemented in the event selection. In order to study the effect of these "fake" photons on the best fit parameters we have repeated the fits by varying the cut on the distance d between 10 and 25 cm in the selection of both data and simulated K ± → π ± π 0 π 0 events. The largest deviations from the results obtained with the default cut value (d=15 cm) are taken as contributions to the systematic errors (see Table 1, "Hadronic showers").
The MC program includes a complete simulation of the beam magnet system and collimators with the purpose of predicting the correlation between the incident K ± momenta and trajectories. However, the absolute beam momentum scale cannot be predicted with the required precision, hence we tune the average value to the measured ones for each continuous data taking period ("run") using K ± → π ± π + π − events which are recorded continuously during data taking, and also simulated by the MC program.
After this adjustement, a residual difference still exists between the measured and simulated K ± momentum distributions. In order to study the sensitivity of the best fit parameters to this distribution, we have corrected the simulated momentum distribution to reproduce the measured one. The corresponding changes of the best fit parameters are included in the contributions to the systematic errors and denoted as 'P K spectrum' in Table 1.
In order to take into account variations of running conditions during data taking, the number of simulated K ± → π ± π 0 π 0 events for each run should be proportional to the corresponding number of events in the data. However, because of small variations of trigger efficiency and acceptance, the ratio between the number of simulated and real events varies by a few percent during the whole data taking period. In order to study the effect of the small mismatch between the two samples on the best fit parameters, we have made them equal run by run by a random rejection of selected events. The corresponding shifts of the best fit parameters are considered as a MC time dependent systematic error, and are listed in Table 1, where they are denoted as "MC(T)".
Correlations between the fit parameters are changed by the systematic uncertainties from the values shown in the Table 2 (purely statistical correlations) to the ones of Table  3.

Conclusion
The square of the K ± → π ± π 0 π 0 matrix element can be written using the empirical approximation where w = 0.00015 (GeV/c 2 ) 2 and U t = (4m 2 π + − s 0 )/m 2 π + with the following parameter values: Near the cusp point U = U t this approximation is only valid if the s 3 distribution is averaged over bins which are wider than the intrinsic width of the peak expected from π + π − bound states and other electromagnetic effects [19], all decaying to π 0 π 0 . This peak is much narrower than the bin width used here, w = 0.00015 (GeV/c 2 ) 2 , which is of the order of the experimental resolution.
The errors are dominated by systematic effects. The systematic errors on the slope parameters g, h are substantially larger than the errors on g, h obtained from our study of the ππ scattering lengths based on the full 2003-2004 data sample [11]. This is mainly because we use here the almost full fit interval in order to give a complete description of the K ± → π ± π 0 π 0 Dalitz plot, while the fitting range used in ref. [11] was optimized to reach the smallest total error for the measured ππ scattering lengths. The wide s 3 fitting range increases the sensitivity of the results to LKr non-linearity and trigger inefficiency.
Finally, we note that there is no model-independent relation between the values of the best fit parameters given above and the S-wave ππ scattering lengths a 0 and a 2 , which are meaningful variables only within a specific formulation of ππ rescattering effects in K ± → π ± π 0 π 0 decay (see [11]). The empirical parameterization proposed here provides a good description of this decay mode, but makes no assumption about the physics mechanisms responsible for the observed cusp structure.