Back to home page

EIC code displayed by LXR
 
 

    

File indexing completed on 2025-09-17 09:18:16 

0001 /// \file Transformation3D.h
0002 /// \author Johannes de Fine Licht (johannes.definelicht@cern.ch)
0003 
0004 #ifndef VECGEOM_BASE_TRANSFORMATION3D_H_
0005 #define VECGEOM_BASE_TRANSFORMATION3D_H_
0006 
0007 #include "VecGeom/base/Config.h"
0008 #include "VecGeom/base/Cuda.h"
0009 #include "VecGeom/base/Global.h"
0010 
0011 #include "VecGeom/base/Vector3D.h"
0012 #include "VecGeom/backend/scalar/Backend.h"
0013 #ifdef VECGEOM_ENABLE_CUDA
0014 #include "VecGeom/backend/cuda/Interface.h"
0015 #endif
0016 #include "VecGeom/management/Logger.h"
0017 
0018 #include <algorithm>
0019 #include <cmath>
0020 #include <cstring>
0021 #include <iostream>
0022 #include <vector>
0023 
0024 #ifdef VECGEOM_ROOT
0025 class TGeoMatrix;
0026 #endif
0027 
0028 typedef int RotationCode;
0029 typedef int TranslationCode;
0030 
0031 namespace vecgeom {
0032 namespace rotation {
0033 enum RotationId { kGeneric = -1, kDiagonal = 0x111, kIdentity = 0x200 };
0034 }
0035 namespace translation {
0036 enum TranslationId { kGeneric = -1, kIdentity = 0 };
0037 }
0038 } // namespace vecgeom
0039 
0040 namespace vecgeom {
0041 
0042 VECGEOM_DEVICE_FORWARD_DECLARE(class Transformation3D;);
0043 
0044 inline namespace VECGEOM_IMPL_NAMESPACE {
0045 
0046 #ifndef VECCORE_CUDA
0047 }
0048 namespace cuda {
0049 class Transformation3D;
0050 }
0051 inline namespace VECGEOM_IMPL_NAMESPACE {
0052 // class vecgeom::cuda::Transformation3D;
0053 #endif
0054 
0055 class Transformation3D {
0056 
0057 private:
0058   // TODO: it might be better to directly store this in terms of Vector3D<Precision> !!
0059   // and would allow for higher level abstraction
0060   Precision fTranslation[3];
0061   Precision fRotation[9];
0062   bool fIdentity;
0063   bool fHasRotation;
0064   bool fHasTranslation;
0065 
0066 public:
0067   VECCORE_ATT_HOST_DEVICE
0068   constexpr Transformation3D()
0069       : fTranslation{0., 0., 0.}, fRotation{1., 0., 0., 0., 1., 0., 0., 0., 1.}, fIdentity(true), fHasRotation(false),
0070         fHasTranslation(false){};
0071 
0072   /**
0073    * Constructor for translation only.
0074    * @param tx Translation in x-coordinate.
0075    * @param ty Translation in y-coordinate.
0076    * @param tz Translation in z-coordinate.
0077    */
0078   VECCORE_ATT_HOST_DEVICE
0079   Transformation3D(const Precision tx, const Precision ty, const Precision tz)
0080       : fTranslation{tx, ty, tz}, fRotation{1., 0., 0., 0., 1., 0., 0., 0., 1.},
0081         fIdentity(tx == 0 && ty == 0 && tz == 0), fHasRotation(false), fHasTranslation(tx != 0 || ty != 0 || tz != 0)
0082   {
0083   }
0084 
0085   /**
0086    * @brief Rotation followed by translation.
0087    * @param tx Translation in x-coordinate.
0088    * @param ty Translation in y-coordinate.
0089    * @param tz Translation in z-coordinate.
0090    * @param phi Rotation angle about z-axis.
0091    * @param theta Rotation angle about new y-axis.
0092    * @param psi Rotation angle about new z-axis.
0093    */
0094   VECCORE_ATT_HOST_DEVICE
0095   Transformation3D(const Precision tx, const Precision ty, const Precision tz, const Precision phi,
0096                    const Precision theta, const Precision psi);
0097 
0098   /**
0099    * @brief RScale followed by rotation followed by translation.
0100    * @param tx Translation in x-coordinate.
0101    * @param ty Translation in y-coordinate.
0102    * @param tz Translation in z-coordinate.
0103    * @param phi Rotation angle about z-axis.
0104    * @param theta Rotation angle about new y-axis.
0105    * @param psi Rotation angle about new z-axis.
0106    */
0107   VECCORE_ATT_HOST_DEVICE
0108   Transformation3D(const Precision tx, const Precision ty, const Precision tz, const Precision phi,
0109                    const Precision theta, const Precision psi, Precision sx, Precision sy, Precision sz);
0110 
0111 
0112   /**
0113    * Constructor to manually set each entry. Used when converting from different
0114    * geometry.
0115    */
0116   VECCORE_ATT_HOST_DEVICE
0117   Transformation3D(const Precision tx, const Precision ty, const Precision tz, const Precision r0, const Precision r1,
0118                    const Precision r2, const Precision r3, const Precision r4, const Precision r5, const Precision r6,
0119                    const Precision r7, const Precision r8);
0120 
0121   /**
0122    * Constructor to manually set each entry. Used when converting from different
0123    * geometry, supporting scaling.
0124    */
0125   VECCORE_ATT_HOST_DEVICE
0126   Transformation3D(const Precision tx, const Precision ty, const Precision tz, const Precision r0, const Precision r1,
0127                    const Precision r2, const Precision r3, const Precision r4, const Precision r5, const Precision r6,
0128                    const Precision r7, const Precision r8, Precision sx, Precision sy, Precision sz);
0129 
0130   /**
0131    * Constructor copying the translation and rotation from memory
0132    * geometry.
0133    */
0134   VECCORE_ATT_HOST_DEVICE
0135   Transformation3D(const Precision *trans, const Precision *rot, bool has_trans, bool has_rot)
0136   {
0137     this->Set(trans, rot, has_trans, has_rot);
0138   }
0139 
0140   /**
0141    * Constructor for a rotation based on a given direction
0142    * @param axis direction of the new z axis
0143    * @param inverse if true the origial axis will be rotated into (0,0,u)
0144                     if false a vector (0,0,u) will be rotated into the original axis
0145    */
0146   VECCORE_ATT_HOST_DEVICE
0147   Transformation3D(const Vector3D<Precision> &axis, bool inverse = true);
0148 
0149   VECCORE_ATT_HOST_DEVICE
0150   VECGEOM_FORCE_INLINE
0151   Transformation3D(Transformation3D const &other);
0152 
0153   VECCORE_ATT_HOST_DEVICE
0154   VECGEOM_FORCE_INLINE
0155   Transformation3D &operator=(Transformation3D const &rhs);
0156 
0157   VECCORE_ATT_HOST_DEVICE
0158   VECGEOM_FORCE_INLINE
0159   bool operator==(Transformation3D const &rhs) const;
0160 
0161   VECCORE_ATT_HOST_DEVICE
0162   ~Transformation3D() {}
0163 
0164   VECCORE_ATT_HOST_DEVICE
0165   void ApplyScale(const Precision sx, const Precision sy, const Precision sz)
0166   {
0167     fRotation[0] *= sx;
0168     fRotation[1] *= sy;
0169     fRotation[2] *= sz;
0170     fRotation[3] *= sx;
0171     fRotation[4] *= sy;
0172     fRotation[5] *= sz;
0173     fRotation[6] *= sx;
0174     fRotation[7] *= sy;
0175     fRotation[8] *= sz;
0176   }
0177 
0178   VECCORE_ATT_HOST_DEVICE
0179   void Clear()
0180   {
0181     fTranslation[0] = 0.;
0182     fTranslation[1] = 0.;
0183     fTranslation[2] = 0.;
0184     fRotation[0]    = 1.;
0185     fRotation[1]    = 0.;
0186     fRotation[2]    = 0.;
0187     fRotation[3]    = 0.;
0188     fRotation[4]    = 1.;
0189     fRotation[5]    = 0.;
0190     fRotation[6]    = 0.;
0191     fRotation[7]    = 0.;
0192     fRotation[8]    = 1.;
0193     fIdentity       = true;
0194     fHasRotation    = false;
0195     fHasTranslation = false;
0196   }
0197 
0198   int MemorySize() const { return sizeof(*this); }
0199 
0200   VECCORE_ATT_HOST_DEVICE
0201   void FixZeroes()
0202   {
0203     for (unsigned int i = 0; i < 9; ++i) {
0204       if (std::abs(fRotation[i]) < vecgeom::kTolerance) fRotation[i] = 0.;
0205     }
0206     for (unsigned int i = 0; i < 3; ++i) {
0207       if (std::abs(fTranslation[i]) < vecgeom::kTolerance) fTranslation[i] = 0.;
0208     }
0209   }
0210 
0211   VECCORE_ATT_HOST_DEVICE
0212   VECGEOM_FORCE_INLINE
0213   Vector3D<Precision> Translation() const
0214   {
0215     return Vector3D<Precision>(fTranslation[0], fTranslation[1], fTranslation[2]);
0216   }
0217 
0218   VECCORE_ATT_HOST_DEVICE
0219   VECGEOM_FORCE_INLINE
0220   Vector3D<Precision> Scale() const
0221   {
0222     Precision sx = std::sqrt(fRotation[0]*fRotation[0] + fRotation[3]*fRotation[3] + fRotation[6]*fRotation[6]);
0223     Precision sy = std::sqrt(fRotation[1]*fRotation[1] + fRotation[4]*fRotation[4] + fRotation[7]*fRotation[7]);
0224     Precision sz = std::sqrt(fRotation[2]*fRotation[2] + fRotation[5]*fRotation[5] + fRotation[8]*fRotation[8]);
0225 
0226     if (Determinant() < 0) sz = -sz;
0227 
0228     return Vector3D<Precision>(sx, sy, sz);
0229   }
0230 
0231   /**
0232    * No safety against faulty indexing.
0233    * @param index Index of translation entry in the range [0-2].
0234    */
0235   VECCORE_ATT_HOST_DEVICE
0236   VECGEOM_FORCE_INLINE
0237   Precision Translation(const int index) const { return fTranslation[index]; }
0238 
0239   VECCORE_ATT_HOST_DEVICE
0240   VECGEOM_FORCE_INLINE
0241   Precision const *Rotation() const { return fRotation; }
0242 
0243   /**
0244    * No safety against faulty indexing.
0245    * \param index Index of rotation entry in the range [0-8].
0246    */
0247   VECCORE_ATT_HOST_DEVICE
0248   VECGEOM_FORCE_INLINE
0249   Precision Rotation(const int index) const { return fRotation[index]; }
0250 
0251   VECCORE_ATT_HOST_DEVICE
0252   VECGEOM_FORCE_INLINE
0253   bool IsIdentity() const { return fIdentity; }
0254 
0255   VECCORE_ATT_HOST_DEVICE
0256   VECGEOM_FORCE_INLINE
0257   bool IsReflected() const { return Determinant() < 0; }
0258 
0259   VECCORE_ATT_HOST_DEVICE
0260   VECGEOM_FORCE_INLINE
0261   bool HasRotation() const { return fHasRotation; }
0262 
0263   VECCORE_ATT_HOST_DEVICE
0264   VECGEOM_FORCE_INLINE
0265   bool HasTranslation() const { return fHasTranslation; }
0266 
0267   VECCORE_ATT_HOST_DEVICE
0268   void Print() const;
0269 
0270   // print to a stream
0271   void Print(std::ostream &) const;
0272 
0273   // Mutators
0274 
0275   VECCORE_ATT_HOST_DEVICE
0276   void SetTranslation(const Precision tx, const Precision ty, const Precision tz);
0277 
0278   VECCORE_ATT_HOST_DEVICE
0279   void SetTranslation(Vector3D<Precision> const &vec);
0280 
0281   VECCORE_ATT_HOST_DEVICE
0282   void SetProperties();
0283 
0284   VECCORE_ATT_HOST_DEVICE
0285   void SetRotation(const Precision phi, const Precision theta, const Precision psi);
0286 
0287   VECCORE_ATT_HOST_DEVICE
0288   void SetRotation(Vector3D<Precision> const &vec);
0289 
0290   VECCORE_ATT_HOST_DEVICE
0291   void SetRotation(const Precision rot0, const Precision rot1, const Precision rot2, const Precision rot3,
0292                    const Precision rot4, const Precision rot5, const Precision rot6, const Precision rot7,
0293                    const Precision rot8);
0294 
0295   /**
0296    * Set transformation and rotation.
0297    * \param trans Pointer to at least 3 values.
0298    * \param rot Pointer to at least 9 values.
0299    */
0300   VECCORE_ATT_HOST_DEVICE
0301   VECGEOM_FORCE_INLINE
0302   void Set(const Precision *trans, const Precision *rot, bool has_trans, bool has_rot)
0303   {
0304     if (has_trans) {
0305       for (size_t i = 0; i < 3; ++i)
0306         fTranslation[i] = trans[i];
0307     }
0308 
0309     if (has_rot) {
0310       for (size_t i = 0; i < 9; ++i)
0311         fRotation[i] = rot[i];
0312     }
0313 
0314     fHasTranslation = has_trans;
0315     fHasRotation    = has_rot;
0316     fIdentity       = !fHasTranslation && !fHasRotation;
0317   }
0318 
0319   // Generation of template parameter codes
0320 
0321   VECCORE_ATT_HOST_DEVICE
0322   RotationCode GenerateRotationCode() const;
0323 
0324   VECCORE_ATT_HOST_DEVICE
0325   TranslationCode GenerateTranslationCode() const;
0326 
0327 private:
0328   // Templated rotation and translation methods which inline and compile to
0329   // optimized versions.
0330 
0331   template <RotationCode code, typename InputType>
0332   VECGEOM_FORCE_INLINE
0333   VECCORE_ATT_HOST_DEVICE
0334   void DoRotation(Vector3D<InputType> const &master, Vector3D<InputType> &local) const;
0335 
0336   template <RotationCode code, typename InputType>
0337   VECGEOM_FORCE_INLINE
0338   VECCORE_ATT_HOST_DEVICE
0339   void DoRotation_new(Vector3D<InputType> const &master, Vector3D<InputType> &local) const;
0340 
0341 private:
0342   template <typename InputType>
0343   VECGEOM_FORCE_INLINE
0344   VECCORE_ATT_HOST_DEVICE
0345   void DoTranslation(Vector3D<InputType> const &master, Vector3D<InputType> &local) const;
0346 
0347   template <bool vectortransform, typename InputType>
0348   VECGEOM_FORCE_INLINE
0349   VECCORE_ATT_HOST_DEVICE
0350   void InverseTransformKernel(Vector3D<InputType> const &local, Vector3D<InputType> &master) const;
0351 
0352 public:
0353   // Transformation interface
0354 
0355   template <TranslationCode trans_code, RotationCode rot_code, typename InputType>
0356   VECGEOM_FORCE_INLINE
0357   VECCORE_ATT_HOST_DEVICE
0358   void Transform(Vector3D<InputType> const &master, Vector3D<InputType> &local) const;
0359 
0360   template <TranslationCode trans_code, RotationCode rot_code, typename InputType>
0361   VECGEOM_FORCE_INLINE
0362   VECCORE_ATT_HOST_DEVICE
0363   Vector3D<InputType> Transform(Vector3D<InputType> const &master) const;
0364 
0365   template <typename InputType>
0366   VECGEOM_FORCE_INLINE
0367   VECCORE_ATT_HOST_DEVICE
0368   void Transform(Vector3D<InputType> const &master, Vector3D<InputType> &local) const;
0369 
0370   template <typename InputType>
0371   VECGEOM_FORCE_INLINE
0372   VECCORE_ATT_HOST_DEVICE
0373   Vector3D<InputType> Transform(Vector3D<InputType> const &master) const;
0374 
0375   template <RotationCode code, typename InputType>
0376   VECGEOM_FORCE_INLINE
0377   VECCORE_ATT_HOST_DEVICE
0378   void TransformDirection(Vector3D<InputType> const &master, Vector3D<InputType> &local) const;
0379 
0380   template <RotationCode code, typename InputType>
0381   VECGEOM_FORCE_INLINE
0382   VECCORE_ATT_HOST_DEVICE
0383   Vector3D<InputType> TransformDirection(Vector3D<InputType> const &master) const;
0384 
0385   template <typename InputType>
0386   VECGEOM_FORCE_INLINE
0387   VECCORE_ATT_HOST_DEVICE
0388   void TransformDirection(Vector3D<InputType> const &master, Vector3D<InputType> &local) const;
0389 
0390   template <typename InputType>
0391   VECGEOM_FORCE_INLINE
0392   VECCORE_ATT_HOST_DEVICE
0393   Vector3D<InputType> TransformDirection(Vector3D<InputType> const &master) const;
0394 
0395   /** The inverse transformation ( aka LocalToMaster ) of an object transform like a point
0396    *  this does not need to currently template on placement since such a transformation is much less used
0397    */
0398   template <typename InputType>
0399   VECGEOM_FORCE_INLINE
0400   VECCORE_ATT_HOST_DEVICE
0401   void InverseTransform(Vector3D<InputType> const &local, Vector3D<InputType> &master) const;
0402 
0403   template <typename InputType>
0404   VECGEOM_FORCE_INLINE
0405   VECCORE_ATT_HOST_DEVICE
0406   Vector3D<InputType> InverseTransform(Vector3D<InputType> const &local) const;
0407 
0408   /** The inverse transformation of an object transforming like a vector */
0409   template <typename InputType>
0410   VECGEOM_FORCE_INLINE
0411   VECCORE_ATT_HOST_DEVICE
0412   void InverseTransformDirection(Vector3D<InputType> const &master, Vector3D<InputType> &local) const;
0413 
0414   template <typename InputType>
0415   VECGEOM_FORCE_INLINE
0416   VECCORE_ATT_HOST_DEVICE
0417   Vector3D<InputType> InverseTransformDirection(Vector3D<InputType> const &master) const;
0418 
0419   /** compose transformations - multiply transformations */
0420   VECCORE_ATT_HOST_DEVICE
0421   VECGEOM_FORCE_INLINE
0422   void MultiplyFromRight(Transformation3D const &rhs);
0423 
0424   /** compose transformations - multiply transformations */
0425   VECCORE_ATT_HOST_DEVICE
0426   VECGEOM_FORCE_INLINE
0427   void CopyFrom(Transformation3D const &rhs)
0428   {
0429     // not sure this compiles under CUDA
0430     copy(&rhs, &rhs + 1, this);
0431   }
0432 
0433   VECCORE_ATT_HOST_DEVICE
0434   double Determinant() const
0435   {
0436     // Computes determinant in double precision
0437     double xx_ = fRotation[0];
0438     double zz_ = fRotation[8];
0439     double yy_ = fRotation[4];
0440     double xy_ = fRotation[1];
0441     double xz_ = fRotation[2];
0442     double yx_ = fRotation[3];
0443     double yz_ = fRotation[5];
0444     double zx_ = fRotation[6];
0445     double zy_ = fRotation[7];
0446 
0447     double detxx = yy_ * zz_ - yz_ * zy_;
0448     double detxy = yx_ * zz_ - yz_ * zx_;
0449     double detxz = yx_ * zy_ - yy_ * zx_;
0450     double det   = xx_ * detxx - xy_ * detxy + xz_ * detxz;
0451     return det;
0452   }
0453 
0454 
0455   // stores the inverse of this matrix into inverse
0456   // taken from CLHEP implementation
0457   VECCORE_ATT_HOST_DEVICE
0458   void Inverse(Transformation3D &inverse) const
0459   {
0460     double xx_ = fRotation[0];
0461     double zz_ = fRotation[8];
0462     double yy_ = fRotation[4];
0463     double xy_ = fRotation[1];
0464     double xz_ = fRotation[2];
0465     double yx_ = fRotation[3];
0466     double yz_ = fRotation[5];
0467     double zx_ = fRotation[6];
0468     double zy_ = fRotation[7];
0469     double dx_ = fTranslation[0];
0470     double dy_ = fTranslation[1];
0471     double dz_ = fTranslation[2];
0472 
0473     double detxx = yy_ * zz_ - yz_ * zy_;
0474     double detxy = yx_ * zz_ - yz_ * zx_;
0475     double detxz = yx_ * zy_ - yy_ * zx_;
0476     double det   = xx_ * detxx - xy_ * detxy + xz_ * detxz;
0477     if (det == 0) {
0478       VECGEOM_LOG(error) << "Transform3D::inverse failed: zero determinant";
0479     }
0480     det = 1. / det;
0481     detxx *= det;
0482     detxy *= det;
0483     detxz *= det;
0484     double detyx            = (xy_ * zz_ - xz_ * zy_) * det;
0485     double detyy            = (xx_ * zz_ - xz_ * zx_) * det;
0486     double detyz            = (xx_ * zy_ - xy_ * zx_) * det;
0487     double detzx            = (xy_ * yz_ - xz_ * yy_) * det;
0488     double detzy            = (xx_ * yz_ - xz_ * yx_) * det;
0489     double detzz            = (xx_ * yy_ - xy_ * yx_) * det;
0490     inverse.fRotation[0]    = detxx;
0491     inverse.fRotation[1]    = -detyx;
0492     inverse.fRotation[2]    = detzx;
0493     inverse.fTranslation[0] = -detxx * dx_ + detyx * dy_ - detzx * dz_;
0494     inverse.fRotation[3] = -detxy, inverse.fRotation[4] = detyy, inverse.fRotation[5] = -detzy,
0495     inverse.fTranslation[1] = detxy * dx_ - detyy * dy_ + detzy * dz_;
0496     inverse.fRotation[6] = detxz, inverse.fRotation[7] = -detyz, inverse.fRotation[8] = detzz,
0497     inverse.fTranslation[2] = -detxz * dx_ + detyz * dy_ - detzz * dz_;
0498 
0499     inverse.fHasTranslation = HasTranslation();
0500     inverse.fHasRotation    = HasRotation();
0501     inverse.fIdentity       = fIdentity;
0502   }
0503 
0504   // Utility and CUDA
0505 
0506 #ifdef VECGEOM_CUDA_INTERFACE
0507   size_t DeviceSizeOf() const { return DevicePtr<cuda::Transformation3D>::SizeOf(); }
0508   DevicePtr<cuda::Transformation3D> CopyToGpu() const;
0509   DevicePtr<cuda::Transformation3D> CopyToGpu(DevicePtr<cuda::Transformation3D> const gpu_ptr) const;
0510   static void CopyManyToGpu(const std::vector<Transformation3D const *> &trafos,
0511                             const std::vector<DevicePtr<cuda::Transformation3D>> &gpu_ptrs);
0512 #endif
0513 
0514 #ifdef VECGEOM_ROOT
0515   // function to convert this transformation to a TGeo transformation
0516   // mainly used for the benchmark comparisons with ROOT
0517   static TGeoMatrix *ConvertToTGeoMatrix(Transformation3D const&);
0518 #endif
0519 
0520 public:
0521   static const Transformation3D kIdentity;
0522 
0523 }; // End class Transformation3D
0524 
0525 VECCORE_ATT_HOST_DEVICE
0526 Transformation3D::Transformation3D(Transformation3D const &other)
0527     : fIdentity(false), fHasRotation(false), fHasTranslation(false)
0528 {
0529   *this = other;
0530 }
0531 
0532 VECCORE_ATT_HOST_DEVICE
0533 Transformation3D &Transformation3D::operator=(Transformation3D const &rhs)
0534 {
0535   copy(rhs.fTranslation, rhs.fTranslation + 3, fTranslation);
0536   copy(rhs.fRotation, rhs.fRotation + 9, fRotation);
0537   fIdentity       = rhs.fIdentity;
0538   fHasTranslation = rhs.fHasTranslation;
0539   fHasRotation    = rhs.fHasRotation;
0540   return *this;
0541 }
0542 
0543 VECCORE_ATT_HOST_DEVICE
0544 bool Transformation3D::operator==(Transformation3D const &rhs) const
0545 {
0546   return equal(fTranslation, fTranslation + 3, rhs.fTranslation) && equal(fRotation, fRotation + 9, rhs.fRotation);
0547 }
0548 
0549 /**
0550  * Rotates a vector to this transformation's frame of reference.
0551  * Templates on the RotationCode generated by GenerateTranslationCode() to
0552  * perform specialized rotation.
0553  * \sa GenerateTranslationCode()
0554  * \param master Vector in original frame of reference.
0555  * \param local Output vector rotated to the new frame of reference.
0556  */
0557 template <RotationCode code, typename InputType>
0558 VECGEOM_FORCE_INLINE
0559 VECCORE_ATT_HOST_DEVICE
0560 void Transformation3D::DoRotation(Vector3D<InputType> const &master, Vector3D<InputType> &local) const
0561 {
0562 
0563   if (code == 0x1B1) {
0564     local[0] = master[0] * fRotation[0];
0565     local[1] = master[1] * fRotation[4] + master[2] * fRotation[7];
0566     local[2] = master[1] * fRotation[5] + master[2] * fRotation[8];
0567     return;
0568   }
0569   if (code == 0x18E) {
0570     local[0] = master[1] * fRotation[3];
0571     local[1] = master[0] * fRotation[1] + master[2] * fRotation[7];
0572     local[2] = master[0] * fRotation[2] + master[2] * fRotation[8];
0573     return;
0574   }
0575   if (code == 0x076) {
0576     local[0] = master[2] * fRotation[6];
0577     local[1] = master[0] * fRotation[1] + master[1] * fRotation[4];
0578     local[2] = master[0] * fRotation[2] + master[1] * fRotation[5];
0579     return;
0580   }
0581   if (code == 0x16A) {
0582     local[0] = master[1] * fRotation[3] + master[2] * fRotation[6];
0583     local[1] = master[0] * fRotation[1];
0584     local[2] = master[2] * fRotation[5] + master[2] * fRotation[8];
0585     return;
0586   }
0587   if (code == 0x155) {
0588     local[0] = master[0] * fRotation[0] + master[2] * fRotation[6];
0589     local[1] = master[1] * fRotation[4];
0590     local[2] = master[0] * fRotation[2] + master[2] * fRotation[8];
0591     return;
0592   }
0593   if (code == 0x0AD) {
0594     local[0] = master[0] * fRotation[0] + master[1] * fRotation[3];
0595     local[1] = master[2] * fRotation[7];
0596     local[2] = master[0] * fRotation[2] + master[1] * fRotation[5];
0597     return;
0598   }
0599   if (code == 0x0DC) {
0600     local[0] = master[1] * fRotation[3] + master[2] * fRotation[6];
0601     local[1] = master[1] * fRotation[4] + master[2] * fRotation[7];
0602     local[2] = master[0] * fRotation[2];
0603     return;
0604   }
0605   if (code == 0x0E3) {
0606     local[0] = master[0] * fRotation[0] + master[2] * fRotation[6];
0607     local[1] = master[0] * fRotation[1] + master[2] * fRotation[7];
0608     local[2] = master[1] * fRotation[5];
0609     return;
0610   }
0611   if (code == 0x11B) {
0612     local[0] = master[0] * fRotation[0] + master[1] * fRotation[3];
0613     local[1] = master[0] * fRotation[1] + master[1] * fRotation[4];
0614     local[2] = master[2] * fRotation[8];
0615     return;
0616   }
0617   if (code == 0x0A1) {
0618     local[0] = master[0] * fRotation[0];
0619     local[1] = master[2] * fRotation[7];
0620     local[2] = master[1] * fRotation[5];
0621     return;
0622   }
0623   if (code == 0x10A) {
0624     local[0] = master[1] * fRotation[3];
0625     local[1] = master[0] * fRotation[1];
0626     local[2] = master[2] * fRotation[8];
0627     return;
0628   }
0629   if (code == 0x046) {
0630     local[0] = master[1] * fRotation[3];
0631     local[1] = master[2] * fRotation[7];
0632     local[2] = master[0] * fRotation[2];
0633     return;
0634   }
0635   if (code == 0x062) {
0636     local[0] = master[2] * fRotation[6];
0637     local[1] = master[0] * fRotation[1];
0638     local[2] = master[1] * fRotation[5];
0639     return;
0640   }
0641   if (code == 0x054) {
0642     local[0] = master[2] * fRotation[6];
0643     local[1] = master[1] * fRotation[4];
0644     local[2] = master[0] * fRotation[2];
0645     return;
0646   }
0647 
0648   // code = 0x111;
0649   if (code == rotation::kDiagonal) {
0650     local[0] = master[0] * fRotation[0];
0651     local[1] = master[1] * fRotation[4];
0652     local[2] = master[2] * fRotation[8];
0653     return;
0654   }
0655 
0656   // code = 0x200;
0657   if (code == rotation::kIdentity) {
0658     local = master;
0659     return;
0660   }
0661 
0662   // General case
0663   local[0] = master[0] * fRotation[0];
0664   local[1] = master[0] * fRotation[1];
0665   local[2] = master[0] * fRotation[2];
0666   local[0] += master[1] * fRotation[3];
0667   local[1] += master[1] * fRotation[4];
0668   local[2] += master[1] * fRotation[5];
0669   local[0] += master[2] * fRotation[6];
0670   local[1] += master[2] * fRotation[7];
0671   local[2] += master[2] * fRotation[8];
0672 }
0673 
0674 /**
0675  * Rotates a vector to this transformation's frame of reference.
0676  * Templates on the RotationCode generated by GenerateTranslationCode() to
0677  * perform specialized rotation.
0678  * \sa GenerateTranslationCode()
0679  * \param master Vector in original frame of reference.
0680  * \param local Output vector rotated to the new frame of reference.
0681  */
0682 template <RotationCode code, typename InputType>
0683 VECGEOM_FORCE_INLINE
0684 VECCORE_ATT_HOST_DEVICE
0685 void Transformation3D::DoRotation_new(Vector3D<InputType> const &master, Vector3D<InputType> &local) const
0686 {
0687 
0688   // code = 0x200;
0689   if (code == rotation::kIdentity) {
0690     local = master;
0691     return;
0692   }
0693 
0694   // General case
0695   local = Vector3D<InputType>(); // reset to zero -- any better way to do this???
0696   if (code & 0x001) {
0697     local[0] += master[0] * fRotation[0];
0698   }
0699   if (code & 0x002) {
0700     local[1] += master[0] * fRotation[1];
0701   }
0702   if (code & 0x004) {
0703     local[2] += master[0] * fRotation[2];
0704   }
0705   if (code & 0x008) {
0706     local[0] += master[1] * fRotation[3];
0707   }
0708   if (code & 0x010) {
0709     local[1] += master[1] * fRotation[4];
0710   }
0711   if (code & 0x020) {
0712     local[2] += master[1] * fRotation[5];
0713   }
0714   if (code & 0x040) {
0715     local[0] += master[2] * fRotation[6];
0716   }
0717   if (code & 0x080) {
0718     local[1] += master[2] * fRotation[7];
0719   }
0720   if (code & 0x100) {
0721     local[2] += master[2] * fRotation[8];
0722   }
0723 }
0724 
0725 template <typename InputType>
0726 VECGEOM_FORCE_INLINE
0727 VECCORE_ATT_HOST_DEVICE
0728 void Transformation3D::DoTranslation(Vector3D<InputType> const &master, Vector3D<InputType> &local) const
0729 {
0730 
0731   local[0] = master[0] - fTranslation[0];
0732   local[1] = master[1] - fTranslation[1];
0733   local[2] = master[2] - fTranslation[2];
0734 }
0735 
0736 /**
0737  * Transform a point to the local reference frame.
0738  * \param master Point to be transformed.
0739  * \param local Output destination. Should never be the same as the input
0740  *              vector!
0741  */
0742 template <TranslationCode trans_code, RotationCode rot_code, typename InputType>
0743 VECGEOM_FORCE_INLINE
0744 VECCORE_ATT_HOST_DEVICE
0745 void Transformation3D::Transform(Vector3D<InputType> const &master, Vector3D<InputType> &local) const
0746 {
0747 
0748   // Identity
0749   if (trans_code == translation::kIdentity && rot_code == rotation::kIdentity) {
0750     local = master;
0751     return;
0752   }
0753 
0754   // Only translation
0755   if (trans_code != translation::kIdentity && rot_code == rotation::kIdentity) {
0756     DoTranslation(master, local);
0757     return;
0758   }
0759 
0760   // Only rotation
0761   if (trans_code == translation::kIdentity && rot_code != rotation::kIdentity) {
0762     DoRotation<rot_code>(master, local);
0763     return;
0764   }
0765 
0766   // General case
0767   Vector3D<InputType> tmp;
0768   DoTranslation(master, tmp);
0769   DoRotation<rot_code>(tmp, local);
0770 }
0771 
0772 /**
0773  * Since transformation cannot be done in place, allows the transformed vector
0774  * to be constructed by Transform directly.
0775  * \param master Point to be transformed.
0776  * \return Newly constructed Vector3D with the transformed coordinates.
0777  */
0778 template <TranslationCode trans_code, RotationCode rot_code, typename InputType>
0779 VECGEOM_FORCE_INLINE
0780 VECCORE_ATT_HOST_DEVICE
0781 Vector3D<InputType> Transformation3D::Transform(Vector3D<InputType> const &master) const
0782 {
0783 
0784   Vector3D<InputType> local;
0785   Transform<trans_code, rot_code>(master, local);
0786   return local;
0787 }
0788 
0789 template <typename InputType>
0790 VECGEOM_FORCE_INLINE
0791 VECCORE_ATT_HOST_DEVICE
0792 void Transformation3D::Transform(Vector3D<InputType> const &master, Vector3D<InputType> &local) const
0793 {
0794   Transform<translation::kGeneric, rotation::kGeneric>(master, local);
0795 }
0796 template <typename InputType>
0797 VECGEOM_FORCE_INLINE
0798 VECCORE_ATT_HOST_DEVICE
0799 Vector3D<InputType> Transformation3D::Transform(Vector3D<InputType> const &master) const
0800 {
0801   return Transform<translation::kGeneric, rotation::kGeneric>(master);
0802 }
0803 
0804 template <bool transform_direction, typename InputType>
0805 VECGEOM_FORCE_INLINE
0806 VECCORE_ATT_HOST_DEVICE
0807 void Transformation3D::InverseTransformKernel(Vector3D<InputType> const &local, Vector3D<InputType> &master) const
0808 {
0809 
0810   // we are just doing the full stuff here ( LocalToMaster is less critical
0811   // than other way round )
0812 
0813   if (transform_direction) {
0814     master[0] = local[0] * fRotation[0];
0815     master[0] += local[1] * fRotation[1];
0816     master[0] += local[2] * fRotation[2];
0817     master[1] = local[0] * fRotation[3];
0818     master[1] += local[1] * fRotation[4];
0819     master[1] += local[2] * fRotation[5];
0820     master[2] = local[0] * fRotation[6];
0821     master[2] += local[1] * fRotation[7];
0822     master[2] += local[2] * fRotation[8];
0823   } else {
0824     master[0] = fTranslation[0];
0825     master[0] += local[0] * fRotation[0];
0826     master[0] += local[1] * fRotation[1];
0827     master[0] += local[2] * fRotation[2];
0828     master[1] = fTranslation[1];
0829     master[1] += local[0] * fRotation[3];
0830     master[1] += local[1] * fRotation[4];
0831     master[1] += local[2] * fRotation[5];
0832     master[2] = fTranslation[2];
0833     master[2] += local[0] * fRotation[6];
0834     master[2] += local[1] * fRotation[7];
0835     master[2] += local[2] * fRotation[8];
0836   }
0837 }
0838 
0839 template <typename InputType>
0840 VECGEOM_FORCE_INLINE
0841 VECCORE_ATT_HOST_DEVICE
0842 void Transformation3D::InverseTransform(Vector3D<InputType> const &local, Vector3D<InputType> &master) const
0843 {
0844   InverseTransformKernel<false, InputType>(local, master);
0845 }
0846 
0847 template <typename InputType>
0848 VECGEOM_FORCE_INLINE
0849 VECCORE_ATT_HOST_DEVICE
0850 Vector3D<InputType> Transformation3D::InverseTransform(Vector3D<InputType> const &local) const
0851 {
0852   Vector3D<InputType> tmp;
0853   InverseTransform(local, tmp);
0854   return tmp;
0855 }
0856 
0857 template <typename InputType>
0858 VECGEOM_FORCE_INLINE
0859 VECCORE_ATT_HOST_DEVICE
0860 void Transformation3D::InverseTransformDirection(Vector3D<InputType> const &local, Vector3D<InputType> &master) const
0861 {
0862   InverseTransformKernel<true, InputType>(local, master);
0863 }
0864 
0865 template <typename InputType>
0866 VECGEOM_FORCE_INLINE
0867 VECCORE_ATT_HOST_DEVICE
0868 Vector3D<InputType> Transformation3D::InverseTransformDirection(Vector3D<InputType> const &local) const
0869 {
0870   Vector3D<InputType> tmp;
0871   InverseTransformDirection(local, tmp);
0872   return tmp;
0873 }
0874 
0875 VECCORE_ATT_HOST_DEVICE
0876 VECGEOM_FORCE_INLINE
0877 void Transformation3D::MultiplyFromRight(Transformation3D const &rhs)
0878 {
0879   // TODO: this code should directly operator on Vector3D and Matrix3D
0880 
0881   if (rhs.fIdentity) return;
0882   fIdentity = false;
0883 
0884   if (rhs.HasTranslation()) {
0885     fHasTranslation = true;
0886     // ideal for fused multiply add
0887     fTranslation[0] += fRotation[0] * rhs.fTranslation[0];
0888     fTranslation[0] += fRotation[1] * rhs.fTranslation[1];
0889     fTranslation[0] += fRotation[2] * rhs.fTranslation[2];
0890 
0891     fTranslation[1] += fRotation[3] * rhs.fTranslation[0];
0892     fTranslation[1] += fRotation[4] * rhs.fTranslation[1];
0893     fTranslation[1] += fRotation[5] * rhs.fTranslation[2];
0894 
0895     fTranslation[2] += fRotation[6] * rhs.fTranslation[0];
0896     fTranslation[2] += fRotation[7] * rhs.fTranslation[1];
0897     fTranslation[2] += fRotation[8] * rhs.fTranslation[2];
0898   }
0899 
0900   if (rhs.HasRotation()) {
0901     fHasRotation   = true;
0902     Precision tmpx = fRotation[0];
0903     Precision tmpy = fRotation[1];
0904     Precision tmpz = fRotation[2];
0905 
0906     // first row of matrix
0907     fRotation[0] = tmpx * rhs.fRotation[0];
0908     fRotation[1] = tmpx * rhs.fRotation[1];
0909     fRotation[2] = tmpx * rhs.fRotation[2];
0910     fRotation[0] += tmpy * rhs.fRotation[3];
0911     fRotation[1] += tmpy * rhs.fRotation[4];
0912     fRotation[2] += tmpy * rhs.fRotation[5];
0913     fRotation[0] += tmpz * rhs.fRotation[6];
0914     fRotation[1] += tmpz * rhs.fRotation[7];
0915     fRotation[2] += tmpz * rhs.fRotation[8];
0916 
0917     tmpx = fRotation[3];
0918     tmpy = fRotation[4];
0919     tmpz = fRotation[5];
0920 
0921     // second row of matrix
0922     fRotation[3] = tmpx * rhs.fRotation[0];
0923     fRotation[4] = tmpx * rhs.fRotation[1];
0924     fRotation[5] = tmpx * rhs.fRotation[2];
0925     fRotation[3] += tmpy * rhs.fRotation[3];
0926     fRotation[4] += tmpy * rhs.fRotation[4];
0927     fRotation[5] += tmpy * rhs.fRotation[5];
0928     fRotation[3] += tmpz * rhs.fRotation[6];
0929     fRotation[4] += tmpz * rhs.fRotation[7];
0930     fRotation[5] += tmpz * rhs.fRotation[8];
0931 
0932     tmpx = fRotation[6];
0933     tmpy = fRotation[7];
0934     tmpz = fRotation[8];
0935 
0936     // third row of matrix
0937     fRotation[6] = tmpx * rhs.fRotation[0];
0938     fRotation[7] = tmpx * rhs.fRotation[1];
0939     fRotation[8] = tmpx * rhs.fRotation[2];
0940     fRotation[6] += tmpy * rhs.fRotation[3];
0941     fRotation[7] += tmpy * rhs.fRotation[4];
0942     fRotation[8] += tmpy * rhs.fRotation[5];
0943     fRotation[6] += tmpz * rhs.fRotation[6];
0944     fRotation[7] += tmpz * rhs.fRotation[7];
0945     fRotation[8] += tmpz * rhs.fRotation[8];
0946   }
0947 }
0948 
0949 /**
0950  * Only transforms by rotation, ignoring the translation part. This is useful
0951  * when transforming directions.
0952  * \param master Point to be transformed.
0953  * \param local Output destination of transformation.
0954  */
0955 template <RotationCode code, typename InputType>
0956 VECGEOM_FORCE_INLINE
0957 VECCORE_ATT_HOST_DEVICE
0958 void Transformation3D::TransformDirection(Vector3D<InputType> const &master, Vector3D<InputType> &local) const
0959 {
0960 
0961   // Rotational fIdentity
0962   if (code == rotation::kIdentity) {
0963     local = master;
0964     return;
0965   }
0966 
0967   // General case
0968   DoRotation<code>(master, local);
0969 }
0970 
0971 /**
0972  * Since transformation cannot be done in place, allows the transformed vector
0973  * to be constructed by TransformDirection directly.
0974  * \param master Point to be transformed.
0975  * \return Newly constructed Vector3D with the transformed coordinates.
0976  */
0977 template <RotationCode code, typename InputType>
0978 VECGEOM_FORCE_INLINE
0979 VECCORE_ATT_HOST_DEVICE
0980 Vector3D<InputType> Transformation3D::TransformDirection(Vector3D<InputType> const &master) const
0981 {
0982 
0983   Vector3D<InputType> local;
0984   TransformDirection<code>(master, local);
0985   return local;
0986 }
0987 
0988 template <typename InputType>
0989 VECGEOM_FORCE_INLINE
0990 VECCORE_ATT_HOST_DEVICE
0991 void Transformation3D::TransformDirection(Vector3D<InputType> const &master, Vector3D<InputType> &local) const
0992 {
0993   TransformDirection<rotation::kGeneric>(master, local);
0994 }
0995 
0996 template <typename InputType>
0997 VECGEOM_FORCE_INLINE
0998 VECCORE_ATT_HOST_DEVICE
0999 Vector3D<InputType> Transformation3D::TransformDirection(Vector3D<InputType> const &master) const
1000 {
1001   return TransformDirection<rotation::kGeneric>(master);
1002 }
1003 
1004 std::ostream &operator<<(std::ostream &os, Transformation3D const &trans);
1005 }
1006 } // namespace vecgeom
1007 
1008 #endif // VECGEOM_BASE_TRANSFORMATION3D_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!