Precision measurement of the $\eta\to\pi^+\pi^-\pi^0$ Dalitz plot distribution with the KLOE detector

Using $1.6$ fb$^{-1}$ of $e^+ e^-\to\phi\to\eta\gamma$ data collected with the KLOE detector at DA$\Phi$NE, the Dalitz plot distribution for the $\eta \to \pi^+ \pi^- \pi^0$ decay is studied with the world's largest sample of $\sim 4.7 \cdot 10^6$ events. The Dalitz plot density is parametrized as a polynomial expansion up to cubic terms in the normalized dimensionless variables $X$ and $Y$. The experiment is sensitive to all charge conjugation conserving terms of the expansion, including a $gX^2Y$ term. The statistical uncertainty of all parameters is improved by a factor two with respect to earlier measurements.


I. INTRODUCTION
The isospin violating η → π + π − π 0 decay can proceed via electromagnetic interactions or via strong interactions due to the difference between the masses of u and d quarks. The electromagnetic part of the decay amplitude is long known to be strongly suppressed [1,2]. The recent calculations performed at next-to-leading order (NLO) of the chiral perturbation theory (ChPT) [3,4] reaffirm that the decay amplitude is dominated by the isospin violating part of the strong interaction.
Defining the quark mass ratio, Q, as the decay width at up to NLO ChPT is proportional to Q −4 [5]. The definition in Eq. (1), neglectingm 2 /m 2 s , gives an ellipse in the m s /m d , m u /m d plane with major semi-axis Q [6]: a determination of Q puts a stringent constraint on the light quark masses.
Using Dashen's theorem [7] to account for the electromagnetic effects, Q can be determined at the lowest order from a combination of kaon and pion masses. With this value of Q = 24.2, the ChPT results for the η → π + π − π 0 decay width at LO, Γ LO = 66 eV, and NLO, Γ N LO = 160 ± 50 eV [8], are far from the experimental value Γ exp = 300 ± 11 eV [9]. The discrepancy could originate from higher order contributions to the decay amplitude or from the corrections to the Q value.
Several theoretical efforts have aimed at a better de-arXiv:1601.06985v1 [hep-ex] 26 Jan 2016 scription of the η → 3π decay amplitude. A full NNLO ChPT calculation has been performed [10]. However, it depends on the values of a large number of the coupling constants of the chiral lagrangian which are not known precisely and therefore is affected by large uncertainties. A calculation in unitarized ChPT using relativistic coupled channels [11] gives better agreement with the data but it is model dependent. With a non-relativistic effective field theory framework, the effects of higher order isospin breaking in the final state interactions have been investigated [12]. The ππ rescattering seems to play an important role in this decay, giving about half of the correction in going from LO to NLO [8]. The rescattering can be accounted for to all orders using dispersive integrals. In this framework, there are three approaches to improve the ChPT predictions [13][14][15]. The reliability of the ChPT calculations could be checked by a comparison with the experimental pion distributions represented by the Dalitz plot. Conversely, precise experimental distributions could be used directly for the dispersive analyses [13][14][15] to determine the Q ratio without relying on the higher order ChPT calculations. For the η → π + π − π 0 Dalitz plot distribution, the normalized variables X and Y are commonly used: T i are kinetic energies of the pions in the η rest frame. The squared amplitude of the decay is parametrized by a polynomial expansion around (X, Y ) = (0, 0): The Dalitz plot distribution can then be fit using this formula to extract the parameters a, b, . . ., usually called the Dalitz plot parameters. Note that coefficients multiplying odd powers of X (c, e, h and l) must be zero assuming charge conjugation invariance.
The experimental values of the Dalitz plot parameters are shown in Tab. I together with the parametrization of theoretical calculations. The most precise previous measurement is from KLOE 2008 analysis which was based on 1.34 · 10 6 events [19]. There is some disagreement among the experiments, specially for the b but also for the a parameter. Looking at the theory, both the b and the f parameters deviate from experiment. The new high statistics measurement presented in this paper can help to clarify the tension among the experimental results, and can be used as a more precise input for the dispersive calculations.

