Determination of the S-wave pi pi scattering lengths from a study of K+- ->pi+- pi0 pi0 decays

We report the results from a study of the full sample of ~ 6.031 10^7 K+- ->pi+- pi0 pi0 decays recorded by the NA48/2 experiment at the CERN SPS. As first observed in this experiment, the pi0 pi0 invariant mass M00 distribution shows a cusp-like anomaly in the region around M00 = 2m+, where m+ is the charged pion mass. This anomaly has been interpreted as an effect due mainly to the final state charge exchange scattering process pi+ pi- ->p0 p0 in K+- ->pi+- pi+ pi- decay. Fits to the M00 distribution using two different theoretical formulations provide the presently most precise determination of (a0 - a2), the difference between the pi pi S-wave scattering lengths in the isospin I=0 and I=2 states. Higher-order pi pi rescattering terms, included in the two formulations, allow also an independent, though less precise, determination of a2.


Introduction
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 [1,2,3]. The experiment used simultaneous K + and K − beams with momenta of 60 GeV/c propagating through the detector along the same beam line. Data were collected in 2003-2004, providing large samples of fully reconstructed K ± → π ± π + π − and K ± → π ± π 0 π 0 decays. From the analysis of the data collected in 2003, we have already reported the observation of a cusp-like anomaly in the π 0 π 0 invariant mass (M 00 ) distribution of K ± → π ± π 0 π 0 decays in the region around M 00 = 2m + , where m + is the charged pion mass [4]. The existence of this threshold anomaly had been first predicted in 1961 by Budini and Fonda [5], as a result of the charge exchange scattering process π + π − → π 0 π 0 in K ± → π ± π + π − decay. These authors had also suggested that the study of this anomaly, once found experimentally, would allow the determination of the cross-section for π + π − → π 0 π 0 at energies very close to threshold. However, samples of K ± → π ± π 0 π 0 decay events available in those years were not sufficient to observe the effect, nor was the M 00 resolution. As a consequence, in the absence of any experimental verification, the article by Budini and Fonda [5] was forgotten.
More recently, Cabibbo [6] has proposed an interpretation of the cusp-like anomaly along the lines proposed by Budini and Fonda [5], but expressing the K ± → π ± π 0 π 0 decay amplitude in terms of the π + π − → π 0 π 0 amplitude at threshold, a x . In the limit of exact isospin symmetry a x can be written as (a 0 − a 2 )/3, where a 0 and a 2 are the S-wave ππ scattering lengths in the isospin I = 0 and I = 2 states, respectively.

Beam and detectors
The layout of the beams and detectors is shown schematically in Fig. 1. The two simultaneous beams are produced by 400 GeV/c 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% (rms) produced at zero angle are selected by two systems of dipole magnets forming "achromats" with null total deflection, focusing quadrupoles, muon sweepers and collimators. With 7×10 11 protons per pulse 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.40 m for the rest.
A liquid Krypton calorimeter (LKr) [11] 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 single electromagnetic showers can be parameterized as σ x = σ y = 0.42/ √ E ⊕ 0.06 cm for each transverse coordinate x, y.
An additional hodoscope consisting of a plane of scintillating fibers is installed in the LKr calorimeter at a depth of ∼ 9.5 X 0 with the purpose of sampling electromagnetic showers. It is divided into four quadrants, each consisting of eight bundles of vertical fibers optically connected to photomultiplier tubes.

