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   VECCORE_ATT_HOST_DEVICE
0036   static void PrintType()
0037   {
0038     //  printf("SpecializedBox<%i, %i>", transCodeT, rotCodeT);
0039   }
0040
0041   template <typename Stream>
0042   static void PrintType(Stream &st, int transCodeT = translation::kGeneric, int rotCodeT = rotation::kGeneric)
0043   {
0044     st << "SpecializedExtruded<" << transCodeT << "," << rotCodeT << ">";
0045   }
0046
0047   template <typename Stream>
0048   static void PrintImplementationType(Stream &st)
0049   {
0050     (void)st;
0051   }
0052
0053   template <typename Stream>
0054   static void PrintUnplacedType(Stream &st)
0055   {
0056     (void)st;
0057   }
0058
0059   template <typename Real_v, typename Bool_v>
0060   VECGEOM_FORCE_INLINE
0061   VECCORE_ATT_HOST_DEVICE
0062   static void Contains(UnplacedStruct_t const &extruded, Vector3D<Real_v> const &point, Bool_v &inside)
0063   {
0064     inside = false;
0065     if (extruded.fIsSxtru) {
0066       SExtruImplementation::Contains<Real_v, Bool_v>(extruded.fSxtruHelper, point, inside);
0067       return;
0068     }
0069
0070 #ifndef VECGEOM_ENABLE_CUDA
0071     if (extruded.fUseTslSections) {
0072       // Find the Z section
0073       int zIndex = extruded.FindZSegment(point[2]);
0074       if ((zIndex < 0) || (zIndex >= (int)extruded.GetNSegments())) return;
0075       inside = extruded.fTslSections[zIndex]->Contains(point);
0076       return;
0077     }
0078 #endif
0079     TessellatedImplementation::Contains<Real_v, Bool_v>(extruded.fTslHelper, point, inside);
0080   }
0081
0082   template <typename Real_v, typename Inside_v>
0083   VECGEOM_FORCE_INLINE
0084   VECCORE_ATT_HOST_DEVICE
0085   static void Inside(UnplacedStruct_t const &extruded, Vector3D<Real_v> const &point, Inside_v &inside)
0086   {
0087     inside = EInside::kOutside;
0088
0089     if (extruded.fIsSxtru) {
0090       SExtruImplementation::Inside<Real_v, Inside_v>(extruded.fSxtruHelper, point, inside);
0091       return;
0092     }
0093
0094 #ifndef VECGEOM_ENABLE_CUDA
0095     if (extruded.fUseTslSections) {
0096       const int nseg = (int)extruded.GetNSegments();
0097       int zIndex     = extruded.FindZSegment(point[2]);
0098       if ((zIndex < 0) || (zIndex > nseg)) return;
0099       inside = extruded.fTslSections[Min(zIndex, nseg - 1)]->Inside(point);
0100       if (inside == EInside::kOutside) return;
0101       if (inside == EInside::kInside) {
0102         // Need to check if point on Z section
0103         if (((zIndex == 0) || (zIndex == nseg)) &&
0104             vecCore::math::Abs(point[2] - extruded.fZPlanes[zIndex]) < kTolerance) {
0105           inside = EInside::kSurface;
0106         }
0107       } else {
0108         inside = EInside::kSurface;
0109       }
0110       return;
0111     }
0112 #endif
0113     TessellatedImplementation::Inside<Real_v, Inside_v>(extruded.fTslHelper, point, inside);
0114   }
0115
0116   template <typename Real_v>
0117   VECGEOM_FORCE_INLINE
0118   VECCORE_ATT_HOST_DEVICE
0119   static void DistanceToIn(UnplacedStruct_t const &extruded, Vector3D<Real_v> const &point,
0120                            Vector3D<Real_v> const &direction, Real_v const &stepMax, Real_v &distance)
0121   {
0122 // Note that Real_v is always double here
0123 #ifdef EFFICIENT_TSL_DISTANCETOIN
0124     if (extruded.fUseTslSections) {
0125       // Check if the bounding box is hit
0126       const Vector3D<Real_v> invdir(Real_v(1.0) / NonZero(direction.x()), Real_v(1.0) / NonZero(direction.y()),
0127                                     Real_v(1.0) / NonZero(direction.z()));
0128       Vector3D<int> sign;
0129       sign[0]  = invdir.x() < 0;
0130       sign[1]  = invdir.y() < 0;
0131       sign[2]  = invdir.z() < 0;
0132       distance = BoxImplementation::IntersectCachedKernel2<Real_v, Real_v>(&extruded.fTslHelper.fMinExtent, point,
0133                                                                            invdir, sign.x(), sign.y(), sign.z(),
0134                                                                            -kTolerance, InfinityLength<Real_v>());
0135       if (distance >= stepMax) return;
0136
0137       // Perform explicit Inside check to detect wrong side points. This impacts
0138       // DistanceToIn performance by about 5% for all topologies
0139       // auto inside = ScalarInsideKernel(unplaced, point);
0140       // if (inside == kInside) return -1.;
0141
0142       int zIndex     = extruded.FindZSegment(point[2]);
0143       const int zMax = extruded.GetNSegments();
0144       // Don't go out of bounds here, as the first/last segment should be checked
0145       // even if the point is outside of Z-bounds
0146       bool fromOutZ =
0147           (point[2] < extruded.fZPlanes[0] + kTolerance) || (point[2] > extruded.fZPlanes[zMax] - kTolerance);
0148       zIndex = zIndex < 0 ? 0 : (zIndex >= zMax ? zMax - 1 : zIndex);
0149
0150       // Traverse Z-segments left or right depending on sign of direction
0151       bool goingRight = direction[2] >= 0;
0152
0153       distance = InfinityLength<Real_v>();
0154       if (goingRight) {
0155         for (int zSegCount = zMax; zIndex < zSegCount; ++zIndex) {
0156           bool skipZ = fromOutZ && (zSegCount == 0);
0157           if (skipZ)
0158             distance = extruded.fTslSections[zIndex]->DistanceToIn<true>(point, direction, invdir.z(), stepMax);
0159           else
0160             distance = extruded.fTslSections[zIndex]->DistanceToIn<false>(point, direction, invdir.z(), stepMax);
0161           // No segment further away can be at a shorter distance to the point, so
0162           // if a valid distance is found, only endcaps remain to be investigated
0163           if (distance >= -kTolerance && distance < InfinityLength<Precision>()) break;
0164         }
0165       } else {
0166         // Going left
0167         for (; zIndex >= 0; --zIndex) {
0168           bool skipZ = fromOutZ && (zIndex == zMax);
0169           if (skipZ)
0170             distance = extruded.fTslSections[zIndex]->DistanceToIn<true>(point, direction, invdir.z(), stepMax);
0171           else
0172             distance = extruded.fTslSections[zIndex]->DistanceToIn<false>(point, direction, invdir.z(), stepMax);
0173           // No segment further away can be at a shorter distance to the point, so
0174           // if a valid distance is found, only endcaps remain to be investigated
0175           if (distance >= -kTolerance && distance < InfinityLength<Precision>()) break;
0176         }
0177       }
0178     }
0179 #endif
0180     if (extruded.fIsSxtru)
0181       SExtruImplementation::DistanceToIn<Real_v>(extruded.fSxtruHelper, point, direction, stepMax, distance);
0182     else
0183       TessellatedImplementation::DistanceToIn<Real_v>(extruded.fTslHelper, point, direction, stepMax, distance);
0184   }
0185
0186   template <typename Real_v>
0187   VECGEOM_FORCE_INLINE
0188   VECCORE_ATT_HOST_DEVICE
0189   static void DistanceToOut(UnplacedStruct_t const &extruded, Vector3D<Real_v> const &point,
0190                             Vector3D<Real_v> const &direction, Real_v const &stepMax, Real_v &distance)
0191   {
0192     if (extruded.fIsSxtru)
0193       SExtruImplementation::DistanceToOut<Real_v>(extruded.fSxtruHelper, point, direction, stepMax, distance);
0194     else
0195       TessellatedImplementation::DistanceToOut<Real_v>(extruded.fTslHelper, point, direction, stepMax, distance);
0196   }
0197
0198   template <typename Real_v>
0199   VECGEOM_FORCE_INLINE
0200   VECCORE_ATT_HOST_DEVICE
0201   static void SafetyToIn(UnplacedStruct_t const &extruded, Vector3D<Real_v> const &point, Real_v &safety)
0202   {
0203     if (extruded.fIsSxtru)
0204       SExtruImplementation::SafetyToIn<Real_v>(extruded.fSxtruHelper, point, safety);
0205     else
0206       TessellatedImplementation::SafetyToIn<Real_v>(extruded.fTslHelper, point, safety);
0207   }
0208
0209   template <typename Real_v>
0210   VECGEOM_FORCE_INLINE
0211   VECCORE_ATT_HOST_DEVICE
0212   static void SafetyToOut(UnplacedStruct_t const &extruded, Vector3D<Real_v> const &point, Real_v &safety)
0213   {
0214     if (extruded.fIsSxtru)
0215       SExtruImplementation::SafetyToOut<Real_v>(extruded.fSxtruHelper, point, safety);
0216     else
0217       TessellatedImplementation::SafetyToOut<Real_v>(extruded.fTslHelper, point, safety);
0218   }
0219
0220   template <typename Real_v>
0221   VECGEOM_FORCE_INLINE
0222   VECCORE_ATT_HOST_DEVICE
0223   static Vector3D<Real_v> NormalKernel(UnplacedStruct_t const &extruded, Vector3D<Real_v> const &point,
0224                                        typename vecCore::Mask_v<Real_v> &valid)
0225   {
0226     // Computes the normal on a surface and returns it as a unit vector
0227     if (extruded.fIsSxtru) return SExtruImplementation::NormalKernel<Real_v>(extruded.fSxtruHelper, point, valid);
0228     return TessellatedImplementation::NormalKernel<Real_v>(extruded.fTslHelper, point, valid);
0229   }
0230
0231 }; // end ExtrudedImplementation
0232
0233 } // namespace VECGEOM_IMPL_NAMESPACE
0234 } // namespace vecgeom
0235
0236 #endif // VECGEOM_VOLUMES_KERNEL_EXTRUDEDIMPLEMENTATION_H_