Measurement of the $e^+e^- \to \pi^+\pi^-\pi^0\eta$ cross section below $\sqrt{s}=2$ GeV

The process $e^+e^- \to \pi^+\pi^-\pi^0\eta$ is studied in the center-of-mass energy region below 2 GeV with the SND detector at the VEPP-2000 $e^+e^-$ collider. The four intermediate states contribute to this process: $\omega\eta$, $\phi\eta$, $a_0(980)\rho$, and a structureless $\pi^+\pi^-\pi^0\eta$ state. We measure the total $e^+e^- \to \pi^+\pi^-\pi^0\eta$ cross section and the cross sections for its components: $\omega\eta$, $\phi\eta$, and a sum of $a_0(980)\rho$ and the structureless state. Our results are in agreement with previous measurements and have comparable or better accuracies.


I. INTRODUCTION
The main goal of experiments at the VEPP-2000 e + e − collider [1] is the measurement of the total cross section of e + e − annihilation into hadrons.The information on the cross section is necessary, in particular, for calculation of the hadronic vacuum polarization contribution into the muon anomalous magnetic moment and the running electromagnetic coupling * e-mail: A.A.Botov@inp.nsk.suconstant.Below 2 GeV the total hadronic cross section is determined as a sum of exclusive cross sections for all possible hadronic modes.At the moment most of these cross sections have been measured.At the same time still there are processes, which cross sections have not been measured yet or have been measured with insufficient accuracy: e + e − → π + π − π 0 η, π + π − 3π 0 , π + π − 2π 0 η, π + π − 4π 0 , etc.
The e + e − → φη process was studied in the BABAR experiment in several φ and η decay modes [7][8][9], and in the SND experiment in the K + K − γγ final state [10].The study of the e + e − → π + π − π 0 η process was performed by the CMD-3 detector [4], which measured for the first time the total e + e − → π + π − π 0 η cross section in the energy region E < 2 GeV and the cross sections for the subprocesses (4) and (5).

II. DATA AND SIMULATION
In this work we analyze a data sample with an integrated luminosity of 27 pb −1 recorded by the SND detector [2] at the VEPP-2000 e + e − collider in 2011 and 2012.In the energy range under study, 1.34-2.00GeV, data were collected in 36 energy points.The c.m. energies for these points were determined with accuracy of 2-6 MeV by the CMD-3 detector, which collected data simultaneously with SND, using momentum measurements in Bhabha and e + e − → pp events [13].Because of the absence of any narrow structures in the cross sections under study, the 36 energy points are merged into 13 energy intervals listed in Table I.The weighted average energy for each interval also listed in the Table I is calculated taking into account the luminosity distribution over the energy points entering into the interval.
Simulation of the processes e + e − → ωη and φη is done with Monte Carlo (MC) event generators based on the e + e − → V P model with V → ρπ → π + π − π 0 , where V is a vector meson and P is the η-meson.In simulation of the a 0 (980)ρ intermediate state, it is assumed that the a 0 (980) and ρ are produced in the S-wave.The simulation of the structureless process e + e − → nres → 3πη is performed in the hypothesis e + e − → ρ(1450)π with ρ(1450) → ρη.The event generators take into account radiative corrections to the initial particles calculated according to Ref. [15].The angular distribution of additional photons radiated by the initial electron and positron is simulated according to Ref. [16].The energy dependencies of Born cross sections needed for calculations of the radiative corrections are obtained from data using an iterative procedure.As a first approximation we take the cross sections from Ref. [6] for the ωη, from Ref. [7] for the φη, and from Ref. [4] for the a 0 ρ and nres intermediate states.In this work we also measure the total e + e − → π + π − π 0 η cross section and the cross section for the sum of the a 0 ρ and nres intermediate states.To determine the energy dependence of the detection efficiency for these processes, simulation in the e + e − → ρ(1450)π, ρ(1450) → ρη model is used.The initial Born cross sections are taken from Refs.[3,4].The iterative procedure is following.From the MC simulation we determinate the detection efficiency and measure the Born cross section from our data as described in Sec.VII.The obtained cross section is used in the second-iteration simulation.
We do not perform full simulation, but instead reweight simulated events using the new Born cross section.The iterations are stopped when the relative difference in the efficiency between two successive iterations is less than 1%.
Interactions of the generated particles with the detector materials are simulated using GEANT4 software [17].The simulation takes into account variation of experimental conditions during data taking, in particular dead detector channels and beam-induced background.The beam background leads to the appearance of spurious photons and charged particles in detected events.To take this effect into account, simulation uses special background events recorded during data taking with a random trigger, which are superimposed on simulated events.
The integrated luminosity is measured using Bhabha scattering events with a 1% systematic uncertainty [6].The luminosity for the 13 energy intervals is listed in Table I.

