Study of the process e + e − → π + π − π 0 using initial state radiation with BABAR

The process e þ e − → π þ π − π 0 γ is studied at a center-of-mass energy near the ϒ ð 4 S Þ resonance using a data sample of 469 fb − 1 collected with the BABAR detector at the PEP-II collider. We have performed a precise measurement of the e þ e − → π þ π − π 0 cross section in the center-of-mass energy range from 0.62 to 3.5 GeV. In the energy regions of the ω and ϕ resonances, the cross section is measured with a systematic uncertainty of 1.3%. The leading-order hadronic contribution to the muon magnetic anomaly calculated using the measured e þ e − → π þ π − π 0 cross section from threshold to 2.0 GeV is ð 45 . 86 (cid:2) 0 . 14 (cid:2) 0 . 58 Þ × 10 − 10 . From the fit to the measured 3 π mass spectrum we have determined the resonance parameters Γ


I. INTRODUCTION
The process e þ e − → π þ π − π 0 1 has the second largest hadronic cross section after e þ e − → π þ π − in the energy region below 1 GeV and is therefore very important for the Standard Model calculation of the anomalous magnetic moment of the muon a μ ≡ ðg − 2Þ μ =2. Currently, the accuracy of the e þ e − → π þ π − π 0 contribution to the muon magnetic anomaly (a 3π μ ) is about 3% [1] and needs to be improved.
The most precise measurements of the e þ e − → π þ π − π 0 cross section in the energy region of the ω and ϕ resonances were performed by the SND and CMD-2 Collaborations at the VEPP-2M e þ e − collider [2][3][4][5]. Above the ϕ meson resonance the latest measurements come from the BABAR experiment [6], which used the initial-state radiation (ISR) technique, and the SND experiment at the VEPP-2000 e þ e − collider [7]. There is also a preliminary result from the BESIII experiment [8], which measured the e þ e − → π þ π − π 0 cross section in the energy range between 0.7 and 3.0 GeV using the ISR technique.
One of the reasons for the relatively low accuracy of a 3π μ is the difference between the cross section measurements in different experiments. For example, the SND cross section near the ω [3] is about 8% (1.8σ) larger than the cross section measured by CMD-2 [4]. BABAR did not measure the cross section in this region, but fitted to the 3π mass spectrum in the e þ e − → π þ π − π 0 γ reaction with the vectormeson-dominance (VMD) model [6] and determined the ω parameters. The BABAR value for the ω peak cross section as well as the BESIII preliminary result [8] support a larger cross section value, as obtained by SND [3]. It is generally accepted that the process e þ e − → π þ π − π 0 proceeds mainly through the ρð770Þπ (ρ þ π − þρ − π þ þρ 0 π 0 ) intermediate state. This assumption has been well tested at * Deceased. the ω and ϕ resonances [9,10]. The dynamics of e þ e − → π þ π − π 0 in the energy range between 1.1 and 2 GeV were recently studied in Ref. [7]. This study confirms the dominance of the ρð770Þπ channel below 1.5 GeV. However, in this region there is a 10-20% contribution from the isovector ωπ 0 mechanism and its interference with the dominant ρð770Þπ amplitude. In the region of the ωð1650Þ resonance (1.55-1.75 GeV), a large contribution of the ρð1450Þπ intermediate state was observed, which is comparable with that of the ρð770Þπ. A relatively large fraction (∼10%) of the ρð1450Þπ channel was also observed in the J=ψ → 3π decay [11].
In this article we update the BABAR e þ e − → π þ π − π 0 measurement [6] using a data set that is 5 times larger. We study the process e þ e − → π þ π − π 0 γ, where the photon emission is caused by initial-state radiation. The Born cross section for this process integrated over the momenta of the hadrons is given by where ffiffi ffi s p is the e þ e − center-of-mass (c.m.) energy, x ≡ 2E γ = ffiffi ffi s p , E γ , and θ γ are the photon energy and polar angle in the c.m. frame, and σ 0 is the Born cross section for e þ e − → π þ π − π 0 . The so-called radiator function (see, for example, Ref. [12]) describes the probability of ISR photon emission for θ γ ≫ m e = ffiffi ffi s p . Here, α is the fine structure constant and m e is the electron mass. The ISR photons are emitted predominantly at small angles relative to the initial electron or positron directions; however, about 10% of the photons have c.m. polar angles in the range 30°< θ γ < 150°. In the present analysis, we require that the ISR photon be detected.
The goal of this analysis is to improve the accuracy of the e þ e − → π þ π − π 0 cross section measurement and the contribution of this process to a μ .

II. THE BABAR DETECTOR AND DATA SAMPLES
In this article a data sample of 469 fb −1 , collected with the BABAR detector [13] at the PEP-II asymmetric-energy storage ring at the SLAC National Accelerator Laboratory, is analyzed. At PEP-II, 9 GeV electrons collide with 3.1 GeV positrons at a center-of-mass energy of 10.58 GeV [ϒð4SÞ resonance]. About 91% of the integrated luminosity was recorded at 10.58 GeV, while 9% was recorded at 10.54 GeV.
Charged-particle tracking for the BABAR detector is provided by a five-layer silicon vertex tracker (SVT) and a 40-layer drift chamber (DCH), operating in a 1.5 T axial magnetic field. The transverse momentum resolution is 0.47% at 1 GeV=c. Energies of photons and electrons are measured with a CsI(Tl) electromagnetic calorimeter (EMC) with a resolution of 3% at 1 GeV. Charged-particle identification is provided by measurements of ionization losses, dE=dx, in the SVT and DCH, and by an internally reflecting ring-imaging Cherenkov detector. Muons are identified in the solenoid's instrumented flux return.
Signal and background ISR processes are simulated by a Monte Carlo (MC) event generator based on the approach suggested in Ref. [14]. A model of the ρð770Þπ intermediate state is used to simulate the signal process e þ e − → 3πγ. The extra-photon radiation from the initial state is implemented with the structure function technique [15], while the final-state radiation (FSR) is simulated using the PHOTOS package [16]. Since the ISR photon is emitted predominantly at small angles relative to the beam directions, the events are generated with the restriction 20°< θ γ < 160°, where θ γ is the ISR photon polar angle in the c.m. frame. We also require that the invariant mass of the hadron system and ISR photon together be greater than 8 GeV=c 2 . This condition restricts the maximum energy of extra photons emitted by the initial particles.
The following background ISR processes are simulated: e þ e − → π þ π − γ, μ þ μ − γ, K þ K − γ, K S K L γ, K þ K − π 0 γ, K S K − π þ γ, π þ π − π 0 π 0 γ, π þ π − ηγ, ωηγ, and ωπ 0 π 0 γ. The backgrounds from non-ISR hadronic processes e þ e − → qq, where q ¼ u, d, s, and from e þ e − → τ þ τ − are simulated with the JETSET [17] and KK2f [18] packages, respectively. The interaction of the generated particles with the BABAR detector and the detector response are simulated using the GEANT4 [19] package. The simulation takes into account the variation of the detector and accelerator conditions, and in particular describes the beam-induced background, which leads to the appearance of spurious photons and tracks in the events of interest.

III. EVENT SELECTION
The selection of e þ e − → π þ π − π 0 γ candidates is based on the requirement that all the final particles be detected and well reconstructed. We select events with exactly two good quality opposite-sign charged tracks, which are considered as π þ and π − candidates, and at least three photons. The "good" tracks are required to have a transverse momentum above 100 MeV=c, originate from the interaction region, and to be not identified as an electron. Their laboratory polar angle must be between 23°and 140°. An event can contain any number of extra tracks not satisfying the above criteria.
The photons must have energies above 100 MeV and be in the well-understood region of the calorimeter 23°< θ < 137.5°. One of the photons (the ISR candidate photon) is required to have a c.m. energy larger than 3 GeV.
The remaining photons must form at least one π 0 candidate, a pair of photons with invariant mass in the range 0.1-0.17 GeV=c 2 .
For events satisfying the selection criteria described above, a kinematic fit is performed with requirements of energy and momentum conservation, and the π 0 mass constraint for the candidate π 0 . The MC simulation does not accurately reproduce the shape of the resolution function for the photon energy. To reduce the effect of the data-MC simulation difference in the energy resolution, the fit uses only the measured direction for the ISR photon candidate; its energy is a free fit parameter. For events with two or more candidate π 0 s, all possible π þ π − π 0 γ combinations are tested and the one with the minimum χ 2 of the kinematic fit (χ 2 3πγ ) is used. As a result of the kinematic fit we obtain the corrected three-pion invariant mass (M 3π ).
The χ 2 3πγ distribution for events from the 3π mass region near the ω mass (0.67-0.87 GeV=c 2 ), where the contribution of background processes is small, is shown in Fig. 1. In further analysis we use two conditions on this parameter: the standard χ 2 3πγ < 40 and the tighter χ 2 3πγ < 20. The latter is applied for the e þ e − → π þ π − π 0 cross section measurement. The χ 2 3πγ distribution for the mass range 1.05 < M 3π < 3.00 GeV=c 2 is shown in Fig. 2. In this region the background contribution is significantly larger.
Events with charged kaons in the final state (e þ e − → K þ K − π 0 γ, e þ e − → K þ K − γ, etc.) are suppressed by the requirement that none of the "good" charged tracks be identified as a kaon.
Two-body ISR events from the processes e þ e − → π þ π − γ and e þ e − → μ þ μ − γ with extra spurious photons are suppressed by the two conditions: E π 0 > 0.4 GeV, where E π 0 is the energy of the candidate π 0 , and M 2 rec > 5 GeV 2 =c 4 , where M rec is the mass recoiling against the π þ π − pair. Some fraction of background e þ e − → π þ π − π 0 π 0 γ events contain additional π 0 candidates. For these events we perform a kinematic fit to the 4πγ hypothesis and apply the condition χ 2 4πγ > 30, which reduces the 4πγ background by a factor of 2.
Another important background source is e þ e − → qq events containing a very energetic π 0 in the final state. A fraction of these events is seen as a peak at the π 0 mass in the M Ã γγ distribution, where M Ã γγ is the invariant mass of two photons, one of which is the most energetic in an event. The second photon is required to have an energy above 100 MeV. Once all possible photon pair combinations are checked, the one with closest invariant mass to the π 0 mass is chosen. Events with 0.10 < M Ã γγ < 0.17 GeV=c 2 are rejected.
The e þ e − → qq background is dominated by e þ e − → π þ π − π 0 π 0 events. Events of this process passing the π þ π − π 0 γ selection criteria have a χ 2 3πγ distribution peaked at low values, similar to the signal. A fraction of these events proceeding via ρ þ ρ − intermediate state is rejected by the condition M πγ > 1.5 GeV=c 2 , where M πγ is the invariant mass of the most energetic photon and one of the charged  pions. This condition also rejects e þ e − → τ þ τ − events, which imitate π þ π − π 0 γ events when both τ's decay into ρν.
In the 3π mass region below 1.1 GeV=c 2 , which is the most important for the calculation of a 3π μ , the background suppression requirements decrease the fraction of background events from 5% to 2%, with loss of signal events of 15%. The χ 2 3πγ distribution for data events rejected by the background suppression requirements in the mass region 1.05 < M 3π < 3 GeV=c 2 is shown as the shaded histogram in Fig. 2. In this region, the background is suppressed by a factor of 2.6 with a loss of signal events of 17%.

