Four-Vectors

  1. Constructors and basic operators
  2. Member methods for input
  3. Member methods for output
  4. Friend methods for output
  5. Operations with four-vectors
  6. Rotations and boosts
The Vec4 class gives a simple implementation of four-vectors. The member function names are based on the assumption that these represent four-momentum vectors. Thus one can get or set p_x, p_y, p_z and e, but not x, y, z or t. This is only a matter of naming, however; a Vec4 can equally well be used to store a space-time four-vector.

The Particle object contains a Vec4 p that stores the particle four-momentum, and another Vec4 vProd for the production vertex. For the latter the input/output method names are adapted to the space-time character rather than the normal energy-momentum one. Thus a user would not normally access the Vec4 classes directly, but only via the methods of the Particle class, see Particle Properties.

Nevertheless you are free to use the PYTHIA four-vectors, e.g. as part of some simple analysis code based directly on the PYTHIA output, say to define the four-vector sum of a set of particles. But note that this class was never set up to allow complete generality, only to provide the operations that are of use inside PYTHIA. There is no separate class for three-vectors, since such can easily be represented by four-vectors where the fourth component is not used.

Four-vectors have the expected functionality: they can be created, copied, added, multiplied, rotated, boosted, and manipulated in other ways. Operator overloading is implemented where reasonable. Properties can be read out, not only the components themselves but also for derived quantities such as absolute momentum and direction angles.

Constructors and basic operators

A few methods are available to create or copy a four-vector:

Vec4::Vec4(double x = 0., double y = 0., double z = 0., double t = 0.)  
creates a four-vector, by default with all components set to 0.

Vec4::Vec4(const Vec4& v)  
creates a four-vector copy of the input four-vector.

Vec4& Vec4::operator=(const Vec4& v)  
copies the input four-vector.

Vec4& Vec4::operator=(double value)  
gives a four-vector with all components set to value.

Member methods for input

The values stored in a four-vector can be modified in a few different ways:

void Vec4::reset()  
sets all components to 0.

void Vec4::p(double pxIn, double pyIn, double pzIn, double eIn)  
sets all components to their input values.

void Vec4::p(Vec4 pIn)  
sets all components equal to those of the input four-vector.

void Vec4::px(double pxIn)  
void Vec4::py(double pyIn)  
void Vec4::pz(double pzIn)  
void Vec4::e(double eIn)  
sets the respective component to the input value.

Member methods for output

A number of methods provides output of basic or derived quantities:

double Vec4::px()  
double Vec4::py()  
double Vec4::pz()  
double Vec4::e()  
gets the respective component.

double& operator[](int i)  
returns component by index, where 1 gives p_x, 2 gives p_y, 3 gives p_z, and anything else gives e.

double Vec4::mCalc()  
double Vec4::m2Calc()  
the (squared) mass, calculated from the four-vectors. If m^2 < 0 the mass is given with a minus sign, -sqrt(-m^2). Note the possible loss of precision in the calculation of E^2 - p^2; for particles the correct mass is stored separately to avoid such problems.

double Vec4::pT()  
double Vec4::pT2()  
the (squared) transverse momentum.

double Vec4::pAbs()  
double Vec4::pAbs2()  
the (squared) absolute momentum.

double Vec4::eT()  
double Vec4::eT2()  
the (squared) transverse energy, eT = e * sin(theta) = e * pT / pAbs.

double Vec4::theta()  
the polar angle, in the range 0 through pi.

double Vec4::phi()  
the azimuthal angle, in the range -pi through pi.

double Vec4::thetaXZ()  
the angle in the xz plane, in the range -pi through pi, with 0 along the +z axis.

double Vec4::pPos()  
double Vec4::pNeg()  
the combinations E+-p_z.

double Vec4::rap()  
double Vec4::eta()  
true rapidity y and pseudorapidity eta. In case of massless or inconsistent input, such as e ≤ |p_z| (within machine precision), a value of +-20 is set, where the sign reflects that of p_z.

