3
<title>Four-Vectors</title>
4
<link rel="stylesheet" type="text/css" href="pythia.css"/>
5
<link rel="shortcut icon" href="pythia32.gif"/>
11
The <code>Vec4</code> class gives a simple implementation of four-vectors.
12
The member function names are based on the assumption that these
13
represent four-momentum vectors. Thus one can get or set
14
<i>p_x, p_y, p_z</i> and <i>e</i>, but not <i>x, y, z</i>
15
or <i>t</i>. This is only a matter of naming, however; a
16
<code>Vec4</code> can equally well be used to store a space-time
20
The <code>Particle</code> object contains a <code>Vec4 p</code> that
21
stores the particle four-momentum, and another <code>Vec4 vProd</code>
22
for the production vertex. For the latter the input/output method
23
names are adapted to the space-time character rather than the normal
24
energy-momentum one. Thus a user would not normally access the
25
<code>Vec4</code> classes directly, but only via the methods of the
26
<code>Particle</code> class,
27
see <a href="ParticleProperties.html" target="page">Particle Properties</a>.
30
Nevertheless you are free to use the PYTHIA four-vectors, e.g. as
31
part of some simple analysis code based directly on the PYTHIA output,
32
say to define the four-vector sum of a set of particles. But note that
33
this class was never set up to allow complete generality, only to
34
provide the operations that are of use inside PYTHIA. There is no
35
separate class for three-vectors, since such can easily be represented
36
by four-vectors where the fourth component is not used.
39
Four-vectors have the expected functionality: they can be created,
40
copied, added, multiplied, rotated, boosted, and manipulated in other
41
ways. Operator overloading is implemented where reasonable. Properties
42
can be read out, not only the components themselves but also for derived
43
quantities such as absolute momentum and direction angles.
45
<h3>Constructors and basic operators</h3>
47
A few methods are available to create or copy a four-vector:
49
<a name="method1"></a>
50
<p/><strong>Vec4::Vec4(double x = 0., double y = 0., double z = 0., double t = 0.) </strong> <br/>
51
creates a four-vector, by default with all components set to 0.
54
<a name="method2"></a>
55
<p/><strong>Vec4::Vec4(const Vec4& v) </strong> <br/>
56
creates a four-vector copy of the input four-vector.
59
<a name="method3"></a>
60
<p/><strong>Vec4& Vec4::operator=(const Vec4& v) </strong> <br/>
61
copies the input four-vector.
64
<a name="method4"></a>
65
<p/><strong>Vec4& Vec4::operator=(double value) </strong> <br/>
66
gives a four-vector with all components set to <i>value</i>.
69
<h3>Member methods for input</h3>
71
The values stored in a four-vector can be modified in a few different
74
<a name="method5"></a>
75
<p/><strong>void Vec4::reset() </strong> <br/>
76
sets all components to 0.
79
<a name="method6"></a>
80
<p/><strong>void Vec4::p(double pxIn, double pyIn, double pzIn, double eIn) </strong> <br/>
81
sets all components to their input values.
84
<a name="method7"></a>
85
<p/><strong>void Vec4::p(Vec4 pIn) </strong> <br/>
86
sets all components equal to those of the input four-vector.
89
<a name="method8"></a>
90
<p/><strong>void Vec4::px(double pxIn) </strong> <br/>
92
<strong>void Vec4::py(double pyIn) </strong> <br/>
94
<strong>void Vec4::pz(double pzIn) </strong> <br/>
96
<strong>void Vec4::e(double eIn) </strong> <br/>
97
sets the respective component to the input value.
100
<h3>Member methods for output</h3>
102
A number of methods provides output of basic or derived quantities:
104
<a name="method9"></a>
105
<p/><strong>double Vec4::px() </strong> <br/>
107
<strong>double Vec4::py() </strong> <br/>
109
<strong>double Vec4::pz() </strong> <br/>
111
<strong>double Vec4::e() </strong> <br/>
112
gets the respective component.
115
<a name="method10"></a>
116
<p/><strong>double& operator[](int i) </strong> <br/>
117
returns component by index, where 1 gives <i>p_x</i>, 2 gives <i>p_y</i>,
118
3 gives <i>p_z</i>, and anything else gives <i>e</i>.
121
<a name="method11"></a>
122
<p/><strong>double Vec4::mCalc() </strong> <br/>
124
<strong>double Vec4::m2Calc() </strong> <br/>
125
the (squared) mass, calculated from the four-vectors.
126
If <i>m^2 < 0</i> the mass is given with a minus sign,
127
<i>-sqrt(-m^2)</i>. Note the possible loss of precision
128
in the calculation of <i>E^2 - p^2</i>; for particles the
129
correct mass is stored separately to avoid such problems.
132
<a name="method12"></a>
133
<p/><strong>double Vec4::pT() </strong> <br/>
135
<strong>double Vec4::pT2() </strong> <br/>
136
the (squared) transverse momentum.
139
<a name="method13"></a>
140
<p/><strong>double Vec4::pAbs() </strong> <br/>
142
<strong>double Vec4::pAbs2() </strong> <br/>
143
the (squared) absolute momentum.
146
<a name="method14"></a>
147
<p/><strong>double Vec4::eT() </strong> <br/>
149
<strong>double Vec4::eT2() </strong> <br/>
150
the (squared) transverse energy,
151
<i>eT = e * sin(theta) = e * pT / pAbs</i>.
154
<a name="method15"></a>
155
<p/><strong>double Vec4::theta() </strong> <br/>
156
the polar angle, in the range 0 through
160
<a name="method16"></a>
161
<p/><strong>double Vec4::phi() </strong> <br/>
162
the azimuthal angle, in the range <i>-pi</i> through <i>pi</i>.
165
<a name="method17"></a>
166
<p/><strong>double Vec4::thetaXZ() </strong> <br/>
167
the angle in the <i>xz</i> plane, in the range <i>-pi</i> through
168
<i>pi</i>, with 0 along the <i>+z</i> axis.
171
<a name="method18"></a>
172
<p/><strong>double Vec4::pPos() </strong> <br/>
174
<strong>double Vec4::pNeg() </strong> <br/>
175
the combinations <i>E+-p_z</i>.
177
<a name="method19"></a>
178
<p/><strong>double Vec4::rap() </strong> <br/>
180
<strong>double Vec4::eta() </strong> <br/>
181
true rapidity <i>y</i> and pseudorapidity <i>eta</i>.
184
<h3>Friend methods for output</h3>
186
There are also some <code>friend</code> methods that take one, two
187
or three four-vectors as argument. Several of them only use the
188
three-vector part of the four-vector.
190
<a name="method20"></a>
191
<p/><strong>friend ostream& operator<<(ostream&, const Vec4& v) </strong> <br/>
192
writes out the values of the four components of a <code>Vec4</code> and,
193
within brackets, a fifth component being the invariant length of the
194
four-vector, as provided by <code>mCalc()</code> above, and it all
195
ended with a newline.
198
<a name="method21"></a>
199
<p/><strong>friend double m(const Vec4& v1, const Vec4& v2) </strong> <br/>
201
<strong>friend double m2(const Vec4& v1, const Vec4& v2) </strong> <br/>
202
the (squared) invariant mass.
205
<a name="method22"></a>
206
<p/><strong>friend double dot3(const Vec4& v1, const Vec4& v2) </strong> <br/>
210
<a name="method23"></a>
211
<p/><strong>friend double cross3(const Vec4& v1, const Vec4& v2) </strong> <br/>
215
<a name="method24"></a>
216
<p/><strong>friend double theta(const Vec4& v1, const Vec4& v2) </strong> <br/>
218
<strong>friend double costheta(const Vec4& v1, const Vec4& v2) </strong> <br/>
219
the (cosine) of the opening angle between the vectors,
220
in the range 0 through <i>pi</i>.
223
<a name="method25"></a>
224
<p/><strong>friend double phi(const Vec4& v1, const Vec4& v2) </strong> <br/>
226
<strong>friend double cosphi(const Vec4& v1, const Vec4& v2) </strong> <br/>
227
the (cosine) of the azimuthal angle between the vectors around the
228
<i>z</i> axis, in the range 0 through <i>pi</i>.
231
<a name="method26"></a>
232
<p/><strong>friend double phi(const Vec4& v1, const Vec4& v2, const Vec4& v3) </strong> <br/>
234
<strong>friend double cosphi(const Vec4& v1, const Vec4& v2, const Vec4& v3) </strong> <br/>
235
the (cosine) of the azimuthal angle between the first two vectors
236
around the direction of the third, in the range 0 through <i>pi</i>.
239
<a name="method27"></a>
240
<p/><strong>friend double RRapPhi(const Vec4& v1, const Vec4& v2) </strong> <br/>
242
<strong>friend double REtaPhi(const Vec4& v1, const Vec4& v2) </strong> <br/>
243
the <i>R</i> distance measure, in <i>(y, phi)</i> or
244
<i>(eta, phi)</i> cylindrical coordinates, i.e.
245
<i>R^2 = (y_1 - y_2)^2 + (phi_1 - phi_2)^2</i> and equivalent.
248
<a name="method28"></a>
249
<p/><strong>friend bool pShift( Vec4& p1Move, Vec4& p2Move, double m1New, double m2New) </strong> <br/>
250
transfer four-momentum between the two four-vectors so that they get
251
the masses <code>m1New</code> and <code>m2New</code>, respectively.
252
Note that <code>p1Move</code> and <code>p2Move</code> act both as
253
input and output arguments. The method will return false if the invariant
254
mass of the four-vectors is too small to accommodate the new masses,
255
and then the four-vectors are not changed.
258
<h3>Operations with four-vectors</h3>
260
Of course one should be able to add, subtract and scale four-vectors,
263
<a name="method29"></a>
264
<p/><strong>Vec4 Vec4::operator-() </strong> <br/>
265
return a vector with flipped sign for all components, while leaving
266
the original vector unchanged.
269
<a name="method30"></a>
270
<p/><strong>Vec4& Vec4::operator+=(const Vec4& v) </strong> <br/>
271
add a four-vector to an existing one.
274
<a name="method31"></a>
275
<p/><strong>Vec4& Vec4::operator-=(const Vec4& v) </strong> <br/>
276
subtract a four-vector from an existing one.
279
<a name="method32"></a>
280
<p/><strong>Vec4& Vec4::operator*=(double f) </strong> <br/>
281
multiply all four-vector components by a real number.
284
<a name="method33"></a>
285
<p/><strong>Vec4& Vec4::operator/=(double f) </strong> <br/>
286
divide all four-vector components by a real number.
289
<a name="method34"></a>
290
<p/><strong>friend Vec4 operator+(const Vec4& v1, const Vec4& v2) </strong> <br/>
291
add two four-vectors.
294
<a name="method35"></a>
295
<p/><strong>friend Vec4 operator-(const Vec4& v1, const Vec4& v2) </strong> <br/>
296
subtract two four-vectors.
299
<a name="method36"></a>
300
<p/><strong>friend Vec4 operator*(double f, const Vec4& v) </strong> <br/>
302
<strong>friend Vec4 operator*(const Vec4& v, double f) </strong> <br/>
303
multiply a four-vector by a real number.
306
<a name="method37"></a>
307
<p/><strong>friend Vec4 operator/(const Vec4& v, double f) </strong> <br/>
308
divide a four-vector by a real number.
311
<a name="method38"></a>
312
<p/><strong>friend double operator*(const Vec4& v1, const Vec4 v2) </strong> <br/>
317
There are also a few related operations that are normal member methods:
319
<a name="method39"></a>
320
<p/><strong>void Vec4::rescale3(double f) </strong> <br/>
321
multiply the three-vector components by <i>f</i>, but keep the
322
fourth component unchanged.
325
<a name="method40"></a>
326
<p/><strong>void Vec4::rescale4(double f) </strong> <br/>
327
multiply all four-vector components by <i>f</i>.
330
<a name="method41"></a>
331
<p/><strong>void Vec4::flip3() </strong> <br/>
332
flip the sign of the three-vector components, but keep the
333
fourth component unchanged.
336
<a name="method42"></a>
337
<p/><strong>void Vec4::flip4() </strong> <br/>
338
flip the sign of all four-vector components.
341
<h3>Rotations and boosts</h3>
343
A common task is to rotate or boost four-vectors. In case only one
344
four-vector is affected the operation may be performed directly on it.
345
However, in case many particles are affected, the helper class
346
<code>RotBstMatrix</code> can be used to speed up operations.
348
<a name="method43"></a>
349
<p/><strong>void Vec4::rot(double theta, double phi) </strong> <br/>
350
rotate the three-momentum with the polar angle <i>theta</i>
351
and the azimuthal angle <i>phi</i>.
354
<a name="method44"></a>
355
<p/><strong>void Vec4::rotaxis(double phi, double nx, double ny, double nz) </strong> <br/>
356
rotate the three-momentum with the azimuthal angle <i>phi</i>
357
around the direction defined by the <i>(n_x, n_y, n_z)</i>
361
<a name="method45"></a>
362
<p/><strong>void Vec4::rotaxis(double phi, Vec4& n) </strong> <br/>
363
rotate the three-momentum with the azimuthal angle <i>phi</i>
364
around the direction defined by the three-vector part of <i>n</i>.
367
<a name="method46"></a>
368
<p/><strong>void Vec4::bst(double betaX, double betaY, double betaZ) </strong> <br/>
369
boost the four-momentum by <i>beta = (beta_x, beta_y, beta_z)</i>.
372
<a name="method47"></a>
373
<p/><strong>void Vec4::bst(double betaX, double betaY, double betaZ, double gamma) </strong> <br/>
374
boost the four-momentum by <i>beta = (beta_x, beta_y, beta_z)</i>,
375
where the <i>gamma = 1/sqrt(1 - beta^2)</i> is also input to allow
376
better precision when <i>beta</i> is close to unity.
379
<a name="method48"></a>
380
<p/><strong>void Vec4::bst(const Vec4& p) </strong> <br/>
381
boost the four-momentum by <i>beta = (p_x/E, p_y/E, p_z/E)</i>.
384
<a name="method49"></a>
385
<p/><strong>void Vec4::bst(const Vec4& p, double m) </strong> <br/>
386
boost the four-momentum by <i>beta = (p_x/E, p_y/E, p_z/E)</i>,
387
where the <i>gamma = E/m</i> is also calculated from input to allow
388
better precision when <i>beta</i> is close to unity.
391
<a name="method50"></a>
392
<p/><strong>void Vec4::bstback(const Vec4& p) </strong> <br/>
393
boost the four-momentum by <i>beta = (-p_x/E, -p_y/E, -p_z/E)</i>.
396
<a name="method51"></a>
397
<p/><strong>void Vec4::bstback(const Vec4& p, double m) </strong> <br/>
398
boost the four-momentum by <i>beta = (-p_x/E, -p_y/E, -p_z/E)</i>,
399
where the <i>gamma = E/m</i> is also calculated from input to allow
400
better precision when <i>beta</i> is close to unity.
403
<a name="method52"></a>
404
<p/><strong>void Vec4::rotbst(const RotBstMatrix& M) </strong> <br/>
405
perform a combined rotation and boost; see below for a description
406
of the <code>RotBstMatrix</code>.
410
For a longer sequence of rotations and boosts, and where several
411
<code>Vec4</code> are to be rotated and boosted in the same way,
412
a more efficient approach is to define a <code>RotBstMatrix</code>,
413
which forms a separate auxiliary class. You can build up this
414
4-by-4 matrix by successive calls to the methods of the class,
415
such that the matrix encodes the full sequence of operations.
416
The order in which you do these calls must agree with the imagined
417
order in which the rotations/boosts should be applied to a
418
four-momentum, since in general the operations do not commute.
419
<br/>(Mathematically you would e.g. define <i>M = M_3 M_2 M_1</i>
420
in that <i>M p = M_3( M_2( M_1 p) ) )</i>. That is, operations
421
on the four-vector <i>p</i> are carried out in the order first
422
<i>M_1</i>, then <i>M_2</i> and finally <i>M_3</i>. Thus
423
<i>M_1, M_2, M_3</i> is also the order in which you should input
424
rotations and boosts to <i>M</i>.)
426
<a name="method53"></a>
427
<p/><strong>RotBstMatrix::RotBstMatrix() </strong> <br/>
428
creates a diagonal unit matrix, i.e. one that leaves a four-vector
432
<a name="method54"></a>
433
<p/><strong>RotBstMatrix::RotBstMatrix(const RotBstMatrix& Min) </strong> <br/>
434
creates a copy of the input matrix.
437
<a name="method55"></a>
438
<p/><strong>RotBstMatrix& RotBstMatrix::operator=(const RotBstMatrix4& Min) </strong> <br/>
439
copies the input matrix.
442
<a name="method56"></a>
443
<p/><strong>void RotBstMatrix::rot(double theta = 0., double phi = 0.) </strong> <br/>
444
rotate by this polar and azimuthal angle.
447
<a name="method57"></a>
448
<p/><strong>void RotBstMatrix::rot(const Vec4& p) </strong> <br/>
449
rotate so that a vector originally along the <i>+z</i> axis becomes
450
parallel with <i>p</i>. More specifically, rotate by <i>-phi</i>,
451
<i>theta</i> and <i>phi</i>, with angles defined by <i>p</i>.
454
<a name="method58"></a>
455
<p/><strong>void RotBstMatrix::bst(double betaX = 0., double betaY = 0., double betaZ = 0.) </strong> <br/>
456
boost by this <i>beta</i> vector.
459
<a name="method59"></a>
460
<p/><strong>void RotBstMatrix::bst(const Vec4&) </strong> <br/>
462
<strong>void RotBstMatrix::bstback(const Vec4&) </strong> <br/>
463
boost with a <i>beta = p/E</i> or <i>beta = -p/E</i>, respectively.
466
<a name="method60"></a>
467
<p/><strong>void RotBstMatrix::bst(const Vec4& p1, const Vec4& p2) </strong> <br/>
468
boost so that <i>p_1</i> is transformed to <i>p_2</i>. It is assumed
469
that the two vectors obey <i>p_1^2 = p_2^2</i>.
472
<a name="method61"></a>
473
<p/><strong>void RotBstMatrix::toCMframe(const Vec4& p1, const Vec4& p2) </strong> <br/>
474
boost and rotate to the rest frame of <i>p_1</i> and <i>p_2</i>,
475
with <i>p_1</i> along the <i>+z</i> axis.
478
<a name="method62"></a>
479
<p/><strong>void RotBstMatrix::fromCMframe(const Vec4& p1, const Vec4& p2) </strong> <br/>
480
rotate and boost from the rest frame of <i>p_1</i> and <i>p_2</i>,
481
with <i>p_1</i> along the <i>+z</i> axis, to the actual frame of
482
<i>p_1</i> and <i>p_2</i>, i.e. the inverse of the above.
485
<a name="method63"></a>
486
<p/><strong>void RotBstMatrix::rotbst(const RotBstMatrix& Min); </strong> <br/>
487
combine the current matrix with another one.
490
<a name="method64"></a>
491
<p/><strong>void RotBstMatrix::invert() </strong> <br/>
492
invert the matrix, which corresponds to an opposite sequence and sign
493
of rotations and boosts.
496
<a name="method65"></a>
497
<p/><strong>void RotBstMatrix::reset() </strong> <br/>
498
reset to no rotation/boost; i.e. the default at creation.
501
<a name="method66"></a>
502
<p/><strong>double RotBstMatrix::deviation() </strong> <br/>
503
crude estimate how much a matrix deviates from the unit matrix:
504
the sum of the absolute values of all non-diagonal matrix elements
505
plus the sum of the absolute deviation of the diagonal matrix
509
<a name="method67"></a>
510
<p/><strong>friend ostream& operator<<(ostream&, const RotBstMatrix& M) </strong> <br/>
511
writes out the values of the sixteen components of a
512
<code>RotBstMatrix</code>, on four consecutive lines and
513
ended with a newline.
519
<!-- Copyright (C) 2017 Torbjorn Sjostrand -->