III. EVENT SELECTION
The process e + e − → π + π − π 0 η is analyzed in the decay mode η → γγ.Therefore, we select events with • two or three charged tracks originated from the interaction region; • at least four photons with energy greater than 20 MeV; • the energy deposition in the calorimeter greater than 300 MeV.
For selected events the vertex fit is performed using the parameters of two charged tracks.
The quality of the fit is characterized by the parameter χ 2 r .If there are three charged tracks in an event, the two of them with the lowest χ 2 r value are selected.The found vertex is used to refine the measured angles of charged particles and photons.Then we select events containing at least one π 0 candidate and one η candidate, which are defined as two-photon pairs with invariant masses in the windows 70 < m 12 < 200 MeV/c 2 and 400 < m 34 < 700 MeV/c 2 , respectively.For these events, a kinematic fit to the π + π − π 0 γγ hypothesis is performed with the four constraints of energy and momentum balance, and the π 0 mass constrain for the π 0 candidate.The χ 2 of the fit (χ 2 3π2γ ) is required to be less than 30.For events containing more than one combination of π 0 and η candidates, the combination with the smallest value of χ 2 3π2γ is chosen.The photon parameters after the kinematic fit are used to recalculate the η-candidate invariant mass (M η ).The event is refitted with the η-mass constraint.The refined energy of the η-meson candidate is used to calculate the invariant mass of the system recoiling against the η meson (M rec η ).The main background source for the process under study is the process For its suppression, a kinematic fit to the e + e − → π + π − π 0 π 0 (γ) hypothesis is performed, in which radiation of an additional photon along the beam axis is allowed.The condition χ 2 4π(γ) < 200 is applied, where χ 2 4π(γ) is χ 2 of the kinematic fit.It decreases the number of background e + e − → π + π − π 0 π 0 events by a factor of 10, rejecting about 40% of the signal events.

IV. DETERMINATION OF THE NUMBER OF SIGNAL EVENTS
The M η spectrum for the selected 13113 data events is shown in Fig. 2. It is seen that only about 35% of events contain an η meson.To determine the number of signal events the M η spectrum in each energy interval is fitted with a sum of signal and background distributions.
The background distribution is obtained using simulation of the main background process e + e − → π + π − π 0 π 0 .A possible inaccuracy in prediction of number background events is taken into account by introducing a scale factor α 4π .For energies below 1.594 GeV, the value of α 4π found in the fit is consistent with unity.At higher energies, there is significant background contribution from other processes, e.g., e + e − → π + π − π 0 π 0 π 0 .In this region α 4π is fixed to unity, and a linear function is added to describe contribution of other background processes.It is worth noting that in the energy region above 1.594GeV the shape of the M η distribution for e + e − → π + π − π 0 π 0 events is close to linear.
The signal distribution is described by a sum of three Gaussian functions with parameters determined from the fit to the simulated M η distribution obtained as a sum of the distributions of the four subprocesses (2)-( 5).To take into account a possible inaccuracy of the signal simulation, two parameters are introduced: mass shift ∆M η and ∆σ Mη .The latter is quadratically added to all Gaussian sigmas.These parameters are determined from the fit to the M η spectrum for events with E ≥ 1.544 GeV and found to be ∆M η = −0.9±1.0MeV/c 2 and ∆σ Mη = 12.0 ± 3 MeV/c 2 .
The number of signal and background events obtained from the fit to the M η spectrum in Fig. 2 is 4643 ± 126 and 8519 ± 185, respectively.The events with η meson belong to the process e + e − → π + π − π 0 η.The same fit is performed in the 13 energy intervals.The obtained numbers of e + e − → π + π − π 0 η events are listed in Table I.
To estimate the systematic uncertainty in the number of signal events due to the imperfect description of the shape of the background distribution, we perform the fit with an additional linear background below 1.594 GeV, and without the linear background but with free α 4π above.The uncertainty associated with the difference between data and simulation in the η-meson line shape is estimated by variation of the parameters ∆M η and ∆σ Mη within their errors.The obtained systematic uncertainties are listed in Table I.

