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 // FMA
0024 template <typename T>
0025 VECCORE_ATT_HOST_DEVICE VECGEOM_FORCE_INLINE T fma_hd(T a, T b, T c)
0026 {
0027 #ifdef VECCORE_CUDA_DEVICE_COMPILATION
0028   return __fma_rn(a, b, c);
0029 #else
0030   return std::fma(a, b, c);
0031 #endif
0032 }
0033
0034 template <typename T>
0035 VECGEOM_FORCE_INLINE VECCORE_ATT_HOST_DEVICE static T Dot_fma(Vector3D<T> const &a, Vector3D<T> const &b)
0036 {
0037   return fma_hd(a[0], b[0], fma_hd(a[1], b[1], a[2] * b[2]));
0038 }
0039
0040 template <typename T>
0041 VECGEOM_FORCE_INLINE VECCORE_ATT_HOST_DEVICE static Vector3D<T> Cross_fma(Vector3D<T> const &left,
0042                                                                           Vector3D<T> const &right)
0043 {
0044   return Vector3D<T>(fma_hd(left[1], right[2], -left[2] * right[1]), fma_hd(left[2], right[0], -left[0] * right[2]),
0045                      fma_hd(left[0], right[1], -left[1] * right[0]));
0046 }
0047
0048 // typesafe sign
0049 template <typename Real_t>
0050 VECGEOM_FORCE_INLINE VECCORE_ATT_HOST_DEVICE constexpr int kSign(Real_t x)
0051 {
0052   return (Real_t(0) < x) - (x < Real_t(0));
0053 }
0054
0055 // relative tolerance specializations
0056 template <typename Real_t>
0057 VECGEOM_FORCE_INLINE VECCORE_ATT_HOST_DEVICE constexpr Real_t kRelTolerance(Real_t x,
0058                                                                             Real_t tolerance = kToleranceStrict<Real_t>)
0059 {
0060   return (vecCore::math::Abs(x) + Real_t(1.)) * tolerance;
0061 }
0062
0063 template <typename Real_t>
0064 VECGEOM_FORCE_INLINE VECCORE_ATT_HOST_DEVICE constexpr Real_t kRelTolerance(Vector3D<Real_t> const &p,
0065                                                                             Real_t tolerance = kToleranceStrict<Real_t>)
0066 {
0067   auto max_abs_coord = vecCore::math::Max(vecCore::math::Abs(p[0]), vecCore::math::Abs(p[1]), vecCore::math::Abs(p[2]));
0068   return kRelTolerance<Real_t>(max_abs_coord, tolerance);
0069 }
0070
0071 template <typename T>
0072 VECGEOM_FORCE_INLINE VECCORE_ATT_HOST_DEVICE constexpr T MakePlusTolerantRel(T const x,
0073                                                                              T tolerance = kToleranceStrict<T>)
0074 {
0075   return (x + kRelTolerance<T>(x, tolerance));
0076 }
0077
0078 template <typename T>
0079 VECGEOM_FORCE_INLINE VECCORE_ATT_HOST_DEVICE constexpr T MakeMinusTolerantRel(T const x,
0080                                                                               T tolerance = kToleranceStrict<T>)
0081 {
0082   return (x - kRelTolerance<T>(x, tolerance));
0083 }
0084
0085 template <bool tolerant, typename T>
0086 VECGEOM_FORCE_INLINE VECCORE_ATT_HOST_DEVICE T MakePlusTolerant(T const &x, vecCore::Scalar<T> tol = kToleranceDist<T>)
0087 {
0088   return (tolerant) ? x + tol : x;
0089 }
0090
0091 template <bool tolerant, typename T>
0092 VECGEOM_FORCE_INLINE VECCORE_ATT_HOST_DEVICE T MakeMinusTolerant(T const &x, vecCore::Scalar<T> tol = kToleranceDist<T>)
0093 {
0094   return (tolerant) ? x - T(tol) : x;
0095 }
0096
0097 /// @brief Utilities to compute tolerance value for cross products. Length should be an overestimate
0098 /// of the point vector length
0099 /// P = point to check if on one side or the other of the AB segment. Cross product computed as:
0100 ///   cross = AP x AB
0101 /// The distance from point P to segment AB is cross/|AB|, hence the tolerance of cross is kTolerance * |AB|
0102 /// One needs to pass length > |AB| to the method.
0103 template <bool tolerant, typename T>
0104 VECGEOM_FORCE_INLINE VECCORE_ATT_HOST_DEVICE T MakePlusTolerantCrossProduct(T const &x, T const &length,
0105                                                                             vecCore::Scalar<T> tol = kToleranceDist<T>)
0106 {
0107   return (tolerant) ? x + length * T(tol) : x;
0108 }
0109
0110 template <bool tolerant, typename T>
0111 VECGEOM_FORCE_INLINE VECCORE_ATT_HOST_DEVICE T MakeMinusTolerantCrossProduct(T const &x, T const &length,
0112                                                                              vecCore::Scalar<T> tol = kToleranceDist<T>)
0113 {
0114   return (tolerant) ? x - length * T(tol) : x;
0115 }
0116
0117 /// @brief Utility to compute (x + tol)^2 for proper account of tolerances when comparing squares.
0118 template <bool tolerant, typename T>
0119 VECGEOM_FORCE_INLINE VECCORE_ATT_HOST_DEVICE T MakePlusTolerantSquare(T const &x,
0120                                                                       vecCore::Scalar<T> tol = kToleranceDist<T>)
0121 {
0122   // calculate (x + halftol) * (x + halftol) which should always >= 0;
0123   // in order to be fast, we neglect the + tol * tol term (since it should be negligible)
0124   return (tolerant) ? x * (x + T(2.0 * tol)) : x * x;
0125 }
0126
0127 /// @brief Utility to compute (x - tol)^2 for proper account of tolerances when comparing squares.
0128 template <bool tolerant, typename T>
0129 VECGEOM_FORCE_INLINE VECCORE_ATT_HOST_DEVICE T MakeMinusTolerantSquare(T const &x,
0130                                                                        vecCore::Scalar<T> tol = kToleranceDist<T>)
0131 {
0132   // calculate (x - halftol) * (x - halftol) which should always >= 0;
0133   // in order to be fast, we neglect the + tol * tol term (since it should be negligible)
0134   // but we make sure that there is never a negative sign (hence the Abs)
0135   return (tolerant) ? Abs(x * (x - T(2.0 * tol))) : x * x;
0136 }
0137
0138 template <typename Real_t>
0139 VECGEOM_FORCE_INLINE VECCORE_ATT_HOST_DEVICE Vector3D<Real_t> NonZeroVector(Vector3D<Real_t> const &vect)
0140 {
0141   return Vector3D<Real_t>(NonZero(vect[0]), NonZero(vect[1]), NonZero(vect[2]));
0142 }
0143
0144 template <typename Real_t>
0145 VECGEOM_FORCE_INLINE VECCORE_ATT_HOST_DEVICE Real_t InvdirNearParallel(Real_t dist)
0146 {
0147   return kInvTolerance * dist;
0148 }
0149
0150 template <bool treatSurfaceT, class Backend>
0151 struct TreatSurfaceTraits;
0152 template <class Backend>
0153 struct TreatSurfaceTraits<true, Backend> {
0154   typedef typename Backend::inside_v Surface_t;
0155   static const Inside_t kInside  = 0;
0156   static const Inside_t kOutside = 2;
0157 };
0158 template <class Backend>
0159 struct TreatSurfaceTraits<false, Backend> {
0160   typedef typename Backend::bool_v Surface_t;
0161   static const bool kInside  = true;
0162   static const bool kOutside = false;
0163 };
0164
0165 /// \brief Flips the sign of an input value depending on the set template
0166 ///        parameter.
0167 template <bool flipT>
0168 struct Flip;
0169
0170 template <>
0171 struct Flip<true> {
0172   template <class T>
0173   VECGEOM_FORCE_INLINE VECCORE_ATT_HOST_DEVICE static T FlipSign(T const &value)
0174   {
0175     return -value;
0176   }
0177   template <class T>
0178   VECGEOM_FORCE_INLINE VECCORE_ATT_HOST_DEVICE static T FlipLogical(T const &value)
0179   {
0180     return !value;
0181   }
0182 };
0183
0184 template <>
0185 struct Flip<false> {
0186   template <class T>
0187   VECGEOM_FORCE_INLINE VECCORE_ATT_HOST_DEVICE static T FlipSign(T const &value)
0188   {
0189     return value;
0190   }
0191   template <class T>
0192   VECGEOM_FORCE_INLINE VECCORE_ATT_HOST_DEVICE static T FlipLogical(T const &value)
0193   {
0194     return value;
0195   }
0196 };
0197
0198 template <class Backend>
0199 VECGEOM_FORCE_INLINE VECCORE_ATT_HOST_DEVICE typename Backend::precision_v NormalizeAngle(
0200     typename Backend::precision_v a)
0201 {
0202   return a + kTwoPi * ((a < 0) - typename Backend::int_v(a / kTwoPi));
0203 }
0204
0205 // \param corner0 First corner of line segment.
0206 // \param corner1 Second corner of line segment.
0207 // \return Shortest distance from the point to the three dimensional line
0208 //         segment represented by the two input points.
0209 template <class Backend>
0210 VECGEOM_FORCE_INLINE VECCORE_ATT_HOST_DEVICE typename Backend::precision_v DistanceToLineSegmentSquared(
0211     Vector3D<Precision> corner0, Vector3D<Precision> corner1, Vector3D<typename Backend::precision_v> const &point)
0212 {
0213
0214   typedef typename Backend::precision_v Float_t;
0215   typedef typename Backend::bool_v Bool_t;
0216
0217   Float_t result(kInfLength);
0218
0219   // Shortest distance is to corner of segment
0220   Vector3D<Precision> line     = corner1 - corner0;
0221   Vector3D<Float_t> projection = point - corner0;
0222   Float_t dot0                 = projection.Dot(line);
0223   Bool_t condition             = dot0 <= 0;
0224   vecCore__MaskedAssignFunc(result, condition, (point - corner0).Mag2());
0225   if (vecCore::MaskFull(condition)) return result;
0226   Precision dot1 = line.Mag2();
0227   condition      = dot1 <= dot0;
0228   vecCore__MaskedAssignFunc(result, condition, (point - corner1).Mag2());
0229   condition = result < kInfLength;
0230   if (vecCore::MaskFull(condition)) return result;
0231
0232   // Shortest distance is to point on segment
0233   vecCore__MaskedAssignFunc(result, !condition, ((corner0 + (dot0 / dot1) * line) - point).Mag2());
0234
0235   return result;
0236 }
0237
0238 // \param corner0 First corner of line segment.
0239 // \param corner1 Second corner of line segment.
0240 // \return Shortest distance from the point to the three dimensional line
0241 //         segment represented by the two input points.
0242 template <typename Real_v>
0243 VECGEOM_FORCE_INLINE VECCORE_ATT_HOST_DEVICE Real_v DistanceToLineSegmentSquared1(Vector3D<Precision> corner0,
0244                                                                                   Vector3D<Precision> corner1,
0245                                                                                   Vector3D<Real_v> const &point)
0246 {
0247
0248   using Bool_v = vecCore::Mask_v<Real_v>;
0249
0250   Real_v result = InfinityLength<Real_v>();
0251
0252   // Shortest distance is to corner of segment
0253   Vector3D<Precision> line    = corner1 - corner0;
0254   Vector3D<Real_v> projection = point - corner0;
0255   Real_v dot0                 = projection.Dot(line);
0256   Bool_v condition            = dot0 <= 0;
0257   vecCore__MaskedAssignFunc(result, condition, (point - corner0).Mag2());
0258   if (vecCore::MaskFull(condition)) return result;
0259   Precision dot1 = line.Mag2();
0260   condition      = dot1 <= dot0;
0261   vecCore__MaskedAssignFunc(result, condition, (point - corner1).Mag2());
0262   condition = result < kInfLength;
0263   if (vecCore::MaskFull(condition)) return result;
0264
0265   // Shortest distance is to point on segment
0266   vecCore__MaskedAssignFunc(result, !condition, Real_v(((corner0 + (dot0 / dot1) * line) - point).Mag2()));
0267
0268   return result;
0269 }
0270
0271 // \param corner0 First corners of line segments.
0272 // \param corner1 Second corners of line segments.
0273 // \return Shortest distance from the point to the three dimensional set of line
0274 //         segments represented by the input corners.
0275 template <typename Real_v>
0276 VECGEOM_FORCE_INLINE VECCORE_ATT_HOST_DEVICE Real_v DistanceToLineSegmentSquared2(Vector3D<Real_v> const &corner0,
0277                                                                                   Vector3D<Real_v> const &corner1,
0278                                                                                   Vector3D<Real_v> const &point,
0279                                                                                   vecCore::Mask_v<Real_v> const &mask)
0280 {
0281
0282   using Bool_v = vecCore::Mask_v<Real_v>;
0283
0284   Real_v result = InfinityLength<Real_v>();
0285
0286   // Shortest distance is to corner of segment
0287   Vector3D<Real_v> line       = corner1 - corner0;
0288   Vector3D<Real_v> projection = point - corner0;
0289   Real_v dot0                 = projection.Dot(line);
0290   Bool_v condition            = dot0 <= 0 && mask;
0291   vecCore__MaskedAssignFunc(result, condition, (point - corner0).Mag2());
0292   if (vecCore::MaskFull(condition && mask)) return result;
0293   Real_v dot1 = line.Mag2();
0294   condition   = dot1 <= dot0 && mask;
0295   vecCore__MaskedAssignFunc(result, condition, (point - corner1).Mag2());
0296   condition = result < kInfLength;
0297   if (vecCore::MaskFull(condition && mask)) return result;
0298
0299   // Shortest distance is to point on segment
0300   vecCore__MaskedAssignFunc(result, !condition && mask, Real_v(((corner0 + (dot0 / dot1) * line) - point).Mag2()));
0301
0302   return result;
0303 }
0304
0305 /// @brief Solver for quadratic equations
0306 /// @tparam Real_t Floating point type
0307 /// @param a, b, c Quadratic equation coefficients
0308 /// @param roots Equation roots
0309 /// @param numroots Number of roots grater than -kTolerance
0310 template <typename Real_t>
0311 VECGEOM_FORCE_INLINE VECCORE_ATT_HOST_DEVICE int QuadraticSolver(Real_t a, Real_t b, Real_t c, Real_t roots[2])
0312 {
0313   Real_t inva  = 1. / NonZero(a);
0314   Real_t phalf = 0.5 * b * inva;
0315   Real_t q     = c * inva;
0316   Real_t delta = phalf * phalf - q;
0317   if (delta < Real_t(0)) return 0;
0318   Real_t r1    = -phalf - Sign(phalf) * Sqrt(delta);
0319   Real_t r2    = q / vecgeom::NonZero(r1);
0320   int numroots = int(r1 > -vecgeom::kToleranceStrict<Real_t>) + int(r2 > -vecgeom::kToleranceStrict<Real_t>);
0321   roots[0]     = Min(r1, r2);
0322   roots[1]     = Max(r1, r2);
0323   if (numroots == 1) roots[0] = roots[1];
0324   return numroots;
0325 }
0326
0327 /**
0328  * Robust unit normal of triangle (a,b,c).
0329  * - Computes three cross products: (b-a)x(c-a), (c-b)x(a-b), (a-c)x(b-c)
0330  * - Picks the one with largest squared length (best-conditioned)
0331  * - Normalizes with scale-safe norm; returns false if degenerate
0332  */
0333 template <typename Real_t>
0334 VECCORE_ATT_HOST_DEVICE bool TriangleUnitNormalRobust(const Vector3D<Real_t> &a, const Vector3D<Real_t> &b,
0335                                                       const Vector3D<Real_t> &c, Vector3D<Real_t> &n_unit, Real_t &dist)
0336 {
0337   using V3d = Vector3D<double>;
0338   V3d ab    = V3d(b) - V3d(a);
0339   V3d ac    = V3d(c) - V3d(a);
0340   V3d bc    = V3d(c) - V3d(b);
0341
0342   // Three candidates (mathematically identical, numerically different)
0343   V3d n0 = Cross_fma(ab, ac);   // base a
0344   V3d n1 = Cross_fma(bc, -ab);  // base b
0345   V3d n2 = Cross_fma(-ac, -bc); // base c
0346
0347   auto s0 = Dot_fma(n0, n0);
0348   auto s1 = Dot_fma(n1, n1);
0349   auto s2 = Dot_fma(n2, n2);
0350
0351   // Pick the one with largest |cross|^2
0352   V3d n          = n0;
0353   double s2_best = s0;
0354   if (s1 > s2_best) {
0355     n       = n1;
0356     s2_best = s1;
0357   }
0358   if (s2 > s2_best) {
0359     n       = n2;
0360     s2_best = s2;
0361   }
0362
0363   // Degenerate?
0364   if (s2_best == 0.) {
0365     n_unit.Set(0.);
0366     dist = 0.;
0367     return false;
0368   }
0369
0370   // Scale-safe normalization
0371   double s  = n.Abs().Max();
0372   double rx = n.x() / s, ry = n.y() / s, rz = n.z() / s;
0373   double len = s * vecCore::math::Sqrt(rx * rx + ry * ry + rz * rz); // = ||n|| = 2*area
0374
0375   double inv = 1. / len;
0376   n *= inv;
0377   n_unit.Set(Real_t(n.x()), Real_t(n.y()), Real_t(n.z()));
0378   dist = -Dot_fma(n, V3d(a));
0379   return true;
0380 }
0381
0382 } // namespace VECGEOM_IMPL_NAMESPACE
0383
0384 } // namespace vecgeom
0385
0386 #endif // VECGEOM_VOLUMES_KERNEL_GENERICKERNELS_H_