Event selection and reconstruction
The K ± → π ± π 0 π 0 decays are selected by a two level trigger. The first level requires a signal in at least one quadrant of the scintillator hodoscope (Q1) in coincidence with the presence of energy depositions in LKr consistent with at least two photons (NUT). At the second level (MBX), an on-line processor receiving the drift chamber information reconstructs the momentum of charged particles and calculates the missing mass under the assumption that the particle is a π ± originating from the decay of a 60 GeV/c K ± travelling along the nominal beam axis. The requirement that the missing mass is not consistent with the π 0 mass rejects most of the main K ± → π ± π 0 background. The typical rate of this trigger is ∼ 15, 000 per burst.
Events with at least one charged particle track having a momentum above 5 GeV/c, measured with a maximum error of 6% (much larger than the magnetic spectrometer resolution), and at least four energy clusters in the LKr, each consistent, in terms of size and energy, with the electromagnetic shower produced by a photon of energy above 3 GeV, are selected for further analysis. In addition, the relative track and photon timings must be consistent with the same event within 10 ns, and the clusters must be in time between each other within 5 ns.
The distance between any two photons in the LKr is 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 must exceed 15 cm. Fiducial cuts on the distance of each photon from the LKr edges and centre are also applied in order to ensure full containment of the electromagnetic showers. In addition, because of the presence of ∼ 100 LKr cells affected by readout problems ("dead cells"), the minimum distance between the photon and the nearest LKr dead cell is required to be at least 2 cm.
At the following step of the analysis we check the consistency of the surviving events with the K ± → π ± π 0 π 0 decay hypothesis. We assume that each possible pair of photons originates from a π 0 → γγ decay and we calculate the distance D ij between the π 0 decay vertex and the LKr front face: where E i ,E j are the energies of the i-th and j-th photon, respectively, R ij is the distance between their impact points on LKr, and m 0 is the π 0 mass. Among all possible π 0 pairs, only those with D ij values differing by less than 500 cm are retained further, and the distance D of the K ± decay vertex from the LKr is taken as the arithmetic average of the two D ij values. This choice gives the best π 0 π 0 invariant mass resolution near threshold: at M 00 = 2m + it is ∼ 0.56 MeV/c 2 , increasing monotonically to ∼ 1.4 MeV/c 2 at the upper edge of the physical region. The reconstructed distance of the decay vertex from the LKr is further required to be at least 2 m downstream of the final beam collimator to exclude π 0mesons produced from beam particles interacting in the collimator material (the downstream end of the final beam collimator is at Z = −18 m).
Because of the long decay volume, a photon emitted at small angle to the beam axis may cross the aluminium vacuum tube in the spectrometer or the DCH1 central flange, and convert to e + e − before reaching the LKr. In such a case the photon must be rejected because its energy cannot be measured precisely. To this purpose, for each photon detected in LKr we require that its distance from the nominal beam axis at the DCH1 plane must be > 11 cm, assuming an origin on axis at D − 400 cm. In this requirement we take into account the resolution of the D measurement (the rms of the difference between D values for the two photon pairs distribution is about 180 cm).
Each surviving π 0 pair is then combined with a charged particle track, assumed to be a π ± . Only those combinations with a total π ± π 0 π 0 energy between 54 and 66 GeV, consistent with the beam energy distribution, are retained, and the π ± π 0 π 0 invariant mass M is calculated, after correcting the charged track momentum vector for the effect of the small measured residual magnetic field in the decay volume (this correction uses the decay vertex position, D, as obtained from LKr information).
For each π ± π 0 π 0 combination, the energy-weighed average coordinates (center-of-gravity, COG) X COG , Y COG are calculated at each DCH plane using the photon impact points on LKr and the track parameters measured before the magnet (so the event COG is a projection of the initial kaon line of flight). Acceptance cuts are then applied on the COG radial position on each DCH plane in order to select only K ± → π ± π 0 π 0 decays originating from the beam axis. 1 In addition, we require a minimal separation between the COG and the charged track coordinates X t , Y t , as measured in each DCH plane: where the limits depend on the COG and track impact point distributions at each drift chamber (see Table 1). The values of R COG−track min take into account both the beam width (the cut is made with respect to each event COG rather than to the nominal beam center) and the area where the track impact point distribution is still sensitive to the detailed features of the beam shape. In this way the effect of these cuts does not depend strongly on the beam shape and on the precise knowledge of the beam position in space (during data taking, the average beam transverse position was observed to move slightly by up to 2 mm). This cut removes about 28% of events, mainly at large M 2 00 , but the statistical precision of the final results on the ππ scattering lengths is not affected.
For events with more than one accepted track-cluster combination (∼ 1.8% of the total), the K ± → π ± π 0 π 0 decay is selected as the π ± π 0 π 0 combination minimizing a quality estimator based on two variables: the difference ∆D of the two D ij values and the difference ∆M between the π ± π 0 π 0 invariant mass and the nominal K ± mass [12]: where the space and mass resolutions rms D , rms M are functions of D, as obtained from the measured ∆D and ∆M distributions. Fig. 2 shows the distribution of ∆M , the difference between the π ± π 0 π 0 invariant mass and the nominal K ± mass for the selected K ± → π ± π 0 π 0 decays (a total of 6.031 × 10 7 events). This distribution is dominated by the gaussian K ± peak, with a resolution σ = 1.3 MeV/c 2 . There are small non Gaussian tails originating from unidentified π ± → µ ± decay in flight or wrong photon pairing. The fraction of events with wrong photon pairing in this sample is 0.19%, as estimated by the Monte Carlo simulation described in the next Section.  Fig. 2. Distribution of the difference between the π ± π 0 π 0 invariant mass and the nominal K ± mass for the selected K ± → π ± π 0 π 0 decays. Fig. 3 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 centred at M 2 00 = (2m + ) 2 (for most of the physical region the bin width is smaller than the M 2 00 resolution, which is 0.00031 (GeV/c 2 ) 2 at M 2 00 = (2m + ) 2 ). The cusp at M 2 00 = (2m + ) 2 = 0.07792 (GeV/c 2 ) 2 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 [13]. This Monte Carlo (MC) program takes into account all detector effects, including the trigger efficiency and the presence of a small number (< 1%) of "dead" LKr cells. It also includes the simulation of the beam line; the beam parameters are tuned for each SPS burst using fully reconstructed K ± → π ± π + π − events, which provide precise information on the average beam angles and positions with respect to the nominal beam axis. Furthermore, the requirement that the average reconstructed π ± π + π − invariant mass is equal to the nominal K ± mass for both K + and K − fixes the absolute momentum scale of the magnetic spectrometer for each charge sign and magnet polarity, and monitors continuously the beam momentum distributions during data taking.
The Dalitz plot distribution of K ± → π ± π 0 π 0 decays has been generated according to a series expansion in the Lorentz-invariant variable u = (s 3 − s 0 )/m 2 + , where s i = (P K − P i ) 2 (i=1,2,3), s 0 = (s 1 + s 2 + s 3 )/3, P K (P i ) is the K(π) four-momentum, and i = 3 corresponds to the π ± [12]. In our case s 3 = M 2 00 , and s 0 = (m 2 K + 2m 2 0 + m 2 + )/3. 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. The sudden change of slope ("cusp" ) observed in the M 2 00 distribution at M 2 00 = (2m + ) 2 (see Fig. 3) can be interpreted [5] [6] as a threshold effect from the decay K ± → π ± π + π − contributing to the K ± → π ± π 0 π 0 amplitude through the charge exchange reaction π + π − → π 0 π 0 . In the formulation by Cabibbo [6] the K ± → π ± π 0 π 0 decay amplitude is described as the sum of two terms: where M 0 is the tree level K ± → π ± π 0 π 0 weak decay amplitude, and M 1 is the contribution from the K ± → π ± π + π − decay amplitude through π + π − → π 0 π 0 charge exchange, with the normalization condition M 1 = 0 at M 2 00 = (2m + ) 2 . The contribution M 1 is given by where a x is the S-wave π + π − charge exchange scattering length (threshold amplitude), and M + is the K ± → π ± π + π − decay amplitude at M 00 = 2m + . M 1 changes from real to imaginary at M 00 = 2m + with the consequence that M 1 interferes destructively with M 0 in the region M 00 < 2m + , while it adds quadratically above it.
In the limit of exact isospin symmetry a x = (a 0 − a 2 )/3, where a 0 and a 2 are the S-wave ππ scattering lengths in the I = 0 and I = 2 states, respectively. However, it was shown in ref. [4] that a fit of this simple formulation to the NA48/2 M 2 00 distribution in the interval 0.074 < M 2 00 < 0.097 (GeV/c 2 ) 2 using a x m + as a free parameter gave only a qualitative description of the data, with all data points lying systematically above the fit in the region near M 2 00 = (2m + ) 2 . It was also shown in ref.
[4] that a good fit could be obtained using a more complete formulation of ππ final state interaction [7] which took into account all rescattering processes at the one-loop and two-loop level.
In the following sections we present the determination of the ππ scattering lengths a 0 and a 2 by fits of the full data set described in Section 2 to two theoretical approaches: the Cabibbo-Isidori (CI) formulation [7], and the more recent Bern-Bonn (BB) formulation [8].
In the CI approach, the structure of the cusp singularity is treated using unitarity, analiticity and cluster decomposition properties of the S-matrix. The decay amplitude is expanded in powers of ππ scattering lengths up to order (scattering length) 2 , and electromagnetic effects are omitted.
The BB approach uses a non-relativistic Lagrangian framework, which automatically satisfies unitarity and analiticity constraints, and allows one to include electromagnetic contributions in a standard way [9].
In all fits we also need information on the K ± → π ± π + π − decay amplitude. To this purpose, we use a sample of 4.709 × 10 8 K ± → π ± π + π − decays which are also measured in this experiment [14].
four-momentum and i = 3 corresponds to the odd pion (π ± from K ± → π ± π 0 π 0 , π ∓ from K ± → π ± π + π − decay), and v = (s 1 −s 2 )/m 2 + . It must be noted that in ref. [7] the v dependence of both amplitudes had been ignored because the coefficients k 0 and k were consistent with zero from previous experiments. Within the very high statistical precision of the present experiment this assumption is no longer valid.
In the fits to the M 2 00 distribution from K ± → π ± π 0 π 0 decay, the free parameters are (a 0 − a 2 )m + , a 2 m + , g 0 , h 0 , and an overall normalization constant. The coefficient k 0 cannot be directly obtained from a fit to the M 2 00 distribution. Its value is determined independently from the Dalitz plot distribution of K ± → π ± π 0 π 0 decays, as described in the Appendix. The value k 0 = 0.0099 is kept fixed in the fits.
All M + parameters are fixed from data: the coefficients g, h, k are obtained from a separate fit to the K ± → π ± π + π − decay Dalitz plot [14], using M + as given by Eq. (4), and taking into account Coulomb effects; and A + is obtained from the measured ratio, R, of the K ± → π ± π + π − and K ± → π ± π 0 π 0 decay rates, R = 3.175 ± 0.050 [12], which is proportional to A 2 + . The fit gives g = −0.2112 ± 0.0002, h = 0.0067 ± 0.0003, k = −0.00477 ± 0.00008; and we obtain A + = 1.925 ± 0.015. These values are kept fixed in the fits to the M 2 00 distribution from K ± → π ± π 0 π 0 decay. As explained in Section 6 all fits are performed over the M 2 00 interval from 0.074094 to 0.104244 (GeV/c 2 ) 2 (bin 26 to 226). The CI formulation [7] does not include radiative corrections, which are particularly important near M 00 = 2m + , and contribute to the formation of π + π − atoms ("pionium"). For this reason we first exclude from the fit a group of seven consecutive bins centred at M 2 00 = 4m 2 + (an interval of ±0.94 MeV/c 2 in M 00 ). The quality of this fit is illustrated in Fig. 4a, which displays the quantity ∆ ≡ (data -fit)/data as a function of M 2 00 . The small excess of events from pionium formation is clearly visible. Pionium formation and its dominating decay to π 0 π 0 are taken into account in the fit by multiplying the content of the bin centred at M 2 00 = 4m 2 + (bin 51) by 1 + f atom , where 1 + f atom describes the contribution from pionium formation and decay. The pionium width is much narrower than the bin width, since its mean lifetime is measured to be ∼ 3 × 10 −15 s [18]; however, the M 2 00 resolution is taken into account in the fits as described in the last paragraph of Section 3. The results of a fit with f atom as a free parameter and with no excluded bins near M 2 00 = 4m 2 + are given in Tables 2 and 3 (fit CI): the quality of this fit is shown in Fig. 4b. The best fit value f atom = 0.0533 ± 0.0091 corresponds to a rate of K ± → π ± + pionium decay, normalized to the K ± → π ± π + π − decay rate, of (1.69 ± 0.29) × 10 −5 , which is larger than the predicted value ∼ 0.8 × 10 −5 [19,20]. As discussed in Section 5, this difference is due to additional radiative effects, which are not taken into account in the CI formulation [7] and, contrary to pionium formation and decay, affect more than one bin. For this reason for the fits without the radiative effects taken into account we prefer to fix f atom = 0.0533 and to exclude from the fit the seven consecutive bins centred at M 2 00 = 4m 2 + . The results of this fit are listed as Fit CI A in Tables 2 and 3.
We have also performed fits using the constraint between a 2 and a 0 predicted by analyticity and chiral symmetry [21] (we refer to this constraint as the ChPT constraint): The results of these fits are shown in Tables 2 and 3 (fits CI χ and CI χ A ). For fit CI χ no bins near the cusp point are excluded and f atom is a free parameter, while for fit CI χ A the seven bins centred at M 2 00 = 4m 2 + are excluded and f atom is kept fixed at the value obtained from fit CI χ .

