Logo ROOT   6.30.04
Reference Guide
 All Namespaces Files Pages
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 
17 #include "Math/GenVector/GenVectorIO.h"
18 
19 #include "Math/GenVector/LorentzRotation.h"
20 #include "Math/GenVector/LorentzVector.h"
21 #include "Math/GenVector/PxPyPzE4D.h"
22 #include "Math/GenVector/GenVector_exception.h"
23 
24 #include <cmath>
25 #include <algorithm>
26 
27 #include "Math/GenVector/Rotation3D.h"
28 #include "Math/GenVector/RotationX.h"
29 #include "Math/GenVector/RotationY.h"
30 #include "Math/GenVector/RotationZ.h"
31 
32 namespace ROOT {
33 
34 namespace Math {
35 
36 LorentzRotation::LorentzRotation() {
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 
44 LorentzRotation::LorentzRotation(Rotation3D const & r) {
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 
55 LorentzRotation::LorentzRotation(AxisAngle const & a) {
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 
67 LorentzRotation::LorentzRotation(EulerAngles const & e) {
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 
79 LorentzRotation::LorentzRotation(Quaternion const & q) {
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 
91 LorentzRotation::LorentzRotation(RotationX const & r) {
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 
101 LorentzRotation::LorentzRotation(RotationY const & r) {
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 
111 LorentzRotation::LorentzRotation(RotationZ const & r) {
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 
121 void
122 LorentzRotation::Rectify() {
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) {
130  GenVector::Throw (
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 ) {
137  GenVector::Throw (
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 ) {
146  GenVector::Throw (
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 ) {
156  GenVector::Throw (
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 ) {
166  GenVector::Throw (
167  "LorentzRotation:Rectify(): Non-spacelike X row projection - "
168  "cannot rectify");
169  return;
170  }
171  x /= std::sqrt(-m2);
172 }
173 
174 
175 void LorentzRotation::Invert() {
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 
186 LorentzRotation LorentzRotation::Inverse() const {
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 
196 LorentzRotation LorentzRotation::operator * (const LorentzRotation & r) const {
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 
219 std::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
GenVector_exception.h