Rivet  1.8.3
LorentzTrans.hh
1 #ifndef RIVET_MATH_LORENTZTRANS
2 #define RIVET_MATH_LORENTZTRANS
3 
4 #include <iostream>
5 
6 #include "Rivet/Math/MathHeader.hh"
7 #include "Rivet/Math/MathUtils.hh"
8 #include "Rivet/Math/MatrixN.hh"
9 #include "Rivet/Math/Matrix3.hh"
10 #include "Rivet/Math/Vector4.hh"
11 
12 namespace Rivet {
13 
14 
15  inline double lorentzGamma(const double beta) {
16  return 1.0 / sqrt(1 - beta*beta);
17  }
18 
19 
22  friend string toString(const LorentzTransform& lt);
23 
24  public:
26  _boostMatrix = Matrix<4>::mkIdentity();
27  }
28 
30  setBoost(boost);
31  }
32 
33  LorentzTransform(const double betaX, const double betaY, const double betaZ) {
34  setBoost(betaX, betaY, betaZ);
35  }
36 
37  LorentzTransform& setBoost(const Vector3& boost) {
38  assert(boost.mod2() < 1);
39  const double beta = boost.mod();
40  const double gamma = lorentzGamma(beta);
41  _boostMatrix = Matrix<4>::mkIdentity();
42  _boostMatrix.set(0, 0, gamma);
43  _boostMatrix.set(1, 1, gamma);
44  // Positive coeff since these are active boosts
45  _boostMatrix.set(0, 1, +beta*gamma);
46  _boostMatrix.set(1, 0, +beta*gamma);
47  _boostMatrix = rotate(Vector3::mkX(), boost)._boostMatrix;
48  return *this;
49  }
50 
51  // LorentzTransform& setBoost(const Vector3& boostdirection, const double beta) {
52  // const Vector3 boost = boostdirection.unit() * beta;
53  // return setBoost(boost);
54  // }
55 
56  LorentzTransform& setBoost(const double betaX, const double betaY, const double betaZ) {
57  return setBoost(Vector3(betaX, betaY, betaZ));
58  }
59 
60  Vector3 boost() const {
61  FourMomentum boost(_boostMatrix.getColumn(0));
62  //cout << "!!!" << boost << endl;
63  if (boost.isZero()) return boost;
64  assert(boost.E() > 0);
65  const double beta = boost.p().mod() / boost.E();
66  return boost.p().unit() * beta;
67  }
68 
69  double beta() const {
70  return boost().mod();
71  }
72 
73  double gamma() const {
74  return lorentzGamma(beta());
75  }
76 
77  LorentzTransform rotate(const Vector3& from, const Vector3& to) const {
78  return rotate(Matrix3(from, to));
79  }
80 
81  LorentzTransform rotate(const Vector3& axis, const double angle) const {
82  return rotate(Matrix3(axis, angle));
83  }
84 
85  LorentzTransform rotate(const Matrix3& rot) const {
86  LorentzTransform lt = *this;
87  const Matrix4 rot4 = mkMatrix4(rot);
88  const Matrix4 newlt = rot4 * _boostMatrix * rot4.inverse();
89  lt._boostMatrix = newlt;
90  return lt;
91  }
92 
93  FourVector transform(const FourVector& v4) const {
94  return multiply(_boostMatrix, v4);
95  }
96 
97  LorentzTransform inverse() const {
98  LorentzTransform rtn;
99  rtn._boostMatrix = _boostMatrix.inverse();
100  return rtn;
101  }
102 
103 
106  LorentzTransform rtn;
107  rtn._boostMatrix = _boostMatrix * lt._boostMatrix;
108  return rtn;
109  }
110 
111  Matrix4 toMatrix() const {
112  return _boostMatrix;
113  }
114 
115 
116  LorentzTransform operator*(const LorentzTransform& lt) const {
117  return combine(lt);
118  }
119 
120  LorentzTransform preMult(const Matrix3 & m3) {
121  _boostMatrix = multiply(mkMatrix4(m3),_boostMatrix);
122  return *this;
123  }
124 
125  LorentzTransform postMult(const Matrix3 & m3) {
126  _boostMatrix *= mkMatrix4(m3);
127  return *this;
128  }
129 
130  private:
131  Matrix4 mkMatrix4(const Matrix3& m3) const {
132  Matrix4 m4 = Matrix4::mkIdentity();
133  for (size_t i = 0; i < 3; ++i) {
134  for (size_t j = 0; j < 3; ++j) {
135  m4.set(i+1, j+1, m3.get(i, j));
136  }
137  }
138  return m4;
139  }
140 
141 
142  private:
143  Matrix4 _boostMatrix;
144 
145  };
146 
147 
148 
149  inline LorentzTransform inverse(const LorentzTransform& lt) {
150  return lt.inverse();
151  }
152 
153  inline LorentzTransform combine(const LorentzTransform& a, const LorentzTransform& b) {
154  return a.combine(b);
155  }
156 
157  inline FourVector transform(const LorentzTransform& lt, const FourVector& v4) {
158  return lt.transform(v4);
159  }
160 
161 
163 
164 
165  inline string toString(const LorentzTransform& lt) {
166  return toString(lt._boostMatrix);
167  }
168 
169  inline ostream& operator<<(std::ostream& out, const LorentzTransform& lt) {
170  out << toString(lt);
171  return out;
172  }
173 
174 
175 }
176 
177 #endif
Definition: MC_JetAnalysis.hh:9
bool isZero(double tolerance=1E-5) const
Check for nullness, allowing for numerical precision.
Definition: VectorN.hh:67
double mod() const
Calculate the modulus of a vector. .
Definition: VectorN.hh:87
Specialisation of MatrixN to aid 3 dimensional rotations.
Definition: Matrix3.hh:13
LorentzTransform combine(const LorentzTransform &lt) const
Combine LTs, treating this as the LH matrix.
Definition: LorentzTrans.hh:105
Matrix< N > inverse() const
Calculate inverse.
Definition: MatrixN.hh:134
Vector3 unit() const
Definition: Vector3.hh:88
Object implementing Lorentz transform calculations and boosts.
Definition: LorentzTrans.hh:21
double angle(const Vector3 &a, const Vector3 &b)
Angle (in radians) between two 3-vectors.
Definition: Vector3.hh:244
General -dimensional mathematical matrix object.
Definition: MatrixN.hh:14
Specialisation of VectorN to a general (non-momentum) Lorentz 4-vector.
Definition: Vector4.hh:20
double mod2() const
Calculate the modulus-squared of a vector. .
Definition: VectorN.hh:76
Three-dimensional specialisation of Vector.
Definition: Vector3.hh:26
std::ostream & operator<<(std::ostream &os, const AnalysisInfo &ai)
Stream an AnalysisInfo as a text description.
Definition: AnalysisInfo.hh:239
Specialized version of the FourVector with momentum/energy functionality.
Definition: Vector4.hh:324