Because of the threshold behavior, the algorithm is slightly different from that for the Compton process. The mesh for q cannot be equally spaced. Select q's such that all the threshold points 1/n (, = integer part of ) are included, that q are equally spaced between successive thresholds, and that the spaces are not very diffrent. Thus, the total number of q's may not exactly equal to MQ.
The functions stored in the array GG(k,n,i,j,l) are (k=1) and (k=2). The integral GINT(k,n,i,j,l) is calculated as in the Compton case. The sum of GINT(k,n,MY,j,l) over n does not make sense because the result would be a quite discontinuos function of and q, which makes the interpolation inaccurate. However, for given q, the sum from n=n(q) to MPH is continuous where n(q) is the minimum integer which does not exceed 1/q. Thus, the sum over n=n to MPH is stored in GALL(k,n,j,l).
For given initial condition, calculate and q. Then, interpolate
GALL(k,n(q),j,l) for j and l and calculate the
total probability P by summing for appropriate polarization:
Generate a uniform random number and reject an event if .
In this case the helicity of the photon should change to
If , decide to create a pair. The number of laser photons to absorb is determined from GINT(*,n,MY,j,l) as in the Compton case.
To determine the electron energy, generate another random number in the range (-1,+1). Then, y () is determined from as in the Compton case. Adopt -y instead of y if .
The helicity of electron is given by
The positron momentum is calculated by the momentum conservation
and the helicity from the above formula with y replaced by -y.