Math/GenVector/RotationZYX.h

0001 // @(#)root/mathcore:$Id$
0002 // Authors: J. Palacios, L. Moneta    2007
0003
0004  /**********************************************************************
0005   *                                                                    *
0006   * Copyright (c) 2007 , LCG ROOT MathLib Team                         *
0007   *                                                                    *
0008   *                                                                    *
0009   **********************************************************************/
0010
0011 // Header file for class Rotation in 3 dimensions, described by 3 Z-Y-X  Euler angles
0012 // representing a rotation along Z, Y and X
0013 //
0014 // Created by: Lorenzo Moneta, Wed. May 22, 2007
0015 //
0016 // Last update: $Id$
0017 //
0018 #ifndef ROOT_Math_GenVector_RotationZYX
0019 #define ROOT_Math_GenVector_RotationZYX  1
0020
0021 #include "Math/Math.h"
0022
0023 #include "Math/GenVector/Rotation3D.h"
0024
0025
0026 #include "Math/GenVector/DisplacementVector3D.h"
0027
0028 #include "Math/GenVector/PositionVector3D.h"
0029
0030 #include "Math/GenVector/LorentzVector.h"
0031
0032 #include "Math/GenVector/3DConversions.h"
0033
0034 #include "TMath.h"
0035
0036 #include <algorithm>
0037 #include <cassert>
0038 #include <iostream>
0039
0040
0041 namespace ROOT {
0042 namespace Math {
0043
0044
0045 //__________________________________________________________________________________________
0046   /**
0047      Rotation class with the (3D) rotation represented by
0048      angles describing first a rotation of
0049      an angle phi (yaw) about the  Z axis,
0050      followed by a rotation of an angle theta (pitch) about the Y axis,
0051      followed by a third rotation of an angle psi (roll) about the X axis.
0052      Note that the rotations are extrinsic rotations happening around a fixed coordinate system.
0053      This is  different than the convention of the ROOT::Math::EulerAngles class, where the rotation are intrinsic.
0054      Also it has not to be confused with the typical Goldstein definition of the Euler Angles
0055      (Z-X-Z or 313 sequence) which is used by the ROOT::Math::EulerAngles class, while the sequence here is Z-Y-X or 321.
0056      Applying a RotationZYX(phi, theta, psi)  to a vector is then equal to applying RotationX(psi) * RotationY(theta) * RotationZ(phi) to the same vector.
0057
0058
0059      @ingroup GenVector
0060
0061      @see GenVector
0062   */
0063
0064 class RotationZYX {
0065
0066 public:
0067
0068    typedef double Scalar;
0069
0070
0071    // ========== Constructors and Assignment =====================
0072
0073    /**
0074       Default constructor
0075    */
0076    RotationZYX() : fPhi(0.0), fTheta(0.0), fPsi(0.0) { }
0077
0078    /**
0079       Constructor from phi, theta and psi
0080    */
0081    RotationZYX( Scalar phi, Scalar theta, Scalar psi ) :
0082       fPhi(phi), fTheta(theta), fPsi(psi)
0083    {Rectify();}  // Added 27 Jan. 06   JMM
0084
0085    /**
0086       Construct given a pair of pointers or iterators defining the
0087       beginning and end of an array of three Scalars, to be treated as
0088       the angles phi, theta and psi.
0089    */
0090    template<class IT>
0091    RotationZYX(IT begin, IT end) { SetComponents(begin,end); }
0092
0093    // The compiler-generated copy ctor, copy assignment, and dtor are OK.
0094
0095    /**
0096       Re-adjust components place angles in canonical ranges
0097    */
0098    void Rectify();
0099
0100
0101    // ======== Construction and Assignment From other Rotation Forms ==================
0102
0103    /**
0104       Construct from another supported rotation type (see gv_detail::convert )
0105    */
0106    template <class OtherRotation>
0107    explicit constexpr RotationZYX(const OtherRotation & r) {gv_detail::convert(r,*this);}
0108
0109
0110    /**
0111       Assign from another supported rotation type (see gv_detail::convert )
0112    */
0113    template <class OtherRotation>
0114    RotationZYX & operator=( OtherRotation const  & r ) {
0115       gv_detail::convert(r,*this);
0116       return *this;
0117    }
0118
0119
0120    // ======== Components ==============
0121
0122    /**
0123       Set the three Euler angles given a pair of pointers or iterators
0124       defining the beginning and end of an array of three Scalars.
0125    */
0126    template<class IT>
0127    void SetComponents(IT begin, IT end) {
0128       fPhi   = *begin++;
0129       fTheta = *begin++;
0130       fPsi   = *begin++;
0131       (void)end;
0132       assert(begin == end);
0133       Rectify();
0134    }
0135
0136    /**
0137       Get the axis and then the angle into data specified by an iterator begin
0138       and another to the end of the desired data (4 past start).
0139    */
0140    template<class IT>
0141    void GetComponents(IT begin, IT end) const {
0142       *begin++ = fPhi;
0143       *begin++ = fTheta;
0144       *begin++ = fPsi;
0145       (void)end;
0146       assert(begin == end);
0147    }
0148
0149    /**
0150       Get the axis and then the angle into data specified by an iterator begin
0151    */
0152    template<class IT>
0153    void GetComponents(IT begin) const {
0154       *begin++ = fPhi;
0155       *begin++ = fTheta;
0156       *begin   = fPsi;
0157    }
0158
0159    /**
0160       Set the components phi, theta, psi based on three Scalars.
0161    */
0162    void SetComponents(Scalar phi, Scalar theta, Scalar psi) {
0163       fPhi=phi; fTheta=theta; fPsi=psi;
0164       Rectify();
0165    }
0166
0167    /**
0168       Get the components phi, theta, psi into three Scalars.
0169    */
0170    void GetComponents(Scalar & phi, Scalar & theta, Scalar & psi) const {
0171       phi=fPhi; theta=fTheta; psi=fPsi;
0172    }
0173
0174    /**
0175       Set Phi angle (Z rotation angle)
0176    */
0177    void SetPhi(Scalar phi) { fPhi=phi; Rectify(); }
0178
0179    /**
0180       Return Phi angle (Z rotation angle)
0181    */
0182    Scalar Phi() const { return fPhi; }
0183
0184    /**
0185       Set Theta angle (Y' rotation angle)
0186    */
0187    void SetTheta(Scalar theta) { fTheta=theta; Rectify(); }
0188
0189    /**
0190       Return Theta angle (Y' rotation angle)
0191    */
0192    Scalar Theta() const { return fTheta; }
0193
0194    /**
0195       Set Psi angle (X'' rotation angle)
0196    */
0197    void SetPsi(Scalar psi) { fPsi=psi; Rectify(); }
0198
0199    /**
0200       Return Psi angle (X'' rotation angle)
0201    */
0202    Scalar Psi() const { return fPsi; }
0203
0204    // =========== operations ==============
0205
0206
0207    /**
0208       Rotation operation on a displacement vector in any coordinate system and tag
0209    */
0210    template <class CoordSystem, class U>
0211    DisplacementVector3D<CoordSystem,U>
0212    operator() (const DisplacementVector3D<CoordSystem,U> & v) const {
0213       return Rotation3D(*this) ( v );
0214    }
0215
0216    /**
0217       Rotation operation on a position vector in any coordinate system
0218    */
0219    template <class CoordSystem, class U>
0220    PositionVector3D<CoordSystem, U>
0221    operator() (const PositionVector3D<CoordSystem,U> & v) const {
0222       DisplacementVector3D< Cartesian3D<double>,U > xyz(v);
0223       DisplacementVector3D< Cartesian3D<double>,U > rxyz = operator()(xyz);
0224       return PositionVector3D<CoordSystem,U> ( rxyz );
0225    }
0226
0227    /**
0228       Rotation operation on a Lorentz vector in any 4D coordinate system
0229    */
0230    template <class CoordSystem>
0231    LorentzVector<CoordSystem>
0232    operator() (const LorentzVector<CoordSystem> & v) const {
0233       DisplacementVector3D< Cartesian3D<double> > xyz(v.Vect());
0234       xyz = operator()(xyz);
0235       LorentzVector< PxPyPzE4D<double> > xyzt (xyz.X(), xyz.Y(), xyz.Z(), v.E());
0236       return LorentzVector<CoordSystem> ( xyzt );
0237    }
0238
0239    /**
0240       Rotation operation on an arbitrary vector v.
0241       Preconditions:  v must implement methods x(), y(), and z()
0242       and the arbitrary vector type must have a constructor taking (x,y,z)
0243    */
0244    template <class ForeignVector>
0245    ForeignVector
0246    operator() (const  ForeignVector & v) const {
0247       DisplacementVector3D< Cartesian3D<double> > xyz(v);
0248       DisplacementVector3D< Cartesian3D<double> > rxyz = operator()(xyz);
0249       return ForeignVector ( rxyz.X(), rxyz.Y(), rxyz.Z() );
0250    }
0251
0252    /**
0253       Overload operator * for rotation on a vector
0254    */
0255    template <class AVector>
0256    inline
0257    AVector operator* (const AVector & v) const
0258    {
0259       return operator()(v);
0260    }
0261
0262    /**
0263       Invert a rotation in place
0264    */
0265    void Invert();
0266
0267    /**
0268       Return inverse of a rotation
0269    */
0270    RotationZYX Inverse() const {
0271       RotationZYX r(*this);
0272       r.Invert();
0273       return r;
0274    }
0275
0276
0277    // ========= Multi-Rotation Operations ===============
0278
0279    /**
0280       Multiply (combine) two rotations
0281    */
0282    RotationZYX operator * (const RotationZYX & e) const;
0283    RotationZYX operator * (const Rotation3D  & r) const;
0284    RotationZYX operator * (const AxisAngle   & a) const;
0285    RotationZYX operator * (const Quaternion  & q) const;
0286    RotationZYX operator * (const EulerAngles & q) const;
0287    RotationZYX operator * (const RotationX  & rx) const;
0288    RotationZYX operator * (const RotationY  & ry) const;
0289    RotationZYX operator * (const RotationZ  & rz) const;
0290
0291    /**
0292       Post-Multiply (on right) by another rotation :  T = T*R
0293    */
0294    template <class R>
0295    RotationZYX & operator *= (const R & r) { return *this = (*this)*r; }
0296
0297    /**
0298       Distance between two rotations
0299    */
0300    template <class R>
0301    Scalar Distance ( const R & r ) const {return gv_detail::dist(*this,r);}
0302
0303    /**
0304       Equality/inequality operators
0305    */
0306    bool operator == (const RotationZYX & rhs) const {
0307       if( fPhi   != rhs.fPhi   ) return false;
0308       if( fTheta != rhs.fTheta ) return false;
0309       if( fPsi   != rhs.fPsi   ) return false;
0310       return true;
0311    }
0312    bool operator != (const RotationZYX & rhs) const {
0313       return ! operator==(rhs);
0314    }
0315
0316 private:
0317
0318    double fPhi;      // Z rotation angle (yaw)    defined in (-PI,PI]
0319    double fTheta;    // Y' rotation angle (pitch) defined in [-PI/2,PI/2]
0320    double fPsi;      // X'' rotation angle (roll) defined in (-PI,PI]
0321
0322    static double Pi() { return TMath::Pi(); }
0323
0324 };  // RotationZYX
0325
0326 /**
0327    Distance between two rotations
0328  */
0329 template <class R>
0330 inline
0331 typename RotationZYX::Scalar
0332 Distance ( const RotationZYX& r1, const R & r2) {return gv_detail::dist(r1,r2);}
0333
0334 /**
0335    Multiplication of an axial rotation by an AxisAngle
0336  */
0337 RotationZYX operator* (RotationX const & r1, RotationZYX const & r2);
0338 RotationZYX operator* (RotationY const & r1, RotationZYX const & r2);
0339 RotationZYX operator* (RotationZ const & r1, RotationZYX const & r2);
0340
0341 /**
0342    Stream Output and Input
0343  */
0344   // TODO - I/O should be put in the manipulator form
0345
0346 std::ostream & operator<< (std::ostream & os, const RotationZYX & e);
0347
0348
0349 } // namespace Math
0350 } // namespace ROOT
0351
0352 #endif // ROOT_Math_GenVector_RotationZYX