Precision calculation of processes used for luminosity measurement at the ZEUS experiment

The process pe ->pe+gamma with the photon emitted along the electron beam axis is used for luminosity measurement at HERA. In this paper the process is calculated including one-loop QED radiative corrections. In the ZEUS experiment, both the electron and the photon can be detected. Therefore both photon and electron spectra with and without the gamma-e coincidence are analyzed. We also calculate the process pe ->pel+l- which contributes to the background in the electron tagger.


Introduction
The luminosity in the ZEUS experiment at the HERA electron-positron collider is measured using the inclusive bremsstrahlung process pe ± → pe ± γ at very small photon angles (Θ γ ∼ 0.1 mrad). A photon spectrometer system is located on the electron beam axis at a distance of about 100 m from the interaction point [1]. Additional detectors are used to measure forward-going electrons 3 with energies between 4−9 GeV. The electron rate is used to calibrate the acceptance of the photon spectrometer. This study is motivated by augmented analysis precision which makes it necessary to go beyond the classical Bethe-Heitler formula [2] used for simulation in ZEUS so far. These simulations assume scattering on a point-like, spin-less proton. No proton recoil and no higher-order effects were taken into account.
In this paper we make an analysis of the processes affecting both the photon detector and the electron tagger in ZEUS. Two experimental configurations are considered, with and without coincidence of electron and photon. We consider the process pe ± → pe ± γ with the photon emitted along the electron-beam axis. The one-loop QED radiative correction ( Fig. 1) to the leptonic current is calculated. The reaction in which a lepton pair is produced, pe ± → pe ± l + l − (Fig. 2) with at least one detectable electron, is also considered. It does not affect the photon spectrometer but is required for the proper simulation of the electron tagger signals. Finally we discuss the systematic error of the luminosity measurement originating from the uncertainty in the cross section calculation.
In previous theoretical estimates for these processes [3,4] good agreement with the classical formula was observed. The contribution of higher-order effects appeared within the 1%-error region. However, these calculations were made for different experimental configurations (coincidence between electron and photon) and for a narrow photon energy range (8 < E γ < 14 GeV). No generator program was published that can be used for detector simulations.
This study is restricted to the QED leptonic current corrections only. Weak effects are suppressed by the Z-boson mass in the considered momentum transfer region (Q 2 < 10 −5 GeV 2 ). The two-photon exchange vanishes in the forward region according to the existing estimate [5]. The photons radiated by the initial or   Figure 1: The lowest order diagrams for pe ± → pe ± γ (A) and pe ± → pe ± γγ (B) processes and QED loop diagrams considered (C).
final state proton are directed along the proton beam axis and escape the calorimeter. The other processes which are discussed in [6,7] may also be neglected.
All results in this paper are compared to the classical Bethe-Heitler formula [2]: Here E e (E ′ e ) is the energy of the initial(final) electron, α is the fine structure constant and r e = α/m e is the classical electron radius. The approximation for the angular distribution [8] is also checked: The paper is organized as follows: The kinematics are discussed in Sec. 2 where it is shown that kinematic variables may differ by up to a factor 10 30 at opposite edges of the phase space. The matrix element calculation is discussed in Sec. 3. Using those results a generator program has been created [9] which is described in Sec. 4. This section also describes the precision-saving algorithm for the phase space generation procedure. This algorithm is necessary since the computation is complicated by large numerical cancellations at small scattering angles. The results are shown in Sec. 5. The paper closes with a conclusion.

