Event Generators for Multi Hadronic Production in e+e− Colli- sions

An event generator for multi hadron production is presented for measuring the R value in theτ-charm energy region withe+e− collisions. The initial state radiation effects are considered up to second order accuracy, and the radi ative correction factor is calculated with hadronic Born cross sections. The establis hed exclusive processes are generated according to their measured cross sections, whil e the missing processes are generated using the LUND Area Law model, and its parameters a r tuned with data collected at √ s = 3.08 GeV. The optimized values are validated with data in the ra nge √ s = 2.2324∼ 3.671 GeV. These optimized parameters are universally valid f or event generation below the DD̄ threshold.


Introduction
The total cross section for multi hadron production in positron-electron (e + e − ) annihilation is one of the most fundamental observables in particle physics. A precise measurement of the hadronic cross section allows us to determine the hadronic contributions to the running of the quantum electrodynamic (QED) fine structure constant α, electroweak parameters, and the strong coupling α s . The R value, defined as the ratio of the total hadronic cross section to that of e + e − → µ + µ − at Born level, have been measured by many collaborations in e + e − scan experiments, over the center-of-mass energy from the two pion mass threshold (M 2π ) to the Z peak [1]. In the tau-charm energy region, the R values measured at BESII [2] were used in the evaluation of the hadronic contribution from the five quark loops at the energy of Z peak, ∆α (5) had (M 2 Z ), with an improved precision by a factor of 2 [3]. A large number of exclusive processes have been measured over the range from M 2π to 5 GeV [4], but most cross sections have large uncertainties. To improve these measurements, a hadronic event generator is needed for us to get better understanding of background events from e + e − → hadrons.
Especially, a precise R-value measurement requires excellent control of radiative correction (RC) and vacuum polarization (VP) in the Monte Carlo (MC) program. We design an event generator for measuring R values and exclusive decays in e + e − collisions. The generator is constructed in the framework of BesEvtGen [5], incorporating both the RC and VP effects. We also present details of the parameter optimization of the Lund Area Law (LUARLW ) model [6] with data, and validations with various distributions within the energy range √ s = 2.2324 ∼ 3.671 GeV.

Framework of event generator
The generator is constructed as a model of the BesEvtGen package. It provides the 4-momentum of each final state particle for detector simulation, and provides the ISR correction factor and VP factors for users to undress the observed cross section. The basic idea of this generator is to decompose the total hadronic cross section into the measured exclusive processes and remaining unknown processes. The latter are generated with the LUARLW model. Figure 1. Feynman diagrams for the process (a) e + e − → X i , and ISR process (b) e + e − → γ ISR X i .

Initial state radiative correction
In an e + e − energy scan experiment, we consider a measurement of the Born cross section (σ 0 ) for a process e + e − → X i , as shown in Fig. 1 (a), where X i denotes the hadron states of i-th process. Due to ISR, the observed cross section (σ) is actually for the process e + e − → γ ISR X i , as shown in Fig. 1 (b). The observed cross section is related to the Born cross section by the quasi-real electron method [7]: where m is the invariant mass of the final states; π(m) is the vacuum polarization function, which will be discussed later; s is the e + e − center-of-mass energy squared; x ≡ 2E * γ / √ s = 1 − m 2 /s, and E * γ is the total energy carried by ISR photons in the e + e − center-of-mass frame; M th is the mass threshold of a given process. W(s, x) is a radiative function, we use the result of QED calculation up to order α 2 [8].
At the leading order of QED calculation, the ISR photon is characterized by soft energy and beam collinear distribution. A more general result is obtained by the method of Bonneau and Martin [11] up to m 2 e /s terms, and the angular distributions is calculated by with (1−2x) sin 2 θ−x 2 cos 4 θ x 2 −2x+2 where E is the beam energy in the center of mass system of the electron and positron.

Vacuum polarization
The vacuum polarization (VP) has been calculated by many groups and is available in the literature. Comparisons between them are given in Ref. [12]. There are notable differences below 1.6 GeV, and above 2.0 GeV; visible differences appear when approaching the charmonium resonances. We use the results from the Fred Jegerlehner group [13]. It provides leptonic and hadronic VPs both in the spaceand time-like region. For the leptonic VP the complete one-and two-loop results and the known high-energy approximation for the three-loop corrections are included. The hadronic contributions are given in tabulated form in the subroutine HADR5N [14].
The narrow vector resonances, such as ψ(3770), ψ(2S ), J/ψ, ρ(1700), and ω(1420), are also included in the calculation for the ISR correction factor. The cross sections for these narrow resonances are represented with the Breit-Wigner function where M, γ, and γ ee are the mass, total width and partial decay width to e + e − final state, respectively. The distribution of cross section versus center-of-mass energy is described by an empirical function, which is parameterized with a multi-Gaussian function. Its parameters are determined by fitting the cross section mode by mode. These empirical functions are used in the generator for the calculation of the ISR correction factor and event type sampling.
The angular distribution for ISR photons is implemented according to Eq. (2). However, angular distributions are implemented only for two-body decays, namely, 1 − cos 2 θ for PP (where P is a pseudoscalar meson) modes, and 1 + α cos 2 θ for the PV (α = 1) and BB modes, where V is a vector meson, and B is a baryon. The angular distribution parameter α for the BB mode is taken as the quark model prediction [27]. The phase space model is used for multi-body decays.

