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, ROOT::Math::XYZVector vertex, ROOT::Math::LorentzRotation 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, ROOT::Math::XYZVector vertex, ROOT::Math::LorentzRotation 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 31 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 155 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,
		ROOT::Math::XYZVector	vertex,
		ROOT::Math::LorentzRotation	boost
	)

Generates one single event.

Parameters

mcGraph	Reference to the MonteCarlo graph into which the generated particles will be stored.
vertex	Production vertex.
boost	Lorentz boost vector.

Returns: Returns the weight of the event.

Definition at line 80 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 101 of file BBBrem.h.

101 { return m_crossSection; }

Belle2::BBBrem::m_crossSection

double m_crossSection

The cross-section in millibarns.

Definition: BBBrem.h:149

◆ 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 106 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 116 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 111 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 121 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 126 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.

Parameters

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.
densitymode	TODO
densityparameter	TODO

Definition at line 21 of file BBBrem.cc.

◆ storeParticle()

void storeParticle	(	MCParticleGraph &	mcGraph,
		const double *	mom,
		int	pdg,
		ROOT::Math::XYZVector	vertex,
		ROOT::Math::LorentzRotation	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].
pdg	The PDG code of the particle.
vertex	The vertex of the particle in [mm].
boost	Lorentz boost vector.
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 382 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 136 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 140 of file BBBrem.h.

◆ m_weightCount

long m_weightCount

protected

Internal weighted event counter.

Used to calculate the cross-section.

Definition at line 138 of file BBBrem.h.

◆ m_weightCountOver

long m_weightCountOver

protected

Internal overweighted event counter.

Used to calculate the cross-section.

Definition at line 139 of file BBBrem.h.

◆ twopi

constexpr double twopi = 2.0 * M_PI

staticconstexprprivate

2*pi.

To keep things short.

Definition at line 178 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

Public Member Functions

Protected Member Functions

Protected Attributes

Private Attributes

Static Private Attributes

Detailed Description

Member Function Documentation

◆ calcOutgoingLeptonsAndWeight()

◆ generateEvent()

◆ getCrossSection()

◆ getCrossSectionError()

◆ getCrossSectionErrorOver()

◆ getCrossSectionOver()

◆ getMaxWeightDelivered()

◆ getSumWeightDelivered()

◆ init()

◆ storeParticle()

Member Data Documentation

◆ m_eventCount

◆ m_maxWeight

◆ m_weightCount

◆ m_weightCountOver

◆ twopi