Kinematics
We denote the initial and final particle momenta as follows: The energies of the final(initial) electron and the photon are denoted by E ′ e (E e ) and E γ . We consider the phase space region with E γ > E cut γ ∼ 5 GeV and the polar angle Θ γ < 0.001. The electrons can be detected  if 4 < E ′ e < 9 GeV. The initial beam energies at HERA are [1]: E p = 460, 575 and 920 GeV, E e = 27.6 GeV.
The momentum transferred by the proton is denoted by q: The variation of Q 2 is very large: The process of the lepton pair creation is also considered (l = e, µ): It affects the electron tagger only and is expected to be a background for the Bethe-Heitler process. The range of Q 2 for this process is:  Events(scaled) Figure 3: The Q 2 -distribution for the processes; pe ± → pe ± γ (solid line, no cuts on E ′ e ), pe ± → pe ± e + e − (dotted line, 4 < E ′ e < 9 GeV, scaled by factor ×10 2 ), pe ± → pe ± µ + µ − (dashed-dotted line, 4 < E ′ e < 9 GeV, scaled by factor ×10 6 ).
The actual distributions of the Q 2 values for the considered processes are shown in Fig. 3. The effective Q 2 region is much smaller than the allowed region due to the matrix element. In this paper, values of Q 2 > 10 −3 GeV 2 are in fact neglected.

Matrix element calculation
For the one-photon exchange process the cross section is written in terms of the leptonic and hadronic tensors: The leptonic tensor is composed of the electromagnetic current vectors as: The plane wave approximation is used in spite of the macroscopic interaction distance at the lowest edge of the Q 2 region (r ∼ 0.1 mm at Q 2 ∼ 10 −24 GeV 2 ). The hadronic tensor for unpolarized protons is written as: Here G E and G M are the Sachs form factors and τ = Q 2 /4m 2 p . We use the dipole parametrization: Here Λ 2 D = 0.71 GeV 2 and µ p is the magnetic moment of proton. The precision of the dipole formula is very poor for Q 2 values greater than 0.1 GeV 2 [10]. We use it here only to check whether the proton size effects are visible in the considered processes. As can be seen from Fig. 3, the effective region of the transferred momentum is Q 2 < 10 −5 GeV 2 for the photon-and e + e − -pair production processes. The dipole form factor (7) differs from unity by 0.003% in that region. On the other hand, modern form factor formulae [10,11] were mostly fitted with Q 2 > 0.1GeV 2 experimental data. They have a similar (1 − G E )value at very low Q 2 . Hence we may either use the simple dipole formula (7) or neglect proton structure effects accepting the value of (1 − G E (Q 2 max )) ∼ 0.003% as the systematic error of the calculation. The error of our estimate for the pe → peµ + µ − process (10 −12 < Q 2 < 1 GeV 2 ) is therefore approximately 20% due to the form factor uncertainty.
The symbolic algebra program ALHEP [12] is used for analytical computations. It creates diagrams using the Standard Model Lagrangian, calculates matrix elements and reduces tensor integrals over the virtual particles' phase space to scalar ones. Finally the matrix element is written in terms of the scalar products of the initial and final momentum vectors. The infrared divergences are regularized using a fixed photon mass parameter. The contributions of the soft radiative corrections are integrated analytically over the photon energy in the region m IR γ < E γ < E SH γ in the laboratory frame. The basics are well described in Ref. [13]. The E SH γ energy cutoff separates the soft and hard bremsstrahlung contributions. We vary it in the range of 10 −8 − 10 −4 GeV to control the generator consistency. The LoopTools [14] package is used for the numerical estimation of loop integrals. It was recompiled with real*16 precision to avoid internal instability at very small Q 2 .
The matrix elements suffer from huge numerical cancellations in the considered phase space region. The first problem appears when calculating momentum couplings by contracting the corresponding vectors. It is solved as described in Sec. 4. Another problem has to do with the cancellation of peaking terms in the matrix element. Let us consider typical values of the scalar products for the Born process: The major cancellation appears due to p i ·k j couplings. The p 2 ·k j and p 1 ·k values are easily expressed via the others. But the two remaining couplings ( p 1 ·k 1 and p 1 ·k 2 ) are independent and may cancel numerically.
We use the substitution: Typical values for δ 1 are: Thus calculating δ 1 precisely as δ 1 = ( k·k 1 k 1α − k·k 2 k 2α )·p α 1 one can avoid one of the two p 1 ·k i couplings. In this way the Born differential cross section can be written as: Here the remaining k 1 ·p 1 -terms appear coupled to the tiny Q 2 factor. The numerical behavior of this formula is stable enough to be used with fast 8-byte variables in the analysis. Alternatively this matrix element can be written in a more compact but still numerically stable way in terms of Levi-Civita tensor couplings [3].
For the higher-order processes there is no simple substitution available, and the straightforward expansion in terms of one major term is not possible. At least 4 independent couplings may occur as major terms in various phase space regions. This problem can be partially solved using the Levi-Civita tensor couplings method [3]. The other possible solution is to split the phase space and to adjust the matrix element separately for every part. This method leads to enormous growth of the code and problems during the debugging step. In this analysis we solve the precision problem using high-precision floating point numbers [15]. The loop and the hard bremsstrahlung corrections (3) and the lepton pair creation process (5) are calculated with 16-byte variables. They allow to use the same matrix element form in the entire phase space region.
For cross-check purposes several matrix element representations (using various sets of independent variables) are compared. The results are similar. However, numerical precision for different representations is quite different. Expressions which depend on the transferred momentum q·k i -couplings are much more efficient. Some representations require at least 32-byte variables for proper evaluation. The splitting parameter and renormalization constants are also varied to check the consistency of the generator.
The higher-order matrix elements are not compact enough to be given in this paper. One can find all the formulae (converted to Mathematica format) together with ALHEP-scripts used for their derivation in Ref. [16].

