VecGeom/base/Vector3D.h

0001 /// \file vector3d.h
0002 /// \author Johannes de Fine Licht (johannes.definelicht@cern.ch)
0003
0004 #ifndef VECGEOM_BASE_VECTOR3D_H_
0005 #define VECGEOM_BASE_VECTOR3D_H_
0006
0007 #include "VecGeom/base/Global.h"
0008 #include "VecGeom/base/AlignedBase.h"
0009
0010 #include <cstdlib>
0011 #include <ostream>
0012 #include <string>
0013
0014 namespace vecgeom {
0015
0016 VECGEOM_DEVICE_FORWARD_DECLARE(template <typename Type> class Vector3D;);
0017 VECGEOM_DEVICE_DECLARE_CONV_TEMPLATE(class, Vector3D, typename);
0018
0019 inline namespace VECGEOM_IMPL_NAMESPACE {
0020
0021 /**
0022  * @brief Three dimensional vector class supporting most arithmetic operations.
0023  * @details If vector acceleration is enabled, the scalar template instantiation
0024  *          will use vector instructions for operations when possible.
0025  */
0026 template <typename Type>
0027 class Vector3D : public AlignedBase {
0028
0029   typedef Vector3D<Type> VecType;
0030
0031 private:
0032   Type vec[3];
0033
0034 public:
0035   VECCORE_ATT_HOST_DEVICE
0036   VECGEOM_FORCE_INLINE
0037   Vector3D(const Type a, const Type b, const Type c)
0038   {
0039     vec[0] = a;
0040     vec[1] = b;
0041     vec[2] = c;
0042   }
0043
0044   VECCORE_ATT_HOST_DEVICE
0045   VECGEOM_FORCE_INLINE
0046   Vector3D()
0047   {
0048     vec[0] = 0;
0049     vec[1] = 0;
0050     vec[2] = 0;
0051   }
0052
0053   VECCORE_ATT_HOST_DEVICE
0054   VECGEOM_FORCE_INLINE
0055   Vector3D(const Type a)
0056   {
0057     vec[0] = a;
0058     vec[1] = a;
0059     vec[2] = a;
0060   }
0061
0062   template <typename TypeOther>
0063   VECGEOM_FORCE_INLINE
0064   VECCORE_ATT_HOST_DEVICE
0065   Vector3D(Vector3D<TypeOther> const &other)
0066   {
0067     vec[0] = other[0];
0068     vec[1] = other[1];
0069     vec[2] = other[2];
0070   }
0071
0072   /**
0073    * Constructs a vector from an std::string of the same format as output by the
0074    * "<<"-operator for outstreams.
0075    * @param str String formatted as "(%d, %d, %d)".
0076    */
0077   VECCORE_ATT_HOST
0078   Vector3D(std::string const &str)
0079   {
0080     int begin = 1, end = str.find(",");
0081     vec[0] = std::atof(str.substr(begin, end - begin).c_str());
0082     begin  = end + 2;
0083     end    = str.find(",", begin);
0084     vec[1] = std::atof(str.substr(begin, end - begin).c_str());
0085     begin  = end + 2;
0086     end    = str.find(")", begin);
0087     vec[2] = std::atof(str.substr(begin, end - begin).c_str());
0088   }
0089
0090   /**
0091    * Contains no check for correct indexing to avoid impairing performance.
0092    * @param index Index of content in the range [0-2].
0093    */
0094   VECCORE_ATT_HOST_DEVICE
0095   VECGEOM_FORCE_INLINE
0096   Type &operator[](const int index) { return vec[index]; }
0097
0098   /**
0099    * Contains no check for correct indexing to avoid impairing performance.
0100    * @param index Index of content in the range [0-2].
0101    */
0102   VECCORE_ATT_HOST_DEVICE
0103   VECGEOM_FORCE_INLINE
0104   Type const &operator[](const int index) const { return vec[index]; }
0105
0106   VECCORE_ATT_HOST_DEVICE
0107   VECGEOM_FORCE_INLINE
0108   Type &x() { return vec[0]; }
0109
0110   VECCORE_ATT_HOST_DEVICE
0111   VECGEOM_FORCE_INLINE
0112   Type const &x() const { return vec[0]; }
0113
0114   VECCORE_ATT_HOST_DEVICE
0115   VECGEOM_FORCE_INLINE
0116   Type &y() { return vec[1]; }
0117
0118   VECCORE_ATT_HOST_DEVICE
0119   VECGEOM_FORCE_INLINE
0120   Type const &y() const { return vec[1]; }
0121
0122   VECCORE_ATT_HOST_DEVICE
0123   VECGEOM_FORCE_INLINE
0124   Type &z() { return vec[2]; }
0125
0126   VECCORE_ATT_HOST_DEVICE
0127   VECGEOM_FORCE_INLINE
0128   Type const &z() const { return vec[2]; }
0129
0130   VECCORE_ATT_HOST_DEVICE
0131   void Set(Type const &a, Type const &b, Type const &c)
0132   {
0133     vec[0] = a;
0134     vec[1] = b;
0135     vec[2] = c;
0136   }
0137
0138   VECCORE_ATT_HOST_DEVICE
0139   void Set(const Type a) { Set(a, a, a); }
0140
0141   /// \return the length squared perpendicular to z direction
0142   VECCORE_ATT_HOST_DEVICE
0143   VECGEOM_FORCE_INLINE
0144   Type Perp2() const { return vec[0] * vec[0] + vec[1] * vec[1]; }
0145
0146   /// \return the length perpendicular to z direction
0147   VECCORE_ATT_HOST_DEVICE
0148   VECGEOM_FORCE_INLINE
0149   Type Perp() const { return Sqrt(Perp2()); }
0150
0151   /// The dot product of two Vector3D<T> objects
0152   /// \return T (where T is float, double, or various SIMD vector types)
0153   template <typename Type2>
0154   VECGEOM_FORCE_INLINE
0155   VECCORE_ATT_HOST_DEVICE
0156   static Type Dot(Vector3D<Type> const &left, Vector3D<Type2> const &right)
0157   {
0158     return left[0] * right[0] + left[1] * right[1] + left[2] * right[2];
0159   }
0160
0161   /// The dot product of two Vector3D<T> objects
0162   /// \return T (where T is float, double, or various SIMD vector types)
0163   template <typename Type2>
0164   VECGEOM_FORCE_INLINE
0165   VECCORE_ATT_HOST_DEVICE
0166   Type Dot(Vector3D<Type2> const &right) const
0167   {
0168     return Dot(*this, right);
0169   }
0170
0171   // For UVector3 compatibility. It is equal to normal multiplication.
0172   VECCORE_ATT_HOST_DEVICE
0173   VECGEOM_FORCE_INLINE
0174   VecType MultiplyByComponents(VecType const &other) const { return *this * other; }
0175
0176   /// \return Squared magnitude of the vector.
0177   VECCORE_ATT_HOST_DEVICE
0178   VECGEOM_FORCE_INLINE
0179   Type Mag2() const { return Dot(*this, *this); }
0180
0181   /// \return Magnitude of the vector.
0182   VECCORE_ATT_HOST_DEVICE
0183   VECGEOM_FORCE_INLINE
0184   Type Mag() const { return Sqrt(Mag2()); }
0185
0186   VECCORE_ATT_HOST_DEVICE
0187   VECGEOM_FORCE_INLINE
0188   Type Length() const { return Mag(); }
0189
0190   VECCORE_ATT_HOST_DEVICE
0191   VECGEOM_FORCE_INLINE
0192   Type Length2() const { return Mag2(); }
0193
0194   /// Normalizes the vector by dividing each entry by the length.
0195   /// \sa Vector3D::Length()
0196   VECCORE_ATT_HOST_DEVICE
0197   VECGEOM_FORCE_INLINE
0198   void Normalize() { *this *= (Type(1.) / Length()); }
0199
0200   VECCORE_ATT_HOST_DEVICE
0201   VECGEOM_FORCE_INLINE
0202   Vector3D<Type> Normalized() const { return Vector3D<Type>(*this) * (Type(1.) / Length()); }
0203
0204   // checks if vector is normalized
0205   // only reasonable to call with standard scalare usage
0206   VECCORE_ATT_HOST_DEVICE
0207   VECGEOM_FORCE_INLINE
0208   bool IsNormalized() const
0209   {
0210     // static_assert here that Type should be primitive type
0211     Precision norm = Mag2();
0212     return Type(1.) - vecgeom::kTolerance < norm && norm < Type(1.) + vecgeom::kTolerance;
0213   }
0214
0215   /// \return Azimuthal angle between -pi and pi.
0216   VECCORE_ATT_HOST_DEVICE
0217   VECGEOM_FORCE_INLINE
0218   Type Phi() const { return ATan2(vec[1], vec[0]); }
0219
0220   /// \return Polar angle between 0 and pi.
0221   VECCORE_ATT_HOST_DEVICE
0222   VECGEOM_FORCE_INLINE
0223   Type Theta() const { return ACos(vec[2] / Mag()); }
0224
0225   /// The cross (vector) product of two Vector3D<T> objects
0226   /// \return Type (where Type is float, double, or various SIMD vector types)
0227   template <class FirstType, class SecondType>
0228   VECGEOM_FORCE_INLINE
0229   VECCORE_ATT_HOST_DEVICE
0230   static Vector3D<Type> Cross(Vector3D<FirstType> const &left, Vector3D<SecondType> const &right)
0231   {
0232     return Vector3D<Type>(left[1] * right[2] - left[2] * right[1], left[2] * right[0] - left[0] * right[2],
0233                           left[0] * right[1] - left[1] * right[0]);
0234   }
0235
0236   /// The cross (vector) product of two Vector3D<T> objects
0237   /// \return Type (where Type is float, double, or various SIMD vector types)
0238   template <class OtherType>
0239   VECGEOM_FORCE_INLINE
0240   VECCORE_ATT_HOST_DEVICE
0241   Vector3D<Type> Cross(Vector3D<OtherType> const &right) const
0242   {
0243     return Cross<Type, OtherType>(*this, right);
0244   }
0245
0246   /// Maps each vector entry to a function that manipulates the entry type.
0247   /// \param f A function of type "Type f(const Type&)" to map over entries.
0248   VECCORE_ATT_HOST_DEVICE
0249   VECGEOM_FORCE_INLINE
0250   void Map(Type (*f)(const Type &))
0251   {
0252     vec[0] = f(vec[0]);
0253     vec[1] = f(vec[1]);
0254     vec[2] = f(vec[2]);
0255   }
0256
0257   VECCORE_ATT_HOST_DEVICE
0258   VECGEOM_FORCE_INLINE
0259   Vector3D<Type> Abs() const
0260   {
0261     return Vector3D<Type>(vecCore::math::Abs(vec[0]), vecCore::math::Abs(vec[1]), vecCore::math::Abs(vec[2]));
0262   }
0263
0264   template <typename BoolType>
0265   VECGEOM_FORCE_INLINE
0266   VECCORE_ATT_HOST_DEVICE
0267   void MaskedAssign(Vector3D<BoolType> const &condition, Vector3D<Type> const &value)
0268   {
0269     vec[0] = (condition[0]) ? value[0] : vec[0];
0270     vec[1] = (condition[1]) ? value[1] : vec[1];
0271     vec[2] = (condition[2]) ? value[2] : vec[2];
0272   }
0273
0274   VECCORE_ATT_HOST_DEVICE
0275   VECGEOM_FORCE_INLINE
0276   Type Min() const { return vecCore::math::Min(vec[0], vec[1], vec[2]); }
0277
0278   VECCORE_ATT_HOST_DEVICE
0279   VECGEOM_FORCE_INLINE
0280   Type Max() const { return vecCore::math::Max(vec[0], vec[1], vec[2]); }
0281
0282   VECCORE_ATT_HOST_DEVICE
0283   VECGEOM_FORCE_INLINE
0284   VecType Unit() const
0285   {
0286     const Type mag2 = Mag2();
0287     VecType output(*this);
0288     output /= Sqrt(mag2 + kMinimum);
0289     return output;
0290   }
0291
0292   VECCORE_ATT_HOST_DEVICE
0293   VECGEOM_FORCE_INLINE
0294   static VecType FromCylindrical(Type r, Type phi, Type z) { return VecType(r * cos(phi), r * sin(phi), z); }
0295
0296   VECCORE_ATT_HOST_DEVICE
0297   VECGEOM_FORCE_INLINE
0298   VecType &FixZeroes()
0299   {
0300     using vecCore::math::Abs;
0301     for (int i = 0; i < 3; ++i) {
0302       vecCore::MaskedAssign(vec[i], Abs(vec[i]) < kTolerance, Type(0.0));
0303     }
0304     return *this;
0305   }
0306
0307   // Inplace binary operators
0308
0309 #define VECTOR3D_TEMPLATE_INPLACE_BINARY_OP(OPERATOR)                                                       \
0310   VECCORE_ATT_HOST_DEVICE                                                                                   \
0311   VECGEOM_FORCE_INLINE                                                                                      \
0312   VecType &operator OPERATOR(const VecType &other)                                                          \
0313   {                                                                                                         \
0314     vec[0] OPERATOR other.vec[0];                                                                           \
0315     vec[1] OPERATOR other.vec[1];                                                                           \
0316     vec[2] OPERATOR other.vec[2];                                                                           \
0317     return *this;                                                                                           \
0318   }                                                                                                         \
0319   template <typename OtherType>                                                                             \
0320   VECCORE_ATT_HOST_DEVICE VECGEOM_FORCE_INLINE VecType &operator OPERATOR(const Vector3D<OtherType> &other) \
0321   {                                                                                                         \
0322     vec[0] OPERATOR other[0];                                                                               \
0323     vec[1] OPERATOR other[1];                                                                               \
0324     vec[2] OPERATOR other[2];                                                                               \
0325     return *this;                                                                                           \
0326   }                                                                                                         \
0327   VECCORE_ATT_HOST_DEVICE                                                                                   \
0328   VECGEOM_FORCE_INLINE                                                                                      \
0329   VecType &operator OPERATOR(const Type &scalar)                                                            \
0330   {                                                                                                         \
0331     vec[0] OPERATOR scalar;                                                                                 \
0332     vec[1] OPERATOR scalar;                                                                                 \
0333     vec[2] OPERATOR scalar;                                                                                 \
0334     return *this;                                                                                           \
0335   }
0336   VECTOR3D_TEMPLATE_INPLACE_BINARY_OP(+=)
0337   VECTOR3D_TEMPLATE_INPLACE_BINARY_OP(-=)
0338   VECTOR3D_TEMPLATE_INPLACE_BINARY_OP(*=)
0339   VECTOR3D_TEMPLATE_INPLACE_BINARY_OP(/=)
0340 #undef VECTOR3D_TEMPLATE_INPLACE_BINARY_OP
0341
0342   VECCORE_ATT_HOST_DEVICE
0343   VECGEOM_FORCE_INLINE
0344   operator bool() const { return vec[0] && vec[1] && vec[2]; }
0345 };
0346
0347 template <typename T>
0348 std::ostream &operator<<(std::ostream &os, Vector3D<T> const &vec)
0349 {
0350   os << "(" << vec[0] << ", " << vec[1] << ", " << vec[2] << ")";
0351   return os;
0352 }
0353
0354 #define VECTOR3D_BINARY_OP(OPERATOR, INPLACE)                                                                   \
0355   template <typename Type, typename OtherType>                                                                  \
0356   VECGEOM_FORCE_INLINE VECCORE_ATT_HOST_DEVICE Vector3D<Type> operator OPERATOR(const Vector3D<Type> &lhs,      \
0357                                                                                 const Vector3D<OtherType> &rhs) \
0358   {                                                                                                             \
0359     Vector3D<Type> result(lhs);                                                                                 \
0360     result INPLACE rhs;                                                                                         \
0361     return result;                                                                                              \
0362   }                                                                                                             \
0363   template <typename Type, typename ScalarType>                                                                 \
0364   VECGEOM_FORCE_INLINE VECCORE_ATT_HOST_DEVICE Vector3D<Type> operator OPERATOR(Vector3D<Type> const &lhs,      \
0365                                                                                 const ScalarType rhs)           \
0366   {                                                                                                             \
0367     Vector3D<Type> result(lhs);                                                                                 \
0368     result INPLACE rhs;                                                                                         \
0369     return result;                                                                                              \
0370   }                                                                                                             \
0371   template <typename Type, typename ScalarType>                                                                 \
0372   VECGEOM_FORCE_INLINE VECCORE_ATT_HOST_DEVICE Vector3D<Type> operator OPERATOR(const ScalarType lhs,           \
0373                                                                                 Vector3D<Type> const &rhs)      \
0374   {                                                                                                             \
0375     Vector3D<Type> result(lhs);                                                                                 \
0376     result INPLACE rhs;                                                                                         \
0377     return result;                                                                                              \
0378   }
0379 VECTOR3D_BINARY_OP(+, +=)
0380 VECTOR3D_BINARY_OP(-, -=)
0381 VECTOR3D_BINARY_OP(*, *=)
0382 VECTOR3D_BINARY_OP(/, /=)
0383 #undef VECTOR3D_BINARY_OP
0384
0385 VECGEOM_FORCE_INLINE
0386 VECCORE_ATT_HOST_DEVICE
0387 bool operator==(Vector3D<Precision> const &lhs, Vector3D<Precision> const &rhs)
0388 {
0389   return Abs(lhs[0] - rhs[0]) < kTolerance && Abs(lhs[1] - rhs[1]) < kTolerance && Abs(lhs[2] - rhs[2]) < kTolerance;
0390 }
0391
0392 VECGEOM_FORCE_INLINE
0393 VECCORE_ATT_HOST_DEVICE
0394 bool operator!=(Vector3D<Precision> const &lhs, Vector3D<Precision> const &rhs)
0395 {
0396   return !(lhs == rhs);
0397 }
0398
0399 template <typename Type>
0400 VECGEOM_FORCE_INLINE
0401 VECCORE_ATT_HOST_DEVICE
0402 Vector3D<Type> operator-(Vector3D<Type> const &vec)
0403 {
0404   return Vector3D<Type>(-vec[0], -vec[1], -vec[2]);
0405 }
0406
0407 VECCORE_ATT_HOST_DEVICE
0408 VECGEOM_FORCE_INLINE
0409 Vector3D<bool> operator!(Vector3D<bool> const &vec)
0410 {
0411   return Vector3D<bool>(!vec[0], !vec[1], !vec[2]);
0412 }
0413
0414 #pragma GCC diagnostic push
0415 #pragma GCC diagnostic ignored "-Weffc++"
0416 #define VECTOR3D_SCALAR_BOOLEAN_LOGICAL_OP(OPERATOR)                                               \
0417   VECCORE_ATT_HOST_DEVICE                                                                          \
0418   VECGEOM_FORCE_INLINE                                                                             \
0419   Vector3D<bool> operator OPERATOR(Vector3D<bool> const &lhs, Vector3D<bool> const &rhs)           \
0420   {                                                                                                \
0421     return Vector3D<bool>(lhs[0] OPERATOR rhs[0], lhs[1] OPERATOR rhs[1], lhs[2] OPERATOR rhs[2]); \
0422   }
0423 VECTOR3D_SCALAR_BOOLEAN_LOGICAL_OP(&&)
0424 VECTOR3D_SCALAR_BOOLEAN_LOGICAL_OP(||)
0425 #undef VECTOR3D_SCALAR_BOOLEAN_LOGICAL_OP
0426 #pragma GCC diagnostic pop
0427 } // namespace VECGEOM_IMPL_NAMESPACE
0428 } // namespace vecgeom
0429
0430 namespace vecCore {
0431
0432 template <typename T>
0433 VECGEOM_FORCE_INLINE
0434 VECCORE_ATT_HOST_DEVICE
0435 void MaskedAssign(vecgeom::Vector3D<T> &v, const vecCore::Mask<T> &mask, const vecgeom::Vector3D<T> &val)
0436 {
0437   vecCore::MaskedAssign(v[0], mask, val[0]);
0438   vecCore::MaskedAssign(v[1], mask, val[1]);
0439   vecCore::MaskedAssign(v[2], mask, val[2]);
0440 }
0441 } // namespace vecCore
0442
0443 // for use in GEANT4
0444 using UVector3 = VECGEOM_NAMESPACE::Vector3D<double>;
0445
0446 #endif // VECGEOM_BASE_VECTOR3D_H_