Back to home page

EIC code displayed by LXR

 
 

    


File indexing completed on 2026-09-26 09:15:49

0001 /*
0002  * PolyconeImplementation.h
0003  *
0004  *  Created on: Dec 8, 2014
0005  *      Author: swenzel
0006  */
0007 
0008 /// History notes:
0009 /// Jan-March 2017: revision + moving to use new Cone Kernels (Raman Sehgal)
0010 /// May-June 2017: revision + moving to new Structure (Raman Sehgal)
0011 
0012 #ifndef VECGEOM_VOLUMES_KERNEL_POLYCONEIMPLEMENTATION_H_
0013 #define VECGEOM_VOLUMES_KERNEL_POLYCONEIMPLEMENTATION_H_
0014 
0015 #include "VecGeom/base/Vector3D.h"
0016 #include "VecGeom/volumes/PolyconeStruct.h"
0017 #include "VecGeom/volumes/kernel/GenericKernels.h"
0018 #include <VecCore/VecCore>
0019 #include "VecGeom/volumes/kernel/ConeImplementation.h"
0020 #include "VecGeom/volumes/kernel/shapetypes/ConeTypes.h"
0021 
0022 #include <cstdio>
0023 
0024 namespace vecgeom {
0025 
0026 VECGEOM_DEVICE_DECLARE_CONV_TEMPLATE(struct, PolyconeImplementation, typename);
0027 
0028 inline namespace VECGEOM_IMPL_NAMESPACE {
0029 
0030 template <typename T>
0031 class SPlacedPolycone;
0032 template <typename T>
0033 class SUnplacedPolycone;
0034 
0035 template <typename polyconeTypeT>
0036 struct PolyconeImplementation {
0037 
0038   using UnplacedStruct_t = PolyconeStruct<Precision>;
0039   using UnplacedVolume_t = SUnplacedPolycone<polyconeTypeT>;
0040   using PlacedShape_t    = SPlacedPolycone<UnplacedVolume_t>;
0041 
0042   template <typename Real_v, bool ForInside>
0043   VECGEOM_FORCE_INLINE VECCORE_ATT_HOST_DEVICE static void GenericKernelForASection(
0044       UnplacedStruct_t const &unplaced, int isect, Vector3D<Real_v> const &polyconePoint,
0045       typename vecCore::Mask_v<Real_v> &secFullyInside, typename vecCore::Mask_v<Real_v> &secFullyOutside)
0046   {
0047 
0048     // using namespace PolyconeTypes;
0049     using namespace ConeTypes;
0050 
0051     if (isect < 0) {
0052       secFullyInside  = false;
0053       secFullyOutside = true;
0054       return;
0055     }
0056 
0057     PolyconeSection const &sec    = unplaced.GetSection(isect);
0058     Vector3D<Precision> secLocalp = polyconePoint - Vector3D<Precision>(0, 0, sec.fShift);
0059 #ifdef POLYCONEDEBUG
0060     std::cerr << " isect=" << isect << "/" << unplaced.GetNSections() << " secLocalP=" << secLocalp
0061               << ", secShift=" << sec.fShift << " sec.fSolid=" << sec.fSolid << std::endl;
0062     if (sec.fSolid) sec.fSolid->Print();
0063 #endif
0064 
0065     ConeHelpers<Real_v, polyconeTypeT>::template GenericKernelForContainsAndInside<ForInside>(
0066         sec.fSolid, secLocalp, secFullyInside, secFullyOutside);
0067   }
0068 
0069   template <typename Real_v, typename Bool_v>
0070   VECGEOM_FORCE_INLINE VECCORE_ATT_HOST_DEVICE static void Contains(UnplacedStruct_t const &polycone,
0071                                                                     Vector3D<Real_v> const &point, Bool_v &inside)
0072   {
0073 
0074     Bool_v unused(false), outside(false);
0075     GenericKernelForContainsAndInside<Real_v, Bool_v, false>(polycone, point, unused, outside);
0076     inside = !outside;
0077   }
0078 
0079   // BIG QUESTION: DO WE WANT TO GIVE ALL 3 TEMPLATE PARAMETERS
0080   // -- OR -- DO WE WANT TO DEDUCE Bool_v, Index_t from Real_v???
0081   template <typename Real_v, typename Inside_t>
0082   VECGEOM_FORCE_INLINE VECCORE_ATT_HOST_DEVICE static void Inside(UnplacedStruct_t const &polycone,
0083                                                                   Vector3D<Real_v> const &point, Inside_t &inside)
0084   {
0085 
0086     using Bool_v       = vecCore::Mask_v<Real_v>;
0087     using InsideBool_v = vecCore::Mask_v<Inside_t>;
0088     Bool_v completelyinside(false), completelyoutside(false);
0089     GenericKernelForContainsAndInside<Real_v, Bool_v, true>(polycone, point, completelyinside, completelyoutside);
0090     inside = EInside::kSurface;
0091     vecCore::MaskedAssign(inside, (InsideBool_v)completelyoutside, Inside_t(EInside::kOutside));
0092     vecCore::MaskedAssign(inside, (InsideBool_v)completelyinside, Inside_t(EInside::kInside));
0093   }
0094 
0095   template <typename Real_v, typename Bool_v, bool ForInside>
0096   VECGEOM_FORCE_INLINE VECCORE_ATT_HOST_DEVICE static void GenericKernelForContainsAndInside(
0097       UnplacedStruct_t const &unplaced, Vector3D<Real_v> const &localPoint, Bool_v &completelyInside,
0098       Bool_v &completelyOutside)
0099   {
0100 
0101     typedef Bool_v Bool_t;
0102 
0103     int indexLow  = unplaced.GetSectionIndex(localPoint.z() - kTolerance);
0104     int indexHigh = unplaced.GetSectionIndex(localPoint.z() + kTolerance);
0105     if (indexLow < 0 && indexHigh < 0) {
0106       completelyOutside = true;
0107       return;
0108     }
0109     if (indexLow < 0 && indexHigh == 0) {
0110       // Check location in section 0 and return
0111       GenericKernelForASection<Real_v, ForInside>(unplaced, 0, localPoint, completelyInside, completelyOutside);
0112       return;
0113     }
0114     if (indexHigh < 0 && indexLow == (unplaced.GetNSections() - 1)) {
0115       // Check location in section N-1 and return
0116       GenericKernelForASection<Real_v, ForInside>(unplaced, (unplaced.GetNSections() - 1), localPoint, completelyInside,
0117                                                   completelyOutside);
0118 
0119       return;
0120     }
0121     if (indexLow >= 0 && indexHigh >= 0) {
0122       if (indexLow == indexHigh) {
0123         // Check location in section indexLow and return
0124         GenericKernelForASection<Real_v, ForInside>(unplaced, indexLow, localPoint, completelyInside,
0125                                                     completelyOutside);
0126 
0127         return;
0128       } else {
0129 
0130         Bool_t secInLow = false, secOutLow = false;
0131         Bool_t secInHigh = false, secOutHigh = false;
0132         GenericKernelForASection<Real_v, ForInside>(unplaced, indexLow, localPoint, secInLow, secOutLow);
0133         GenericKernelForASection<Real_v, ForInside>(unplaced, indexHigh, localPoint, secInHigh, secOutHigh);
0134         Bool_t surfLow  = !secInLow && !secOutLow;
0135         Bool_t surfHigh = !secInHigh && !secOutHigh;
0136 
0137         if (surfLow && surfHigh) {
0138           completelyInside = true;
0139           return;
0140         } else {
0141           // else if point is on surface of only one of the two sections then point is actually on surface , the default
0142           // case,
0143           // so no need to check
0144 
0145           // What needs to check is if it is outside both ie. Outside indexLow section and Outside indexHigh section
0146           // then it is certainly outside
0147           if (secOutLow && secOutHigh) {
0148             completelyOutside = true;
0149             return;
0150           }
0151         }
0152       }
0153     }
0154   }
0155 
0156   template <typename Real_v>
0157   VECGEOM_FORCE_INLINE VECCORE_ATT_HOST_DEVICE static void DistanceToIn(UnplacedStruct_t const &polycone,
0158                                                                         Vector3D<Real_v> const &point,
0159                                                                         Vector3D<Real_v> const &direction,
0160                                                                         Real_v const &stepMax, Real_v &distance)
0161   {
0162     // using namespace PolyconeTypes;
0163     Vector3D<Real_v> p = point;
0164     Vector3D<Real_v> v = direction;
0165 
0166 #ifdef POLYCONEDEBUG
0167     std::cerr << "Polycone::DistToIn() (spot 1): point=" << point << ", dir=" << direction << ", localPoint=" << p
0168               << ", localDir=" << v << "\n";
0169 #endif
0170 
0171     // TODO: add bounding box check maybe??
0172 
0173     distance      = kInfLength;
0174     int increment = (v.z() > 0) ? 1 : -1;
0175     if (std::fabs(v.z()) < kTolerance) increment = 0;
0176     int index = polycone.GetSectionIndex(p.z());
0177     if (index == -1) index = 0;
0178     if (index == -2) index = polycone.GetNSections() - 1;
0179 
0180     do {
0181       // now we have to find a section
0182       PolyconeSection const &sec = polycone.GetSection(index);
0183 
0184 #ifdef POLYCONEDEBUG
0185       std::cerr << "Polycone::DistToIn() (spot 2):"
0186                 << " index=" << index << " NSec=" << polycone.GetNSections() << " &sec=" << &sec << " - secPars:"
0187                 << " secOffset=" << sec.fShift << " Dz=" << sec.fSolid.GetDz() << " Rmin1=" << sec.fSolid.GetRmin1()
0188                 << " Rmin2=" << sec.fSolid.GetRmin2() << " Rmax1=" << sec.fSolid.GetRmax1()
0189                 << " Rmax2=" << sec.fSolid.GetRmax2() << " -- calling Cone::DistToIn()...\n";
0190 #endif
0191 
0192       ConeImplementation<polyconeTypeT>::template DistanceToIn<Real_v>(
0193           sec.fSolid, p - Vector3D<Precision>(0, 0, sec.fShift), v, stepMax, distance);
0194 
0195 #ifdef POLYCONEDEBUG
0196       std::cerr << "Polycone::DistToIn() (spot 3):"
0197                 << " distToIn() = " << distance << "\n";
0198 #endif
0199 
0200       if (distance < kInfLength || !increment) break;
0201       index += increment;
0202     } while (index >= 0 && index < polycone.GetNSections());
0203     return;
0204   }
0205 
0206   template <typename Real_v>
0207   VECGEOM_FORCE_INLINE VECCORE_ATT_HOST_DEVICE static void DistanceToOut(UnplacedStruct_t const &polycone,
0208                                                                          Vector3D<Real_v> const &point,
0209                                                                          Vector3D<Real_v> const &dir,
0210                                                                          Real_v const &stepMax, Real_v &distance)
0211   {
0212     // using namespace PolyconeTypes;
0213     distance            = kInfLength;
0214     Vector3D<Real_v> pn = point;
0215 
0216     // specialization for N==1??? It should be a cone in the first place
0217     if (polycone.GetNSections() == 1) {
0218       const PolyconeSection &section = polycone.GetSection(0);
0219 
0220       ConeImplementation<polyconeTypeT>::template DistanceToOut<Real_v>(
0221           section.fSolid, point - Vector3D<Precision>(0, 0, section.fShift), dir, stepMax, distance);
0222 
0223       return;
0224     }
0225 
0226     int indexLow  = polycone.GetSectionIndex(point.z() - kTolerance);
0227     int indexHigh = polycone.GetSectionIndex(point.z() + kTolerance);
0228     int index     = 0;
0229 
0230     // section index is -1 when out of left-end
0231     // section index is -2 when beyond right-end
0232 
0233     if (indexLow < 0 && indexHigh < 0) {
0234       distance = -1;
0235       return;
0236     } else if (indexLow < 0 && indexHigh >= 0) {
0237       index                          = indexHigh;
0238       const PolyconeSection &section = polycone.GetSection(index);
0239 
0240       Inside_t inside;
0241       //      ConeImplementation<ConeTypes::UniversalCone>::Inside<Real_v>(
0242       //          section.fSolid, point - Vector3D<Precision>(0, 0, section.fShift), inside);
0243       ConeImplementation<polyconeTypeT>::template Inside<Real_v>(
0244           section.fSolid, point - Vector3D<Precision>(0, 0, section.fShift), inside);
0245       if (inside == EInside::kOutside) {
0246         distance = -1;
0247         return;
0248       }
0249     } else if (indexLow != indexHigh && (indexLow >= 0)) {
0250       // we are close to an intermediate Surface, section has to be identified
0251       const PolyconeSection &section = polycone.GetSection(indexLow);
0252 
0253       Inside_t inside;
0254       //      ConeImplementation<ConeTypes::UniversalCone>::Inside<Real_v>(
0255       //        section.fSolid, point - Vector3D<Precision>(0, 0, section.fShift), inside);
0256 
0257       ConeImplementation<polyconeTypeT>::template Inside<Real_v>(
0258           section.fSolid, point - Vector3D<Precision>(0, 0, section.fShift), inside);
0259 
0260       if (inside == EInside::kOutside) {
0261         index = indexHigh;
0262       } else {
0263         index = indexLow;
0264       }
0265     } else {
0266       index = indexLow;
0267       if (index < 0) index = polycone.GetSectionIndex(point.z());
0268     }
0269     if (index < 0) {
0270       distance = 0.;
0271       return;
0272     }
0273     // Added
0274     else {
0275       const PolyconeSection &section = polycone.GetSection(index);
0276 
0277       Inside_t inside;
0278       //      ConeImplementation<ConeTypes::UniversalCone>::Inside<Real_v>(
0279       //          section.fSolid, point - Vector3D<Precision>(0, 0, section.fShift), inside);
0280       ConeImplementation<polyconeTypeT>::template Inside<Real_v>(
0281           section.fSolid, point - Vector3D<Precision>(0, 0, section.fShift), inside);
0282       if (inside == EInside::kOutside) {
0283         distance = -1;
0284         return;
0285       }
0286     }
0287 
0288     Precision totalDistance = 0.;
0289     Precision dist;
0290     int increment = (dir.z() > 0) ? 1 : -1;
0291     if (std::fabs(dir.z()) < kTolerance) increment = 0;
0292 
0293     // What is the relevance of istep?
0294     int istep = 0;
0295     do {
0296       const PolyconeSection &section = polycone.GetSection(index);
0297 
0298       if ((totalDistance != 0) || (istep < 2)) {
0299         pn = point + totalDistance * dir; // point must be shifted, so it could eventually get into another solid
0300         pn.z() -= section.fShift;
0301         Inside_t inside;
0302         //        ConeImplementation<ConeTypes::UniversalCone>::Inside<Real_v>(section.fSolid, pn, inside);
0303         ConeImplementation<polyconeTypeT>::template Inside<Real_v>(section.fSolid, pn, inside);
0304 
0305         if (inside == EInside::kOutside) {
0306           break;
0307         }
0308       } else
0309         pn.z() -= section.fShift;
0310 
0311       istep++;
0312 
0313       // ConeImplementation<ConeTypes::UniversalCone>::DistanceToOut<Real_v>(section.fSolid, pn, dir, stepMax, dist);
0314       ConeImplementation<polyconeTypeT>::template DistanceToOut<Real_v>(section.fSolid, pn, dir, stepMax, dist);
0315       if (dist == -1) return;
0316 
0317       // Section Surface case
0318       if (std::fabs(dist) < 0.5 * kTolerance) {
0319         int index1 = index;
0320         if ((index > 0) && (index < polycone.GetNSections() - 1)) {
0321           index1 += increment;
0322         } else {
0323           if ((index == 0) && (increment > 0)) index1 += increment;
0324           if ((index == polycone.GetNSections() - 1) && (increment < 0)) index1 += increment;
0325         }
0326 
0327         Vector3D<Precision> pte         = point + (totalDistance + dist) * dir;
0328         const PolyconeSection &section1 = polycone.GetSection(index1);
0329         pte.z() -= section1.fShift;
0330         Vector3D<Precision> localp;
0331         Inside_t inside22;
0332         // ConeImplementation<ConeTypes::UniversalCone>::Inside<Real_v>(section1.fSolid, pte, inside22);
0333         ConeImplementation<polyconeTypeT>::template Inside<Real_v>(section1.fSolid, pte, inside22);
0334         if (inside22 == 3 || (increment == 0)) {
0335           break;
0336         }
0337       } // end if surface case
0338 
0339       totalDistance += dist;
0340       index += increment;
0341     } while (increment != 0 && index >= 0 && index < polycone.GetNSections());
0342 
0343     distance = totalDistance;
0344 
0345     return;
0346   }
0347 
0348   template <typename Real_v>
0349   VECGEOM_FORCE_INLINE VECCORE_ATT_HOST_DEVICE static void SafetyToIn(UnplacedStruct_t const &polycone,
0350                                                                       Vector3D<Real_v> const &point, Real_v &safety)
0351   {
0352 
0353     Vector3D<Real_v> p = point;
0354     int index          = polycone.GetSectionIndex(p.z());
0355 
0356     bool needZ = false;
0357     if (index < 0) {
0358       needZ = true;
0359       if (index == -1) index = 0;
0360       if (index == -2) index = polycone.GetNSections() - 1;
0361     }
0362     Precision minSafety        = 0; //= SafetyFromOutsideSection(index, p);
0363     PolyconeSection const &sec = polycone.GetSection(index);
0364     // safety to current segment
0365     if (needZ) {
0366       //      ConeImplementation<ConeTypes::UniversalCone>::SafetyToIn<Real_v>(
0367       //          sec.fSolid, p - Vector3D<Precision>(0, 0, sec.fShift), safety);
0368       ConeImplementation<polyconeTypeT>::template SafetyToIn<Real_v>(sec.fSolid,
0369                                                                      p - Vector3D<Precision>(0, 0, sec.fShift), safety);
0370     } else {
0371 
0372       //      ConeImplementation<ConeTypes::UniversalCone>::SafetyToIn<Real_v>(
0373       //          sec.fSolid, p - Vector3D<Precision>(0, 0, sec.fShift), safety);
0374       ConeImplementation<polyconeTypeT>::template SafetyToIn<Real_v>(sec.fSolid,
0375                                                                      p - Vector3D<Precision>(0, 0, sec.fShift), safety);
0376 
0377       if (safety < kTolerance) return;
0378       minSafety       = safety;
0379       Precision zbase = polycone.fZs[index + 1];
0380       // going right
0381       for (int i = index + 1; i < polycone.GetNSections(); ++i) {
0382         Precision dz = polycone.fZs[i] - zbase;
0383         if (dz >= minSafety) break;
0384 
0385         PolyconeSection const &sect = polycone.GetSection(i);
0386 
0387         //      ConeImplementation<ConeTypes::UniversalCone>::SafetyToIn<Real_v>(
0388         //          sect.fSolid, p - Vector3D<Precision>(0, 0, sect.fShift), safety);
0389 
0390         ConeImplementation<polyconeTypeT>::template SafetyToIn<Real_v>(
0391             sect.fSolid, p - Vector3D<Precision>(0, 0, sect.fShift), safety);
0392 
0393         if (safety < minSafety) minSafety = safety;
0394       }
0395 
0396       // going left if this is possible
0397       if (index > 0) {
0398         zbase = polycone.fZs[index - 1];
0399         for (int i = index - 1; i >= 0; --i) {
0400           Precision dz = zbase - polycone.fZs[i];
0401           if (dz >= minSafety) break;
0402           PolyconeSection const &sect = polycone.GetSection(i);
0403 
0404           //        ConeImplementation<ConeTypes::UniversalCone>::SafetyToIn<Real_v>(
0405           //            sect.fSolid, p - Vector3D<Precision>(0, 0, sect.fShift), safety);
0406 
0407           ConeImplementation<polyconeTypeT>::template SafetyToIn<Real_v>(
0408               sect.fSolid, p - Vector3D<Precision>(0, 0, sect.fShift), safety);
0409 
0410           if (safety < minSafety) minSafety = safety;
0411         }
0412       }
0413       safety = minSafety;
0414     }
0415     return;
0416   }
0417   /* A function to calculation the shortest distance of a point from a segment */
0418   template <typename Real_v>
0419   VECGEOM_FORCE_INLINE VECCORE_ATT_HOST_DEVICE static Real_v DistanceToSeg(Vector3D<Real_v> const &point,
0420                                                                            Vector3D<Precision> segment_start,
0421                                                                            Vector3D<Precision> segment_end)
0422   {
0423     Vector3D<Real_v> segment_vector   = segment_end - segment_start;
0424     Vector3D<Real_v> point_vector     = point - segment_start;
0425     Real_v projection_scalar          = point_vector.Dot(segment_vector) / segment_vector.Mag2();
0426     Vector3D<Real_v> projection_point = segment_start + projection_scalar * segment_vector;
0427     vecCore__MaskedAssignFunc(projection_point, projection_scalar < Real_v(0.), Vector3D<Real_v>(segment_start));
0428     vecCore__MaskedAssignFunc(projection_point, projection_scalar > Real_v(1.), Vector3D<Real_v>(segment_end));
0429     Vector3D<Real_v> distVec = point - projection_point;
0430     return distVec.Mag();
0431   }
0432 
0433   /*
0434   ** New definition that uses Segments and Single Wedge, instead of SafetyToOut from ConeImplementation
0435   */
0436   template <typename Real_v>
0437   VECGEOM_FORCE_INLINE VECCORE_ATT_HOST_DEVICE static void SafetyToOut(UnplacedStruct_t const &polycone,
0438                                                                        Vector3D<Real_v> const &point, Real_v &safety)
0439   {
0440     typedef typename vecCore::Mask_v<Real_v> Bool_v;
0441     safety = Real_v(kInfLength);
0442     Bool_v compIn(false), compOut(false), done(false);
0443     GenericKernelForContainsAndInside<Real_v, Bool_v, true>(polycone, point, compIn, compOut);
0444     vecCore__MaskedAssignFunc(safety, compOut && !done, Real_v(-1.));
0445     done |= compOut;
0446     if (vecCore::MaskFull(done)) return;
0447 
0448     Vector3D<Real_v> twoDPoint = Vector3D<Real_v>(point.Perp(), point.z(), 0.);
0449     for (unsigned int currSegIndex = 0; currSegIndex < polycone.fTwoDVec.size(); currSegIndex++) {
0450       unsigned int nextSegIndex = currSegIndex + 1;
0451       if (currSegIndex == polycone.fTwoDVec.size() - 1) {
0452         nextSegIndex = 0;
0453       }
0454 
0455       Real_v distance(kInfLength);
0456       if ((polycone.fTwoDVec[currSegIndex].x() == 0 && polycone.fTwoDVec[nextSegIndex].x() == 0) ||
0457           (polycone.fTwoDVec[currSegIndex].x() == polycone.fTwoDVec[nextSegIndex].x() &&
0458            polycone.fTwoDVec[currSegIndex].y() == polycone.fTwoDVec[nextSegIndex].y()))
0459         distance = Real_v(kInfLength);
0460       else
0461         distance = DistanceToSeg(twoDPoint, polycone.fTwoDVec[currSegIndex], polycone.fTwoDVec[nextSegIndex]);
0462 
0463       vecCore__MaskedAssignFunc(safety, (distance < safety) && !done, distance);
0464     }
0465 
0466     if (polycone.fDeltaPhi < 2 * kPi) {
0467       Real_v safetyPhi = polycone.fPhiWedge.SafetyToOut<Real_v>(point);
0468       vecCore__MaskedAssignFunc(safety, safetyPhi < safety, safetyPhi);
0469     }
0470   }
0471 
0472 #if (0)
0473   /* Retaining the old definition for the time being */
0474   template <typename Real_v>
0475   VECGEOM_FORCE_INLINE VECCORE_ATT_HOST_DEVICE static void SafetyToOut_Old(UnplacedStruct_t const &polycone,
0476                                                                            Vector3D<Real_v> const &point,
0477                                                                            Real_v &safety)
0478   {
0479     typedef typename vecCore::Mask_v<Real_v> Bool_v;
0480     Bool_v compIn(false), compOut(false);
0481     GenericKernelForContainsAndInside<Real_v, Bool_v, true>(polycone, point, compIn, compOut);
0482     if (compOut) {
0483       safety = -1;
0484       return;
0485     }
0486 
0487     int index = polycone.GetSectionIndex(point.z());
0488     if (index < 0) {
0489       safety = -1;
0490       return;
0491     }
0492 
0493     PolyconeSection const &sec = polycone.GetSection(index);
0494 
0495     Vector3D<Real_v> p = point - Vector3D<Precision>(0, 0, sec.fShift);
0496     // ConeImplementation<ConeTypes::UniversalCone>::SafetyToOut<Real_v>(sec.fSolid, p, safety);
0497     ConeImplementation<polyconeTypeT>::template SafetyToOut<Real_v>(sec.fSolid, p, safety);
0498 
0499     Precision minSafety = safety;
0500     if (minSafety == kInfLength) {
0501       safety = 0.;
0502       return;
0503     }
0504     if (minSafety < kTolerance) {
0505       safety = 0.;
0506       return;
0507     }
0508 
0509     Precision zbase = polycone.fZs[index + 1];
0510     for (int i = index + 1; i < polycone.GetNSections(); ++i) {
0511       Precision dz = polycone.fZs[i] - zbase;
0512       if (dz >= minSafety) break;
0513       PolyconeSection const &sect = polycone.GetSection(i);
0514       p                           = point - Vector3D<Precision>(0, 0, sect.fShift);
0515 
0516       // ConeImplementation<ConeTypes::UniversalCone>::SafetyToIn<Real_v>(sect.fSolid, p, safety);
0517       ConeImplementation<polyconeTypeT>::template SafetyToIn<Real_v>(sect.fSolid, p, safety);
0518 
0519       if (safety < minSafety) minSafety = safety;
0520     }
0521 
0522     if (index > 0) {
0523       zbase = polycone.fZs[index - 1];
0524       for (int i = index - 1; i >= 0; --i) {
0525         Precision dz = zbase - polycone.fZs[i];
0526         if (dz >= minSafety) break;
0527         PolyconeSection const &sect = polycone.GetSection(i);
0528         p                           = point - Vector3D<Precision>(0, 0, sect.fShift);
0529 
0530         // ConeImplementation<ConeTypes::UniversalCone>::SafetyToIn<Real_v>(sect.fSolid, p, safety);
0531         ConeImplementation<polyconeTypeT>::template SafetyToIn<Real_v>(sect.fSolid, p, safety);
0532 
0533         if (safety < minSafety) minSafety = safety;
0534       }
0535     }
0536 
0537     safety = minSafety;
0538     return;
0539   }
0540 #endif
0541 };
0542 } // namespace VECGEOM_IMPL_NAMESPACE
0543 } // namespace vecgeom
0544 
0545 #endif // VECGEOM_VOLUMES_KERNEL_polyconeIMPLEMENTATION_H_