IV. BACKGROUND ESTIMATION AND SUBTRACTION
To estimate background, the samples of simulated events listed in Sec. II are normalized to the collected integrated luminosity. Before normalization, the hadron mass spectrum for a particular simulated process is reweighted using Eq. (1) and the existing data on its Born cross section. For the most important background ISR processes e þ e − → K þ K − π 0 γ, e þ e − → π þ π − γ, and e þ e − → π þ π − π 0 π 0 γ, data samples selected with special criteria are used to determine additional scale factors.
The mass distribution for events with two charged kaons surviving our selection (dN 0K =dm) is obtained from the distribution of events with two identified kaons: simulation corrected for datasimulation differences in the charged-kaon identification efficiency. The observed spectrum of two-kaon background events is almost completely saturated by the e þ e − → K þ K − π 0 γ process.
The scale factor for the e þ e − → π þ π − γ process is estimated using events with 40 < χ 2 3πγ < 250 and M 2 rec < 10 GeV 2 =c 4 . The latter condition suppresses contributions of all processes except e þ e − → π þ π − γ. The scale factor is found to be 1.6 AE 0.2. The quoted systematic uncertainty is estimated by variation of the conditions on χ 2 3πγ and M 2 rec . The large difference between the fitted and expected numbers of e þ e − → π þ π − γ background events may be the result of an inaccurate simulation of the nuclear interactions of charged π mesons in the calorimeter. In particular, the number of fake photons due to nuclear interactions may be different in data and simulation.
The process e þ e − → π þ π − π 0 π 0 γ is the main source of background for the process under study. Several intermediate states (ωπ 0 , a 1 π, ρ þ ρ − , etc.) contribute to this process. Our MC event generator incorrectly reproduces both the 4π mass spectrum for e þ e − → π þ π − π 0 π 0 γ events and the relation between intermediate states, in particular, the fraction of ωπ 0 events. Therefore, the normalization for this process is performed in two stages. In the first stage, we select events with two charged particles and at least five photons, perform a kinematic fit to the e þ e − → π þ π − π 0 π 0 γ hypothesis, and select events with χ 2 4πγ < 30. We measure the 4π mass spectrum and reweight the e þ e − → π þ π − π 0 π 0 γ simulation using the ratio of the data and simulated spectra as a weight function. The reweighting is performed separately for ωπ 0 and non-ωπ 0 events. In the second stage, we analyze the 3π mass spectrum below 1.1 GeV=c 2 for events with 50 < χ 2 3πγ < 500 and M 2 rec > 10 GeV 2 =c 4 . The latter condition is applied to suppress the e þ e − → π þ π − γ background. The spectrum shown in Fig. 3 is fitted with a sum of simulated signal and background distributions. The fitted parameters are scale factors for the e þ e − → π þ π − π 0 γ and e þ e − → π þ π − π 0 π 0 γ distributions. The difference in the line shape of the ω peak between data and simulation seen in Fig. 3 is attributed to inaccurate simulation of the tails of the M 3π resolution function at large χ 2 3πγ values. The π þ π − π 0 π 0 γ scale factor is found to be 1.30 AE 0.15. The quoted uncertainty is systematic. It is estimated by variation of the conditions on χ 2 3πγ and M 2 rec . The total contribution to the background from other ISR processes at M 3π < 1.1 GeV is calculated to be less than 1=50 of the e þ e − → π þ π − π 0 π 0 γ background.
The e þ e − → qq background events can be divided into two classes. The first class (4π) contains events from the e þ e − → π þ π − π 0 π 0 process. The second (non-4π) contains events from all other processes. The 4π events has a χ 2 3πγ distribution peaked at small values similar to the signal The 3π invariant mass spectrum for data events with 50 < χ 2 3πγ < 500 and M 2 rec > 10 GeV 2 =c 4 . The solid histogram represents the result of the fit with a sum of signal and background distributions. The light-shaded (green) area represents the fitted e þ e − → π þ π − π 0 π 0 γ contribution, while the dark-shaded (blue) histograms are the spectrum for all other background processes.
process e þ e − → π þ π − π 0 γ. For the second class, the χ 2 3πγ distribution has a wide maximum near χ 2 3πγ ¼ 300. The ISR photon in 4π and most non-4πqq events is imitated by a photon from the π 0 decay. Therefore, to estimate these backgrounds we study the M Ã γγ distribution. The M Ã γγ distribution for data events with 0.6 < M 3π < 3.5 GeV=c 2 selected using our standard selection criteria except for the condition on M Ã γγ is shown in Fig. 4. The events in the π 0 peak originate mainly from the 4π class, while the nearly flat distribution is dominated by e þ e − → π þ π − π 0 γ events. The distribution is fitted with the sum of a Gaussian and a linear function. The non-4π background is also estimated from the number of events in the π 0 peak in the M Ã γγ distribution, but for events with 40 < χ 2 3πγ < 200. The 3π mass region 0.6-3.5 GeV=c 2 is divided into 29 intervals with 0.1 GeV=c 2 width. For each M 3π interval, we determine the numbers of 4π and non-4π events in data and e þ e − → qq simulation from the fit to the M Ã γγ distribution. The obtained data spectrum for 4π events is compared with the same spectrum for simulated events in Fig. 5. It is seen that the e þ e − → qq simulation reproduces reasonably well the total number of selected e þ e − → π þ π − π 0 π 0 events. The overall scale factor for the simulation is 0.83 AE 0.05. However the shapes of the M 3π spectra for data and simulation are different, especially in the region 1.3-1.8 GeV=c 2 . At M 3π > 0.9 GeV=c 2 the ratio of the data and simulated spectra shown in Fig. 5 is used to reweight the yield of simulated e þ e − → π þ π − π 0 π 0 events. It should be noted that the ratio of the number of 4π events selected with our standard criteria to the number of events shown in Fig. 5 is about five. The uncertainty in the number of 4π background events obtained using the reweighted simulation is dominated by the uncertainty in the number of events in each mass bin in Fig. 5.
An excess of data over simulation is seen in Fig. 5 in the mass region 0.7-0.9 GeV=c 2 . This excess may be an indication of a contribution from the e þ e − → ωπ 0 → π þ π − π 0 π 0 process, which is absent in our e þ e − → qq simulation. This process produces events peaked at the ω mass. We repeat the fitting procedure described above with finer binning. The result is shown in Fig. 6. This spectrum is used to estimate the e þ e − → π þ π − π 0 π 0 background in the π þ π − π 0 mass region from 0.7 to 0.9 GeV=c 2 . To do this, the data spectrum in Fig. 6 is multiplied by a scale factor of five, obtained in the region M 3π > 0.9 GeV. The systematic uncertainty in this estimation is taken to be 100%. The same scale factor is used for the interval 0.6 < M 3π < 0.7 GeV=c 2 , where the number of fitted 4π events in Fig. 5 is 0.7 AE 1.

