BBBrem Class Reference

Generator for low scattering angle radiative Bhabha events (Beam-Beam Bremsstrahlung). More...

#include <BBBrem.h>

Public Member Functions
	BBBrem ()
	Constructor.
	~BBBrem ()
	Destructor.
void	init (double cmsEnergy, double minPhotonEFrac, bool unweighted=true, double maxWeight=2000.0, int densitymode=1, double densityparameter=1.68e-17)
	Initializes the generator. More...
double	generateEvent (MCParticleGraph &mcGraph, TVector3 vertex, TLorentzRotation boost)
	Generates one single event. More...
double	getCrossSection ()
	Returns the total cross section of the generated process in millibarn. More...
double	getCrossSectionError ()
	Returns the error on the total cross section of the generated process in millibarn. More...
double	getCrossSectionOver ()
	Returns the total overweight bias cross section of the generated process in millibarn. More...
double	getCrossSectionErrorOver ()
	Returns the error on the total overweight bias cross section of the generated process in millibarn. More...
double	getMaxWeightDelivered ()
	Returns the maximum weight given by the BBBrem generation. More...
double	getSumWeightDelivered ()
	Returns the sum of all weights returned by the BBBrem generation. More...
void	term ()
	Terminates the generator.

Protected Member Functions
void	calcOutgoingLeptonsAndWeight ()
	Calculate the outgoing leptons and the event weight for one single radiative Bhabha scattering. More...
void	storeParticle (MCParticleGraph &mcGraph, const double *mom, int pdg, TVector3 vertex, TLorentzRotation boost, bool isVirtual=false, bool isInitial=false)
	Store a single generated particle into the MonteCarlo graph. More...

Protected Attributes
int	m_eventCount
	Internal event counter. More...
bool	m_unweighted
	True if BBBrem should produce unweighted events.
long	m_weightCount
	Internal weighted event counter. More...
long	m_weightCountOver
	Internal overweighted event counter. More...
double	m_maxWeight
	The maximum weight. More...
double	m_maxWeightDelivered
	The maximum weight given by the BBBrem generation.
double	m_sumWeightDelivered
	The sum of all weights returned by the BBBrem generation.
double	m_sumWeightDeliveredSqr
	The square of the sum of all weights returned by the BBBrem generation.
double	m_sumWeightDeliveredOver
	The sum of all overweights.
double	m_sumWeightDeliveredSqrOver
	The square of the sum of all overweights.
double	m_cmsEnergy
	Center of mass energy (sqrt(s)).
double	m_photonEFrac
	Minimum photon energy fraction.
double	m_crossSection
	The cross-section in millibarns.
double	m_crossSectionError
	The error on the cross-section in millibarns.
double	m_crossSectionOver
	The overweight bias cross-section in millibarns.
double	m_crossSectionErrorOver
	The overweight bias error on the cross-section in millibarns.
int	m_densityCorrectionMode
	Mode for bunch density correction.
double	m_densityCorrectionParameter
	Density correction parameter tc.

Private Attributes
double	alpha
	variable
double	rme
	in MeV
double	s
	variable
double	rme2
	variable
double	rme2s
	variable
double	rls
	variable
double	z0
	variable
double	a1
	variable
double	a2
	variable
double	ac
	variable
double	sigapp
	variable
double	eb
	variable
double	pb
	variable
double	rin2pb
	variable
double	p1 [4]
	variable
double	p2 [4]
	variable
double	q1 [4]
	variable
double	q2 [4]
	variable
double	qk [4]
	variable
double	weight
	variable

Static Private Attributes
static constexpr double	tomb = 3.8937966e5 / 1e6
	Conversion factor (hc)^2.
static constexpr double	twopi = 2.0 * M_PI
	2*pi. More...

Detailed Description

Generator for low scattering angle radiative Bhabha events (Beam-Beam Bremsstrahlung).