Friend methods for output

There are also some friend methods that take one, two or three four-vectors as argument. Several of them only use the three-vector part of the four-vector.

friend ostream& operator<<(ostream&, const Vec4& v)  
writes out the values of the four components of a Vec4 and, within brackets, a fifth component being the invariant length of the four-vector, as provided by mCalc() above, and it all ended with a newline.

friend bool isnan(const Vec4& v)  
friend double isinf(const Vec4& v)  
returns true if any of the four elements is not a number or infinite. friend bool isnan(const Vec4& v)  
returns true if all four elements are finite.

friend double m(const Vec4& v1, const Vec4& v2)  
friend double m2(const Vec4& v1, const Vec4& v2)  
the (squared) invariant mass.

friend double dot3(const Vec4& v1, const Vec4& v2)  
the three-product.

friend double cross3(const Vec4& v1, const Vec4& v2)  
the cross-product.

friend double cross4(const Vec4& v1, const Vec4& v2, const Vec4& v3)  
the cross-product of three four-vectors: v_i = epsilon_{iabc} v1_a v2_b v3_c.

friend double theta(const Vec4& v1, const Vec4& v2)  
friend double costheta(const Vec4& v1, const Vec4& v2)  
the (cosine) of the opening angle between the vectors, in the range 0 through pi.

friend double phi(const Vec4& v1, const Vec4& v2)  
friend double cosphi(const Vec4& v1, const Vec4& v2)  
the (cosine) of the azimuthal angle between the vectors around the z axis, in the range 0 through pi.

friend double phi(const Vec4& v1, const Vec4& v2, const Vec4& v3)  
friend double cosphi(const Vec4& v1, const Vec4& v2, const Vec4& v3)  
the (cosine) of the azimuthal angle between the first two vectors around the direction of the third, in the range 0 through pi.

friend double RRapPhi(const Vec4& v1, const Vec4& v2)  
friend double REtaPhi(const Vec4& v1, const Vec4& v2)  
the R distance measure, in (y, phi) or (eta, phi) cylindrical coordinates, i.e. R^2 = (y_1 - y_2)^2 + (phi_1 - phi_2)^2 and equivalent.

friend bool pShift( Vec4& p1Move, Vec4& p2Move, double m1New, double m2New)  
transfer four-momentum between the two four-vectors so that they get the masses m1New and m2New, respectively. Note that p1Move and p2Move act both as input and output arguments. The method will return false if the invariant mass of the four-vectors is too small to accommodate the new masses, and then the four-vectors are not changed.

friend pair<Vec4,Vec4> getTwoPerpendicular(const Vec4& v1, const Vec4& v2)  
create a pair of four-vectors that are perpendicular to both input vectors and to each other, and have the squared norm -1.

Operations with four-vectors

Of course one should be able to add, subtract and scale four-vectors, and more:

Vec4 Vec4::operator-()  
return a vector with flipped sign for all components, while leaving the original vector unchanged.

Vec4& Vec4::operator+=(const Vec4& v)  
add a four-vector to an existing one.

Vec4& Vec4::operator-=(const Vec4& v)  
subtract a four-vector from an existing one.

Vec4& Vec4::operator*=(double f)  
multiply all four-vector components by a real number.

Vec4& Vec4::operator/=(double f)  
divide all four-vector components by a real number.

friend Vec4 operator+(const Vec4& v1, const Vec4& v2)  
add two four-vectors.

friend Vec4 operator-(const Vec4& v1, const Vec4& v2)  
subtract two four-vectors.

friend Vec4 operator*(double f, const Vec4& v)  
friend Vec4 operator*(const Vec4& v, double f)  
multiply a four-vector by a real number.

friend Vec4 operator/(const Vec4& v, double f)  
divide a four-vector by a real number.

friend double operator*(const Vec4& v1, const Vec4 v2)  
four-vector product.

There are also a few related operations that are normal member methods:

void Vec4::rescale3(double f)  
multiply the three-vector components by f, but keep the fourth component unchanged.

