volumes/kernel/ExtrudedImplementation.h

0001 //===-- kernel/ExtrudedImplementation.h ----------------------------------*- C++ -*-===//
0002 //===--------------------------------------------------------------------------===//
0003 /// @file ExtrudedImplementation.h
0004 /// @author mihaela.gheata@cern.ch
0005
0006 #ifndef VECGEOM_VOLUMES_KERNEL_EXTRUDEDIMPLEMENTATION_H_
0007 #define VECGEOM_VOLUMES_KERNEL_EXTRUDEDIMPLEMENTATION_H_
0008
0009 #include <cstdio>
0010 #include <VecCore/VecCore>
0011 #include "VecGeom/base/Config.h"
0012 #include "VecGeom/volumes/kernel/GenericKernels.h"
0013 #include "VecGeom/base/Vector3D.h"
0014
0015 #include "TessellatedImplementation.h"
0016 #include "SExtruImplementation.h"
0017
0018 namespace vecgeom {
0019
0020 VECGEOM_DEVICE_FORWARD_DECLARE(struct ExtrudedImplementation;);
0021 VECGEOM_DEVICE_DECLARE_CONV(struct, ExtrudedImplementation);
0022
0023 inline namespace VECGEOM_IMPL_NAMESPACE {
0024
0025 class PlacedExtruded;
0026 class ExtrudedStruct;
0027 class UnplacedExtruded;
0028
0029 struct ExtrudedImplementation {
0030
0031   using PlacedShape_t    = PlacedExtruded;
0032   using UnplacedStruct_t = ExtrudedStruct;
0033   using UnplacedVolume_t = UnplacedExtruded;
0034
0035   template <typename Real_v, typename Bool_v>
0036   VECGEOM_FORCE_INLINE VECCORE_ATT_HOST_DEVICE static void Contains(UnplacedStruct_t const &extruded,
0037                                                                     Vector3D<Real_v> const &point, Bool_v &inside)
0038   {
0039     inside = false;
0040     if (extruded.fIsSxtru) {
0041       SExtruImplementation::Contains<Real_v, Bool_v>(extruded.fSxtruHelper, point, inside);
0042       return;
0043     }
0044
0045 #ifndef VECGEOM_ENABLE_CUDA
0046     if (extruded.fUseTslSections) {
0047       // Find the Z section
0048       int zIndex = extruded.FindZSegment(point[2]);
0049       if ((zIndex < 0) || (zIndex >= (int)extruded.GetNSegments())) return;
0050       inside = extruded.fTslSections[zIndex]->Contains(point);
0051       return;
0052     }
0053 #endif
0054     TessellatedImplementation::Contains<Real_v, Bool_v>(extruded.fTslHelper, point, inside);
0055   }
0056
0057   template <typename Real_v, typename Inside_v>
0058   VECGEOM_FORCE_INLINE VECCORE_ATT_HOST_DEVICE static void Inside(UnplacedStruct_t const &extruded,
0059                                                                   Vector3D<Real_v> const &point, Inside_v &inside)
0060   {
0061     inside = EInside::kOutside;
0062
0063     if (extruded.fIsSxtru) {
0064       SExtruImplementation::Inside<Real_v, Inside_v>(extruded.fSxtruHelper, point, inside);
0065       return;
0066     }
0067
0068 #ifndef VECGEOM_ENABLE_CUDA
0069     if (extruded.fUseTslSections) {
0070       const int nseg = (int)extruded.GetNSegments();
0071       int zIndex     = extruded.FindZSegment(point[2]);
0072       if ((zIndex < 0) || (zIndex > nseg)) return;
0073       inside = extruded.fTslSections[Min(zIndex, nseg - 1)]->Inside(point);
0074       if (inside == EInside::kOutside) return;
0075       if (inside == EInside::kInside) {
0076         // Need to check if point on Z section
0077         if (((zIndex == 0) || (zIndex == nseg)) &&
0078             vecCore::math::Abs(point[2] - extruded.fZPlanes[zIndex]) < kTolerance) {
0079           inside = EInside::kSurface;
0080         }
0081       } else {
0082         inside = EInside::kSurface;
0083       }
0084       return;
0085     }
0086 #endif
0087     TessellatedImplementation::Inside<Real_v, Inside_v>(extruded.fTslHelper, point, inside);
0088   }
0089
0090   template <typename Real_v>
0091   VECGEOM_FORCE_INLINE VECCORE_ATT_HOST_DEVICE static void DistanceToIn(UnplacedStruct_t const &extruded,
0092                                                                         Vector3D<Real_v> const &point,
0093                                                                         Vector3D<Real_v> const &direction,
0094                                                                         Real_v const &stepMax, Real_v &distance)
0095   {
0096 // Note that Real_v is always double here
0097 #ifdef EFFICIENT_TSL_DISTANCETOIN
0098     if (extruded.fUseTslSections) {
0099       // Check if the bounding box is hit
0100       const Vector3D<Real_v> invdir(Real_v(1.0) / NonZero(direction.x()), Real_v(1.0) / NonZero(direction.y()),
0101                                     Real_v(1.0) / NonZero(direction.z()));
0102       Vector3D<int> sign;
0103       sign[0]  = invdir.x() < 0;
0104       sign[1]  = invdir.y() < 0;
0105       sign[2]  = invdir.z() < 0;
0106       distance = BoxImplementation::IntersectCachedKernel2<Real_v, Real_v>(&extruded.fTslHelper.fMinExtent, point,
0107                                                                            invdir, sign.x(), sign.y(), sign.z(),
0108                                                                            -kTolerance, InfinityLength<Real_v>());
0109       if (distance >= stepMax) return;
0110
0111       // Perform explicit Inside check to detect wrong side points. This impacts
0112       // DistanceToIn performance by about 5% for all topologies
0113       // auto inside = ScalarInsideKernel(unplaced, point);
0114       // if (inside == kInside) return -1.;
0115
0116       int zIndex     = extruded.FindZSegment(point[2]);
0117       const int zMax = extruded.GetNSegments();
0118       // Don't go out of bounds here, as the first/last segment should be checked
0119       // even if the point is outside of Z-bounds
0120       bool fromOutZ =
0121           (point[2] < extruded.fZPlanes[0] + kTolerance) || (point[2] > extruded.fZPlanes[zMax] - kTolerance);
0122       zIndex = zIndex < 0 ? 0 : (zIndex >= zMax ? zMax - 1 : zIndex);
0123
0124       // Traverse Z-segments left or right depending on sign of direction
0125       bool goingRight = direction[2] >= 0;
0126
0127       distance = InfinityLength<Real_v>();
0128       if (goingRight) {
0129         for (int zSegCount = zMax; zIndex < zSegCount; ++zIndex) {
0130           bool skipZ = fromOutZ && (zSegCount == 0);
0131           if (skipZ)
0132             distance = extruded.fTslSections[zIndex]->DistanceToIn<true>(point, direction, invdir.z(), stepMax);
0133           else
0134             distance = extruded.fTslSections[zIndex]->DistanceToIn<false>(point, direction, invdir.z(), stepMax);
0135           // No segment further away can be at a shorter distance to the point, so
0136           // if a valid distance is found, only endcaps remain to be investigated
0137           if (distance >= -kTolerance && distance < InfinityLength<Precision>()) break;
0138         }
0139       } else {
0140         // Going left
0141         for (; zIndex >= 0; --zIndex) {
0142           bool skipZ = fromOutZ && (zIndex == zMax);
0143           if (skipZ)
0144             distance = extruded.fTslSections[zIndex]->DistanceToIn<true>(point, direction, invdir.z(), stepMax);
0145           else
0146             distance = extruded.fTslSections[zIndex]->DistanceToIn<false>(point, direction, invdir.z(), stepMax);
0147           // No segment further away can be at a shorter distance to the point, so
0148           // if a valid distance is found, only endcaps remain to be investigated
0149           if (distance >= -kTolerance && distance < InfinityLength<Precision>()) break;
0150         }
0151       }
0152     }
0153 #endif
0154     if (extruded.fIsSxtru)
0155       SExtruImplementation::DistanceToIn<Real_v>(extruded.fSxtruHelper, point, direction, stepMax, distance);
0156     else
0157       TessellatedImplementation::DistanceToIn<Real_v>(extruded.fTslHelper, point, direction, stepMax, distance);
0158   }
0159
0160   template <typename Real_v>
0161   VECGEOM_FORCE_INLINE VECCORE_ATT_HOST_DEVICE static void DistanceToOut(UnplacedStruct_t const &extruded,
0162                                                                          Vector3D<Real_v> const &point,
0163                                                                          Vector3D<Real_v> const &direction,
0164                                                                          Real_v const &stepMax, Real_v &distance)
0165   {
0166     if (extruded.fIsSxtru)
0167       SExtruImplementation::DistanceToOut<Real_v>(extruded.fSxtruHelper, point, direction, stepMax, distance);
0168     else
0169       TessellatedImplementation::DistanceToOut<Real_v>(extruded.fTslHelper, point, direction, stepMax, distance);
0170   }
0171
0172   template <typename Real_v>
0173   VECGEOM_FORCE_INLINE VECCORE_ATT_HOST_DEVICE static void SafetyToIn(UnplacedStruct_t const &extruded,
0174                                                                       Vector3D<Real_v> const &point, Real_v &safety)
0175   {
0176     if (extruded.fIsSxtru)
0177       SExtruImplementation::SafetyToIn<Real_v>(extruded.fSxtruHelper, point, safety);
0178     else
0179       TessellatedImplementation::SafetyToIn<Real_v>(extruded.fTslHelper, point, safety);
0180   }
0181
0182   template <typename Real_v>
0183   VECGEOM_FORCE_INLINE VECCORE_ATT_HOST_DEVICE static void SafetyToOut(UnplacedStruct_t const &extruded,
0184                                                                        Vector3D<Real_v> const &point, Real_v &safety)
0185   {
0186     if (extruded.fIsSxtru)
0187       SExtruImplementation::SafetyToOut<Real_v>(extruded.fSxtruHelper, point, safety);
0188     else
0189       TessellatedImplementation::SafetyToOut<Real_v>(extruded.fTslHelper, point, safety);
0190   }
0191
0192   template <typename Real_v>
0193   VECGEOM_FORCE_INLINE VECCORE_ATT_HOST_DEVICE static Vector3D<Real_v> NormalKernel(
0194       UnplacedStruct_t const &extruded, Vector3D<Real_v> const &point, typename vecCore::Mask_v<Real_v> &valid)
0195   {
0196     // Computes the normal on a surface and returns it as a unit vector
0197     if (extruded.fIsSxtru) return SExtruImplementation::NormalKernel<Real_v>(extruded.fSxtruHelper, point, valid);
0198     return TessellatedImplementation::NormalKernel<Real_v>(extruded.fTslHelper, point, valid);
0199   }
0200
0201 }; // end ExtrudedImplementation
0202
0203 } // namespace VECGEOM_IMPL_NAMESPACE
0204 } // namespace vecgeom
0205
0206 #endif // VECGEOM_VOLUMES_KERNEL_EXTRUDEDIMPLEMENTATION_H_