Back to home page

EIC code displayed by LXR

 
 

    


File indexing completed on 2026-09-16 09:25:15

0001 // @(#)root/mathcore:$Id$
0002 // Authors: W. Brown, M. Fischler, L. Moneta    2005
0003 
0004 /**********************************************************************
0005  *                                                                    *
0006  * Copyright (c) 2005 , LCG ROOT MathLib Team                         *
0007  *                                                                    *
0008  *                                                                    *
0009  **********************************************************************/
0010 
0011 // Header file for class Transform3D
0012 //
0013 // Created by: Lorenzo Moneta  October 21 2005
0014 //
0015 //
0016 #ifndef ROOT_Math_GenVector_Transform3D
0017 #define ROOT_Math_GenVector_Transform3D  1
0018 
0019 
0020 
0021 #include "Math/GenVector/DisplacementVector3D.h"
0022 
0023 #include "Math/GenVector/PositionVector3D.h"
0024 
0025 #include "Math/GenVector/Rotation3D.h"
0026 
0027 #include "Math/GenVector/Translation3D.h"
0028 
0029 
0030 #include "Math/GenVector/AxisAnglefwd.h"
0031 #include "Math/GenVector/EulerAnglesfwd.h"
0032 #include "Math/GenVector/Quaternionfwd.h"
0033 #include "Math/GenVector/RotationZYXfwd.h"
0034 #include "Math/GenVector/RotationXfwd.h"
0035 #include "Math/GenVector/RotationYfwd.h"
0036 #include "Math/GenVector/RotationZfwd.h"
0037 
0038 #include <iostream>
0039 #include <type_traits>
0040 #include <cmath>
0041 
0042 //#include "Math/Vector3Dfwd.h"
0043 
0044 
0045 
0046 namespace ROOT {
0047 
0048 namespace Math {
0049 
0050 namespace Impl {
0051 
0052 //_________________________________________________________________________________________
0053 /**
0054     Basic 3D Transformation class describing  a rotation and then a translation
0055     The internal data are a 3D rotation data (represented as a 3x3 matrix) and a 3D vector data.
0056     They are represented and held in this class like a 3x4 matrix (a simple array of 12 numbers).
0057 
0058     The class can be constructed from any 3D rotation object
0059     (ROOT::Math::Rotation3D, ROOT::Math::AxisAngle, ROOT::Math::Quaternion, etc...) and/or
0060     a 3D Vector (ROOT::Math::DislacementVector3D or via ROOT::Math::Translation ) representing a Translation.
0061     The Transformation is defined by applying first the rotation and then the translation.
0062     A transformation defined by applying first a translation and then a rotation is equivalent to the
0063     transformation obtained applying first the rotation and then a translation equivalent to the rotated vector.
0064     The operator * can be used to obtain directly such transformations, in addition to combine various
0065     transformations.
0066     Keep in mind that the operator * (like in the case of rotations ) is not commutative.
0067     The operator * is used (in addition to operator() ) to apply a transformations on the vector
0068     (DisplacementVector3D and LorentzVector classes) and point (PositionVector3D)  classes.
0069     In the case of Vector objects the transformation only rotates them and does not translate them.
0070     Only Point objects are able to be both rotated and translated.
0071 
0072 
0073     @ingroup GenVector
0074 
0075     @see GenVector
0076 
0077 */
0078 
0079 template <typename T = double>
0080 class Transform3D {
0081 
0082 public:
0083    typedef T Scalar;
0084 
0085    typedef DisplacementVector3D<Cartesian3D<T>, DefaultCoordinateSystemTag> Vector;
0086    typedef PositionVector3D<Cartesian3D<T>, DefaultCoordinateSystemTag>     Point;
0087 
0088    enum ETransform3DMatrixIndex {
0089       kXX = 0, kXY = 1, kXZ = 2, kDX = 3,
0090       kYX = 4, kYY = 5, kYZ = 6, kDY = 7,
0091       kZX = 8, kZY = 9, kZZ =10, kDZ = 11
0092    };
0093 
0094 
0095 
0096    /**
0097        Default constructor (identy rotation) + zero translation
0098    */
0099    Transform3D()
0100    {
0101       SetIdentity();
0102    }
0103 
0104    /**
0105       Construct given a pair of pointers or iterators defining the
0106       beginning and end of an array of 12 Scalars
0107    */
0108    template<class IT>
0109    Transform3D(IT begin, IT end)
0110    {
0111       SetComponents(begin,end);
0112    }
0113 
0114    /**
0115       Construct from a rotation and then a translation described by a Vector
0116    */
0117    Transform3D( const Rotation3D & r, const Vector & v)
0118    {
0119       AssignFrom( r, v );
0120    }
0121    /**
0122       Construct from a rotation and then a translation described by a Translation3D class
0123    */
0124    Transform3D(const Rotation3D &r, const Translation3D<T> &t) { AssignFrom(r, t.Vect()); }
0125 
0126    /**
0127       Construct from a rotation (any rotation object)  and then a translation
0128       (represented by any DisplacementVector)
0129       The requirements on the rotation and vector objects are that they can be transformed in a
0130       Rotation3D class and in a Cartesian3D Vector
0131    */
0132    template <class ARotation, class CoordSystem, class Tag>
0133    Transform3D( const ARotation & r, const DisplacementVector3D<CoordSystem,Tag> & v)
0134    {
0135       AssignFrom( Rotation3D(r), Vector (v.X(),v.Y(),v.Z()) );
0136    }
0137 
0138    /**
0139       Construct from a rotation (any rotation object)  and then a translation
0140       represented by a Translation3D class
0141       The requirements on the rotation is that it can be transformed in a
0142       Rotation3D class
0143    */
0144    template <class ARotation>
0145    Transform3D(const ARotation &r, const Translation3D<T> &t)
0146    {
0147       AssignFrom( Rotation3D(r), t.Vect() );
0148    }
0149 
0150 
0151 #ifdef OLD_VERSION
0152    /**
0153       Construct from a translation and then a rotation (inverse assignment)
0154    */
0155    Transform3D( const Vector & v, const Rotation3D & r)
0156    {
0157       // is equivalent from having first the rotation and then the translation vector rotated
0158       AssignFrom( r, r(v) );
0159    }
0160 #endif
0161 
0162    /**
0163       Construct from a 3D Rotation only with zero translation
0164    */
0165    explicit constexpr Transform3D( const Rotation3D & r) {
0166       AssignFrom(r);
0167    }
0168 
0169    // convenience methods for constructing a Transform3D from all the 3D rotations classes
0170    // (cannot use templates for conflict with LA)
0171 
0172    explicit constexpr Transform3D( const AxisAngle & r) {
0173       AssignFrom(Rotation3D(r));
0174    }
0175    explicit constexpr Transform3D( const EulerAngles & r) {
0176       AssignFrom(Rotation3D(r));
0177    }
0178    explicit constexpr Transform3D( const Quaternion & r) {
0179       AssignFrom(Rotation3D(r));
0180    }
0181    explicit constexpr Transform3D( const RotationZYX & r) {
0182       AssignFrom(Rotation3D(r));
0183    }
0184 
0185    // Constructors from axial rotations
0186    // TO DO: implement direct methods for axial rotations without going through Rotation3D
0187    explicit constexpr Transform3D( const RotationX & r) {
0188       AssignFrom(Rotation3D(r));
0189    }
0190    explicit constexpr Transform3D( const RotationY & r) {
0191       AssignFrom(Rotation3D(r));
0192    }
0193    explicit constexpr Transform3D( const RotationZ & r) {
0194       AssignFrom(Rotation3D(r));
0195    }
0196 
0197    /**
0198       Construct from a translation only, represented by any DisplacementVector3D
0199       and with an identity rotation
0200    */
0201    template<class CoordSystem, class Tag>
0202    explicit constexpr Transform3D( const DisplacementVector3D<CoordSystem,Tag> & v) {
0203       AssignFrom(Vector(v.X(),v.Y(),v.Z()));
0204    }
0205    /**
0206       Construct from a translation only, represented by a Cartesian 3D Vector,
0207       and with an identity rotation
0208    */
0209    explicit constexpr Transform3D( const Vector & v) {
0210       AssignFrom(v);
0211    }
0212    /**
0213       Construct from a translation only, represented by a Translation3D class
0214       and with an identity rotation
0215    */
0216    explicit constexpr Transform3D(const Translation3D<T> &t) { AssignFrom(t.Vect()); }
0217 
0218 #ifdef OLD_VERSION
0219    /**
0220       Construct from a translation (using any type of DisplacementVector )
0221       and then a rotation (any rotation object).
0222       Requirement on the rotation and vector objects are that they can be transformed in a
0223       Rotation3D class and in a Vector
0224    */
0225    template <class ARotation, class CoordSystem, class Tag>
0226    Transform3D(const DisplacementVector3D<CoordSystem,Tag> & v , const ARotation & r)
0227    {
0228       // is equivalent from having first the rotation and then the translation vector rotated
0229       Rotation3D r3d(r);
0230       AssignFrom( r3d, r3d( Vector(v.X(),v.Y(),v.Z()) ) );
0231    }
0232 #endif
0233 
0234 public:
0235    /**
0236       Construct transformation from one coordinate system defined by three
0237       points (origin + two axis) to
0238       a new coordinate system defined by other three points (origin + axis)
0239       Scalar version.
0240       @param fr0  point defining origin of original reference system
0241       @param fr1  point defining first axis of original reference system
0242       @param fr2  point defining second axis of original reference system
0243       @param to0  point defining origin of transformed reference system
0244       @param to1  point defining first axis transformed reference system
0245       @param to2  point defining second axis transformed reference system
0246    */
0247    template <typename SCALAR = T, typename std::enable_if<std::is_arithmetic<SCALAR>::value>::type * = nullptr>
0248    Transform3D(const Point &fr0, const Point &fr1, const Point &fr2, const Point &to0, const Point &to1,
0249                const Point &to2)
0250    {
0251       // takes impl. from CLHEP ( E.Chernyaev). To be checked
0252 
0253       Vector x1 = (fr1 - fr0).Unit();
0254       Vector y1 = (fr2 - fr0).Unit();
0255       Vector x2 = (to1 - to0).Unit();
0256       Vector y2 = (to2 - to0).Unit();
0257 
0258       //   C H E C K   A N G L E S
0259 
0260       const T cos1 = x1.Dot(y1);
0261       const T cos2 = x2.Dot(y2);
0262 
0263       if (std::fabs(T(1) - cos1) <= T(0.000001) || std::fabs(T(1) - cos2) <= T(0.000001)) {
0264          std::cerr << "Transform3D: Error : zero angle between axes" << std::endl;
0265          SetIdentity();
0266       } else {
0267          if (std::fabs(cos1 - cos2) > T(0.000001)) {
0268             std::cerr << "Transform3D: Warning: angles between axes are not equal" << std::endl;
0269          }
0270 
0271          //   F I N D   R O T A T I O N   M A T R I X
0272 
0273          Vector z1 = (x1.Cross(y1)).Unit();
0274          y1        = z1.Cross(x1);
0275 
0276          Vector z2 = (x2.Cross(y2)).Unit();
0277          y2        = z2.Cross(x2);
0278 
0279          T x1x = x1.x();
0280          T x1y = x1.y();
0281          T x1z = x1.z();
0282          T y1x = y1.x();
0283          T y1y = y1.y();
0284          T y1z = y1.z();
0285          T z1x = z1.x();
0286          T z1y = z1.y();
0287          T z1z = z1.z();
0288 
0289          T x2x = x2.x();
0290          T x2y = x2.y();
0291          T x2z = x2.z();
0292          T y2x = y2.x();
0293          T y2y = y2.y();
0294          T y2z = y2.z();
0295          T z2x = z2.x();
0296          T z2y = z2.y();
0297          T z2z = z2.z();
0298 
0299          T detxx = (y1y * z1z - z1y * y1z);
0300          T detxy = -(y1x * z1z - z1x * y1z);
0301          T detxz = (y1x * z1y - z1x * y1y);
0302          T detyx = -(x1y * z1z - z1y * x1z);
0303          T detyy = (x1x * z1z - z1x * x1z);
0304          T detyz = -(x1x * z1y - z1x * x1y);
0305          T detzx = (x1y * y1z - y1y * x1z);
0306          T detzy = -(x1x * y1z - y1x * x1z);
0307          T detzz = (x1x * y1y - y1x * x1y);
0308 
0309          T txx = x2x * detxx + y2x * detyx + z2x * detzx;
0310          T txy = x2x * detxy + y2x * detyy + z2x * detzy;
0311          T txz = x2x * detxz + y2x * detyz + z2x * detzz;
0312          T tyx = x2y * detxx + y2y * detyx + z2y * detzx;
0313          T tyy = x2y * detxy + y2y * detyy + z2y * detzy;
0314          T tyz = x2y * detxz + y2y * detyz + z2y * detzz;
0315          T tzx = x2z * detxx + y2z * detyx + z2z * detzx;
0316          T tzy = x2z * detxy + y2z * detyy + z2z * detzy;
0317          T tzz = x2z * detxz + y2z * detyz + z2z * detzz;
0318 
0319          //   S E T    T R A N S F O R M A T I O N
0320 
0321          T dx1 = fr0.x(), dy1 = fr0.y(), dz1 = fr0.z();
0322          T dx2 = to0.x(), dy2 = to0.y(), dz2 = to0.z();
0323 
0324          SetComponents(txx, txy, txz, dx2 - txx * dx1 - txy * dy1 - txz * dz1, tyx, tyy, tyz,
0325                        dy2 - tyx * dx1 - tyy * dy1 - tyz * dz1, tzx, tzy, tzz, dz2 - tzx * dx1 - tzy * dy1 - tzz * dz1);
0326       }
0327    }
0328 
0329    /**
0330       Construct transformation from one coordinate system defined by three
0331       points (origin + two axis) to
0332       a new coordinate system defined by other three points (origin + axis)
0333       Vectorised version.
0334       @param fr0  point defining origin of original reference system
0335       @param fr1  point defining first axis of original reference system
0336       @param fr2  point defining second axis of original reference system
0337       @param to0  point defining origin of transformed reference system
0338       @param to1  point defining first axis transformed reference system
0339       @param to2  point defining second axis transformed reference system
0340    */
0341    template <typename SCALAR = T, typename std::enable_if<!std::is_arithmetic<SCALAR>::value>::type * = nullptr>
0342    Transform3D(const Point &fr0, const Point &fr1, const Point &fr2, const Point &to0, const Point &to1,
0343                const Point &to2)
0344    {
0345       // takes impl. from CLHEP ( E.Chernyaev). To be checked
0346 
0347       Vector x1 = (fr1 - fr0).Unit();
0348       Vector y1 = (fr2 - fr0).Unit();
0349       Vector x2 = (to1 - to0).Unit();
0350       Vector y2 = (to2 - to0).Unit();
0351 
0352       //   C H E C K   A N G L E S
0353 
0354       const T cos1 = x1.Dot(y1);
0355       const T cos2 = x2.Dot(y2);
0356 
0357       const auto m1 = (abs(T(1) - cos1) <= T(0.000001) || abs(T(1) - cos2) <= T(0.000001));
0358 
0359       const auto m2 = (abs(cos1 - cos2) > T(0.000001));
0360       if (any_of(m2)) {
0361          std::cerr << "Transform3D: Warning: angles between axes are not equal" << std::endl;
0362       }
0363 
0364       //   F I N D   R O T A T I O N   M A T R I X
0365 
0366       Vector z1 = (x1.Cross(y1)).Unit();
0367       y1        = z1.Cross(x1);
0368 
0369       Vector z2 = (x2.Cross(y2)).Unit();
0370       y2        = z2.Cross(x2);
0371 
0372       T x1x = x1.x();
0373       T x1y = x1.y();
0374       T x1z = x1.z();
0375       T y1x = y1.x();
0376       T y1y = y1.y();
0377       T y1z = y1.z();
0378       T z1x = z1.x();
0379       T z1y = z1.y();
0380       T z1z = z1.z();
0381 
0382       T x2x = x2.x();
0383       T x2y = x2.y();
0384       T x2z = x2.z();
0385       T y2x = y2.x();
0386       T y2y = y2.y();
0387       T y2z = y2.z();
0388       T z2x = z2.x();
0389       T z2y = z2.y();
0390       T z2z = z2.z();
0391 
0392       T detxx = (y1y * z1z - z1y * y1z);
0393       T detxy = -(y1x * z1z - z1x * y1z);
0394       T detxz = (y1x * z1y - z1x * y1y);
0395       T detyx = -(x1y * z1z - z1y * x1z);
0396       T detyy = (x1x * z1z - z1x * x1z);
0397       T detyz = -(x1x * z1y - z1x * x1y);
0398       T detzx = (x1y * y1z - y1y * x1z);
0399       T detzy = -(x1x * y1z - y1x * x1z);
0400       T detzz = (x1x * y1y - y1x * x1y);
0401 
0402       T txx = x2x * detxx + y2x * detyx + z2x * detzx;
0403       T txy = x2x * detxy + y2x * detyy + z2x * detzy;
0404       T txz = x2x * detxz + y2x * detyz + z2x * detzz;
0405       T tyx = x2y * detxx + y2y * detyx + z2y * detzx;
0406       T tyy = x2y * detxy + y2y * detyy + z2y * detzy;
0407       T tyz = x2y * detxz + y2y * detyz + z2y * detzz;
0408       T tzx = x2z * detxx + y2z * detyx + z2z * detzx;
0409       T tzy = x2z * detxy + y2z * detyy + z2z * detzy;
0410       T tzz = x2z * detxz + y2z * detyz + z2z * detzz;
0411 
0412       //   S E T    T R A N S F O R M A T I O N
0413 
0414       T dx1 = fr0.x(), dy1 = fr0.y(), dz1 = fr0.z();
0415       T dx2 = to0.x(), dy2 = to0.y(), dz2 = to0.z();
0416 
0417       SetComponents(txx, txy, txz, dx2 - txx * dx1 - txy * dy1 - txz * dz1, tyx, tyy, tyz,
0418                     dy2 - tyx * dx1 - tyy * dy1 - tyz * dz1, tzx, tzy, tzz, dz2 - tzx * dx1 - tzy * dy1 - tzz * dz1);
0419 
0420       if (any_of(m1)) {
0421          std::cerr << "Transform3D: Error : zero angle between axes" << std::endl;
0422          SetIdentity(m1);
0423       }
0424    }
0425 
0426    // use compiler generated copy ctor, copy assignment and destructor
0427 
0428    /**
0429       Construct from a linear algebra matrix of size at least 3x4,
0430       which must support operator()(i,j) to obtain elements (0,0) thru (2,3).
0431       The 3x3 sub-block is assumed to be the rotation part and the translations vector
0432       are described by the 4-th column
0433    */
0434    template<class ForeignMatrix>
0435    explicit constexpr Transform3D(const ForeignMatrix & m) {
0436       SetComponents(m);
0437    }
0438 
0439    /**
0440       Raw constructor from 12 Scalar components
0441    */
0442    Transform3D(T xx, T xy, T xz, T dx, T yx, T yy, T yz, T dy, T zx, T zy, T zz, T dz)
0443    {
0444       SetComponents (xx, xy, xz, dx, yx, yy, yz, dy, zx, zy, zz, dz);
0445    }
0446 
0447 
0448    /**
0449       Construct from a linear algebra matrix of size at least 3x4,
0450       which must support operator()(i,j) to obtain elements (0,0) thru (2,3).
0451       The 3x3 sub-block is assumed to be the rotation part and the translations vector
0452       are described by the 4-th column
0453    */
0454    template <class ForeignMatrix>
0455    Transform3D<T> &operator=(const ForeignMatrix &m)
0456    {
0457       SetComponents(m);
0458       return *this;
0459    }
0460 
0461 
0462    // ======== Components ==============
0463 
0464 
0465    /**
0466       Set the 12 matrix components given an iterator to the start of
0467       the desired data, and another to the end (12 past start).
0468    */
0469    template<class IT>
0470    void SetComponents(IT begin, IT end) {
0471       for (int i = 0; i <12; ++i) {
0472          fM[i] = *begin;
0473          ++begin;
0474       }
0475       (void)end;
0476       assert (end==begin);
0477    }
0478 
0479    /**
0480       Get the 12 matrix components into data specified by an iterator begin
0481       and another to the end of the desired data (12 past start).
0482    */
0483    template<class IT>
0484    void GetComponents(IT begin, IT end) const {
0485       for (int i = 0; i <12; ++i) {
0486          *begin = fM[i];
0487          ++begin;
0488       }
0489       (void)end;
0490       assert (end==begin);
0491    }
0492 
0493    /**
0494       Get the 12 matrix components into data specified by an iterator begin
0495    */
0496    template<class IT>
0497    void GetComponents(IT begin) const {
0498       std::copy(fM, fM + 12, begin);
0499    }
0500 
0501    /**
0502       Set components from a linear algebra matrix of size at least 3x4,
0503       which must support operator()(i,j) to obtain elements (0,0) thru (2,3).
0504       The 3x3 sub-block is assumed to be the rotation part and the translations vector
0505       are described by the 4-th column
0506    */
0507    template<class ForeignMatrix>
0508    void
0509    SetTransformMatrix (const ForeignMatrix & m) {
0510       fM[kXX]=m(0,0);  fM[kXY]=m(0,1);  fM[kXZ]=m(0,2); fM[kDX]=m(0,3);
0511       fM[kYX]=m(1,0);  fM[kYY]=m(1,1);  fM[kYZ]=m(1,2); fM[kDY]=m(1,3);
0512       fM[kZX]=m(2,0);  fM[kZY]=m(2,1);  fM[kZZ]=m(2,2); fM[kDZ]=m(2,3);
0513    }
0514 
0515    /**
0516       Get components into a linear algebra matrix of size at least 3x4,
0517       which must support operator()(i,j) for write access to elements
0518       (0,0) thru (2,3).
0519    */
0520    template<class ForeignMatrix>
0521    void
0522    GetTransformMatrix (ForeignMatrix & m) const {
0523       m(0,0)=fM[kXX];  m(0,1)=fM[kXY];  m(0,2)=fM[kXZ];  m(0,3)=fM[kDX];
0524       m(1,0)=fM[kYX];  m(1,1)=fM[kYY];  m(1,2)=fM[kYZ];  m(1,3)=fM[kDY];
0525       m(2,0)=fM[kZX];  m(2,1)=fM[kZY];  m(2,2)=fM[kZZ];  m(2,3)=fM[kDZ];
0526    }
0527 
0528 
0529    /**
0530       Set the components from 12 scalars
0531    */
0532    void SetComponents(T xx, T xy, T xz, T dx, T yx, T yy, T yz, T dy, T zx, T zy, T zz, T dz)
0533    {
0534       fM[kXX]=xx;  fM[kXY]=xy;  fM[kXZ]=xz;  fM[kDX]=dx;
0535       fM[kYX]=yx;  fM[kYY]=yy;  fM[kYZ]=yz;  fM[kDY]=dy;
0536       fM[kZX]=zx;  fM[kZY]=zy;  fM[kZZ]=zz;  fM[kDZ]=dz;
0537    }
0538 
0539    /**
0540       Get the components into 12 scalars
0541    */
0542    void GetComponents(T &xx, T &xy, T &xz, T &dx, T &yx, T &yy, T &yz, T &dy, T &zx, T &zy, T &zz, T &dz) const
0543    {
0544       xx=fM[kXX];  xy=fM[kXY];  xz=fM[kXZ];  dx=fM[kDX];
0545       yx=fM[kYX];  yy=fM[kYY];  yz=fM[kYZ];  dy=fM[kDY];
0546       zx=fM[kZX];  zy=fM[kZY];  zz=fM[kZZ];  dz=fM[kDZ];
0547    }
0548 
0549 
0550    /**
0551       Get the rotation and translation vector representing the 3D transformation
0552       in any rotation and any vector (the Translation class could also be used)
0553    */
0554    template<class AnyRotation, class V>
0555    void GetDecomposition(AnyRotation &r, V &v) const {
0556       GetRotation(r);
0557       GetTranslation(v);
0558    }
0559 
0560 
0561    /**
0562       Get the rotation and translation vector representing the 3D transformation
0563    */
0564    void GetDecomposition(Rotation3D &r, Vector &v) const {
0565       GetRotation(r);
0566       GetTranslation(v);
0567    }
0568 
0569    /**
0570       Get the 3D rotation representing the 3D transformation
0571    */
0572    Rotation3D Rotation() const {
0573       return Rotation3D( fM[kXX], fM[kXY], fM[kXZ],
0574                          fM[kYX], fM[kYY], fM[kYZ],
0575                          fM[kZX], fM[kZY], fM[kZZ] );
0576    }
0577 
0578    /**
0579       Get the rotation representing the 3D transformation
0580    */
0581    template <class AnyRotation>
0582    AnyRotation Rotation() const {
0583       return AnyRotation(Rotation3D(fM[kXX], fM[kXY], fM[kXZ], fM[kYX], fM[kYY], fM[kYZ], fM[kZX], fM[kZY], fM[kZZ]));
0584    }
0585 
0586    /**
0587       Get the  rotation (any type) representing the 3D transformation
0588    */
0589    template <class AnyRotation>
0590    void GetRotation(AnyRotation &r) const {
0591       r = Rotation();
0592    }
0593 
0594    /**
0595       Get the translation representing the 3D transformation in a Cartesian vector
0596    */
0597    Translation3D<T> Translation() const { return Translation3D<T>(fM[kDX], fM[kDY], fM[kDZ]); }
0598 
0599    /**
0600       Get the translation representing the 3D transformation in any vector
0601       which implements the SetXYZ method
0602    */
0603    template <class AnyVector>
0604    void GetTranslation(AnyVector &v) const {
0605       v.SetXYZ(fM[kDX], fM[kDY], fM[kDZ]);
0606    }
0607 
0608 
0609 
0610    // operations on points and vectors
0611 
0612    /**
0613       Transformation operation for Position Vector in Cartesian coordinate
0614       For a Position Vector first a rotation and then a translation is applied
0615    */
0616    Point operator() (const Point & p) const {
0617       return Point ( fM[kXX]*p.X() + fM[kXY]*p.Y() + fM[kXZ]*p.Z() + fM[kDX],
0618                      fM[kYX]*p.X() + fM[kYY]*p.Y() + fM[kYZ]*p.Z() + fM[kDY],
0619                      fM[kZX]*p.X() + fM[kZY]*p.Y() + fM[kZZ]*p.Z() + fM[kDZ] );
0620    }
0621 
0622 
0623    /**
0624       Transformation operation for Displacement Vectors in Cartesian coordinate
0625       For the Displacement Vectors only the rotation applies - no translations
0626    */
0627    Vector operator() (const Vector & v) const {
0628       return Vector( fM[kXX]*v.X() + fM[kXY]*v.Y() + fM[kXZ]*v.Z() ,
0629                      fM[kYX]*v.X() + fM[kYY]*v.Y() + fM[kYZ]*v.Z() ,
0630                      fM[kZX]*v.X() + fM[kZY]*v.Y() + fM[kZZ]*v.Z()  );
0631    }
0632 
0633 
0634    /**
0635       Transformation operation for Position Vector in any coordinate system
0636    */
0637    template <class CoordSystem>
0638    PositionVector3D<CoordSystem> operator()(const PositionVector3D<CoordSystem> &p) const
0639    {
0640       return PositionVector3D<CoordSystem>(operator()(Point(p)));
0641    }
0642    /**
0643       Transformation operation for Position Vector in any coordinate system
0644    */
0645    template <class CoordSystem>
0646    PositionVector3D<CoordSystem> operator*(const PositionVector3D<CoordSystem> &v) const
0647    {
0648       return operator()(v);
0649    }
0650 
0651    /**
0652       Transformation operation for Displacement Vector in any coordinate system
0653    */
0654    template<class CoordSystem >
0655    DisplacementVector3D<CoordSystem> operator() (const DisplacementVector3D <CoordSystem> & v) const {
0656       return DisplacementVector3D<CoordSystem>(operator()(Vector(v)));
0657    }
0658    /**
0659       Transformation operation for Displacement Vector in any coordinate system
0660    */
0661    template <class CoordSystem>
0662    DisplacementVector3D<CoordSystem> operator*(const DisplacementVector3D<CoordSystem> &v) const
0663    {
0664       return operator()(v);
0665    }
0666 
0667    /**
0668       Directly apply the inverse affine transformation on vectors.
0669       Avoids having to calculate the inverse as an intermediate result.
0670       This is possible since the inverse of a rotation is its transpose.
0671    */
0672    Vector ApplyInverse(const Vector &v) const
0673    {
0674       return Vector(fM[kXX] * v.X() + fM[kYX] * v.Y() + fM[kZX] * v.Z(),
0675                     fM[kXY] * v.X() + fM[kYY] * v.Y() + fM[kZY] * v.Z(),
0676                     fM[kXZ] * v.X() + fM[kYZ] * v.Y() + fM[kZZ] * v.Z());
0677    }
0678 
0679    /**
0680       Directly apply the inverse affine transformation on points
0681       (first inverse translation then inverse rotation).
0682       Avoids having to calculate the inverse as an intermediate result.
0683       This is possible since the inverse of a rotation is its transpose.
0684    */
0685    Point ApplyInverse(const Point &p) const
0686    {
0687       Point tmp(p.X() - fM[kDX], p.Y() - fM[kDY], p.Z() - fM[kDZ]);
0688       return Point(fM[kXX] * tmp.X() + fM[kYX] * tmp.Y() + fM[kZX] * tmp.Z(),
0689                    fM[kXY] * tmp.X() + fM[kYY] * tmp.Y() + fM[kZY] * tmp.Z(),
0690                    fM[kXZ] * tmp.X() + fM[kYZ] * tmp.Y() + fM[kZZ] * tmp.Z());
0691    }
0692 
0693    /**
0694       Directly apply the inverse affine transformation on an arbitrary
0695       coordinate-system point.
0696       Involves casting to Point(p) type.
0697    */
0698    template <class CoordSystem>
0699    PositionVector3D<CoordSystem> ApplyInverse(const PositionVector3D<CoordSystem> &p) const
0700    {
0701       return PositionVector3D<CoordSystem>(ApplyInverse(Point(p)));
0702    }
0703 
0704    /**
0705       Directly apply the inverse affine transformation on an arbitrary
0706       coordinate-system vector.
0707       Involves casting to Vector(p) type.
0708    */
0709    template <class CoordSystem>
0710    DisplacementVector3D<CoordSystem> ApplyInverse(const DisplacementVector3D<CoordSystem> &p) const
0711    {
0712       return DisplacementVector3D<CoordSystem>(ApplyInverse(Vector(p)));
0713    }
0714 
0715    /**
0716       Transformation operation for points between different coordinate system tags
0717    */
0718    template <class CoordSystem, class Tag1, class Tag2>
0719    void Transform(const PositionVector3D<CoordSystem, Tag1> &p1, PositionVector3D<CoordSystem, Tag2> &p2) const
0720    {
0721       const Point xyzNew = operator()(Point(p1.X(), p1.Y(), p1.Z()));
0722       p2.SetXYZ( xyzNew.X(), xyzNew.Y(), xyzNew.Z() );
0723    }
0724 
0725 
0726    /**
0727       Transformation operation for Displacement Vector of different coordinate systems
0728    */
0729    template <class CoordSystem, class Tag1, class Tag2>
0730    void Transform(const DisplacementVector3D<CoordSystem, Tag1> &v1, DisplacementVector3D<CoordSystem, Tag2> &v2) const
0731    {
0732       const Vector xyzNew = operator()(Vector(v1.X(), v1.Y(), v1.Z()));
0733       v2.SetXYZ( xyzNew.X(), xyzNew.Y(), xyzNew.Z() );
0734    }
0735 
0736    /**
0737       Transformation operation for a Lorentz Vector in any  coordinate system
0738    */
0739    template <class CoordSystem >
0740    LorentzVector<CoordSystem> operator() (const LorentzVector<CoordSystem> & q) const {
0741       const Vector xyzNew = operator()(Vector(q.Vect()));
0742       return LorentzVector<CoordSystem>(xyzNew.X(), xyzNew.Y(), xyzNew.Z(), q.E());
0743    }
0744    /**
0745       Transformation operation for a Lorentz Vector in any  coordinate system
0746    */
0747    template <class CoordSystem>
0748    LorentzVector<CoordSystem> operator*(const LorentzVector<CoordSystem> &q) const
0749    {
0750       return operator()(q);
0751    }
0752 
0753    /**
0754       Transformation on a 3D plane
0755    */
0756    template <typename TYPE>
0757    Plane3D<TYPE> operator()(const Plane3D<TYPE> &plane) const
0758    {
0759       // transformations on a 3D plane
0760       const auto n = plane.Normal();
0761       // take a point on the plane. Use origin projection on the plane
0762       // ( -ad, -bd, -cd) if (a**2 + b**2 + c**2 ) = 1
0763       const auto d = plane.HesseDistance();
0764       Point p(-d * n.X(), -d * n.Y(), -d * n.Z());
0765       return Plane3D<TYPE>(operator()(n), operator()(p));
0766    }
0767 
0768    /// Multiplication operator for 3D plane
0769    template <typename TYPE>
0770    Plane3D<TYPE> operator*(const Plane3D<TYPE> &plane) const
0771    {
0772       return operator()(plane);
0773    }
0774 
0775    // skip transformation for arbitrary vectors - not really defined if point or displacement vectors
0776 
0777    /**
0778       multiply (combine) with another transformation in place
0779    */
0780    inline Transform3D<T> &operator*=(const Transform3D<T> &t);
0781 
0782    /**
0783       multiply (combine) two transformations
0784    */
0785    inline Transform3D<T> operator*(const Transform3D<T> &t) const;
0786 
0787    /**
0788        Invert the transformation in place (scalar)
0789    */
0790    template <typename SCALAR = T, typename std::enable_if<std::is_arithmetic<SCALAR>::value>::type * = nullptr>
0791    void Invert()
0792    {
0793       //
0794       // Name: Transform3D::inverse                     Date:    24.09.96
0795       // Author: E.Chernyaev (IHEP/Protvino)            Revised:
0796       //
0797       // Function: Find inverse affine transformation.
0798 
0799       T detxx = fM[kYY] * fM[kZZ] - fM[kYZ] * fM[kZY];
0800       T detxy = fM[kYX] * fM[kZZ] - fM[kYZ] * fM[kZX];
0801       T detxz = fM[kYX] * fM[kZY] - fM[kYY] * fM[kZX];
0802       T det   = fM[kXX] * detxx - fM[kXY] * detxy + fM[kXZ] * detxz;
0803       if (det == T(0)) {
0804          std::cerr << "Transform3D::inverse error: zero determinant" << std::endl;
0805          return;
0806       }
0807       det = T(1) / det;
0808       detxx *= det;
0809       detxy *= det;
0810       detxz *= det;
0811       T detyx = (fM[kXY] * fM[kZZ] - fM[kXZ] * fM[kZY]) * det;
0812       T detyy = (fM[kXX] * fM[kZZ] - fM[kXZ] * fM[kZX]) * det;
0813       T detyz = (fM[kXX] * fM[kZY] - fM[kXY] * fM[kZX]) * det;
0814       T detzx = (fM[kXY] * fM[kYZ] - fM[kXZ] * fM[kYY]) * det;
0815       T detzy = (fM[kXX] * fM[kYZ] - fM[kXZ] * fM[kYX]) * det;
0816       T detzz = (fM[kXX] * fM[kYY] - fM[kXY] * fM[kYX]) * det;
0817       SetComponents(detxx, -detyx, detzx, -detxx * fM[kDX] + detyx * fM[kDY] - detzx * fM[kDZ], -detxy, detyy, -detzy,
0818                     detxy * fM[kDX] - detyy * fM[kDY] + detzy * fM[kDZ], detxz, -detyz, detzz,
0819                     -detxz * fM[kDX] + detyz * fM[kDY] - detzz * fM[kDZ]);
0820    }
0821 
0822    /**
0823        Invert the transformation in place (vectorised)
0824    */
0825    template <typename SCALAR = T, typename std::enable_if<!std::is_arithmetic<SCALAR>::value>::type * = nullptr>
0826    void Invert()
0827    {
0828       //
0829       // Name: Transform3D::inverse                     Date:    24.09.96
0830       // Author: E.Chernyaev (IHEP/Protvino)            Revised:
0831       //
0832       // Function: Find inverse affine transformation.
0833 
0834       T          detxx    = fM[kYY] * fM[kZZ] - fM[kYZ] * fM[kZY];
0835       T          detxy    = fM[kYX] * fM[kZZ] - fM[kYZ] * fM[kZX];
0836       T          detxz    = fM[kYX] * fM[kZY] - fM[kYY] * fM[kZX];
0837       T          det      = fM[kXX] * detxx - fM[kXY] * detxy + fM[kXZ] * detxz;
0838       const auto detZmask = (det == T(0));
0839       if (any_of(detZmask)) {
0840          std::cerr << "Transform3D::inverse error: zero determinant" << std::endl;
0841          where(detZmask, det) = T(1);
0842       }
0843       det = T(1) / det;
0844       detxx *= det;
0845       detxy *= det;
0846       detxz *= det;
0847       T detyx = (fM[kXY] * fM[kZZ] - fM[kXZ] * fM[kZY]) * det;
0848       T detyy = (fM[kXX] * fM[kZZ] - fM[kXZ] * fM[kZX]) * det;
0849       T detyz = (fM[kXX] * fM[kZY] - fM[kXY] * fM[kZX]) * det;
0850       T detzx = (fM[kXY] * fM[kYZ] - fM[kXZ] * fM[kYY]) * det;
0851       T detzy = (fM[kXX] * fM[kYZ] - fM[kXZ] * fM[kYX]) * det;
0852       T detzz = (fM[kXX] * fM[kYY] - fM[kXY] * fM[kYX]) * det;
0853       // Set det=0 cases to 0
0854       if (any_of(detZmask)) {
0855          where(detZmask, detxx) = T(0);
0856          where(detZmask, detxy) = T(0);
0857          where(detZmask, detxz) = T(0);
0858          where(detZmask, detyx) = T(0);
0859          where(detZmask, detyy) = T(0);
0860          where(detZmask, detyz) = T(0);
0861          where(detZmask, detzx) = T(0);
0862          where(detZmask, detzy) = T(0);
0863          where(detZmask, detzz) = T(0);
0864       }
0865       // set final components
0866       SetComponents(detxx, -detyx, detzx, -detxx * fM[kDX] + detyx * fM[kDY] - detzx * fM[kDZ], -detxy, detyy, -detzy,
0867                     detxy * fM[kDX] - detyy * fM[kDY] + detzy * fM[kDZ], detxz, -detyz, detzz,
0868                     -detxz * fM[kDX] + detyz * fM[kDY] - detzz * fM[kDZ]);
0869    }
0870 
0871    /**
0872       Return the inverse of the transformation.
0873    */
0874    Transform3D<T> Inverse() const
0875    {
0876       Transform3D<T> t(*this);
0877       t.Invert();
0878       return t;
0879    }
0880 
0881    /**
0882       Equality operator. Check equality for each element
0883       To do: use T tolerance
0884    */
0885    bool operator==(const Transform3D<T> &rhs) const
0886    {
0887       return (fM[0] == rhs.fM[0] && fM[1] == rhs.fM[1] && fM[2] == rhs.fM[2] && fM[3] == rhs.fM[3] &&
0888               fM[4] == rhs.fM[4] && fM[5] == rhs.fM[5] && fM[6] == rhs.fM[6] && fM[7] == rhs.fM[7] &&
0889               fM[8] == rhs.fM[8] && fM[9] == rhs.fM[9] && fM[10] == rhs.fM[10] && fM[11] == rhs.fM[11]);
0890    }
0891 
0892    /**
0893       Inequality operator. Check equality for each element
0894       To do: use T tolerance
0895    */
0896    bool operator!=(const Transform3D<T> &rhs) const { return !operator==(rhs); }
0897 
0898 protected:
0899 
0900    /**
0901       make transformation from first a rotation then a translation
0902    */
0903    void AssignFrom(const Rotation3D &r, const Vector &v)
0904    {
0905       // assignment  from rotation + translation
0906 
0907       T rotData[9];
0908       r.GetComponents(rotData, rotData + 9);
0909       // first raw
0910       for (int i = 0; i < 3; ++i) fM[i] = rotData[i];
0911       // second raw
0912       for (int i = 0; i < 3; ++i) fM[kYX + i] = rotData[3 + i];
0913       // third raw
0914       for (int i = 0; i < 3; ++i) fM[kZX + i] = rotData[6 + i];
0915 
0916       // translation data
0917       T vecData[3];
0918       v.GetCoordinates(vecData, vecData + 3);
0919       fM[kDX] = vecData[0];
0920       fM[kDY] = vecData[1];
0921       fM[kDZ] = vecData[2];
0922    }
0923 
0924    /**
0925       make transformation from only rotations (zero translation)
0926    */
0927    void AssignFrom(const Rotation3D &r)
0928    {
0929       // assign from only a rotation  (null translation)
0930       T rotData[9];
0931       r.GetComponents(rotData, rotData + 9);
0932       for (int i = 0; i < 3; ++i) {
0933          for (int j = 0; j < 3; ++j) fM[4 * i + j] = rotData[3 * i + j];
0934          // empty vector data
0935          fM[4 * i + 3] = T(0);
0936       }
0937    }
0938 
0939    /**
0940       make transformation from only translation (identity rotations)
0941    */
0942    void AssignFrom(const Vector &v)
0943    {
0944       // assign from a translation only (identity rotations)
0945       fM[kXX] = T(1);
0946       fM[kXY] = T(0);
0947       fM[kXZ] = T(0);
0948       fM[kDX] = v.X();
0949       fM[kYX] = T(0);
0950       fM[kYY] = T(1);
0951       fM[kYZ] = T(0);
0952       fM[kDY] = v.Y();
0953       fM[kZX] = T(0);
0954       fM[kZY] = T(0);
0955       fM[kZZ] = T(1);
0956       fM[kDZ] = v.Z();
0957    }
0958 
0959    /**
0960       Set identity transformation (identity rotation , zero translation)
0961    */
0962    void SetIdentity()
0963    {
0964       // set identity ( identity rotation and zero translation)
0965       fM[kXX] = T(1);
0966       fM[kXY] = T(0);
0967       fM[kXZ] = T(0);
0968       fM[kDX] = T(0);
0969       fM[kYX] = T(0);
0970       fM[kYY] = T(1);
0971       fM[kYZ] = T(0);
0972       fM[kDY] = T(0);
0973       fM[kZX] = T(0);
0974       fM[kZY] = T(0);
0975       fM[kZZ] = T(1);
0976       fM[kDZ] = T(0);
0977    }
0978 
0979    /**
0980       Set identity transformation (identity rotation , zero translation)
0981       vectorised version that sets using a mask
0982    */
0983    template <typename SCALAR = T, typename std::enable_if<!std::is_arithmetic<SCALAR>::value>::type * = nullptr>
0984    void SetIdentity(const typename SCALAR::mask_type m)
0985    {
0986       // set identity ( identity rotation and zero translation)
0987       where(m, fM[kXX]) = T(1);
0988       where(m, fM[kXY]) = T(0);
0989       where(m, fM[kXZ]) = T(0);
0990       where(m, fM[kDX]) = T(0);
0991       where(m, fM[kYX]) = T(0);
0992       where(m, fM[kYY]) = T(1);
0993       where(m, fM[kYZ]) = T(0);
0994       where(m, fM[kDY]) = T(0);
0995       where(m, fM[kZX]) = T(0);
0996       where(m, fM[kZY]) = T(0);
0997       where(m, fM[kZZ]) = T(1);
0998       where(m, fM[kDZ]) = T(0);
0999    }
1000 
1001 private:
1002    T fM[12]; // transformation elements (3x4 matrix)
1003 };
1004 
1005 
1006 
1007 
1008 // inline functions (combination of transformations)
1009 
1010 template <class T>
1011 inline Transform3D<T> &Transform3D<T>::operator*=(const Transform3D<T> &t)
1012 {
1013    // combination of transformations
1014 
1015    SetComponents(fM[kXX]*t.fM[kXX]+fM[kXY]*t.fM[kYX]+fM[kXZ]*t.fM[kZX],
1016                  fM[kXX]*t.fM[kXY]+fM[kXY]*t.fM[kYY]+fM[kXZ]*t.fM[kZY],
1017                  fM[kXX]*t.fM[kXZ]+fM[kXY]*t.fM[kYZ]+fM[kXZ]*t.fM[kZZ],
1018                  fM[kXX]*t.fM[kDX]+fM[kXY]*t.fM[kDY]+fM[kXZ]*t.fM[kDZ]+fM[kDX],
1019 
1020                  fM[kYX]*t.fM[kXX]+fM[kYY]*t.fM[kYX]+fM[kYZ]*t.fM[kZX],
1021                  fM[kYX]*t.fM[kXY]+fM[kYY]*t.fM[kYY]+fM[kYZ]*t.fM[kZY],
1022                  fM[kYX]*t.fM[kXZ]+fM[kYY]*t.fM[kYZ]+fM[kYZ]*t.fM[kZZ],
1023                  fM[kYX]*t.fM[kDX]+fM[kYY]*t.fM[kDY]+fM[kYZ]*t.fM[kDZ]+fM[kDY],
1024 
1025                  fM[kZX]*t.fM[kXX]+fM[kZY]*t.fM[kYX]+fM[kZZ]*t.fM[kZX],
1026                  fM[kZX]*t.fM[kXY]+fM[kZY]*t.fM[kYY]+fM[kZZ]*t.fM[kZY],
1027                  fM[kZX]*t.fM[kXZ]+fM[kZY]*t.fM[kYZ]+fM[kZZ]*t.fM[kZZ],
1028                  fM[kZX]*t.fM[kDX]+fM[kZY]*t.fM[kDY]+fM[kZZ]*t.fM[kDZ]+fM[kDZ]);
1029 
1030    return *this;
1031 }
1032 
1033 template <class T>
1034 inline Transform3D<T> Transform3D<T>::operator*(const Transform3D<T> &t) const
1035 {
1036    // combination of transformations
1037 
1038    return Transform3D<T>(fM[kXX] * t.fM[kXX] + fM[kXY] * t.fM[kYX] + fM[kXZ] * t.fM[kZX],
1039                          fM[kXX] * t.fM[kXY] + fM[kXY] * t.fM[kYY] + fM[kXZ] * t.fM[kZY],
1040                          fM[kXX] * t.fM[kXZ] + fM[kXY] * t.fM[kYZ] + fM[kXZ] * t.fM[kZZ],
1041                          fM[kXX] * t.fM[kDX] + fM[kXY] * t.fM[kDY] + fM[kXZ] * t.fM[kDZ] + fM[kDX],
1042 
1043                          fM[kYX] * t.fM[kXX] + fM[kYY] * t.fM[kYX] + fM[kYZ] * t.fM[kZX],
1044                          fM[kYX] * t.fM[kXY] + fM[kYY] * t.fM[kYY] + fM[kYZ] * t.fM[kZY],
1045                          fM[kYX] * t.fM[kXZ] + fM[kYY] * t.fM[kYZ] + fM[kYZ] * t.fM[kZZ],
1046                          fM[kYX] * t.fM[kDX] + fM[kYY] * t.fM[kDY] + fM[kYZ] * t.fM[kDZ] + fM[kDY],
1047 
1048                          fM[kZX] * t.fM[kXX] + fM[kZY] * t.fM[kYX] + fM[kZZ] * t.fM[kZX],
1049                          fM[kZX] * t.fM[kXY] + fM[kZY] * t.fM[kYY] + fM[kZZ] * t.fM[kZY],
1050                          fM[kZX] * t.fM[kXZ] + fM[kZY] * t.fM[kYZ] + fM[kZZ] * t.fM[kZZ],
1051                          fM[kZX] * t.fM[kDX] + fM[kZY] * t.fM[kDY] + fM[kZZ] * t.fM[kDZ] + fM[kDZ]);
1052 }
1053 
1054 
1055 
1056 
1057 //--- global functions resulting in Transform3D -------
1058 
1059 
1060 // ------ combination of a  translation (first)  and a rotation ------
1061 
1062 
1063 /**
1064    combine a translation and a rotation to give a transform3d
1065    First the translation then the rotation
1066  */
1067 template <class T>
1068 inline Transform3D<T> operator*(const Rotation3D &r, const Translation3D<T> &t)
1069 {
1070    return Transform3D<T>(r, r(t.Vect()));
1071 }
1072 template <class T>
1073 inline Transform3D<T> operator*(const RotationX &r, const Translation3D<T> &t)
1074 {
1075    Rotation3D r3(r);
1076    return Transform3D<T>(r3, r3(t.Vect()));
1077 }
1078 template <class T>
1079 inline Transform3D<T> operator*(const RotationY &r, const Translation3D<T> &t)
1080 {
1081    Rotation3D r3(r);
1082    return Transform3D<T>(r3, r3(t.Vect()));
1083 }
1084 template <class T>
1085 inline Transform3D<T> operator*(const RotationZ &r, const Translation3D<T> &t)
1086 {
1087    Rotation3D r3(r);
1088    return Transform3D<T>(r3, r3(t.Vect()));
1089 }
1090 template <class T>
1091 inline Transform3D<T> operator*(const RotationZYX &r, const Translation3D<T> &t)
1092 {
1093    Rotation3D r3(r);
1094    return Transform3D<T>(r3, r3(t.Vect()));
1095 }
1096 template <class T>
1097 inline Transform3D<T> operator*(const AxisAngle &r, const Translation3D<T> &t)
1098 {
1099    Rotation3D r3(r);
1100    return Transform3D<T>(r3, r3(t.Vect()));
1101 }
1102 template <class T>
1103 inline Transform3D<T> operator*(const EulerAngles &r, const Translation3D<T> &t)
1104 {
1105    Rotation3D r3(r);
1106    return Transform3D<T>(r3, r3(t.Vect()));
1107 }
1108 template <class T>
1109 inline Transform3D<T> operator*(const Quaternion &r, const Translation3D<T> &t)
1110 {
1111    Rotation3D r3(r);
1112    return Transform3D<T>(r3, r3(t.Vect()));
1113 }
1114 
1115 // ------ combination of a  rotation (first)  and then a translation ------
1116 
1117 /**
1118    combine a rotation and a translation to give a transform3d
1119    First a rotation then the translation
1120  */
1121 template <class T>
1122 inline Transform3D<T> operator*(const Translation3D<T> &t, const Rotation3D &r)
1123 {
1124    return Transform3D<T>(r, t.Vect());
1125 }
1126 template <class T>
1127 inline Transform3D<T> operator*(const Translation3D<T> &t, const RotationX &r)
1128 {
1129    return Transform3D<T>(Rotation3D(r), t.Vect());
1130 }
1131 template <class T>
1132 inline Transform3D<T> operator*(const Translation3D<T> &t, const RotationY &r)
1133 {
1134    return Transform3D<T>(Rotation3D(r), t.Vect());
1135 }
1136 template <class T>
1137 inline Transform3D<T> operator*(const Translation3D<T> &t, const RotationZ &r)
1138 {
1139    return Transform3D<T>(Rotation3D(r), t.Vect());
1140 }
1141 template <class T>
1142 inline Transform3D<T> operator*(const Translation3D<T> &t, const RotationZYX &r)
1143 {
1144    return Transform3D<T>(Rotation3D(r), t.Vect());
1145 }
1146 template <class T>
1147 inline Transform3D<T> operator*(const Translation3D<T> &t, const EulerAngles &r)
1148 {
1149    return Transform3D<T>(Rotation3D(r), t.Vect());
1150 }
1151 template <class T>
1152 inline Transform3D<T> operator*(const Translation3D<T> &t, const Quaternion &r)
1153 {
1154    return Transform3D<T>(Rotation3D(r), t.Vect());
1155 }
1156 template <class T>
1157 inline Transform3D<T> operator*(const Translation3D<T> &t, const AxisAngle &r)
1158 {
1159    return Transform3D<T>(Rotation3D(r), t.Vect());
1160 }
1161 
1162 // ------ combination of a Transform3D and a pure translation------
1163 
1164 /**
1165    combine a transformation and a translation to give a transform3d
1166    First the translation then the transform3D
1167  */
1168 template <class T>
1169 inline Transform3D<T> operator*(const Transform3D<T> &t, const Translation3D<T> &d)
1170 {
1171    Rotation3D r = t.Rotation();
1172    return Transform3D<T>(r, r(d.Vect()) + t.Translation().Vect());
1173 }
1174 
1175 /**
1176    combine a translation and a transformation to give a transform3d
1177    First the transformation then the translation
1178  */
1179 template <class T>
1180 inline Transform3D<T> operator*(const Translation3D<T> &d, const Transform3D<T> &t)
1181 {
1182    return Transform3D<T>(t.Rotation(), t.Translation().Vect() + d.Vect());
1183 }
1184 
1185 // ------ combination of a Transform3D and any rotation------
1186 
1187 
1188 /**
1189    combine a transformation and a rotation to give a transform3d
1190    First the rotation then the transform3D
1191  */
1192 template <class T>
1193 inline Transform3D<T> operator*(const Transform3D<T> &t, const Rotation3D &r)
1194 {
1195    return Transform3D<T>(t.Rotation() * r, t.Translation());
1196 }
1197 template <class T>
1198 inline Transform3D<T> operator*(const Transform3D<T> &t, const RotationX &r)
1199 {
1200    return Transform3D<T>(t.Rotation() * r, t.Translation());
1201 }
1202 template <class T>
1203 inline Transform3D<T> operator*(const Transform3D<T> &t, const RotationY &r)
1204 {
1205    return Transform3D<T>(t.Rotation() * r, t.Translation());
1206 }
1207 template <class T>
1208 inline Transform3D<T> operator*(const Transform3D<T> &t, const RotationZ &r)
1209 {
1210    return Transform3D<T>(t.Rotation() * r, t.Translation());
1211 }
1212 template <class T>
1213 inline Transform3D<T> operator*(const Transform3D<T> &t, const RotationZYX &r)
1214 {
1215    return Transform3D<T>(t.Rotation() * r, t.Translation());
1216 }
1217 template <class T>
1218 inline Transform3D<T> operator*(const Transform3D<T> &t, const EulerAngles &r)
1219 {
1220    return Transform3D<T>(t.Rotation() * r, t.Translation());
1221 }
1222 template <class T>
1223 inline Transform3D<T> operator*(const Transform3D<T> &t, const AxisAngle &r)
1224 {
1225    return Transform3D<T>(t.Rotation() * r, t.Translation());
1226 }
1227 template <class T>
1228 inline Transform3D<T> operator*(const Transform3D<T> &t, const Quaternion &r)
1229 {
1230    return Transform3D<T>(t.Rotation() * r, t.Translation());
1231 }
1232 
1233 
1234 
1235 /**
1236    combine a rotation and a transformation to give a transform3d
1237    First the transformation then the rotation
1238  */
1239 template <class T>
1240 inline Transform3D<T> operator*(const Rotation3D &r, const Transform3D<T> &t)
1241 {
1242    return Transform3D<T>(r * t.Rotation(), r * t.Translation().Vect());
1243 }
1244 template <class T>
1245 inline Transform3D<T> operator*(const RotationX &r, const Transform3D<T> &t)
1246 {
1247    Rotation3D r3d(r);
1248    return Transform3D<T>(r3d * t.Rotation(), r3d * t.Translation().Vect());
1249 }
1250 template <class T>
1251 inline Transform3D<T> operator*(const RotationY &r, const Transform3D<T> &t)
1252 {
1253    Rotation3D r3d(r);
1254    return Transform3D<T>(r3d * t.Rotation(), r3d * t.Translation().Vect());
1255 }
1256 template <class T>
1257 inline Transform3D<T> operator*(const RotationZ &r, const Transform3D<T> &t)
1258 {
1259    Rotation3D r3d(r);
1260    return Transform3D<T>(r3d * t.Rotation(), r3d * t.Translation().Vect());
1261 }
1262 template <class T>
1263 inline Transform3D<T> operator*(const RotationZYX &r, const Transform3D<T> &t)
1264 {
1265    Rotation3D r3d(r);
1266    return Transform3D<T>(r3d * t.Rotation(), r3d * t.Translation().Vect());
1267 }
1268 template <class T>
1269 inline Transform3D<T> operator*(const EulerAngles &r, const Transform3D<T> &t)
1270 {
1271    Rotation3D r3d(r);
1272    return Transform3D<T>(r3d * t.Rotation(), r3d * t.Translation().Vect());
1273 }
1274 template <class T>
1275 inline Transform3D<T> operator*(const AxisAngle &r, const Transform3D<T> &t)
1276 {
1277    Rotation3D r3d(r);
1278    return Transform3D<T>(r3d * t.Rotation(), r3d * t.Translation().Vect());
1279 }
1280 template <class T>
1281 inline Transform3D<T> operator*(const Quaternion &r, const Transform3D<T> &t)
1282 {
1283    Rotation3D r3d(r);
1284    return Transform3D<T>(r3d * t.Rotation(), r3d * t.Translation().Vect());
1285 }
1286 
1287 
1288 //---I/O functions
1289 // TODO - I/O should be put in the manipulator form
1290 
1291 /**
1292    print the 12 components of the Transform3D
1293  */
1294 template <class T>
1295 std::ostream &operator<<(std::ostream &os, const Transform3D<T> &t)
1296 {
1297    // TODO - this will need changing for machine-readable issues
1298    //        and even the human readable form needs formatting improvements
1299 
1300    T m[12];
1301    t.GetComponents(m, m + 12);
1302    os << "\n" << m[0] << "  " << m[1] << "  " << m[2] << "  " << m[3];
1303    os << "\n" << m[4] << "  " << m[5] << "  " << m[6] << "  " << m[7];
1304    os << "\n" << m[8] << "  " << m[9] << "  " << m[10] << "  " << m[11] << "\n";
1305    return os;
1306 }
1307 
1308 } // end namespace Impl
1309 
1310 // typedefs for double and float versions
1311 typedef Impl::Transform3D<double> Transform3D;
1312 typedef Impl::Transform3D<float>  Transform3DF;
1313 
1314 } // end namespace Math
1315 
1316 } // end namespace ROOT
1317 
1318 
1319 #endif /* ROOT_Math_GenVector_Transform3D */