II. THE KLOE DETECTOR
The KLOE detector at the DAΦNE e + e − collider in Frascati consists of a large cylindrical Drift chamber (DC) and an electromagnetic calorimeter (EMC) in a 0.52 T axial magnetic field. The DC [23] is 4 m in diameter and 3.3 m long and is operated with a helium -isobutane gas mixture (90% -10%). Charged particles are reconstructed with a momentum resolution of σ(p ⊥ )/p ⊥ 0.4%.
The EMC [24] consists of alternating layers of lead and scintillating fibers covering 98% of the solid angle. The lead-fiber layers are arranged in ∼ (4.4 × 4.4) cm 2 cells, five in depth, and these are read out at both ends. Hits in cells close in time and space are grouped together in clusters. Cluster energy is obtained from the signal amplitude and has a resolution of σ(E)/E = 5.7%/ E(GeV). Cluster time, t cluster , and position are energy weighted averages, with time resolution σ(t) = (57 ps)/ E(GeV) ⊕ 100 ps. The cluster position along the fibers is obtained from time differences of the signals.
The KLOE trigger [25] uses both EMC and DC information. The trigger conditions are chosen to minimize beam background. In this analysis, events are selected with the calorimeter trigger, requiring two energy deposits with E > 50 MeV for the barrel and E > 150 MeV for the endcaps.
The analysis is performed using data collected at the φ meson peak with the KLOE detector in 2004-2005, and corresponds to an integrated luminosity of ∼ 1.6 fb −1 . Due to DAΦNE crossing angle φ mesons have a small horizontal momentum, p φ of about 13 MeV/c. The η mesons are produced in the radiative decay φ → ηγ φ . The photon from the φ radiative decay, γ φ , has an energy E ∼ 363 MeV. The data sample used for this analysis is independent and about four times larger than the one used in the previous KLOE(08) η → π + π − π 0 Dalitz plot analysis [19].
The reconstructed data are sorted by an event classification procedure which rejects beam and cosmic ray backgrounds and splits the events into separate streams according to their topology [26]. The beam and background conditions are monitored. The corresponding parameters are stored for each run and included in the GEANT3 based Monte Carlo (MC) simulation of the detector. The event generators for the production and decays of the φ-meson include simulation of initial state radiation. The final state radiation is included for the simulation of the signal process. The simulation of e + e − → ωπ 0 process (an important background in this analysis) assumes a cross section of 8 nb. The simulations of the background channels used in this analysis correspond to the integrated luminosity of the experimental data set, while the signal simulation corresponds to ten times larger luminosity.
Selection steps are listed below: • A candidate event has at least three prompt neutral clusters in the EMC. The clusters are required to have energy at least 10 MeV and polar angles 23 • < θ < 157 • , where θ is calculated from the distance of the cluster to the beam crossing point (R cluster ). The time of the prompt clusters should be within the time window for massless particles, |t cluster − R cluster /c| < 5σ(t), while neutral clusters do not have an associated track in the DC.
• At least one of the prompt neutral clusters has energy greater than 250 MeV. The highest energy cluster is assumed to originate from the γ φ photon.
• The two tracks with smallest distance from the beam crossing, and with opposite curvature, are chosen. In the following these tracks are assumed to be due to charged pions. Discrimination against electron contamination from Bhabha scattering is achieved by means of Time Of Flight as discussed in the following.
• P φ , the four-momentum of the φ meson, is determined using the beam-beam energy √ s and the φ transverse momentum measured in Bhabha scattering events for each run.
• The γ φ direction is obtained from the position of the EMC cluster while its energy/momentum is calculated from the two body kinematics of the φ → ηγ φ decay: where θ φ,γ is the angle between the φ and the γ φ momenta. The four-momentum of the η meson is then: • The π 0 four-momentum is calculated from the missing four-momentum to η and the charged pions: • To reduce the Bhabha scattering background, the following two cuts are applied: a cut in the (θ +γ ,θ −γ ) plane as shown in Fig. 1, where θ +γ (θ −γ ) is the angle between the π + (π − ) and the closest photon from π 0 decay. a cut in the (∆t e , ∆t π ) plane as shown in Fig. 2, to discriminate electrons from pions, where ∆t e , ∆t π are calculated for tracks which have an associated cluster, ∆t e/π ≡ t track e/π − t cluster , where t track e/π , is the expected arrival time to EMC for e/π with the measured momentum, and t cluster the measured time of the EMC cluster.
• To improve the agreement between simulation and data, a correction for the relative yields of: (i) e + e − → ωπ 0 , and (ii) sum of all other backgrounds, with respect to the signal is applied. The correction factors are obtained from a fit to the distribution of the azimuthal angle between the π 0 decay photons, in the π 0 rest frame, θ * γγ (Fig. 3). The uncertainties of the correction factors are evaluated by comparing the corresponding fit to the distribution of the missing mass squared, P 2 π 0 (Fig. 4).
• To further reduce the background contamination, two more cuts are applied: The overall signal efficiency is 37.6% at the end of the analysis chain and the signal to background ratio is 133.
As can be seen in Figs. 3, 4 and 5 the agreement of simulation with the experimental data is good.

