Logo ROOT   6.30.04
Reference Guide
 All Namespaces Files Pages
LorentzRotation.h
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 ROOT MathLib Team *
7  * *
8  * *
9  **********************************************************************/
10 
11 // Header file for LorentzRotation
12 //
13 // Created by: Mark Fischler Mon Aug 8 2005
14 //
15 // Last update: $Id$
16 //
17 #ifndef ROOT_Math_GenVector_LorentzRotation
18 #define ROOT_Math_GenVector_LorentzRotation 1
19 
20 #include "Math/GenVector/LorentzRotationfwd.h"
21 
22 #include "Math/GenVector/LorentzVector.h"
23 #include "Math/GenVector/PxPyPzE4D.h"
24 
25 #include "Math/GenVector/Rotation3Dfwd.h"
26 #include "Math/GenVector/AxisAnglefwd.h"
27 #include "Math/GenVector/EulerAnglesfwd.h"
28 #include "Math/GenVector/Quaternionfwd.h"
29 #include "Math/GenVector/RotationXfwd.h"
30 #include "Math/GenVector/RotationYfwd.h"
31 #include "Math/GenVector/RotationZfwd.h"
32 #include "Math/GenVector/Boost.h"
33 #include "Math/GenVector/BoostX.h"
34 #include "Math/GenVector/BoostY.h"
35 #include "Math/GenVector/BoostZ.h"
36 
37 namespace ROOT {
38 
39  namespace Math {
40 
41 //__________________________________________________________________________________________
42  /**
43  Lorentz transformation class with the (4D) transformation represented by
44  a 4x4 orthosymplectic matrix.
45  See also Boost, BoostX, BoostY and BoostZ for classes representing
46  specialized Lorentz transformations.
47  Also, the 3-D rotation classes can be considered to be special Lorentz
48  transformations which do not mix space and time components.
49 
50  @ingroup GenVector
51 
52  */
53 
54 class LorentzRotation {
55 
56 public:
57 
58  typedef double Scalar;
59 
60  enum ELorentzRotationMatrixIndex {
61  kXX = 0, kXY = 1, kXZ = 2, kXT = 3
62  , kYX = 4, kYY = 5, kYZ = 6, kYT = 7
63  , kZX = 8, kZY = 9, kZZ = 10, kZT = 11
64  , kTX = 12, kTY = 13, kTZ = 14, kTT = 15
65  };
66 
67  // ========== Constructors and Assignment =====================
68 
69  /**
70  Default constructor (identity transformation)
71  */
72  LorentzRotation();
73 
74  /**
75  Construct given a pair of pointers or iterators defining the
76  beginning and end of an array of sixteen Scalars
77  */
78  template<class IT>
79  LorentzRotation(IT begin, IT end) { SetComponents(begin,end); }
80 
81  // The compiler-generated and dtor are OK but we have implementwd the copy-ctor and
82  // assignment operators since we have a template assignment
83 
84  /**
85  Copy constructor
86  */
87  LorentzRotation( LorentzRotation const & r ) {
88  *this = r;
89  }
90 
91  /**
92  Construct from a pure boost
93  */
94  explicit LorentzRotation( Boost const & b ) { b.GetLorentzRotation( fM+0 ); }
95  explicit LorentzRotation( BoostX const & bx ) { bx.GetLorentzRotation( fM+0 ); }
96  explicit LorentzRotation( BoostY const & by ) { by.GetLorentzRotation( fM+0 ); }
97  explicit LorentzRotation( BoostZ const & bz ) { bz.GetLorentzRotation( fM+0 ); }
98 
99  /**
100  Construct from a 3-D rotation (no space-time mixing)
101  */
102  explicit LorentzRotation( Rotation3D const & r );
103  explicit LorentzRotation( AxisAngle const & a );
104  explicit LorentzRotation( EulerAngles const & e );
105  explicit LorentzRotation( Quaternion const & q );
106  explicit LorentzRotation( RotationX const & r );
107  explicit LorentzRotation( RotationY const & r );
108  explicit LorentzRotation( RotationZ const & r );
109 
110  /**
111  Construct from a linear algebra matrix of size at least 4x4,
112  which must support operator()(i,j) to obtain elements (0,3) thru (3,3).
113  Precondition: The matrix is assumed to be orthosymplectic. NO checking
114  or re-adjusting is performed.
115  Note: (0,0) refers to the XX component; (3,3) refers to the TT component.
116  */
117  template<class ForeignMatrix>
118  explicit LorentzRotation(const ForeignMatrix & m) { SetComponents(m); }
119 
120  /**
121  Construct from four orthosymplectic vectors (which must have methods
122  x(), y(), z() and t()) which will be used as the columns of the Lorentz
123  rotation matrix. The orthosymplectic conditions will be checked, and
124  values adjusted so that the result will always be a good Lorentz rotation
125  matrix.
126  */
127  template<class Foreign4Vector>
128  LorentzRotation(const Foreign4Vector& v1,
129  const Foreign4Vector& v2,
130  const Foreign4Vector& v3,
131  const Foreign4Vector& v4 ) { SetComponents(v1, v2, v3, v4); }
132 
133 
134  /**
135  Raw constructor from sixteen Scalar components (without any checking)
136  */
137  LorentzRotation(Scalar xx, Scalar xy, Scalar xz, Scalar xt,
138  Scalar yx, Scalar yy, Scalar yz, Scalar yt,
139  Scalar zx, Scalar zy, Scalar zz, Scalar zt,
140  Scalar tx, Scalar ty, Scalar tz, Scalar tt)
141  {
142  SetComponents (xx, xy, xz, xt,
143  yx, yy, yz, yt,
144  zx, zy, zz, zt,
145  tx, ty, tz, tt);
146  }
147 
148  /**
149  Assign from another LorentzRotation
150  */
151  LorentzRotation &
152  operator=( LorentzRotation const & rhs ) {
153  SetComponents( rhs.fM[0], rhs.fM[1], rhs.fM[2], rhs.fM[3],
154  rhs.fM[4], rhs.fM[5], rhs.fM[6], rhs.fM[7],
155  rhs.fM[8], rhs.fM[9], rhs.fM[10], rhs.fM[11],
156  rhs.fM[12], rhs.fM[13], rhs.fM[14], rhs.fM[15] );
157  return *this;
158  }
159 
160  /**
161  Assign from a pure boost
162  */
163  LorentzRotation &
164  operator=( Boost const & b ) { return operator=(LorentzRotation(b)); }
165  LorentzRotation &
166  operator=( BoostX const & b ) { return operator=(LorentzRotation(b)); }
167  LorentzRotation &
168  operator=( BoostY const & b ) { return operator=(LorentzRotation(b)); }
169  LorentzRotation &
170  operator=( BoostZ const & b ) { return operator=(LorentzRotation(b)); }
171 
172  /**
173  Assign from a 3-D rotation
174  */
175  LorentzRotation &
176  operator=( Rotation3D const & r ) { return operator=(LorentzRotation(r)); }
177  LorentzRotation &
178  operator=( AxisAngle const & a ) { return operator=(LorentzRotation(a)); }
179  LorentzRotation &
180  operator=( EulerAngles const & e ) { return operator=(LorentzRotation(e)); }
181  LorentzRotation &
182  operator=( Quaternion const & q ) { return operator=(LorentzRotation(q)); }
183  LorentzRotation &
184  operator=( RotationZ const & r ) { return operator=(LorentzRotation(r)); }
185  LorentzRotation &
186  operator=( RotationY const & r ) { return operator=(LorentzRotation(r)); }
187  LorentzRotation &
188  operator=( RotationX const & r ) { return operator=(LorentzRotation(r)); }
189 
190  /**
191  Assign from a linear algebra matrix of size at least 4x4,
192  which must support operator()(i,j) to obtain elements (0,3) thru (3,3).
193  Precondition: The matrix is assumed to be orthosymplectic. NO checking
194  or re-adjusting is performed.
195  */
196  template<class ForeignMatrix>
197  LorentzRotation &
198  operator=(const ForeignMatrix & m) {
199  SetComponents( m(0,0), m(0,1), m(0,2), m(0,3),
200  m(1,0), m(1,1), m(1,2), m(1,3),
201  m(2,0), m(2,1), m(2,2), m(2,3),
202  m(3,0), m(3,1), m(3,2), m(3,3) );
203  return *this;
204  }
205 
206  /**
207  Re-adjust components to eliminate small deviations from a perfect
208  orthosyplectic matrix.
209  */
210  void Rectify();
211 
212  // ======== Components ==============
213 
214  /**
215  Set components from four orthosymplectic vectors (which must have methods
216  x(), y(), z(), and t()) which will be used as the columns of the
217  Lorentz rotation matrix. The values will be adjusted
218  so that the result will always be a good Lorentz rotation matrix.
219  */
220  template<class Foreign4Vector>
221  void
222  SetComponents (const Foreign4Vector& v1,
223  const Foreign4Vector& v2,
224  const Foreign4Vector& v3,
225  const Foreign4Vector& v4 ) {
226  fM[kXX]=v1.x(); fM[kXY]=v2.x(); fM[kXZ]=v3.x(); fM[kXT]=v4.x();
227  fM[kYX]=v1.y(); fM[kYY]=v2.y(); fM[kYZ]=v3.y(); fM[kYT]=v4.y();
228  fM[kZX]=v1.z(); fM[kZY]=v2.z(); fM[kZZ]=v3.z(); fM[kZT]=v4.z();
229  fM[kTX]=v1.t(); fM[kTY]=v2.t(); fM[kTZ]=v3.t(); fM[kTT]=v4.t();
230  Rectify();
231  }
232 
233  /**
234  Get components into four 4-vectors which will be the (orthosymplectic)
235  columns of the rotation matrix. (The 4-vector class must have a
236  constructor from 4 Scalars used as x, y, z, t)
237  */
238  template<class Foreign4Vector>
239  void
240  GetComponents ( Foreign4Vector& v1,
241  Foreign4Vector& v2,
242  Foreign4Vector& v3,
243  Foreign4Vector& v4 ) const {
244  v1 = Foreign4Vector ( fM[kXX], fM[kYX], fM[kZX], fM[kTX] );
245  v2 = Foreign4Vector ( fM[kXY], fM[kYY], fM[kZY], fM[kTY] );
246  v3 = Foreign4Vector ( fM[kXZ], fM[kYZ], fM[kZZ], fM[kTZ] );
247  v4 = Foreign4Vector ( fM[kXT], fM[kYT], fM[kZT], fM[kTT] );
248  }
249 
250  /**
251  Set the 16 matrix components given an iterator to the start of
252  the desired data, and another to the end (16 past start).
253  */
254  template<class IT>
255 #ifndef NDEBUG
256  void SetComponents(IT begin, IT end) {
257 #else
258  void SetComponents(IT begin, IT ) {
259 #endif
260  for (int i = 0; i <16; ++i) {
261  fM[i] = *begin;
262  ++begin;
263  }
264  assert (end==begin);
265  }
266 
267  /**
268  Get the 16 matrix components into data specified by an iterator begin
269  and another to the end of the desired data (16 past start).
270  */
271  template<class IT>
272 #ifndef NDEBUG
273  void GetComponents(IT begin, IT end) const {
274 #else
275  void GetComponents(IT begin, IT ) const {
276 #endif
277  for (int i = 0; i <16; ++i) {
278  *begin = fM[i];
279  ++begin;
280  }
281  assert (end==begin);
282  }
283 
284  /**
285  Get the 16 matrix components into data specified by an iterator begin
286  */
287  template<class IT>
288  void GetComponents(IT begin) const {
289  std::copy ( fM+0, fM+16, begin );
290  }
291 
292  /**
293  Set components from a linear algebra matrix of size at least 4x4,
294  which must support operator()(i,j) to obtain elements (0,0) thru (3,3).
295  Precondition: The matrix is assumed to be orthosymplectic. NO checking
296  or re-adjusting is performed.
297  */
298  template<class ForeignMatrix>
299  void
300  SetRotationMatrix (const ForeignMatrix & m) {
301  fM[kXX]=m(0,0); fM[kXY]=m(0,1); fM[kXZ]=m(0,2); fM[kXT]=m(0,3);
302  fM[kYX]=m(1,0); fM[kYY]=m(1,1); fM[kYZ]=m(1,2); fM[kYT]=m(1,3);
303  fM[kZX]=m(2,0); fM[kZY]=m(2,1); fM[kZZ]=m(2,2); fM[kZT]=m(2,3);
304  fM[kTX]=m(3,0); fM[kTY]=m(3,1); fM[kTZ]=m(3,2); fM[kTT]=m(3,3);
305  }
306 
307  /**
308  Get components into a linear algebra matrix of size at least 4x4,
309  which must support operator()(i,j) for write access to elements
310  (0,0) thru (3,3).
311  */
312  template<class ForeignMatrix>
313  void
314  GetRotationMatrix (ForeignMatrix & m) const {
315  m(0,0)=fM[kXX]; m(0,1)=fM[kXY]; m(0,2)=fM[kXZ]; m(0,3)=fM[kXT];
316  m(1,0)=fM[kYX]; m(1,1)=fM[kYY]; m(1,2)=fM[kYZ]; m(1,3)=fM[kYT];
317  m(2,0)=fM[kZX]; m(2,1)=fM[kZY]; m(2,2)=fM[kZZ]; m(2,3)=fM[kZT];
318  m(3,0)=fM[kTX]; m(3,1)=fM[kTY]; m(3,2)=fM[kTZ]; m(3,3)=fM[kTT];
319  }
320 
321  /**
322  Set the components from sixteen scalars -- UNCHECKED for orthosymplectic
323  */
324  void
325  SetComponents (Scalar xx, Scalar xy, Scalar xz, Scalar xt,
326  Scalar yx, Scalar yy, Scalar yz, Scalar yt,
327  Scalar zx, Scalar zy, Scalar zz, Scalar zt,
328  Scalar tx, Scalar ty, Scalar tz, Scalar tt) {
329  fM[kXX]=xx; fM[kXY]=xy; fM[kXZ]=xz; fM[kXT]=xt;
330  fM[kYX]=yx; fM[kYY]=yy; fM[kYZ]=yz; fM[kYT]=yt;
331  fM[kZX]=zx; fM[kZY]=zy; fM[kZZ]=zz; fM[kZT]=zt;
332  fM[kTX]=tx; fM[kTY]=ty; fM[kTZ]=tz; fM[kTT]=tt;
333  }
334 
335  /**
336  Get the sixteen components into sixteen scalars
337  */
338  void
339  GetComponents (Scalar &xx, Scalar &xy, Scalar &xz, Scalar &xt,
340  Scalar &yx, Scalar &yy, Scalar &yz, Scalar &yt,
341  Scalar &zx, Scalar &zy, Scalar &zz, Scalar &zt,
342  Scalar &tx, Scalar &ty, Scalar &tz, Scalar &tt) const {
343  xx=fM[kXX]; xy=fM[kXY]; xz=fM[kXZ]; xt=fM[kXT];
344  yx=fM[kYX]; yy=fM[kYY]; yz=fM[kYZ]; yt=fM[kYT];
345  zx=fM[kZX]; zy=fM[kZY]; zz=fM[kZZ]; zt=fM[kZT];
346  tx=fM[kTX]; ty=fM[kTY]; tz=fM[kTZ]; tt=fM[kTT];
347  }
348 
349  // =========== operations ==============
350 
351  /**
352  Lorentz transformation operation on a Minkowski ('Cartesian')
353  LorentzVector
354  */
355  LorentzVector< ROOT::Math::PxPyPzE4D<double> >
356  operator() (const LorentzVector< ROOT::Math::PxPyPzE4D<double> > & v) const {
357  Scalar x = v.Px();
358  Scalar y = v.Py();
359  Scalar z = v.Pz();
360  Scalar t = v.E();
361  return LorentzVector< PxPyPzE4D<double> >
362  ( fM[kXX]*x + fM[kXY]*y + fM[kXZ]*z + fM[kXT]*t
363  , fM[kYX]*x + fM[kYY]*y + fM[kYZ]*z + fM[kYT]*t
364  , fM[kZX]*x + fM[kZY]*y + fM[kZZ]*z + fM[kZT]*t
365  , fM[kTX]*x + fM[kTY]*y + fM[kTZ]*z + fM[kTT]*t );
366  }
367 
368  /**
369  Lorentz transformation operation on a LorentzVector in any
370  coordinate system
371  */
372  template <class CoordSystem>
373  LorentzVector<CoordSystem>
374  operator() (const LorentzVector<CoordSystem> & v) const {
375  LorentzVector< PxPyPzE4D<double> > xyzt(v);
376  LorentzVector< PxPyPzE4D<double> > r_xyzt = operator()(xyzt);
377  return LorentzVector<CoordSystem> ( r_xyzt );
378  }
379 
380  /**
381  Lorentz transformation operation on an arbitrary 4-vector v.
382  Preconditions: v must implement methods x(), y(), z(), and t()
383  and the arbitrary vector type must have a constructor taking (x,y,z,t)
384  */
385  template <class Foreign4Vector>
386  Foreign4Vector
387  operator() (const Foreign4Vector & v) const {
388  LorentzVector< PxPyPzE4D<double> > xyzt(v);
389  LorentzVector< PxPyPzE4D<double> > r_xyzt = operator()(xyzt);
390  return Foreign4Vector ( r_xyzt.X(), r_xyzt.Y(), r_xyzt.Z(), r_xyzt.T() );
391  }
392 
393  /**
394  Overload operator * for rotation on a vector
395  */
396  template <class A4Vector>
397  inline
398  A4Vector operator* (const A4Vector & v) const
399  {
400  return operator()(v);
401  }
402 
403  /**
404  Invert a Lorentz rotation in place
405  */
406  void Invert();
407 
408  /**
409  Return inverse of a rotation
410  */
411  LorentzRotation Inverse() const;
412 
413  // ========= Multi-Rotation Operations ===============
414 
415  /**
416  Multiply (combine) this Lorentz rotation by another LorentzRotation
417  */
418  LorentzRotation operator * (const LorentzRotation & r) const;
419 
420  //#ifdef TODO_LATER
421  /**
422  Multiply (combine) this Lorentz rotation by a pure Lorentz boost
423  */
424  //TODO: implement directly in a more efficient way. Now are implemented
425  // going through another LorentzRotation
426  LorentzRotation operator * (const Boost & b) const { LorentzRotation tmp(b); return (*this)*tmp; }
427  LorentzRotation operator * (const BoostX & b) const { LorentzRotation tmp(b); return (*this)*tmp; }
428  LorentzRotation operator * (const BoostY & b) const { LorentzRotation tmp(b); return (*this)*tmp; }
429  LorentzRotation operator * (const BoostZ & b) const { LorentzRotation tmp(b); return (*this)*tmp; }
430 
431  /**
432  Multiply (combine) this Lorentz rotation by a 3-D Rotation
433  */
434  LorentzRotation operator * (const Rotation3D & r) const { LorentzRotation tmp(r); return (*this)*tmp; }
435  LorentzRotation operator * (const AxisAngle & a) const { LorentzRotation tmp(a); return (*this)*tmp; }
436  LorentzRotation operator * (const EulerAngles & e) const { LorentzRotation tmp(e); return (*this)*tmp; }
437  LorentzRotation operator * (const Quaternion & q) const { LorentzRotation tmp(q); return (*this)*tmp; }
438  LorentzRotation operator * (const RotationX & rx) const { LorentzRotation tmp(rx); return (*this)*tmp; }
439  LorentzRotation operator * (const RotationY & ry) const { LorentzRotation tmp(ry); return (*this)*tmp; }
440  LorentzRotation operator * (const RotationZ & rz) const { LorentzRotation tmp(rz); return (*this)*tmp; }
441  //#endif
442 
443  /**
444  Post-Multiply (on right) by another LorentzRotation, Boost, or
445  rotation : T = T*R
446  */
447  template <class R>
448  LorentzRotation & operator *= (const R & r) { return *this = (*this)*r; }
449 
450  /**
451  Equality/inequality operators
452  */
453  bool operator == (const LorentzRotation & rhs) const {
454  for (unsigned int i=0; i < 16; ++i) {
455  if( fM[i] != rhs.fM[i] ) return false;
456  }
457  return true;
458  }
459  bool operator != (const LorentzRotation & rhs) const {
460  return ! operator==(rhs);
461  }
462 
463 private:
464 
465  Scalar fM[16];
466 
467 }; // LorentzRotation
468 
469 // ============ Class LorentzRotation ends here ============
470 
471 
472 /**
473  Stream Output and Input
474  */
475  // TODO - I/O should be put in the manipulator form
476 
477 std::ostream & operator<< (std::ostream & os, const LorentzRotation & r);
478 
479 // ============================================ vetted to here ============
480 
481 #ifdef NOTYET
482 /**
483  Distance between two Lorentz rotations
484  */
485 template <class R>
486 inline
487 typename Rotation3D::Scalar
488 Distance ( const Rotation3D& r1, const R & r2) {return gv_detail::dist(r1,r2);}
489 #endif
490 
491 } //namespace Math
492 } //namespace ROOT
493 
494 
495 
496 
497 
498 
499 
500 #endif /* ROOT_Math_GenVector_LorentzRotation */