void Vec4::rescale4(double f)  
multiply all four-vector components by f.

void Vec4::flip3()  
flip the sign of the three-vector components, but keep the fourth component unchanged.

void Vec4::flip4()  
flip the sign of all four-vector components.

Rotations and boosts

A common task is to rotate or boost four-vectors. In case only one four-vector is affected the operation may be performed directly on it. However, in case many particles are affected, the helper class RotBstMatrix can be used to speed up operations.

void Vec4::rot(double theta, double phi)  
rotate the three-momentum with the polar angle theta and the azimuthal angle phi.

void Vec4::rotaxis(double phi, double nx, double ny, double nz)  
rotate the three-momentum with the azimuthal angle phi around the direction defined by the (n_x, n_y, n_z) three-vector.

void Vec4::rotaxis(double phi, Vec4& n)  
rotate the three-momentum with the azimuthal angle phi around the direction defined by the three-vector part of n.

void Vec4::bst(double betaX, double betaY, double betaZ)  
boost the four-momentum by beta = (beta_x, beta_y, beta_z).

void Vec4::bst(double betaX, double betaY, double betaZ, double gamma)  
boost the four-momentum by beta = (beta_x, beta_y, beta_z), where the gamma = 1/sqrt(1 - beta^2) is also input to allow better precision when beta is close to unity.

void Vec4::bst(const Vec4& p)  
boost the four-momentum by beta = (p_x/E, p_y/E, p_z/E).

void Vec4::bst(const Vec4& p, double m)  
boost the four-momentum by beta = (p_x/E, p_y/E, p_z/E), where the gamma = E/m is also calculated from input to allow better precision when beta is close to unity.

void Vec4::bstback(const Vec4& p)  
boost the four-momentum by beta = (-p_x/E, -p_y/E, -p_z/E).

void Vec4::bstback(const Vec4& p, double m)  
boost the four-momentum by beta = (-p_x/E, -p_y/E, -p_z/E), where the gamma = E/m is also calculated from input to allow better precision when beta is close to unity.

void Vec4::rotbst(const RotBstMatrix& M)  
perform a combined rotation and boost; see below for a description of the RotBstMatrix.

void Vec4::eInFrame(const Vec4& pIn)  
calculate the energy of the four-momentum in the frame of the four-momentum pIn.

For a longer sequence of rotations and boosts, and where several Vec4 are to be rotated and boosted in the same way, a more efficient approach is to define a RotBstMatrix, which forms a separate auxiliary class. You can build up this 4-by-4 matrix by successive calls to the methods of the class, such that the matrix encodes the full sequence of operations. The order in which you do these calls must agree with the imagined order in which the rotations/boosts should be applied to a four-momentum, since in general the operations do not commute.
(Mathematically you would e.g. define M = M_3 M_2 M_1 in that M p = M_3( M_2( M_1 p) ) ). That is, operations on the four-vector p are carried out in the order first M_1, then M_2 and finally M_3. Thus M_1, M_2, M_3 is also the order in which you should input rotations and boosts to M.)

RotBstMatrix::RotBstMatrix()  
creates a diagonal unit matrix, i.e. one that leaves a four-vector unchanged.

RotBstMatrix::RotBstMatrix(const RotBstMatrix& Min)  
creates a copy of the input matrix.

RotBstMatrix& RotBstMatrix::operator=(const RotBstMatrix4& Min)  
copies the input matrix.

void RotBstMatrix::rot(double theta = 0., double phi = 0.)  
rotate by this polar and azimuthal angle.

void RotBstMatrix::rot(const Vec4& p)  
rotate so that a vector originally along the +z axis becomes parallel with p. More specifically, rotate by -phi, theta and phi, with angles defined by p.

void RotBstMatrix::bst(double betaX = 0., double betaY = 0., double betaZ = 0.)  
boost by this beta vector.

