An extension of BDK event generator

. An extension of the well known ee → 4 l event generator of F. A. Berends, P. H. Daverveldt, and R. Kleiss is presented. This extension accounts for production of narrow resonances, as well as the initial state radiation. In addition, it allows one to simulate the production of l + l − π + π − with the lepton pair and the pion pair in 1 −− states.


Introduction
BDK is the popular nickname for ee → 4l generator of F. A. Berends, P. H. Daverveldt, and R. Kleiss [1]. The authors named their generator DIAG36 according to number of QED diagrams accounted in e + e − → e + e − + e + e − simulation. Our work on its extension is inspired by needs of some analyses in KEDR [2] and Belle [3].

BDK generator
The generator is suitable in the energy region of C− and B− factories where electroweak corrections can be neglected. Electroweak processes are accounted in NEXTCALIBUR generator by F.A. Berends et al. [4].
DIAG36 includes four subgenerators for the four types of diagrams: multiperipheral, bremsstrahlung, annihilation and conversion. Generated events are reweighted according to the square of the exact lowest-order QED matrix element so that the interference effects are accounted.
The best experimental tests of generator for e + e − → e + e − +e + e − and e + e − → e + e − +µ + µ − has been performed by L3 experiment at LEP [5], the accuracy of about 2% has been reached.
3 Identity effect in e + e − → e + e − + e + e − process The BDK generator takes into account all diagrams corresponding to permutation of identical particle momenta. It is known that in the total cross section the identity effect is negative and is suppressed by the factor ln(s/m 2 e ) [6]. However, BDK generator predicts that identity becomes important when the cut on the invariant mass is applied. This statement is illustrated by Figure 1. Worth note that identity effect is not considered in some popular two-photon event generators.

Alternative matrix element expression
The matrix element expressions in the BDK extension have been obtained by P.A. Krachkov and A.I. Milstein. Using these expressions one can • check the identity effect accounting • simulate e + e − → l + l − π + π − for the pairs in 1 − − states • perform calculations, using C++, with the multiple precision, which is important for simulation of e + e − → e + e − + e + e − at high energies.
The two sets of matrix elements for leptonic processes give identical numerical results despite to very different approaches to calculations.
Among the e + e − → l + l − π + π − diagrams only those are considered which are not suppressed with the pion QED form-factor at the collision energy. They can not be treated accurately, however, they can be neglected at collision energy above 3 GeV. The code for the form-factor calculation was taken from the MCGPJ generator [7].

Accounting of narrow resonances
The narrow resonances φ, J/ψ, ψ(2S ), Υ(1S ), Υ(2S ) and Υ(3S ) were accounted via the vacuum polarization operator. In the original BDK codes and in the alternative matrix element calculations the photon q 2 was replaced with (1− Π R ) q 2 , where the resonance contributions in the vacuum polarization is given by where Γ (0) ee is the electron partial width, Γ 0 is the total resonance width and M 0 is its mass. The naked parameters employed differ from those in PDG tables. For discussion see [8] The cross section of e + e − → µ + µ − π + π − process as a function of the muon and pion pairs invariant masses are shown in Figure 2. The resonance-continuum interference is clearly seen. The total cross section of 2µ 2π and 4µ productions at 10 GeV accounting for ISR are about 550 and 380 fb, respectively. The characteristic Υ(1S ) cross section is 0.15 fb.

Accounting of initial state radiation
For the two-photon processes, the initial state radiation (ISR) is important if the multiperipheral graphs do not dominate. For the process e + e − → µ + µ − µ + µ − at 10 GeV the ISR correction to the total cross section is about 10%. The ISR is accounted using the result of E.A. Kuraev and V.S. Fadin [9] for the cross section σ(s), where s is the c.m. energy squared, F (s, x) is the radiative correction function, and σ 0 (s) is the cross section without account of radiation. To calculate it, the preparatory run of the generator is required. The cross sections of e + e − → µ + µ − µ + µ − reaches the maximum of about 950 fb at the energy of 2.5 GeV approximately. For the e + e − → µ + µ − π + π − processes correspondent numbers are 2500 fb and 2 GeV. Eq. (2) is used to generate the amount of radiated energy in the event and to determine the collision energy at which the event should be generated. The produced particles are boosted to the laboratory frame under the assumption that only one photon is radiated along or opposite to the initial electron.  Accounting of ISR becomes essential when the recoil mass method is employed in the data analysis. In this method it is assumed that the energy and the momentum of the produced system equal to those of colliding beams. The machine energy spread and ISR violate this assumption, however, the recoil mass analysis is useful in many cases. Figure 3 illustrates the applications of this method to studying of e + e − → Υ(1S )π + π − process with the dimuon decay of Υ(1S ) at Belle conditions. The two-photon production e + e − → π + π − µ + µ − gives the peaking background to the single photon process. Because of the large mass of Υ(1S ) the shape of recoil mass signal is rather similar to that of the e + e − → Υ(1S )(γ) cross section.

Comparison with data
The prediction of the BDK extension are compared with the experimental data at the c.m. energy of 10.52 GeV in Figure 4. The distribution in the pion invariant mass were obtained for the muon pair mass in the range 8.7-9.4 GeV.  Figure 4. The recoil mass for pion pairs (left), their invariant mass (right). Points with error bar represent preliminary Belle data [3], red lines show simulated contributions of e + e − → 2π2µ, green and blue lines are for contributions of misidentified muon and electron pairs, respectively. The black lines show the sum. The peaks at 9.6 and 10.2 GeV are reflections of the transition between narrow Υ states, which was not simulated.