IV. DALITZ PLOT
For the Dalitz plot, a two dimensional histogram representation is used. The bin width is determined both by the resolution in the X and Y variables and the number of events in each bin, which should be large enough to justify χ 2 fitting. The resolution of the X and Y variables is evaluated with MC signal simulation. The distribution of the difference between the true and reconstructed values is fit with a double Gaussian. The standard deviations of the narrower Gaussians are δ X = 0.021 and δ Y = 0.032. The range (−1, 1) for the X and Y variables was divided into 31 and 20 bins, respectively. Therefore the bin widths correspond to approximately three standard deviations. The minimum bin content is 3.3 · 10 3 events. Fig. 6 shows the distributions of the θ * γγ and the P 2 π 0 variables for two bins in the Dalitz plot, one with the largest content and one with the smallest. As can be seen, the signal and the background are well reproduced by the simulation. Fig. 7 shows the experimental Dalitz plot distribution after background subtraction, which is fit to the amplitude expansion from Eq. (5) to extract the Dalitz plot parameters. Only n = 371 bins which are fully inside the kinematic boundaries are used and there are ∼ 4.7 · 10 6 entries in the background subtracted Dalitz plot.
The fit is performed by minimizing the χ 2 like function where: given by Eq. (5). The integral is over X and Y in the allowed phase space for bin j. The sum over j bins includes all Dalitz plot bins at least partly inside the physical border, n T .
is the background subtracted content of Dalitz plot bin i, where β 1,2 are the scaling factors, B i1 is the ωπ 0 background in the bin i and B i2 is the same for the remaining background.
• S ij is the acceptance and smearing matrix from bin j to bin i in the Dalitz plot. It is determined from signal MC by S ij = N rec,i;gen,j /N gen,j , where N rec,i;gen,j denotes the number of events reconstructed in bin i which were generated in bin j and N gen,j denotes the total number of events generated in bin j.
The input-output test of the fit procedure was performed using signal MC generated with the same statistics as the experimental data. The extracted values for the parameters were within one standard deviation with respect to the input.
The fit has been performed using different choices of the free parameters in Eq. (5), with the normalization N and the parameters a, b and d always let free. The main fit results are summarized in Tab. II. The first row (set #1) includes all parameters of the cubic expansion, Eq. (5). The fit values of the charge conjugation violating parameters c, e, h and l are consistent with zero (c = (4.3 ± 3.4) · 10 −3 , e = (2.5 ± 3.2) · 10 −3 , h = (1.1 ± 0.9) · 10 −2 , l = (1.1 ± 6.5) · 10 −3 ) and are omitted from the table. Therefore in all remaining fits the charge conjugation violating parameters c, e, h and l are set to zero. The result #2 with f = g = 0 demonstrates that it is not possible to describe the distribution with only quadratic terms. The fits including in addition the f parameter improves the accuracy largely providing a good description. A complementary test with f fixed to zero and free g parameter leads to much worse χ 2 value. This observation is consistent with the recent experimental results that used the same free parameters, in particular the KLOE(08) analysis. In the set #4 which includes both f and g parameters in the fit, the g parameter differs from zero at the 4.9σ level and the fit probability is as good as for the fit #1. For completeness in the further discussions we include also set #3 with g = 0, since it enables direct comparison to the previous experiments (KLOE(08), WASA (14) and BESIII (15)). The correlation matrices for fits #3 and #4 are: The fit #4 is compared to the background subtracted Dalitz plot data, N i , in Fig. 8. The red lines represent the fit result and correspond to separate slices in the Y variable. Fig. 9 shows the distribution of the normalized residuals for the fit #4: r i = (N i − n j=1 S ij N T,j )/σ i . The location of the residuals r i > 1 and r i < −1 on

V. ASYMMETRIES
While the extracted Dalitz plot parameters are consistent with charge conjugation symmetry, the unbinned integrated charge asymmetries provide a more sensitive test. The left-right (A LR ), quadrant (A Q ) and sextant (A S ) asymmetries are defined in Ref. [27]. The same background subtraction is applied as for the Dalitz plot parameter analysis. For each region in the Dalitz plot used in the calculation of the asymmetries, the acceptance is calculated from the signal MC as the ratio between the number of the reconstructed and the generated events. The yields are then corrected for the corresponding efficiency. The procedure was tested using signal MC generated with the same statistics as the experimental data.

VI. SYSTEMATIC CHECKS
To quantify and account for systematic effects in the results, several checks have been made.
• Minimum photon energy cut (EGmin) is changed from 10 MeV to 20 MeV (for comparison the EMC energy resolution varies from 60% to 40% for this energy range). The systematic error is taken as half of the difference.
• Background subtraction (BkgSub) is checked by determining the background scaling factors for each bin (or region for the asymmetries) of the Dalitz plot separately. With the same method as for the whole data sample, using the θ * γγ and P 2 π 0 distributions, background scaling factors are determined for each bin (or region). The systematic error is taken as half the difference with the standard result.
• Choice of binning (BIN) is tested by varying number of bins of the Dalitz plot. For X and Y simultaneously, the bin width is varied from ∼ 2δ X,Y to ∼ 5δ X,Y , in total 10 configurations. The systematic uncertainty is given by the standard deviation of the results.
• ∆t e , ∆t π cut: the offsets of the horizontal and diagonal lines shown in Fig. 2 were varied by ±0.22 ns and ±0.21 ns, respectively.
• Missing mass cut (MM) is tested by varying the cut by ±2.0 MeV, ∼ 1σ. For this cut a stronger dependence of the parameters on the cut was noted. This has been further investigated by performing the Dalitz plot parameter fit for one parameter at a time, for each step, and keeping the other parameters fixed at the value for the standard result. Since the dependence was reduced when varying just one parameter, we conclude that it is mostly due to the correlations between parameters.  • Event classification procedure (ECL) is investigated by using a prescaled data sample without the event classification bias (collected with prescaling factor 1/20). The fraction of events remaining in each Dalitz plot bin after the event classification conditions varies between 94% and 80% for different bins and it is very well described by the MC within the errors. The analysis of the prescaled data follows the standard chain. The systematic error is extracted as half the difference between the results of the analysis with and without the event classification procedure.
Unless stated otherwise the systematic error is calculated as the difference between the two tests and the standard result. If both differences have the same sign, the asymmetric error is taken with one boundary set at zero and the other at the largest of the differences. The resulting systematic error contributions for the Dalitz plot parameters for the sets #4 and #3 are summarized in Tab. III  and Tab. IV, respectively. The results for the charge asymmetries are summarized in Tab. V.    These results confirm the tension with the theoretical calculations on the b parameter, and also the need for the f parameter. In comparison to the previous measurements shown in Tab. I, the present results are the most precise and the first including the g parameter. The improvement over KLOE(08) analysis comes from four times larger statistics and improvement in the systematic uncertainties which are in some cases reduced by factor 2 − 3. The major improvement in the systematic uncertainties comes from the analysis of the effect of the Event classification with an unbiased prescaled data sample.
The final values of the charge asymmetries are all con- The systematic and statistical uncertainties are of the same size except for the A LR which is dominated by the systematic uncertainty due to the description of the Bhabha background.