Belle II Software  release-05-01-25
BBBrem.cc
1 /**************************************************************************
2  * BASF2 (Belle Analysis Framework 2) *
3  * Copyright(C) 2012 Belle II Collaboration *
4  * *
5  * Author: The Belle II Collaboration *
6  * Contributors: Andreas Moll, Hiroyuki Nakayama *
7  * *
8  * This software is provided "as is" without any warranty. *
9  **************************************************************************/
10 
11 #include <generators/bbbrem/BBBrem.h>
12 #include <framework/logging/Logger.h>
13 #include <framework/gearbox/Const.h>
14 #include <mdst/dataobjects/MCParticleGraph.h>
15 
16 #include <TRandom3.h>
17 //#include <cmath>
18 
19 using namespace std;
20 using namespace Belle2;
21 
22 
23 void BBBrem::init(double cmsEnergy, double minPhotonEFrac, bool unweighted, double maxWeight, int densitymode,
24  double densityparameter)
25 {
26  m_cmsEnergy = cmsEnergy;
27  m_photonEFrac = minPhotonEFrac;
28 
29  m_unweighted = unweighted;
30  m_maxWeight = maxWeight;
31  m_densityCorrectionMode = densitymode;
32  m_densityCorrectionParameter = densityparameter;
33 
34  m_maxWeightDelivered = 0.0;
35  m_sumWeightDelivered = 0.0;
36  m_sumWeightDeliveredSqr = 0.0;
37  m_weightCount = 0;
38 
39  if ((minPhotonEFrac <= 0.0) || (minPhotonEFrac >= 1.0)) {
40  B2ERROR("BBBrem: The minimum photon energy fraction has to be in the range ]0,1[ !");
41  return;
42  }
43 
44  B2DEBUG(100, "BBBrem: Center of mass energy: " << cmsEnergy);
45  B2DEBUG(100, "BBBrem: Minimum photon energy: " << minPhotonEFrac << " * beam energy");
46 
47  //Initialize the constants (in order to be consistent with the FORTRAN source code)
48  alpha = Const::fineStrConst;
49  rme = Const::electronMass;
50 
51  //Initialize the derived constants
52  s = cmsEnergy * cmsEnergy;
53  rme2 = rme * rme;
54  rme2s = rme2 / s;
55  rls = -log(rme2s);
56  z0 = minPhotonEFrac / (1.0 - minPhotonEFrac);
57 
58  //Approximate total cross section
59  a1 = log((1.0 + z0) / z0);
60 // a2 = (log(1.0 + z0)) / z0;
61 // Expression 'log(1 + x)' is replaced by 'log1p(x)' to avoid loss of precision.
62  a2 = (log1p(z0)) / z0;
63  ac = a1 / (a1 + a2);
64  sigapp = 8.0 * (alpha * alpha * alpha) / rme2 * (-log(rme2s)) * (a1 + a2) * tomb;
65  B2DEBUG(100, "BBBrem: Approximate cross section: " << sigapp << " millibarn");
66 
67  //Initial-state momenta
68  eb = 0.5 * cmsEnergy;
69  pb = sqrt(eb * eb - rme2);
70  rin2pb = 0.5 / pb;
71  p1[0] = 0.0;
72  p1[1] = 0.0;
73  p1[2] = -pb;
74  p1[3] = eb;
75  q1[0] = 0.0;
76  q1[1] = 0.0;
77  q1[2] = pb;
78  q1[3] = eb;
79 }
80 
81 
82 double BBBrem::generateEvent(MCParticleGraph& mcGraph, TVector3 vertex, TLorentzRotation boost)
83 {
84 
85  if (m_unweighted) {
86  double ran = 0.0;
87  do {
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);
92  m_weightCountOver++;
93  m_sumWeightDeliveredOver = m_sumWeightDeliveredOver + (weight - m_maxWeight);
94  m_sumWeightDeliveredSqrOver = m_sumWeightDeliveredSqrOver + (weight - m_maxWeight) * (weight - m_maxWeight);
95  }
96  m_weightCount++;
97  m_sumWeightDelivered += weight;
98  m_sumWeightDeliveredSqr += weight * weight;
99 
100  ran = gRandom->Uniform();
101  } while (weight <= ran * m_maxWeight);
102  } else {
103  calcOutgoingLeptonsAndWeight();
104  if (weight > m_maxWeightDelivered) m_maxWeightDelivered = weight;
105  m_weightCount++;
106  m_sumWeightDelivered += weight;
107  m_sumWeightDeliveredSqr += weight * weight;
108  }
109 
110  //Store the incoming particles as virtual particles, the outgoing particles as real particles
111  storeParticle(mcGraph, p1, 11, vertex, boost, false, true);
112  storeParticle(mcGraph, q1, -11, vertex, boost, false, true);
113 
114  // BBBrem emits gammma only from electron(p1)
115  // To emit gamma from positron we need to swap electron and positron here
116  bool swapflag = (gRandom->Uniform() > 0.5) ? true : false;
117  if (swapflag) {
118  double tmp[4];
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];
122  p2[2] = -p2[2];
123  q2[2] = -q2[2];
124  qk[2] = -qk[2];
125  }
126 
127  //Store the outgoing particles as real particles
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);
131 
132  if (!m_unweighted) return weight;
133  return 1.0;
134 }
135 
136 
137 void BBBrem::term()
138 {
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) /
142  m_weightCount;
143 
144  if (m_unweighted) {
145  m_crossSectionOver = m_sumWeightDeliveredOver / m_weightCount;
146  m_crossSectionErrorOver = sqrt(abs(((m_sumWeightDeliveredSqrOver) / m_weightCount - m_crossSectionOver * m_crossSectionOver) /
147  m_weightCount));
148  }
149  }
150 }
151 
152 
153 //=========================================================================
154 // Protected methods
155 //=========================================================================
156 
157 void BBBrem::calcOutgoingLeptonsAndWeight()
158 {
159  //Generate z
160  double z;
161  if (gRandom->Uniform() < ac) {
162  double temp1 = gRandom->Uniform();
163 // z = 1.0 / (temp1 * (exp(a1 * gRandom->Uniform()) - 1.0));
164 // Expression 'exp(x) - 1' is replaced by 'expm1(x)' to avoid loss of precision.
165  z = 1.0 / (temp1 * (expm1(a1 * gRandom->Uniform())));
166  } else {
167  z = z0 / gRandom->Uniform();
168  }
169 
170  //Bounds on t
171  double y = rme2s * z;
172  double q0 = eb * y;
173  double temp1 = pb * pb - eb * q0;
174  double temp2 = temp1 * temp1 - rme2 - q0 * q0;
175 
176  //If temp2<0 (very very rare): set weight to 0
177  if (temp2 < 0.0) {
178  B2WARNING("BBBrem: y too large: delta_t^2 = " << temp2 << " !!!");
179  weight = 0.0;
180  } else {
181  double tmin = -2.0 * (temp1 + sqrt(temp2));
182  double tmax = rme2 * s * y * y / tmin;
183 
184  //Generate t
185  double sy = s * y;
186  double w2 = sy + rme2;
187  temp1 = sy + tmax;
188  double rlamx = sqrt(temp1 * temp1 - 4.0 * w2 * tmax);
189 
190  if (temp1 <= 0.0) {
191  temp1 = rlamx - temp1;
192  } else {
193  temp1 = -4.0 * w2 * tmax / (rlamx + temp1);
194  }
195 
196  double t = 0.0;
197  do {
198  double b;
199 // b = exp(gRandom->Uniform() * log(1.0 + 2.0 * sy / temp1));
200 // Expression 'log(1 + x)' is replaced by 'log1p(x)' to avoid loss of precision.
201  b = exp(gRandom->Uniform() * log1p(2.0 * sy / temp1));
202  t = -b * z * z * rme2 / ((b - 1) * (b * z + b - 1));
203  } while (t <= tmin);
204 
205  //Generate cgam
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);
209 // double vgam = eps * (exp(gRandom->Uniform() * rl) - 1.0);
210 // Expression 'exp(x) - 1' is replaced by 'expm1(x)' to avoid loss of precision.
211  double vgam = eps * (expm1(gRandom->Uniform() * rl));
212  double cgam = 1.0 - vgam;
213  double sgam = sqrt(vgam * (2 - vgam));
214 
215  //Generate azimuthal angles
216  double phi = twopi * gRandom->Uniform();
217  double phig = twopi * gRandom->Uniform();
218 
219  //Construct momentum transfer q(mu)
220  double ql = (2.0 * eb * q0 - t) * rin2pb;
221  double qt = sqrt((tmax - t) * (t - tmin)) * rin2pb;
222  double q[4];
223  q[0] = qt * sin(phi);
224  q[1] = qt * cos(phi);
225  q[2] = ql;
226  q[3] = q0;
227 
228  //Construct momentum of outgoing positron in lab frame
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];
233 
234  //Find euler angles of p1(mu) in cm frame
235  double r0 = eb + q0;
236  double w = sqrt(w2);
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;
248  double vthat;
249  if (phat3 > 0.0) {
250  vthat = sthat * sthat / (1 - sqrt(1 - sthat * sthat));
251  } else {
252  vthat = sthat * sthat / (1 + sqrt(1 - sthat * sthat));
253  }
254  double cthat = vthat - 1.0;
255 
256  //Rotate using these euler angles to get the qk direction in the cm
257  double sfg = sin(phig);
258  double cfg = cos(phig);
259  temp1 = sgam * sfg;
260  temp2 = cthat * sgam * cfg + sthat * cgam;
261  double veg = vthat + vgam - vthat * vgam - sthat * sgam * cfg;
262  double qkhat[4];
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);
267 
268  //Boost the photon momentum to the lab frame
269  temp1 = pb * qkhat[2];
270  if (temp1 > 0.0) {
271  temp2 = (rme2 * qkhat[3] * qkhat[3] + pb * pb * (qkhat[0] * qkhat[0] + qkhat[1] * qkhat[1])) / (eb * qkhat[3] + temp1);
272  } else {
273  temp2 = eb * qkhat[3] - temp1;
274  }
275 
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);
281 
282  //Construct p2 by momentum conservation
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;
287 
288  //Impose cut on the photon energy: qk[3]>eb*rk0
289  if (qk[3] < eb * m_photonEFrac) {
290  weight = 0.0;
291  } else {
292  //The event is now accepted: compute matrix element and weight
293  //Compute fudge factor c1
294 // double c1 = log(1 + z) / log((2 + eps) / eps);
295 // Expression 'log(1 + x)' is replaced by 'log1p(x)' to avoid loss of precision.
296  double c1 = log1p(z) / log((2.0 + eps) / eps);
297 
298  //Compute fudge factor c2
299  temp1 = sy - tmax;
300  double vnumx = sqrt(temp1 * temp1 - 4.0 * rme2 * tmax) + temp1;
301  temp1 = sy + tmax;
302 
303  double vdenx;
304  if (temp1 < 0.0) {
305  vdenx = sqrt(temp1 * temp1 - 4.0 * w2 * tmax) - temp1;
306  } else {
307  vdenx = -4.0 * w2 * tmax / (sqrt(temp1 * temp1 - 4.0 * w2 * tmax) + temp1);
308  }
309  temp1 = sy - tmin;
310  double vnumn = sqrt(temp1 * temp1 - 4.0 * rme2 * tmin) + temp1;
311  temp1 = sy + tmin;
312 
313  double vdenn;
314  if (temp1 < 0.0) {
315  vdenn = sqrt(temp1 * temp1 - 4.0 * w2 * tmin) - temp1;
316  } else {
317  vdenn = -4.0 * w2 * tmin / (sqrt(temp1 * temp1 - 4.0 * w2 * tmin) + temp1);
318  }
319  double c2 = 2.0 * rls / log((vnumx * vdenn) / (vdenx * vnumn));
320 
321  //Compute vector (small) r in cm frame, and (big) z
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]);
329 
330  //The other invariants
331  double s1 = 4.0 * eb * (eb - qk[3]);
332  double d1 = sy * rlam * (eps + vgam) * rin2w * rin2w;
333  double d2 = 0.5 * sy;
334 
335  //The exact matrix element
336  //Kleiss-burkhardt cross section multiplied by t^2
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;
349 
350  //The approximate matrix element without c1,2, multiplied by t^2
351  double rmap = 4.0 * s * s * rind1 * rind2 * (-t) * c1 * c2;
352 
353  //The weight
354  weight = rmex / rmap * sigapp;
355 
356  //Isnan check (not sure if this is ok)
357  if (std::isnan(weight)) {
358  B2WARNING("BBBrem: Weight is nan! Setting the weight to zero.");
359  weight = 0.0;
360  }
361 
362 
363  //========================================================
364  // beam size effect (or density effect)
365  // ref: arxiv:hep-ph/9401333
366 
367  // tc = (hbarc/simga_y)^2
368 // double tc = 1.68e-17; //SuperKEKB LER (sigma_y*=48nm)
369  //double tc = 9.81e-18; //SuperKEKB HER (sigma_y*=63nm)
370  double tc = m_densityCorrectionParameter;
371 // int cutflag = 0 ; // 0: no cut, 1: hard cut, 2: soft cut
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); // t<0<tc, always t!=tc
376  }
377  //========================================================
378 
379  }
380  }
381 }
382 
383 
384 void BBBrem::storeParticle(MCParticleGraph& mcGraph, const double* mom, int pdg, TVector3 vertex, TLorentzRotation boost,
385  bool isVirtual, bool isInitial)
386 {
387 
388  // RG 6/25/14 Add new flag for ISR "c_Initial"
390  if (isVirtual) {
391  part.setStatus(MCParticle::c_IsVirtual);
392  } else if (isInitial) {
393  part.setStatus(MCParticle::c_Initial);
394  }
395 
396  //all particles from a generator are primary
397  part.addStatus(MCParticle::c_PrimaryParticle);
398 
399  //all non virtual photons from BBBREM are ISR or FSR
400  if (pdg == 22 && !isVirtual) {
401  part.addStatus(MCParticle::c_IsISRPhoton);
402  part.addStatus(MCParticle::c_IsFSRPhoton);
403  }
404 
405  part.setPDG(pdg);
406  part.setFirstDaughter(0);
407  part.setLastDaughter(0);
408  part.setMomentum(TVector3(mom[0], mom[1], -mom[2])); //Switch direction, because electrons fly into the +z direction at Belle II
409  part.setEnergy(mom[3]);
410  part.setMassFromPDG();
411 // part.setDecayVertex(0.0, 0.0, 0.0);
412 
413  //set the production vertex of non initial particles
414  if (!isInitial) {
415  TVector3 v3 = part.getProductionVertex();
416  v3 = v3 + vertex;
417  part.setProductionVertex(v3);
418  part.setValidVertex(true);
419  }
420 
421  //If boosting is enable boost the particles to the lab frame
422  TLorentzVector p4 = part.get4Vector();
423  p4 = boost * p4;
424  part.set4Vector(p4);
425 }
Belle2::MCParticleGraph
Class to build, validate and sort a particle decay chain.
Definition: MCParticleGraph.h:48
Belle2
Abstract base class for different kinds of events.
Definition: MillepedeAlgorithm.h:19
Belle2::MCParticleGraph::addParticle
GraphParticle & addParticle()
Add new particle to the graph.
Definition: MCParticleGraph.h:305
Belle2::MCParticleGraph::GraphParticle
Class to represent Particle data in graph.
Definition: MCParticleGraph.h:86