volumes/kernel/GenericKernels.h

0001 /// \file GenericKernels.h
0002 /// \author Johannes de Fine Licht (johannes.definelicht@cern.ch)
0003
0004 #ifndef VECGEOM_VOLUMES_KERNEL_GENERICKERNELS_H_
0005 #define VECGEOM_VOLUMES_KERNEL_GENERICKERNELS_H_
0006
0007 #include "VecGeom/base/Global.h"
0008 #include "VecGeom/base/Transformation3D.h"
0009 #include "VecGeom/base/Vector3D.h"
0010
0011 namespace vecgeom {
0012 inline namespace VECGEOM_IMPL_NAMESPACE {
0013
0014 template <class Backend>
0015 struct GenericKernels {
0016
0017   typedef typename Backend::precision_v Float_t;
0018   typedef typename Backend::int_v Int_t;
0019   typedef typename Backend::bool_v Bool_t;
0020
0021 }; // End struct GenericKernels
0022
0023 template <bool tolerant, typename T>
0024 VECGEOM_FORCE_INLINE
0025 VECCORE_ATT_HOST_DEVICE
0026 T MakePlusTolerant(T const &x, vecCore::Scalar<T> halftol = kHalfTolerance)
0027 {
0028   return (tolerant) ? x + halftol : x;
0029 }
0030
0031 template <bool tolerant, typename T>
0032 VECGEOM_FORCE_INLINE
0033 VECCORE_ATT_HOST_DEVICE
0034 T MakeMinusTolerant(T const &x, vecCore::Scalar<T> halftol = kHalfTolerance)
0035 {
0036   return (tolerant) ? x - T(halftol) : x;
0037 }
0038
0039 /// @brief Utilities to compute tolerance value for cross products. Length should be an overestimate
0040 /// of the point vector length
0041 /// P = point to check if on one side or the other of the AB segment. Cross product computed as:
0042 ///   cross = AP x AB
0043 /// The distance from point P to segment AB is cross/|AB|, hence the tolerance of cross is kTolerance * |AB|
0044 /// One needs to pass length > |AB| to the method.
0045 template <bool tolerant, typename T>
0046 VECGEOM_FORCE_INLINE
0047 VECCORE_ATT_HOST_DEVICE
0048 T MakePlusTolerantCrossProduct(T const &x, T const &length, vecCore::Scalar<T> halftol = kHalfTolerance)
0049 {
0050   return (tolerant) ? x + length * T(halftol) : x;
0051 }
0052
0053 template <bool tolerant, typename T>
0054 VECGEOM_FORCE_INLINE
0055 VECCORE_ATT_HOST_DEVICE
0056 T MakeMinusTolerantCrossProduct(T const &x,  T const &length, vecCore::Scalar<T> halftol = kHalfTolerance)
0057 {
0058   return (tolerant) ? x - length * T(halftol) : x;
0059 }
0060
0061 /// @brief Utility to compute (x + tol)^2 for proper account of tolerances when comparing squares.
0062 template <bool tolerant, typename T>
0063 VECGEOM_FORCE_INLINE
0064 VECCORE_ATT_HOST_DEVICE
0065 T MakePlusTolerantSquare(T const &x, vecCore::Scalar<T> tol = kTolerance)
0066 {
0067   // calculate (x + halftol) * (x + halftol) which should always >= 0;
0068   // in order to be fast, we neglect the + tol * tol term (since it should be negligible)
0069   return (tolerant) ? x * (x + T(2.0 * tol)) : x * x;
0070 }
0071
0072 /// @brief Utility to compute (x - tol)^2 for proper account of tolerances when comparing squares.
0073 template <bool tolerant, typename T>
0074 VECGEOM_FORCE_INLINE
0075 VECCORE_ATT_HOST_DEVICE
0076 T MakeMinusTolerantSquare(T const &x, vecCore::Scalar<T> tol = kTolerance)
0077 {
0078   // calculate (x - halftol) * (x - halftol) which should always >= 0;
0079   // in order to be fast, we neglect the + tol * tol term (since it should be negligible)
0080   // but we make sure that there is never a negative sign (hence the Abs)
0081   return (tolerant) ? Abs(x * (x - T(2.0 * tol))) : x * x;
0082 }
0083
0084 template <bool treatSurfaceT, class Backend>
0085 struct TreatSurfaceTraits;
0086 template <class Backend>
0087 struct TreatSurfaceTraits<true, Backend> {
0088   typedef typename Backend::inside_v Surface_t;
0089   static const Inside_t kInside  = 0;
0090   static const Inside_t kOutside = 2;
0091 };
0092 template <class Backend>
0093 struct TreatSurfaceTraits<false, Backend> {
0094   typedef typename Backend::bool_v Surface_t;
0095   static const bool kInside  = true;
0096   static const bool kOutside = false;
0097 };
0098
0099 /// \brief Flips the sign of an input value depending on the set template
0100 ///        parameter.
0101 template <bool flipT>
0102 struct Flip;
0103
0104 template <>
0105 struct Flip<true> {
0106   template <class T>
0107   VECGEOM_FORCE_INLINE
0108   VECCORE_ATT_HOST_DEVICE
0109   static T FlipSign(T const &value)
0110   {
0111     return -value;
0112   }
0113   template <class T>
0114   VECGEOM_FORCE_INLINE
0115   VECCORE_ATT_HOST_DEVICE
0116   static T FlipLogical(T const &value)
0117   {
0118     return !value;
0119   }
0120 };
0121
0122 template <>
0123 struct Flip<false> {
0124   template <class T>
0125   VECGEOM_FORCE_INLINE
0126   VECCORE_ATT_HOST_DEVICE
0127   static T FlipSign(T const &value)
0128   {
0129     return value;
0130   }
0131   template <class T>
0132   VECGEOM_FORCE_INLINE
0133   VECCORE_ATT_HOST_DEVICE
0134   static T FlipLogical(T const &value)
0135   {
0136     return value;
0137   }
0138 };
0139
0140 template <class Backend>
0141 VECGEOM_FORCE_INLINE
0142 VECCORE_ATT_HOST_DEVICE
0143 typename Backend::precision_v NormalizeAngle(typename Backend::precision_v a)
0144 {
0145   return a + kTwoPi * ((a < 0) - typename Backend::int_v(a / kTwoPi));
0146 }
0147
0148 // \param corner0 First corner of line segment.
0149 // \param corner1 Second corner of line segment.
0150 // \return Shortest distance from the point to the three dimensional line
0151 //         segment represented by the two input points.
0152 template <class Backend>
0153 VECGEOM_FORCE_INLINE
0154 VECCORE_ATT_HOST_DEVICE
0155 typename Backend::precision_v DistanceToLineSegmentSquared(Vector3D<Precision> corner0, Vector3D<Precision> corner1,
0156                                                            Vector3D<typename Backend::precision_v> const &point)
0157 {
0158
0159   typedef typename Backend::precision_v Float_t;
0160   typedef typename Backend::bool_v Bool_t;
0161
0162   Float_t result(kInfLength);
0163
0164   // Shortest distance is to corner of segment
0165   Vector3D<Precision> line     = corner1 - corner0;
0166   Vector3D<Float_t> projection = point - corner0;
0167   Float_t dot0                 = projection.Dot(line);
0168   Bool_t condition             = dot0 <= 0;
0169   vecCore__MaskedAssignFunc(result, condition, (point - corner0).Mag2());
0170   if (vecCore::MaskFull(condition)) return result;
0171   Precision dot1 = line.Mag2();
0172   condition      = dot1 <= dot0;
0173   vecCore__MaskedAssignFunc(result, condition, (point - corner1).Mag2());
0174   condition = result < kInfLength;
0175   if (vecCore::MaskFull(condition)) return result;
0176
0177   // Shortest distance is to point on segment
0178   vecCore__MaskedAssignFunc(result, !condition, ((corner0 + (dot0 / dot1) * line) - point).Mag2());
0179
0180   return result;
0181 }
0182
0183 // \param corner0 First corner of line segment.
0184 // \param corner1 Second corner of line segment.
0185 // \return Shortest distance from the point to the three dimensional line
0186 //         segment represented by the two input points.
0187 template <typename Real_v>
0188 VECGEOM_FORCE_INLINE
0189 VECCORE_ATT_HOST_DEVICE
0190 Real_v DistanceToLineSegmentSquared1(Vector3D<Precision> corner0, Vector3D<Precision> corner1,
0191                                      Vector3D<Real_v> const &point)
0192 {
0193
0194   using Bool_v = vecCore::Mask_v<Real_v>;
0195
0196   Real_v result = InfinityLength<Real_v>();
0197
0198   // Shortest distance is to corner of segment
0199   Vector3D<Precision> line    = corner1 - corner0;
0200   Vector3D<Real_v> projection = point - corner0;
0201   Real_v dot0                 = projection.Dot(line);
0202   Bool_v condition            = dot0 <= 0;
0203   vecCore__MaskedAssignFunc(result, condition, (point - corner0).Mag2());
0204   if (vecCore::MaskFull(condition)) return result;
0205   Precision dot1 = line.Mag2();
0206   condition      = dot1 <= dot0;
0207   vecCore__MaskedAssignFunc(result, condition, (point - corner1).Mag2());
0208   condition = result < kInfLength;
0209   if (vecCore::MaskFull(condition)) return result;
0210
0211   // Shortest distance is to point on segment
0212   vecCore__MaskedAssignFunc(result, !condition, Real_v(((corner0 + (dot0 / dot1) * line) - point).Mag2()));
0213
0214   return result;
0215 }
0216
0217 // \param corner0 First corners of line segments.
0218 // \param corner1 Second corners of line segments.
0219 // \return Shortest distance from the point to the three dimensional set of line
0220 //         segments represented by the input corners.
0221 template <typename Real_v>
0222 VECGEOM_FORCE_INLINE
0223 VECCORE_ATT_HOST_DEVICE
0224 Real_v DistanceToLineSegmentSquared2(Vector3D<Real_v> const &corner0, Vector3D<Real_v> const &corner1,
0225                                      Vector3D<Real_v> const &point, vecCore::Mask_v<Real_v> const &mask)
0226 {
0227
0228   using Bool_v = vecCore::Mask_v<Real_v>;
0229
0230   Real_v result = InfinityLength<Real_v>();
0231
0232   // Shortest distance is to corner of segment
0233   Vector3D<Real_v> line       = corner1 - corner0;
0234   Vector3D<Real_v> projection = point - corner0;
0235   Real_v dot0                 = projection.Dot(line);
0236   Bool_v condition            = dot0 <= 0 && mask;
0237   vecCore__MaskedAssignFunc(result, condition, (point - corner0).Mag2());
0238   if (vecCore::MaskFull(condition && mask)) return result;
0239   Real_v dot1 = line.Mag2();
0240   condition   = dot1 <= dot0 && mask;
0241   vecCore__MaskedAssignFunc(result, condition, (point - corner1).Mag2());
0242   condition = result < kInfLength;
0243   if (vecCore::MaskFull(condition && mask)) return result;
0244
0245   // Shortest distance is to point on segment
0246   vecCore__MaskedAssignFunc(result, !condition && mask, Real_v(((corner0 + (dot0 / dot1) * line) - point).Mag2()));
0247
0248   return result;
0249 }
0250
0251 } // End inline namespace
0252
0253 } // End global namespace
0254
0255 #endif // VECGEOM_VOLUMES_KERNEL_GENERICKERNELS_H_