Fits using the Bern-Bonn theoretical formulation
The Bern-Bonn (BB) formulation [8] describes the K → 3π decay amplitudes using two expansion parameters: a, the generic ππ scattering amplitude at threshold; and a formal parameter ǫ such that in the K-meson rest frame the pion momentum is of order ǫ, and its kinetic energy T is of order ǫ 2 . In the formulation of ref. [8] the K → 3π decay amplitudes include terms up to O(ǫ 2 , aǫ 3 , a 2 ǫ 2 ). However, in the formulae used in the fits described below these amplitudes include terms up to O(ǫ 4 , aǫ 5 , a 2 ǫ 2 ). In the BB formulation the description of the K → 3π decay amplitudes is valid over the full physical region 2 .
At tree level the K → 3π decay amplitudes are expressed as polynomials containing terms in T 3 , T 2 3 , and (T 1 − T 2 ) 2 , where T 3 is the kinetic energy of the "odd" pion (π ± from K ± → π ± π 0 π 0 , π ∓ from K ± → π ± π + π − decay) in the K ± rest frame, while T 1 and T 2 are the kinetic energies of the two same-sign pions. Since these variables can be expressed as functions of the relativistic invariants u and v defined previously, for consistency with the fits described in the previous subsection we prefer to use the same forms as given in Eqs. (3) and (4). It must be noted, however, that the best fit polynomial coefficients are not expected to be equal to those obtained from the fits to the CI formulation [7] because the loop diagram contributions are different in the two formulations.
As for CI, also in the BB formulation rescattering effects are much smaller in K ± → π ± π + π − than in the K ± → π ± π 0 π 0 decay, and a good fit to the M 2 ±± distribution alone can be obtained with or without the addition of rescattering terms to the tree-level weak amplitude of K ± → π ± π + π − decay. However, contrary to CI, the coefficients of the tree-level K ± → π ± π + π − amplitudes enter into the K ± → π ± π 0 π 0 rescattering terms in different combinations. Therefore, the use of a phenomenological description of the K ± → π ± π + π − decay amplitude extracted from a fit to K ± → π ± π + π − data alone is not justified in this case. Thus, in order to obtain a precision on the fit parameters which matches the BB approximation level, the value of each coefficient of the K ± → π ± π + π − tree-level amplitude is obtained from the fit. 3 We perform simultaneous fits to two distributions: the M 2 00 distribution described in Section 2 and the M 2 ±± distribution from K ± → π ± π + π − decay, obtained as a projection of the Dalitz plot described in ref. [14]. This latter distribution is made with the same binning as for the M 2 00 distribution from K ± → π ± π 0 π 0 decay and consists of 4.709 × 10 8 events.
All fits are performed over the M 2 00 interval from 0.074094 to 0.104244 (GeV/c 2 ) 2 (bin 26 to 226), and from 0.080694 to 0.119844 (GeV/c 2 ) 2 (bin 70 to 330) for the M 2 ±± distribution from K ± → π ± π + π − decay. As for the M 2 00 distribution from K ± → π ± π 0 π 0 decay, a very large sample of simulated K ± → π ± π + π − decays (see ref. [14]) is used to obtain the detection probability and the distribution function for the reconstructed value M 2 ±± for any generated value of M 2 ±± . In all fits the free parameters are (a 0 − a 2 )m + and a 2 m + (or only a 0 m + for the fit using the ChPT constraint given by Eq. (5)), the coefficients of the tree-level weak amplitudes g 0 , h 0 , g, h, k (see Eqs. (3,4)), and two overall normalization constants (one for each distribution). The coefficient k 0 (see Eq. (3)) is determined independently from a separate fit to the Dalitz plot distribution of K ± → π ± π 0 π 0 decays (see the Appendix). The fixed value k 0 = 0.0085 is used in the fits. In some of the fits the contribution from pionium formation, described by f atom , is also a free parameter.
Since the detection of K ± → π ± π 0 π 0 and K ± → π ± π + π − decays involves different detector components and different triggers (no use of LKr information is made to select K ± → π ± π + π − decays), the ratio of the detection efficiencies for the two decay modes is not known with the precision needed to extract the value of A + (see Eq. (4)) from the fit. Therefore, as for the CI fits, also for the BB fits A + is obtained from the ratio of the K ± → π ± π + π − and K ± → π ± π 0 π 0 decay rates, measured by other experiments, R = 3.175 ± 0.050 [12]. Tables 2 and 3 show the results of a fit (fit BB) using f atom as a free parameter and including all bins around the cusp point in the fit; for fit BB A the value of f atom is fixed and seven bins centred at M 2 00 = 4m 2 + are excluded. A comparison with the results of the corresponding CI fits (fits CI and CI A , respectively) shows that the difference between the best fit values of (a 0 − a 2 )m + is rather small (about 3%), while the difference between the two a 2 m + values is much larger. We note that in the BB fits a 2 m + has a stronger correlation with other fit parameters than in the CI fits (see Tables 4 and 5). Fits BB χ and BB χ A (see Tables 2 and 3) are similar to BB and BB A , respectively, but the ChPT constraint given by Eq. (5) is used. Here the best fit value of a 0 m + agrees well with the value obtained from the CI fit (fit CI χ A ).

