CaloGeometryTools/interface/Transform3DPJ.h

0001 // @(#)root/mathcore:$Name: V02-02-29 $:$Id: Transform3DPJ.h,v 1.1 2007/05/27 20:12:30 pjanot Exp $
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 Transform3DPJ
0012 //
0013 // Created by: Lorenzo Moneta  October 21 2005
0014 //
0015 //
0016 #ifndef ROOT_Math_GenVector_Transform3DPJ
0017 #define ROOT_Math_GenVector_Transform3DPJ 1
0018
0019 #include "Math/GenVector/Cartesian3D.h"
0020 #include "Math/GenVector/DisplacementVector3D.h"
0021 #include "Math/GenVector/PositionVector3D.h"
0022 #include "Math/GenVector/LorentzVector.h"
0023 #include "Math/GenVector/Rotation3D.h"
0024 #include "Math/GenVector/AxisAnglefwd.h"
0025 #include "Math/GenVector/EulerAnglesfwd.h"
0026 #include "Math/GenVector/Quaternionfwd.h"
0027 #include "Math/GenVector/RotationXfwd.h"
0028 #include "Math/GenVector/RotationYfwd.h"
0029 #include "Math/GenVector/RotationZfwd.h"
0030 #include "Math/GenVector/Plane3D.h"
0031
0032 #include <iostream>
0033
0034 //#include "Math/Vector3Dfwd.h"
0035
0036 namespace ROOT {
0037
0038   namespace Math {
0039
0040     using ROOT::Math::Plane3D;
0041
0042     /**
0043      Basic 3D Transformation class describing  a rotation and then a translation
0044      The internal data are a rotation data and a 3D vector data and they can be represented
0045      like a 3x4 matrix
0046      The class has a template parameter the coordinate system tag of the reference system
0047      to which the transformatioon will be applied. For example for transforming from
0048      global to local coordinate systems, the transfrom3D has to be instantiated with the
0049      coordinate of the traget system
0050
0051      @ingroup GenVector
0052
0053   */
0054
0055     class Transform3DPJ {
0056     public:
0057       typedef DisplacementVector3D<Cartesian3D<double>, DefaultCoordinateSystemTag> Vector;
0058       typedef PositionVector3D<Cartesian3D<double>, DefaultCoordinateSystemTag> Point;
0059
0060       enum ETransform3DMatrixIndex {
0061         kXX = 0,
0062         kXY = 1,
0063         kXZ = 2,
0064         kDX = 3,
0065         kYX = 4,
0066         kYY = 5,
0067         kYZ = 6,
0068         kDY = 7,
0069         kZX = 8,
0070         kZY = 9,
0071         kZZ = 10,
0072         kDZ = 11
0073       };
0074
0075       /**
0076     Default constructor (identy rotation) + zero translation
0077     */
0078       Transform3DPJ() { SetIdentity(); }
0079
0080       /**
0081        Construct given a pair of pointers or iterators defining the
0082        beginning and end of an array of 12 Scalars
0083     */
0084       template <class IT>
0085       Transform3DPJ(IT begin, IT end) {
0086         SetComponents(begin, end);
0087       }
0088
0089       /**
0090        Construct from a rotation and then a translation described by a Vector
0091     */
0092       Transform3DPJ(const Rotation3D &r, const Vector &v) { AssignFrom(r, v); }
0093       /**
0094        Construct from a translation and then a rotation (inverse assignment)
0095     */
0096       Transform3DPJ(const Vector &v, const Rotation3D &r) {
0097         // is equivalent from having first the rotation and then the translation vector rotated
0098         AssignFrom(r, r(v));
0099       }
0100
0101       /**
0102        Construct from a 3D Rotation only with zero translation
0103     */
0104       explicit Transform3DPJ(const Rotation3D &r) { AssignFrom(r); }
0105       // convenience methods for the other rotations (cannot use templates for conflict with LA)
0106       explicit Transform3DPJ(const AxisAngle &r) { AssignFrom(Rotation3D(r)); }
0107       explicit Transform3DPJ(const EulerAngles &r) { AssignFrom(Rotation3D(r)); }
0108       explicit Transform3DPJ(const Quaternion &r) { AssignFrom(Rotation3D(r)); }
0109       // TO DO: implement direct methods for axial rotations without going through Rotation3D
0110       explicit Transform3DPJ(const RotationX &r) { AssignFrom(Rotation3D(r)); }
0111       explicit Transform3DPJ(const RotationY &r) { AssignFrom(Rotation3D(r)); }
0112       explicit Transform3DPJ(const RotationZ &r) { AssignFrom(Rotation3D(r)); }
0113
0114       /**
0115        Construct from a translation only, represented by any DisplacementVector3D
0116        and with an identity rotation
0117     */
0118       template <class CoordSystem, class Tag>
0119       explicit Transform3DPJ(const DisplacementVector3D<CoordSystem, Tag> &v) {
0120         AssignFrom(Vector(v.X(), v.Y(), v.Z()));
0121       }
0122       /**
0123        Construct from a translation only, represented by a Cartesian 3D Vector,
0124        and with an identity rotation
0125     */
0126       explicit Transform3DPJ(const Vector &v) { AssignFrom(v); }
0127
0128       //#if !defined(__MAKECINT__) && !defined(G__DICTIONARY)  // this is ambigous with double * , double *
0129       /**
0130        Construct from a rotation (any rotation object)  and then a translation
0131        (represented by any DisplacementVector)
0132        The requirements on the rotation and vector objects are that they can be transformed in a
0133        Rotation3D class and in a Vector
0134     */
0135       // to do : change to displacement vector3D
0136       template <class ARotation, class CoordSystem, class Tag>
0137       Transform3DPJ(const ARotation &r, const DisplacementVector3D<CoordSystem, Tag> &v) {
0138         AssignFrom(Rotation3D(r), Vector(v.X(), v.Y(), v.Z()));
0139       }
0140       /**
0141        Construct from a translation (using any type of DisplacementVector )
0142        and then a rotation (any rotation object).
0143        Requirement on the rotation and vector objects are that they can be transformed in a
0144        Rotation3D class and in a Vector
0145     */
0146       template <class ARotation, class CoordSystem, class Tag>
0147       Transform3DPJ(const DisplacementVector3D<CoordSystem, Tag> &v, const ARotation &r) {
0148         // is equivalent from having first the rotation and then the translation vector rotated
0149         Rotation3D r3d(r);
0150         AssignFrom(r3d, r3d(Vector(v.X(), v.Y(), v.Z())));
0151       }
0152
0153       //#endif
0154
0155       /**
0156        Construct transformation from one coordinate system defined by three
0157        points (origin + two axis) to
0158        a new coordinate system defined by other three points (origin + axis)
0159        @param fr0  point defining origin of original reference system
0160        @param fr1  point defining first axis of original reference system
0161        @param fr2  point defining second axis of original reference system
0162        @param to0  point defining origin of transformed reference system
0163        @param to1  point defining first axis transformed reference system
0164        @param to2  point defining second axis transformed reference system
0165
0166      */
0167       Transform3DPJ(
0168           const Point &fr0, const Point &fr1, const Point &fr2, const Point &to0, const Point &to1, const Point &to2);
0169
0170       // use compiler generated copy ctor, copy assignmet and dtor
0171
0172       /**
0173        Construct from a linear algebra matrix of size at least 3x4,
0174        which must support operator()(i,j) to obtain elements (0,0) thru (2,3).
0175        The 3x3 sub-block is assumed to be the rotation part and the translations vector
0176        are described by the 4-th column
0177     */
0178       template <class ForeignMatrix>
0179       explicit Transform3DPJ(const ForeignMatrix &m) {
0180         SetComponents(m);
0181       }
0182
0183       /**
0184        Raw constructor from 12 Scalar components
0185     */
0186       Transform3DPJ(double xx,
0187                     double xy,
0188                     double xz,
0189                     double dx,
0190                     double yx,
0191                     double yy,
0192                     double yz,
0193                     double dy,
0194                     double zx,
0195                     double zy,
0196                     double zz,
0197                     double dz) {
0198         SetComponents(xx, xy, xz, dx, yx, yy, yz, dy, zx, zy, zz, dz);
0199       }
0200
0201       /**
0202        Construct from a linear algebra matrix of size at least 3x4,
0203        which must support operator()(i,j) to obtain elements (0,0) thru (2,3).
0204        The 3x3 sub-block is assumed to be the rotation part and the translations vector
0205        are described by the 4-th column
0206     */
0207       template <class ForeignMatrix>
0208       Transform3DPJ &operator=(const ForeignMatrix &m) {
0209         SetComponents(m);
0210         return *this;
0211       }
0212
0213       // ======== Components ==============
0214
0215       /**
0216        Set the 12 matrix components given an iterator to the start of
0217        the desired data, and another to the end (12 past start).
0218     */
0219       template <class IT>
0220       void SetComponents(IT begin, IT end) {
0221         for (int i = 0; i < 12; ++i) {
0222           fM[i] = *begin;
0223           ++begin;
0224         }
0225         assert(end == begin);
0226       }
0227
0228       /**
0229        Get the 12 matrix components into data specified by an iterator begin
0230        and another to the end of the desired data (12 past start).
0231     */
0232       template <class IT>
0233       void GetComponents(IT begin, IT end) const {
0234         for (int i = 0; i < 12; ++i) {
0235           *begin = fM[i];
0236           ++begin;
0237         }
0238         assert(end == begin);
0239       }
0240
0241       /**
0242        Get the 12 matrix components into data specified by an iterator begin
0243     */
0244       template <class IT>
0245       void GetComponents(IT begin) const {
0246         std::copy(fM, fM + 12, begin);
0247       }
0248
0249       /**
0250        Set components from a linear algebra matrix of size at least 3x4,
0251        which must support operator()(i,j) to obtain elements (0,0) thru (2,3).
0252        The 3x3 sub-block is assumed to be the rotation part and the translations vector
0253        are described by the 4-th column
0254     */
0255       template <class ForeignMatrix>
0256       void SetTransformMatrix(const ForeignMatrix &m) {
0257         fM[kXX] = m(0, 0);
0258         fM[kXY] = m(0, 1);
0259         fM[kXZ] = m(0, 2);
0260         fM[kDX] = m(0, 3);
0261         fM[kYX] = m(1, 0);
0262         fM[kYY] = m(1, 1);
0263         fM[kYZ] = m(1, 2);
0264         fM[kDY] = m(1, 3);
0265         fM[kZX] = m(2, 0);
0266         fM[kZY] = m(2, 1);
0267         fM[kZZ] = m(2, 2);
0268         fM[kDZ] = m(2, 3);
0269       }
0270
0271       /**
0272        Get components into a linear algebra matrix of size at least 3x4,
0273        which must support operator()(i,j) for write access to elements
0274        (0,0) thru (2,3).
0275     */
0276       template <class ForeignMatrix>
0277       void GetTransformMatrix(ForeignMatrix &m) const {
0278         m(0, 0) = fM[kXX];
0279         m(0, 1) = fM[kXY];
0280         m(0, 2) = fM[kXZ];
0281         m(0, 3) = fM[kDX];
0282         m(1, 0) = fM[kYX];
0283         m(1, 1) = fM[kYY];
0284         m(1, 2) = fM[kYZ];
0285         m(1, 3) = fM[kDY];
0286         m(2, 0) = fM[kZX];
0287         m(2, 1) = fM[kZY];
0288         m(2, 2) = fM[kZZ];
0289         m(2, 3) = fM[kDZ];
0290       }
0291
0292       /**
0293        Set the components from 12 scalars
0294     */
0295       void SetComponents(double xx,
0296                          double xy,
0297                          double xz,
0298                          double dx,
0299                          double yx,
0300                          double yy,
0301                          double yz,
0302                          double dy,
0303                          double zx,
0304                          double zy,
0305                          double zz,
0306                          double dz) {
0307         fM[kXX] = xx;
0308         fM[kXY] = xy;
0309         fM[kXZ] = xz;
0310         fM[kDX] = dx;
0311         fM[kYX] = yx;
0312         fM[kYY] = yy;
0313         fM[kYZ] = yz;
0314         fM[kDY] = dy;
0315         fM[kZX] = zx;
0316         fM[kZY] = zy;
0317         fM[kZZ] = zz;
0318         fM[kDZ] = dz;
0319       }
0320
0321       /**
0322        Get the nine components into 12 scalars
0323     */
0324       void GetComponents(double &xx,
0325                          double &xy,
0326                          double &xz,
0327                          double &dx,
0328                          double &yx,
0329                          double &yy,
0330                          double &yz,
0331                          double &dy,
0332                          double &zx,
0333                          double &zy,
0334                          double &zz,
0335                          double &dz) const {
0336         xx = fM[kXX];
0337         xy = fM[kXY];
0338         xz = fM[kXZ];
0339         dx = fM[kDX];
0340         yx = fM[kYX];
0341         yy = fM[kYY];
0342         yz = fM[kYZ];
0343         dy = fM[kDY];
0344         zx = fM[kZX];
0345         zy = fM[kZY];
0346         zz = fM[kZZ];
0347         dz = fM[kDZ];
0348       }
0349
0350       /**
0351        Get the rotation and translation vector representing the 3D transformation
0352     */
0353       void GetDecomposition(Rotation3D &r, Vector &v) const;
0354
0355       // operations on points and vectors
0356
0357       /**
0358        Transformation operation for Position Vector in Cartesian coordinate
0359     */
0360       Point operator()(const Point &p) const;
0361
0362       /**
0363        Transformation operation for Displacement Vectors in Cartesian coordinate
0364        For the Displacement Vectors only the rotation applies - no translations
0365     */
0366       Vector operator()(const Vector &v) const;
0367
0368       /**
0369        Transformation operation for Position Vector in any coordinate system
0370     */
0371       template <class CoordSystem>
0372       PositionVector3D<CoordSystem> operator()(const PositionVector3D<CoordSystem> &p) const {
0373         Point xyzNew = operator()(Point(p));
0374         return PositionVector3D<CoordSystem>(xyzNew);
0375       }
0376
0377       /**
0378        Transformation operation for Displacement Vector in any coordinate system
0379     */
0380       template <class CoordSystem>
0381       DisplacementVector3D<CoordSystem> operator()(const DisplacementVector3D<CoordSystem> &v) const {
0382         Vector xyzNew = operator()(Vector(v));
0383         return DisplacementVector3D<CoordSystem>(xyzNew);
0384       }
0385
0386       /**
0387        Transformation operation for points between different coordinate system tags
0388     */
0389       template <class CoordSystem, class Tag1, class Tag2>
0390       void Transform(const PositionVector3D<CoordSystem, Tag1> &p1, PositionVector3D<CoordSystem, Tag2> &p2) const {
0391         Point xyzNew = operator()(Point(p1.X(), p1.Y(), p1.Z()));
0392         p2.SetXYZ(xyzNew.X(), xyzNew.Y(), xyzNew.Z());
0393       }
0394
0395       /**
0396        Transformation operation for Displacement Vector of different coordinate systems
0397     */
0398       template <class CoordSystem, class Tag1, class Tag2>
0399       void Transform(const DisplacementVector3D<CoordSystem, Tag1> &v1,
0400                      DisplacementVector3D<CoordSystem, Tag2> &v2) const {
0401         Vector xyzNew = operator()(Vector(v1.X(), v1.Y(), v1.Z()));
0402         v2.SetXYZ(xyzNew.X(), xyzNew.Y(), xyzNew.Z());
0403       }
0404
0405       /**
0406        Transformation operation for a Lorentz Vector in any  coordinate system
0407     */
0408       template <class CoordSystem>
0409       LorentzVector<CoordSystem> operator()(const LorentzVector<CoordSystem> &q) const {
0410         Vector xyzNew = operator()(Vector(q.Vect()));
0411         return LorentzVector<CoordSystem>(xyzNew.X(), xyzNew.Y(), xyzNew.Z(), q.E());
0412       }
0413
0414       /**
0415        Transformation on a 3D plane
0416     */
0417       Plane3D operator()(const Plane3D &plane) const;
0418
0419       // skip transformation for arbitrary vectors - not really defined if point or displacement vectors
0420
0421       // same but with operator *
0422       /**
0423        Transformation operation for Vectors. Apply same rules as operator()
0424        depending on type of vector.
0425        Will work only for DisplacementVector3D, PositionVector3D and LorentzVector
0426     */
0427       template <class AVector>
0428       AVector operator*(const AVector &v) const {
0429         return operator()(v);
0430       }
0431
0432       /**
0433        multiply (combine) with another transformation in place
0434      */
0435       Transform3DPJ &operator*=(const Transform3DPJ &t);
0436
0437       /**
0438        multiply (combine) two transformations
0439      */
0440       Transform3DPJ operator*(const Transform3DPJ &t) const {
0441         Transform3DPJ tmp(*this);
0442         tmp *= t;
0443         return tmp;
0444       }
0445
0446       /**
0447     Invert the transformation in place
0448     */
0449       void Invert();
0450
0451       /**
0452        Return the inverse of the transformation.
0453     */
0454       Transform3DPJ Inverse() const {
0455         Transform3DPJ t(*this);
0456         t.Invert();
0457         return t;
0458       }
0459
0460       /**
0461        Equality/inequality operators
0462     */
0463       bool operator==(const Transform3DPJ &rhs) const {
0464         if (fM[0] != rhs.fM[0])
0465           return false;
0466         if (fM[1] != rhs.fM[1])
0467           return false;
0468         if (fM[2] != rhs.fM[2])
0469           return false;
0470         if (fM[3] != rhs.fM[3])
0471           return false;
0472         if (fM[4] != rhs.fM[4])
0473           return false;
0474         if (fM[5] != rhs.fM[5])
0475           return false;
0476         if (fM[6] != rhs.fM[6])
0477           return false;
0478         if (fM[7] != rhs.fM[7])
0479           return false;
0480         if (fM[8] != rhs.fM[8])
0481           return false;
0482         if (fM[9] != rhs.fM[9])
0483           return false;
0484         if (fM[10] != rhs.fM[10])
0485           return false;
0486         if (fM[11] != rhs.fM[11])
0487           return false;
0488         return true;
0489       }
0490
0491       bool operator!=(const Transform3DPJ &rhs) const { return !operator==(rhs); }
0492
0493     protected:
0494       /**
0495        make transformation from first a rotation then a translation
0496      */
0497       void AssignFrom(const Rotation3D &r, const Vector &v);
0498
0499       /**
0500        make transformation from only rotations (zero translation)
0501      */
0502       void AssignFrom(const Rotation3D &r);
0503
0504       /**
0505        make transformation from only translation (identity rotations)
0506      */
0507       void AssignFrom(const Vector &v);
0508
0509       /**
0510        Set identity transformation (identity rotation , zero translation)
0511      */
0512       void SetIdentity();
0513
0514     private:
0515       double fM[12];
0516     };
0517
0518     // global functions
0519
0520     // TODO - I/O should be put in the manipulator form
0521
0522     std::ostream &operator<<(std::ostream &os, const Transform3DPJ &t);
0523
0524     // need a function Transform = Translation * Rotation ???
0525
0526   }  // end namespace Math
0527
0528 }  // end namespace ROOT
0529
0530 #endif /* MATHCORE_BASIC3DTRANSFORMATION */