Based on the FORTRAN source code from R. Kleiss. In order to be consistent with the original source code, the variable names are kept and might therefore violate the coding conventions.

Original authors: r. Kleiss and h. Burkhardt Reference: [arXiv:hep-ph/9401333v1]

Definition at line 40 of file BBBrem.h.

Member Function Documentation

calcOutgoingLeptonsAndWeight()

void calcOutgoingLeptonsAndWeight ( )

protected

Calculate the outgoing leptons and the event weight for one single radiative Bhabha scattering.

The main method. A direct translation from the BBBrem Fortran source code.

Definition at line 157 of file BBBrem.cc.

{
  //Generate z
  double z;
  if (gRandom->Uniform() < ac) {
    double temp1 = gRandom->Uniform();
//     z = 1.0 / (temp1 * (exp(a1 * gRandom->Uniform()) - 1.0));
// Expression 'exp(x) - 1' is replaced by 'expm1(x)' to avoid loss of precision.
    z = 1.0 / (temp1 * (expm1(a1 * gRandom->Uniform())));
  } else {
    z = z0 / gRandom->Uniform();
  }
 
  //Bounds on t
  double y = rme2s * z;
  double q0 = eb * y;
  double temp1 = pb * pb - eb * q0;
  double temp2 = temp1 * temp1 - rme2 - q0 * q0;
 
  //If temp2<0 (very very rare): set weight to 0
  if (temp2 < 0.0) {
    B2WARNING("BBBrem: y too large: delta_t^2 = " << temp2 << " !!!");
    weight = 0.0;
  } else {
    double tmin = -2.0 * (temp1 + sqrt(temp2));
    double tmax = rme2 * s * y * y / tmin;
 
    //Generate t
    double sy = s * y;
    double w2 = sy + rme2;
    temp1 = sy + tmax;
    double rlamx = sqrt(temp1 * temp1 - 4.0 * w2 * tmax);
 
    if (temp1 <= 0.0) {
      temp1 = rlamx - temp1;
    } else {
      temp1 = -4.0 * w2 * tmax / (rlamx + temp1);
    }
 
    double t = 0.0;
    do {
      double b;
// b = exp(gRandom->Uniform() * log(1.0 + 2.0 * sy / temp1));
// Expression 'log(1 + x)' is replaced by 'log1p(x)' to avoid loss of precision.
      b = exp(gRandom->Uniform() * log1p(2.0 * sy / temp1));
      t = -b * z * z * rme2 / ((b - 1) * (b * z + b - 1));
    } while (t <= tmin);
 
    //Generate cgam
    double rlam = sqrt((sy - t) * (sy - t) - 4 * rme2 * t);
    double eps = 4 * rme2 * w2 / (rlam * (rlam + w2 + rme2 - t));
    double rl = log((2 + eps) / eps);
//     double vgam = eps * (exp(gRandom->Uniform() * rl) - 1.0);
// Expression 'exp(x) - 1' is replaced by 'expm1(x)' to avoid loss of precision.
    double vgam = eps * (expm1(gRandom->Uniform() * rl));
    double cgam = 1.0 - vgam;
    double sgam = sqrt(vgam * (2 - vgam));
 
    //Generate azimuthal angles
    double phi  = twopi * gRandom->Uniform();
    double phig = twopi * gRandom->Uniform();
 
    //Construct momentum transfer q(mu)
    double ql = (2.0 * eb * q0 - t) * rin2pb;
    double qt = sqrt((tmax - t) * (t - tmin)) * rin2pb;
    double q[4];
    q[0] = qt * sin(phi);
    q[1] = qt * cos(phi);
    q[2] = ql;
    q[3] = q0;
 
    //Construct momentum of outgoing positron in lab frame
    q2[0] = q1[0] - q[0];
    q2[1] = q1[1] - q[1];
    q2[2] = q1[2] - q[2];
    q2[3] = q1[3] - q[3];
 
    //Find euler angles of p1(mu) in cm frame
    double r0 = eb + q0;
    double w = sqrt(w2);
    double rin2w = 0.5 / w;
    double rinr0w = 1.0 / (r0 + w);
    double eta = -(sy + 2 * w * q0 + t) * rin2w * rinr0w;
    double phat1 = -q[0] * (1 + eta);
    double phat2 = -q[1] * (1 + eta);
    double phat3 = pb * eta - ql * (1 + eta);
    double phatl = rlam * rin2w;
    double phatt = sqrt(phat1 * phat1 + phat2 * phat2);
    double sfhat = phat1 / phatt;
    double cfhat = phat2 / phatt;
    double sthat = phatt / phatl;
    double vthat;
    if (phat3 > 0.0) {
      vthat = sthat * sthat / (1 - sqrt(1 - sthat * sthat));
    } else {
      vthat = sthat * sthat / (1 + sqrt(1 - sthat * sthat));
    }
    double cthat = vthat - 1.0;
 
    //Rotate using these euler angles to get the qk direction in the cm
    double sfg = sin(phig);
    double cfg = cos(phig);
    temp1 = sgam * sfg;
    temp2 = cthat * sgam * cfg + sthat * cgam;
    double veg = vthat + vgam - vthat * vgam - sthat * sgam * cfg;
    double qkhat[4];
    qkhat[3] = sy * rin2w;
    qkhat[0] = qkhat[3] * (cfhat * temp1 + sfhat * temp2);
    qkhat[1] = qkhat[3] * (-sfhat * temp1 + cfhat * temp2);
    qkhat[2] = qkhat[3] * (veg - 1.0);
 
    //Boost the photon momentum to the lab frame
    temp1 = pb * qkhat[2];
    if (temp1 > 0.0) {
      temp2 = (rme2 * qkhat[3] * qkhat[3] + pb * pb * (qkhat[0] * qkhat[0] + qkhat[1] * qkhat[1])) / (eb * qkhat[3] + temp1);
    } else {
      temp2 = eb * qkhat[3] - temp1;
    }
 
    qk[3] = (temp2 + qkhat[3] * q[3] + qkhat[0] * q[0] + qkhat[1] * q[1] + qkhat[2] * q[2]) / w;
    temp1 = (qk[3] + qkhat[3]) * rinr0w;
    qk[0] = qkhat[0] + temp1 * q[0];
    qk[1] = qkhat[1] + temp1 * q[1];
    qk[2] = qkhat[2] + temp1 * (-pb + ql);
 
    //Construct p2 by momentum conservation
    p2[0] = -q2[0] - qk[0];
    p2[1] = -q2[1] - qk[1];
    p2[2] = -q2[2] - qk[2];
    p2[3] = -q2[3] - qk[3] + m_cmsEnergy;
 
    //Impose cut on the photon energy: qk[3]>eb*rk0
    if (qk[3] < eb * m_photonEFrac) {
      weight = 0.0;
    } else {
      //The event is now accepted: compute matrix element and weight
      //Compute fudge factor c1
//       double c1 = log(1 + z) / log((2 + eps) / eps);
//       Expression 'log(1 + x)' is replaced by 'log1p(x)' to avoid loss of precision.
      double c1 = log1p(z) / log((2.0 + eps) / eps);
 
      //Compute fudge factor c2
      temp1 = sy - tmax;
      double vnumx = sqrt(temp1 * temp1 - 4.0 * rme2 * tmax) + temp1;
      temp1 = sy + tmax;
 
      double vdenx;
      if (temp1 < 0.0) {
        vdenx = sqrt(temp1 * temp1 - 4.0 * w2 * tmax) - temp1;
      } else {
        vdenx = -4.0 * w2 * tmax / (sqrt(temp1 * temp1 - 4.0 * w2 * tmax) + temp1);
      }
      temp1 = sy - tmin;
      double vnumn = sqrt(temp1 * temp1 - 4.0 * rme2 * tmin) + temp1;
      temp1 = sy + tmin;
 
      double vdenn;
      if (temp1 < 0.0) {
        vdenn = sqrt(temp1 * temp1 - 4.0 * w2 * tmin) - temp1;
      } else {
        vdenn = -4.0 * w2 * tmin / (sqrt(temp1 * temp1 - 4.0 * w2 * tmin) + temp1);
      }
      double c2 = 2.0 * rls / log((vnumx * vdenn) / (vdenx * vnumn));
 
      //Compute vector (small) r in cm frame, and (big) z
      double rlabl = (t - 2.0 * rme2 * y) * rin2pb;
      double rhat4 = -(2.0 * rme2 * y + (1 - y) * t) * rin2w;
      double etar = rhat4 * rinr0w;
      double rhat1 = -q[0] * (1 + etar);
      double rhat2 = -q[1] * (1 + etar);
      double rhat3 = rlabl + (pb - ql) * etar;
      double zz = s * (rhat4 * qkhat[3] - rhat1 * qkhat[0] - rhat2 * qkhat[1] - rhat3 * qkhat[2]);
 
      //The other invariants
      double s1 = 4.0 * eb * (eb - qk[3]);
      double d1 = sy * rlam * (eps + vgam) * rin2w * rin2w;
      double d2 = 0.5 * sy;
 
      //The exact matrix element
      //Kleiss-burkhardt cross section multiplied by t^2
      double rind1 = 1.0 / d1;
      double rind12 = rind1 * rind1;
      double rind2 = 1.0 / d2;
      double rind22 = rind2 * rind2;
      temp1 = s + t - 2 * d2;
      temp2 = s1 + t + 2 * d1;
      double aa0 = (s * s + s1 * s1 + temp1 * temp1 + temp2 * temp2) * rind1 * rind2 * (-t);
      double aa1 = -4.0 * rme2 * zz * zz * rind12 * rind22;
      double aa2 = -8.0 * rme2 * (d1 * d1 + d2 * d2) * rind1 * rind2;
      double aa3 = 16.0 * rme2 * rme2 * (d1 - d2) * zz * rind12 * rind22;
      double aa4 = -16.0 * rme2 * rme2 * rme2 * (d1 - d2) * (d1 - d2) * rind12 * rind22;
      double rmex = aa0 + aa1 + aa2 + aa3 + aa4;
 
      //The approximate matrix element without c1,2, multiplied by t^2
      double rmap = 4.0 * s * s * rind1 * rind2 * (-t) * c1 * c2;
 
      //The weight
      weight = rmex / rmap * sigapp;
 
      //Isnan check (not sure if this is ok)
      if (std::isnan(weight)) {
        B2WARNING("BBBrem: Weight is nan! Setting the weight to zero.");
        weight = 0.0;
      }
 
 
      //========================================================
      // beam size effect (or density effect)
      // ref: arxiv:hep-ph/9401333
 
      // tc = (hbarc/simga_y)^2
//       double tc = 1.68e-17;  //SuperKEKB LER (sigma_y*=48nm)
      //double tc = 9.81e-18;  //SuperKEKB HER (sigma_y*=63nm)
      double tc = m_densityCorrectionParameter;
//       int cutflag = 0 ;  // 0: no cut, 1: hard cut, 2: soft cut
      if (m_densityCorrectionMode == 1) {
        if (abs(t) < tc) weight = 0.0;
      } else if (m_densityCorrectionMode == 2) {
        if (t != tc) weight *= t * t / (t - tc) / (t - tc); // t<0<tc, always t!=tc
      }
      //========================================================
 
    }
  }
}

