Logo ROOT  
Reference Guide
LorentzRotation.cxx
Go to the documentation of this file.
1// @(#)root/mathcore:$Id$
2// Authors: W. Brown, M. Fischler, L. Moneta 2005
3
4 /**********************************************************************
5 * *
6 * Copyright (c) 2005 , LCG ROOT FNAL MathLib Team *
7 * *
8 * *
9 **********************************************************************/
10
11// Header file for class LorentzRotation, a 4x4 matrix representation of
12// a general Lorentz transformation
13//
14// Created by: Mark Fischler Mon Aug 8 2005
15//
16
18
23
24#include <cmath>
25#include <algorithm>
26
31
32namespace ROOT {
33
34namespace Math {
35
37 // constructor of an identity LR
38 fM[kXX] = 1.0; fM[kXY] = 0.0; fM[kXZ] = 0.0; fM[kXT] = 0.0;
39 fM[kYX] = 0.0; fM[kYY] = 1.0; fM[kYZ] = 0.0; fM[kYT] = 0.0;
40 fM[kZX] = 0.0; fM[kZY] = 0.0; fM[kZZ] = 1.0; fM[kZT] = 0.0;
41 fM[kTX] = 0.0; fM[kTY] = 0.0; fM[kTZ] = 0.0; fM[kTT] = 1.0;
42}
43
45 // construct from Rotation3D
46 r.GetComponents ( fM[kXX], fM[kXY], fM[kXZ],
47 fM[kYX], fM[kYY], fM[kYZ],
48 fM[kZX], fM[kZY], fM[kZZ] );
49 fM[kXT] = 0.0;
50 fM[kYT] = 0.0;
51 fM[kZT] = 0.0;
52 fM[kTX] = 0.0; fM[kTY] = 0.0; fM[kTZ] = 0.0; fM[kTT] = 1.0;
53}
54
56 // construct from AxisAngle
57 const Rotation3D r(a);
58 r.GetComponents ( fM[kXX], fM[kXY], fM[kXZ],
59 fM[kYX], fM[kYY], fM[kYZ],
60 fM[kZX], fM[kZY], fM[kZZ] );
61 fM[kXT] = 0.0;
62 fM[kYT] = 0.0;
63 fM[kZT] = 0.0;
64 fM[kTX] = 0.0; fM[kTY] = 0.0; fM[kTZ] = 0.0; fM[kTT] = 1.0;
65}
66
68 // construct from EulerAngles
69 const Rotation3D r(e);
70 r.GetComponents ( fM[kXX], fM[kXY], fM[kXZ],
71 fM[kYX], fM[kYY], fM[kYZ],
72 fM[kZX], fM[kZY], fM[kZZ] );
73 fM[kXT] = 0.0;
74 fM[kYT] = 0.0;
75 fM[kZT] = 0.0;
76 fM[kTX] = 0.0; fM[kTY] = 0.0; fM[kTZ] = 0.0; fM[kTT] = 1.0;
77}
78
80 // construct from Quaternion
81 const Rotation3D r(q);
82 r.GetComponents ( fM[kXX], fM[kXY], fM[kXZ],
83 fM[kYX], fM[kYY], fM[kYZ],
84 fM[kZX], fM[kZY], fM[kZZ] );
85 fM[kXT] = 0.0;
86 fM[kYT] = 0.0;
87 fM[kZT] = 0.0;
88 fM[kTX] = 0.0; fM[kTY] = 0.0; fM[kTZ] = 0.0; fM[kTT] = 1.0;
89}
90
92 // construct from RotationX
93 Scalar s = r.SinAngle();
94 Scalar c = r.CosAngle();
95 fM[kXX] = 1.0; fM[kXY] = 0.0; fM[kXZ] = 0.0; fM[kXT] = 0.0;
96 fM[kYX] = 0.0; fM[kYY] = c ; fM[kYZ] = -s ; fM[kYT] = 0.0;
97 fM[kZX] = 0.0; fM[kZY] = s ; fM[kZZ] = c ; fM[kZT] = 0.0;
98 fM[kTX] = 0.0; fM[kTY] = 0.0; fM[kTZ] = 0.0; fM[kTT] = 1.0;
99}
100
102 // construct from RotationY
103 Scalar s = r.SinAngle();
104 Scalar c = r.CosAngle();
105 fM[kXX] = c ; fM[kXY] = 0.0; fM[kXZ] = s ; fM[kXT] = 0.0;
106 fM[kYX] = 0.0; fM[kYY] = 1.0; fM[kYZ] = 0.0; fM[kYT] = 0.0;
107 fM[kZX] = -s ; fM[kZY] = 0.0; fM[kZZ] = c ; fM[kZT] = 0.0;
108 fM[kTX] = 0.0; fM[kTY] = 0.0; fM[kTZ] = 0.0; fM[kTT] = 1.0;
109}
110
112 // construct from RotationX
113 Scalar s = r.SinAngle();
114 Scalar c = r.CosAngle();
115 fM[kXX] = c ; fM[kXY] = -s ; fM[kXZ] = 0.0; fM[kXT] = 0.0;
116 fM[kYX] = s ; fM[kYY] = c ; fM[kYZ] = 0.0; fM[kYT] = 0.0;
117 fM[kZX] = 0.0; fM[kZY] = 0.0; fM[kZZ] = 1.0; fM[kZT] = 0.0;
118 fM[kTX] = 0.0; fM[kTY] = 0.0; fM[kTZ] = 0.0; fM[kTT] = 1.0;
119}
120
121void
123 // Assuming the representation of this is close to a true Lorentz Rotation,
124 // but may have drifted due to round-off error from many operations,
125 // this forms an "exact" orthosymplectic matrix for the Lorentz Rotation
126 // again.
127
128 typedef LorentzVector< PxPyPzE4D<Scalar> > FourVector;
129 if (fM[kTT] <= 0) {
131 "LorentzRotation:Rectify(): Non-positive TT component - cannot rectify");
132 return;
133 }
134 FourVector t ( fM[kTX], fM[kTY], fM[kTZ], fM[kTT] );
135 Scalar m2 = t.M2();
136 if ( m2 <= 0 ) {
138 "LorentzRotation:Rectify(): Non-timelike time row - cannot rectify");
139 return;
140 }
141 t /= std::sqrt(m2);
142 FourVector z ( fM[kZX], fM[kZY], fM[kZZ], fM[kZT] );
143 z = z - z.Dot(t)*t;
144 m2 = z.M2();
145 if ( m2 >= 0 ) {
147 "LorentzRotation:Rectify(): Non-spacelike Z row projection - "
148 "cannot rectify");
149 return;
150 }
151 z /= std::sqrt(-m2);
152 FourVector y ( fM[kYX], fM[kYY], fM[kYZ], fM[kYT] );
153 y = y - y.Dot(t)*t - y.Dot(z)*z;
154 m2 = y.M2();
155 if ( m2 >= 0 ) {
157 "LorentzRotation:Rectify(): Non-spacelike Y row projection - "
158 "cannot rectify");
159 return;
160 }
161 y /= std::sqrt(-m2);
162 FourVector x ( fM[kXX], fM[kXY], fM[kXZ], fM[kXT] );
163 x = x - x.Dot(t)*t - x.Dot(z)*z - x.Dot(y)*y;
164 m2 = x.M2();
165 if ( m2 >= 0 ) {
167 "LorentzRotation:Rectify(): Non-spacelike X row projection - "
168 "cannot rectify");
169 return;
170 }
171 x /= std::sqrt(-m2);
172}
173
174
176 // invert modifying current content
177 Scalar temp;
178 temp = fM[kXY]; fM[kXY] = fM[kYX]; fM[kYX] = temp;
179 temp = fM[kXZ]; fM[kXZ] = fM[kZX]; fM[kZX] = temp;
180 temp = fM[kYZ]; fM[kYZ] = fM[kZY]; fM[kZY] = temp;
181 temp = fM[kXT]; fM[kXT] = -fM[kTX]; fM[kTX] = -temp;
182 temp = fM[kYT]; fM[kYT] = -fM[kTY]; fM[kTY] = -temp;
183 temp = fM[kZT]; fM[kZT] = -fM[kTZ]; fM[kTZ] = -temp;
184}
185
187 // return an inverse LR
188 return LorentzRotation
189 ( fM[kXX], fM[kYX], fM[kZX], -fM[kTX]
190 , fM[kXY], fM[kYY], fM[kZY], -fM[kTY]
191 , fM[kXZ], fM[kYZ], fM[kZZ], -fM[kTZ]
192 , -fM[kXT], -fM[kYT], -fM[kZT], fM[kTT]
193 );
194}
195
197 // combination with another LR
198 return LorentzRotation
199 ( fM[kXX]*r.fM[kXX] + fM[kXY]*r.fM[kYX] + fM[kXZ]*r.fM[kZX] + fM[kXT]*r.fM[kTX]
200 , fM[kXX]*r.fM[kXY] + fM[kXY]*r.fM[kYY] + fM[kXZ]*r.fM[kZY] + fM[kXT]*r.fM[kTY]
201 , fM[kXX]*r.fM[kXZ] + fM[kXY]*r.fM[kYZ] + fM[kXZ]*r.fM[kZZ] + fM[kXT]*r.fM[kTZ]
202 , fM[kXX]*r.fM[kXT] + fM[kXY]*r.fM[kYT] + fM[kXZ]*r.fM[kZT] + fM[kXT]*r.fM[kTT]
203 , fM[kYX]*r.fM[kXX] + fM[kYY]*r.fM[kYX] + fM[kYZ]*r.fM[kZX] + fM[kYT]*r.fM[kTX]
204 , fM[kYX]*r.fM[kXY] + fM[kYY]*r.fM[kYY] + fM[kYZ]*r.fM[kZY] + fM[kYT]*r.fM[kTY]
205 , fM[kYX]*r.fM[kXZ] + fM[kYY]*r.fM[kYZ] + fM[kYZ]*r.fM[kZZ] + fM[kYT]*r.fM[kTZ]
206 , fM[kYX]*r.fM[kXT] + fM[kYY]*r.fM[kYT] + fM[kYZ]*r.fM[kZT] + fM[kYT]*r.fM[kTT]
207 , fM[kZX]*r.fM[kXX] + fM[kZY]*r.fM[kYX] + fM[kZZ]*r.fM[kZX] + fM[kZT]*r.fM[kTX]
208 , fM[kZX]*r.fM[kXY] + fM[kZY]*r.fM[kYY] + fM[kZZ]*r.fM[kZY] + fM[kZT]*r.fM[kTY]
209 , fM[kZX]*r.fM[kXZ] + fM[kZY]*r.fM[kYZ] + fM[kZZ]*r.fM[kZZ] + fM[kZT]*r.fM[kTZ]
210 , fM[kZX]*r.fM[kXT] + fM[kZY]*r.fM[kYT] + fM[kZZ]*r.fM[kZT] + fM[kZT]*r.fM[kTT]
211 , fM[kTX]*r.fM[kXX] + fM[kTY]*r.fM[kYX] + fM[kTZ]*r.fM[kZX] + fM[kTT]*r.fM[kTX]
212 , fM[kTX]*r.fM[kXY] + fM[kTY]*r.fM[kYY] + fM[kTZ]*r.fM[kZY] + fM[kTT]*r.fM[kTY]
213 , fM[kTX]*r.fM[kXZ] + fM[kTY]*r.fM[kYZ] + fM[kTZ]*r.fM[kZZ] + fM[kTT]*r.fM[kTZ]
214 , fM[kTX]*r.fM[kXT] + fM[kTY]*r.fM[kYT] + fM[kTZ]*r.fM[kZT] + fM[kTT]*r.fM[kTT]
215 );
216}
217
218
219std::ostream & operator<< (std::ostream & os, const LorentzRotation & r) {
220 // TODO - this will need changing for machine-readable issues
221 // and even the human readable form needs formatiing improvements
222 double m[16];
223 r.GetComponents(m, m+16);
224 os << "\n" << m[0] << " " << m[1] << " " << m[2] << " " << m[3];
225 os << "\n" << m[4] << " " << m[5] << " " << m[6] << " " << m[7];
226 os << "\n" << m[8] << " " << m[9] << " " << m[10] << " " << m[11];
227 os << "\n" << m[12] << " " << m[13] << " " << m[14] << " " << m[15] << "\n";
228 return os;
229}
230
231
232} //namespace Math
233} //namespace ROOT
ROOT::R::TRInterface & r
Definition: Object.C:4
#define c(i)
Definition: RSha256.hxx:101
#define e(i)
Definition: RSha256.hxx:103
float * q
Definition: THbookFile.cxx:87
double sqrt(double)
AxisAngle class describing rotation represented with direction axis (3D Vector) and an angle of rotat...
Definition: AxisAngle.h:41
EulerAngles class describing rotation as three angles (Euler Angles).
Definition: EulerAngles.h:43
Lorentz transformation class with the (4D) transformation represented by a 4x4 orthosymplectic matrix...
A4Vector operator*(const A4Vector &v) const
Overload operator * for rotation on a vector.
LorentzRotation Inverse() const
Return inverse of a rotation.
void Rectify()
Re-adjust components to eliminate small deviations from a perfect orthosyplectic matrix.
void Invert()
Invert a Lorentz rotation in place.
LorentzRotation()
Default constructor (identity transformation)
Rotation class with the (3D) rotation represented by a unit quaternion (u, i, j, k).
Definition: Quaternion.h:47
Rotation class with the (3D) rotation represented by a 3x3 orthogonal matrix.
Definition: Rotation3D.h:65
Rotation class representing a 3D rotation about the X axis by the angle of rotation.
Definition: RotationX.h:43
Rotation class representing a 3D rotation about the Y axis by the angle of rotation.
Definition: RotationY.h:43
Rotation class representing a 3D rotation about the Z axis by the angle of rotation.
Definition: RotationZ.h:43
Double_t y[n]
Definition: legend1.C:17
Double_t x[n]
Definition: legend1.C:17
Namespace for new Math classes and functions.
void Throw(const char *)
function throwing exception, by creating internally a GenVector_exception only when needed
std::ostream & operator<<(std::ostream &os, const AxisAngle &a)
Stream Output and Input.
Definition: AxisAngle.cxx:91
VSD Structures.
Definition: StringConv.hxx:21
static constexpr double s
static constexpr double m2
auto * m
Definition: textangle.C:8
auto * a
Definition: textangle.C:12