2.
A similar procedure is used to reweight the non-4πqq simulation. For this background we also analyze events with M Ã γγ near the η meson mass. The M 3π spectra for the 4π and non-4πqq background events selected with the standard criteria are shown in Fig. 7. It is seen that the fraction of non-4πqq events is relatively small.
The mass region 0.6 < M 3π < 1. standard criteria is shown in Fig. 8. The points with error bars in Fig. 8 represent the estimated total background contribution from the sources described above. The background M 3π spectrum on a linear scale is displayed in Fig. 9 (left). The filled histogram represents the contribution of all background sources except 2πγ and 4πγ. About two-thirds of events in this histogram come from the e þ e − → π þ π − π 0 π 0 process. The open histogram is a sum of the filled histogram and the e þ e − → π þ π − γ background spectrum. It is seen from Fig. 9 (left) that the background in this M 3π region is dominated by the processes e þ e − → π þ π − π 0 π 0 γ and e þ e − → π þ π − γ. The ratio of the background spectrum to the data spectrum is shown in Fig. 9 (right). The background fraction decreases from ð25 AE 15Þ% at 0.65 GeV=c 2 to ð7 AE 3Þ% at 0.7 GeV=c 2 and to ð0.5 AE 0.1Þ% in the ω region, then increases to ð9 AE 2Þ% at 0.9 GeV=c 2 and decreases again to ð0.5 AE 0.1Þ% at the ϕ. Near 1.05 GeV=c 2 , where the e þ e − → π þ π − π 0 cross section has a minimum, the background fraction is ð27 AE 5Þ%. With the tighter selection χ 2 3πγ < 20 the background fraction decreases by a factor of about two.
In the region 0.6 < M 3π < 1.1 GeV=c 2 the estimated background is subtracted from the number of selected data events in each M 3π bin. It should be noted that the numbers of background events in different mass bins are correlated. This correlation arises from the uncertainties in the scale factors for e þ e − → π þ π − π 0 π 0 γ and e þ e − → π þ π − γ events, which are equal to 10.5% and 12.5%, respectively.
The M 3π spectrum for selected data events with 1.1 < M 3π < 3.5 GeV=c 2 is shown in Fig. 10 (left). The points with error bars in Fig. 10 (left) represent the calculated spectrum for background events, while the filled histogram shows the background spectrum with the e þ e − → π þ π − π 0 π 0 contribution subtracted. It is seen that the process e þ e − → π þ π − π 0 π 0 becomes the dominant background source above 1.5 GeV=c 2 . This background has a χ 2 3πγ distribution similar to that for signal events. It is estimated as described above and subtracted from the data M 3π spectrum. The e þ e − → K þ K − π 0 γ background is also estimated from data. It is found to be relatively small, about 4% of the e þ e − → π þ π − π 0 π 0 contribution. Figure 10 (right) displays the calculated background from all other sources. Here, the dominant contribution arises from the e þ e − → π þ π − π 0 π 0 γ process. The next largest contribution comes from non-4πqq events.  The mass region 1.1 < M 3π < 3.5 GeV=c 2 is divided into 72 bins. The bin width is 25 MeV=c 2 below 2.7 GeV=c 2 and 100 MeV=c 2 above. In this region, the background from ISR processes, even from e þ e − → π þ π − π 0 π 0 γ, cannot be estimated with the same precision as at low masses, because the MC event generator does not include many intermediate states contributing to the ISR processes. Therefore, a procedure of background subtraction based on the difference in χ 2 3πγ distributions for signal and background events is used. In each mass bin, we subtract background events of the The points with error bars show the calculated spectrum for background events. The filled histogram represents the background spectrum with the e þ e − → π þ π − π 0 π 0 contribution subtracted. Right panel: The M 3π spectrum for background events from all sources except e þ e − → π þ π − π 0 π 0 and e þ e − → K þ K − π 0 γ (points with error bars). The histogram represents the same spectrum with the e þ e − → π þ π − π 0 π 0 γ contribution subtracted.
e þ e − → K þ K − π 0 γ and e þ e − → π þ π − π 0 π 0 processes and determine the numbers of events with χ 2 3πγ ≤ 20 (N 1 ) and 20 < χ 2 3πγ < 40 (N 2 ). The numbers of signal (N sig ) and remaining background (N bkg ) events are then determined from the system of linear equations: The coefficients α ¼ N 1 =ðN 1 þ N 2 Þ for pure signal and background events are determined from simulation. The mass dependence of the coefficient α sig is shown in Fig. 11. The values of α sig at the ϕ and J=ψ masses can be extracted from data. In the ϕ mass region, we determine N 1 and N 2 for pure signal events by subtracting the calculated background. In the J=ψ mass region, the same numbers are obtained using a fit to the M 3π spectrum by a sum of a J=ψ line shape and a linear function (see Sec. XI). The resulting values of α sig are 0.859 AE 0.003 at the ϕ mass and 0.890 AE 0.005 at the J=ψ mass. Their ratios to the corresponding values obtained from simulation are R ϕ ¼ 1.004 AE 0.004 and R J=ψ ¼ 1.018 AE 0.007, respectively. In Eqs. (3), we use for α sig the fitting function shown in Fig. 11 multiplied by a linear function interpolating between R ϕ and R J=ψ .
The α bkg coefficient is determined using a mixture of background simulated events shown in Fig. 10 (right). The coefficient is practically independent of mass and equal to 0.316 AE 0.007. To estimate the systematic uncertainty in α bkg , we vary the fraction of non-π þ π − π 0 π 0 γ events in the mixture of simulated background events by a factor of two. The variation in the α bkg value is taken as a measure of the systematic uncertainty. It is less than 5% below 2 GeV=c 2 , 8% between 2 and 3 GeV=c 2 , and 15% above 3 GeV=c 2 .
The M 3π spectrum for background events obtained by the solution of the system of equations (3) is shown in Fig. 12 in comparison with the same spectrum obtained using simulation. It is seen that the simulation reproduces the data spectrum reasonably well up to 3 GeV. The M 3π spectrum for signal events is shown in Fig. 13.

V. FINAL-STATE RADIATION
A high-energy photon can be also emitted from the final state. Since the 3π system in the ISR and final-state radiation processes has different C-parity, the contribution of the interference between them to the total cross section vanishes when integrating over the final hadron momenta. We analyze two FSR mechanisms. The first is emission of the photon by charged pions. Its cross section is calculated as σ 3π ð10.58 GeVÞf FSR , where f FSR is the FSR probability. The e þ e − → π þ π − π 0 cross section at 10.58 GeV can be estimated from the CLEO measurement at 3.67 GeV σ 3π ð3.67 1Þ pb [20]. Perturbative QCD (pQCD) predicts the same asymptotic energy dependence 1=E 8 , where E is the c.m. energy, for all vectorpseudoscalar (e þ e − → VP) cross sections [21,22]. This prediction can be tested experimentally using the CLEO [20] and Belle [23,24] results for e þ e − → VP cross sections at 3.67 GeV and 10.58 GeV, respectively. For the most accurately measured cross sections for e þ e − → ρη, ωπ 0 , and K Ã K, the ratio σð3.67 GeVÞ=σð10.58 GeVÞ ≈ 3000, which corresponds to the dependence 1=E 7.6 . With this dependence, σ 3π ð10.58 GeVÞ is expected to be about 4.4 fb.
The mass region under study M 3π < 3.5 GeV=c 2 corresponds to the FSR photon c.m. energy E Ã γ > 4.7 GeV. Such a photon can be radiated only by the most energetic pion in the e þ e − → π þ π − π 0 process. For the dominant mechanism e þ e − → ρð770Þπ → π þ π − π 0 , the c.m. energy of the most energetic pion is 5.26 GeV. To estimate the FSR probability, we use the formula for the FSR e þ e − → π þ π − γ cross section from Ref. [25] obtained for pointlike pions. The FSR probability for the π þ π − final state at 10.52 GeV [f 2π ðE Ã γ > 4.7 GeVÞ ¼ 0.26α=π] must be multiplied by a factor of 1=3 (only the most energetic pion in the 3π final state can emit such a photon and this pion must be charged). Thus, the FSR contribution to the e þ e − → π þ π − π 0 γ cross section under the assumption that the photon is emitted by charged pions is estimated to be about 0.001 fb and is negligible.
The second FSR mechanism is photon emission from the quarks, which then hadronize into π þ π − π 0 . In the 3π mass region under study, this process is expected to be dominated by production of C ¼ þ1 resonances decaying to π þ π − π 0 , e.g., the processes e þ e − → ηγ, a 1 ð1260Þγ, The process e þ e − → ηγ has a 3π invariant mass well below the mass range under study. This process and the process e þ e − → η 0 γ were studied by BABAR in Ref. [26]. The measured e þ e − → ηγ and e þ e − → η 0 γ cross sections are 4.5 þ1.2 −1.1 AE 0.3 fb and 5.4 AE 0.8 AE 0.3 fb, respectively. In Ref. [26] they are compared with the pQCD prediction obtained with asymptotic η and η 0 distribution amplitudes, 2.2 fb and 5.5 fb, respectively.
The cross section for the processes e þ e − → a 1 ð1260Þγ, a 2 ð1320Þγ at large c.m. energy is given by [27] dσðe þ e − → MγÞ d cos where F Mγγ is a meson-photon transition form factor for the helicity-zero state, which dominates at large momentum transfers, and where I M is an integral depending on the shape of the meson distribution amplitude. For the asymptotic distribution amplitude, I a 1 ¼ 6 and I a 2 ¼ 10. With the meson decay constants, f a 1 ≈ 200 MeV [28], and f a 2 ≈ f f 2 ≈ 110 MeV [29], the cross sections for the processes e þ e − → a 1 ð1260Þγ and a 2 ð1320Þγ are estimated to be 6.4 fb and 5.4 fb, respectively. There are no experimental data for these cross sections. There is, however, a measurement of the e þ e − → f 2 ð1260Þγ cross section at 10.58 GeV performed by BABAR [30]: ð37 þ24 −18 Þ fb, which is in reasonable agreement with the prediction σ f 2 γ ≈ ð25=9Þσ a 2 γ ≈ 15 fb [27]. The radiative process with an excited pion e þ e − → πð1300Þγ is expected to be small because of the suppression of the πð1300Þ leptonic decay constant [31].
To estimate the detection efficiency for the FSR processes we assume that the efficiency is weakly dependent on the internal structure of the 3π state and reweight simulated ISR e þ e − → 3πγ events to reproduce the photon angular distribution given by Eq. (4). The obtained detection efficiency at the a 2 ð1320Þ mass is 17.9% for the standard selection criteria. The 3π mass distribution for the e þ e − → Mγ process has a resonance shape. The expected mass spectrum for the FSR processes, calculated as a sum of the a 1 ð1260Þ, a 2 ð1320Þ, a 1 ð1640Þ, and a 2 ð1700Þ Breit-Wigner functions, is shown in Fig. 14 by the solid histogram. Interference between amplitudes of different resonances may strongly modify this spectrum. The effect of interference is demonstrated in Fig. 14. We take into account the interference between the a 1 ð1260Þ and a 1 ð1640Þ amplitudes, and a 2 ð1320Þ and a 2 ð1700Þ amplitudes, but neglect the interference between a 1 and a 2 states. The dotted (dashed) histogram represents the result for relative phases between resonances equal to 0 (π). We subtract the spectrum without interference from the spectra for the selected data events shown in Figs. 8 and 13. The systematic uncertainty in the FSR contribution, which takes into account the uncertainty in the theoretical prediction and the effect of interference, is estimated to be 100%. The fraction of the FSR background is maximal (7-8%) in the region 1.05-1.08 GeV=c 2 , where the measured M 3π spectrum has a minimum, and near M 3π ¼ 1.32 GeV=c 2 . Near 1.7 GeV=c 2 , the background fraction is about 6%.
In the mass region near 2 GeV=c 2 , there are several poorly established excited a 1 and a 2 states [35]. We model their contribution by a sum of the a 1 ð1930Þ and a 2 ð2030Þ resonances assuming that f a 1 ð1930Þ Bða 1 ð1930Þ → π þ π − π 0 Þ ≈ 0.2f 2 and f a 2 ð2030Þ Bða 2 ð2030Þ → π þ π − π 0 Þ ≈ 0.2f 2 a 2 ð1320Þ . The latter relation is based on the results of the measurement of the γγ → π þ π − π 0 cross section in Ref. [34]. We find that the radiative production of the excited a 1 and a 2 states with mass near 2 GeV=c 2 may give a 10% contribution to the measured M 3π spectrum above 1.8 GeV. This value is taken as an estimate of the systematic uncertainty associated with FSR at M 3π > 1.8 GeV=c 2 .