generateEvent()

double generateEvent	(	MCParticleGraph &	mcGraph,
		TVector3	vertex,
		TLorentzRotation	boost
	)

Generates one single event.

Parameters

mcGraph

Reference to the MonteCarlo graph into which the generated particles will be stored.

Returns

Returns the weight of the event.

Definition at line 82 of file BBBrem.cc.

getCrossSection()

double getCrossSection ( )

inline

Returns the total cross section of the generated process in millibarn.

Returns

The total cross section.

Definition at line 106 of file BBBrem.h.

{ return m_crossSectionOver; }

getCrossSectionError()

double getCrossSectionError ( )

inline

Returns the error on the total cross section of the generated process in millibarn.

Returns

The error on the total cross section.

Definition at line 111 of file BBBrem.h.

getCrossSectionErrorOver()

double getCrossSectionErrorOver ( )

inline

Returns the error on the total overweight bias cross section of the generated process in millibarn.

Returns

The error on the total overweight bias cross section.

Definition at line 121 of file BBBrem.h.

getCrossSectionOver()

double getCrossSectionOver ( )

inline

Returns the total overweight bias cross section of the generated process in millibarn.

Returns

The total overweight bias cross section.

Definition at line 116 of file BBBrem.h.

getMaxWeightDelivered()

double getMaxWeightDelivered ( )