Generator
The typical adaptive Monte-Carlo generator is based on iterative phase space splitting according to the distribution of the integrated function (or its derivative). At the first (integration) step, the appropriate splitting grid is created and at the second (generation) step, the events are sampled. However, the narrow peaks of the function may lead to significant integration errors. Since only a finite number of events is generated in every cell, the peak may remain undetected. The integration in the forward region is complicated by the fact that the effective region of the phase space is negligible as compared to the total range of the variables. For example, the FOAM [17] generator (with default parameters) being efficient for other purposes failed to integrate the Born cross section (9) in the whole phase space region (4).
Cuts may also lead to integration errors. In our case several cuts must be applied at the integration step: • Minimum energy of the detected photon, E min γ ∼ 5 GeV; • Minimum energy of the undetected photon, E SH γ ∼ 10 −8 · · · 10 −5 GeV (soft and hard bremsstrahlung separator); • Energy cut on the detected electron in the l + l − -production process: 4 < E e < 9 GeV. The cross section with and without this cut is about 0.1 and 4.0 mb respectively. This cut needs to be applied at the integration level. Otherwise the event generation is slowed by a factor 40.
Therefore a new adaptive Monte-Carlo integration program was created. It uses the photon energy as integration variable, thus avoiding E γ cuts. The main idea is to carefully tune the sub-space splitting algorithm to minimize the integration time and the amount of memory used. The distribution of the derivative is examined to detect narrow peaks. A separate splitting algorithm is used if a function cut is detected inside the cell. Generated function values are re-used for further sub-cell integration wherever possible. The maximum number of steps of the phase space splitting is about 110 for 2 → 3 processes and 130 for 2 → 4 processes. The statistical error of the integration is distributed according to the function value: δI i ∝ I i . For cross-check purposes, the Born process is partially integrated with the FOAM [17] generator. The FOAM results are stable in the sub region Q 2 > 10 −9 GeV 2 and are in good agreement with our results.
Another problem is to compose the appropriate phase space reconstruction procedure. This procedure is used to map a hypercube point r i ∈ (0, 1) N produced by the generator to the phase space of final particles. The standard method of factorizing the phase space to 2-body decay sub spaces [18] is used: dΓ n (p n ) = (2π) −1 dp 2 n−1 dΓ 2 (p n → k n + p n−1 ) dΓ n−1 (p n−1 ).
However, the vectors of the 2-particle sub space are not calculated in the rest frame of the total momentum p n since Lorentz-transformations in the forward region lead to about 6-digit precision loss. The procedure is described in the following: We generate the recoiling proton first: Here φ z is the free axial rotation angle for the event, p n−1 and Γ n−1 are the momentum and the phase space of the remaining particles (e.g. p n−1 = k + k 2 for the Born process), and the term Q 4 cancels one in the denominator of the matrix element. The limits for Q 2 are defined in Sec. 2 for each process. The limits on p 2 n−1 are: For the pe → peγ process the most convenient order of the remaining integration is: Here φ k+ k2 is the angle of the plane spanned by the electron and the photon with respect to the direction of their sum vector k + k 2 in the lab. frame. The integration limits on the energy are defined by the general formulae: For the Born process here: p = p n−1 = k + k 2 , k i = k γ , m i = 0, (p − k i ) 2 = m 2 e and the E cut max is unused. Finally, the limits on the photon energy are: For the two photon bremsstrahlung process p n−1 = k + k 2 + k ′ and the phase space integral is given by: Here the limits on E γ are defined by Eq. (14) with p n−1 = k + k 2 and the additional condition Similar formulae are used for the lepton pair production process (p n−1 = k 2 + k 3 + k 4 ): Here the limits on E ′ e are given by Eq. (13) with p = p n−1 , m i = m e , (p − k i ) 2 = (k 3 + k 4 ) 2 and the cuts on the final electron are the E cut min/max parameters. The integration limits for the second lepton energy E l − are calculated using the same formula with p = k 3 + k 4 , m i = m l , (p − k i ) 2 = m 2 l and no additional cuts for the µ + µ − -production process. For the case l = e we set E max l − = E ′ e in Eq. (13). Alternatively, one may apply the detector cuts on the second electron together with the factor 1/2 to avoid double-counting of identical particles.
Calculating the scalar products by contracting previously reconstructed momenta also leads to large numerical errors. Therefore we define products in the generator using the original r i -values.
The new code saves up to 30 decimal digits of precision when compared to the typical algorithm that uses momenta in the two-body rest frames [18]. The Born process converges with the standard 8-byte floating point variables. Sixteen-byte floating point variables (QD package [15]) are used for the higher-order process calculation.
A different set of integration variables is implemented to cross check the generator. The recoil proton is still reconstructed according to Eqs. (10)(11) and the sequential approach from Ref. [18] is used for the remaining particles. Here the energy cuts are applied on the integrating function and the center-of-mass frames are used for momentum reconstruction. The result of the integration agrees with the previous one but the generator requires 32-byte variables (and huge computer resources) for numerical stability.
Due to the use of high-precision floating point numbers, the performance of the generator must be carefully tuned. Both integration and generation steps require sufficient computer resources: up to 2Gb of memory and about a day of the CPU time. To optimize the generator performance we make parallel runs wherever possible. Thus the integration of the sub processes is performed separately and events are generated in parallel runs using independent sets of random numbers. We also integrate the loopand factorisable soft bremsstrahlung corrections using the previously created Born-process pattern. At