VI. DETECTION EFFICIENCY
The detection efficiency is determined using MC simulation as the ratio of the true 3π mass spectra computed after and before applying the selection criteria. The detection efficiency calculated in this way is shown in Fig. 15. Its mass dependence is fitted by a combination of a third-order polynomial in the range 0.62-2.3 GeV=c 2 , a linear function in the range 2.3-2.9 GeV=c 2 , and a constant above 2.9 GeV=c 2 . For the tighter requirement χ 2 3πγ < 20 the detection efficiency is smaller by 12-17%. The statistical uncertainty of the fitted detection efficiency is about  0.1% at the ω, 0.2% at the ϕ, and then increases to 1% at 2 GeV=c 2 and up to 2.2% at 2.5 GeV and above.
The decrease in efficiency below 0.62 GeV=c 2 is due to the merging of clusters from photons and charged pions in the calorimeter. This effect leads to π 0 loss, which increases as the 3π mass decreases. To avoid a possible systematic uncertainty due to imperfect simulation of this effect, we perform the measurement of the e þ e − → π þ π − π 0 cross section at masses above 0.62 GeV=c 2 .
The efficiency (ε MC ) found using MC simulation must be corrected to account for data-MC simulation differences in detector response: where δ i are efficiency corrections for the different effects discussed below.

A. ISR photon inefficiency
A correction is applied to the ISR photon detection efficiency. There are two sources of this correction: data-MC simulation differences in the probability of photon conversion in the detector material before the DCH, and dead calorimeter channels. A sample of e þ e − → μ þ μ − γ events is used to determine the calorimeter photon inefficiency in data. Events with exactly two charged tracks identified as muons are selected, and a one-constraint kinematic fit is performed with the requirement that the recoil mass against the muon pair be zero. A tight condition on the χ 2 of the kinematic fit selects events with only one photon in the final state. The photon direction is determined from the fit. The detection inefficiency is calculated as the ratio of the number of events not satisfying the condition E Ã γ > 3 GeV, to the total number of selected events. The same procedure is applied to simulated e þ e − → μ þ μ − γ events. The efficiency correction is determined from the data-MC simulation ratio as a function of the photon polar angle and the μ þ μ − invariant mass. The data-MC simulation difference in the probability of photon conversion is also studied using e þ e − → μ þ μ − γ events. In addition to two identified muons, we require that an event contain a converted-photon candidate, i.e., a pair of oppositely charged tracks with e þ e − invariant mass close to zero, momentum directed along the expected photon direction, and forming a secondary vertex well separated from the interaction region. The data-MC difference in the probability of photon conversion is measured as a function of the photon polar angle. Then we calculate the total correction to the ISR photon efficiency due to calorimeter inefficiency and photon conversion.
The measured angular dependence of the correction is used to reweight the simulated e þ e − → π þ π − π 0 γ events and calculate the efficiency correction. It is found to be −ð1.0 AE 0.2Þ% for 3π masses below 1.1 GeV=c 2 , −ð1.2 AE 0.2Þ% in the mass range 1.1-2.0 GeV=c 2 , and −ð1.4 AE 0.2Þ% in the range 2.0-3.5 GeV=c 2 . The contribution to this correction from photon conversion is about −0.2%.

B. π 0 efficiency and kinematic-fit χ 2 distribution
From the study of the ISR photon inefficiency it is expected that the difference between data and simulation in the π 0 detection efficiency is at least −2%. To study the π 0 losses more accurately, we perform a kinematic fit for data and simulated events to the e þ e − → π þ π − π 0 γ hypothesis using the measured parameters for only the two charged tracks and the ISR photon. The π 0 energy and angles are determined as a result of the fit. We apply a very tight condition on the fit quality and the background suppression conditions described in Sec. III. Because of the high level of remaining background, we restrict our study to the ω mass region.
The π 0 detection efficiency is determined as the fraction of selected signal events with a detected π 0 . The result depends on the definition of the π 0 candidate. For the simple π 0 definition as a pair of photons with invariant mass near the π 0 mass, for example, in the range 0.1-0.17 GeV=c 2 , there is a substantial probability to observe a false π 0 candidate due to a large number of spurious photons in an event. To avoid difficulties with false π 0 's, we require that an event containing the π 0 candidate satisfy our standard kinematic-fit condition χ 2 3πγ < 40. The 3π mass spectra for selected events with χ 2 3πγ < 40 and χ 2 3πγ > 40 are shown in Fig. 16. The mass spectra are fitted with a sum of distributions for signal and background events. The signal distribution is extracted from the simulation. The background spectrum is a sum of the simulated distribution for e þ e − → π þ π − π 0 π 0 γ events and a second-order polynomial with free coefficients. The efficiency correction due to π 0 losses is determined to be 5Þ%. Here f is the fraction of selected events with χ 2 3πγ < 40. In a similar way, we determine the efficiency correction in different ranges of the π 0 energy. At the current level of statistical precision the correction is found to be independent of the π 0 energy. Therefore, the efficiency correction due to π 0 losses determined at the ω region is also used for higher 3π masses.
The π 0 correction includes a part of the efficiency correction due to the χ 2 3πγ < 40 requirement related to the photons from the π 0 decay. To understand the influence of the data-simulation difference in the parameters of the charged tracks and the ISR photon, we study e þ e − → μ þ μ − γ events. We select events with two charged particles identified as muons and a photon with c.m. energy larger than 3 GeV. As mentioned in Sec. II, the simulation uses the requirement that the invariant mass of the muon pair and ISR photon be greater than 8 GeV=c 2 . To ensure compliance with this requirement in the data, we apply an additional condition that the invariant mass of the muon and the ISR photon candidates be greater than 9 GeV=c 2 .
For such selected events, a kinematic fit is performed with the requirements of energy and momentum balance. The fit uses measured momenta and angles of the muons and only angles of the ISR photon. The χ 2 distributions for selected data and simulated e þ e − → μ þ μ − γ events with μ þ μ − invariant mass M μμ < 1.1 GeV=c 2 are compared in Fig. 17. It is seen that the data and simulated distributions are in agreement. To estimate the difference between them numerically, we calculate the double ratio R 2 χ ¼ ½Nðχ 2 < cÞ=N 0 data = ½Nðχ 2 < cÞ=N 0 MC , where N 0 is the total number of selected μ þ μ − γ events, and Nðχ 2 < cÞ is the number of events satisfying the condition χ 2 < c. This ratio is practically independent of the c value in the range 20 < c < 40. Its deviation from unity, R 2 χ − 1, in the invariant mass range 0.6 < M μμ < 1.1 GeV=c 2 , equal to −ð0.4 AE 0.2Þ%, can be used as an estimation of the efficiency correction for e þ e − → π þ π − π 0 γ events.
We take into account the difference in the chargedparticle momentum distributions for the processes e þ e − → μ þ μ − γ and e þ e − → π þ π − π 0 γ. To understand a possible effect of this difference, we study the dependence of R 2 χ on the minimum muon momentum in an event and do not observe any statistically significant dependence. However, since the phase space distribution of charged pions from the e þ e − → π þ π − π 0 γ reaction cannot be fully reproduced using the e þ e − → μ þ μ − γ events, we assign a 100% systematic uncertainty to this correction.
In summary, the efficiency correction associated with the difference in the χ 2 distribution between data and simulation is estimated to be −ð0.4 AE 0.4Þ% in the 3π mass region 0.6-1.1 GeV=c 2 . For higher masses the correction is larger. Its average value in the mass range 1.1-3.5 GeV=c 2 is ð1 AE 1Þ%.  C. Efficiency correction due to the selection criteria Our preliminary selection contains the requirement of exactly two good quality charged tracks in an event. The definition of a good charged track is given in Sec. III. To determine an efficiency correction due to this requirement, we analyze events with three good tracks. Two of them with opposite charge having closest distance to the beam axis are selected as candidates for charged pions from the e þ e − → π þ π − π 0 γ reaction. The fraction of three-track events determined in the 3π mass regions near the ω and ϕ resonances is about 0.4% both in data and simulation. No efficiency correction due to the requirement of exactly two charged tracks is needed.
Radiative Bhabha events are rejected by the requirement that none of the good charged tracks be identified as an electron. The rejected events are prescaled by a factor of 40.
We study a sample of prescaled events passing our standard selection criteria, except for the electron identification requirement, and find that the efficiency correction is −ð0.01 AE 0.12Þ%.
The efficiency correction for the background suppression requirements described in Sec. III is determined near the ω, ϕ, and J=ψ resonances from ratios of the number of events selected with and without these requirements, in data and MC simulation. The fraction of signal events rejected by the background suppression requirements varies from 15% in the ω and ϕ mass region to 25% at the J=ψ. This dependence is reproduced by the simulation. The efficiency correction is ð0.4 AE 0.2Þ% at the ω and ϕ, and ð0.6 AE 0.8Þ% at the J=ψ. The latter correction is used in the energy region above 1.1 GeV.