inline

Returns the maximum weight given by the BBBrem generation.

Returns

The maximum weight given by the BBBrem generation.

Definition at line 126 of file BBBrem.h.

getSumWeightDelivered()

double getSumWeightDelivered ( )

inline

Returns the sum of all weights returned by the BBBrem generation.

Returns

The sum of all weights returned by the BBBrem generation.

Definition at line 131 of file BBBrem.h.

init()

void init	(	double	cmsEnergy,
		double	minPhotonEFrac,
		bool	unweighted = true,
		double	maxWeight = 2000.0,
		int	densitymode = 1,
		double	densityparameter = 1.68e-17
	)

Initializes the generator.

@cmsEnergy The center of mass energy in [GeV]. @minPhotonEFrac The minimum photon energy fraction. The min. photon energy is this fraction * beam energy. @unweighted Set to true to generate unweighted events (e.g. for detector simulation), set to false for weighted events. @maxWeight The maximum weight. Only required in the case of unweighted events.

Definition at line 23 of file BBBrem.cc.

storeParticle()

void storeParticle ( MCParticleGraph &  mcGraph, const double *  mom, int  pdg, TVector3  vertex, TLorentzRotation  boost, bool  isVirtual = false, bool  isInitial = false  )

protected

Store a single generated particle into the MonteCarlo graph.