V. SEPARATION OF INTERMEDIATE STATES
The distribution of the invariant mass of the system recoiling against the η meson M rec η for selected data events containing η meson is shown in Fig.  noted that the a 0 0 peak contains only one-third of a 0 ρ events.Unfortunately, our ηπ-mass resolution is not sufficiently good to unambiguously separate between a 0 ρ and nres contributions below 1.844 GeV, where the fraction of a 0 ρ events is lower.Therefore, in the further analysis we do not separate these two mechanisms.Instead, we analyze M rec η distributions at different c.m. energy intervals and measure the ωη and φη contributions, and the total contribution of the a 0 ρ and nres intermediate states.
The M rec η distributions are fitted with a sum of simulated distributions for the intermediate states ωη, φη, a 0 ρ, and nres.The ωη and φη distributions are described by triple-Gaussian functions, the parameters of which are determined from the fit to the simulated distributions.The data-simulation difference in the ω-meson line shape was studied in Ref. [6] using e + e − → ωπ 0 events.The difference is parametrized by the two parameters ∆M ω and ∆σ Mω .The latter parameter was found consistent with zero, while ∆M ω = 7.5 ± 1.9 MeV.
The same correction is applied to the φ-meson line shape.The histograms obtained using MC simulation are used to describe the a 0 ρ and nres distributions.The free fit parameters are the numbers of ωη and φη events, and the total number of the a 0 ρ and nres events.The ratio of a 0 ρ to nres events is fixed at the value measured in Ref. [4] and is allowed to vary within its uncertainty during the fit.The result of the fit to the M rec η for all selected data events containing η meson is represented in Fig. 3.
The e + e − → ωη process was studied previously in the SND experiment in Ref. [6], using the same data sample.The number of ωη events in the current analysis is about 15% larger than that in Ref. [6].One-third of this difference is due to the difference between the two analyses in the total number of selected events containing η meson.In the current analysis we use more accurate η-meson line shape, which take into account all four intermediate states contributing into the e + e − → π + π − π 0 η reaction, while in the previous the line shape was extracted from e + e − → ωη simulation.The remaining 10% is associated with the shape of the non-ωη distribution, which was described in Ref. [6] by a linear function.The a 0 ρ + nres model used in this analysis is certainly more realistic.In particular, it describes well the M rec η spectrum in the full range of M rec η variation (see Fig. 3).On the other hand, we observe that the measurement of the e + e − → ωη cross section is very sensitive to the shape of the non-ωη M rec η distribution.Therefore, an additional systematic uncertainty on the number of fitted ωη events is introduced to account this sensitivity, which is estimated to be 10%, the difference between results obtained with the linear background shape and the a 0 ρ + nres model.The same uncertainty is assigned to the number of fitted φη events.
The ωη and φη systematic uncertainties are translated to the uncertainty on the number of a 0 ρ + nres events as ∆N 2  ωη + ∆N 2 φη .

VI. DETECTION EFFICIENCY
The MC simulation used in this analysis takes into account radiative corrections, in particular photon emission from the initial state.Therefore, the detection efficiency (ε MC ) depends on the Born cross section for the simulated process and is determined using the iterative procedure described in Sec.II.In our case of the resonance energy region, the detection efficiency ε MC may strongly deviate, due to the initial state radiation, from the efficiency ε 0 determined at E γ = 0, where E γ is the energy of the photon emitted from the initial state.
The dependence of the efficiency on E γ may be parametrized as ε 0 (E)r(E, E γ ).The shape of the function r(E, E γ ) is defined by the cut on χ 2 3π2γ and is practically the same for all processes under study.The energy dependences of ε 0 for the four intermediate states contributing to the process e + e − → π + π − π 0 η are shown in Fig. 5.It is seen that in the energy range under study the efficiencies change from 10 to 14%.However, for the given energy the differences between the four efficiencies are within 5%.This difference is taken as an estimate of the additional model uncertainty on the detection efficiency in the measurements of the total e + e − → π + π − π 0 η cross section and the cross section for the sum of the a 0 ρ and nres intermediate states.For these measurements the simulation in the e + e − → ρ(1450)π, ρ(1450) → ρη model is used.The efficiency obtained in this model is corrected with a scale factor of ε0 (E)/ε nres 0 (E), where ε nres 0 is the efficiency at E γ = 0 obtained in the e + e − → ρ(1450)π model, and ε0 is the efficiency averaged over the four (two) intermediate states for the total (a 0 ρ + nres) cross section measurement with the weights equal to their expected fractions in the number of selected events.
The detection efficiencies obtained using MC simulation are shown in Fig. 6.The steep decrease of the e + e − → ωη efficiency is explained by a very small value of the Born cross section above 1.8 GeV.In this energy region practically all detected ωη events contain a photon emitted from the initial state, which distorts the event kinematics.tion efficiency ε and the efficiency determined using MC simulation