D. Efficiency correction due to track losses
The data-MC simulation difference in track losses for isolated tracks is studied using e þ e − → τ þ τ − events, with one τ decaying leptonically and the other τ hadronically with three charged particles. No difference between data and simulation in the tracking efficiency is observed within an uncertainty of 0.24% per track. In e þ e − → π þ π − π 0 γ events, especially at small M 3π , the angle between charged tracks may be small, and the effect of track overlap in the DCH should be taken into account. To study this effect we analyze the distribution of the azimuthal angle difference between the positive and negative tracks Δφ In the BABAR magnetic field, events with Δφ AE > 0 exhibit a "fishtail" two-track configuration in which the tracks tend to overlap. The Δφ AE distribution for simulated signal events from the ω mass region is shown in Fig. 18. The track overlap leads to an asymmetry in the distribution. It should be noted that larger Δφ AE values correspond to larger differences between charged pion momenta. Therefore, the asymmetry in the distribution is seen even at relatively large Δφ AE ∼ 0.5. With larger values of M 3π , the Δφ AE distribution becomes wider and more symmetric. We estimate the fraction of events lost because of track overlap as This fraction is 11% at the ω, 8% at the ϕ, and about 1% at the J=ψ. We do not observe any significant difference in this fraction between data and simulation. The difference calculated over the 3π mass range 0.6-1.1 GeV=c 2 is −ð0.03 AE 0.23Þ%. The uncertainty in this difference is used as an estimate of the systematic uncertainty associated with track overlap in the DCH. The total systematic uncertainty in the detector efficiency due to data-MC simulation differences in the tracking efficiency and track overlap is estimated to be 0.5%.

E. Efficiency correction due to trigger and background filters
We also studied the quality of the simulation of the trigger and background filters [36] used in event reconstruction. In the analysis we use events passing through the two trigger lines L3OutDCH and L3OutEMC, which are based on signals from the DCH and EMC, respectively. The inefficiency of these lines in the simulation is 6.1% for L3OutDCH and ð2.3 AE 0.2Þ × 10 −4 for L3OutEMC. A logical OR of the L3OutDCH and L3OutEMC lines has a very small inefficiency, ð1.2 AE 0.1Þ × 10 −4 . The inefficiencies for the trigger lines in data can be estimated using the overlap of the samples of events passing through them. The simulation shows that these estimates are very close to the true inefficiencies of the trigger lines. This method applied to data results in an inefficiency of ð6.6 AE 0.1Þ% for L3OutDCH and ð3.8 AE 0.6Þ × 10 −4 for L3OutEMC. Although the efficiency in data is lower, the efficiency itself is very close to 100%. Therefore, no correction is applied for the trigger inefficiency. The inefficiency in the background filters in simulation is about 1.8% at the ω and ϕ mass regions and then decreases to 0.5% at 2 GeV=c 2 and to 0.3% at the J=ψ. To measure this inefficiency in data we use a subsample of prescaled events that does not pass through the background filters. The prescale factor is 200. The filter inefficiency for 3π masses below 1.1 GeV is measured to be ð3.2 AE 0.7Þ%. The efficiency correction in this mass region is −ð1.4 AE 0.7Þ%. For M 3π > 1.1 GeV=c 2 , insufficient statistical precision and large background do not allow us to determine the inefficiency with acceptable accuracy. Therefore, in this region we use the correction −ð1 AE 1Þ%, which covers the range of its possible variations as a function of mass.
The efficiency corrections δ i are summarized in Table I. The total efficiency correction is about −6%.
In Sec. VII, we also analyze the mass spectrum for events with χ 2 3πγ < 20. The additional correction related to this requirement is 0.001 AE 0.001 at the ω, 0.004 AE 0.004 at the ϕ, and 0.018 AE 0.007 at the J=ψ. A linear interpolation is used between the resonances.

F. Model uncertainty
The signal simulation uses the model of the ρð770Þπ intermediate state. This model works reasonably well for the ω and ϕ decays [9,10]. A comparison of the data and simulated two-pion distributions in different mass regions for BABAR 3π data was performed in our previous work [6]. Data and simulation agree well below 1.1 GeV, in the ω and ϕ regions. For higher masses, the difference was observed associated with additional intermediate mechanisms ωπ 0 and ρð1450Þπ.
To estimate the model dependence of the detection efficiency in the mass range 1.1-2 GeV=c 2 , the simulated signal events are reweighted using the model with a sum of the ρð770Þπ, ωπ 0 , and ρð1450Þπ mechanisms with coefficients and relative phases taken from the SND measurement [7]. The difference in the detection efficiencies between the two models depends on energy but does not exceed 1.5%. This number is taken as an estimate of the model uncertainty in the detection efficiency in the region 1.1-3.5 GeV=c 2 .
A similar procedure is used to find the correction to the detection efficiency at the J=ψ. Here we use the result of the Dalitz plot analysis of Ref. [11]. An ∼ 10% contribution from the ρð1450Þπ channel leads to a shift in the detection efficiency of −ð0.5 AE 0.1Þ%.

VII. FIT TO THE π + π − π 0 INVARIANT MASS DISTRIBUTION
To measure the e þ e − → π þ π − π 0 cross section, detector resolution effects need to be unfolded from the measured 3π invariant-mass spectrum. In Fig. 19 (left) the simulated distribution of the true 3π mass in the energy regions of the ω and ϕ resonances is compared with the distribution of the reconstructed 3π mass. The true spectrum varies by four orders of magnitude and has two narrow peaks. The reconstructed spectrum strongly differs from the true one. For such a spectrum, the result of the unfolding procedure is very sensitive to the quality of simulation used to obtain the resolution function. To study the difference between data and simulation in resolution, we fit the measured 3π mass spectrum with the vector-meson-dominance model including several resonances. The ω and ϕ masses and widths are known with relatively high accuracy. Therefore, from the fit we can extract the mass shift and standard deviation of an additional smearing Gaussian function needed to describe the data-MC simulation difference in the mass resolution. These parameters are determined separately for the ω and ϕ resonances.
The detector resolution function has long non-Gaussian tails as seen in Fig. 19 (right), where the distribution of the difference between the reconstructed and true mass (ΔM 3π ) is shown for events from the ω peak. To increase the fraction of events in the non-Gaussian tails, events are selected with the condition 20 < χ 2 3πγ < 40. The distribution is fitted by a sum of three Gaussians and a Lorentzian function LðxÞ ¼ ðγ=πÞ=ððx − x 0 Þ 2 þ γ 2 Þ. The latter is shown in Fig. 19 (right) by the dashed histogram. Because of the asymmetry in the ΔM 3π distribution, the maximum of the Lorentzian function is shifted from zero by about −30 MeV. The same shift is observed in the ΔM 3π distribution at 750, 900 GeV=c 2 , and at the ϕ resonance.
To describe a possible difference between data and simulation in the tails of the resolution function, we introduce to the fit to the M 3π data spectrum a smearing Lorentzian function with x 0 ¼ −30 MeV.
The following probability density function is used in the fit to the measured M 3π spectrum: dN dm where the theoretical spectrum of true 3π mass (dN=dm) is convolved with the smearing Gaussian (G) and Lorentzian (L) functions. The spectra of true and measured masses are presented as histograms with the same binning. The folding matrix P ij obtained using simulation gives the probability that an event with true mass in bin j is actually reconstructed in bin i. From the fit we determine the standard deviations of the smearing Gaussian and mass shifts at the ω and ϕ, ϵ, and γ of the smearing Lorentzian. The width of the mass bins near the ω and ϕ resonances is chosen to be 2.5 MeV=c 2 . This width is not much smaller than the resonance widths. Therefore, the elements of the matrix P ij depend on the values of the resonance parameters used in simulation. We correct the folding matrix using an iterative procedure. The procedure uses results of the fit without Lorentzian smearing (Model 4 in Table II) for events with χ 2 3πγ < 20. Simulated events are reweighted by the ratio of the fitted spectrum ðdN=dmÞ Ã G to the true simulated spectrum. The reweighting is performed with a bin width of 0.5 MeV=c 2 . Then a new matrix P ij is obtained, and the fit is repeated. We iterate until the change in ðdN=dmÞ true Ã G between two successive iterations is less than 0.1%.
The true mass spectrum in the fit is described by the following function:

Model Lorentzian smearing
Bðρ → 3πÞ where σ 3π ðmÞ is the Born cross section for e þ e − → 3π, dL=dm is the so-called ISR differential luminosity, ε is the detection efficiency as a function of mass, and R is a radiative correction factor accounting for the Born mass spectrum distortion due to emission of several photons by the initial electron and positron. The ISR luminosity is calculated using the total integrated luminosity L and the probability density function for ISR photon emission [Eq. (2)]: Here, x ¼ 1 − m 2 =s, ffiffi ffi s p is the e þ e − c.m. energy, C ¼ cos θ 0 , and θ 0 determines the range of polar angles in the c.m. frame: θ 0 < θ γ < 180°− θ 0 for the ISR photon. In our case θ 0 is equal to 20°, since we determine the detector efficiency using the simulation with 20°< θ γ < 160°. The total integrated luminosity (L ¼ 468.6 fb −1 ) is measured with an accuracy of 0.43% [37].
The Born cross section for e þ e − → 3π can be written as the sum of the contributions of five resonances ρ ≡ ρð770Þ, ω ≡ ωð782Þ, ϕ ≡ ϕð1020Þ, ω 0 ≡ ωð1420Þ, and ω 00 ≡ ωð1650Þ: where m V and Γ V are the mass and width of the resonance V, φ V is its phase, and BðV → e þ e − Þ and BðV → 3πÞ are the branching fractions of V into e þ e − and 3π, Here Γ f ðmÞ is the mass-dependent partial width of the resonance decay into the final state f, and Γ f ðm V Þ ¼ Γ V BðV → fÞ. The mass-dependent width for the ω and ϕ mesons has been calculated taking into account all significant decay modes. The corresponding formulas can be found, for example, in Ref. [3]. We assume that the V → 3π decay proceeds via the ρπ intermediate state, and F ρπ ðmÞ is the 3π phase space volume calculated under this hypothesis. The formula for the F ρπ calculation can be found in Ref. [3]. The radiative correction factor R is determined using Monte Carlo simulation (at the generator level, with no detector simulation) with the PHOKHARA event generator [38]. This generator includes the next-to-leading order (NLO) ISR contributions. The accuracy of the cross section calculation for ISR processes with the ISR photon emitted at large angle is estimated to be 0.5% [39]. Since the radiative correction is independent of process, we generate e þ e − → μ þ μ − γ events with no FSR included in the simulation, and calculate a ratio of the mass spectra obtained in the NLO and LO generator modes. With the requirement on the invariant mass of the μ þ μ − γ system M μμγ > 8 GeV=c 2 , this ratio is weakly dependent on mass and is equal to R ¼ 1.0077 AE 0.0005 below 1.1 GeV=c 2 , 1.0086 AE 0.0004 between 1.1 and 2 GeV=c 2 , and 1.0091 AE 0.0004 in the range 2-3.5 GeV=c 2 . The quoted uncertainty reflects the observed R variation in the specified mass range. The radiative correction factor does not include the corrections due to leptonic and hadronic vacuum polarization. Here we follow the generally accepted practice [40] of including the vacuum polarization correction to the resonance electronic width.
The free parameters in the fit are the scale factors for products of the branching fractions BðV → e þ e − ÞBðV → 3πÞ, and the masses and widths of the ω 0 and ω 00 . The masses and widths of the ω and ϕ mesons are fixed at the Particle Data Group (PDG) values [35].
The phase φ ω is set to zero. The relative phase between the ω and ϕ amplitudes, φ ϕ ¼ ð163 AE 7Þ°, is taken from Ref. [3]. The phases of the ω 0 and ω 00 are fixed at values of 180°and 0° [41] with an uncertainty of 20°. This uncertainty is estimated from the deviation of φ ϕ from 180°. Our fitting function does not take into account the isovector e þ e − → ωπ 0 γ → 3πγ contribution and the presence of the ω 00 → ρð1450Þπ decay. Therefore, we do not expect that the parameters of excited ω states are determined correctly. Their inclusion into the fit is needed to study the effect of interference of the ω and ϕ amplitudes with the contributions of the excited states. The fitted mass region is restricted to masses below 1.8 GeV.
The branching fraction of the ρ → π þ π − π 0 decay can be estimated assuming that the dominant mechanism of the ω → π þ π − and ρ → π þ π − π 0 decay is ρ − ω mixing. Under this assumption, the coupling g ρ→3π ¼ ξg ω→3π , where the mixing parameter jξj 2 ≈ Γðω → 2πÞ=Γðρ → 2πÞ, and Bðρ → 3πÞ ≈ 0.4 × 10 −4 . The phase φ ρ is expected to be close to −90°. A significantly larger value of ξ is obtained in Ref. [42], where data on the pion electromagnetic form factor are analyzed in the model including both ρ − ω mixing and direct isospin-breaking ω → 2π decay. Using the result of Ref. [42], we obtain Bðρ → 3πÞ ≈ ð2.5 AE 1.0Þ × 10 −4 and φ ρ ¼ −ð40 AE 13Þ°. The values of the branching fraction and phase measured in the SND experiment [3] Table II, respectively. current experimental accuracy, ∼1% at the ω resonance, the contribution of the ρ → π þ π − π 0 decay with a branching fraction of about 10 −4 must be taken into account in the fit. We fit the 3π mass spectrum with a series of 4 models, allowing or not for Lorentzian smearing and leaving the ρ → 3π branching fraction as a free parameter or forced to zero, as shown in Table II. The results of the fit in Models 1,2,3 are shown in Fig. 20. The smearing Gaussian standard deviation obtained in Model 1 is 1.5 AE 0.2 GeV=c 2 at the ω and 1.5 AE 0.4 GeV=c 2 at the ϕ. Therefore, in what follows we use a mass-independent smearing standard deviation. The fitted parameters of the smearing Lorentzian function are the following: ϵ ¼ 0.007 AE 0.002 and γ ¼ 63 AE 35 GeV=c 2 . The physical fit parameters will be discussed below.
It is seen from Fig. 20 that the fit in Model 2 [Bðρ → 3πÞ ¼ 0, no smearing Lorentzian] cannot describe data well below 0.73 GeV=c 2 and in the region 0.82-0.9 GeV=c 2 . Including the smearing Lorentzian function (Model 3) improves the fit in the energy region below the ω. The region 0.82-0.9 GeV=c 2 cannot be described reasonably well without the ρ → 3π decay. Models 2 and 3 also have a worse fit quality in the mass range 1.05-1.8 GeV=c 2 . This is because the fit tries to compensate the absence of the ρ → 3π decay by increasing the contribution from the tails of the ωð1420Þ and ωð1650Þ resonances.
The difference between Model 1 and Model 4 [free Bðρ → 3πÞ, no Lorentzian smearing] is maximal in the 3π mass region 0.62-0.72 GeV=c 2 (see Fig. 21). To decrease the difference between the measured and predicted spectrum below the ω peak in absence of the Lorentzian smearing, the fit in Model 4 increases Bðρ → 3πÞ by 14%.
To reduce the influence on the fitted parameters of the data-simulation difference in resolution, we tighten the condition on χ 2 3πγ from 40 to 20. In the last four rows of Table II, we compare the quality of the fit with Models 1 − 4 to the mass spectrum obtained with the tighter χ 2 3πγ requirement. It is seen that inclusion of the Lorentzian smearing in this case improves the fit quality insignificantly. The obtained parameters of the Lorentzian smearing function are ϵ ¼ 0.0022 AE 0.0016 and γ ¼ 59 þ54 −31 MeV=c 2 . Therefore, below, we quote the fit parameters for Model 4. The standard deviation of the smearing Gaussian function and the mass shifts for the ω and ϕ mesons are found to be  Table II, respectively.   TABLE III. Contributions to the systematic uncertainties in fit parameters from different effects ( Effect and m ϕ − m PDG ϕ ¼ 0.08 AE 0.08 MeV=c 2 . The latter two parameters are consistent with zero. Since the mass resolution (full width at half-maximum is about 13 MeV=c 2 at the ω and 15 MeV=c 2 at the ϕ) is larger than the ω and ϕ widths, the parameter with least sensitivity to resolution effects is the area under the resonance curve, i.e., P V ≡ ΓðV → e þ e − ÞBðV → π þ π − π 0 Þ. From the fit we obtain P ω ¼ ð0.5698 AE 0.0031 AE 0.0082Þ keV; For the ρ meson, we determine Bðρ → 3πÞ ¼ ð0.88 AE 0.23 AE 0.30Þ × 10 −4 ;  [4,5] data on the e þ e − → π þ π − π 0 cross section and the cross section calculated using Eq. (11) with parameters obtained from the fit to the 3π mass spectrum. The uncertainties shown for the SND and CMD2-2 data are statistical. The systematic uncertainty is 3.4% for the SND data at the ω [3], 1.3% for the CMD-2 data at the ω [4], 5% for the SND data at the ϕ [3], and 2.5% for the CMD-2 data at the ϕ [5]. The systematic uncertainty in the BABAR cross section is about 1.5%.
The significance of the ρ → 3π decay estimated from the difference between the χ 2 values for Models 4 and 2 is greater than 6σ. The first uncertainties in Eqs. (13) and (14) are statistical, while the second are systematic. The latter include the systematic uncertainties in the ISR luminosity, radiative correction, and detection efficiency. The uncertainty due to the data/MC difference in the mass resolution line shape is estimated as a difference between results of the fits with Models 4 and 1. We also vary within the uncertainties the values of the parameters φ ϕ , Γ ω , and Γ ϕ , and the scale factors for the background processes. These contributions to the systematic uncertainties are listed in Table III. Our fitting model given by Eq. (11) assumes that the e þ e − → π þ π − π 0 process proceeds via the ρð770Þπ intermediate state. Actually, due to a sizable ωð1650Þ → ρð1450Þπ → 3π transition [7] and the existence of the e þ e − → ωπ 0 → π þ π − π 0 process, the e þ e − → π þ π − π 0 amplitude above the ϕ cannot be presented as a simple coherent sum of the ωð1420Þ and ωð1650Þ amplitudes. To study the effect of non-ρð770Þπ mechanisms, we substitute the BABAR data above 1.1 GeV=c 2 by the SND data on the ρð770Þπ and ρð1450Þπ cross sections [7]. The data on the phase difference between the ρð770Þπ and ρð1450Þπ amplitudes measured in Ref. [7] is also included in the fit. The new fitting function takes into account ωð1650Þ transitions to ρð770Þπ and ρð1450Þπ, and interference between ρð770Þπ and ρð1450Þπ amplitudes. This new approach modifies the contribution of the ωð1420Þ and ωð1650Þ resonances in the 3π mass region below 1.1 GeV=c 2 and shifts the parameters of the ρ, ω and ϕ resonances. In particular, for the ρ decay we obtain Bðρ → 3πÞ ¼ ð1.14 þ0. 32 −0.28 Þ × 10 −4 : ð15Þ The difference between the results of this new fit and our nominal fit is used as an estimate of systematic uncertainty due to the ωð1650Þ → ρð1450Þπ decay (see Table III.) The fitted values of ΓðV → e þ e − ÞBðV → 3πÞ for the ω and ϕ mesons given by Eq. (13) are in reasonable agreement with the corresponding world average values [35]: 0.557 AE 0.011 keV and 0.1925 AE 0.0043 keV, respectively. For the ω meson the accuracy of our result is comparable with the accuracy of the PDG value. For the ϕ meson, we have a large systematic uncertainty related to the interference between ϕ-meson amplitude and amplitudes of the resonances of the ω family. The fitted values of Bðρ → 3πÞ and φ ρ given by Eq. (14) are in agreement with the SND results: ð1.01 þ0.54 −0.36 AE 0.34Þ × 10 −4 and −ð135 þ17 −13 AE 9Þ° [3]. It is instructive to compare the cross section calculated using Eq. (11) with the SND and CMD-2 data [2][3][4][5]. Such a comparison is presented in Fig. 22, where the difference between SND and CMD-2 data and the BABAR fit is shown in the energy region of the ω and ϕ resonances. A shift in the energy (3π mass) scale between different sets of data leads to the appearance of wiggles in the relative difference between them near the resonance maximum. To eliminate these wiggles we shift the SND (CMD-2) data by −0.18 (0.09) MeV at the ω region, and 0.09 (0.13) MeV at the ϕ region. It is seen that the BABAR cross section is in reasonable agreement with the SND data below the ϕ. At the ω the difference between the SND and BABAR cross sections is about 2%, well below the systematic uncertainty (3.4% for SND and 1.3% for BABAR). The CMD-2 data in the vicinity of the ω lie about 7% below zero. With the CMD-2 statistical and systematic uncertainties of 1.8% and 1.3%, respectively, the difference between CMD-2 and BABAR is 2.7σ. Near the maximum of the ϕ-meson resonance the CMD-2 and SND data with systematic uncertainties of 2.5% and 5%, respectively, lie about 4% and 11% higher than the BABAR cross section.
VIII. MEASUREMENT OF THE e + e − → π + π − π 0 CROSS SECTION BELOW 1.1 GeV=c 2 In the M 3π region below 1.1 GeV=c 2 , the detector resolution strongly distorts the 3π mass spectrum as shown in Fig. 19 (left). To obtain the true mass (M true 3π ) spectrum, unfolding must be applied to the measured M 3π spectrum. Similar to the previous BABAR analyses [43,44], we use a simplified version of the iterative unfolding method developed in Ref. [45].
In Sec. VII we reweight the signal MC simulation using the results of the fit to the measured M 3π spectrum and obtain the folding matrix P ij . This matrix must be corrected to take into account the data-MC difference in the mass resolution. This difference is described by the smearing Gaussian (G) and Lorentzian (L) functions, the parameters of which are determined from the fit described in Sec. VII. The corrected folding matrix is calculated as where the matrices G ij and L ij are obtained using the fitted theoretical mass spectrum ðdN=dmÞ true FIT and its convolution with the smearing functions G and L, respectively. In the unfolding procedure described below we use P Ã ij with ϵ ¼ 0. Unfolding with nonzero ϵ ¼ 0.002 (Model 1 in Table II) is performed to estimate a systematic uncertainty due to possible unaccounted Lorentzian smearing. Figure 23 represents the transfer matrix A ij ¼ P Ã ij T j , where the vector T i is obtained by integration of the theoretical mass spectrum ðdN=dmÞ true FIT over bin j. The unfolding matrix can also be obtained asP ij T j is the reconstructed spectrum corresponding to the true spectrum T i . For small bin size the folding matrix describes detector resolution and FSR effects and does not depend on the true spectrum T i , while A ij andP ij depend on it. The unfolding method used is based on the idea that if T i is close to the true spectrum and the folding matrix describes resolution and FSR effects well, the matrixP ij can be applied to the measured spectrum to obtain the true spectrum.
The unfolding process consists of several iteration steps. At each step, differences between the unfolded data spectrum and T i are used to correct T i and the unfolding matrixP ij , keeping the folding probabilities unchanged. A regularization function is used to suppress unfolding large statistical fluctuations in the data and guarantee the stability of the method.
Since the transfer matrix shown in Fig. 23 is nondiagonal, the values obtained for the true data spectrum are correlated. The covariance matrix containing the statistical uncertainties and their bin-to-bin correlations is obtained from pseudoexperiments (toy MC), where both the spectrum and the transfer matrix are statistically fluctuated. In this analysis, we generate 1000 toy-MC samples.
To test the unfolding procedure and choose parameters of the regularization functions, we examine two model spectra, each with the number of events equal to the number of events in data. The first spectrum is the true MC spectrum T i , while the second is based on the fit with zero ρð770Þ → 3π amplitude (Model 2 in Table II). Both spectra are convolved with the folding matrix P Ã ij , statistically fluctuated, and then subjected to the same unfolding procedure. In Fig. 24 the relative difference between the two unfolded spectra is compared with the same difference in the true spectra. The regularization parameter is chosen to minimize the difference between the points and the curve. The unfolded data shown in Fig. 24 are obtained after the first iteration step, and further iterations do not improve the result. Another test is performed to assess a systematic uncertainty in the unfolding method. A set of 100 spectra are generated as described above, using the true spectrum T i . They are unfolded and averaged. The deviation of the  IV. Measured e þ e − → π þ π − π 0 cross section below 1.1 GeV=c 2 . The first uncertainty is the square root of the diagonal element of the statistical covariance matrix, the second is correlated systematic, and the third is the square root of the diagonal element of the systematic covariance matrix.  (Table continued) average unfolded spectrum from T i is taken as a measure of the systematic uncertainty in the unfolding method. Figure 25 shows the difference between the unfolded data spectrum and the result of fit to the measured 3π mass spectrum (vector T i ). The error bars correspond to the diagonal elements of the covariance matrix for the unfolded spectrum given in Ref. [46]. The comparison demonstrates good agreement of fit results and unfolding and establishes the adequacy of the model used in the fit.
Using the unfolded 3π mass spectrum and Eq. (9), we calculate the Born cross section listed in Table IV, where the first uncertainty is the square root of the diagonal element of the statistical covariance matrix [46]. Systematic uncertainty is divided into two parts. The second error in Table IV represents a correlated uncertainty that includes the uncertainties in the luminosity, radiative correction, detection efficiency, and the uncertainty due to the unfolding procedure. For the remaining part of the systematic uncertainty associated with background subtraction and data-simulation difference in the mass resolution, we provide the covariance matrix [46]. The square root of the diagonal element of this matrix is listed in Table IV as the third error.
The mass dependence of the total systematic uncertainty is compared with the uncertainties from the different sources in Fig. 26. It is seen that in the mass region between 0.73 and 1.03 GeV=c 2 the systematic uncertainty is dominated by the uncertainties in the luminosity, radiative correction, and detection efficiency, whose total contribution (1.3%) is independent of mass.
IX. MEASUREMENT OF THE e + e − → π + π − π 0 CROSS SECTION ABOVE 1.1 GeV=c 2 Above 1.1 GeV=c 2 , the resolution effects distort the 3π mass spectrum insignificantly. We test this by a convolution of the theoretical mass spectrum (9) in the mass range 1-2 GeV=c 2 with the resolution function obtained using simulation. The observed difference between the true and measured spectra does not exceed 1%. Therefore, the e þ e − → π þ π − π 0 in the mass region 1.1-3.5 GeV=c 2 is determined as  The cross section thus obtained is listed in Table V. The quoted uncertainties are statistical and systematic. The latter includes uncertainties in the integrated luminosity (0.4%) and radiative correction (0.5%), the statistical (0.3-2.4%), systematic (1.7-1.8%), and model (1.5%) uncertainties in the detection efficiency, and the uncertainty associated with background subtraction (3-15%).
In Fig. 27 (left) the measured cross section is compared with the SND measurement [7] in the mass range 1.1-2 GeV=c 2 . A sizable difference between the two measurements is observed near 1.25 GeV=c 2 and 1.5 GeV=c 2 . The cross section above 2 GeV is shown in Fig. 27 (right).
The leading-order hadronic contribution to the muon anomalous magnetic moment is calculated using the measured total hadronic cross section via the dispersion integral (see, for example, Ref. [47]) where the kernel function KðsÞ can be found in Ref. [47] and TABLE V. Measured e þ e − → π þ π − π 0 cross section above 1.1 GeV=c 2 . The first uncertainty is statistical, the second is systematic. In the M 3π intervals 3.0-3.1 GeV=c 2 and 3.1-3.2 GeV=c 2 , the values of the nonresonant cross section are listed, which are obtained by subtraction of the J=ψ contribution (see Sec. XI).
Here σ 0 is the bare cross section, excluding effects from vacuum polarization.
To calculate a 3π μ we substitute σ 0 ðe þ e − → hadronsÞðsÞ in Eq. (19) by where σ 3π ðsÞ is the "dressed" cross section measured in Secs. VIII and IX. The vacuum polarization operator ΠðsÞ is tabulated in Ref. [48]. The integral (18) is substituted by a sum over mass bins with s ¼ M 2 3π;i . In the sum, the values of the functions K, dL=dm, ε, and j1 − ΠðsÞj 2 are taken at the center of the bin. To estimate the uncertainty due to the substitution of the integral by the sum, we perform a 3π μ calculations using the theoretical cross section [Eq. (11)] and mass spectrum [Eq. (9)] with parameters [Eq. (14)]. The difference between the sum and integral is found to be 0.03% for the mass range 0.62-1.1 GeV=c 2 and 0.007% for the mass range 1.1-2.0 GeV=c 2 . The main reason for the larger difference in the lower-mass region is the strong s dependence of j1 − ΠðsÞj 2 near the ω and ϕ resonances.
It should be noted that the exclusion/inclusion of the factor j1 − ΠðsÞj 2 in the a 3π μ calculation changes its value by about 3.5%. The theoretical cross section [Eq. (11)] is used to estimate the uncertainty in a 3π μ associated with ΠðsÞ uncertainties [ΔΠðsÞ] given in Ref. [48]. Assuming that the ΔΠðsÞ at different s are fully correlated, we calculate a 3π μ with ΠðsÞ substituted by ΠðsÞ AE ΔΠðsÞ. The resulting uncertainty in a 3π μ is found to be 0.06% for the mass range 0.62-1.1 GeV=c 2 and 0.03% for the mass range 1.1-2.0 GeV=c 2 .
The above-mentioned strong s dependence of the j1 − ΠðsÞj 2 term near the resonances leads to the a 3π μ systematic uncertainty associated with the 3π mass scale calibration. The e þ e − → hadrons cross section near the ω and ϕ resonances used for the ΠðsÞ calculation in Ref. [48] is based mainly on data obtained in the SND and CMD-2 experiments at the VEPP-2M collider. In Sec. VII, where the SND and CMD-2 measurements are compared with the BABAR fit, we observe 0.1-0.2 MeV=c 2 shifts between the energy/mass scales of the BABAR and VEPP-2M experiments. To estimate the associated systematic uncertainty we introduce a mass shift ΔM ¼ AE0.2 MeV=c 2 in the theoretical cross section [Eq. (11)] and calculate a 3π μ . The relative difference with the zero ΔM value is found to be 0.2% for the mass range 0.62-1.1 GeV=c 2 and 0.03% for the mass range 1.1-2.0 GeV=c 2 .
Combining the three systematic uncertainties described above in quadrature, we find the systematic uncertainty in a 3π μ associated with the vacuum polarization factor to be 0.21% at M 3π ¼ 0.62-1.1 GeV=c 2 and 0.04% at M 3π ¼ 1.1-2.0 GeV=c 2 .
The a 3π μ values for different mass intervals obtained using the e þ e − → π þ π − π 0 cross section measured in this work are listed in Table VI. For the mass range 0.62-1.10 GeV=c 2 the quoted uncertainties are statistical, systematic due to the cross section measurement, and For the BABAR data, the error bar represents the statistical uncertainty, while the shaded box shows the systematic uncertainty. For the SND data, only the statistical uncertainty is shown; the systematic uncertainty is 4.4%.
systematic due to the vacuum polarization. The statistical uncertainty in a 3π μ is calculated using the toy MC simulation as described in Sec. VIII.
The contributions to the systematic uncertainty in a 3π μ ð0.62 < M 3π < 1.1 GeV=c 2 Þ from different effects are listed in Table VII. The uncertainties in the detection efficiency, luminosity, and radiative correction dominate. These three contributions are common for the mass intervals below and above 1.1 GeV=c 2 . However, in the mass range 1.10-2.00 GeV=c 2 the largest contribution to the systematic uncertainty comes from the FSR background. In the a 3π μ ð1.1 < M 3π < 2.0 GeV=c 2 Þ calculation, the systematic uncertainties listed in Table V are conservatively taken to be 100% correlated.
XI. MEASUREMENT OF THE J=ψ → π + π − π 0 DECAY The 3π mass spectrum in the J=ψ mass region for data events selected with the standard criteria is shown in Fig. 28. The small width of the J=ψ resonance leads to negligible peaking background. In particular, e þ e − → J=ψγ → K þ K − π 0 γ events reconstructed under the 3πγ hypothesis have the 3π invariant mass in the range 2.8-3.0 GeV=c 2 . To determine the number of J=ψ events, the spectrum is fitted with a sum of a resonance distribution and a linear background. The resonance line shape is a Breit-Wigner function convolved with a triple-Gaussian function describing detector resolution. The Breit-Wigner width is fixed at its PDG value [35]. The parameters of the resolution function are determined from simulation. To account for possible differences in detector response between data and simulation, the simulated resolution function is modified by adding a smearing variance σ 2 s to each of the three variances of the triple-Gaussian function. The free parameters in the fit are the number of resonance events (N J=ψ ), the number of nonresonant background events, the slope of the background, σ s , and the resonance mass.
The differential cross section for ISR production of a narrow resonance, such as J=ψ, can be calculated using [12] TABLE VI. Values of a 3π μ for different mass intervals. The first three rows represent the BABAR result, while the last three are the calculations [1,[49][50][51] based on previous e þ e − → π þ π − π 0 measurements.

XII. SUMMARY
The cross section for the process e þ e − → π þ π − π 0 has been measured by the BABAR experiment in the c.m. energy range from 0.62 to 3.5 GeV, using the ISR method. The cross section is dominated by the ω and ϕ resonances. Near the maxima of these resonances it is measured with a systematic uncertainty of 1.3%. The leading-order hadronic contribution to the muon magnetic anomaly, calculated using the measured e þ e − → π þ π − π 0 cross section from threshold to 2.0 GeV, is ð45.86 AE 0.14 AE 0.58Þ × 10 −10 . Our a 3π μ value is in reasonable agreement with the calculations [1,[49][50][51] based on previous e þ e − → π þ π − π 0 measurements but is more precise by about a factor of about 2. From the fit to the measured 3π mass spectrum in the process e þ e − → π þ π − π 0 γ we have determined the resonance parameters P ω ¼ ð0.5698 AE 0.0031 AE 0.0082Þ keV; P ϕ ¼ ð0.1841 AE 0.0021 AE 0.0080Þ keV; Bðρ → 3πÞ ¼ ð0.88 AE 0.23 AE 0.30Þ × 10 −4 ; where P V ¼ ΓðV → e þ e − ÞBðV → π þ π − π 0 Þ. The significance of the ρ → 3π decay is found to be greater than 6σ.
The measured values of ΓðV → e þ e − ÞBðV → 3πÞ for the ω and ϕ mesons are in agreement with the world average values [35]. For the ω meson, the accuracy of our result is comparable with the accuracy of the PDG value. For the ϕ meson we have a large systematic uncertainty related to the interference between the ϕ-meson amplitude and amplitudes of the resonances of the ω family. The measured values of Bðρ → 3πÞ and φ ρ are in agreement with the SND results [3].