Event Generation
Performance: ~10 6 evts/hr Parallel runs are available 16-byte variables Parallel runs Figure 4: The generator program structure. The cross section values correspond to E cut γ = 10 GeV for the photoproduction processes and 4 < Ee < 9 GeV for the lepton pair production process. every point in phase space the Born matrix element is within 10 − 30% of the corrected one, |M Born | 2 ≈ |M Born | 2 + 2ℜ(M * Born M loop ) + |M R sof t | 2 . Therefore we use the Born process integration grid and simply update the integral value in every sub cell (using the distributed computation again). The scheme of optimal usage of the generator is outlined in Fig. 4.

Results
We consider the process pe ± → pe ± γ with at least one high energy photon in the final state (E > E cut γ ). The value of E cut γ is set to 10 GeV. We also simulate the energy acceptance of the photon spectrometer [1] using an approximate formula. For the simulation of the lepton-photon coincidence, an additional cut on the lepton 4 < E e < 9 GeV is used. The background lepton pair creation process (5) is also considered. When the two identical particles (e or γ) appear in the detectable region only the hardest one is observed. The PDG'09 values are used for physical constants. The results for the cross sections are given in Table 1 for different beam energies. All other results are for E p = 920 GeV.
The uncertainty of the the Born process integration is less than 0.01%. The generator error for the process with the one-loop radiative correction (RC) is about 0.2%. This uncertainty is dominated by the precision of the sub process pe → peγγ in the neighborhood of the soft and hard bremsstrahlung separation border E γ ′ ∼ E SH γ . The large peak of the radiative process here results in the systematic error of ∼ 0.1 − 0.2% in the total cross section. Better precision can be achieved by increasing the integration time and computer resources used. However, since the total error of the ZEUS luminosity measurement is about 1 − 2% the theoretical error of about 0.2% is acceptable for the cross section calculation. The uncertainty for e + e − -pair production is about 0.1% of the total electron tagger signal.
The Born cross section and the photon energy spectrum coincide with the classical Bethe-Heitler formula (1) predictions. When the radiative corrections are applied, the total cross section is still equal to the classical one within the integration error. But the energy spectrum is slightly different at very high E γ values (Fig. 5). This is exactly the region of the high systematic error of the radiative sub-process integration.
However, this region is suppressed by the small acceptance of the photon spectrometer [1]. Therefore there is no significant difference in the energy spectrum after the detector acceptance is applied (Fig. 6). The angular spectrum of the photons differs from the expected classical shape (2) (see Fig. 7). However, the divergence of the electron beam in ZEUS (∼ 0.09 − 0.21 mrad [1]) is much greater than the mean photon emission angle. Thus the actual angular spectrum shape should not affect the luminosity measurement. To estimate the expected effect we consider the size of the photon detector. We estimate the fraction of events outside the spectrometer depending on the angular distribution used. The exit window of the ZEUS spectrometer system has a diameter of about 10 cm and is located at 92.5 m from the interaction point [1]. To simulate the angular spread of the incoming electron beam we used a Gaussian distribution with σ x = 21 · 10 −4 , σ y = 9 · 10 −4 . The results, ignoring the angular acceptance of the spectrometer, are given in Table 2. The effect of the photon angular spectrum is less than 0.05%. This uncertainty is below the required generator error and can be neglected. If one ignores the scattering angle and only considers Θ γ = 0 for all photons the difference is about 0.4%.
The electron energy spectrum obtained from the generator is slightly above the classical formula (see Fig. 8). This arises from the additional e + e − -pair production process (5). The difference increases in the low-E e region up to a few percent. This contribution may affect the calibration procedure of the ZEUS photon spectrometer.
The photon energy spectrum for e − γ coincidence is presented in Fig. 9. The deviation from the classical formula at the edges of the plot is due to the additional momentum carried away by the undetected photon.

Conclusions
In this paper we have considered the processes that affect the luminosity measurement in the ZEUS experiment. The basic process pe ± → pe ± γ has been calculated including one-loop QED radiative corrections. The contribution of the lepton pair production process has also been estimated.
The following conclusions can be drawn: • The Born cross section and the photon energy spectrum coincide with the classical Bethe-Heitler formula values (1) within 0.01%; • The angular spectrum is significantly different from the classical estimate (2). However, for the ZEUS photon spectrometer the angular effect is less than 0.05%; Detector radius, mm Θ γ by generator Θ γ by Eq. (2) Table 2: Fractions of events (in %) outside the photon detector acceptance depending on the angular distribution used and the size of the ZEUS spectrometer exit window. The angular spread of initial electron beam is simulated with σx = 21 · 10 −4 , σy = 9 · 10 −4 .  • The contribution of the higher-order effects to the total cross section of the inclusive pe ± → pe ± γ process is within the 0.2% value; • The rate of the electrons in the detector is about 1−2% above the the Bethe-Heitler formula prediction due to the additional e + e − -pair production process.