volumes/kernel/PolyhedronImplementation.h

0001 /// \file PolyhedronImplementation.h
0002 /// \author Johannes de Fine Licht (johannes.definelicht@cern.ch)
0003
0004 #ifndef VECGEOM_VOLUMES_KERNEL_POLYHEDRONIMPLEMENTATION_H_
0005 #define VECGEOM_VOLUMES_KERNEL_POLYHEDRONIMPLEMENTATION_H_
0006
0007 #include <cstdio>
0008
0009 #include "VecGeom/base/Vector3D.h"
0010 #include "VecGeom/volumes/kernel/GenericKernels.h"
0011 #include "VecGeom/volumes/kernel/TubeImplementation.h"
0012 #include "VecGeom/volumes/Quadrilaterals.h"
0013 #include "VecGeom/volumes/PolyhedronStruct.h"
0014
0015 namespace vecgeom {
0016
0017 // VECGEOM_DEVICE_FORWARD_DECLARE(struct PolyhedronImplementation;);
0018
0019 VECGEOM_DEVICE_DECLARE_CONV_TEMPLATE_2v(struct, PolyhedronImplementation, Polyhedron::EInnerRadii,
0020                                         Polyhedron::EInnerRadii::kGeneric, Polyhedron::EPhiCutout,
0021                                         Polyhedron::EPhiCutout::kGeneric);
0022
0023 inline namespace VECGEOM_IMPL_NAMESPACE {
0024
0025 class PlacedPolyhedron;
0026 class UnplacedPolyhedron;
0027
0028 template <Polyhedron::EInnerRadii innerRadiiT, Polyhedron::EPhiCutout phiCutoutT>
0029 struct PolyhedronImplementation {
0030
0031   using PlacedShape_t    = PlacedPolyhedron;
0032   using UnplacedStruct_t = PolyhedronStruct<Precision>;
0033   using UnplacedVolume_t = UnplacedPolyhedron;
0034
0035   VECCORE_ATT_HOST_DEVICE
0036   static void PrintType() {}
0037
0038   template <typename Stream>
0039   static void PrintType(Stream &s, int transCodeT = translation::kGeneric, int rotCodeT = rotation::kGeneric)
0040   {
0041     s << "SpecializedPolyhedron<" << transCodeT << "," << rotCodeT << ">";
0042   }
0043
0044   template <typename Stream>
0045   static void PrintImplementationType(Stream & /*s*/)
0046   {
0047   }
0048
0049   template <typename Stream>
0050   static void PrintUnplacedType(Stream & /*s*/)
0051   {
0052   }
0053
0054   /// \param pointZ Z-coordinate of a point.
0055   /// \return Index of the Z-segment in which the passed point is located. If
0056   ///         point is outside the polyhedron, -1 will be returned for Z smaller
0057   ///         than the first Z-plane, or N for Z larger than the last Z-plane,
0058   ///         where N is the amount of segments.
0059   template <typename Real_v>
0060   VECGEOM_FORCE_INLINE
0061   VECCORE_ATT_HOST_DEVICE
0062   static int FindZSegment(UnplacedStruct_t const &unplaced, Real_v const &pointZ);
0063
0064   /// \return Index of the phi-segment in which the passed point is located.
0065   ///         Assuming the polyhedron has been constructed properly, this should
0066   ///         always be a valid index.
0067   template <typename Real_v>
0068   VECGEOM_FORCE_INLINE
0069   VECCORE_ATT_HOST_DEVICE
0070   static int FindPhiSegment(UnplacedStruct_t const &unplaced, Vector3D<Real_v> const &point);
0071
0072   /// \param segmentIndex Index to the Z-segment to which the distance should be
0073   ///                     computed.
0074   /// \return Distance to the closest quadrilateral intersection by the passed
0075   ///         ray. Only intersections from the correct direction are accepted,
0076   ///         so value is always positive.
0077   template <typename Real_v>
0078   VECGEOM_FORCE_INLINE
0079   VECCORE_ATT_HOST_DEVICE
0080   static Real_v DistanceToInZSegment(UnplacedStruct_t const &unplaced, int segmentIndex, Vector3D<Real_v> const &point,
0081                                      Vector3D<Real_v> const &direction);
0082
0083   /// \param segmentIndex Index to the Z-segment to which the distance should be
0084   ///                     computed.
0085   /// \return Distance to the closest quadrilateral intersection by the passed
0086   ///         ray. Only intersections from the correct direction are accepted,
0087   ///         so value is always positive.
0088   template <typename Real_v>
0089   VECGEOM_FORCE_INLINE
0090   VECCORE_ATT_HOST_DEVICE
0091   static Real_v DistanceToOutZSegment(UnplacedStruct_t const &unplaced, int segmentIndex, Precision zMin,
0092                                       Precision zMax, Vector3D<Real_v> const &point, Vector3D<Real_v> const &direction);
0093
0094   /// \param segmentIndex Index to the Z-segment for which the safety should be
0095   ///        computed.
0096   /// \param phiIndex Index to the phi-segment for which the safety should be
0097   ///                 computed.
0098   /// \return Exact squared distance from the passed point to the quadrilateral
0099   ///         at the Z-segment and phi indices passed.
0100   VECCORE_ATT_HOST_DEVICE
0101   VECGEOM_FORCE_INLINE
0102   static Precision ScalarSafetyToZSegmentSquared(UnplacedStruct_t const &unplaced, int segmentIndex, int &phiIndex,
0103                                                  Vector3D<Precision> const &point, bool pt_inside, int &iSurf);
0104
0105   /// \param goingRight Whether the point is travelling along the Z-axis (true)
0106   ///        or opposite of the Z-axis (false).
0107   /// \param distance Output argument which will be minimized with the found
0108   ///                 distance.
0109   template <bool pointInsideT>
0110   VECGEOM_FORCE_INLINE
0111   VECCORE_ATT_HOST_DEVICE
0112   static void ScalarDistanceToEndcaps(UnplacedStruct_t const &unplaced, bool goingRight,
0113                                       Vector3D<Precision> const &point, Vector3D<Precision> const &direction,
0114                                       Precision &distance);
0115
0116   /// \brief Computes the exact distance to the closest endcap and minimizes it
0117   ///        with the output argument.
0118   /// \param distance Output argument which will be minimized with the found
0119   ///                 distance.
0120   VECCORE_ATT_HOST_DEVICE
0121   VECGEOM_FORCE_INLINE
0122   static void ScalarSafetyToEndcapsSquared(UnplacedStruct_t const &unplaced, Vector3D<Precision> const &point,
0123                                            Precision &distance, int &iz);
0124
0125   /// \param largePhiCutout Whether the phi cutout angle is larger than pi.
0126   /// \return Whether a point is within the infinite phi wedge formed from
0127   ///         origin in the cutout angle between the first and last vector.
0128   template <typename Real_v>
0129   VECGEOM_FORCE_INLINE
0130   VECCORE_ATT_HOST_DEVICE
0131   static vecCore::Mask_v<Real_v> InPhiCutoutWedge(ZSegment const &segment, bool largePhiCutout,
0132                                                   Vector3D<Real_v> const &point);
0133
0134   VECCORE_ATT_HOST_DEVICE
0135   VECGEOM_FORCE_INLINE
0136   static bool ScalarContainsKernel(UnplacedStruct_t const &unplaced, Vector3D<Precision> const &point);
0137
0138   VECCORE_ATT_HOST_DEVICE
0139   VECGEOM_FORCE_INLINE
0140   static Inside_t ScalarInsideKernel(UnplacedStruct_t const &unplaced, Vector3D<Precision> const &point);
0141
0142   VECCORE_ATT_HOST_DEVICE
0143   VECGEOM_FORCE_INLINE
0144   static Inside_t ScalarInsideSegPhi(UnplacedStruct_t const &unplaced, Vector3D<Precision> const &point, int zIndex,
0145                                      int phiIndex);
0146
0147   VECCORE_ATT_HOST_DEVICE
0148   VECGEOM_FORCE_INLINE
0149   static Inside_t ScalarInsideSegBorder(UnplacedStruct_t const &unplaced, Vector3D<Precision> const &point, int zIndex);
0150
0151   VECCORE_ATT_HOST_DEVICE
0152   VECGEOM_FORCE_INLINE
0153   static Precision ScalarDistanceToInKernel(UnplacedStruct_t const &unplaced, Vector3D<Precision> const &point,
0154                                             Vector3D<Precision> const &direction, const Precision stepMax);
0155
0156   VECCORE_ATT_HOST_DEVICE
0157   VECGEOM_FORCE_INLINE
0158   static Precision ScalarDistanceToOutKernel(UnplacedStruct_t const &unplaced, Vector3D<Precision> const &point,
0159                                              Vector3D<Precision> const &direction, const Precision stepMax);
0160
0161   VECCORE_ATT_HOST_DEVICE
0162   VECGEOM_FORCE_INLINE
0163   static Precision ScalarSafetyKernel(UnplacedStruct_t const &unplaced, Vector3D<Precision> const &point,
0164                                       bool pt_inside);
0165
0166   VECCORE_ATT_HOST_DEVICE
0167   VECGEOM_FORCE_INLINE
0168   static bool ScalarNormalKernel(UnplacedStruct_t const &unplaced, Vector3D<Precision> const &point,
0169                                  Vector3D<Precision> &normal);
0170
0171   /// Not implemented. Scalar version is called from SpecializedPolyhedron.
0172   template <typename Real_v, typename Bool_v>
0173   VECGEOM_FORCE_INLINE
0174   VECCORE_ATT_HOST_DEVICE
0175   static void UnplacedContains(UnplacedStruct_t const &unplaced, Vector3D<Real_v> const &point, Bool_v &inside);
0176
0177   /// Not implemented. Scalar version is called from SpecializedPolyhedron.
0178   template <typename Real_v, typename Bool_v>
0179   VECGEOM_FORCE_INLINE
0180   VECCORE_ATT_HOST_DEVICE
0181   static void Contains(UnplacedStruct_t const &unplaced, Vector3D<Real_v> const &point, Bool_v &inside);
0182
0183   /// Not implemented. Scalar version is called from Specializedunplaced.
0184   template <typename Real_v, typename Inside_v>
0185   VECGEOM_FORCE_INLINE
0186   VECCORE_ATT_HOST_DEVICE
0187   static void Inside(UnplacedStruct_t const &unplaced, Vector3D<Real_v> const &point, Inside_v &inside);
0188
0189   /// Not implemented. Scalar version is called from SpecializedPolyhedron.
0190   template <typename Real_v>
0191   VECGEOM_FORCE_INLINE
0192   VECCORE_ATT_HOST_DEVICE
0193   static void DistanceToIn(UnplacedStruct_t const &unplaced, Vector3D<Real_v> const &point,
0194                            Vector3D<Real_v> const &direction, Real_v const &stepMax, Real_v &distance);
0195
0196   /// Not implemented. Scalar version is called from SpecializedPolyhedron.
0197   template <typename Real_v>
0198   VECGEOM_FORCE_INLINE
0199   VECCORE_ATT_HOST_DEVICE
0200   static void DistanceToOut(UnplacedStruct_t const &unplaced, Vector3D<Real_v> const &point,
0201                             Vector3D<Real_v> const &direction, Real_v const &stepMax, Real_v &distance);
0202
0203   /// Not implemented. Scalar version is called from SpecializedPolyhedron.
0204   template <typename Real_v>
0205   VECGEOM_FORCE_INLINE
0206   VECCORE_ATT_HOST_DEVICE
0207   static void SafetyToIn(UnplacedStruct_t const &unplaced, Vector3D<Real_v> const &point, Real_v &safety);
0208
0209   /// Not implemented. Scalar version is called from SpecializedPolyhedron.
0210   template <typename Real_v>
0211   VECGEOM_FORCE_INLINE
0212   VECCORE_ATT_HOST_DEVICE
0213   static void SafetyToOut(UnplacedStruct_t const &unplaced, Vector3D<Real_v> const &point, Real_v &safety);
0214
0215 }; // End struct PolyhedronImplementation
0216
0217 namespace {
0218
0219 /// Polyhedron-specific trait class typedef'ing the tube specialization that
0220 /// should be called as a bounds check in Contains, Inside and DistanceToIn.
0221
0222 // SW (19.6.2015): switching to UniversalTube as Phi section was not
0223 // correctly treated with a hollow tube
0224 // TODO: this could be CORRECTLY put back for optimization
0225 template <Polyhedron::EInnerRadii innerRadiiT>
0226 struct HasInnerRadiiTraits {
0227   /// If polyhedron has inner radii, use a hollow tube
0228   typedef TubeImplementation<TubeTypes::UniversalTube> TubeKernels;
0229 };
0230
0231 template <>
0232 struct HasInnerRadiiTraits<Polyhedron::EInnerRadii::kFalse> {
0233   /// If polyhedron has no inner radii, use a non-hollow tube
0234   typedef TubeImplementation<TubeTypes::UniversalTube> TubeKernels;
0235 };
0236
0237 template <Polyhedron::EInnerRadii innerRadiiT>
0238 VECGEOM_FORCE_INLINE
0239 VECCORE_ATT_HOST_DEVICE
0240 bool TreatInner(bool hasInnerRadius)
0241 {
0242   return hasInnerRadius;
0243 }
0244
0245 template <>
0246 VECGEOM_FORCE_INLINE
0247 VECCORE_ATT_HOST_DEVICE
0248 bool TreatInner<Polyhedron::EInnerRadii::kFalse>(bool /*hasInnerRadius*/)
0249 {
0250   return false;
0251 }
0252
0253 template <Polyhedron::EPhiCutout phiCutoutT>
0254 VECGEOM_FORCE_INLINE
0255 VECCORE_ATT_HOST_DEVICE
0256 bool TreatPhi(bool /*hasPhiCutout*/)
0257 {
0258   return true;
0259 }
0260
0261 template <>
0262 VECGEOM_FORCE_INLINE
0263 VECCORE_ATT_HOST_DEVICE
0264 bool TreatPhi<Polyhedron::EPhiCutout::kFalse>(bool /*hasPhiCutout*/)
0265 {
0266   return false;
0267 }
0268
0269 template <>
0270 VECGEOM_FORCE_INLINE
0271 VECCORE_ATT_HOST_DEVICE
0272 bool TreatPhi<Polyhedron::EPhiCutout::kGeneric>(bool hasPhiCutout)
0273 {
0274   return hasPhiCutout;
0275 }
0276
0277 template <Polyhedron::EPhiCutout phiCutoutT>
0278 VECGEOM_FORCE_INLINE
0279 VECCORE_ATT_HOST_DEVICE
0280 bool LargePhiCutout(bool largePhiCutout)
0281 {
0282   return largePhiCutout;
0283 }
0284
0285 template <>
0286 VECGEOM_FORCE_INLINE
0287 VECCORE_ATT_HOST_DEVICE
0288 bool LargePhiCutout<Polyhedron::EPhiCutout::kTrue>(bool /*largePhiCutout*/)
0289 {
0290   return false;
0291 }
0292
0293 template <>
0294 VECGEOM_FORCE_INLINE
0295 VECCORE_ATT_HOST_DEVICE
0296 bool LargePhiCutout<Polyhedron::EPhiCutout::kLarge>(bool /*largePhiCutout*/)
0297 {
0298   return true;
0299 }
0300
0301 } // End anonymous namespace
0302
0303 namespace {
0304
0305 template <typename Real_v>
0306 VECGEOM_FORCE_INLINE
0307 VECCORE_ATT_HOST_DEVICE
0308 int FindZSegmentKernel(Precision const *begin, Precision const *end, Real_v const &pointZ);
0309
0310 template <>
0311 VECGEOM_FORCE_INLINE
0312 VECCORE_ATT_HOST_DEVICE
0313 int FindZSegmentKernel<Precision>(Precision const *begin, Precision const *end, Precision const &pointZ)
0314 {
0315   // TODO: vectorize this and move the brute-force algorithm to the CUDA
0316   //       implementation. Inspiration can be found at:
0317   //       http://schani.wordpress.com/2010/04/30/linear-vs-binary-search/
0318   int index = -1;
0319   // Modified algorithm to select the first section the position is close to
0320   // within boundary tolerance. This is important for degenerated Z polyhedra
0321   while (begin < end - 1 && pointZ - kTolerance > *begin) {
0322     ++index;
0323     ++begin;
0324   }
0325   if (pointZ + kTolerance > *begin) return (index + 1);
0326   return index;
0327 }
0328 } // End anonymous namespace
0329
0330 template <Polyhedron::EInnerRadii innerRadiiT, Polyhedron::EPhiCutout phiCutoutT>
0331 template <typename Real_v>
0332 VECGEOM_FORCE_INLINE
0333 VECCORE_ATT_HOST_DEVICE
0334 int PolyhedronImplementation<innerRadiiT, phiCutoutT>::FindZSegment(UnplacedStruct_t const &unplaced,
0335                                                                     Real_v const &pointZ)
0336 {
0337   return FindZSegmentKernel<Real_v>(&unplaced.fZPlanes[0], &unplaced.fZPlanes[0] + unplaced.fZPlanes.size(), pointZ);
0338 }
0339
0340 template <Polyhedron::EInnerRadii innerRadiiT, Polyhedron::EPhiCutout phiCutoutT>
0341 template <typename Real_v>
0342 VECGEOM_FORCE_INLINE
0343 VECCORE_ATT_HOST_DEVICE
0344 int PolyhedronImplementation<innerRadiiT, phiCutoutT>::FindPhiSegment(UnplacedStruct_t const &unplaced,
0345                                                                       Vector3D<Real_v> const &point)
0346 {
0347
0348   // Bounds between phi sections are represented as planes through the origin,
0349   // with the normal pointing along the phi direction.
0350   // To find the correct section, the point is projected onto each plane. If the
0351   // point is in front of a plane, but behind the subsequent plane, it must be
0352   // between them.
0353
0354   int index                           = -1;
0355   SOA3D<Precision> const &phiSections = unplaced.fPhiSections;
0356   Real_v projectionFirst, projectionSecond;
0357   projectionFirst = point[0] * phiSections.x(0) + point[1] * phiSections.y(0) + point[2] * phiSections.z(0);
0358   for (int i = 1, iMax = unplaced.fSideCount + 1; i < iMax; ++i) {
0359     projectionSecond = point[0] * phiSections.x(i) + point[1] * phiSections.y(i) + point[2] * phiSections.z(i);
0360     vecCore__MaskedAssignFunc(index, projectionFirst > -kTolerance && projectionSecond < kTolerance, i - 1);
0361     if (vecCore::MaskFull(index >= 0)) break;
0362     projectionFirst = projectionSecond;
0363   }
0364
0365   return index;
0366 }
0367
0368 template <Polyhedron::EInnerRadii innerRadiiT, Polyhedron::EPhiCutout phiCutoutT>
0369 template <typename Real_v>
0370 VECCORE_ATT_HOST_DEVICE
0371 Real_v PolyhedronImplementation<innerRadiiT, phiCutoutT>::DistanceToInZSegment(UnplacedStruct_t const &unplaced,
0372                                                                                int segmentIndex,
0373                                                                                Vector3D<Real_v> const &point,
0374                                                                                Vector3D<Real_v> const &direction)
0375 {
0376
0377   using Bool_v = vecCore::Mask_v<Real_v>;
0378
0379   Real_v distance;
0380   Bool_v done;
0381
0382   ZSegment const &segment = unplaced.fZSegments[segmentIndex];
0383
0384   // If the outer shell is hit, this will always be the correct result
0385   distance = segment.outer.DistanceToIn<Real_v, false>(point, direction);
0386   done     = distance < InfinityLength<Real_v>();
0387   if (vecCore::MaskFull(done)) return distance;
0388
0389   // If the outer shell is not hit and the phi cutout sides are hit, this will
0390   // always be the correct result
0391   if (TreatPhi<phiCutoutT>(unplaced.fHasPhiCutout)) {
0392     vecCore__MaskedAssignFunc(distance, !done, (segment.phi.DistanceToIn<Real_v, false>(point, direction)));
0393     if (unplaced.fHasLargePhiCutout) {
0394       // NOTE: The statement above is NOT always true: if fHasLargePhiCutout is false there can be a first hit of the
0395       // inner surface coming from the endcap holes
0396       done |= distance < InfinityLength<Real_v>();
0397       if (vecCore::MaskFull(done)) return distance;
0398     }
0399   }
0400
0401   // Finally treat inner shell
0402   if (TreatInner<innerRadiiT>(segment.hasInnerRadius())) {
0403     Real_v distrmin = segment.inner.DistanceToIn<Real_v, true>(point, direction);
0404     vecCore__MaskedAssignFunc(distance, !done && distance > distrmin, distrmin);
0405   }
0406
0407   return distance;
0408 }
0409
0410 template <Polyhedron::EInnerRadii innerRadiiT, Polyhedron::EPhiCutout phiCutoutT>
0411 template <typename Real_v>
0412 VECCORE_ATT_HOST_DEVICE
0413 Real_v PolyhedronImplementation<innerRadiiT, phiCutoutT>::DistanceToOutZSegment(UnplacedStruct_t const &unplaced,
0414                                                                                 int segmentIndex, Precision zMin,
0415                                                                                 Precision zMax,
0416                                                                                 Vector3D<Real_v> const &point,
0417                                                                                 Vector3D<Real_v> const &direction)
0418 {
0419
0420   using Bool_v = vecCore::Mask_v<Real_v>;
0421
0422   Bool_v done(false);
0423   Real_v distance = InfinityLength<Real_v>();
0424
0425   ZSegment const &segment = unplaced.fZSegments[segmentIndex];
0426
0427   // Check inner shell first, as it would always be the correct result
0428   if (TreatInner<innerRadiiT>(segment.hasInnerRadius())) {
0429     distance = segment.inner.DistanceToIn<Real_v, false>(point, direction);
0430     // Even if an inner surface is hit, there may be a phi hit before if there is no large phi cut
0431     if (unplaced.fHasLargePhiCutout) {
0432       done = distance < InfinityLength<Real_v>();
0433       if (vecCore::MaskFull(done)) return distance;
0434     }
0435   }
0436
0437   // Check phi cutout if necessary. It is also possible to return here if a
0438   // result is found
0439   if (TreatPhi<phiCutoutT>(unplaced.fHasPhiCutout)) {
0440     Real_v distphi = segment.phi.DistanceToIn<Real_v, true>(point, direction);
0441     vecCore::MaskedAssign(distance, distance > distphi, distphi);
0442   }
0443
0444   done = distance > -kTolerance && distance < InfinityLength<Real_v>();
0445   if (vecCore::MaskFull(done)) return distance;
0446
0447   // Finally check outer shell
0448   Real_v distout = segment.outer.DistanceToOut<Real_v>(point, direction, zMin, zMax);
0449   vecCore::MaskedAssign(distance, !done, distout);
0450
0451   return distance;
0452 }
0453
0454 template <Polyhedron::EInnerRadii innerRadiiT, Polyhedron::EPhiCutout phiCutoutT>
0455 VECGEOM_FORCE_INLINE
0456 VECCORE_ATT_HOST_DEVICE
0457 Precision PolyhedronImplementation<innerRadiiT, phiCutoutT>::ScalarSafetyToZSegmentSquared(
0458     UnplacedStruct_t const &unplaced, int segmentIndex, int &phiIndex, Vector3D<Precision> const &point, bool pt_inside,
0459     int &iSurf)
0460 {
0461
0462   ZSegment const &segment = unplaced.fZSegments[segmentIndex];
0463   bool in_cutout          = phiIndex < 0;
0464
0465   Precision safetySquared = InfinityLength<Precision>();
0466   if (TreatPhi<phiCutoutT>(unplaced.fHasPhiCutout) && segment.phi.size() == 2) {
0467     //  Check if points is in the cutout wedge first.
0468     if (pt_inside || in_cutout) {
0469       // If point is in the cutout or if the call comes from SafetyToOut we need to check both phi planes
0470       iSurf         = 0;
0471       safetySquared = segment.phi.ScalarDistanceSquared(0, point);
0472       Precision saf = segment.phi.ScalarDistanceSquared(1, point);
0473       if (saf < safetySquared) {
0474         safetySquared = saf;
0475         iSurf         = 1;
0476       }
0477       // If the point is within the phi cutout wedge, we still need to check the
0478       // inner part
0479       if (in_cutout) {
0480         if (TreatInner<innerRadiiT>(segment.hasInnerRadius())) {
0481           if (segment.inner.size() > 0) {
0482             Precision safetySquaredInner = segment.inner.ScalarDistanceSquared(0, point);
0483             if (safetySquaredInner < safetySquared) {
0484               iSurf         = 2;
0485               phiIndex      = 0;
0486               safetySquared = safetySquaredInner;
0487             }
0488             if (segment.inner.size() > 1) {
0489               safetySquaredInner = segment.inner.ScalarDistanceSquared(segment.inner.size() - 1, point);
0490               if (safetySquaredInner < safetySquared) {
0491                 iSurf         = 2;
0492                 phiIndex      = segment.inner.size() - 1;
0493                 safetySquared = safetySquaredInner;
0494               }
0495             }
0496           }
0497         }
0498         return safetySquared;
0499       }
0500     }
0501   }
0502
0503   if (in_cutout && segmentIndex > 0 && segmentIndex < unplaced.fZSegments.size() - 1 &&
0504       unplaced.fZPlanes[segmentIndex] == unplaced.fZPlanes[segmentIndex + 1]) {
0505     // We are checking a segment at same Z. We have to check the inner and outer
0506     // quadrilaterals for first and last phi
0507     Precision safetySquaredOuter = InfinityLength<Precision>();
0508     if (segment.outer.size() > 0) {
0509       safetySquaredOuter = segment.outer.ScalarDistanceSquared(0, point);
0510       if (safetySquaredOuter < safetySquared) {
0511         iSurf         = 3;
0512         phiIndex      = 0;
0513         safetySquared = safetySquaredOuter;
0514       }
0515       if (segment.outer.size() > 1) {
0516         safetySquaredOuter = segment.outer.ScalarDistanceSquared(segment.outer.size() - 1, point);
0517         if (safetySquaredOuter < safetySquared) {
0518           iSurf         = 3;
0519           phiIndex      = segment.outer.size() - 1;
0520           safetySquared = safetySquaredOuter;
0521         }
0522       }
0523     }
0524     Precision safetySquaredInner = InfinityLength<Precision>();
0525     if (TreatInner<innerRadiiT>(segment.hasInnerRadius())) {
0526       if (segment.inner.size() > 0) {
0527         safetySquaredInner = segment.inner.ScalarDistanceSquared(0, point);
0528         if (safetySquaredInner < safetySquared) {
0529           iSurf         = 2;
0530           phiIndex      = 0;
0531           safetySquared = safetySquaredInner;
0532         }
0533         if (segment.inner.size() > 1) {
0534           safetySquaredInner = segment.inner.ScalarDistanceSquared(segment.inner.size() - 1, point);
0535           if (safetySquaredInner < safetySquared) {
0536             iSurf         = 2;
0537             phiIndex      = segment.inner.size() - 1;
0538             safetySquared = safetySquaredInner;
0539           }
0540         }
0541       }
0542     }
0543     return safetySquared;
0544   }
0545
0546   // Otherwise check the outer shell
0547   // TODO: we need to check segment.outer.size() > 0
0548   Precision safetySquaredOuter = InfinityLength<Precision>();
0549   if (segment.outer.size() > 0) safetySquaredOuter = segment.outer.ScalarDistanceSquared(phiIndex, point);
0550
0551   // And finally the inner
0552   Precision safetySquaredInner = InfinityLength<Precision>();
0553   if (TreatInner<innerRadiiT>(segment.hasInnerRadius())) {
0554     if (segment.inner.size() > 0) safetySquaredInner = segment.inner.ScalarDistanceSquared(phiIndex, point);
0555   }
0556   if (safetySquaredInner < safetySquared) {
0557     iSurf         = 2;
0558     safetySquared = safetySquaredInner;
0559   }
0560   if (safetySquaredOuter < safetySquared) {
0561     iSurf         = 3;
0562     safetySquared = safetySquaredOuter;
0563   }
0564   return safetySquared;
0565 }
0566
0567 template <Polyhedron::EInnerRadii innerRadiiT, Polyhedron::EPhiCutout phiCutoutT>
0568 template <bool pointInsideT>
0569 VECGEOM_FORCE_INLINE
0570 VECCORE_ATT_HOST_DEVICE
0571 void PolyhedronImplementation<innerRadiiT, phiCutoutT>::ScalarDistanceToEndcaps(UnplacedStruct_t const &unplaced,
0572                                                                                 bool /*goingRight*/,
0573                                                                                 Vector3D<Precision> const &point,
0574                                                                                 Vector3D<Precision> const &direction,
0575                                                                                 Precision &distance)
0576 {
0577
0578   ZSegment const *segment;
0579   Precision zPlane;
0580
0581   // Determine whether to use first segment/first endcap or last segment/second
0582   // endcap
0583   // NOTE: might make this more elegant
0584   if (pointInsideT) // inside version
0585   {
0586     if (direction[2] < 0) {
0587       segment = &unplaced.fZSegments[0];
0588       zPlane  = unplaced.fZPlanes[0];
0589     } else {
0590       segment = &unplaced.fZSegments[unplaced.fZSegments.size() - 1];
0591       zPlane  = unplaced.fZPlanes[unplaced.fZSegments.size()];
0592     }
0593   } else // outside version
0594   {
0595     if (direction[2] < 0) {
0596       segment = &unplaced.fZSegments[unplaced.fZSegments.size() - 1];
0597       zPlane  = unplaced.fZPlanes[unplaced.fZSegments.size()];
0598     } else {
0599       segment = &unplaced.fZSegments[0];
0600       zPlane  = unplaced.fZPlanes[0];
0601     }
0602   }
0603
0604   Precision distanceTest = (zPlane - point[2]) / NonZero(direction[2]);
0605   // If the distance is not better there's no reason to check for validity
0606   if (distanceTest < -kTolerance || distanceTest >= distance) return;
0607
0608   Vector3D<Precision> intersection = point + distanceTest * direction;
0609   // Intersection point must be inside outer shell and outside inner shell
0610   if (!segment->outer.Contains<Precision>(intersection)) return;
0611   if (TreatInner<innerRadiiT>(segment->hasInnerRadius())) {
0612     if (segment->inner.Contains<Precision>(intersection)) return;
0613   }
0614   // Intersection point must not be in phi cutout wedge
0615   if (TreatPhi<phiCutoutT>(unplaced.fHasPhiCutout)) {
0616     if (InPhiCutoutWedge<Precision>(*segment, unplaced.fHasLargePhiCutout, intersection)) {
0617       return;
0618     }
0619   }
0620
0621   distance = distanceTest;
0622 }
0623
0624 template <Polyhedron::EInnerRadii innerRadiiT, Polyhedron::EPhiCutout phiCutoutT>
0625 VECGEOM_FORCE_INLINE
0626 VECCORE_ATT_HOST_DEVICE
0627 void PolyhedronImplementation<innerRadiiT, phiCutoutT>::ScalarSafetyToEndcapsSquared(UnplacedStruct_t const &unplaced,
0628                                                                                      Vector3D<Precision> const &point,
0629                                                                                      Precision &distanceSquared,
0630                                                                                      int &iz)
0631 {
0632
0633   // Compute both distances (simple subtractions) to determine which is closer
0634   Precision firstDistance = unplaced.fZPlanes[0] - point[2];
0635   Precision lastDistance  = unplaced.fZPlanes[unplaced.fZSegments.size()] - point[2];
0636
0637   // Only treat the closest endcap
0638   bool isFirst            = Abs(firstDistance) < Abs(lastDistance);
0639   iz                      = 0;
0640   ZSegment const &segment = isFirst ? unplaced.fZSegments[0] : unplaced.fZSegments[unplaced.fZSegments.size() - 1];
0641
0642   Precision distanceTest        = isFirst ? firstDistance : lastDistance;
0643   Precision distanceTestSquared = distanceTest * distanceTest;
0644   // No need to investigate further if distance is larger anyway
0645   if (distanceTestSquared >= distanceSquared) return;
0646
0647   // Check if projection is within the endcap bounds
0648   Vector3D<Precision> intersection(point[0], point[1], point[2] + distanceTest);
0649   if (!segment.outer.Contains<Precision>(intersection)) return;
0650   if (TreatInner<innerRadiiT>(segment.hasInnerRadius())) {
0651     if (segment.inner.Contains<Precision>(intersection)) return;
0652   }
0653   if (TreatPhi<phiCutoutT>(unplaced.fHasPhiCutout)) {
0654     if (InPhiCutoutWedge<Precision>(segment, unplaced.fHasLargePhiCutout, intersection)) {
0655       return;
0656     }
0657   }
0658
0659   iz              = (isFirst) ? -1 : 1;
0660   distanceSquared = distanceTestSquared;
0661 }
0662
0663 template <Polyhedron::EInnerRadii innerRadiiT, Polyhedron::EPhiCutout phiCutoutT>
0664 template <typename Real_v>
0665 VECGEOM_FORCE_INLINE
0666 VECCORE_ATT_HOST_DEVICE
0667 vecCore::Mask_v<Real_v> PolyhedronImplementation<innerRadiiT, phiCutoutT>::InPhiCutoutWedge(
0668     ZSegment const &segment, bool largePhiCutout, Vector3D<Real_v> const &point)
0669 {
0670   using Bool_v     = vecCore::Mask_v<Real_v>;
0671   Bool_v pointSeg0 = point.Dot(segment.phi.GetNormal(0)) + segment.phi.GetDistance(0) >= 0;
0672   Bool_v pointSeg1 = point.Dot(segment.phi.GetNormal(1)) + segment.phi.GetDistance(1) >= 0;
0673   // For a cutout larger than 180 degrees, the point is in the wedge if it is
0674   // in front of at least one plane.
0675   if (LargePhiCutout<phiCutoutT>(largePhiCutout)) {
0676     return pointSeg0 || pointSeg1;
0677   }
0678   // Otherwise it should be in front of both planes
0679   return pointSeg0 && pointSeg1;
0680 }
0681
0682 template <Polyhedron::EInnerRadii innerRadiiT, Polyhedron::EPhiCutout phiCutoutT>
0683 VECCORE_ATT_HOST_DEVICE
0684 bool PolyhedronImplementation<innerRadiiT, phiCutoutT>::ScalarContainsKernel(UnplacedStruct_t const &unplaced,
0685                                                                              Vector3D<Precision> const &point)
0686 {
0687
0688   // First check if in bounding tube
0689   {
0690     bool inBounds;
0691     // Correct tube algorithm obtained from trait class
0692     HasInnerRadiiTraits<innerRadiiT>::TubeKernels::template Contains(
0693         unplaced.fBoundingTube, Vector3D<Precision>(point[0], point[1], point[2] - unplaced.fBoundingTubeOffset),
0694         inBounds);
0695     if (!inBounds) return false;
0696   }
0697
0698   // Find correct segment by checking Z-bounds
0699   int zIndex = FindZSegment<Precision>(unplaced, point[2]);
0700   if (!((zIndex >= 0) && (zIndex < unplaced.fZSegments.size()))) return false;
0701
0702   ZSegment const &segment = unplaced.fZSegments[zIndex];
0703
0704   // In case the point lies at the same Z as 2 consecutive planes, the lesser
0705   // index is selected. The Quadrilaterals algorithm for Contains in this case
0706   // does not work.
0707   if (unplaced.fSameZ[zIndex]) {
0708     // Identify phi index
0709     int phiIndex = FindPhiSegment<Precision>(unplaced, point);
0710     if (phiIndex < 0) return false;
0711     // Get the vector perpendicular to the rmax edge of the outer quadrilateral
0712     Vector3D<Precision> const &vout = (segment.outer.size()) ? segment.outer.GetSideVectors()[0].GetNormals()[phiIndex]
0713                                                              : segment.inner.GetSideVectors()[0].GetNormals()[phiIndex];
0714     // Compute the projection of the point vectoron the vout vector. This
0715     // corresponds to a "radius" or the point.
0716     Precision rdotvout = vecCore::math::Abs<Precision>(point.Dot(vout));
0717     // Now compare the point radius with the ranges corresponding to the lower
0718     // and upper segments
0719     bool in1 = (rdotvout >= unplaced.fRMin[zIndex]) && (rdotvout <= unplaced.fRMax[zIndex]);
0720     bool in2 = (rdotvout >= unplaced.fRMin[zIndex + 1]) && (rdotvout <= unplaced.fRMax[zIndex + 1]);
0721     return (in1 | in2);
0722   }
0723
0724   // Check that the point is in the outer shell
0725   if (!segment.outer.Contains<Precision>(point)) return false;
0726
0727   // Check that the point is not in the inner shell
0728   if (TreatInner<innerRadiiT>(segment.hasInnerRadius())) {
0729     if (segment.inner.Contains<Precision>(point)) return false;
0730   }
0731
0732   // In principle, handling of phi should not be needed here since it is
0733   // contained in the bounding tube. However, we need to check again due
0734   // to different handling of tolerances.
0735   if (TreatPhi<phiCutoutT>(unplaced.fHasPhiCutout)) {
0736     if (!segment.phi.Contains<Precision>(point)) return false;
0737   }
0738
0739   return true;
0740 }
0741
0742 // TODO: check this code -- maybe unify with previous function
0743 template <Polyhedron::EInnerRadii innerRadiiT, Polyhedron::EPhiCutout phiCutoutT>
0744 VECCORE_ATT_HOST_DEVICE
0745 Inside_t PolyhedronImplementation<innerRadiiT, phiCutoutT>::ScalarInsideKernel(UnplacedStruct_t const &unplaced,
0746                                                                                Vector3D<Precision> const &point)
0747 {
0748
0749   // First check if in bounding tube
0750   {
0751     bool inBounds;
0752     // Correct tube algorithm obtained from trait class
0753     // FIX: the bounding tube was wrong. Since the fast UnplacedContains is
0754     // used for early return, the bounding tube has to be larger than the
0755     // ideal bounding tube to account for the tolerance (offset was wrong)
0756     HasInnerRadiiTraits<innerRadiiT>::TubeKernels::template Contains(
0757         unplaced.fBoundingTube, Vector3D<Precision>(point[0], point[1], point[2] - unplaced.fBoundingTubeOffset),
0758         inBounds);
0759     if (!inBounds) return EInside::kOutside;
0760   }
0761
0762   // Find correct segment by checking Z-bounds
0763   // The FindZSegment was fixed for the degenerated Z case when 2 planes
0764   // have identical Z. In this case, if the point is close within tolerance
0765   // to such section, the returned index has to be the first of the 2, so that
0766   // all navigation functions start by checking the degenerated segment.
0767   int zIndex = FindZSegment<Precision>(unplaced, point[2]);
0768   if (zIndex > (unplaced.fZSegments.size() - 1)) zIndex = unplaced.fZSegments.size() - 1;
0769   if (zIndex < 0) zIndex = 0;
0770
0771   ZSegment const &segment = unplaced.fZSegments[zIndex];
0772
0773   // Point in between 2 planes at same Z
0774   if (unplaced.fSameZ[zIndex]) return ScalarInsideSegBorder(unplaced, point, zIndex);
0775
0776   // Check that the point is in the outer shell
0777   {
0778     Inside_t insideOuter = segment.outer.Inside<Precision, Inside_t>(point);
0779     if (insideOuter != EInside::kInside) return insideOuter;
0780   }
0781
0782   // Check that the point is not in the inner shell
0783   if (TreatInner<innerRadiiT>(segment.hasInnerRadius())) {
0784     Inside_t insideInner = segment.inner.Inside<Precision, Inside_t>(point);
0785     if (insideInner == EInside::kInside) return EInside::kOutside;
0786     if (insideInner == EInside::kSurface) return EInside::kSurface;
0787   }
0788
0789   // Check that the point is not in the phi cutout wedge
0790   if (TreatPhi<phiCutoutT>(unplaced.fHasPhiCutout)) {
0791     // Inside_t insidePhi = unplaced.fPhiWedge.Inside<Precision, Inside_t>(point);
0792     Inside_t insidePhi = segment.phi.Inside<Precision, Inside_t>(point);
0793     if (insidePhi != EInside::kInside) return insidePhi;
0794   }
0795
0796   // FIX: Still need to check if not on one of the Z boundaries.
0797   Precision dz = vecCore::math::Abs(vecCore::math::Abs(point[2] - unplaced.fBoundingTubeOffset) -
0798                                     0.5 * (unplaced.fZPlanes[unplaced.fZSegments.size()] - unplaced.fZPlanes[0]));
0799   if (dz < kTolerance) return EInside::kSurface;
0800   return EInside::kInside;
0801 }
0802
0803 template <Polyhedron::EInnerRadii innerRadiiT, Polyhedron::EPhiCutout phiCutoutT>
0804 VECCORE_ATT_HOST_DEVICE
0805 Inside_t PolyhedronImplementation<innerRadiiT, phiCutoutT>::ScalarInsideSegBorder(UnplacedStruct_t const &unplaced,
0806                                                                                   Vector3D<Precision> const &point,
0807                                                                                   int zIndex)
0808 {
0809   // Check Inside if the point is in between two non-continuous "border-like"
0810   // segments. The zIndex corresponds to the lesser index of the 2 planes having the same Z.
0811   // The Quadrilaterals algorithm for Inside in this case does not work.
0812
0813   ZSegment const &segment = unplaced.fZSegments[zIndex];
0814   // Identify phi index
0815   int phiIndex = FindPhiSegment<Precision>(unplaced, point);
0816   if (phiIndex < 0) return EInside::kOutside;
0817   // Get the vector perpendicular to the rmax edge of the outer quadrilateral
0818   Vector3D<Precision> const &vout = (segment.outer.size()) ? segment.outer.GetSideVectors()[0].GetNormals()[phiIndex]
0819                                                            : segment.inner.GetSideVectors()[0].GetNormals()[phiIndex];
0820   // Compute the projection of the point vectoron the vout vector. This
0821   // corresponds to a "radius" or the point.
0822   Precision rdotvout = vecCore::math::Abs<Precision>(point.Dot(vout));
0823   // Now compare the point radius with the ranges corresponding to the lower
0824   // and upper segments
0825   bool in1 = (rdotvout > unplaced.fRMin[zIndex] - kTolerance) && (rdotvout < unplaced.fRMax[zIndex] + kTolerance);
0826   bool in2 =
0827       (rdotvout > unplaced.fRMin[zIndex + 1] - kTolerance) && (rdotvout < unplaced.fRMax[zIndex + 1] + kTolerance);
0828   if (in1 && in2) {
0829     if ((rdotvout < unplaced.fRMin[zIndex] + kTolerance) || (rdotvout > unplaced.fRMax[zIndex] - kTolerance) ||
0830         (rdotvout < unplaced.fRMin[zIndex + 1] + kTolerance) || (rdotvout > unplaced.fRMax[zIndex + 1] - kTolerance))
0831       return EInside::kSurface;
0832     // Need to check phi surface
0833     if (TreatPhi<phiCutoutT>(unplaced.fHasPhiCutout)) {
0834       Inside_t insidePhi = unplaced.fPhiWedge.Inside<Precision, Inside_t>(point);
0835       return insidePhi;
0836     }
0837     return EInside::kInside;
0838   }
0839   if (!in1 && !in2) return EInside::kOutside;
0840   return EInside::kSurface;
0841 }
0842
0843 template <Polyhedron::EInnerRadii innerRadiiT, Polyhedron::EPhiCutout phiCutoutT>
0844 VECCORE_ATT_HOST_DEVICE
0845 Inside_t PolyhedronImplementation<innerRadiiT, phiCutoutT>::ScalarInsideSegPhi(UnplacedStruct_t const &unplaced,
0846                                                                                Vector3D<Precision> const &point,
0847                                                                                int zIndex, int phiIndex)
0848 {
0849   // Check inside for a specified z segment and phi edge
0850   if (phiIndex < 0) return EInside::kOutside;
0851
0852   // Z range
0853   Precision dz = vecCore::math::Abs(point[2] - unplaced.fBoundingTubeOffset) -
0854                  0.5 * (unplaced.fZPlanes[unplaced.fZSegments.size()] - unplaced.fZPlanes[0]);
0855   //  if (vecCore::math::Abs(dz) < kHalfTolerance) return EInside::kSurface;
0856   if (dz > kHalfTolerance) return EInside::kOutside;
0857
0858   if (unplaced.fSameZ[zIndex]) return ScalarInsideSegBorder(unplaced, point, zIndex);
0859
0860   ZSegment const &segment = unplaced.fZSegments[zIndex];
0861
0862   // Check that the point is in the outer shell
0863   {
0864     Inside_t insideOuter = segment.outer.Inside<Precision, Inside_t>(point, phiIndex);
0865     if (insideOuter != EInside::kInside) return insideOuter;
0866   }
0867
0868   // Check that the point is not in the inner shell
0869   if (TreatInner<innerRadiiT>(segment.hasInnerRadius())) {
0870     Inside_t insideInner = segment.inner.Inside<Precision, Inside_t>(point, phiIndex);
0871     if (insideInner == EInside::kInside) return EInside::kOutside;
0872     if (insideInner == EInside::kSurface) return EInside::kSurface;
0873   }
0874
0875   // Check that the point is not in the phi cutout wedge
0876   if (TreatPhi<phiCutoutT>(unplaced.fHasPhiCutout)) {
0877     Inside_t insidePhi = unplaced.fPhiWedge.Inside<Precision, Inside_t>(point);
0878     if (insidePhi != EInside::kInside) return insidePhi;
0879   }
0880
0881   if (vecCore::math::Abs(dz) < kHalfTolerance) return EInside::kSurface;
0882   return EInside::kInside;
0883 }
0884
0885 template <Polyhedron::EInnerRadii innerRadiiT, Polyhedron::EPhiCutout phiCutoutT>
0886 VECCORE_ATT_HOST_DEVICE
0887 Precision PolyhedronImplementation<innerRadiiT, phiCutoutT>::ScalarDistanceToInKernel(
0888     UnplacedStruct_t const &unplaced, Vector3D<Precision> const &point, Vector3D<Precision> const &direction,
0889     const Precision stepMax)
0890 {
0891
0892   // Fast exclude points beyond endcaps moving on same side as endcap normal
0893   if ((point[2] < unplaced.fZPlanes[0] + kTolerance) && direction[2] <= 0) return InfinityLength<Precision>();
0894   if ((point[2] > unplaced.fZPlanes[unplaced.fZSegments.size()] - kTolerance) && direction[2] >= 0)
0895     return InfinityLength<Precision>();
0896
0897   // Perform explicit Inside check to detect wrong side points. This impacts
0898   // DistanceToIn performance by about 5% for all topologies
0899   auto inside = ScalarInsideKernel(unplaced, point);
0900   if (inside == kInside) return -1.;
0901
0902   // Check if the point is within the bounding tube
0903   bool inBounds;
0904   Precision tubeDistance = 0.;
0905   {
0906     Vector3D<Precision> boundsPoint(point[0], point[1], point[2] - unplaced.fBoundingTubeOffset);
0907     HasInnerRadiiTraits<innerRadiiT>::TubeKernels::template Contains(unplaced.fBoundingTube, boundsPoint, inBounds);
0908     // If the point is inside the bounding tube, the result of DistanceToIn is
0909     // unreliable and cannot be used to reject rays.
0910     // TODO: adjust tube DistanceToIn function to correctly return a negative
0911     //       value for points inside the tube. This will allow the removal of
0912     //       the contains check here.
0913     if (!inBounds) {
0914       // If the point is outside the bounding tube, check if the ray misses
0915       // the bounds
0916       HasInnerRadiiTraits<innerRadiiT>::TubeKernels::template DistanceToIn(unplaced.fBoundingTube, boundsPoint,
0917                                                                            direction, stepMax, tubeDistance);
0918       if (tubeDistance == InfinityLength<Precision>()) {
0919         return InfinityLength<Precision>();
0920       }
0921     }
0922   }
0923
0924   int zIndex     = FindZSegment<Precision>(unplaced, point[2]);
0925   const int zMax = unplaced.fZSegments.size();
0926   // Don't go out of bounds here, as the first/last segment should be checked
0927   // even if the point is outside of Z-bounds
0928   zIndex = zIndex < 0 ? 0 : (zIndex >= zMax ? zMax - 1 : zIndex);
0929
0930   // Traverse Z-segments left or right depending on sign of direction
0931   bool goingRight = direction[2] >= 0;
0932
0933   Precision distance = InfinityLength<Precision>();
0934   if (goingRight) {
0935     for (int zSegCount = unplaced.fZSegments.size(); zIndex < zSegCount; ++zIndex) {
0936       distance = DistanceToInZSegment<Precision>(unplaced, zIndex, point, direction);
0937       // No segment further away can be at a shorter distance to the point, so
0938       // if a valid distance is found, only endcaps remain to be investigated
0939       if (distance >= 0 && distance < InfinityLength<Precision>()) break;
0940     }
0941   } else {
0942     // Going left
0943     for (; zIndex >= 0; --zIndex) {
0944       distance = DistanceToInZSegment<Precision>(unplaced, zIndex, point, direction);
0945       // No segment further away can be at a shorter distance to the point, so
0946       // if a valid distance is found, only endcaps remain to be investigated
0947       if (distance >= 0 && distance < InfinityLength<Precision>()) break;
0948     }
0949   }
0950
0951   // Minimize with distance to endcaps
0952   ScalarDistanceToEndcaps<false>(unplaced, goingRight, point, direction, distance);
0953
0954   // last sanity check: distance should be larger than estimate from bounding tube
0955   return (distance >= tubeDistance - 1E-6) ? distance : vecgeom::InfinityLength<Precision>();
0956 }
0957
0958 template <Polyhedron::EInnerRadii innerRadiiT, Polyhedron::EPhiCutout phiCutoutT>
0959 VECCORE_ATT_HOST_DEVICE
0960 Precision PolyhedronImplementation<innerRadiiT, phiCutoutT>::ScalarSafetyKernel(UnplacedStruct_t const &unplaced,
0961                                                                                 Vector3D<Precision> const &point,
0962                                                                                 bool pt_inside)
0963 {
0964
0965   Precision safety = InfinityLength<Precision>();
0966   Precision dz;
0967   int iSurf, iz;
0968
0969   const int zMax = unplaced.fZSegments.size();
0970   int zIndex     = FindZSegment<Precision>(unplaced, point[2]);
0971   zIndex         = zIndex < 0 ? 0 : (zIndex >= zMax ? zMax - 1 : zIndex);
0972
0973   int phiIndex = FindPhiSegment<Precision>(unplaced, point);
0974
0975   // Check if point is on the 'pt_inside' side
0976   // Perform explicit Inside check to detect wrong side points. This impacts
0977   // Safety performance by 5-10% for all topologies
0978   Inside_t inside = ScalarInsideSegPhi(unplaced, point, zIndex, phiIndex);
0979   if (inside == EInside::kSurface) return 0.;
0980   bool contains = (inside == EInside::kInside);
0981   if (contains ^ pt_inside) return -1.;
0982
0983   // Right
0984   for (int z = zIndex; z < zMax;) {
0985     safety = Min(safety, ScalarSafetyToZSegmentSquared(unplaced, z, phiIndex, point, pt_inside, iSurf));
0986     ++z;
0987     dz = unplaced.fZPlanes[z] - point[2];
0988     // Fixed bug: dz was compared directly to safety to stop the search, while safety is a squared
0989     if (dz * dz > safety) break;
0990   }
0991   // Left
0992   for (int z = zIndex - 1; z >= 0; --z) {
0993     safety = Min(safety, ScalarSafetyToZSegmentSquared(unplaced, z, phiIndex, point, pt_inside, iSurf));
0994     dz     = point[2] - unplaced.fZPlanes[z];
0995     if (dz * dz > safety) break;
0996   }
0997
0998   // Endcap
0999   ScalarSafetyToEndcapsSquared(unplaced, point, safety, iz);
1000
1001   safety = vecCore::math::Sqrt(safety);
1002   return safety;
1003 }
1004
1005 template <Polyhedron::EInnerRadii innerRadiiT, Polyhedron::EPhiCutout phiCutoutT>
1006 VECCORE_ATT_HOST_DEVICE
1007 bool PolyhedronImplementation<innerRadiiT, phiCutoutT>::ScalarNormalKernel(UnplacedStruct_t const &unplaced,
1008                                                                            Vector3D<Precision> const &point,
1009                                                                            Vector3D<Precision> &normal)
1010 {
1011
1012   Precision safety = InfinityLength<Precision>();
1013   const int zMax   = unplaced.fZSegments.size();
1014   int zIndex       = FindZSegment<Precision>(unplaced, point[2]);
1015   if (zIndex < 0) {
1016     normal = Vector3D<Precision>(0, 0, -1);
1017     return true;
1018   }
1019
1020   if (zIndex >= zMax) {
1021     normal = Vector3D<Precision>(0, 0, 1);
1022     return true;
1023   }
1024
1025   int iSeg = zIndex;
1026   Precision dz;
1027   int iSurf    = -1;
1028   int iz       = 0;
1029   int phiIndex = FindPhiSegment<Precision>(unplaced, point);
1030
1031   // Right
1032   for (int z = zIndex; z < zMax;) {
1033     int iSurfCrt        = -1;
1034     Precision safetySeg = ScalarSafetyToZSegmentSquared(unplaced, z, phiIndex, point, true, iSurfCrt);
1035     if (safetySeg < safety) {
1036       safety = safetySeg;
1037       iSeg   = z;
1038       iSurf  = iSurfCrt;
1039     }
1040     ++z;
1041     dz = unplaced.fZPlanes[z] - point[2];
1042     if (dz * dz > safety) break;
1043   }
1044   // Left
1045   for (int z = zIndex - 1; z >= 0; --z) {
1046     int iSurfCrt        = -1;
1047     Precision safetySeg = ScalarSafetyToZSegmentSquared(unplaced, z, phiIndex, point, true, iSurfCrt);
1048     if (safetySeg < safety) {
1049       safety = safetySeg;
1050       iSeg   = z;
1051       iSurf  = iSurfCrt;
1052     }
1053     dz = point[2] - unplaced.fZPlanes[z];
1054     if (dz * dz > safety) break;
1055   }
1056
1057   // Endcap
1058   ScalarSafetyToEndcapsSquared(unplaced, point, safety, iz);
1059   if (iz != 0) {
1060     normal = Vector3D<Precision>(0, 0, iz);
1061     return true;
1062   }
1063
1064   // Retrieve the segment the point is closest to.
1065   ZSegment const &segment = unplaced.fZSegments[iSeg];
1066   if (iSurf >= 0 && iSurf < 2) {
1067     normal = segment.phi.GetNormal(iSurf);
1068   } else {
1069     if (iSurf == 2)
1070       normal = -1. * segment.inner.GetNormal(phiIndex);
1071     else
1072       normal = segment.outer.GetNormal(phiIndex);
1073   }
1074   return true;
1075 }
1076
1077 template <Polyhedron::EInnerRadii innerRadiiT, Polyhedron::EPhiCutout phiCutoutT>
1078 VECCORE_ATT_HOST_DEVICE
1079 Precision PolyhedronImplementation<innerRadiiT, phiCutoutT>::ScalarDistanceToOutKernel(
1080     UnplacedStruct_t const &unplaced, Vector3D<Precision> const &point, Vector3D<Precision> const &direction,
1081     const Precision /*stepMax*/)
1082 {
1083   // Fast exclusion if out of Z range
1084   const int zMax = unplaced.fZSegments.size();
1085   if ((point[2] < unplaced.fZPlanes[0] - kTolerance) || (point[2] > unplaced.fZPlanes[zMax] + kTolerance)) return -1.;
1086
1087   // Perform explicit Inside check to detect wrong side points. This impacts
1088   // DistanceToOut performance by about 20% for all topologies
1089   auto inside = ScalarInsideKernel(unplaced, point);
1090   if (inside == kOutside) return -1.;
1091
1092   int zIndex = FindZSegment<Precision>(unplaced, point[2]);
1093   // Don't go out of bounds
1094   zIndex = zIndex < 0 ? 0 : (zIndex >= zMax ? zMax - 1 : zIndex);
1095
1096   // Traverse Z-segments left or right depending on sign of direction
1097   bool goingRight = direction[2] >= 0;
1098
1099   Precision distance = InfinityLength<Precision>();
1100   if (goingRight) {
1101     for (; zIndex < zMax; ++zIndex) {
1102       distance = DistanceToOutZSegment<Precision>(unplaced, zIndex, unplaced.fZPlanes[zIndex],
1103                                                   unplaced.fZPlanes[zIndex + 1], point, direction);
1104       if (distance >= 0 && distance < InfinityLength<Precision>()) break;
1105       if (unplaced.fZPlanes[zIndex] - point[2] > distance) break;
1106     }
1107   } else {
1108     // Going left
1109     for (; zIndex >= 0; --zIndex) {
1110       distance = DistanceToOutZSegment<Precision>(unplaced, zIndex, unplaced.fZPlanes[zIndex],
1111                                                   unplaced.fZPlanes[zIndex + 1], point, direction);
1112       if (distance >= 0 && distance < InfinityLength<Precision>()) break;
1113       if (point[2] - unplaced.fZPlanes[zIndex] > distance) break;
1114     }
1115   }
1116
1117   // Endcaps
1118   ScalarDistanceToEndcaps<true>(unplaced, goingRight, point, direction, distance);
1119
1120   // disabling stepMax until convention revised and clear
1121   // there is a problem when distance = infinity due to some error condition but stepMax finite
1122   // return distance < stepMax ? distance : stepMax;
1123   // If not hitting anything, we must be on an edge since point is not outside
1124   if (distance >= InfinityLength<Precision>()) distance = 0.;
1125   return distance;
1126 }
1127
1128 template <Polyhedron::EInnerRadii innerRadiiT, Polyhedron::EPhiCutout phiCutoutT>
1129 template <typename Real_v, typename Bool_v>
1130 VECCORE_ATT_HOST_DEVICE
1131 void PolyhedronImplementation<innerRadiiT, phiCutoutT>::UnplacedContains(UnplacedStruct_t const &unplaced,
1132                                                                          Vector3D<Real_v> const &point, Bool_v &inside)
1133 {
1134
1135   inside = ScalarContainsKernel(unplaced, point);
1136 }
1137
1138 template <Polyhedron::EInnerRadii innerRadiiT, Polyhedron::EPhiCutout phiCutoutT>
1139 template <typename Real_v, typename Bool_v>
1140 VECGEOM_FORCE_INLINE
1141 VECCORE_ATT_HOST_DEVICE
1142 void PolyhedronImplementation<innerRadiiT, phiCutoutT>::Contains(UnplacedStruct_t const &unplaced,
1143                                                                  Vector3D<Real_v> const &point, Bool_v &inside)
1144 {
1145
1146   // we should assert if Backend != scalar
1147   inside = ScalarContainsKernel(unplaced, point);
1148 }
1149
1150 template <Polyhedron::EInnerRadii innerRadiiT, Polyhedron::EPhiCutout phiCutoutT>
1151 template <typename Real_v, typename Inside_t>
1152 VECCORE_ATT_HOST_DEVICE
1153 void PolyhedronImplementation<innerRadiiT, phiCutoutT>::Inside(UnplacedStruct_t const &unplaced,
1154                                                                Vector3D<Real_v> const &point, Inside_t &inside)
1155 {
1156
1157   // we should assert if Backend != scalar
1158   inside = ScalarInsideKernel(unplaced, point);
1159 }
1160
1161 template <Polyhedron::EInnerRadii innerRadiiT, Polyhedron::EPhiCutout phiCutoutT>
1162 template <typename Real_v>
1163 VECCORE_ATT_HOST_DEVICE
1164 void PolyhedronImplementation<innerRadiiT, phiCutoutT>::DistanceToIn(UnplacedStruct_t const &unplaced,
1165                                                                      Vector3D<Real_v> const &point,
1166                                                                      Vector3D<Real_v> const &direction,
1167                                                                      Real_v const &stepMax, Real_v &distance)
1168 {
1169   distance = ScalarDistanceToInKernel(unplaced, point, direction, stepMax);
1170 }
1171
1172 template <Polyhedron::EInnerRadii innerRadiiT, Polyhedron::EPhiCutout phiCutoutT>
1173 template <typename Real_v>
1174 VECCORE_ATT_HOST_DEVICE
1175 void PolyhedronImplementation<innerRadiiT, phiCutoutT>::DistanceToOut(UnplacedStruct_t const &unplaced,
1176                                                                       Vector3D<Real_v> const &point,
1177                                                                       Vector3D<Real_v> const &direction,
1178                                                                       Real_v const &stepMax, Real_v &distance)
1179 {
1180
1181   distance = ScalarDistanceToOutKernel(unplaced, point, direction, stepMax);
1182 }
1183
1184 template <Polyhedron::EInnerRadii innerRadiiT, Polyhedron::EPhiCutout phiCutoutT>
1185 template <typename Real_v>
1186 VECCORE_ATT_HOST_DEVICE
1187 void PolyhedronImplementation<innerRadiiT, phiCutoutT>::SafetyToIn(UnplacedStruct_t const &unplaced,
1188                                                                    Vector3D<Real_v> const &point, Real_v &safety)
1189 {
1190
1191   safety = ScalarSafetyKernel(unplaced, point, false);
1192 }
1193
1194 template <Polyhedron::EInnerRadii innerRadiiT, Polyhedron::EPhiCutout phiCutoutT>
1195 template <typename Real_v>
1196 VECGEOM_FORCE_INLINE
1197 VECCORE_ATT_HOST_DEVICE
1198 void PolyhedronImplementation<innerRadiiT, phiCutoutT>::SafetyToOut(UnplacedStruct_t const &unplaced,
1199                                                                     Vector3D<Real_v> const &point, Real_v &safety)
1200 {
1201
1202   safety = ScalarSafetyKernel(unplaced, point, true);
1203 }
1204
1205 } // namespace VECGEOM_IMPL_NAMESPACE
1206 } // namespace vecgeom
1207
1208 #endif // VECGEOM_VOLUMES_KERNEL_POLYHEDRONIMPLEMENTATION_H_