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