11 #include <generators/bbbrem/BBBrem.h>
12 #include <framework/logging/Logger.h>
13 #include <framework/gearbox/Const.h>
14 #include <mdst/dataobjects/MCParticleGraph.h>
23 void BBBrem::init(
double cmsEnergy,
double minPhotonEFrac,
bool unweighted,
double maxWeight,
int densitymode,
24 double densityparameter)
26 m_cmsEnergy = cmsEnergy;
27 m_photonEFrac = minPhotonEFrac;
29 m_unweighted = unweighted;
30 m_maxWeight = maxWeight;
31 m_densityCorrectionMode = densitymode;
32 m_densityCorrectionParameter = densityparameter;
34 m_maxWeightDelivered = 0.0;
35 m_sumWeightDelivered = 0.0;
36 m_sumWeightDeliveredSqr = 0.0;
39 if ((minPhotonEFrac <= 0.0) || (minPhotonEFrac >= 1.0)) {
40 B2ERROR(
"BBBrem: The minimum photon energy fraction has to be in the range ]0,1[ !");
44 B2DEBUG(100,
"BBBrem: Center of mass energy: " << cmsEnergy);
45 B2DEBUG(100,
"BBBrem: Minimum photon energy: " << minPhotonEFrac <<
" * beam energy");
48 alpha = Const::fineStrConst;
49 rme = Const::electronMass;
52 s = cmsEnergy * cmsEnergy;
56 z0 = minPhotonEFrac / (1.0 - minPhotonEFrac);
59 a1 = log((1.0 + z0) / z0);
62 a2 = (log1p(z0)) / z0;
64 sigapp = 8.0 * (alpha * alpha * alpha) / rme2 * (-log(rme2s)) * (a1 + a2) * tomb;
65 B2DEBUG(100,
"BBBrem: Approximate cross section: " << sigapp <<
" millibarn");
69 pb = sqrt(eb * eb - rme2);
82 double BBBrem::generateEvent(
MCParticleGraph& mcGraph, TVector3 vertex, TLorentzRotation boost)
88 calcOutgoingLeptonsAndWeight();
89 if (weight > m_maxWeightDelivered) m_maxWeightDelivered = weight;
90 if (weight > m_maxWeight) {
91 B2INFO(
"BBBrem: OVERWEIGHT, w=" << weight <<
", nw=" << m_weightCountOver <<
", sumw=" << m_sumWeightDeliveredOver);
93 m_sumWeightDeliveredOver = m_sumWeightDeliveredOver + (weight - m_maxWeight);
94 m_sumWeightDeliveredSqrOver = m_sumWeightDeliveredSqrOver + (weight - m_maxWeight) * (weight - m_maxWeight);
97 m_sumWeightDelivered += weight;
98 m_sumWeightDeliveredSqr += weight * weight;
100 ran = gRandom->Uniform();
101 }
while (weight <= ran * m_maxWeight);
103 calcOutgoingLeptonsAndWeight();
104 if (weight > m_maxWeightDelivered) m_maxWeightDelivered = weight;
106 m_sumWeightDelivered += weight;
107 m_sumWeightDeliveredSqr += weight * weight;
111 storeParticle(mcGraph, p1, 11, vertex, boost,
false,
true);
112 storeParticle(mcGraph, q1, -11, vertex, boost,
false,
true);
116 bool swapflag = (gRandom->Uniform() > 0.5) ? true :
false;
119 for (
int i = 0; i < 4; i++) tmp[i] = p2[i];
120 for (
int i = 0; i < 4; i++) p2[i] = q2[i];
121 for (
int i = 0; i < 4; i++) q2[i] = tmp[i];
128 storeParticle(mcGraph, p2, 11, vertex, boost,
false,
false);
129 storeParticle(mcGraph, q2, -11, vertex, boost,
false,
false);
130 storeParticle(mcGraph, qk, 22, vertex, boost,
false,
false);
132 if (!m_unweighted)
return weight;
139 if (m_weightCount > 0.0) {
140 m_crossSection = m_sumWeightDelivered / m_weightCount;
141 m_crossSectionError = sqrt(m_sumWeightDeliveredSqr - m_sumWeightDelivered * m_sumWeightDelivered / m_weightCount) /
145 m_crossSectionOver = m_sumWeightDeliveredOver / m_weightCount;
146 m_crossSectionErrorOver = sqrt(abs(((m_sumWeightDeliveredSqrOver) / m_weightCount - m_crossSectionOver * m_crossSectionOver) /
157 void BBBrem::calcOutgoingLeptonsAndWeight()
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);
283 p2[0] = -q2[0] - qk[0];
284 p2[1] = -q2[1] - qk[1];
285 p2[2] = -q2[2] - qk[2];
286 p2[3] = -q2[3] - qk[3] + m_cmsEnergy;
289 if (qk[3] < eb * m_photonEFrac) {
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));
322 double rlabl = (t - 2.0 * rme2 * y) * rin2pb;
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;
354 weight = rmex / rmap * sigapp;
357 if (std::isnan(weight)) {
358 B2WARNING(
"BBBrem: Weight is nan! Setting the weight to zero.");
370 double tc = m_densityCorrectionParameter;
372 if (m_densityCorrectionMode == 1) {
373 if (abs(t) < tc) weight = 0.0;
374 }
else if (m_densityCorrectionMode == 2) {
375 if (t != tc) weight *= t * t / (t - tc) / (t - tc);
384 void BBBrem::storeParticle(
MCParticleGraph& mcGraph,
const double* mom,
int pdg, TVector3 vertex, TLorentzRotation boost,
385 bool isVirtual,
bool isInitial)
391 part.setStatus(MCParticle::c_IsVirtual);
392 }
else if (isInitial) {
393 part.setStatus(MCParticle::c_Initial);
397 part.addStatus(MCParticle::c_PrimaryParticle);
400 if (pdg == 22 && !isVirtual) {
401 part.addStatus(MCParticle::c_IsISRPhoton);
402 part.addStatus(MCParticle::c_IsFSRPhoton);
406 part.setFirstDaughter(0);
407 part.setLastDaughter(0);
408 part.setMomentum(TVector3(mom[0], mom[1], -mom[2]));
409 part.setEnergy(mom[3]);
410 part.setMassFromPDG();
415 TVector3 v3 = part.getProductionVertex();
417 part.setProductionVertex(v3);
418 part.setValidVertex(
true);
422 TLorentzVector p4 = part.get4Vector();