VecGeom/base/Transformation2DMP.h

0001 /// \file Transformation2DMP.h This provides the precision templated Transformation2D
0002
0003 #ifndef VECGEOM_BASE_TRANSFORMATION2DMP_H_
0004 #define VECGEOM_BASE_TRANSFORMATION2DMP_H_
0005
0006 #include "VecGeom/base/Config.h"
0007 #include "VecGeom/base/Cuda.h"
0008 #include "VecGeom/base/Global.h"
0009
0010 #include <VecGeom/base/Transformation3D.h>
0011
0012 #include "VecGeom/base/Vector3D.h"
0013 #include "VecGeom/base/Vector2D.h"
0014 #include "VecGeom/backend/scalar/Backend.h"
0015 #ifdef VECGEOM_ENABLE_CUDA
0016 #include "VecGeom/backend/cuda/Interface.h"
0017 #endif
0018
0019 #include <algorithm>
0020 #include <cmath>
0021 #include <cstring>
0022 #include <iostream>
0023 #include <vector>
0024
0025 #include <initializer_list>
0026 // #ifdef VECGEOM_ROOT
0027 // class TGeoMatrix;
0028 // #endif
0029
0030 namespace vecgeom {
0031
0032 VECGEOM_DEVICE_FORWARD_DECLARE(class Transformation2DMP;);
0033
0034 inline namespace VECGEOM_IMPL_NAMESPACE {
0035
0036 #ifndef VECCORE_CUDA
0037 }
0038 namespace cuda {
0039 class Transformation2DMP;
0040 }
0041 inline namespace VECGEOM_IMPL_NAMESPACE {
0042 // class vecgeom::cuda::Transformation2D;
0043 #endif
0044
0045 template <typename Real_s>
0046 class Transformation2DMP {
0047
0048   // On the precision templating: Real_s is the precision used for the stored variables,
0049   // Real_i the one for the input variables.
0050
0051 private:
0052   // TODO: it might be better to directly store this in terms of Vector3D<Real_t> !!
0053   // and would allow for higher level abstraction
0054   Real_s tx_{0.}, ty_{0.}, tz_{0.};
0055   Real_s rxx_{1.}, ryx_{0.};
0056   Real_s rxy_{0.}, ryy_{1.};
0057   bool fIdentity;
0058   bool fHasRotation;
0059   bool fHasTranslation;
0060
0061 public:
0062   VECCORE_ATT_HOST_DEVICE
0063   constexpr Transformation2DMP() : fIdentity(true), fHasRotation(false), fHasTranslation(false){};
0064
0065   /**
0066    * Constructor from other transformation in other precision.
0067    * @tparam Real_i precision of other transformation
0068    * @param other other transformation
0069    */
0070   template <typename Real_i>
0071   Transformation2DMP(const Transformation2DMP<Real_i> &other)
0072   {
0073     // Assuming getTranslation returns a pointer or array of 3 elements
0074     auto translation = other.Translation();
0075     tx_              = static_cast<Real_s>(translation[0]);
0076     ty_              = static_cast<Real_s>(translation[1]);
0077     tz_              = static_cast<Real_s>(translation[2]);
0078
0079     // Assuming getRotation returns a pointer or array of 9 elements
0080     auto rotation = other.Rotation();
0081     rxx_          = static_cast<Real_s>(rotation[0]);
0082     ryx_          = static_cast<Real_s>(rotation[1]);
0083     rxy_          = static_cast<Real_s>(rotation[2]);
0084     ryy_          = static_cast<Real_s>(rotation[3]);
0085
0086     fIdentity       = other.IsIdentity();
0087     fHasRotation    = other.HasRotation();
0088     fHasTranslation = other.HasTranslation();
0089   }
0090
0091   /**
0092    * Constructor from 3D transformation in other precision.
0093    * @tparam Real_i precision of other transformation
0094    * @param other other transformation
0095    */
0096   template <typename Real_i>
0097   Transformation2DMP(const Transformation3DMP<Real_i> &other)
0098   {
0099     // Assuming getTranslation returns a pointer or array of 3 elements
0100     auto translation = other.Translation();
0101     tx_              = static_cast<Real_s>(translation[0]);
0102     ty_              = static_cast<Real_s>(translation[1]);
0103     tz_              = static_cast<Real_s>(translation[2]);
0104
0105     // Assuming getRotation returns a pointer or array of 9 elements
0106     auto rotation = other.Rotation();
0107
0108     // 2D transformation should only be created for 2D rotations, so rzx, rxz, rzy, ryz should be 0 and ryy should be
0109     // +-1
0110     VECGEOM_ASSERT(Abs(rotation[2]) < vecgeom::kTolerance && Abs(rotation[5]) < vecgeom::kTolerance &&
0111                    Abs(rotation[6]) < vecgeom::kTolerance && Abs(rotation[7]) < vecgeom::kTolerance &&
0112                    Abs(rotation[8]) - 1. < vecgeom::kTolerance);
0113
0114     rxx_ = static_cast<Real_s>(rotation[0]);
0115     ryx_ = static_cast<Real_s>(rotation[1]);
0116     rxy_ = static_cast<Real_s>(rotation[3]);
0117     ryy_ = static_cast<Real_s>(rotation[4]);
0118
0119     fIdentity       = other.IsIdentity();
0120     fHasRotation    = other.HasRotation();
0121     fHasTranslation = other.HasTranslation();
0122   }
0123
0124   /**
0125    * Constructor for translation only.
0126    * @param tx Translation in x-coordinate.
0127    * @param ty Translation in y-coordinate.
0128    * @param tz Translation in z-coordinate.
0129    */
0130   template <typename Real_i>
0131   VECCORE_ATT_HOST_DEVICE Transformation2DMP(const Real_i tx, const Real_i ty, const Real_i tz)
0132       : tx_{static_cast<Real_s>(tx)}, ty_{static_cast<Real_s>(ty)}, tz_{static_cast<Real_s>(tz)},
0133         fIdentity(tx == 0 && ty == 0 && tz == 0), fHasRotation(false), fHasTranslation(tx != 0 || ty != 0 || tz != 0)
0134   {
0135   }
0136
0137   /**
0138    * Constructor to manually set each entry. Used when converting from different
0139    * geometry.
0140    */
0141   template <typename Real_i>
0142   VECCORE_ATT_HOST_DEVICE Transformation2DMP(const Real_i tx, const Real_i ty, const Real_i tz, const Real_i rxx,
0143                                              const Real_i ryx, const Real_i rxy, const Real_i ryy)
0144       : fIdentity(false), fHasRotation(true), fHasTranslation(true)
0145   {
0146     SetTranslation(tx, ty, tz);
0147     SetRotation(rxx, ryx, rxy, ryy);
0148     SetProperties();
0149   }
0150   /**
0151    * Constructor using the rotation from a different transformation
0152    */
0153   template <typename Real_i>
0154   VECCORE_ATT_HOST_DEVICE Transformation2DMP(const Real_i tx, const Real_i ty, const Real_i tz,
0155                                              Transformation2DMP<Real_i> const &trot)
0156   {
0157     SetTranslation(tx, ty, tz);
0158     SetRotation(trot.Rotation()[0], trot.Rotation()[1], trot.Rotation()[2], trot.Rotation()[3]);
0159     SetProperties();
0160   }
0161   /**
0162    * Constructor to manually set each entry. Used when converting from different
0163    * geometry, supporting scaling.
0164    */
0165   template <typename Real_i>
0166   VECCORE_ATT_HOST_DEVICE Transformation2DMP(const Real_i tx, const Real_i ty, const Real_i tz, const Real_i rxx,
0167                                              const Real_i ryx, const Real_i rxy, const Real_i ryy, Real_i sx, Real_i sy)
0168   {
0169     SetTranslation(tx, ty, tz);
0170     SetRotation(rxx, ryx, rxy, ryy);
0171     ApplyScale(sx, sy);
0172     SetProperties();
0173   }
0174
0175   /**
0176    * Constructor copying the translation and rotation from memory
0177    * geometry.
0178    */
0179   template <typename Real_i>
0180   VECCORE_ATT_HOST_DEVICE Transformation2DMP(const Real_i *trans, const Real_i *rot, bool has_trans, bool has_rot)
0181   {
0182     this->Set(trans, rot, has_trans, has_rot);
0183   }
0184
0185   VECCORE_ATT_HOST_DEVICE
0186   VECGEOM_FORCE_INLINE
0187   bool operator==(Transformation2DMP const &rhs) const;
0188
0189   VECCORE_ATT_HOST_DEVICE
0190   VECGEOM_FORCE_INLINE
0191   Transformation2DMP operator*(Transformation2DMP const &tf) const;
0192
0193   VECCORE_ATT_HOST_DEVICE
0194   VECGEOM_FORCE_INLINE
0195   Transformation2DMP const &operator*=(Transformation2DMP const &rhs);
0196
0197   VECCORE_ATT_HOST_DEVICE
0198   ~Transformation2DMP() {}
0199
0200   template <typename Real_i>
0201   VECCORE_ATT_HOST_DEVICE void ApplyScale(const Real_i sx, const Real_i sy)
0202   {
0203     rxx_ *= static_cast<Real_s>(sx);
0204     ryx_ *= static_cast<Real_s>(sy);
0205     rxy_ *= static_cast<Real_s>(sx);
0206     ryy_ *= static_cast<Real_s>(sy);
0207   }
0208
0209   template <typename Real_i>
0210   VECCORE_ATT_HOST_DEVICE bool ApproxEqual(Transformation2DMP<Real_i> const &rhs,
0211                                            Real_i tolerance = kToleranceDist<Real_i>) const
0212   {
0213     auto const data1 = &tx_;
0214     auto const data2 = &rhs.tx_;
0215     for (int i = 0; i < 7; ++i)
0216       if (vecCore::math::Abs(data1[i] - data2[i]) > tolerance) return false;
0217     return true;
0218   }
0219
0220   VECCORE_ATT_HOST_DEVICE
0221   void Clear()
0222   {
0223     tx_             = 0.;
0224     ty_             = 0.;
0225     tz_             = 0.;
0226     rxx_            = 1.;
0227     ryx_            = 0.;
0228     rxy_            = 0.;
0229     ryy_            = 1.;
0230     fIdentity       = true;
0231     fHasRotation    = false;
0232     fHasTranslation = false;
0233   }
0234   VECCORE_ATT_HOST_DEVICE
0235   VECGEOM_FORCE_INLINE
0236   Real_s Determinant() const { return (rxx_ * ryy_ - rxy_ * ryx_); }
0237
0238   VECCORE_ATT_HOST_DEVICE
0239   VECGEOM_FORCE_INLINE
0240   bool IsReflected() const { return (Determinant() < Real_s(0)); }
0241
0242   int MemorySize() const { return sizeof(*this); }
0243
0244   VECCORE_ATT_HOST_DEVICE
0245   void FixZeroes()
0246   {
0247     if (std::abs(tx_) < vecgeom::kTolerance) tx_ = 0.;
0248     if (std::abs(ty_) < vecgeom::kTolerance) ty_ = 0.;
0249     if (std::abs(tz_) < vecgeom::kTolerance) tz_ = 0.;
0250     if (std::abs(rxx_) < vecgeom::kTolerance) rxx_ = 0.;
0251     if (std::abs(ryx_) < vecgeom::kTolerance) ryx_ = 0.;
0252     if (std::abs(rxy_) < vecgeom::kTolerance) rxy_ = 0.;
0253     if (std::abs(ryy_) < vecgeom::kTolerance) ryy_ = 0.;
0254   }
0255
0256   VECCORE_ATT_HOST_DEVICE
0257   VECGEOM_FORCE_INLINE
0258   Vector3D<Real_s> Translation() const { return Vector3D<Real_s>(tx_, ty_, tz_); }
0259
0260   /**
0261    * No safety against faulty indexing.
0262    * @param index Index of translation entry in the range [0-2].
0263    */
0264   VECCORE_ATT_HOST_DEVICE
0265   VECGEOM_FORCE_INLINE
0266   Real_s Translation(const int index) const { return *(&tx_ + index); }
0267
0268   VECCORE_ATT_HOST_DEVICE
0269   VECGEOM_FORCE_INLINE
0270   Real_s const *Rotation() const { return &rxx_; }
0271
0272   /**
0273    * No safety against faulty indexing.
0274    * \param index Index of rotation entry in the range [0-3].
0275    */
0276   VECCORE_ATT_HOST_DEVICE
0277   VECGEOM_FORCE_INLINE
0278   Real_s Rotation(const int index) const { return *(&rxx_ + index); }
0279
0280   VECCORE_ATT_HOST_DEVICE
0281   VECGEOM_FORCE_INLINE
0282   bool IsIdentity() const { return fIdentity; }
0283
0284   VECCORE_ATT_HOST_DEVICE
0285   VECGEOM_FORCE_INLINE
0286   bool IsXYRotation() const { return (fHasRotation); }
0287
0288   VECCORE_ATT_HOST_DEVICE
0289   VECGEOM_FORCE_INLINE
0290   bool HasRotation() const { return fHasRotation; }
0291
0292   VECCORE_ATT_HOST_DEVICE
0293   VECGEOM_FORCE_INLINE
0294   bool HasTranslation() const { return fHasTranslation; }
0295
0296   void Print(std::ostream &s) const
0297   {
0298     s << "Transformation2DMP {{" << tx_ << "," << ty_ << "," << tz_ << "}";
0299     s << "{" << rxx_ << "," << ryx_ << "," << rxy_ << "," << ryy_ << "}}\n";
0300   }
0301
0302   VECCORE_ATT_HOST_DEVICE
0303   void Print() const
0304   {
0305     printf("Transformation2D {{%.2f, %.2f, %.2f}, ", tx_, ty_, tz_);
0306     printf("{%.2f, %.2f, %.2f, %.2f}}", rxx_, ryx_, rxy_, ryy_);
0307   }
0308
0309   // Mutators
0310
0311   template <typename Real_i>
0312   VECCORE_ATT_HOST_DEVICE void SetTranslation(const Real_i tx, const Real_i ty, const Real_i tz)
0313   {
0314     tx_ = static_cast<Real_s>(tx);
0315     ty_ = static_cast<Real_s>(ty);
0316     tz_ = static_cast<Real_s>(tz);
0317   }
0318
0319   template <typename Real_i>
0320   VECCORE_ATT_HOST_DEVICE void SetTranslation(Vector3D<Real_i> const &vec)
0321   {
0322     SetTranslation(vec[0], vec[1], vec[2]);
0323   }
0324
0325   // FIXME tolerance still compile option
0326   VECCORE_ATT_HOST_DEVICE
0327   void SetProperties()
0328   {
0329     fHasTranslation = (fabs(tx_) > kTolerance || fabs(ty_) > kTolerance || fabs(tz_) > kTolerance) ? true : false;
0330     fHasRotation    = (fabs(rxx_ - 1.) > kTolerance) || (fabs(ryx_) > kTolerance) || (fabs(rxy_) > kTolerance) ||
0331                            (fabs(ryy_ - 1.) > kTolerance)
0332                           ? true
0333                           : false;
0334     fIdentity       = !fHasTranslation && !fHasRotation;
0335   }
0336
0337   template <typename Real_i>
0338   VECCORE_ATT_HOST_DEVICE void SetRotation(const Real_i xx, const Real_i yx, const Real_i xy, const Real_i yy)
0339   {
0340     rxx_ = static_cast<Real_s>(xx);
0341     ryx_ = static_cast<Real_s>(yx);
0342     rxy_ = static_cast<Real_s>(xy);
0343     ryy_ = static_cast<Real_s>(yy);
0344   }
0345
0346   // /**
0347   //  * Set rotation given an arbitrary axis and an angle.
0348   //  * \param aaxis Rotation axis. No need to be a unit vector.
0349   //  * \param ddelta Rotation angle in radians.
0350   //  */
0351   // VECCORE_ATT_HOST_DEVICE
0352   // void Set(Vector3D<double> const &aaxis, double ddelta);
0353
0354   /**
0355    * Set transformation and rotation.
0356    * \param trans Pointer to at least 3 values.
0357    * \param rot Pointer to at least 4 values.
0358    */
0359   template <typename Real_i>
0360   VECCORE_ATT_HOST_DEVICE VECGEOM_FORCE_INLINE void Set(const Real_i *trans, const Real_i *rot, bool has_trans,
0361                                                         bool has_rot)
0362   {
0363     if (has_trans) {
0364       tx_ = static_cast<Real_s>(trans[0]);
0365       ty_ = static_cast<Real_s>(trans[1]);
0366       tz_ = static_cast<Real_s>(trans[2]);
0367     }
0368
0369     if (has_rot) {
0370       rxx_ = static_cast<Real_s>(rot[0]);
0371       ryx_ = static_cast<Real_s>(rot[1]);
0372       rxy_ = static_cast<Real_s>(rot[2]);
0373       ryy_ = static_cast<Real_s>(rot[3]);
0374     }
0375
0376     fHasTranslation = has_trans;
0377     fHasRotation    = has_rot;
0378     fIdentity       = !fHasTranslation && !fHasRotation;
0379   }
0380
0381 private:
0382   // Templated rotation and translation methods which inline and compile to
0383   // optimized versions.
0384
0385   VECGEOM_FORCE_INLINE VECCORE_ATT_HOST_DEVICE void DoRotation(Vector3D<Real_s> const &master,
0386                                                                Vector3D<Real_s> &local) const;
0387
0388 private:
0389   VECGEOM_FORCE_INLINE VECCORE_ATT_HOST_DEVICE void DoTranslation(Vector3D<Real_s> const &master,
0390                                                                   Vector3D<Real_s> &local) const;
0391
0392   template <bool vectortransform>
0393   VECGEOM_FORCE_INLINE VECCORE_ATT_HOST_DEVICE void InverseTransformKernel(Vector3D<Real_s> const &local,
0394                                                                            Vector3D<Real_s> &master) const;
0395
0396 public:
0397   // Transformation interface
0398
0399   VECGEOM_FORCE_INLINE VECCORE_ATT_HOST_DEVICE void Transform(Vector3D<Real_s> const &master,
0400                                                               Vector3D<Real_s> &local) const;
0401
0402   VECGEOM_FORCE_INLINE VECCORE_ATT_HOST_DEVICE Vector3D<Real_s> Transform(Vector3D<Real_s> const &master) const;
0403
0404   VECGEOM_FORCE_INLINE VECCORE_ATT_HOST_DEVICE void TransformDirection(Vector3D<Real_s> const &master,
0405                                                                        Vector3D<Real_s> &local) const;
0406
0407   VECGEOM_FORCE_INLINE VECCORE_ATT_HOST_DEVICE Vector3D<Real_s> TransformDirection(
0408       Vector3D<Real_s> const &master) const;
0409
0410   /** The inverse transformation ( aka LocalToMaster ) of an object transform like a point
0411    *  this does not need to currently template on placement since such a transformation is much less used
0412    */
0413   VECGEOM_FORCE_INLINE VECCORE_ATT_HOST_DEVICE void InverseTransform(Vector3D<Real_s> const &local,
0414                                                                      Vector3D<Real_s> &master) const;
0415
0416   VECGEOM_FORCE_INLINE VECCORE_ATT_HOST_DEVICE Vector3D<Real_s> InverseTransform(Vector3D<Real_s> const &local) const;
0417
0418   /** The inverse transformation of an object transforming like a vector */
0419   VECGEOM_FORCE_INLINE VECCORE_ATT_HOST_DEVICE void InverseTransformDirection(Vector3D<Real_s> const &master,
0420                                                                               Vector3D<Real_s> &local) const;
0421
0422   VECGEOM_FORCE_INLINE VECCORE_ATT_HOST_DEVICE Vector3D<Real_s> InverseTransformDirection(
0423       Vector3D<Real_s> const &master) const;
0424
0425   /** compose transformations - multiply transformations */
0426   VECCORE_ATT_HOST_DEVICE
0427   VECGEOM_FORCE_INLINE
0428   void MultiplyFromRight(Transformation2DMP const &rhs);
0429
0430   /**
0431    * @brief Returns the inverse of this matrix. Taken from G4AffineTransformation
0432    */
0433   VECCORE_ATT_HOST_DEVICE
0434   Transformation2DMP<Real_s> Inverse() const
0435   {
0436     Real_s ttx = -tx_, tty = -ty_, ttz = -tz_;
0437     return Transformation2DMP<Real_s>(ttx * rxx_ + tty * rxy_, ttx * ryx_ + tty * ryy_, ttz, rxx_, rxy_, ryx_, ryy_);
0438   }
0439
0440   // Utility and CUDA
0441
0442 #ifdef VECGEOM_CUDA_INTERFACE
0443   size_t DeviceSizeOf() const { return DevicePtr<cuda::Transformation2DMP>::SizeOf(); }
0444   DevicePtr<cuda::Transformation2DMP> CopyToGpu() const;
0445   DevicePtr<cuda::Transformation2DMP> CopyToGpu(DevicePtr<cuda::Transformation2DMP> const gpu_ptr) const;
0446   static void CopyManyToGpu(const std::vector<Transformation2DMP const *> &trafos,
0447                             const std::vector<DevicePtr<cuda::Transformation2DMP>> &gpu_ptrs);
0448 #endif
0449
0450 public:
0451   static const Transformation2DMP kIdentity;
0452
0453 }; // End class Transformation2D
0454
0455 template <typename Real_t>
0456 VECCORE_ATT_HOST_DEVICE bool Transformation2DMP<Real_t>::operator==(Transformation2DMP<Real_t> const &rhs) const
0457 {
0458   return equal(&tx_, &tx_ + 7, &rhs.tx_);
0459 }
0460
0461 /**
0462  * Rotates a vector to this transformation's frame of reference.
0463  * \param master Vector in original frame of reference.
0464  * \param local Output vector rotated to the new frame of reference.
0465  */
0466 template <typename Real_s>
0467 VECGEOM_FORCE_INLINE VECCORE_ATT_HOST_DEVICE void Transformation2DMP<Real_s>::DoRotation(Vector3D<Real_s> const &master,
0468                                                                                          Vector3D<Real_s> &local) const
0469 {
0470   local[0] = master[0] * rxx_;
0471   local[1] = master[0] * ryx_;
0472   local[2] = master[2]; // this must be set since the non-existing rzz_ = 1.
0473   local[0] += master[1] * rxy_;
0474   local[1] += master[1] * ryy_;
0475 }
0476
0477 template <typename Real_s>
0478 VECGEOM_FORCE_INLINE VECCORE_ATT_HOST_DEVICE void Transformation2DMP<Real_s>::DoTranslation(
0479     Vector3D<Real_s> const &master, Vector3D<Real_s> &local) const
0480 {
0481
0482   local[0] = master[0] - tx_;
0483   local[1] = master[1] - ty_;
0484   local[2] = master[2] - tz_;
0485 }
0486
0487 /**
0488  * Transform a point to the local reference frame.
0489  * \param master Point to be transformed.
0490  * \param local Output destination. Should never be the same as the input
0491  *              vector!
0492  */
0493 template <typename Real_s>
0494 VECGEOM_FORCE_INLINE VECCORE_ATT_HOST_DEVICE void Transformation2DMP<Real_s>::Transform(Vector3D<Real_s> const &master,
0495                                                                                         Vector3D<Real_s> &local) const
0496 {
0497   Vector3D<Real_s> tmp;
0498   DoTranslation(master, tmp);
0499   DoRotation(tmp, local);
0500 }
0501
0502 /**
0503  * Since transformation cannot be done in place, allows the transformed vector
0504  * to be constructed by Transform directly.
0505  * \param master Point to be transformed.
0506  * \return Newly constructed Vector3D with the transformed coordinates.
0507  */
0508 template <typename Real_s>
0509 VECGEOM_FORCE_INLINE VECCORE_ATT_HOST_DEVICE Vector3D<Real_s> Transformation2DMP<Real_s>::Transform(
0510     Vector3D<Real_s> const &master) const
0511 {
0512
0513   Vector3D<Real_s> local;
0514   Transform(master, local);
0515   return local;
0516 }
0517
0518 template <typename Real_s>
0519 template <bool transform_direction>
0520 VECGEOM_FORCE_INLINE VECCORE_ATT_HOST_DEVICE void Transformation2DMP<Real_s>::InverseTransformKernel(
0521     Vector3D<Real_s> const &local, Vector3D<Real_s> &master) const
0522 {
0523
0524   // we are just doing the full stuff here ( LocalToMaster is less critical
0525   // than other way round )
0526
0527   if (transform_direction) {
0528     master[0] = local[0] * rxx_;
0529     master[0] += local[1] * ryx_;
0530     master[1] = local[0] * rxy_;
0531     master[1] += local[1] * ryy_;
0532     master[2] = local[2]; // assuming rzz_ = 1
0533
0534   } else {
0535     master[0] = tx_;
0536     master[0] += local[0] * rxx_;
0537     master[0] += local[1] * ryx_;
0538     master[1] = ty_;
0539     master[1] += local[0] * rxy_;
0540     master[1] += local[1] * ryy_;
0541     master[2] = tz_;
0542     master[2] += local[2];
0543   }
0544 }
0545
0546 template <typename Real_s>
0547 VECGEOM_FORCE_INLINE VECCORE_ATT_HOST_DEVICE void Transformation2DMP<Real_s>::InverseTransform(
0548     Vector3D<Real_s> const &local, Vector3D<Real_s> &master) const
0549 {
0550   InverseTransformKernel<false>(local, master);
0551 }
0552
0553 template <typename Real_s>
0554 VECGEOM_FORCE_INLINE VECCORE_ATT_HOST_DEVICE Vector3D<Real_s> Transformation2DMP<Real_s>::InverseTransform(
0555     Vector3D<Real_s> const &local) const
0556 {
0557   Vector3D<Real_s> tmp;
0558   InverseTransform(local, tmp);
0559   return tmp;
0560 }
0561
0562 template <typename Real_s>
0563 VECGEOM_FORCE_INLINE VECCORE_ATT_HOST_DEVICE void Transformation2DMP<Real_s>::InverseTransformDirection(
0564     Vector3D<Real_s> const &local, Vector3D<Real_s> &master) const
0565 {
0566   InverseTransformKernel<true>(local, master);
0567 }
0568
0569 template <typename Real_s>
0570 VECGEOM_FORCE_INLINE VECCORE_ATT_HOST_DEVICE Vector3D<Real_s> Transformation2DMP<Real_s>::InverseTransformDirection(
0571     Vector3D<Real_s> const &local) const
0572 {
0573   Vector3D<Real_s> tmp;
0574   InverseTransformDirection(local, tmp);
0575   return tmp;
0576 }
0577
0578 template <typename Real_s>
0579 VECCORE_ATT_HOST_DEVICE VECGEOM_FORCE_INLINE Transformation2DMP<Real_s> Transformation2DMP<Real_s>::operator*(
0580     Transformation2DMP<Real_s> const &rhs) const
0581 {
0582   if (rhs.fIdentity) return Transformation2DMP<Real_s>(*this);
0583   return Transformation2DMP<Real_s>(tx_ * rhs.rxx_ + ty_ * rhs.ryx_ + tz_, // tx
0584                                     tx_ * rhs.rxy_ + ty_ * rhs.ryy_ + tz_, // ty
0585                                     tz_ + rhs.tz_,                         // tz
0586                                     rxx_ * rhs.rxx_ + rxy_ * rhs.ryx_,     // rxx
0587                                     ryx_ * rhs.rxx_ + ryy_ * rhs.ryx_,     // ryx
0588                                     rxx_ * rhs.rxy_ + rxy_ * rhs.ryy_,     // rxy
0589                                     ryx_ * rhs.rxy_ + ryy_ * rhs.ryy_);    // ryy
0590 }
0591
0592 template <typename Real_s>
0593 VECCORE_ATT_HOST_DEVICE VECGEOM_FORCE_INLINE void Transformation2DMP<Real_s>::MultiplyFromRight(
0594     Transformation2DMP<Real_s> const &rhs)
0595 {
0596   // TODO: this code should directly operator on Vector3D and Matrix3D
0597
0598   if (rhs.fIdentity) return;
0599   fIdentity = false;
0600
0601   if (rhs.HasTranslation()) {
0602     fHasTranslation = true;
0603     // ideal for fused multiply add
0604     tx_ += rxx_ * rhs.tx_;
0605     tx_ += ryx_ * rhs.ty_;
0606
0607     ty_ += rxy_ * rhs.tx_;
0608     ty_ += ryy_ * rhs.ty_;
0609
0610     tz_ += rhs.tz_;
0611   }
0612
0613   if (rhs.HasRotation()) {
0614     fHasRotation   = true;
0615     Precision tmpx = rxx_;
0616     Precision tmpy = ryx_;
0617
0618     // first row of matrix
0619     rxx_ = tmpx * rhs.rxx_;
0620     ryx_ = tmpx * rhs.ryx_;
0621     rxx_ += tmpy * rhs.rxy_;
0622     ryx_ += tmpy * rhs.ryy_;
0623
0624     tmpx = rxy_;
0625     tmpy = ryy_;
0626
0627     // second row of matrix
0628     rxy_ = tmpx * rhs.rxx_;
0629     ryy_ = tmpx * rhs.ryx_;
0630     rxy_ += tmpy * rhs.rxy_;
0631     ryy_ += tmpy * rhs.ryy_;
0632   }
0633 }
0634
0635 /**
0636  * Only transforms by rotation, ignoring the translation part. This is useful
0637  * when transforming directions.
0638  * \param master Point to be transformed.
0639  * \param local Output destination of transformation.
0640  */
0641 template <typename Real_s>
0642 VECGEOM_FORCE_INLINE VECCORE_ATT_HOST_DEVICE void Transformation2DMP<Real_s>::TransformDirection(
0643     Vector3D<Real_s> const &master, Vector3D<Real_s> &local) const
0644 {
0645   DoRotation(master, local);
0646 }
0647
0648 /**
0649  * Since transformation cannot be done in place, allows the transformed vector
0650  * to be constructed by TransformDirection directly.
0651  * \param master Point to be transformed.
0652  * \return Newly constructed Vector3D with the transformed coordinates.
0653  */
0654 template <typename Real_s>
0655 VECGEOM_FORCE_INLINE VECCORE_ATT_HOST_DEVICE Vector3D<Real_s> Transformation2DMP<Real_s>::TransformDirection(
0656     Vector3D<Real_s> const &master) const
0657 {
0658
0659   Vector3D<Real_s> local;
0660   TransformDirection(master, local);
0661   return local;
0662 }
0663
0664 template <typename Real_s>
0665 std::ostream &operator<<(std::ostream &os, Transformation2DMP<Real_s> const &trans)
0666 {
0667   os << "TransformationMP {" << trans.Translation() << ", " << "(" << trans.Rotation(0) << ", " << trans.Rotation(1)
0668      << ", " << trans.Rotation(2) << ", " << trans.Rotation(3) << ")}" << "; identity(" << trans.IsIdentity()
0669      << "); rotation(" << trans.HasRotation() << ")";
0670   return os;
0671 }
0672 }
0673 } // namespace vecgeom
0674
0675 #endif // VECGEOM_BASE_TRANSFORMATION2DMP_H_