Back to home page

EIC code displayed by LXR
 
 

    

File indexing completed on 2025-01-18 10:13:55 

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