Radiative correction outside the cusp point
Radiative corrections to both K ± → π ± π 0 π 0 and K ± → π ± π + π − decay channels have been recently studied by extending the BB formulation [8] to include real and virtual photons [9]. In the K ± rest frame the emission of real photons is allowed only for photon energies E < E cut . We have performed simultaneous fits to the M 2 00 distribution from K ± → π ± π 0 π 0 and to the M 2 ±± distribution from K ± → π ± π + π − decays using the formulation  of ref. [9]. Our event selection does not exclude the presence of additional photons; however, energetic photons emitted in K ± decays result in a reconstructed π ± π 0 π 0 invariant mass lower than the K mass. We set E cut = 0.010 GeV in order to be consistent with the measured π ± π 0 π 0 invariant mass distribution shown in Fig. 2 (the same is true for the π ± π + π − invariant mass distribution from K ± → π ± π + π − decay measured in this experiment [14]). For each fit we adjust the value of A + (see Eq. (4)) so that the ratio of the K ± → π ± π + π − and K ± → π ± π 0 π 0 decay rates is consistent with the measured one [12]. The formulation of ref. [9] does not include pionium formation, and the K ± → π ± π 0 π 0 amplitude, A rad 00+ , has a non-physical singularity at M 2 00 = (2m + ) 2 . To avoid problems in the fits, the square of decay amplitude at the center of bin 51, where the singularity occurs, is replaced by |A 00+ | 2 (1 + f atom ), where A 00+ is the decay amplitude of the BB formulation without radiative corrections [8], and f atom is again a free parameter.
The results of simultaneous fits to the M 2 00 distribution from K ± → π ± π 0 π 0 decays, and to the M 2 ±± distribution from K ± → π ± π + π − decay are shown in Tables 6 and 7. In all these fits the M 2 00 and M 2 ±± intervals are equal to those of the fits described in Sections 4.1 and 4.2 (see Tables 2 and 3). In fit BB all bins around the cusp point are included and f atom is a free parameter, while in fit BB A seven consecutive bins centred at M 2 00 = (2m + ) 2 are excluded and f atom is fixed to the value given by fit BB. A comparison of fit BB or BB A with radiative corrections taken into account (Table 6) with the corresponding fits without radiative corrections (fits BB, BB A of Table  2) shows that radiative corrections reduce (a 0 − a 2 )m + by ∼ 9%. However, the change in the best fit value of a 2 m + is much larger, possibly suggesting again that the determination of this scattering length is affected by large theoretical uncertainties.
Fits BB χ and BB χ A in Tables 6 and 7 are similar to BB and BB A , respectively, but the constraint between a 2 and a 0 predicted by analyticity and chiral symmetry [21] (see Eq. (5)) is used. A comparison of fits BB χ and BB χ A with the corresponding fits obtained without radiative corrections (fits BB χ , BB χ A of Table 2) shows that radiative corrections reduce a 0 m + by ∼ 6%.
For all fits BB χ to BB χ A in Tables 6 and 7 the effect of changing the maximum allowed photon energy E cut from 0.005 to 0.020 GeV is found to be negligible.
No study of radiative corrections has been performed in the framework of the CI approach [7]. However, the dominating radiative effects (Coulomb interaction and photon emission) are independent of the specific approximation. Therefore, extracting the relative effect of radiative corrections from the BB calculation and using it for the fit to the CI formula is justified. In order to obtain an approximate estimate of radiative effects in this case, we have corrected the fit procedure by multiplying the absolute value of the K ± → π ± π 0 π 0 decay amplitude given in ref. [7] by |A rad 00+ /A 00+ | [22], as obtained in the framework of the BB formulation [8,9]. Because of the non-physical singularity of A rad 00+ at M 2 00 = (2m + ) 2 in the BB formulation, in the calculation of the K ± → π ± π 0 π 0 decay amplitude for the 51 st bin we also multiply the squared amplitude of ref. [7] by 1 + f atom .
The results of these radiative-corrected fits to the M 2 00 distribution from K ± → π ± π 0 π 0 decay performed using   Table 8. Fit parameter correlations for the CI formulation with radiative correction (fit CI in Table 6).  Table 9. Fit parameter correlations for the BB formulation with radiative correction (fit BB in Table 6).  Tables 8 and  9. Fig. 5 illustrates the fit results for the fits CI and BB with and without radiative corrections. All the fits are performed using the same K ± → π ± π 0 π 0 data sample.