LUND Area Law model
The hadronic events produced in the e + e − annihilation are evolved as follows. As the first step, a quark-antiquark (qq) pair is produced from a virtual photon, coupled to the e + e − pair. Then the qq branching proceeds via emitting gluons, and further develops into hadrons. In the high energy region, the cluster model (e.g. HERWIG [28]) and LUND string model (e.g. JETSET/PYTHIA [29]) are available and precise enough to describe the hadronic fragmentation with parameters optimized at boson Z peak. However, in the intermediate and low energy region, parameters need to be optimized or a new model is desirable to describe the light quark fragmentation.
In the tau-charm energy region, the LUARLW model [6] has been proposed to estimate the multiplicity distribution for primary hadrons produced from the string fragmentation. The probability distribution reads: where c 0 , c 1 , c 2 , α, β and γ are parameters to be tuned with data. An interface to access the LUARLW model is designed in the BesEvtGen [5] framework, and is only used to generate the primary hadrons. The further decays into light hadrons are realized with BesEvtGen [5].

Monte Carlo algorithm
The event sampling proceeds via two steps. Firstly, the mass of the hadron system, M hadrons , is sampled according to the distribution of the observed cross section, i.e. dσ(s)/dm, for the process e + e − → γ IS R X i according to Eq. (1). For simplicity, the ISR energy, √ s − M hadrons , is imposed on a single photon. The second step is to sample the event type topology according to the ratios of individual cross sections at the energy point M hadrons .

Sampling of M hadrons
To sample the M hadrons , we split the region M th ∼ √ s into a few hundred intervals. The cumulative cross section up to the i-th interval, m i , iŝ The M hadrons is sampled according to theσ(m i ) distribution with the discrete MC sampling technique.

Sampling of event type
Using the discrete MC sampling technique, the final states for exclusive modes are sampled according to the ratios of their cross sections (σ m ) to the total cross section (σ tot ), i.e., where m is an index for exclusive precess, and events for the remainder part, 1 − m c m , are generated with the LUARLW model.

Strategy to optimize the LUARLW parameters
The LUARLW model parameters are optimized with the parameterized response function method. The optimal values are obtained by simultaneously fitting this function to data distributions. The idea for this method is borrowed from that implemented in the event generator tuning tool Professor and Rivet [30] system, which was introduced by TASSO, and later used by ALEPH, DELPHI [31][32][33][34][35][36], and recently by the LHC [30]. This method has the advantage of eliminating the problem from the socalled manual and brute-force tunings, such as the slow tuning procedure and the sub-optimal results. An ensemble of MC samples was produced within the framework of the BesEvtGen [5] event generator, and then it is subject to detector simulation with BOSS software [40]. 91 independent MC samples were prepared, each one generated with a different set of LUARLW parameters, which were randomly chosen in the parameter space around a given central point p 0 . All MC samples were produced with equal statistics, and were large enough so that the overall statistical uncertainties are negligible.
By including the correlations among the model parameters, the dependence of physical observable is expanded up to the quadratic term as done in Ref. [37], and the response function reads where n is the number of parameters to be fitted, and MC(p 0 + δp, x) denotes the distribution of physical observable x predicted for a given set of parameter values p 0 + δp, where p 0 is the central value and δp i is the deviation of the i-th parameter. The quadratic term in the expansion accounts for the possible correlations between the model parameters. The number of coefficients a (0,1,2) , L, in the expansion is calculated with L = 1 + n + n(n + 1)/2, and the coefficients are determined by fitting Eq. (5) to the L reference simulation distributions. This fit is equivalent to solving a system of linear equations of Eq. (5). Then the optimal values of the parameters p i , their errors σ i , and their correlation coefficients ρ i j will be determined with a standard χ 2 fit to data using package MINUIT [38]. The fit is done simultaneously for all distributions and for all bins.
To minimize statistical uncertainties, the model parameters should be fitted to the distributions that show strong dependence on the parameters under consideration and least dependence on the others. For each distribution, a quality to measure the sensitivity to the model i-th parameter is calculated, i.e.
where δMC(x) is the change of the distribution MC(x) when the model parameter p i is changed by δp i from its central value. Sensitivity values for charged track distributions and event shapes vary within the range from -0.3 to 0.3, but the polar angle and azimuthal distributions for charged tracks are not sensitive to the change of model parameters. This is because the inclusive charged tracks are distributed isotropically over the whole phase space. Taking the sensitivity into consideration, only 12 observable distributions are kept for the model parameter fit. They are the number of photons (N γ ), the number of charged tracks (N track ), momentum of tracks (P track ),  [29], where W is the total reconstructed energy of an event, and P ⊥ is the transverse momentum. We have 12 parameters to be optimized. According to Eq. (6), there are 91 coefficients, a (0,1,2) in Eq. (5) to be determined. Hence we need at least 91 MC samples to determine these coefficients. These were prepared with 0.5 million events for each sample. Then the dependence of response function on model parameters is established, and this analytical expression is used to simultaneously fit to the data distributions after QED background events are subtracted. In the optimization procedure, the χ 2 function is defined over each bin, ie., χ 2 → χ 2 /N, where χ 2 values are calculated over nonempty N Table 1. Optimized parameters at √ s =3.08 GeV. The statistical errors are negligible. (2S +1) P J denotes a meson has spin S , orbital angular momentum (L) and total spin J.

