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