Pionium formation and other electromagnetic effects at the cusp point
Pionium formation in particle decay and in charged particle scattering was studied in early theoretical work [20,23], but a unified description of its production together with other electromagnetic effects near threshold was missing.
In a more recent approach [24], electromagnetic effects in K ± → π ± π 0 π 0 decay have been studied in the framework of nonrelativistic quantum mechanics using a potential model to describe the electromagnetic interaction between the π + π − pair in loop diagrams. This model is equivalent to a perturbative one, in which all simple sequential π + π − loops with electromagnetic interactions between the two charged pions are taken into account to all orders (including the formation of electromagnetically bound final states), but there is no emission of real photons and the electromagnetic interaction with the other π ± from the K ± → π ± π + π − decay is ignored. Because of these limitations, the model of ref. [24] cannot be directly applied to the full physical region of the K ± → π ± π 0 π 0 decay; however, contrary to the BB formulation [9], its integral effect over a narrow region which includes the cusp point (M 2 00 = 4m 2 + ) can be calculated. We have implemented the electromagnetic effects predicted by the model of ref. [24] in the parameterization of the CI formulation [7] (the detailed procedure is described in Eqs. (6,7,8) of ref. [25]). In the theoretical M 2 00 distribution the electromagnetic correction for the bin centred at 4m 2 + (bin 51), averaged over the bin, depends on the bin width, as it includes contributions from both pionium bound states with negligible widths and a very narrow peak of unbound π + π − states annihilating to π 0 π 0 . For the bin width of 0.00015 (GeV/c 2 ) 2 used in the fits, these effects increase the content of bin 51 by 5.8%, in agreement with the results of the fits performed using f atom as a free parameter (see Tables 2, 6). Thus the model of ref. [24] explains why the typical fit result for f atom is nearly twice as large as the prediction for pionium contribution only, as calculated in refs. [19,20].
Near the cusp point the two calculations of electromagnetic effects [9] and [24,25] are very similar numerically, thus increasing the confidence in the central cusp bin radiative effect calculated using Eq. (8) of ref. [25]. However, at larger distances from the cusp the approach of refs. [24,25] leads to deviations from the electromagnetic corrections of ref. [9]. This can be explained by the fact that the model of ref. [24] takes into account only processes that dominate near the cusp point. For this reason we do not use this model in the fits, but we consider it as a complementary calculation limited to a region very close to the cusp point, providing a finite result for the bin centred at M 2 00 = 4m 2 + which the formulation of ref. [9] does not provide.

Systematic uncertainties
As shown below, all systematic corrections affecting the best fit values of the coefficients describing the K ± → π ± π 0 π 0 weak amplitude at tree level, g 0 and h 0 (see Eq. (3)), are found to be much smaller than the statistical errors. We use these corrections as additional contributions to the systematic uncertainties instead of correcting the central values of these parameters.
For a given fit, we find that the systematic uncertainties affecting the best fit parameters do not change appreciably if the fit is performed with or without electromagnetic corrections. In addition, we find that, with the exception of f atom , the systematic uncertainties affecting all other parameters are practically the same if in the fit the seven consecutive bins centred at M 2 00 = 4m 2 + are included (and f atom is used as a free parameter), or if they are excluded (and the value of f atom is fixed).
For these reasons, we give detailed estimates of the systematic uncertainties only for fits CI, CI χ , BB, BB χ performed with the decay amplitude corrected for electromagnetic effects.
The parameters g, h, k which describe the K ± → π ± π + π − weak amplitude at tree level are used as free parameters when fitting the data to the BB formulation [8,9]. However, they enter into the K ± → π ± π 0 π 0 decay amplitude only through rescattering terms, thus we do not consider the best fit values of these parameters as a measurement of physically important values. Here we do not estimate the systematic uncertainties affecting them and we discuss the uncertainties associated with K ± → π ± π + π − decay in Section 7. In the study of the systematic uncertainties affecting the K ± → π ± π 0 π 0 decay parameters we fix the values of the K ± → π ± π + π − decay parameters g, h, k in the BB formulation to their best fit values shown in Table 7.
The fit interval for the presentation of the final results (bins 26-226 of width 0.00015 (GeV/c 2 ) 2 , with bin 51 centred at 4m 2 π + ) has been chosen to minimize the total experimental error of the measured a 0 − a 2 . If the upper limit of the fit region, s max 3 , is increased, the statistical error decreases. All our fits give good χ 2 up to rather high s max 3 values where the acceptance is small 4 . However, the systematic error increases with s max 3 , especially the contributions from trigger inefficiency and non-linearity of the LKr response. The total experimental error on a 0 − a 2 , obtained by adding quadratically the statistical and systematic error, has a minimum when the upper limit of the fit interval corresponds to bin 226.

