Math/GenVector/LorentzRotation.h

0001 // @(#)root/mathcore:$Id$
0002 // Authors: W. Brown, M. Fischler, L. Moneta    2005
0003
0004  /**********************************************************************
0005   *                                                                    *
0006   * Copyright (c) 2005 ROOT MathLib Team                               *
0007   *                                                                    *
0008   *                                                                    *
0009   **********************************************************************/
0010
0011 // Header file for LorentzRotation
0012 //
0013 // Created by: Mark Fischler  Mon Aug 8  2005
0014 //
0015 // Last update: $Id$
0016 //
0017 #ifndef ROOT_Math_GenVector_LorentzRotation
0018 #define ROOT_Math_GenVector_LorentzRotation  1
0019
0020 #include "Math/GenVector/LorentzRotationfwd.h"
0021
0022 #include "Math/GenVector/LorentzVector.h"
0023 #include "Math/GenVector/PxPyPzE4D.h"
0024
0025 #include "Math/GenVector/Rotation3Dfwd.h"
0026 #include "Math/GenVector/AxisAnglefwd.h"
0027 #include "Math/GenVector/EulerAnglesfwd.h"
0028 #include "Math/GenVector/Quaternionfwd.h"
0029 #include "Math/GenVector/RotationXfwd.h"
0030 #include "Math/GenVector/RotationYfwd.h"
0031 #include "Math/GenVector/RotationZfwd.h"
0032 #include "Math/GenVector/Boost.h"
0033 #include "Math/GenVector/BoostX.h"
0034 #include "Math/GenVector/BoostY.h"
0035 #include "Math/GenVector/BoostZ.h"
0036
0037 #include <algorithm>
0038
0039 namespace ROOT {
0040
0041   namespace Math {
0042
0043 //__________________________________________________________________________________________
0044   /**
0045      Lorentz transformation class with the (4D) transformation represented by
0046      a 4x4 orthosymplectic matrix.
0047      See also Boost, BoostX, BoostY and BoostZ for classes representing
0048      specialized Lorentz transformations.
0049      Also, the 3-D rotation classes can be considered to be special Lorentz
0050      transformations which do not mix space and time components.
0051
0052      @ingroup GenVector
0053
0054      @see GenVector
0055   */
0056
0057 class LorentzRotation {
0058
0059 public:
0060
0061   typedef double Scalar;
0062
0063   enum ELorentzRotationMatrixIndex {
0064       kXX =  0, kXY =  1, kXZ =  2, kXT =  3
0065     , kYX =  4, kYY =  5, kYZ =  6, kYT =  7
0066     , kZX =  8, kZY =  9, kZZ = 10, kZT = 11
0067     , kTX = 12, kTY = 13, kTZ = 14, kTT = 15
0068   };
0069
0070   // ========== Constructors and Assignment =====================
0071
0072   /**
0073       Default constructor (identity transformation)
0074   */
0075   LorentzRotation();
0076
0077   /**
0078      Construct given a pair of pointers or iterators defining the
0079      beginning and end of an array of sixteen Scalars
0080    */
0081   template<class IT>
0082   LorentzRotation(IT begin, IT end) { SetComponents(begin,end); }
0083
0084   // The compiler-generated and dtor are OK but we have implementwd the copy-ctor and
0085   // assignment operators since we have a template assignment
0086
0087   /**
0088      Copy constructor
0089    */
0090    LorentzRotation( LorentzRotation const & r ) {
0091       *this = r;
0092    }
0093
0094   /**
0095      Construct from a pure boost
0096   */
0097   explicit LorentzRotation( Boost  const & b  ) {  b.GetLorentzRotation( fM+0 ); }
0098   explicit LorentzRotation( BoostX const & bx ) { bx.GetLorentzRotation( fM+0 ); }
0099   explicit LorentzRotation( BoostY const & by ) { by.GetLorentzRotation( fM+0 ); }
0100   explicit LorentzRotation( BoostZ const & bz ) { bz.GetLorentzRotation( fM+0 ); }
0101
0102   /**
0103      Construct from a 3-D rotation (no space-time mixing)
0104   */
0105   explicit LorentzRotation( Rotation3D  const & r );
0106   explicit LorentzRotation( AxisAngle   const & a );
0107   explicit LorentzRotation( EulerAngles const & e );
0108   explicit LorentzRotation( Quaternion  const & q );
0109   explicit LorentzRotation( RotationX   const & r );
0110   explicit LorentzRotation( RotationY   const & r );
0111   explicit LorentzRotation( RotationZ   const & r );
0112
0113   /**
0114      Construct from a linear algebra matrix of size at least 4x4,
0115      which must support operator()(i,j) to obtain elements (0,3) thru (3,3).
0116      Precondition:  The matrix is assumed to be orthosymplectic.  NO checking
0117      or re-adjusting is performed.
0118      Note:  (0,0) refers to the XX component; (3,3) refers to the TT component.
0119   */
0120   template<class ForeignMatrix>
0121   explicit constexpr LorentzRotation(const ForeignMatrix & m) { SetComponents(m); }
0122
0123   /**
0124      Construct from four orthosymplectic vectors (which must have methods
0125      x(), y(), z() and t()) which will be used as the columns of the Lorentz
0126      rotation matrix.  The orthosymplectic conditions will be checked, and
0127      values adjusted so that the result will always be a good Lorentz rotation
0128      matrix.
0129   */
0130   template<class Foreign4Vector>
0131   LorentzRotation(const Foreign4Vector& v1,
0132                   const Foreign4Vector& v2,
0133                   const Foreign4Vector& v3,
0134                   const Foreign4Vector& v4 ) { SetComponents(v1, v2, v3, v4); }
0135
0136
0137   /**
0138      Raw constructor from sixteen Scalar components (without any checking)
0139   */
0140   LorentzRotation(Scalar  xx, Scalar  xy, Scalar  xz, Scalar xt,
0141                   Scalar  yx, Scalar  yy, Scalar  yz, Scalar yt,
0142                   Scalar  zx, Scalar  zy, Scalar  zz, Scalar zt,
0143                   Scalar  tx, Scalar  ty, Scalar  tz, Scalar tt)
0144  {
0145     SetComponents (xx, xy, xz, xt,
0146                    yx, yy, yz, yt,
0147                    zx, zy, zz, zt,
0148                    tx, ty, tz, tt);
0149  }
0150
0151   /**
0152       Assign from another LorentzRotation
0153   */
0154    LorentzRotation &
0155   operator=( LorentzRotation  const & rhs ) {
0156       SetComponents( rhs.fM[0],  rhs.fM[1],  rhs.fM[2],  rhs.fM[3],
0157                      rhs.fM[4],  rhs.fM[5],  rhs.fM[6],  rhs.fM[7],
0158                      rhs.fM[8],  rhs.fM[9],  rhs.fM[10], rhs.fM[11],
0159                      rhs.fM[12], rhs.fM[13], rhs.fM[14], rhs.fM[15] );
0160       return *this;
0161    }
0162
0163   /**
0164      Assign from a pure boost
0165   */
0166   LorentzRotation &
0167   operator=( Boost  const & b ) { return operator=(LorentzRotation(b)); }
0168   LorentzRotation &
0169   operator=( BoostX const & b ) { return operator=(LorentzRotation(b)); }
0170   LorentzRotation &
0171   operator=( BoostY const & b ) { return operator=(LorentzRotation(b)); }
0172   LorentzRotation &
0173   operator=( BoostZ const & b ) { return operator=(LorentzRotation(b)); }
0174
0175   /**
0176      Assign from a 3-D rotation
0177   */
0178   LorentzRotation &
0179   operator=( Rotation3D  const & r ) { return operator=(LorentzRotation(r)); }
0180   LorentzRotation &
0181   operator=( AxisAngle   const & a ) { return operator=(LorentzRotation(a)); }
0182   LorentzRotation &
0183   operator=( EulerAngles const & e ) { return operator=(LorentzRotation(e)); }
0184   LorentzRotation &
0185   operator=( Quaternion  const & q ) { return operator=(LorentzRotation(q)); }
0186   LorentzRotation &
0187   operator=( RotationZ   const & r ) { return operator=(LorentzRotation(r)); }
0188   LorentzRotation &
0189   operator=( RotationY   const & r ) { return operator=(LorentzRotation(r)); }
0190   LorentzRotation &
0191   operator=( RotationX   const & r ) { return operator=(LorentzRotation(r)); }
0192
0193   /**
0194      Assign from a linear algebra matrix of size at least 4x4,
0195      which must support operator()(i,j) to obtain elements (0,3) thru (3,3).
0196      Precondition:  The matrix is assumed to be orthosymplectic.  NO checking
0197      or re-adjusting is performed.
0198   */
0199   template<class ForeignMatrix>
0200   LorentzRotation &
0201   operator=(const ForeignMatrix & m) {
0202      SetComponents( m(0,0), m(0,1), m(0,2), m(0,3),
0203                     m(1,0), m(1,1), m(1,2), m(1,3),
0204                     m(2,0), m(2,1), m(2,2), m(2,3),
0205                     m(3,0), m(3,1), m(3,2), m(3,3) );
0206      return *this;
0207   }
0208
0209   /**
0210      Re-adjust components to eliminate small deviations from a perfect
0211      orthosyplectic matrix.
0212    */
0213   void Rectify();
0214
0215   // ======== Components ==============
0216
0217   /**
0218      Set components from four orthosymplectic vectors (which must have methods
0219      x(), y(), z(), and t()) which will be used as the columns of the
0220      Lorentz rotation matrix.  The values will be adjusted
0221      so that the result will always be a good Lorentz rotation matrix.
0222   */
0223   template<class Foreign4Vector>
0224   void
0225   SetComponents (const Foreign4Vector& v1,
0226                  const Foreign4Vector& v2,
0227                  const Foreign4Vector& v3,
0228                  const Foreign4Vector& v4 ) {
0229     fM[kXX]=v1.x();  fM[kXY]=v2.x();  fM[kXZ]=v3.x();  fM[kXT]=v4.x();
0230     fM[kYX]=v1.y();  fM[kYY]=v2.y();  fM[kYZ]=v3.y();  fM[kYT]=v4.y();
0231     fM[kZX]=v1.z();  fM[kZY]=v2.z();  fM[kZZ]=v3.z();  fM[kZT]=v4.z();
0232     fM[kTX]=v1.t();  fM[kTY]=v2.t();  fM[kTZ]=v3.t();  fM[kTT]=v4.t();
0233     Rectify();
0234   }
0235
0236   /**
0237      Get components into four 4-vectors which will be the (orthosymplectic)
0238      columns of the rotation matrix.  (The 4-vector class must have a
0239      constructor from 4 Scalars used as x, y, z, t)
0240   */
0241   template<class Foreign4Vector>
0242   void
0243   GetComponents ( Foreign4Vector& v1,
0244                   Foreign4Vector& v2,
0245                   Foreign4Vector& v3,
0246                   Foreign4Vector& v4 ) const {
0247     v1 = Foreign4Vector ( fM[kXX], fM[kYX], fM[kZX], fM[kTX] );
0248     v2 = Foreign4Vector ( fM[kXY], fM[kYY], fM[kZY], fM[kTY] );
0249     v3 = Foreign4Vector ( fM[kXZ], fM[kYZ], fM[kZZ], fM[kTZ] );
0250     v4 = Foreign4Vector ( fM[kXT], fM[kYT], fM[kZT], fM[kTT] );
0251   }
0252
0253   /**
0254      Set the 16 matrix components given an iterator to the start of
0255      the desired data, and another to the end (16 past start).
0256    */
0257   template<class IT>
0258   void SetComponents(IT begin, IT end) {
0259      for (int i = 0; i <16; ++i) {
0260         fM[i] = *begin;
0261         ++begin;
0262      }
0263      (void)end;
0264      assert (end==begin);
0265   }
0266
0267   /**
0268      Get the 16 matrix components into data specified by an iterator begin
0269      and another to the end of the desired data (16 past start).
0270    */
0271   template<class IT>
0272   void GetComponents(IT begin, IT end) const {
0273      for (int i = 0; i <16; ++i) {
0274         *begin = fM[i];
0275         ++begin;
0276      }
0277      (void)end;
0278      assert (end==begin);
0279   }
0280
0281   /**
0282      Get the 16 matrix components into data specified by an iterator begin
0283    */
0284   template<class IT>
0285   void GetComponents(IT begin) const {
0286     std::copy ( fM+0, fM+16, begin );
0287   }
0288
0289   /**
0290      Set components from a linear algebra matrix of size at least 4x4,
0291      which must support operator()(i,j) to obtain elements (0,0) thru (3,3).
0292      Precondition:  The matrix is assumed to be orthosymplectic.  NO checking
0293      or re-adjusting is performed.
0294   */
0295   template<class ForeignMatrix>
0296   void
0297   SetRotationMatrix (const ForeignMatrix & m) {
0298     fM[kXX]=m(0,0);  fM[kXY]=m(0,1);  fM[kXZ]=m(0,2);  fM[kXT]=m(0,3);
0299     fM[kYX]=m(1,0);  fM[kYY]=m(1,1);  fM[kYZ]=m(1,2);  fM[kYT]=m(1,3);
0300     fM[kZX]=m(2,0);  fM[kZY]=m(2,1);  fM[kZZ]=m(2,2);  fM[kZT]=m(2,3);
0301     fM[kTX]=m(3,0);  fM[kTY]=m(3,1);  fM[kTZ]=m(3,2);  fM[kTT]=m(3,3);
0302   }
0303
0304   /**
0305      Get components into a linear algebra matrix of size at least 4x4,
0306      which must support operator()(i,j) for write access to elements
0307      (0,0) thru (3,3).
0308   */
0309   template<class ForeignMatrix>
0310   void
0311   GetRotationMatrix (ForeignMatrix & m) const {
0312     m(0,0)=fM[kXX];  m(0,1)=fM[kXY];  m(0,2)=fM[kXZ]; m(0,3)=fM[kXT];
0313     m(1,0)=fM[kYX];  m(1,1)=fM[kYY];  m(1,2)=fM[kYZ]; m(1,3)=fM[kYT];
0314     m(2,0)=fM[kZX];  m(2,1)=fM[kZY];  m(2,2)=fM[kZZ]; m(2,3)=fM[kZT];
0315     m(3,0)=fM[kTX];  m(3,1)=fM[kTY];  m(3,2)=fM[kTZ]; m(3,3)=fM[kTT];
0316   }
0317
0318   /**
0319      Set the components from sixteen scalars -- UNCHECKED for orthosymplectic
0320    */
0321   void
0322   SetComponents (Scalar  xx, Scalar  xy, Scalar  xz, Scalar  xt,
0323                  Scalar  yx, Scalar  yy, Scalar  yz, Scalar  yt,
0324                  Scalar  zx, Scalar  zy, Scalar  zz, Scalar  zt,
0325                  Scalar  tx, Scalar  ty, Scalar  tz, Scalar  tt) {
0326                  fM[kXX]=xx;  fM[kXY]=xy;  fM[kXZ]=xz;  fM[kXT]=xt;
0327                  fM[kYX]=yx;  fM[kYY]=yy;  fM[kYZ]=yz;  fM[kYT]=yt;
0328                  fM[kZX]=zx;  fM[kZY]=zy;  fM[kZZ]=zz;  fM[kZT]=zt;
0329                  fM[kTX]=tx;  fM[kTY]=ty;  fM[kTZ]=tz;  fM[kTT]=tt;
0330   }
0331
0332   /**
0333      Get the sixteen components into sixteen scalars
0334    */
0335   void
0336   GetComponents (Scalar &xx, Scalar &xy, Scalar &xz, Scalar &xt,
0337                  Scalar &yx, Scalar &yy, Scalar &yz, Scalar &yt,
0338                  Scalar &zx, Scalar &zy, Scalar &zz, Scalar &zt,
0339                  Scalar &tx, Scalar &ty, Scalar &tz, Scalar &tt) const {
0340                  xx=fM[kXX];  xy=fM[kXY];  xz=fM[kXZ];  xt=fM[kXT];
0341                  yx=fM[kYX];  yy=fM[kYY];  yz=fM[kYZ];  yt=fM[kYT];
0342                  zx=fM[kZX];  zy=fM[kZY];  zz=fM[kZZ];  zt=fM[kZT];
0343                  tx=fM[kTX];  ty=fM[kTY];  tz=fM[kTZ];  tt=fM[kTT];
0344   }
0345
0346   // =========== operations ==============
0347
0348   /**
0349      Lorentz transformation operation on a Minkowski ('Cartesian')
0350      LorentzVector
0351   */
0352   LorentzVector< ROOT::Math::PxPyPzE4D<double> >
0353   operator() (const LorentzVector< ROOT::Math::PxPyPzE4D<double> > & v) const {
0354         Scalar x = v.Px();
0355         Scalar y = v.Py();
0356         Scalar z = v.Pz();
0357         Scalar t = v.E();
0358         return LorentzVector< PxPyPzE4D<double> >
0359            ( fM[kXX]*x + fM[kXY]*y + fM[kXZ]*z + fM[kXT]*t
0360              , fM[kYX]*x + fM[kYY]*y + fM[kYZ]*z + fM[kYT]*t
0361              , fM[kZX]*x + fM[kZY]*y + fM[kZZ]*z + fM[kZT]*t
0362              , fM[kTX]*x + fM[kTY]*y + fM[kTZ]*z + fM[kTT]*t );
0363   }
0364
0365   /**
0366      Lorentz transformation operation on a LorentzVector in any
0367      coordinate system
0368    */
0369   template <class CoordSystem>
0370   LorentzVector<CoordSystem>
0371   operator() (const LorentzVector<CoordSystem> & v) const {
0372     LorentzVector< PxPyPzE4D<double> > xyzt(v);
0373     LorentzVector< PxPyPzE4D<double> > r_xyzt = operator()(xyzt);
0374     return LorentzVector<CoordSystem> ( r_xyzt );
0375   }
0376
0377   /**
0378      Lorentz transformation operation on an arbitrary 4-vector v.
0379      Preconditions:  v must implement methods x(), y(), z(), and t()
0380      and the arbitrary vector type must have a constructor taking (x,y,z,t)
0381    */
0382   template <class Foreign4Vector>
0383   Foreign4Vector
0384   operator() (const Foreign4Vector & v) const {
0385     LorentzVector< PxPyPzE4D<double> > xyzt(v);
0386     LorentzVector< PxPyPzE4D<double> > r_xyzt = operator()(xyzt);
0387     return Foreign4Vector ( r_xyzt.X(), r_xyzt.Y(), r_xyzt.Z(), r_xyzt.T() );
0388   }
0389
0390   /**
0391      Overload operator * for rotation on a vector
0392    */
0393   template <class A4Vector>
0394   inline
0395   A4Vector operator* (const A4Vector & v) const
0396   {
0397     return operator()(v);
0398   }
0399
0400   /**
0401       Invert a Lorentz rotation in place
0402    */
0403   void Invert();
0404
0405   /**
0406       Return inverse of  a rotation
0407    */
0408   LorentzRotation Inverse() const;
0409
0410   // ========= Multi-Rotation Operations ===============
0411
0412   /**
0413      Multiply (combine) this Lorentz rotation by another LorentzRotation
0414    */
0415   LorentzRotation operator * (const LorentzRotation & r) const;
0416
0417   //#ifdef TODO_LATER
0418   /**
0419      Multiply (combine) this Lorentz rotation by a pure Lorentz boost
0420    */
0421   //TODO: implement directly in a more efficient way. Now are implemented
0422   // going through another LorentzRotation
0423   LorentzRotation operator * (const Boost  & b) const  { LorentzRotation tmp(b); return (*this)*tmp; }
0424   LorentzRotation operator * (const BoostX & b) const  { LorentzRotation tmp(b); return (*this)*tmp; }
0425   LorentzRotation operator * (const BoostY & b) const  { LorentzRotation tmp(b); return (*this)*tmp; }
0426   LorentzRotation operator * (const BoostZ & b) const  { LorentzRotation tmp(b); return (*this)*tmp; }
0427
0428   /**
0429      Multiply (combine) this Lorentz rotation by a 3-D Rotation
0430    */
0431   LorentzRotation operator * (const Rotation3D  & r) const { LorentzRotation tmp(r); return (*this)*tmp; }
0432   LorentzRotation operator * (const AxisAngle   & a) const { LorentzRotation tmp(a); return (*this)*tmp; }
0433   LorentzRotation operator * (const EulerAngles & e) const { LorentzRotation tmp(e); return (*this)*tmp; }
0434   LorentzRotation operator * (const Quaternion  & q) const { LorentzRotation tmp(q); return (*this)*tmp; }
0435   LorentzRotation operator * (const RotationX  & rx) const { LorentzRotation tmp(rx); return (*this)*tmp; }
0436   LorentzRotation operator * (const RotationY  & ry) const { LorentzRotation tmp(ry); return (*this)*tmp; }
0437   LorentzRotation operator * (const RotationZ  & rz) const { LorentzRotation tmp(rz); return (*this)*tmp; }
0438   //#endif
0439
0440   /**
0441      Post-Multiply (on right) by another LorentzRotation, Boost, or
0442      rotation :  T = T*R
0443    */
0444   template <class R>
0445   LorentzRotation & operator *= (const R & r) { return *this = (*this)*r; }
0446
0447   /**
0448      Equality/inequality operators
0449    */
0450   bool operator == (const LorentzRotation & rhs) const {
0451     for (unsigned int i=0; i < 16; ++i) {
0452       if( fM[i] != rhs.fM[i] )  return false;
0453     }
0454     return true;
0455   }
0456   bool operator != (const LorentzRotation & rhs) const {
0457     return ! operator==(rhs);
0458   }
0459
0460 private:
0461
0462   Scalar fM[16];
0463
0464 };  // LorentzRotation
0465
0466 // ============ Class LorentzRotation ends here ============
0467
0468
0469 /**
0470    Stream Output and Input
0471  */
0472   // TODO - I/O should be put in the manipulator form
0473
0474 std::ostream & operator<< (std::ostream & os, const LorentzRotation & r);
0475
0476 // ============================================ vetted to here  ============
0477
0478 #ifdef NOTYET
0479 /**
0480    Distance between two Lorentz rotations
0481  */
0482 template <class R>
0483 inline
0484 typename Rotation3D::Scalar
0485 Distance ( const Rotation3D& r1, const R & r2) {return gv_detail::dist(r1,r2);}
0486 #endif
0487
0488 } //namespace Math
0489 } //namespace ROOT
0490
0491
0492
0493
0494
0495
0496
0497 #endif /* ROOT_Math_GenVector_LorentzRotation  */