void RotBstMatrix::bst(const Vec4&)  
void RotBstMatrix::bstback(const Vec4&)  
boost with a beta = p/E or beta = -p/E, respectively.

void RotBstMatrix::bst(const Vec4& p1, const Vec4& p2)  
boost so that p_1 is transformed to p_2. It is assumed that the two vectors obey p_1^2 = p_2^2.

void RotBstMatrix::toCMframe(const Vec4& p1, const Vec4& p2)  
boost and rotate to the rest frame of p_1 and p_2, with p_1 along the +z axis.

void RotBstMatrix::fromCMframe(const Vec4& p1, const Vec4& p2, bool flip = false)  
rotate and boost from the rest frame of p_1 and p_2, with p_1 along the +z axis, to the actual frame of p_1 and p_2, i.e. the inverse of the above. The option flip handles the case when p_1 is along the -z axis in the rest frame. This is accomplished by performing the rotation for p_2 and negating the rotational part.

void RotBstMatrix::toSameVframe(const Vec4& p1, const Vec4& p2)  
void RotBstMatrix::fromSameVframe(const Vec4& p1, const Vec4& p2)  
similar to toCMframe and fromCMframe above, but to/from a frame with equal-size and opposite velocities instead of ditto three-momenta.

void RotBstMatrix::rotbst(const RotBstMatrix& Min);  
combine the current matrix with another one.

RotBstMatrix& RotBstMatrix::operator*(const RotBstMatrix& Min)  
return the combination of this and another matrix (corresponding to first making the transformation M_{in} and then this transformation.

Vec4 RotBstMatrix::operator*(Vec4 p)  
return a rotated and boosted version of p.

void RotBstMatrix::invert()  
invert the matrix, which corresponds to an opposite sequence and sign of rotations and boosts.

RotBstMatrix RotBstMatrix::inverse()  
return the inverse matrix, which corresponds to an opposite sequence and sign of rotations and boosts, without overwriting the original matrix.

void RotBstMatrix::reset()  
reset to no rotation/boost; i.e. the default at creation.

double RotBstMatrix::value(int i, int j)  
return the value of the (i,j) = [i][j] matrix element.

double RotBstMatrix::deviation()  
crude estimate how much a matrix deviates from the unit matrix: the sum of the absolute values of all non-diagonal matrix elements plus the sum of the absolute deviation of the diagonal matrix elements from unity.

friend ostream& operator<<(ostream&, const RotBstMatrix& M)  
writes out the values of the sixteen components of a RotBstMatrix, on four consecutive lines and ended with a newline.

RotBstMatrix toCMframe(const Vec4& p)  
Return a RotBstMatrix corresponding to a boost to the rest frame of p.

RotBstMatrix fromCMframe(const Vec4& p)  
Return a RotBstMatrix corresponding to a boost from the rest frame of p .

RotBstMatrix toCMframe(const Vec4& p1, const Vec4& p2)  
Return a RotBstMatrix corresponding to a boost and rotation to the rest frame of p_1 and p_2, with p_1 along the +z axis.

RotBstMatrix fromCMframe(const Vec4& p1, const Vec4& p2, bool flip = false)  
Return a RotBstMatrix corresponding to a boost and rotation from the rest frame of p_1 and p_2, with p_1 along the +z axis. The option flip handles the case when p_1 is along the -z axis in the rest frame. This is accomplished by performing the rotation for p_2 and negating the rotational part.

RotBstMatrix toCMframe(const Vec4& ptot, const Vec4& pz, const Vec4& pxz)  
Return a RotBstMatrix corresponding to a boost and rotation to the rest frame of p_{tot}, where additionally the p_z vector is put along the +z axis and the p_{xz} on put in the xz plane with positive x.

RotBstMatrix fromCMframe(const Vec4& ptot, const Vec4& pz, const Vec4& pxz)  
Return a RotBstMatrix corresponding to a boost and rotation from the rest frame of p_{tot}, where additionally the p_z vector is along the +z axis and the p_{xz} is in the xz plane with positive x when in the rest frame.