Acceptance
The detector acceptance to K ± → π ± π 0 π 0 decays depends strongly on the position of the K ± decay vertex along the nominal beam axis, Z, so the Z distribution provides a sensitive tool to control the quality of the acceptance simulation. Fig. 6 shows the comparison between the data and Monte-Carlo simulated Z distributions. The small difference between the shapes of the two distributions in the region Z < 0 disappears when the trigger efficiency correction is applied, so this difference is taken into account in the contribution to the systematic uncertainties from the trigger efficiency (see . A small difference between the shapes of the two distributions is also present in the large Z region in the area 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 acceptance difference has been checked by introducing a small mismatch in the track radius cuts between real and simulated data, and also by applying small changes to the LKr energy scale (equivalent to shifts of the event Z position similar to the effect observed in the acceptance). The corresponding small changes of the fit results are considered as the acceptance related contribution to the systematic uncertainties (quoted as Acceptance(Z) in Tables  11-14).  The Monte Carlo 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, M, depends only on u. We have stud-ied the sensitivity of the fit results to the presence of a v-dependent term by adding to |M| 2 a term of the form k 0 v 2 or k ′ Re(M)v 2 , consistent with the observed v dependence in the data. The largest variations of the fit results are shown in Tables 11-14 as the contributions to the systematic uncertainties arising from the simplified matrix element used in the Monte Carlo (they are quoted as Acceptance(V)).

Trigger efficiency
During data taking in 2003 and 2004 some changes to the trigger conditions were introduced following improvements in detector and electronics performance. In addition, different minimum bias triggers with different downscaling factors were used. As a consequence, trigger effects have been studied separately for the data samples taken during seven periods of uniform trigger conditions. Details of the trigger efficiency for the K ± → π ± π 0 π 0 decay events are given in [1,3].
As described in Section 2, 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 as described in section 2, it was possible to measure separately two efficiencies: 1. the efficiency of the minimum bias Q1*MBX trigger using the event sample recorded by the minimum bias NUT trigger; 2. 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.
The measured efficiencies for seven different periods are shown in Fig. 7 as a function of the reconstructed M 2 00 . In the initial data taking periods the samples of minimum bias events were rather small, resulting in relatively large statistical errors. However, we can improve the estimate of the trigger efficiency for these periods under the additional assumption that it is a smooth function of M 2 00 (this assumption is justified by the fact that no anomaly is expected nor observed in its behaviour). We find that a 2-nd degree polynomial describes well the trigger efficiency over the M 2 00 fit interval. Moreover, over this interval the dependence is almost linear, so we expect a negligible effect on the determination of the scattering lengths.  Fits are made separately for each of the data taking periods shown in Fig. 7. In a first fit, the M 2 00 distribution from the data and the corresponding trigger efficiency are fitted simultaneously, and the theoretical M 2 00 distribution, distorted by the acceptance and resolution effects, is multiplied by the corresponding trigger efficiency, as parameterized using Eq. (6). The fit to the M 2 00 distribution alone is then repeated under the assumption of a fully efficient trigger, and the results of the two fits are compared to obtain the trigger efficiency correction and its effective error. As an example, Table 10 lists the trigger corrections to the best fit parameters of fits CI and CI χ (see Table 6).
The trigger corrections are all in agreement with zero within their statistical uncertainties. For a conservative estimate, we combine in quadrature the corrections and their errors to obtain the trigger efficiency contribution to the systematic uncertainties of the best fit results (see Tables 11-14).

LKr resolution
As described in Section 2, the π 0 π 0 invariant mass M 00 is determined using only information from the LKr calorimeter (photon energies and coordinates of their impact points). The measurement of the scattering lengths relies, therefore, on the correct description of the M 00 resolution in the Monte Carlo simulation.
In order to check the quality of the LKr energy resolution we cannot use the π 0 mass peak in the two-photon invariant mass distribution, because the nominal π 0 mass [12] is used in the reconstruction of the two-photon decay vertex (see Section 2). 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 twophoton invariant masses for the more and less energetic π 0 , respectively, in the same K ± → π ± π 0 π 0 decay. The distributions of this ratio for real and simulated events are shown in Fig. 8. One can see that 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 centred at zero and with σ = 0.06 GeV (see Fig. 8). 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 fit, we have repeated the fit eleven times using for each fit a data sample obtained by smearing the original photon energies with a different series of random numbers, as described in the previous paragraph. The shifts of the fit parameters, averaged over the eleven fits, represent the systematic effects, while the errors on those average values are the corresponding uncertainties. Conservatively, the quadratic sum of the shifts and their errors is quoted as "LKr resolution" in Tables 11-14.

LKr non-linearity
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 using the following criteria: 1. both π 0 → γγ decays must be close to symmetrical (0.45 < Eγ E π 0 < 0.55); 2. the more energetic π 0 (denoted as π 0 1 ) must fulfil the requirement 22 GeV < E π 0 1 < 26 GeV. For the π 0 pairs selected in such way 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 . Fig. 9 shows the average ratio r as a function of E π 0 2 /2 for both data and simulated events (for symmetric π 0 → γγ decays E π 0 2 /2 is the photon energy).
Because of the resolution effects discussed in the previous subsection 5 , r depends on the lowest pion energy even in the case of perfect LKr linearity. However, as shown in Fig. 9, 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 MC , where r Data and r MC are the average ratios for data and simulated events, respectively. As shown in Fig. 9, 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 uncer- 5 The small resolution mismatch between data and simulated events introduces a negligible effect here. Tables 11-14, where they are quoted as "LKr non-linearity". versus E π 0 2 /2 for π 0 pairs from K ± → π ± π 0 π 0 decays selected as described in the text. Solid circles: data; crosses: simulated events; open circles: simulated events corrected for non-linearity (see text). The π 0 2 energy is divided by 2 to compare with the γ energy for symmetric π 0 decays.

Hadronic showers in LKR
The π ± interaction in the 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 cut on the distance d between each photon and the impact point of any charged particle track at the LKr front face is implemented in the event selection, as described in Section 2. 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 uncertainties (see Tables 11-14).

Other sources
The Monte Carlo program includes a complete simulation of the beam magnet system and collimators with the purpose of reproducing the correlation between the incident K ± momenta and trajectories. However, the absolute beam momentum scale cannot be modelled 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 during data taking, and also simulated by the Monte Carlo program.
After this adjustment, a residual systematic difference still exists between the measured and simulated K ± momentum distributions, as shown in Fig. 10. In order to study the sensitivity of the best fit parameters to this distribution, we have corrected the width of the simulated K ± momentum distribution to reproduce the measured distribution (see Fig. 10) using a method based on the rejection of simulated events. To minimize the random effect of this rejection, a fraction of events has also been removed from the uncorrected MC sample in such a way that the corrected and uncorrected MC samples have a maximum overlap of events and the same statistics. The corresponding changes of the best fit parameters are included in the contributions to the systematic uncertainties and quoted as "P K spectrum" in  In order to take into account changes 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 changes in the trigger efficiency and in acceptance related to minor hardware problems, 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 Monte Carlo time dependent systematic error, and are listed in Tables 11-14, where they are quoted as "MC(T)". Table 11. Fit parameter systematic uncertainties in units of 10 −4 for the CI formulation with electromagnetic corrections (fit CI in Table 6). The factor m+ which should multiply the scattering lengths is omitted for simplicity.

Source
g0 h0 a0 a2 a0 − a2 fatom Acceptance(Z) 22 Table 12. Fit parameter systematic uncertainties in units of 10 −4 for the CI formulation with electromagnetic corrections and with the ChPT constraint (fit CI χ in Table 6). The factor m+ which should multiply the scattering lengths is omitted for simplicity.

External uncertainties
The most important source of external error is the value of |A + |, obtained from the measured ratio of the K ± → π ± π + π − and K ± → π ± π 0 π 0 decay rates, R = 3.175 ± 0.050 [12]. This ratio is proportional to |A + | 2 , so δ|A + |/|A + | = 0.5(δR)/R.  Table 6). The factor m+ which should multiply the scattering lengths is omitted for simplicity.  Table 6). The factor m+ which should multiply the scattering lengths is omitted for simplicity.
Source g0 h0 a0 a2 a0 − a2 fatom Acceptance(Z) 24  The typical |A + | uncertainty is, therefore, δ|A + | ≈ 0.015. We have checked the shifts of the fit results due to the variation of |A + | within its uncertainty. Each fit is redone twice changing the |A + | value by +δ|A + | and −δ|A + |. One half of the variation of the fit parameters corresponding to these two fits is listed in Table 15, and is taken as the external contribution to the full parameter uncertainty. The BB formulation with radiative corrections [9] provides presently the most complete description of rescattering effects in K → 3π decay. For this reason we use the results from the fits to this formulation to present our final results on the ππ scattering lengths: (a 0 − a 2 )m + = 0.2571 ± 0.0048(stat.) ±0.0025(syst.) ± 0.0014(ext.); a 2 m + = −0.024 ± 0.013(stat.) ±0.009(syst.) ± 0.002(ext.).
The values of the ππ scattering lengths, (a 0 − a 2 )m + and a 2 m + , are obtained from fit BB of Table 6. In addition to the statistical, systematic and external errors discussed in the previous sections, these values are affected by a theoretical uncertainty. We note that, at the level of approximation of the BB and CI amplitude expression used in the fits, a difference of 0.0088(3.4%)is found between the values of (a 0 − a 2 )m + and of 0.015(62%) for a 2 m + . For the sake of comparison with other independent results on the ππ scattering lengths we take into account these differences as theoretical uncertainty. From the measurement of the lifetime of pionium by the DIRAC experiment at the CERN PS [18] a value of |a 0 − a 2 |m + = 0.264 +0.033 −0.020 was deduced which agrees, within its quoted uncertainty, with our result (it should be noted that this measurement provides only a determination of |a 0 − a 2 |, while our measurement of K ± → π ± π 0 π 0 decay is also sensitive to the sign).
Previous determinations of the ππ scattering lengths have also relied on the measurement of K ± → π + π − e ± ν e (K e4 ) decay. Fig. 11 compares our results (Eqs. (7,8)) with the results from the most recent analysis of a large sample of K e4 decays, also collected by the NA48/2 collaboration [26].
For this fit the theoretical uncertainty affecting the value of a 0 − a 2 is estimated to be ±2% (±0.0053) from a recent study of the effect of adding three-loop diagrams to the K ± → π ± π 0 π 0 decay amplitude [27] in the frame of the CI formulation [7] (the goals of this study included a more precise estimate of the theoretical uncertainties affecting the ππ scattering lengths). This theoretical uncertainty is smaller than that affecting the result of the fit with a 0 − a 2 and a 2 as free parameters, because the theoretical uncertainty on a 2 becomes negligible when using the ChPT constraint. The 68% confidence level ellipse corresponding to the result given by Eq. (9) is also shown in Fig. 11, together with a fit to the K e4 data which uses the same ChPT constraint. The a 0 −a 2 vs a 2 correlation coefficient for this figure has been calculated taking into account statistical, systematic and external covariances. Its value is −0.774, while the statistical correlation alone is −0.839 (see Table  9).