Imperfect simulation of detector response leads to a difference between the actual detec-
where κ i are the efficiency corrections for different effects.The main selection criterion for signal events χ 2 3π2γ < 30 is mostly affected by uncertainties in simulation.The quality of the simulation of the χ 2 3π2γ distribution is studied in Ref. [6] using e + e − → ωπ 0 → π + π − π 0 π 0 events, which has a large cross section in the energy region under study and the same number of the final particles as the process under study.The correction is found to be κ 1 = (2.5 ± 1.1)%.The correction on the χ 2 4π(γ) > 200 condition was also studied in Ref. [6] and was found to be consistent with zero with a systematic uncertainty of 4.6%.
The difference between data and simulation in photon conversion in detector material before the tracking system is studied using events of the process e + e − → γγ.The corresponding efficiency correction is κ 2 = (−1.35± 0.05)%.The largest part of the systematic uncertainties associated with data-MC simulation difference in track reconstruction cancels as a result of luminosity normalization.The difference in the track reconstruction for electrons and pions was studied in Ref. [18].The corresponding correction κ 3 = (−0.3± 0.2)%.
The total correction is κ = (0.9 ± 4.7)%.The corrected values of the detection efficiency for the process e + e − → π + π − π 0 η are listed in Table I.For the processes e + e − → ωη and φη, TABLE I: The weighted average energy (E ) for the c.m. energy interval specified, integrated luminosity (L), number of selected π + π − π 0 η data events (N ), the radiative correction factor (1+δ) and detection efficiency (ε) for the process e + e − → π + π − π 0 η.For N and ε the statistical and systematic errors are quoted.For 1 + δ the systematic errors are quoted.the systematic error is 4.8% and includes the uncertainty of the efficiency correction (4.7%), the uncertainty of the iterative procedure described in Sec.II.For e + e − → π + π − π 0 η and e + e − → a 0 ρ + nres, the 5% uncertainty related to inaccurately known process dynamics is added, and the total uncertainty is 6.9%.

VII. DETERMINATION OF THE BORN CROSS SECTIONS
The experimental values of the visible cross section for the processes under study are calculated as follows: where N i , L i , and ε i are the number of selected data events, integrated luminosity, and detection efficiency for the ith energy interval.The parameter B is equal to the branching fraction B(ω → π + π − π 0 ) = 0.892 ± 0.007 [14] for e + e − → ωη, B(φ → π + π − π 0 ) = 0.1524 ± 0.0033 for e + e − → φη, and unity for others.
The visible cross section (σ vis ) is related to the Born cross section (σ) by the following expression [15]: where x is the beam-energy fraction carried away by photons emitted from the initial state, the function F (x, E) describes the probability of radiation of photons with total energy xE/2, and The right side of Eq. ( 9) can be rewritten in the more conventional form where δ(E) is the radiative correction.
Experimental values of the Born cross section are obtained as follows.The energy dependence of the measured visible cross section is fitted with Eq. ( 9), in which the Born cross section is given by a theoretical model describing data well.The model parameters obtained in the fit are used to calculate δ(E i ) using Eq.(10).Here E i is the weighted average c.m. energy for ith energy interval.The values of the Born cross section are then obtained using the equation The two-resonance model is used to parametrize the Born cross sections where m V and Γ V are the mass and width of the resonance ) is the product of the branching fractions for the V decay to e + e − and the final state f , and P f (E) is the phase space factor.
For e + e − → ωη, the first term in Eq. ( 12) is associated with the ω(1420) resonance, the second is a sum of contributions of the ω(1650) and φ(1680) resonances, and P f (E) = q 3 ω (E), where q ω (E) is the ω momentum in the reaction e + e − → ωη.The phase between the first and second terms in Eq. ( 12) is chosen to be equal to π [6].In the fit to the e + e − → ωη cross-section data, the free parameters are B V ′ , B V ′′ , m V ′′ , and Γ V ′′ .The V ′ mass and width are fixed at the Particle Data Group (PDG) values for ω(1420) [14].The obtained fit parameters are listed in Table II.The fitted V ′′ mass is in agreement with the PDG mass of both ω(1650) and φ(1680) resonances [14], while the fitted width is smaller than the PDG estimate for the ω(1650) width, 315 ± 35 MeV [14], but agrees with the PDG value, 150 ± 50 MeV [14], for the φ(1680) resonance.
The φη data are well described by a single-resonance model with B V ′ = 0 and P f (E) = q 3 φ (E), where q φ (E) is the φ momentum in the reaction e + e − → φη.The obtained V ′′ mass and width listed in Table II are in agreement with the PDG values for the φ(1680).
The e + e − → a 0 ρ + nres → π + π − π 0 η and e + e − → π + π − π 0 η cross sections are described by the model (12) with seven free parameters (B and ϕ) and . This model has no physical sense, but describes data well, and therefore can be used to calculate radiative corrections.
The obtained values of the Born cross sections are listed in Table III and shown in Fig. 7 together with the fitted curves.The numerical values of radiative correction for the total e + e − → π + π − π 0 η are listed in Table I.The uncertainty on the radiative correction is estimated by varying the fitted parameters within their errors.For the cross sections we quote the statistical and systematic errors.The latter includes the uncertainties on the detection efficiency (systematic and model), number of selected events, luminosity, and radiative correction.The energy-independent (correlated) part of the uncertainty for the total e + e − → π + π − π 0 η cross section is 7%.
In Fig. 7 our results on the cross sections are compared with the previous measure-  ments.The obtained e + e − → ωη cross section agrees with the CMD-3 measurement [4].
Both the SND and CMD-3 results lie below the BABAR data [5].The SND and BABAR [7] measurements of the e + e − → φη cross sections are in reasonable agreement.The significant difference between the SND and CMD-3 measurements is observed for the e + e − → a 0 ρ+nres cross sections.The total e + e − → π + π − π 0 η cross section measured by SND is, in general, consistent with the CMD-3 result [4].The ∼ 15% difference in the cross section maximum is within the systematic uncertainties, which are 7% for SND and 11% for CMD-3.

