Calculate the outgoing leptons and the event weight for one single radiative Bhabha scattering.
161 if (gRandom->Uniform() <
ac) {
162 double temp1 = gRandom->Uniform();
165 z = 1.0 / (temp1 * (expm1(
a1 * gRandom->Uniform())));
167 z =
z0 / gRandom->Uniform();
171 double y =
rme2s * z;
173 double temp1 =
pb *
pb -
eb * q0;
174 double temp2 = temp1 * temp1 -
rme2 - q0 * q0;
178 B2WARNING(
"BBBrem: y too large: delta_t^2 = " << temp2 <<
" !!!");
181 double tmin = -2.0 * (temp1 + sqrt(temp2));
182 double tmax =
rme2 *
s * y * y / tmin;
186 double w2 = sy +
rme2;
188 double rlamx = sqrt(temp1 * temp1 - 4.0 * w2 * tmax);
191 temp1 = rlamx - temp1;
193 temp1 = -4.0 * w2 * tmax / (rlamx + temp1);
201 b = exp(gRandom->Uniform() * log1p(2.0 * sy / temp1));
202 t = -b * z * z *
rme2 / ((b - 1) * (b * z + b - 1));
206 double rlam = sqrt((sy - t) * (sy - t) - 4 *
rme2 * t);
207 double eps = 4 *
rme2 * w2 / (rlam * (rlam + w2 +
rme2 - t));
208 double rl = log((2 + eps) / eps);
211 double vgam = eps * (expm1(gRandom->Uniform() * rl));
212 double cgam = 1.0 - vgam;
213 double sgam = sqrt(vgam * (2 - vgam));
216 double phi =
twopi * gRandom->Uniform();
217 double phig =
twopi * gRandom->Uniform();
220 double ql = (2.0 *
eb * q0 - t) *
rin2pb;
221 double qt = sqrt((tmax - t) * (t - tmin)) *
rin2pb;
223 q[0] = qt * sin(phi);
224 q[1] = qt * cos(phi);
229 q2[0] =
q1[0] - q[0];
230 q2[1] =
q1[1] - q[1];
231 q2[2] =
q1[2] - q[2];
232 q2[3] =
q1[3] - q[3];
237 double rin2w = 0.5 / w;
238 double rinr0w = 1.0 / (r0 + w);
239 double eta = -(sy + 2 * w * q0 + t) * rin2w * rinr0w;
240 double phat1 = -q[0] * (1 + eta);
241 double phat2 = -q[1] * (1 + eta);
242 double phat3 =
pb * eta - ql * (1 + eta);
243 double phatl = rlam * rin2w;
244 double phatt = sqrt(phat1 * phat1 + phat2 * phat2);
245 double sfhat = phat1 / phatt;
246 double cfhat = phat2 / phatt;
247 double sthat = phatt / phatl;
250 vthat = sthat * sthat / (1 - sqrt(1 - sthat * sthat));
252 vthat = sthat * sthat / (1 + sqrt(1 - sthat * sthat));
254 double cthat = vthat - 1.0;
257 double sfg = sin(phig);
258 double cfg = cos(phig);
260 temp2 = cthat * sgam * cfg + sthat * cgam;
261 double veg = vthat + vgam - vthat * vgam - sthat * sgam * cfg;
263 qkhat[3] = sy * rin2w;
264 qkhat[0] = qkhat[3] * (cfhat * temp1 + sfhat * temp2);
265 qkhat[1] = qkhat[3] * (-sfhat * temp1 + cfhat * temp2);
266 qkhat[2] = qkhat[3] * (veg - 1.0);
269 temp1 =
pb * qkhat[2];
271 temp2 = (
rme2 * qkhat[3] * qkhat[3] +
pb *
pb * (qkhat[0] * qkhat[0] + qkhat[1] * qkhat[1])) / (
eb * qkhat[3] + temp1);
273 temp2 =
eb * qkhat[3] - temp1;
276 qk[3] = (temp2 + qkhat[3] * q[3] + qkhat[0] * q[0] + qkhat[1] * q[1] + qkhat[2] * q[2]) / w;
277 temp1 = (
qk[3] + qkhat[3]) * rinr0w;
278 qk[0] = qkhat[0] + temp1 * q[0];
279 qk[1] = qkhat[1] + temp1 * q[1];
280 qk[2] = qkhat[2] + temp1 * (-
pb + ql);
296 double c1 = log1p(z) / log((2.0 + eps) / eps);
300 double vnumx = sqrt(temp1 * temp1 - 4.0 *
rme2 * tmax) + temp1;
305 vdenx = sqrt(temp1 * temp1 - 4.0 * w2 * tmax) - temp1;
307 vdenx = -4.0 * w2 * tmax / (sqrt(temp1 * temp1 - 4.0 * w2 * tmax) + temp1);
310 double vnumn = sqrt(temp1 * temp1 - 4.0 *
rme2 * tmin) + temp1;
315 vdenn = sqrt(temp1 * temp1 - 4.0 * w2 * tmin) - temp1;
317 vdenn = -4.0 * w2 * tmin / (sqrt(temp1 * temp1 - 4.0 * w2 * tmin) + temp1);
319 double c2 = 2.0 *
rls / log((vnumx * vdenn) / (vdenx * vnumn));
323 double rhat4 = -(2.0 *
rme2 * y + (1 - y) * t) * rin2w;
324 double etar = rhat4 * rinr0w;
325 double rhat1 = -q[0] * (1 + etar);
326 double rhat2 = -q[1] * (1 + etar);
327 double rhat3 = rlabl + (
pb - ql) * etar;
328 double zz =
s * (rhat4 * qkhat[3] - rhat1 * qkhat[0] - rhat2 * qkhat[1] - rhat3 * qkhat[2]);
331 double s1 = 4.0 *
eb * (
eb -
qk[3]);
332 double d1 = sy * rlam * (eps + vgam) * rin2w * rin2w;
333 double d2 = 0.5 * sy;
337 double rind1 = 1.0 / d1;
338 double rind12 = rind1 * rind1;
339 double rind2 = 1.0 / d2;
340 double rind22 = rind2 * rind2;
341 temp1 =
s + t - 2 * d2;
342 temp2 = s1 + t + 2 * d1;
343 double aa0 = (
s *
s + s1 * s1 + temp1 * temp1 + temp2 * temp2) * rind1 * rind2 * (-t);
344 double aa1 = -4.0 *
rme2 * zz * zz * rind12 * rind22;
345 double aa2 = -8.0 *
rme2 * (d1 * d1 + d2 * d2) * rind1 * rind2;
346 double aa3 = 16.0 *
rme2 *
rme2 * (d1 - d2) * zz * rind12 * rind22;
347 double aa4 = -16.0 *
rme2 *
rme2 *
rme2 * (d1 - d2) * (d1 - d2) * rind12 * rind22;
348 double rmex = aa0 + aa1 + aa2 + aa3 + aa4;
351 double rmap = 4.0 *
s *
s * rind1 * rind2 * (-t) * c1 * c2;
358 B2WARNING(
"BBBrem: Weight is nan! Setting the weight to zero.");
373 if (abs(t) < tc)
weight = 0.0;
375 if (t != tc)
weight *= t * t / (t - tc) / (t - tc);