Summary and conclusions
We have studied the π 0 π 0 invariant mass distribution measured from the final sample of 6.031 × 10 7 K ± → π ± π 0 π 0 fully reconstructed decays collected by the NA48/2 experiment at the CERN SPS. As first observed in this experiment [4], this distribution shows a cusp-like anomaly at M 00 = 2m + which is interpreted as an effect due mainly to the final state charge-exchange scattering process π + π − → π 0 π 0 in K ± → π ± π + π − decay [5,6]. Good fits to the M 2 00 distribution have been obtained using two different theoretical formulations [7] and [8,9], all including next-to-leading order rescattering terms. We use the results of the fit to the formulation which includes radiative corrections [9] to determine the difference a 0 −a 2 , which enters in the leading-order rescattering term, and a 2 , which enters in the higher-order rescattering terms, where a 0 and a 2 are the I = 0 and I = 2 S-wave ππ scattering lengths, respectively. These values are given in Eqs. (7) and (8), while Eq. (9) gives the result from a fit that uses the constraint between a 2 and a 0 predicted by analyticity and chiral symmetry [21] (see Eq. (5)).
As discussed in Section 8, our results agree with the values of the ππ scattering lengths obtained from the study of K e4 decay [26], which have errors of comparable magnitude. The value of a 0 −a 2 as quoted in Eqs. (7) and (9) are also in agreement with theoretical calculation performed in the framework of Chiral Perturbation Theory [28,29], which predict (a 0 − a 2 )m + = 0.265 ± 0.004.
We finally note a major difference between K ± → π ± π + π − and K ± → π ± π 0 π 0 decays. In the case of K ± → π ± π + π − decay there is no cusp singularity in the physical region because the invariant mass of any pion pair is always ≥ 2m + . As a consequence, rescattering effects can be reabsorbed in the values of the Dalitz plot parameters g, h, k obtained from fits without rescattering, such as those discussed in ref. [14]. On the contrary, a correct description of the K ± → π ± π 0 π 0 Dalitz plot is only possible if rescattering effects are taken into account to the next-to-leading order. Furthermore, the values of the parameters g 0 , h 0 , k 0 which describe the weak K ± → π ± π 0 π 0 amplitude at tree level depend on the specific theoretical formulation of rescattering effects used to fit the data.
In a forthcoming paper we propose an empirical parameterization capable of giving a description of the K ± → π ± π 0 π 0 Dalitz plot, which does not rely on any ππ rescattering mechanisms, but nevertheless reproduces the cusp anomaly at M 00 = 2m + . This parameterization is useful for computer simulations of K ± → π ± π 0 π 0 decay requiring a precise description of all Dalitz plot details.
M 2 00 and cos(θ), where θ is the angle between the momentum vectors of the π ± and one of the two π 0 in the rest frame of the π 0 pair (with this choice of variables the Dalitz plot has a rectangular physical boundary). The M 2 00 fit interval is identical to the one used for the onedimensional fits described in Sections 4.1, 4.2, but the bin width is increased from 0.00015 to 0.0003 (GeV/c 2 ) 2 , and four consecutive bins around M 2 00 = 4m 2 + are excluded. The cos(θ) variable is divided into 21 equal bins from −1.05 to 1.05, but only the interval −0.85 < cos(θ) < 0.85 (17 bins) is used in the fits.
In order to take into account the distortions of the theoretical Dalitz plot due to acceptance and resolution effects, a four-dimensional matrix (with dimensions 210 × 21×210×21) is obtained from the Monte Carlo simulation described in Section 3. This matrix is used to transform the true simulated Dalitz plot into an expected one which can be directly compared with the measured Dalitz plot at each step of the χ 2 minimization.
Fits to the CI formulation [7] are performed with a fixed value a 2 = −0.044. If the k 0 parameter is kept fixed at zero, the fit quality is very poor (χ 2 = 4784.4 for 1237 degrees of freedom); however, if k 0 is used as a free parameter in the fit, the best fit value is k 0 = 0.00974 ± 0.00016, and χ 2 = 1223.5 for 1236 degrees of freedom. The results of these two fits are shown in Fig. 12, where the data and best fit Dalitz plots are projected onto the cos(θ) axis.
A simultaneous fit to the Dalitz plot from K ± → π ± π 0 π 0 decay and to the M 2 ±± distribution from K ± → π ± π + π − decay is performed in the frame of the BB formulation [8] using the constraint between a 2 and a 0 predicted by analyticity and chiral symmetry (see Eq. (5)). The best fit gives k 0 = 0.00850±0.00014, with χ 2 = 1975.5 for 1901 degrees of freedom. The difference between the k 0 value so obtained and that obtained from a fit to the CI formulation [7] is due to the rescattering contributions which are different in the two formulations. When radiative corrections are included in the fit [9], k 0 is practically unchanged (its best fit value is 0.008495), demonstrating that electromagnetic corrections have a negligible effect on its determination.
The second fitting method is based on the event weighting technique. In order to study the size of the trigger effect on the fit parameters, we use a fraction of the data taken with uniform trigger conditions and associated with a large minimum bias event sample which allows a precise evaluation of the trigger efficiency.
The Dalitz plot is described by the u and |v| variables (see Eq. (3)), and the intervals −1.45 < u < 1.35 and |v| < 2.8 are each sudivided into 50 equal size bins. The fits are performed using the CI formulation [7] over a wide region which excludes only the tails of the distribution (0 < |v| < 0.9 v max , u < 0.9). All bins around the cusp point are included, and pionium formation is taken into account by multiplying the theoretical K ± → π ± π 0 π 0 decay probability by the factor 1.055 in the interval |M 2 00 − 4m 2 + | < 0.000075 (GeV/c 2 ) 2 . The fits are performed with a fixed value a 2 = −0.044.
In the fits we use the Dalitz plots distributions of the selected events, corrected (or not corrected) for the trigger efficiency, and of a corresponding subsample of ∼ 2.8×10 7 simulated events generated with a simple matrix element M sim without rescattering effects and with fixed values of g 0 , h 0 and k 0 . At every iteration in the χ 2 minimization, each simulated event is reweighted by the ratio |M| 2 |Msim| 2 , where M is the matrix element which includes rescattering and is calculated with the new fitting parameters, and both M and M sim are calculated at the generated u, |v| values. The simulated events so weighted are then rebinned, and their two-dimensional u, |v| distribution is compared with that of the data.
A good fit (χ 2 = 1166 for 1257 degrees of freedom) is obtained when the trigger efficiency is taken into account, giving k 0 = 0.00966 ± 0.00018. If the trigger effect is ignored, the χ 2 value is somewhat worse (χ 2 = 1276) and we obtain k 0 = 0.01010 ± 0.00017. This result demonstrates that the trigger effect is important for the wide region of the Dalitz plot used in the fit, increasing the measured k 0 by ≈ 0.0004.
The data used in these fits overlap only partially with the data used in the fit to the CI formulation [7] performed using the first method and discussed above, but the results have almost equal statistical errors. We average the two results from the fits without trigger correction, obtaining k 0 = (0.00974 + 0.01010)/2 = 0.0099. We take the statistical error of one of them as the statistical error of the measured k 0 value, and conservatively take one half of the difference between them as the contribution to the systematic error due to the different fitting techniques. As mentioned above, the trigger correction shifts the k 0 central value by −0.0004. Because this effect is measured only with a partial data sample, we also add it in quadrature to the systematic error. So our measurement of k 0 in the frame of the CI rescattering formulation [7] gives k 0 = 0.0095 ± 0.00017(stat.) ± 0.00048(syst.) = 0.0095 ± 0.0005.
For most of the one-dimensional fits discussed in the present paper we do not apply any trigger correction, so here we use the effective value k 0 = 0.0099 for the fits to the CI formulation [7], and k 0 = 0.0085 for the fits to the BB formulation [8,9]. Since k 0 is kept fixed in those fits, we check the variations of all the best fit parameters by varying k 0 within the limits defined by its full error. These variations are listed in Tables 11-14, where they are denoted as "k 0 error".