Back to home page

EIC code displayed by LXR
 
 

    

File indexing completed on 2026-04-17 08:35:29 

0001 /// \file Transformation3DMP.h This provides the precision templated Transformation3D, which is required for mixed precision
0002 
0003 #ifndef VECGEOM_BASE_TRANSFORMATION3DMP_H_
0004 #define VECGEOM_BASE_TRANSFORMATION3DMP_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/backend/scalar/Backend.h"
0014 #ifdef VECGEOM_ENABLE_CUDA
0015 #include "VecGeom/backend/cuda/Interface.h"
0016 #endif
0017 
0018 #include <algorithm>
0019 #include <cmath>
0020 #include <cstring>
0021 #include <iostream>
0022 #include <vector>
0023 
0024 #include <initializer_list>
0025 // #ifdef VECGEOM_ROOT
0026 // class TGeoMatrix;
0027 // #endif
0028 
0029 namespace vecgeom {
0030 
0031 VECGEOM_DEVICE_FORWARD_DECLARE(class Transformation3DMP;);
0032 
0033 inline namespace VECGEOM_IMPL_NAMESPACE {
0034 
0035 #ifndef VECCORE_CUDA
0036 }
0037 namespace cuda {
0038 class Transformation3DMP;
0039 }
0040 inline namespace VECGEOM_IMPL_NAMESPACE {
0041 // class vecgeom::cuda::Transformation3D;
0042 #endif
0043 
0044 template <typename Real_s>
0045 class Transformation3DMP {
0046 
0047   // On the precision templating: Real_s is the precision used for the stored variables,
0048   // Real_i the one for the input variables.
0049 
0050 private:
0051   // TODO: it might be better to directly store this in terms of Vector3D<Real_t> !!
0052   // and would allow for higher level abstraction
0053   Real_s tx_{0.}, ty_{0.}, tz_{0.};
0054   Real_s rxx_{1.}, ryx_{0.}, rzx_{0.};
0055   Real_s rxy_{0.}, ryy_{1.}, rzy_{0.};
0056   Real_s rxz_{0.}, ryz_{0.}, rzz_{1.};
0057   bool fIdentity;
0058   bool fHasRotation;
0059   bool fHasTranslation;
0060 
0061 public:
0062   VECCORE_ATT_HOST_DEVICE
0063   constexpr Transformation3DMP() : 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   Transformation3DMP(const Transformation3DMP<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     rzx_          = static_cast<Real_s>(rotation[2]);
0084     rxy_          = static_cast<Real_s>(rotation[3]);
0085     ryy_          = static_cast<Real_s>(rotation[4]);
0086     rzy_          = static_cast<Real_s>(rotation[5]);
0087     rxz_          = static_cast<Real_s>(rotation[6]);
0088     ryz_          = static_cast<Real_s>(rotation[7]);
0089     rzz_          = static_cast<Real_s>(rotation[8]);
0090 
0091     fIdentity       = other.IsIdentity();
0092     fHasRotation    = other.HasRotation();
0093     fHasTranslation = other.HasTranslation();
0094   }
0095 
0096   /**
0097    * Constructor for translation only.
0098    * @param tx Translation in x-coordinate.
0099    * @param ty Translation in y-coordinate.
0100    * @param tz Translation in z-coordinate.
0101    */
0102   template <typename Real_i>
0103   VECCORE_ATT_HOST_DEVICE Transformation3DMP(const Real_i tx, const Real_i ty, const Real_i tz)
0104       : tx_{static_cast<Real_s>(tx)}, ty_{static_cast<Real_s>(ty)}, tz_{static_cast<Real_s>(tz)},
0105         fIdentity(tx == 0 && ty == 0 && tz == 0), fHasRotation(false), fHasTranslation(tx != 0 || ty != 0 || tz != 0)
0106   {
0107   }
0108 
0109   /**
0110    * @brief Rotation followed by translation.
0111    * @param tx Translation in x-coordinate.
0112    * @param ty Translation in y-coordinate.
0113    * @param tz Translation in z-coordinate.
0114    * @param phi Rotation angle about z-axis.
0115    * @param theta Rotation angle about new y-axis.
0116    * @param psi Rotation angle about new z-axis.
0117    */
0118   template <typename Real_i>
0119   VECCORE_ATT_HOST_DEVICE Transformation3DMP(const Real_i tx, const Real_i ty, const Real_i tz, const Real_i phi,
0120                                              const Real_i theta, const Real_i psi)
0121       : fIdentity(false), fHasRotation(true), fHasTranslation(true)
0122   {
0123     SetTranslation(tx, ty, tz);
0124     SetRotation(phi, theta, psi);
0125     SetProperties();
0126   }
0127 
0128   // enabling initilizer lists
0129   template <typename Real_i>
0130   Transformation3DMP(std::initializer_list<Real_i> values)
0131   {
0132     auto it      = values.begin();
0133     tx_          = *it++;
0134     ty_          = *it++;
0135     tz_          = *it++;
0136     Real_i phi   = *it++;
0137     Real_i theta = *it++;
0138     Real_i psi   = *it;
0139     SetTranslation(tx_, ty_, tz_);
0140     SetRotation(phi, theta, psi);
0141     SetProperties();
0142   }
0143 
0144   /**
0145    * @brief RScale followed by rotation followed by translation.
0146    * @param tx Translation in x-coordinate.
0147    * @param ty Translation in y-coordinate.
0148    * @param tz Translation in z-coordinate.
0149    * @param phi Rotation angle about z-axis.
0150    * @param theta Rotation angle about new y-axis.
0151    * @param psi Rotation angle about new z-axis.
0152    * @param sx X scaling factor
0153    * @param sy Y scaling factor
0154    * @param sz Z scaling factor
0155    */
0156   template <typename Real_i>
0157   VECCORE_ATT_HOST_DEVICE Transformation3DMP(const Real_i tx, const Real_i ty, const Real_i tz, const Real_i phi,
0158                                              const Real_i theta, const Real_i psi, Real_i sx, Real_i sy, Real_i sz)
0159       : fIdentity(false), fHasRotation(true), fHasTranslation(true)
0160   {
0161     SetTranslation(tx, ty, tz);
0162     SetRotation(phi, theta, psi);
0163     ApplyScale(sx, sy, sz);
0164     SetProperties();
0165   }
0166 
0167   /**
0168    * Constructor to manually set each entry. Used when converting from different
0169    * geometry.
0170    */
0171   template <typename Real_i>
0172   VECCORE_ATT_HOST_DEVICE Transformation3DMP(const Real_i tx, const Real_i ty, const Real_i tz, const Real_i rxx,
0173                                              const Real_i ryx, const Real_i rzx, const Real_i rxy, const Real_i ryy,
0174                                              const Real_i rzy, const Real_i rxz, const Real_i ryz, const Real_i rzz)
0175       : fIdentity(false), fHasRotation(true), fHasTranslation(true)
0176   {
0177     SetTranslation(tx, ty, tz);
0178     SetRotation(rxx, ryx, rzx, rxy, ryy, rzy, rxz, ryz, rzz);
0179     SetProperties();
0180   }
0181   /**
0182    * Constructor using the rotation from a different transformation
0183    */
0184   template <typename Real_i>
0185   VECCORE_ATT_HOST_DEVICE Transformation3DMP(const Real_i tx, const Real_i ty, const Real_i tz,
0186                                              Transformation3DMP<Real_i> const &trot)
0187   {
0188     SetTranslation(tx, ty, tz);
0189     SetRotation(trot.Rotation()[0], trot.Rotation()[1], trot.Rotation()[2], trot.Rotation()[3], trot.Rotation()[4],
0190                 trot.Rotation()[5], trot.Rotation()[6], trot.Rotation()[7], trot.Rotation()[8]);
0191     SetProperties();
0192   }
0193   /**
0194    * Constructor to manually set each entry. Used when converting from different
0195    * geometry, supporting scaling.
0196    */
0197   template <typename Real_i>
0198   VECCORE_ATT_HOST_DEVICE Transformation3DMP(const Real_i tx, const Real_i ty, const Real_i tz, const Real_i rxx,
0199                                              const Real_i ryx, const Real_i rzx, const Real_i rxy, const Real_i ryy,
0200                                              const Real_i rzy, const Real_i rxz, const Real_i ryz, const Real_i rzz,
0201                                              Real_i sx, Real_i sy, Real_i sz)
0202   {
0203     SetTranslation(tx, ty, tz);
0204     SetRotation(rxx, ryx, rzx, rxy, ryy, rzy, rxz, ryz, rzz);
0205     ApplyScale(sx, sy, sz);
0206     SetProperties();
0207   }
0208 
0209   /**
0210    * Constructor copying the translation and rotation from memory
0211    * geometry.
0212    */
0213   template <typename Real_i>
0214   VECCORE_ATT_HOST_DEVICE Transformation3DMP(const Real_i *trans, const Real_i *rot, bool has_trans, bool has_rot)
0215   {
0216     this->Set(trans, rot, has_trans, has_rot);
0217   }
0218 
0219   /**
0220    * @brief Construct from a Transformation3D holding the elements in vecgeom::Precition
0221    */
0222   VECCORE_ATT_HOST_DEVICE Transformation3DMP(Transformation3D const &other)
0223   {
0224     auto rot = other.Rotation();
0225     SetTranslation(other.Translation(0), other.Translation(1), other.Translation(2));
0226     SetRotation(rot[0], rot[1], rot[2], rot[3], rot[4], rot[5], rot[6], rot[7], rot[8]);
0227     SetProperties();
0228   }
0229 
0230   // /**
0231   //  * Constructor for a rotation based on a given direction
0232   //  * @param axis direction of the new z axis
0233   //  * @param inverse if true the origial axis will be rotated into (0,0,u)
0234   //                   if false a vector (0,0,u) will be rotated into the original axis
0235   //  */
0236   // VECCORE_ATT_HOST_DEVICE
0237   // Transformation3D(const Vector3D<Precision> &axis, bool inverse = true);
0238 
0239   VECCORE_ATT_HOST_DEVICE
0240   VECGEOM_FORCE_INLINE
0241   bool operator==(Transformation3DMP const &rhs) const;
0242 
0243   VECCORE_ATT_HOST_DEVICE
0244   VECGEOM_FORCE_INLINE
0245   Transformation3DMP operator*(Transformation3DMP const &tf) const;
0246 
0247   VECCORE_ATT_HOST_DEVICE
0248   VECGEOM_FORCE_INLINE
0249   Transformation3DMP const &operator*=(Transformation3DMP const &rhs);
0250 
0251   // operator for multiplication with Transformation3D
0252   VECCORE_ATT_HOST_DEVICE
0253   VECGEOM_FORCE_INLINE
0254   Transformation3DMP operator*(Transformation3D const &rhs) const;
0255 
0256   VECCORE_ATT_HOST_DEVICE
0257   VECGEOM_FORCE_INLINE
0258   Transformation3DMP const &operator*=(Transformation3D const &rhs);
0259 
0260   VECCORE_ATT_HOST_DEVICE
0261   ~Transformation3DMP() {}
0262 
0263   template <typename Real_i>
0264   VECCORE_ATT_HOST_DEVICE void ApplyScale(const Real_i sx, const Real_i sy, const Real_i sz)
0265   {
0266     rxx_ *= static_cast<Real_s>(sx);
0267     ryx_ *= static_cast<Real_s>(sy);
0268     rzx_ *= static_cast<Real_s>(sz);
0269     rxy_ *= static_cast<Real_s>(sx);
0270     ryy_ *= static_cast<Real_s>(sy);
0271     rzy_ *= static_cast<Real_s>(sz);
0272     rxz_ *= static_cast<Real_s>(sx);
0273     ryz_ *= static_cast<Real_s>(sy);
0274     rzz_ *= static_cast<Real_s>(sz);
0275   }
0276 
0277   // CHECK THIS IF WE NEED BOTH WITH MP and WITHOUT MP INPUT OR ONLY THIS ONE FIXME
0278   template <typename Real_i>
0279   VECCORE_ATT_HOST_DEVICE bool ApproxEqual(Transformation3DMP<Real_i> const &rhs,
0280                                            Real_i tolerance = kToleranceDist<Real_i>) const
0281   {
0282     auto const data1 = &tx_;
0283     auto const data2 = &rhs.tx_;
0284     for (int i = 0; i < 12; ++i)
0285       if (vecCore::math::Abs(data1[i] - data2[i]) > tolerance) return false;
0286     return true;
0287   }
0288 
0289   VECCORE_ATT_HOST_DEVICE
0290   void Clear()
0291   {
0292     tx_             = 0.;
0293     ty_             = 0.;
0294     tz_             = 0.;
0295     rxx_            = 1.;
0296     ryx_            = 0.;
0297     rzx_            = 0.;
0298     rxy_            = 0.;
0299     ryy_            = 1.;
0300     rzy_            = 0.;
0301     rxz_            = 0.;
0302     ryz_            = 0.;
0303     rzz_            = 1.;
0304     fIdentity       = true;
0305     fHasRotation    = false;
0306     fHasTranslation = false;
0307   }
0308 
0309   int MemorySize() const { return sizeof(*this); }
0310 
0311   VECCORE_ATT_HOST_DEVICE
0312   void FixZeroes()
0313   {
0314     if (std::abs(tx_) < vecgeom::kTolerance) tx_ = 0.;
0315     if (std::abs(ty_) < vecgeom::kTolerance) ty_ = 0.;
0316     if (std::abs(tz_) < vecgeom::kTolerance) tz_ = 0.;
0317     if (std::abs(rxx_) < vecgeom::kTolerance) rxx_ = 0.;
0318     if (std::abs(ryx_) < vecgeom::kTolerance) ryx_ = 0.;
0319     if (std::abs(rzx_) < vecgeom::kTolerance) rzx_ = 0.;
0320     if (std::abs(rxy_) < vecgeom::kTolerance) rxy_ = 0.;
0321     if (std::abs(ryy_) < vecgeom::kTolerance) ryy_ = 0.;
0322     if (std::abs(rzy_) < vecgeom::kTolerance) rzy_ = 0.;
0323     if (std::abs(rxz_) < vecgeom::kTolerance) rxz_ = 0.;
0324     if (std::abs(ryz_) < vecgeom::kTolerance) ryz_ = 0.;
0325     if (std::abs(rzz_) < vecgeom::kTolerance) rzz_ = 0.;
0326   }
0327 
0328   VECCORE_ATT_HOST_DEVICE
0329   VECGEOM_FORCE_INLINE
0330   Vector3D<Real_s> Translation() const { return Vector3D<Real_s>(tx_, ty_, tz_); }
0331 
0332   // VECCORE_ATT_HOST_DEVICE
0333   // VECGEOM_FORCE_INLINE
0334   // Vector3D<Real_s> Scale() const
0335   // {
0336   //   Real_s sx = std::sqrt(rxx_ * rxx_ + rxy_ * rxy_ + rxz_ * rxz_);
0337   //   Real_s sy = std::sqrt(ryx_ * ryx_ + ryy_ * ryy_ + ryz_ * ryz_);
0338   //   Real_s sz = std::sqrt(rzx_ * rzx_ + rzy_ * rzy_ + rzz_ * rzz_);
0339 
0340   //   if (Determinant() < 0) sz = -sz;
0341 
0342   //   return Vector3D<Real_s>(sx, sy, sz);
0343   // }
0344 
0345   /**
0346    * No safety against faulty indexing.
0347    * @param index Index of translation entry in the range [0-2].
0348    */
0349   VECCORE_ATT_HOST_DEVICE
0350   VECGEOM_FORCE_INLINE
0351   Real_s Translation(const int index) const { return *(&tx_ + index); }
0352 
0353   VECCORE_ATT_HOST_DEVICE
0354   VECGEOM_FORCE_INLINE
0355   Real_s const *Rotation() const { return &rxx_; }
0356 
0357   /**
0358    * No safety against faulty indexing.
0359    * \param index Index of rotation entry in the range [0-8].
0360    */
0361   VECCORE_ATT_HOST_DEVICE
0362   VECGEOM_FORCE_INLINE
0363   Real_s Rotation(const int index) const { return *(&rxx_ + index); }
0364 
0365   VECCORE_ATT_HOST_DEVICE
0366   VECGEOM_FORCE_INLINE
0367   bool IsIdentity() const { return fIdentity; }
0368 
0369   VECCORE_ATT_HOST_DEVICE
0370   VECGEOM_FORCE_INLINE
0371   bool IsXYRotation() const
0372   {
0373     return (fHasRotation && (std::abs(rzx_) < vecgeom::kTolerance) && (std::abs(rzy_) < vecgeom::kTolerance) &&
0374             (std::abs(rxz_) < vecgeom::kTolerance) && (std::abs(ryz_) < vecgeom::kTolerance) &&
0375             (std::abs(rzz_ - Real_s(1.)) < vecgeom::kTolerance));
0376   }
0377 
0378   VECCORE_ATT_HOST_DEVICE
0379   VECGEOM_FORCE_INLINE
0380   bool IsXYFlipRotation() const
0381   {
0382     return (fHasRotation && (std::abs(rzx_) < vecgeom::kTolerance) && (std::abs(rzy_) < vecgeom::kTolerance) &&
0383             (std::abs(rxz_) < vecgeom::kTolerance) && (std::abs(ryz_) < vecgeom::kTolerance) &&
0384             (std::abs(rzz_ + Real_s(1.)) < vecgeom::kTolerance));
0385   }
0386 
0387   VECCORE_ATT_HOST_DEVICE
0388   VECGEOM_FORCE_INLINE
0389   bool IsReflected() const { return Determinant() < 0; }
0390 
0391   VECCORE_ATT_HOST_DEVICE
0392   VECGEOM_FORCE_INLINE
0393   bool HasRotation() const { return fHasRotation; }
0394 
0395   VECCORE_ATT_HOST_DEVICE
0396   VECGEOM_FORCE_INLINE
0397   bool HasTranslation() const { return fHasTranslation; }
0398 
0399   // VECCORE_ATT_HOST_DEVICE
0400   // void Print() const;
0401 
0402   // // print to a stream
0403   // void Print(std::ostream &) const;
0404 
0405   // VECCORE_ATT_HOST_DEVICE
0406   // void PrintG4() const;
0407 
0408   void Print(std::ostream &s) const
0409   {
0410     s << "Transformation3DMP {{" << tx_ << "," << ty_ << "," << tz_ << "}";
0411     s << "{" << rxx_ << "," << ryx_ << "," << rzx_ << "," << rxy_ << "," << ryy_ << "," << rzy_ << "," << rxz_ << ","
0412       << ryz_ << "," << rzz_ << "}}\n";
0413   }
0414 
0415   VECCORE_ATT_HOST_DEVICE
0416   void Print() const
0417   {
0418     printf("Transformation3D {{%.2f, %.2f, %.2f}, ", tx_, ty_, tz_);
0419     printf("{%.2f, %.2f, %.2f, %.2f, %.2f, %.2f, %.2f, %.2f, %.2f}}", rxx_, ryx_, rzx_, rxy_, ryy_, rzy_, rxz_, ryz_,
0420            rzz_);
0421   }
0422 
0423   VECCORE_ATT_HOST_DEVICE
0424   void PrintG4() const
0425   {
0426     using vecCore::math::Abs;
0427     using vecCore::math::Max;
0428     constexpr double deviationTolerance = 1.0e-05;
0429     printf("  TransformationMP: \n");
0430 
0431     bool UnitTr            = !fHasRotation;
0432     double diagDeviation   = Max(Abs(rxx_ - 1.0), Abs(ryy_ - 1.0), Abs(rzz_ - 1.0));
0433     double offdDeviationUL = Max(Abs(rxy_), Abs(rxz_), Abs(ryx_));
0434     double offdDeviationDR = Max(Abs(ryz_), Abs(rzx_), Abs(rzy_));
0435     double offdDeviation   = Max(offdDeviationUL, offdDeviationDR);
0436 
0437     if (UnitTr || Max(diagDeviation, offdDeviation) < deviationTolerance) {
0438       printf("    UNIT  Rotation \n");
0439     } else {
0440       printf("rx/x,y,z: %.6g %.6g %.6g\nry/x,y,z: %.6g %.6g %.6g\nrz/x,y,z: %.6g %.6g %.6g\n", rxx_, rxy_, rxz_, ryx_,
0441              ryy_, ryz_, rzx_, rzy_, rzz_);
0442     }
0443 
0444     printf("tr/x,y,z: %.6g %.6g %.6g\n", tx_, ty_, tz_);
0445   }
0446 
0447   // Mutators
0448 
0449   template <typename Real_i>
0450   VECCORE_ATT_HOST_DEVICE void SetTranslation(const Real_i tx, const Real_i ty, const Real_i tz)
0451   {
0452     tx_ = static_cast<Real_s>(tx);
0453     ty_ = static_cast<Real_s>(ty);
0454     tz_ = static_cast<Real_s>(tz);
0455   }
0456 
0457   template <typename Real_i>
0458   VECCORE_ATT_HOST_DEVICE void SetTranslation(Vector3D<Real_i> const &vec)
0459   {
0460     SetTranslation(vec[0], vec[1], vec[2]);
0461   }
0462 
0463   // FIXME tolerance still compile option
0464   VECCORE_ATT_HOST_DEVICE
0465   void SetProperties()
0466   {
0467     fHasTranslation = (fabs(tx_) > kTolerance || fabs(ty_) > kTolerance || fabs(tz_) > kTolerance) ? true : false;
0468     fHasRotation    = (fabs(rxx_ - 1.) > kTolerance) || (fabs(ryx_) > kTolerance) || (fabs(rzx_) > kTolerance) ||
0469                            (fabs(rxy_) > kTolerance) || (fabs(ryy_ - 1.) > kTolerance) || (fabs(rzy_) > kTolerance) ||
0470                            (fabs(rxz_) > kTolerance) || (fabs(ryz_) > kTolerance) || (fabs(rzz_ - 1.) > kTolerance)
0471                           ? true
0472                           : false;
0473     fIdentity       = !fHasTranslation && !fHasRotation;
0474   }
0475 
0476   template <typename Real_i>
0477   VECCORE_ATT_HOST_DEVICE void SetRotation(const Real_i phi, const Real_i theta, const Real_i psi)
0478   {
0479 
0480     const auto sinphi = sin(kDegToRad * phi);
0481     const auto cosphi = cos(kDegToRad * phi);
0482     const auto sinthe = sin(kDegToRad * theta);
0483     const auto costhe = cos(kDegToRad * theta);
0484     const auto sinpsi = sin(kDegToRad * psi);
0485     const auto cospsi = cos(kDegToRad * psi);
0486 
0487     rxx_ = static_cast<Real_s>(cospsi * cosphi - costhe * sinphi * sinpsi);
0488     ryx_ = static_cast<Real_s>(-sinpsi * cosphi - costhe * sinphi * cospsi);
0489     rzx_ = static_cast<Real_s>(sinthe * sinphi);
0490     rxy_ = static_cast<Real_s>(cospsi * sinphi + costhe * cosphi * sinpsi);
0491     ryy_ = static_cast<Real_s>(-sinpsi * sinphi + costhe * cosphi * cospsi);
0492     rzy_ = static_cast<Real_s>(-sinthe * cosphi);
0493     rxz_ = static_cast<Real_s>(sinpsi * sinthe);
0494     ryz_ = static_cast<Real_s>(cospsi * sinthe);
0495     rzz_ = static_cast<Real_s>(costhe);
0496   }
0497 
0498   template <typename Real_i>
0499   VECCORE_ATT_HOST_DEVICE void SetRotation(Vector3D<Real_i> const &vec)
0500   {
0501     SetRotation(vec[0], vec[1], vec[2]);
0502   }
0503 
0504   template <typename Real_i>
0505   VECCORE_ATT_HOST_DEVICE void SetRotation(const Real_i xx, const Real_i yx, const Real_i zx, const Real_i xy,
0506                                            const Real_i yy, const Real_i zy, const Real_i xz, const Real_i yz,
0507                                            const Real_i zz)
0508   {
0509     rxx_ = static_cast<Real_s>(xx);
0510     ryx_ = static_cast<Real_s>(yx);
0511     rzx_ = static_cast<Real_s>(zx);
0512     rxy_ = static_cast<Real_s>(xy);
0513     ryy_ = static_cast<Real_s>(yy);
0514     rzy_ = static_cast<Real_s>(zy);
0515     rxz_ = static_cast<Real_s>(xz);
0516     ryz_ = static_cast<Real_s>(yz);
0517     rzz_ = static_cast<Real_s>(zz);
0518   }
0519 
0520   // /**
0521   //  * Set rotation given an arbitrary axis and an angle.
0522   //  * \param aaxis Rotation axis. No need to be a unit vector.
0523   //  * \param ddelta Rotation angle in radians.
0524   //  */
0525   // VECCORE_ATT_HOST_DEVICE
0526   // void Set(Vector3D<double> const &aaxis, double ddelta);
0527 
0528   /**
0529    * Set transformation and rotation.
0530    * \param trans Pointer to at least 3 values.
0531    * \param rot Pointer to at least 9 values.
0532    */
0533   template <typename Real_i>
0534   VECCORE_ATT_HOST_DEVICE VECGEOM_FORCE_INLINE void Set(const Real_i *trans, const Real_i *rot, bool has_trans,
0535                                                         bool has_rot)
0536   {
0537     if (has_trans) {
0538       tx_ = static_cast<Real_s>(trans[0]);
0539       ty_ = static_cast<Real_s>(trans[1]);
0540       tz_ = static_cast<Real_s>(trans[2]);
0541     }
0542 
0543     if (has_rot) {
0544       rxx_ = static_cast<Real_s>(rot[0]);
0545       ryx_ = static_cast<Real_s>(rot[1]);
0546       rzx_ = static_cast<Real_s>(rot[2]);
0547       rxy_ = static_cast<Real_s>(rot[3]);
0548       ryy_ = static_cast<Real_s>(rot[4]);
0549       rzy_ = static_cast<Real_s>(rot[5]);
0550       rxz_ = static_cast<Real_s>(rot[6]);
0551       ryz_ = static_cast<Real_s>(rot[7]);
0552       rzz_ = static_cast<Real_s>(rot[8]);
0553     }
0554 
0555     fHasTranslation = has_trans;
0556     fHasRotation    = has_rot;
0557     fIdentity       = !fHasTranslation && !fHasRotation;
0558   }
0559 
0560 private:
0561   // Templated rotation and translation methods which inline and compile to
0562   // optimized versions.
0563 
0564   VECGEOM_FORCE_INLINE VECCORE_ATT_HOST_DEVICE void DoRotation(Vector3D<Real_s> const &master,
0565                                                                Vector3D<Real_s> &local) const;
0566 
0567 private:
0568   VECGEOM_FORCE_INLINE VECCORE_ATT_HOST_DEVICE void DoTranslation(Vector3D<Real_s> const &master,
0569                                                                   Vector3D<Real_s> &local) const;
0570 
0571   template <bool vectortransform>
0572   VECGEOM_FORCE_INLINE VECCORE_ATT_HOST_DEVICE void InverseTransformKernel(Vector3D<Real_s> const &local,
0573                                                                            Vector3D<Real_s> &master) const;
0574 
0575 public:
0576   // Transformation interface
0577 
0578   VECGEOM_FORCE_INLINE VECCORE_ATT_HOST_DEVICE void Transform(Vector3D<Real_s> const &master,
0579                                                               Vector3D<Real_s> &local) const;
0580 
0581   VECGEOM_FORCE_INLINE VECCORE_ATT_HOST_DEVICE Vector3D<Real_s> Transform(Vector3D<Real_s> const &master) const;
0582 
0583   VECGEOM_FORCE_INLINE VECCORE_ATT_HOST_DEVICE void TransformDirection(Vector3D<Real_s> const &master,
0584                                                                        Vector3D<Real_s> &local) const;
0585 
0586   VECGEOM_FORCE_INLINE VECCORE_ATT_HOST_DEVICE Vector3D<Real_s> TransformDirection(
0587       Vector3D<Real_s> const &master) const;
0588 
0589   /** The inverse transformation ( aka LocalToMaster ) of an object transform like a point
0590    *  this does not need to currently template on placement since such a transformation is much less used
0591    */
0592   VECGEOM_FORCE_INLINE VECCORE_ATT_HOST_DEVICE void InverseTransform(Vector3D<Real_s> const &local,
0593                                                                      Vector3D<Real_s> &master) const;
0594 
0595   VECGEOM_FORCE_INLINE VECCORE_ATT_HOST_DEVICE Vector3D<Real_s> InverseTransform(Vector3D<Real_s> const &local) const;
0596 
0597   /** The inverse transformation of an object transforming like a vector */
0598   VECGEOM_FORCE_INLINE VECCORE_ATT_HOST_DEVICE void InverseTransformDirection(Vector3D<Real_s> const &master,
0599                                                                               Vector3D<Real_s> &local) const;
0600 
0601   VECGEOM_FORCE_INLINE VECCORE_ATT_HOST_DEVICE Vector3D<Real_s> InverseTransformDirection(
0602       Vector3D<Real_s> const &master) const;
0603 
0604   /** compose transformations - multiply transformations */
0605   VECCORE_ATT_HOST_DEVICE
0606   VECGEOM_FORCE_INLINE
0607   void MultiplyFromRight(Transformation3DMP const &rhs);
0608 
0609   VECCORE_ATT_HOST_DEVICE
0610   VECGEOM_FORCE_INLINE
0611   void MultiplyFromRight(Transformation3D const &rhs);
0612 
0613   // FIXME
0614   // /** compose transformations - multiply transformations */
0615   // VECCORE_ATT_HOST_DEVICE
0616   // VECGEOM_FORCE_INLINE
0617   // void CopyFrom(Transformation3D const &rhs)
0618   // {
0619   //   // not sure this compiles under CUDA
0620   //   copy(&rhs, &rhs + 1, this);
0621   // }
0622 
0623   // VECCORE_ATT_HOST_DEVICE
0624   // Transformation3D &RotateX(double a);
0625 
0626   // VECCORE_ATT_HOST_DEVICE
0627   // Transformation3D &RotateY(double a);
0628 
0629   // VECCORE_ATT_HOST_DEVICE
0630   // Transformation3D &RotateZ(double a);
0631 
0632   /**
0633    * @brief Returns determinant for the rotation.
0634    */
0635   VECCORE_ATT_HOST_DEVICE
0636   double Determinant() const
0637   {
0638     // Computes determinant in double precision
0639     double xx_ = rxx_;
0640     double xy_ = rxy_;
0641     double xz_ = rxz_;
0642     double yx_ = ryx_;
0643     double yy_ = ryy_;
0644     double yz_ = ryz_;
0645     double zx_ = rzx_;
0646     double zy_ = rzy_;
0647     double zz_ = rzz_;
0648 
0649     double detxx = yy_ * zz_ - yz_ * zy_;
0650     double detxy = yx_ * zz_ - yz_ * zx_;
0651     double detxz = yx_ * zy_ - yy_ * zx_;
0652     double det   = xx_ * detxx - xy_ * detxy + xz_ * detxz;
0653     return det;
0654   }
0655 
0656   // /**
0657   //  * @brief Returns rotation axis as a unit vector
0658   //  */
0659   // VECCORE_ATT_HOST_DEVICE
0660   // Vector3D<double> Axis() const;
0661 
0662   /**
0663    * @brief Rectifies the round-off making the rotation true.
0664    * @details  From clhep/src/RotationC.cc
0665      Assuming the representation of this is close to a true Rotation,
0666      but may have drifted due to round-off error from many operations,
0667      this forms an "exact" orthonormal matrix for the rotation again.
0668 
0669      The first step is to average with the transposed inverse.  This
0670      will correct for small errors such as those occuring when decomposing
0671      a LorentzTransformation.  Then we take the bull by the horns and
0672      formally extract the axis and delta (assuming the Rotation were true)
0673      and re-setting the rotation according to those.
0674    */
0675   // VECCORE_ATT_HOST_DEVICE
0676   // Transformation3D &Rectify();
0677 
0678   // /**
0679   //  * @brief Stores the inverse of this matrix into inverse. Taken from G4AffineTransformation
0680   //  */
0681   // VECCORE_ATT_HOST_DEVICE
0682   // Transformation3D &Invert()
0683   // {
0684   //   double ttx = -tx_, tty = -ty_, ttz = -tz_;
0685   //   tx_ = ttx * rxx_ + tty * rxy_ + ttz * rxz_;
0686   //   ty_ = ttx * ryx_ + tty * ryy_ + ttz * ryz_;
0687   //   tz_ = ttx * rzx_ + tty * rzy_ + ttz * rzz_;
0688 
0689   //   auto tmp1 = ryx_;
0690   //   ryx_      = rxy_;
0691   //   rxy_      = tmp1;
0692   //   auto tmp2 = rzx_;
0693   //   rzx_      = rxz_;
0694   //   rxz_      = tmp2;
0695   //   auto tmp3 = rzy_;
0696   //   rzy_      = ryz_;
0697   //   ryz_      = tmp3;
0698   //   return *this;
0699   // }
0700 
0701   /**
0702    * @brief Returns the inverse of this matrix. Taken from G4AffineTransformation
0703    */
0704   VECCORE_ATT_HOST_DEVICE
0705   Transformation3DMP<Real_s> Inverse() const
0706   {
0707     Real_s ttx = -tx_, tty = -ty_, ttz = -tz_;
0708     return Transformation3DMP<Real_s>(ttx * rxx_ + tty * rxy_ + ttz * rxz_, ttx * ryx_ + tty * ryy_ + ttz * ryz_,
0709                                       ttx * rzx_ + tty * rzy_ + ttz * rzz_, rxx_, rxy_, rxz_, ryx_, ryy_, ryz_, rzx_,
0710                                       rzy_, rzz_);
0711   }
0712 
0713   // Utility and CUDA
0714 
0715 #ifdef VECGEOM_CUDA_INTERFACE
0716   size_t DeviceSizeOf() const { return DevicePtr<cuda::Transformation3DMP>::SizeOf(); }
0717   DevicePtr<cuda::Transformation3DMP> CopyToGpu() const;
0718   DevicePtr<cuda::Transformation3DMP> CopyToGpu(DevicePtr<cuda::Transformation3DMP> const gpu_ptr) const;
0719   static void CopyManyToGpu(const std::vector<Transformation3DMP const *> &trafos,
0720                             const std::vector<DevicePtr<cuda::Transformation3DMP>> &gpu_ptrs);
0721 #endif
0722 
0723   // #ifdef VECGEOM_ROOT
0724   //   // function to convert this transformation to a TGeo transformation
0725   //   // mainly used for the benchmark comparisons with ROOT
0726   //   static TGeoMatrix *ConvertToTGeoMatrix(Transformation3D const &);
0727   // #endif
0728 
0729 public:
0730   static const Transformation3DMP kIdentity;
0731 
0732 }; // End class Transformation3D
0733 
0734 template <typename Real_t>
0735 VECCORE_ATT_HOST_DEVICE bool Transformation3DMP<Real_t>::operator==(Transformation3DMP<Real_t> const &rhs) const
0736 {
0737   return equal(&tx_, &tx_ + 12, &rhs.tx_);
0738 }
0739 
0740 /**
0741  * Rotates a vector to this transformation's frame of reference.
0742  * \param master Vector in original frame of reference.
0743  * \param local Output vector rotated to the new frame of reference.
0744  */
0745 template <typename Real_s>
0746 VECGEOM_FORCE_INLINE VECCORE_ATT_HOST_DEVICE void Transformation3DMP<Real_s>::DoRotation(Vector3D<Real_s> const &master,
0747                                                                                          Vector3D<Real_s> &local) const
0748 {
0749   local[0] = master[0] * rxx_;
0750   local[1] = master[0] * ryx_;
0751   local[2] = master[0] * rzx_;
0752   local[0] += master[1] * rxy_;
0753   local[1] += master[1] * ryy_;
0754   local[2] += master[1] * rzy_;
0755   local[0] += master[2] * rxz_;
0756   local[1] += master[2] * ryz_;
0757   local[2] += master[2] * rzz_;
0758 }
0759 
0760 template <typename Real_s>
0761 VECGEOM_FORCE_INLINE VECCORE_ATT_HOST_DEVICE void Transformation3DMP<Real_s>::DoTranslation(
0762     Vector3D<Real_s> const &master, Vector3D<Real_s> &local) const
0763 {
0764 
0765   local[0] = master[0] - tx_;
0766   local[1] = master[1] - ty_;
0767   local[2] = master[2] - tz_;
0768 }
0769 
0770 /**
0771  * Transform a point to the local reference frame.
0772  * \param master Point to be transformed.
0773  * \param local Output destination. Should never be the same as the input
0774  *              vector!
0775  */
0776 template <typename Real_s>
0777 VECGEOM_FORCE_INLINE VECCORE_ATT_HOST_DEVICE void Transformation3DMP<Real_s>::Transform(Vector3D<Real_s> const &master,
0778                                                                                         Vector3D<Real_s> &local) const
0779 {
0780   Vector3D<Real_s> tmp;
0781   DoTranslation(master, tmp);
0782   DoRotation(tmp, local);
0783 }
0784 
0785 /**
0786  * Since transformation cannot be done in place, allows the transformed vector
0787  * to be constructed by Transform directly.
0788  * \param master Point to be transformed.
0789  * \return Newly constructed Vector3D with the transformed coordinates.
0790  */
0791 template <typename Real_s>
0792 VECGEOM_FORCE_INLINE VECCORE_ATT_HOST_DEVICE Vector3D<Real_s> Transformation3DMP<Real_s>::Transform(
0793     Vector3D<Real_s> const &master) const
0794 {
0795 
0796   Vector3D<Real_s> local;
0797   Transform(master, local);
0798   return local;
0799 }
0800 
0801 template <typename Real_s>
0802 template <bool transform_direction>
0803 VECGEOM_FORCE_INLINE VECCORE_ATT_HOST_DEVICE void Transformation3DMP<Real_s>::InverseTransformKernel(
0804     Vector3D<Real_s> const &local, Vector3D<Real_s> &master) const
0805 {
0806 
0807   // we are just doing the full stuff here ( LocalToMaster is less critical
0808   // than other way round )
0809 
0810   if (transform_direction) {
0811     master[0] = local[0] * rxx_;
0812     master[0] += local[1] * ryx_;
0813     master[0] += local[2] * rzx_;
0814     master[1] = local[0] * rxy_;
0815     master[1] += local[1] * ryy_;
0816     master[1] += local[2] * rzy_;
0817     master[2] = local[0] * rxz_;
0818     master[2] += local[1] * ryz_;
0819     master[2] += local[2] * rzz_;
0820   } else {
0821     master[0] = tx_;
0822     master[0] += local[0] * rxx_;
0823     master[0] += local[1] * ryx_;
0824     master[0] += local[2] * rzx_;
0825     master[1] = ty_;
0826     master[1] += local[0] * rxy_;
0827     master[1] += local[1] * ryy_;
0828     master[1] += local[2] * rzy_;
0829     master[2] = tz_;
0830     master[2] += local[0] * rxz_;
0831     master[2] += local[1] * ryz_;
0832     master[2] += local[2] * rzz_;
0833   }
0834 }
0835 
0836 template <typename Real_s>
0837 VECGEOM_FORCE_INLINE VECCORE_ATT_HOST_DEVICE void Transformation3DMP<Real_s>::InverseTransform(
0838     Vector3D<Real_s> const &local, Vector3D<Real_s> &master) const
0839 {
0840   InverseTransformKernel<false>(local, master);
0841 }
0842 
0843 template <typename Real_s>
0844 VECGEOM_FORCE_INLINE VECCORE_ATT_HOST_DEVICE Vector3D<Real_s> Transformation3DMP<Real_s>::InverseTransform(
0845     Vector3D<Real_s> const &local) const
0846 {
0847   Vector3D<Real_s> tmp;
0848   InverseTransform(local, tmp);
0849   return tmp;
0850 }
0851 
0852 template <typename Real_s>
0853 VECGEOM_FORCE_INLINE VECCORE_ATT_HOST_DEVICE void Transformation3DMP<Real_s>::InverseTransformDirection(
0854     Vector3D<Real_s> const &local, Vector3D<Real_s> &master) const
0855 {
0856   InverseTransformKernel<true>(local, master);
0857 }
0858 
0859 template <typename Real_s>
0860 VECGEOM_FORCE_INLINE VECCORE_ATT_HOST_DEVICE Vector3D<Real_s> Transformation3DMP<Real_s>::InverseTransformDirection(
0861     Vector3D<Real_s> const &local) const
0862 {
0863   Vector3D<Real_s> tmp;
0864   InverseTransformDirection(local, tmp);
0865   return tmp;
0866 }
0867 
0868 template <typename Real_s>
0869 VECCORE_ATT_HOST_DEVICE VECGEOM_FORCE_INLINE Transformation3DMP<Real_s> Transformation3DMP<Real_s>::operator*(
0870     Transformation3DMP<Real_s> const &rhs) const
0871 {
0872   if (rhs.fIdentity) return Transformation3DMP<Real_s>(*this);
0873   return Transformation3DMP<Real_s>(tx_ * rhs.rxx_ + ty_ * rhs.ryx_ + tz_ * rhs.rzx_ + rhs.tx_, // tx
0874                                     tx_ * rhs.rxy_ + ty_ * rhs.ryy_ + tz_ * rhs.rzy_ + rhs.ty_, // ty
0875                                     tx_ * rhs.rxz_ + ty_ * rhs.ryz_ + tz_ * rhs.rzz_ + rhs.tz_, // tz
0876                                     rxx_ * rhs.rxx_ + rxy_ * rhs.ryx_ + rxz_ * rhs.rzx_,        // rxx
0877                                     ryx_ * rhs.rxx_ + ryy_ * rhs.ryx_ + ryz_ * rhs.rzx_,        // ryx
0878                                     rzx_ * rhs.rxx_ + rzy_ * rhs.ryx_ + rzz_ * rhs.rzx_,        // rzx
0879                                     rxx_ * rhs.rxy_ + rxy_ * rhs.ryy_ + rxz_ * rhs.rzy_,        // rxy
0880                                     ryx_ * rhs.rxy_ + ryy_ * rhs.ryy_ + ryz_ * rhs.rzy_,        // ryy
0881                                     rzx_ * rhs.rxy_ + rzy_ * rhs.ryy_ + rzz_ * rhs.rzy_,        // rzy
0882                                     rxx_ * rhs.rxz_ + rxy_ * rhs.ryz_ + rxz_ * rhs.rzz_,        // rxz
0883                                     ryx_ * rhs.rxz_ + ryy_ * rhs.ryz_ + ryz_ * rhs.rzz_,        // ryz
0884                                     rzx_ * rhs.rxz_ + rzy_ * rhs.ryz_ + rzz_ * rhs.rzz_);       // rzz
0885 }
0886 
0887 template <typename Real_s>
0888 VECCORE_ATT_HOST_DEVICE VECGEOM_FORCE_INLINE Transformation3DMP<Real_s> Transformation3DMP<Real_s>::operator*(
0889     Transformation3D const &rhs) const
0890 {
0891   if (rhs.IsIdentity()) return Transformation3DMP<Real_s>(*this);
0892   return Transformation3DMP<Real_s>(
0893       tx_ * rhs.Rotation()[0] + ty_ * rhs.Rotation()[3] + tz_ * rhs.Rotation()[6] + rhs.Translation(0),
0894       tx_ * rhs.Rotation()[1] + ty_ * rhs.Rotation()[4] + tz_ * rhs.Rotation()[7] + rhs.Translation(1),
0895       tx_ * rhs.Rotation()[2] + ty_ * rhs.Rotation()[5] + tz_ * rhs.Rotation()[8] + rhs.Translation(2),
0896       rxx_ * rhs.Rotation()[0] + rxy_ * rhs.Rotation()[3] + rxz_ * rhs.Rotation()[6],
0897       ryx_ * rhs.Rotation()[0] + ryy_ * rhs.Rotation()[3] + ryz_ * rhs.Rotation()[6],
0898       rzx_ * rhs.Rotation()[0] + rzy_ * rhs.Rotation()[3] + rzz_ * rhs.Rotation()[6],
0899       rxx_ * rhs.Rotation()[1] + rxy_ * rhs.Rotation()[4] + rxz_ * rhs.Rotation()[7],
0900       ryx_ * rhs.Rotation()[1] + ryy_ * rhs.Rotation()[4] + ryz_ * rhs.Rotation()[7],
0901       rzx_ * rhs.Rotation()[1] + rzy_ * rhs.Rotation()[4] + rzz_ * rhs.Rotation()[7],
0902       rxx_ * rhs.Rotation()[2] + rxy_ * rhs.Rotation()[5] + rxz_ * rhs.Rotation()[8],
0903       ryx_ * rhs.Rotation()[2] + ryy_ * rhs.Rotation()[5] + ryz_ * rhs.Rotation()[8],
0904       rzx_ * rhs.Rotation()[2] + rzy_ * rhs.Rotation()[5] + rzz_ * rhs.Rotation()[8]);
0905 }
0906 
0907 template <typename Real_s>
0908 VECCORE_ATT_HOST_DEVICE VECGEOM_FORCE_INLINE Transformation3DMP<Real_s> const &Transformation3DMP<Real_s>::operator*=(
0909     Transformation3DMP<Real_s> const &rhs)
0910 {
0911   if (rhs.fIdentity) return *this;
0912   fIdentity = false;
0913 
0914   fHasTranslation |= rhs.HasTranslation();
0915   if (fHasTranslation) {
0916     auto tx = tx_ * rhs.rxx_ + ty_ * rhs.ryx_ + tz_ * rhs.rzx_ + rhs.tx_;
0917     auto ty = tx_ * rhs.rxy_ + ty_ * rhs.ryy_ + tz_ * rhs.rzy_ + rhs.ty_;
0918     auto tz = tx_ * rhs.rxz_ + ty_ * rhs.ryz_ + tz_ * rhs.rzz_ + rhs.tz_;
0919     tx_     = tx;
0920     ty_     = ty;
0921     tz_     = tz;
0922   }
0923 
0924   fHasRotation |= rhs.HasRotation();
0925   if (rhs.HasRotation()) {
0926     auto rxx = rxx_ * rhs.rxx_ + rxy_ * rhs.ryx_ + rxz_ * rhs.rzx_;
0927     auto ryx = ryx_ * rhs.rxx_ + ryy_ * rhs.ryx_ + ryz_ * rhs.rzx_;
0928     auto rzx = rzx_ * rhs.rxx_ + rzy_ * rhs.ryx_ + rzz_ * rhs.rzx_;
0929     auto rxy = rxx_ * rhs.rxy_ + rxy_ * rhs.ryy_ + rxz_ * rhs.rzy_;
0930     auto ryy = ryx_ * rhs.rxy_ + ryy_ * rhs.ryy_ + ryz_ * rhs.rzy_;
0931     auto rzy = rzx_ * rhs.rxy_ + rzy_ * rhs.ryy_ + rzz_ * rhs.rzy_;
0932     auto rxz = rxx_ * rhs.rxz_ + rxy_ * rhs.ryz_ + rxz_ * rhs.rzz_;
0933     auto ryz = ryx_ * rhs.rxz_ + ryy_ * rhs.ryz_ + ryz_ * rhs.rzz_;
0934     auto rzz = rzx_ * rhs.rxz_ + rzy_ * rhs.ryz_ + rzz_ * rhs.rzz_;
0935 
0936     rxx_ = rxx;
0937     rxy_ = rxy;
0938     rxz_ = rxz;
0939     ryx_ = ryx;
0940     ryy_ = ryy;
0941     ryz_ = ryz;
0942     rzx_ = rzx;
0943     rzy_ = rzy;
0944     rzz_ = rzz;
0945   }
0946 
0947   return *this;
0948 }
0949 
0950 // This function is slow due to static cast and must only be used at construction but never at for propagation of rays!
0951 // It is only used for interfacing Transformation3D's which might be obtained from solids directly
0952 // This is needed e.g., in TopMatrixImpl in NavStateIndex.h
0953 template <typename Real_s>
0954 VECCORE_ATT_HOST_DEVICE VECGEOM_FORCE_INLINE Transformation3DMP<Real_s> const &Transformation3DMP<Real_s>::operator*=(
0955     Transformation3D const &rhs)
0956 {
0957   if (rhs.IsIdentity()) return *this;
0958   fIdentity = false;
0959 
0960   fHasTranslation |= rhs.HasTranslation();
0961   if (fHasTranslation) {
0962     auto tx = tx_ * rhs.Rotation()[0] + ty_ * rhs.Rotation()[1] + tz_ * rhs.Rotation()[2] + rhs.Translation(0);
0963     auto ty = tx_ * rhs.Rotation()[3] + ty_ * rhs.Rotation()[4] + tz_ * rhs.Rotation()[5] + rhs.Translation(1);
0964     auto tz = tx_ * rhs.Rotation()[6] + ty_ * rhs.Rotation()[7] + tz_ * rhs.Rotation()[8] + rhs.Translation(2);
0965     tx_     = static_cast<Real_s>(tx);
0966     ty_     = static_cast<Real_s>(ty);
0967     tz_     = static_cast<Real_s>(tz);
0968   }
0969 
0970   fHasRotation |= rhs.HasRotation();
0971   if (rhs.HasRotation()) {
0972     auto rxx = rxx_ * rhs.Rotation()[0] + rxy_ * rhs.Rotation()[1] + rxz_ * rhs.Rotation()[2];
0973     auto ryx = ryx_ * rhs.Rotation()[0] + ryy_ * rhs.Rotation()[1] + ryz_ * rhs.Rotation()[2];
0974     auto rzx = rzx_ * rhs.Rotation()[0] + rzy_ * rhs.Rotation()[1] + rzz_ * rhs.Rotation()[2];
0975     auto rxy = rxx_ * rhs.Rotation()[3] + rxy_ * rhs.Rotation()[4] + rxz_ * rhs.Rotation()[5];
0976     auto ryy = ryx_ * rhs.Rotation()[3] + ryy_ * rhs.Rotation()[4] + ryz_ * rhs.Rotation()[5];
0977     auto rzy = rzx_ * rhs.Rotation()[3] + rzy_ * rhs.Rotation()[4] + rzz_ * rhs.Rotation()[5];
0978     auto rxz = rxx_ * rhs.Rotation()[6] + rxy_ * rhs.Rotation()[7] + rxz_ * rhs.Rotation()[8];
0979     auto ryz = ryx_ * rhs.Rotation()[6] + ryy_ * rhs.Rotation()[7] + ryz_ * rhs.Rotation()[8];
0980     auto rzz = rzx_ * rhs.Rotation()[6] + rzy_ * rhs.Rotation()[7] + rzz_ * rhs.Rotation()[8];
0981 
0982     rxx_ = static_cast<Real_s>(rxx);
0983     rxy_ = static_cast<Real_s>(rxy);
0984     rxz_ = static_cast<Real_s>(rxz);
0985     ryx_ = static_cast<Real_s>(ryx);
0986     ryy_ = static_cast<Real_s>(ryy);
0987     ryz_ = static_cast<Real_s>(ryz);
0988     rzx_ = static_cast<Real_s>(rzx);
0989     rzy_ = static_cast<Real_s>(rzy);
0990     rzz_ = static_cast<Real_s>(rzz);
0991   }
0992 
0993   return *this;
0994 }
0995 
0996 template <typename Real_s>
0997 VECCORE_ATT_HOST_DEVICE VECGEOM_FORCE_INLINE void Transformation3DMP<Real_s>::MultiplyFromRight(
0998     Transformation3DMP<Real_s> const &rhs)
0999 {
1000   // TODO: this code should directly operator on Vector3D and Matrix3D
1001 
1002   if (rhs.fIdentity) return;
1003   fIdentity = false;
1004 
1005   if (rhs.HasTranslation()) {
1006     fHasTranslation = true;
1007     // ideal for fused multiply add
1008     tx_ += rxx_ * rhs.tx_;
1009     tx_ += ryx_ * rhs.ty_;
1010     tx_ += rzx_ * rhs.tz_;
1011 
1012     ty_ += rxy_ * rhs.tx_;
1013     ty_ += ryy_ * rhs.ty_;
1014     ty_ += rzy_ * rhs.tz_;
1015 
1016     tz_ += rxz_ * rhs.tx_;
1017     tz_ += ryz_ * rhs.ty_;
1018     tz_ += rzz_ * rhs.tz_;
1019   }
1020 
1021   if (rhs.HasRotation()) {
1022     fHasRotation   = true;
1023     Precision tmpx = rxx_;
1024     Precision tmpy = ryx_;
1025     Precision tmpz = rzx_;
1026 
1027     // first row of matrix
1028     rxx_ = tmpx * rhs.rxx_;
1029     ryx_ = tmpx * rhs.ryx_;
1030     rzx_ = tmpx * rhs.rzx_;
1031     rxx_ += tmpy * rhs.rxy_;
1032     ryx_ += tmpy * rhs.ryy_;
1033     rzx_ += tmpy * rhs.rzy_;
1034     rxx_ += tmpz * rhs.rxz_;
1035     ryx_ += tmpz * rhs.ryz_;
1036     rzx_ += tmpz * rhs.rzz_;
1037 
1038     tmpx = rxy_;
1039     tmpy = ryy_;
1040     tmpz = rzy_;
1041 
1042     // second row of matrix
1043     rxy_ = tmpx * rhs.rxx_;
1044     ryy_ = tmpx * rhs.ryx_;
1045     rzy_ = tmpx * rhs.rzx_;
1046     rxy_ += tmpy * rhs.rxy_;
1047     ryy_ += tmpy * rhs.ryy_;
1048     rzy_ += tmpy * rhs.rzy_;
1049     rxy_ += tmpz * rhs.rxz_;
1050     ryy_ += tmpz * rhs.ryz_;
1051     rzy_ += tmpz * rhs.rzz_;
1052 
1053     tmpx = rxz_;
1054     tmpy = ryz_;
1055     tmpz = rzz_;
1056 
1057     // third row of matrix
1058     rxz_ = tmpx * rhs.rxx_;
1059     ryz_ = tmpx * rhs.ryx_;
1060     rzz_ = tmpx * rhs.rzx_;
1061     rxz_ += tmpy * rhs.rxy_;
1062     ryz_ += tmpy * rhs.ryy_;
1063     rzz_ += tmpy * rhs.rzy_;
1064     rxz_ += tmpz * rhs.rxz_;
1065     ryz_ += tmpz * rhs.ryz_;
1066     rzz_ += tmpz * rhs.rzz_;
1067   }
1068 }
1069 
1070 // multiply from right with Transformation3D. This will be slow and should only be called at initialization, not during
1071 // the computation!
1072 template <typename Real_s>
1073 VECCORE_ATT_HOST_DEVICE VECGEOM_FORCE_INLINE void Transformation3DMP<Real_s>::MultiplyFromRight(
1074     Transformation3D const &rhs)
1075 {
1076   // TODO: this code should directly operator on Vector3D and Matrix3D
1077 
1078   if (rhs.IsIdentity()) return;
1079   fIdentity = false;
1080 
1081   if (rhs.HasTranslation()) {
1082     fHasTranslation = true;
1083     // ideal for fused multiply add
1084     tx_ += rxx_ * static_cast<Real_s>(rhs.Translation(0));
1085     tx_ += ryx_ * static_cast<Real_s>(rhs.Translation(1));
1086     tx_ += rzx_ * static_cast<Real_s>(rhs.Translation(2));
1087 
1088     ty_ += rxy_ * static_cast<Real_s>(rhs.Translation(0));
1089     ty_ += ryy_ * static_cast<Real_s>(rhs.Translation(1));
1090     ty_ += rzy_ * static_cast<Real_s>(rhs.Translation(2));
1091 
1092     tz_ += rxz_ * static_cast<Real_s>(rhs.Translation(0));
1093     tz_ += ryz_ * static_cast<Real_s>(rhs.Translation(1));
1094     tz_ += rzz_ * static_cast<Real_s>(rhs.Translation(2));
1095   }
1096 
1097   if (rhs.HasRotation()) {
1098     fHasRotation = true;
1099     Real_s tmpx  = rxx_;
1100     Real_s tmpy  = ryx_;
1101     Real_s tmpz  = rzx_;
1102 
1103     // first row of matrix
1104     rxx_ = tmpx * static_cast<Real_s>(rhs.Rotation()[0]);
1105     ryx_ = tmpx * static_cast<Real_s>(rhs.Rotation()[1]);
1106     rzx_ = tmpx * static_cast<Real_s>(rhs.Rotation()[2]);
1107     rxx_ += tmpy * static_cast<Real_s>(rhs.Rotation()[3]);
1108     ryx_ += tmpy * static_cast<Real_s>(rhs.Rotation()[4]);
1109     rzx_ += tmpy * static_cast<Real_s>(rhs.Rotation()[5]);
1110     rxx_ += tmpz * static_cast<Real_s>(rhs.Rotation()[6]);
1111     ryx_ += tmpz * static_cast<Real_s>(rhs.Rotation()[7]);
1112     rzx_ += tmpz * static_cast<Real_s>(rhs.Rotation()[8]);
1113 
1114     tmpx = rxy_;
1115     tmpy = ryy_;
1116     tmpz = rzy_;
1117 
1118     // second row of matrix
1119     rxy_ = tmpx * static_cast<Real_s>(rhs.Rotation()[0]);
1120     ryy_ = tmpx * static_cast<Real_s>(rhs.Rotation()[1]);
1121     rzy_ = tmpx * static_cast<Real_s>(rhs.Rotation()[2]);
1122     rxy_ += tmpy * static_cast<Real_s>(rhs.Rotation()[3]);
1123     ryy_ += tmpy * static_cast<Real_s>(rhs.Rotation()[4]);
1124     rzy_ += tmpy * static_cast<Real_s>(rhs.Rotation()[5]);
1125     rxy_ += tmpz * static_cast<Real_s>(rhs.Rotation()[6]);
1126     ryy_ += tmpz * static_cast<Real_s>(rhs.Rotation()[7]);
1127     rzy_ += tmpz * static_cast<Real_s>(rhs.Rotation()[8]);
1128 
1129     tmpx = rxz_;
1130     tmpy = ryz_;
1131     tmpz = rzz_;
1132 
1133     // third row of matrix
1134     rxz_ = tmpx * static_cast<Real_s>(rhs.Rotation()[0]);
1135     ryz_ = tmpx * static_cast<Real_s>(rhs.Rotation()[1]);
1136     rzz_ = tmpx * static_cast<Real_s>(rhs.Rotation()[2]);
1137     rxz_ += tmpy * static_cast<Real_s>(rhs.Rotation()[3]);
1138     ryz_ += tmpy * static_cast<Real_s>(rhs.Rotation()[4]);
1139     rzz_ += tmpy * static_cast<Real_s>(rhs.Rotation()[5]);
1140     rxz_ += tmpz * static_cast<Real_s>(rhs.Rotation()[6]);
1141     ryz_ += tmpz * static_cast<Real_s>(rhs.Rotation()[7]);
1142     rzz_ += tmpz * static_cast<Real_s>(rhs.Rotation()[8]);
1143   }
1144 }
1145 
1146 /**
1147  * Only transforms by rotation, ignoring the translation part. This is useful
1148  * when transforming directions.
1149  * \param master Point to be transformed.
1150  * \param local Output destination of transformation.
1151  */
1152 template <typename Real_s>
1153 VECGEOM_FORCE_INLINE VECCORE_ATT_HOST_DEVICE void Transformation3DMP<Real_s>::TransformDirection(
1154     Vector3D<Real_s> const &master, Vector3D<Real_s> &local) const
1155 {
1156   DoRotation(master, local);
1157 }
1158 
1159 /**
1160  * Since transformation cannot be done in place, allows the transformed vector
1161  * to be constructed by TransformDirection directly.
1162  * \param master Point to be transformed.
1163  * \return Newly constructed Vector3D with the transformed coordinates.
1164  */
1165 template <typename Real_s>
1166 VECGEOM_FORCE_INLINE VECCORE_ATT_HOST_DEVICE Vector3D<Real_s> Transformation3DMP<Real_s>::TransformDirection(
1167     Vector3D<Real_s> const &master) const
1168 {
1169 
1170   Vector3D<Real_s> local;
1171   TransformDirection(master, local);
1172   return local;
1173 }
1174 
1175 template <typename Real_s>
1176 std::ostream &operator<<(std::ostream &os, Transformation3DMP<Real_s> const &trans)
1177 {
1178   os << "TransformationMP {" << trans.Translation() << ", " << "(" << trans.Rotation(0) << ", " << trans.Rotation(1)
1179      << ", " << trans.Rotation(2) << ", " << trans.Rotation(3) << ", " << trans.Rotation(4) << ", " << trans.Rotation(5)
1180      << ", " << trans.Rotation(6) << ", " << trans.Rotation(7) << ", " << trans.Rotation(8) << ")}" << "; identity("
1181      << trans.IsIdentity() << "); rotation(" << trans.HasRotation() << ")";
1182   return os;
1183 }
1184 }
1185 } // namespace vecgeom
1186 
1187 #endif // VECGEOM_BASE_TRANSFORMATION3DMP_H_
[ Source navigation ] [ Diff markup ] [ Identifier search ] [ general search ]

This page was automatically generated by the 2.3.7 LXR engine.

The LXR team
Valid CSS 2.1! Valid HTML 4.01!