VecGeom/management/ABBoxImplementation.h

0001 //===-- management/ABBoxImplementation.h -------------------------*- C++ -*-===//
0002 //===----------------------------------------------------------------------===//
0003 /// @file ABBoxImplementation.h
0004 /// @brief Aligned bounding-box helper kernels for navigation acceleration.
0005
0006 #ifndef VECGEOM_MANAGEMENT_ABBOXIMPLEMENTATION_H_
0007 #define VECGEOM_MANAGEMENT_ABBOXIMPLEMENTATION_H_
0008
0009 #include "VecGeom/base/Vector3D.h"
0010 #include "VecGeom/volumes/kernel/GenericKernels.h"
0011 #include <VecCore/VecCore>
0012
0013 namespace vecgeom {
0014 inline namespace VECGEOM_IMPL_NAMESPACE {
0015
0016 /// @brief Aligned bounding-box helper kernels used by acceleration structures.
0017 /// @details These methods intentionally keep vector mask support because they
0018 /// operate on bounding-box batches in navigation and voxel managers, unlike the
0019 /// scalar `UnplacedBox` navigation kernels.
0020 struct ABBoxImplementation {
0021
0022   /// @brief Intersect one ray with one scalar or lane-vector group of aligned boxes.
0023   /// @details The caller caches inverse direction and sign bits once per ray;
0024   /// vector callers use this to prune ABBox batches without touching the
0025   /// ordinary `UnplacedBox` navigation kernels.
0026   template <typename Real_v, typename basep>
0027   VECGEOM_FORCE_INLINE VECCORE_ATT_HOST_DEVICE static Real_v IntersectCachedKernel2(Vector3D<Real_v> const *corners,
0028                                                                                     Vector3D<basep> const &point,
0029                                                                                     Vector3D<basep> const &inverseray,
0030                                                                                     int signx, int signy, int signz,
0031                                                                                     basep t0, basep t1)
0032   {
0033
0034     using Bool_v = vecCore::Mask<Real_v>;
0035
0036     Real_v tmin  = (corners[signx].x() - Real_v(point.x())) * inverseray.x();
0037     Real_v tymax = (corners[1 - signy].y() - Real_v(point.y())) * inverseray.y();
0038     Bool_v done  = tmin > tymax;
0039     if (vecCore::MaskFull(done)) return InfinityLength<Real_v>();
0040
0041     Real_v tmax  = (corners[1 - signx].x() - Real_v(point.x())) * inverseray.x();
0042     Real_v tymin = (corners[signy].y() - Real_v(point.y())) * inverseray.y();
0043
0044     done |= (tymin > tmax);
0045     if (vecCore::MaskFull(done)) return InfinityLength<Real_v>();
0046
0047     tmin = Max(tmin, tymin);
0048     tmax = Min(tmax, tymax);
0049
0050     Real_v tzmin = (corners[signz].z() - point.z()) * inverseray.z();
0051     Real_v tzmax = (corners[1 - signz].z() - point.z()) * inverseray.z();
0052
0053     done |= (Real_v(tmin) > Real_v(tzmax)) || (Real_v(tzmin) > Real_v(tmax));
0054     if (vecCore::MaskFull(done)) return InfinityLength<Real_v>();
0055
0056     tmin = Max(tmin, tzmin);
0057     tmax = Min(tmax, tzmax);
0058
0059     done |= !((tmin <= Real_v(t1 + kTolerance)) && (tmax > Real_v(t0 - kTolerance)));
0060     vecCore__MaskedAssignFunc(tmin, done, InfinityLength<Real_v>());
0061     return tmin;
0062   }
0063
0064   // a contains kernel to be used with aligned bounding boxes
0065   // scalar and vector modes (aka backend) for boxes but only single points
0066   // should be useful to test one point against many bounding boxes
0067   // TODO: check if this can be unified with the normal generic box kernel
0068   template <typename Real_v, typename Bool_v = vecCore::Mask<Real_v>>
0069   VECGEOM_FORCE_INLINE VECCORE_ATT_HOST_DEVICE static void ABBoxContainsKernel(Vector3D<Real_v> const &lowercorner,
0070                                                                                Vector3D<Real_v> const &uppercorner,
0071                                                                                Vector3D<Precision> const &point,
0072                                                                                Bool_v &inside)
0073   {
0074
0075     inside = lowercorner.x() < Real_v(point.x());
0076     inside &= uppercorner.x() > Real_v(point.x());
0077     if (vecCore::MaskEmpty(inside)) return;
0078
0079     inside &= lowercorner.y() < Real_v(point.y());
0080     inside &= uppercorner.y() > Real_v(point.y());
0081     if (vecCore::MaskEmpty(inside)) return;
0082
0083     inside &= lowercorner.z() < Real_v(point.z());
0084     inside &= uppercorner.z() > Real_v(point.z());
0085   }
0086
0087   // playing with a kernel that can do multi-box - single particle; multi-box -- multi-particle, single-box --
0088   // multi-particle
0089   template <typename T1, typename T2, typename Bool_v>
0090   VECGEOM_FORCE_INLINE VECCORE_ATT_HOST_DEVICE static void ABBoxContainsKernelGeneric(Vector3D<T1> const &lowercorner,
0091                                                                                       Vector3D<T1> const &uppercorner,
0092                                                                                       Vector3D<T2> const &point,
0093                                                                                       Bool_v &inside)
0094   {
0095     inside = lowercorner.x() < T1(point.x());
0096     inside &= uppercorner.x() > T1(point.x());
0097     if (vecCore::MaskEmpty(inside)) return;
0098
0099     inside &= lowercorner.y() < T1(point.y());
0100     inside &= uppercorner.y() > T1(point.y());
0101     if (vecCore::MaskEmpty(inside)) return;
0102
0103     inside &= lowercorner.z() < T1(point.z());
0104     inside &= uppercorner.z() > T1(point.z());
0105   }
0106
0107   // safety square for Bounding boxes
0108   // generic kernel treating one track and one or multiple boxes
0109   // in case a point is inside a box a squared value
0110   // is returned but given an overall negative sign
0111   template <typename Real_v, typename Real_s = typename vecCore::TypeTraits<Real_v>::ScalarType>
0112   VECGEOM_FORCE_INLINE VECCORE_ATT_HOST_DEVICE static Real_v ABBoxSafetySqr(Vector3D<Real_v> const &lowercorner,
0113                                                                             Vector3D<Real_v> const &uppercorner,
0114                                                                             Vector3D<Real_s> const &point)
0115   {
0116
0117     using Vector3D_v = Vector3D<Real_v>;
0118     using Bool_v     = vecCore::Mask<Real_v>;
0119
0120     const Vector3D_v kHalf(Real_v(static_cast<Real_s>(0.5)));
0121     const Vector3D_v origin((uppercorner + lowercorner) * kHalf);
0122     const Vector3D_v delta((uppercorner - lowercorner) * kHalf);
0123     // promote scalar point to vector point
0124     const Vector3D_v promotedpoint(Real_v(point.x()), Real_v(point.y()), Real_v(point.z()));
0125
0126     // it would be nicer to have a standalone Abs function taking Vector3D as input
0127     const Vector3D_v safety = ((promotedpoint - origin).Abs()) - delta;
0128     const Bool_v outsidex   = safety.x() > Real_s(0.);
0129     const Bool_v outsidey   = safety.y() > Real_s(0.);
0130     const Bool_v outsidez   = safety.z() > Real_s(0.);
0131
0132     Real_v runningsafetysqr(0.);                  // safety squared from outside
0133     Real_v runningmax(-InfinityLength<Real_v>()); // relevant for safety when we are inside
0134
0135     // loop over dimensions manually unrolled
0136     // treat x dim
0137     {
0138       // this will be much simplified with operator notation
0139       Real_v tmp(0.);
0140       vecCore__MaskedAssignFunc(tmp, outsidex, safety.x() * safety.x());
0141       runningsafetysqr += tmp;
0142       runningmax = Max(runningmax, safety.x());
0143     }
0144
0145     // treat y dim
0146     {
0147       Real_v tmp(0.);
0148       vecCore__MaskedAssignFunc(tmp, outsidey, safety.y() * safety.y());
0149       runningsafetysqr += tmp;
0150       runningmax = Max(runningmax, safety.y());
0151     }
0152
0153     // treat z dim
0154     {
0155       Real_v tmp(0.);
0156       vecCore__MaskedAssignFunc(tmp, outsidez, safety.z() * safety.z());
0157       runningsafetysqr += tmp;
0158       runningmax = Max(runningmax, safety.z());
0159     }
0160
0161     Bool_v inside = !(outsidex || outsidey || outsidez);
0162     if (!vecCore::MaskEmpty(inside)) vecCore__MaskedAssignFunc(runningsafetysqr, inside, -runningmax * runningmax);
0163     return runningsafetysqr;
0164   }
0165
0166   // safety square for Bounding boxes, returning the squared range for any point in the box
0167   // generic kernel treating one track and one or multiple boxes
0168   // in case a point is inside a box a squared value
0169   // is returned but given an overall negative sign
0170   template <typename Real_v, typename Real_s = typename vecCore::TypeTraits<Real_v>::ScalarType>
0171   VECGEOM_FORCE_INLINE VECCORE_ATT_HOST_DEVICE static Real_v ABBoxSafetyRangeSqr(Vector3D<Real_v> const &lowercorner,
0172                                                                                  Vector3D<Real_v> const &uppercorner,
0173                                                                                  Vector3D<Real_s> const &point,
0174                                                                                  Real_v &safetymaxsqr)
0175   {
0176
0177     using Vector3D_v = Vector3D<Real_v>;
0178     using Bool_v     = vecCore::Mask<Real_v>;
0179
0180     const Vector3D_v kHalf(Real_v(static_cast<Real_s>(0.5)));
0181     const Vector3D_v origin((uppercorner + lowercorner) * kHalf);
0182     const Vector3D_v delta((uppercorner - lowercorner) * kHalf);
0183     // promote scalar point to vector point
0184     const Vector3D_v promotedpoint(Real_v(point.x()), Real_v(point.y()), Real_v(point.z()));
0185
0186     // it would be nicer to have a standalone Abs function taking Vector3D as input
0187     const Vector3D_v safety  = ((promotedpoint - origin).Abs()) - delta;
0188     const Vector3D_v safetyp = ((promotedpoint - origin).Abs()) + delta;
0189     const Bool_v outsidex    = safety.x() > Real_s(0.);
0190     const Bool_v outsidey    = safety.y() > Real_s(0.);
0191     const Bool_v outsidez    = safety.z() > Real_s(0.);
0192
0193     Real_v runningsafetysqr(0.);                  // safety squared from outside
0194     safetymaxsqr = safetyp.Mag2();                // safetymax squared from outside
0195     Real_v runningmax(-InfinityLength<Real_v>()); // relevant for safety when we are inside
0196
0197     // loop over dimensions manually unrolled
0198     // treat x dim
0199     {
0200       // this will be much simplified with operator notation
0201       Real_v tmp(0.);
0202       vecCore__MaskedAssignFunc(tmp, outsidex, safety.x() * safety.x());
0203       runningsafetysqr += tmp;
0204       runningmax = Max(runningmax, safety.x());
0205       //      vecCore__MaskedAssignFunc(tmp, outsidex, safetyp.x() * safetyp.x());
0206       //      safetymaxsqr += tmp;
0207     }
0208
0209     // treat y dim
0210     {
0211       Real_v tmp(0.);
0212       vecCore__MaskedAssignFunc(tmp, outsidey, safety.y() * safety.y());
0213       runningsafetysqr += tmp;
0214       runningmax = Max(runningmax, safety.y());
0215       //      vecCore__MaskedAssignFunc(tmp, outsidey, safetyp.y() * safetyp.y());
0216       //      safetymaxsqr += tmp;
0217     }
0218
0219     // treat z dim
0220     {
0221       Real_v tmp(0.);
0222       vecCore__MaskedAssignFunc(tmp, outsidez, safety.z() * safety.z());
0223       runningsafetysqr += tmp;
0224       runningmax = Max(runningmax, safety.z());
0225       //      vecCore__MaskedAssignFunc(tmp, outsidez, safetyp.z() * safetyp.z());
0226       //      safetymaxsqr += tmp;
0227     }
0228
0229     Bool_v inside = !(outsidex || outsidey || outsidez);
0230     if (!vecCore::MaskEmpty(inside)) vecCore__MaskedAssignFunc(runningsafetysqr, inside, -runningmax * runningmax);
0231     return runningsafetysqr;
0232   }
0233
0234 }; // end aligned bounding box struct
0235
0236 } // namespace VECGEOM_IMPL_NAMESPACE
0237 } // namespace vecgeom
0238
0239 #endif // VECGEOM_MANAGEMENT_ABBOXIMPLEMENTATION_H_