VIII. SUMMARY
In the experiment with the SND detector at the VEPP-2000 e + e − collider in the energy range 1.34-2.00GeV the cross sections for the process e + e − → π + π − π 0 η and its subprocesses e + e − → ωη, e + e − → φη and e + e − → a 0 ρ + nres, where nres is the structureless π + π − π 0 η state, are measured.The cross sections have a peak near the energy ≃ 1650 MeV.The maximum value of the e + e − → π + π − π 0 η cross section is 5 nb, about 8% of the total hadron cross section at this energy.The obtained cross sections are rather consistent with the previous experiments and have comparable or better accuracy.The result on the e + e − → ωη cross section supersedes the previous SND measurement [6].The obtained e + e − → ωη and

FIG. 2 :
FIG.2:The M η spectrum for selected data events (points with error bars).The solid histogram is the result of the fit to the data spectrum with a sum of signal and background distributions.The fitted background contribution is shown by the dashed histogram.

FIG. 3 :
FIG. 3: The M rec η distribution for selected e + e − → π + π − π 0 η data events (points with error bars).The solid histogram represents the result of the fit to the M rec η distribution with the sum of simulated distributions for the processes (2)-(5).The dashed histogram shows the fitted distribution for the sum of the processes (4) and (5).

3 .
The number of events in each M rec η bin is determined from the fit to the M η distribution as described in Sec.IV.The peaks in the M rec η spectrum at the ω and φ masses correspond to ωη and φη events.The a 0 ρ and nres intermediate states have wide M rec η distributions with similar shapes.The contribution of the a 0 ρ mechanism can be observed in the ηπ invariant mass spectrum.The ηπ 0 mass spectrum for data events containing η meson from the energy region E ≥ 1.844 GeV, where the a 0 ρ mechanism dominates[4], is shown in Fig.4.The spectrum is fitted by the sum of the a 0 ρ, nres, and φη distribution.The number of φη events and the absence of ωη events are determined from the fit to the M rec η distribution for the energy region E ≥ 1.844 GeV.The clear signal of the a 0 0 meson is seen in Fig.4.It should be

FIG. 4 :
FIG. 4: The distribution of ηπ 0 invariant mass for data events (points with error bars) with energy E ≥ 1.844 GeV.The solid histogram represents the result of the fit with the sum of simulated distributions for the a 0 ρ, φη, and nres intermediate states.The dashed histogram shows the fitted sum of the φη and nres distributions.

FIG. 6 :
FIG.6:The energy dependence of the detection efficiency, corrected as described in text, for the process e + e − → 3πη and the intermediate states ωη, φη, and a 0 ρ + nres.

TABLE II :
The fit parameters for the processes e + e − → ωη and e + e − → φη, and χ 2 /ν, where ν is the number degrees of freedom.