Parameters
Tuned Description PARJ (1) 0.118 Suppression of diquark-antidiquark pair production PARJ (2) 0.670 Suppression of s quark pair production PARJ (11) 0.868 Probability that a light meson has spin 1 PARJ (12) 0.644 Probability that a strange meson has spin 1 PARJ (14) 0.188 Probability for a 1 P 1 meson production PARJ (15) 0.232 Probability for a 3 P 0 meson production PARJ (16) 0.518 Probability for a 3 P 1 meson production PARJ (17) 0.320 Probability for a 3 P 2 meson production PARJ (21) 0.201 Width of Gaussian for transverse momentum RALPA(67) 0.191 LUARLW model parameter RALPA (16) 1.000 LUARLW model parameter RALPA(17) -0.537 LUARLW model parameter bins. To consider the requirement of fit goodness on the multiplicity of charged tracks, this distribution is weighted with a factor of 10, while other distributions are weighted with a unitary factor. This weighted factor is chosen by requiring that the fit quality of all distributions are satisfactory.

Event selection and fit results
We use the data taken at √ s =3.08 GeV to optimize the parameters. To validate the parameters, we compare the MC distributions to the data distribution within the energy region 2.0 -4.26 GeV. The QED backgrounds, e.g. e + e − → e + e − , γγ, γ * γ * , µ + µ − , and τ + τ − are subtracted using MC samples, and they are normalized according to their cross sections to the luminosity of data sets. The event selection criteria for light hadrons are similar to those applied to the R value measurements [2,39].
The selected candidates are characterized by the distributions of charged track multiplicity (N track ), track energy (E track ) and momentum (p track ), polar angle (cos θ), azimuthal angle (φ), rapidity, peseudorapidity, and a set of event shapes. These distributions are normalized to one and the errors are scaled for all bins.
To consider the possible correlations between these observable quantities, different observable combinations were tried. In each combination, track observables, N γ , N track , E track , x f and x ⊥ , must be included, while the p track distribution or event shapes are partly included in the simultaneous fit. Generally speaking, the more observable distributions are involved in the fit, the worse fit quality one gets. To validate the resulted parameters, they are reused to generate MC samples, and compared to data.

Validation of tuned parameters
To select the most optimal values, we compare the data to the MC distributions, which are generated with optimized parameters for all sets. We require that the parameters can produce MC distributions having the best fit goodness quality χ 2 /N, where N is the total number of bins for calculating the χ 2 values. The optimal values are given in Table 1. To note that these values are only responsible for unknown processes other than the exclusive modes.
To validate this set of parameters for the MC generation below the DD threshold, we compare the multiplicity distributions for charged tracks at 14 energy points from √ s = 2.2324 to 3.671 GeV, as shown in Fig. 2. When extending these parameters from 3.65 GeV to low energy points, the agreement between the data and MC gets better. This is due to the fact that the total cross section equals the sum of the exclusive ones when approaching the energy 2.0 GeV.

Discussion and summary
To summarize, we have developed an event generator for R measurement at energy scan experiments, incorporating the initial state radiation effects up to the second order correction. The ISR correction factor is calculated using the measured Born cross sections. The established exclusive processes are generated according to their measured cross sections, while missing processes are generated using the LUARLW model, with tuned parameters at 3.08 GeV. To validate the optimized parameters, various MC distributions are compared to the data distributions with data from energy √ s = 2.2324 to 3.671 GeV. We conclude that the optimized parameters are valid for MC generation below the DD threshold. Above the DD threshold, the parameters should be optimized with the charm meson decays.