Parameters

mcGraph	Reference to the MonteCarlo graph into which the particle should be stored.
mom	The 3-momentum of the particle in [GeV].
vtx	The vertex of the particle in [mm].
pdg	The PDG code of the particle.
isVirtual	If the particle is a virtual particle, such as the incoming particles, set this to true.
isInitial	If the particle is a initial particle for ISR, set this to true.

Definition at line 384 of file BBBrem.cc.

Member Data Documentation

m_eventCount

int m_eventCount

protected

Internal event counter.

Used to calculate the cross-section.

Definition at line 141 of file BBBrem.h.

m_maxWeight

double m_maxWeight

protected

The maximum weight.

Used for the event rejection procedure to produce unweighted events.

Definition at line 145 of file BBBrem.h.

m_weightCount

long m_weightCount

protected

Internal weighted event counter.

Used to calculate the cross-section.

Definition at line 143 of file BBBrem.h.

m_weightCountOver

long m_weightCountOver

protected

Internal overweighted event counter.

Used to calculate the cross-section.

Definition at line 144 of file BBBrem.h.

twopi

constexpr double twopi = 2.0 * M_PI

staticconstexprprivate

2*pi.

To keep things short.

Definition at line 182 of file BBBrem.h.

---

The documentation for this class was generated from the following files:

• generators/bbbrem/include/BBBrem.h
• generators/bbbrem/src/BBBrem.cc