volumes/kernel/ConeImplementation.h

0001 /*
0002  * ConeImplementation.h
0003  *
0004  *  Created on: May 14, 2014
0005  *      Author: swenzel
0006  */
0007
0008 /// History notes:
0009 /// revision + moving to Vectorized Cone Kernels (Raman Sehgal)
0010 /// May-June 2017: revision + moving to new Structure (Raman Sehgal)
0011 /// 20180323 Guilherme Lima  Adapted to new UnplacedVolume factory
0012
0013 #ifndef VECGEOM_VOLUMES_KERNEL_CONEIMPLEMENTATION_H_
0014 #define VECGEOM_VOLUMES_KERNEL_CONEIMPLEMENTATION_H_
0015
0016 #include "VecGeom/base/Vector3D.h"
0017 #include "VecGeom/volumes/kernel/GenericKernels.h"
0018 #include "VecGeom/volumes/kernel/shapetypes/ConeTypes.h"
0019 #include "VecGeom/volumes/ConeStruct.h"
0020 #include <cstdio>
0021 #include "VecGeom/volumes/ConeUtilities.h"
0022
0023 namespace vecgeom {
0024
0025 VECGEOM_DEVICE_DECLARE_CONV_TEMPLATE(struct, ConeImplementation, typename);
0026
0027 inline namespace VECGEOM_IMPL_NAMESPACE {
0028
0029 template <typename T>
0030 class SPlacedCone;
0031 template <typename T>
0032 class SUnplacedCone;
0033
0034 template <typename coneTypeT>
0035 struct ConeImplementation {
0036
0037   using UnplacedStruct_t = ConeStruct<Precision>;
0038   using UnplacedVolume_t = SUnplacedCone<coneTypeT>;
0039   using PlacedShape_t    = SPlacedCone<UnplacedVolume_t>;
0040
0041   VECCORE_ATT_HOST_DEVICE
0042   static void PrintType() {}
0043
0044   template <typename Stream>
0045   static void PrintType(Stream &s, int transCodeT = translation::kGeneric, int rotCodeT = rotation::kGeneric)
0046   {
0047     s << "SpecializedCone<" << transCodeT << "," << rotCodeT << ">";
0048   }
0049
0050   template <typename Stream>
0051   static void PrintImplementationType(Stream & /*s*/)
0052   {
0053   }
0054
0055   template <typename Stream>
0056   static void PrintUnplacedType(Stream & /*s*/)
0057   {
0058   }
0059
0060   /* A Function that will just check if the point is on the CONICAL (circle) edge
0061    * assuming that it is on either lowerZ or upperZ
0062    *
0063    * Beware : It will not do any checks on Z
0064    */
0065   template <typename Real_v, bool ForInnerSurface, bool ForLowerZ>
0066   VECGEOM_FORCE_INLINE
0067   VECCORE_ATT_HOST_DEVICE
0068   static typename vecCore::Mask_v<Real_v> IsOnRing(UnplacedStruct_t const &cone, Vector3D<Real_v> const &point)
0069   {
0070     using Bool_v = typename vecCore::Mask_v<Real_v>;
0071
0072     Real_v rad2 = point.Perp2();
0073     Bool_v onRing(false);
0074
0075     if (ForLowerZ) {
0076       if (ForInnerSurface) {
0077         onRing = (rad2 <= MakePlusTolerantSquare<true>(cone.fRmin1)) &&
0078                  (rad2 >= MakeMinusTolerantSquare<true>(cone.fRmin1));
0079       } else {
0080     onRing = (rad2 <= MakePlusTolerantSquare<true>(cone.fRmax1)) &&
0081                  (rad2 >= MakeMinusTolerantSquare<true>(cone.fRmax1));
0082       }
0083     } else {
0084       if (ForInnerSurface) {
0085         onRing = (rad2 <= MakePlusTolerantSquare<true>(cone.fRmin2)) &&
0086                  (rad2 >= MakeMinusTolerantSquare<true>(cone.fRmin2));
0087       } else {
0088         onRing = (rad2 <= MakePlusTolerantSquare<true>(cone.fRmax2)) &&
0089                  (rad2 >= MakeMinusTolerantSquare<true>(cone.fRmax2));
0090       }
0091     }
0092
0093     return onRing;
0094   }
0095
0096   template <typename Real_v>
0097   VECGEOM_FORCE_INLINE
0098   VECCORE_ATT_HOST_DEVICE
0099   static void Contains(UnplacedStruct_t const &cone, Vector3D<Real_v> const &point,
0100                        typename vecCore::Mask_v<Real_v> &inside)
0101   {
0102     typedef typename vecCore::Mask_v<Real_v> Bool_v;
0103     Bool_v unused(false);
0104     Bool_v outside(false);
0105     ConeHelpers<Real_v, coneTypeT>::template GenericKernelForContainsAndInside<false>(cone, point, unused, outside);
0106     inside = !outside;
0107   }
0108
0109   template <typename Real_v, typename Inside_v>
0110   VECGEOM_FORCE_INLINE
0111   VECCORE_ATT_HOST_DEVICE
0112   static void Inside(UnplacedStruct_t const &cone, Vector3D<Real_v> const &point, Inside_v &inside)
0113   {
0114     ConeHelpers<Real_v, coneTypeT>::template Inside<Inside_v>(cone, point, inside);
0115   }
0116
0117   template <typename Real_v>
0118   VECGEOM_FORCE_INLINE
0119   VECCORE_ATT_HOST_DEVICE
0120   static void DistanceToIn(UnplacedStruct_t const &cone, Vector3D<Real_v> const &point, Vector3D<Real_v> const &dir,
0121                            Real_v const & /*stepMax*/, Real_v &distance)
0122   {
0123     using namespace ConeUtilities;
0124     using namespace ConeTypes;
0125     typedef Real_v Float_t;
0126     typedef typename vecCore::Mask_v<Real_v> Bool_t;
0127
0128     Bool_t done(false);
0129     const Real_v zero(0.);
0130
0131     //=== First, for points outside and moving away --> return infinity
0132     distance = kInfLength;
0133
0134     // outside of Z range and going away?
0135     Float_t distz          = Abs(point.z()) - cone.fDz; // avoid a division for now
0136     Bool_t outZAndGoingOut = (distz > kConeTolerance && (point.z() * dir.z()) >= zero) ||
0137                              (Abs(distz) < kConeTolerance && (point.z() * dir.z()) > zero);
0138     done |= outZAndGoingOut;
0139     if (vecCore::MaskFull(done)) return;
0140
0141     // outside or *on* outer cone and going away?
0142     Float_t outerRad  = GetRadiusOfConeAtPoint<Real_v, false>(cone, point.z());
0143     Float_t rsq            = point.Perp2(); // point.x()*point.x() + point.y()*point.y();
0144     done |= (rsq > MakeMinusTolerantSquare<true>(outerRad, cone.fOuterTolerance))
0145          && (dir.Dot(GetNormal<Real_v, false>(cone, point)) >= zero);
0146     if (vecCore::MaskFull(done)) return;
0147
0148     //=== Next, check all dimensions of the cone: for points inside --> return -1
0149     vecCore__MaskedAssignFunc(distance, !done, Float_t(-1.0));
0150
0151     // For points inside z-range, return -1
0152     Bool_t inside = distz < -kConeTolerance;
0153
0154     inside &= rsq < MakeMinusTolerantSquare<true>(outerRad, cone.fOuterTolerance);
0155
0156     if (checkRminTreatment<coneTypeT>(cone)) {
0157       Float_t innerRad = GetRadiusOfConeAtPoint<Real_v, true>(cone, point.z());
0158       inside &= rsq > MakePlusTolerantSquare<true>(innerRad, cone.fInnerTolerance);
0159     }
0160     if (checkPhiTreatment<coneTypeT>(cone)) { // && !vecCore::MaskEmpty(inside)) {
0161       Bool_t insector;
0162       PointInCyclicalSector<Real_v, coneTypeT, false, false>(cone, point.x(), point.y(), insector);
0163       inside &= insector;
0164     }
0165     done |= inside;
0166     if (vecCore::MaskFull(done)) return;
0167
0168     //=== Next step: check if z-plane is the right entry point (both r,phi
0169     // should be valid at z-plane crossing)
0170     vecCore__MaskedAssignFunc(distance, !done, Float_t(kInfLength));
0171
0172     distz /= NonZero(Abs(dir.z()));
0173
0174
0175 #ifdef EDGE_POINTS
0176     Bool_t onZsurf  = (Abs(point.z()) - cone.fDz) < Real_v(kConeTolerance);
0177     Bool_t onLoZSrf = onZsurf && point.z() < zero;
0178     Bool_t onHiZSrf = onZsurf && point.z() > zero;
0179     Bool_t loZcond  = onLoZSrf && (IsOnRing<Real_v, false, true>(cone, point));
0180     Bool_t hiZcond  = onHiZSrf && (IsOnRing<Real_v, false,false>(cone, point));
0181     if (checkRminTreatment<coneTypeT>(cone)) {
0182        loZcond |= onLoZSrf && IsOnRing<Real_v, true, true>(cone, point);
0183        hiZcond |= onHiZSrf && IsOnRing<Real_v, true,false>(cone, point);
0184     }
0185     vecCore::MaskedAssign(distz, loZcond || hiZcond, zero);
0186 #endif
0187
0188
0189     Float_t hitx = point.x() + distz * dir.x();
0190     Float_t hity = point.y() + distz * dir.y();
0191
0192     Float_t r2 = (hitx * hitx) + (hity * hity);
0193
0194     Precision innerZTol         = cone.fTolIz;
0195     Bool_t isHittingTopPlane    = (point.z() >= innerZTol) && (r2 <= cone.fSqRmax2 + kTolerance);  // cone.fSqRmax2Tol
0196     Bool_t isHittingBottomPlane = (point.z() <= -innerZTol) && (r2 <= cone.fSqRmax1 + kTolerance); // cone.fSqRmax1Tol
0197     Bool_t okz                  = (isHittingTopPlane || isHittingBottomPlane);
0198
0199     if (checkRminTreatment<coneTypeT>(cone)) {
0200       isHittingTopPlane &= (r2 >= cone.fSqRmin2 - kTolerance);    // cone.fSqRmin2Tol
0201       isHittingBottomPlane &= (r2 >= cone.fSqRmin1 - kTolerance); // cone.fSqRmin1Tol
0202       okz &= ((isHittingTopPlane || isHittingBottomPlane));
0203     }
0204
0205     if (checkPhiTreatment<coneTypeT>(cone)) {
0206       Bool_t insector;
0207       PointInCyclicalSector<Real_v, coneTypeT, false>(cone, hitx, hity, insector);
0208       okz &= insector;
0209     }
0210     vecCore::MaskedAssign(distance, !done && okz, distz);
0211     done |= okz;
0212     if (vecCore::MaskFull(done)) return;
0213
0214     Float_t dist_rOuter(kInfLength);
0215     Bool_t ok_outerCone =
0216         ConeHelpers<Real_v, coneTypeT>::template DetectIntersectionAndCalculateDistanceToConicalSurface<true, false>(
0217             cone, point, dir, dist_rOuter);
0218     ok_outerCone &= dist_rOuter < distance;
0219     vecCore::MaskedAssign(distance, !done && ok_outerCone, dist_rOuter);
0220     done |= ok_outerCone;
0221     if (vecCore::MaskFull(done)) return;
0222
0223     Float_t dist_rInner(kInfLength);
0224     if (checkRminTreatment<coneTypeT>(cone)) {
0225
0226       Bool_t ok_innerCone =
0227           ConeHelpers<Real_v, coneTypeT>::template DetectIntersectionAndCalculateDistanceToConicalSurface<true, true>(
0228               cone, point, dir, dist_rInner);
0229       ok_innerCone &= dist_rInner < distance;
0230       vecCore::MaskedAssign(distance, !done && ok_innerCone, dist_rInner);
0231     }
0232
0233     if (checkPhiTreatment<coneTypeT>(cone)) {
0234
0235       evolution::Wedge const &w = cone.fPhiWedge;
0236
0237       Float_t dist_phi;
0238       Bool_t ok_phi;
0239       PhiPlaneTrajectoryIntersection<Real_v, coneTypeT, SectorType<coneTypeT>::value != kOnePi, true>(
0240           cone.fAlongPhi1x, cone.fAlongPhi1y, w.GetNormal1().x(), w.GetNormal1().y(), cone, point, dir, dist_phi,
0241           ok_phi);
0242       ok_phi &= dist_phi < distance;
0243       vecCore::MaskedAssign(distance, !done && ok_phi, dist_phi);
0244       done |= ok_phi;
0245
0246       if (SectorType<coneTypeT>::value != kOnePi) {
0247
0248         PhiPlaneTrajectoryIntersection<Real_v, coneTypeT, true, true>(cone.fAlongPhi2x, cone.fAlongPhi2y,
0249                                                                       w.GetNormal2().x(), w.GetNormal2().y(), cone,
0250                                                                       point, dir, dist_phi, ok_phi);
0251
0252         vecCore::MaskedAssign(distance, ok_phi && dist_phi < distance, dist_phi);
0253       }
0254     }
0255   }
0256
0257   template <typename Real_v>
0258   VECGEOM_FORCE_INLINE
0259   VECCORE_ATT_HOST_DEVICE
0260   static void DistanceToOut(UnplacedStruct_t const &cone, Vector3D<Real_v> const &point,
0261                             Vector3D<Real_v> const &direction, Real_v const & /*stepMax*/, Real_v &distance)
0262   {
0263
0264     distance = kInfLength;
0265     using namespace ConeUtilities;
0266     using namespace ConeTypes;
0267
0268     typedef typename vecCore::Mask_v<Real_v> Bool_t;
0269
0270     Bool_t done(false);
0271     const Real_v zero(0.0);
0272
0273     // Using this logic will improve performance of Scalar code
0274     Real_v distz         = Abs(point.z()) - cone.fDz;
0275
0276     //=== Next, check all dimensions of the cone: for points outside --> return -1
0277     vecCore__MaskedAssignFunc(distance, !done, Real_v(-1.0));
0278
0279     Bool_t outside = distz > Real_v(kConeTolerance);
0280
0281     Real_v rsq           = point.Perp2();
0282     Real_v outerRad = ConeUtilities::GetRadiusOfConeAtPoint<Real_v, false>(cone, point.z());
0283     outside |=  rsq > MakePlusTolerantSquare<true>(outerRad, cone.fOuterTolerance);
0284     done |= outside;
0285     if (vecCore::MaskFull(done)) return;
0286
0287     if (checkRminTreatment<coneTypeT>(cone) && !vecCore::MaskFull(outside)) {
0288       Real_v innerRad = ConeUtilities::GetRadiusOfConeAtPoint<Real_v, true>(cone, point.z());
0289       outside |= rsq < MakeMinusTolerantSquare<true>(innerRad, cone.fInnerTolerance);
0290       done |= outside;
0291       if (vecCore::MaskFull(done)) return;
0292     }
0293     if (checkPhiTreatment<coneTypeT>(cone) && !vecCore::MaskEmpty(outside)) {
0294       Bool_t insector;
0295       ConeUtilities::PointInCyclicalSector<Real_v, coneTypeT, false, false>(cone, point.x(), point.y(), insector);
0296       outside |= !insector;
0297     }
0298     done |= outside;
0299     if (vecCore::MaskFull(done)) return;
0300
0301
0302     Bool_t isGoingUp   = direction.z() > zero;
0303     Bool_t isGoingDown = direction.z() < zero;
0304     Bool_t isOnZPlaneAndMovingOutside(false);
0305     isOnZPlaneAndMovingOutside = !outside && ((isGoingUp && point.z() > zero && Abs(distz) < kConeTolerance) ||
0306                                               (isGoingDown && point.z() < zero && Abs(distz) < kConeTolerance));
0307     vecCore__MaskedAssignFunc(distance, !done && isOnZPlaneAndMovingOutside, distz);
0308     done |= isOnZPlaneAndMovingOutside;
0309     if (vecCore::MaskFull(done)) return;
0310
0311
0312     //=== Next step: check if z-plane is the right entry point (both r,phi
0313     // should be valid at z-plane crossing)
0314     vecCore__MaskedAssignFunc(distance, !done, Real_v(kInfLength));
0315
0316     Precision fDz  = cone.fDz;
0317     Real_v dirZInv = Real_v(1.) / NonZero(direction.z());
0318     vecCore__MaskedAssignFunc(distance, isGoingUp, (fDz - point.z()) * dirZInv);
0319     vecCore__MaskedAssignFunc(distance, isGoingDown, (-fDz - point.z()) * dirZInv);
0320
0321     Real_v dist_rOuter(kInfLength);
0322     Bool_t ok_outerCone =
0323         ConeHelpers<Real_v, coneTypeT>::template DetectIntersectionAndCalculateDistanceToConicalSurface<false, false>(
0324             cone, point, direction, dist_rOuter);
0325
0326     vecCore::MaskedAssign(distance, !done && ok_outerCone && dist_rOuter < distance, dist_rOuter);
0327
0328     Real_v dist_rInner(kInfLength);
0329     if (checkRminTreatment<coneTypeT>(cone)) {
0330       Bool_t ok_innerCone =
0331           ConeHelpers<Real_v, coneTypeT>::template DetectIntersectionAndCalculateDistanceToConicalSurface<false, true>(
0332               cone, point, direction, dist_rInner);
0333       vecCore::MaskedAssign(distance, !done && ok_innerCone && dist_rInner < distance, dist_rInner);
0334     }
0335
0336     if (checkPhiTreatment<coneTypeT>(cone)) {
0337
0338       Bool_t isOnStartPhi      = ConeUtilities::IsOnStartPhi<Real_v>(cone, point);
0339       Bool_t isOnEndPhi        = ConeUtilities::IsOnEndPhi<Real_v>(cone, point);
0340       Vector3D<Real_v> normal1 = cone.fPhiWedge.GetNormal1();
0341       Vector3D<Real_v> normal2 = cone.fPhiWedge.GetNormal2();
0342       Bool_t cond              = (isOnStartPhi && direction.Dot(-normal1) > zero) ||
0343                     (isOnEndPhi && direction.Dot(-normal2) > zero);
0344       vecCore__MaskedAssignFunc(distance, !done && cond, zero);
0345       done |= cond;
0346       if (vecCore::MaskFull(done)) return;
0347
0348       Real_v dist_phi;
0349       Bool_t ok_phi;
0350       evolution::Wedge const &w = cone.fPhiWedge;
0351       PhiPlaneTrajectoryIntersection<Real_v, coneTypeT, SectorType<coneTypeT>::value != kOnePi, false>(
0352           cone.fAlongPhi1x, cone.fAlongPhi1y, w.GetNormal1().x(), w.GetNormal1().y(), cone, point, direction, dist_phi,
0353           ok_phi);
0354       ok_phi &= dist_phi < distance;
0355       vecCore::MaskedAssign(distance, !done && ok_phi, dist_phi);
0356       done |= ok_phi;
0357
0358       if (SectorType<coneTypeT>::value != kOnePi) {
0359         ConeUtilities::PhiPlaneTrajectoryIntersection<Real_v, coneTypeT, true, false>(
0360             cone.fAlongPhi2x, cone.fAlongPhi2y, w.GetNormal2().x(), w.GetNormal2().y(), cone, point, direction,
0361             dist_phi, ok_phi);
0362         vecCore::MaskedAssign(distance, ok_phi && dist_phi < distance, dist_phi);
0363       }
0364     }
0365     vecCore__MaskedAssignFunc(distance, distance < zero && Abs(distance) < kTolerance, zero);
0366   }
0367
0368   template <typename Real_v>
0369   VECGEOM_FORCE_INLINE
0370   VECCORE_ATT_HOST_DEVICE
0371   static void SafetyToIn(UnplacedStruct_t const &cone, Vector3D<Real_v> const &point, Real_v &safety)
0372   {
0373     using namespace ConeUtilities;
0374     using namespace ConeTypes;
0375     safety = -kInfLength;
0376     typedef typename vecCore::Mask_v<Real_v> Bool_t;
0377     typedef Real_v Float_t;
0378
0379     Bool_t done(false);
0380     Precision fDz = cone.fDz;
0381     Float_t distz = Abs(point.z()) - fDz;
0382
0383     // Next, check all dimensions of the cone, whether points are inside -->
0384     // return -1
0385     vecCore__MaskedAssignFunc(safety, !done, Float_t(-1.0));
0386
0387     // For points inside z-range, return -1
0388     Bool_t inside = distz < -kConeTolerance;
0389
0390     // This logic to check if the point is inside is far better than
0391     // using GenericKernel and will improve performance.
0392     Float_t outerRad = GetRadiusOfConeAtPoint<Real_v, false>(cone, point.z());
0393     Float_t rsq      = point.Perp2();
0394     inside &= rsq < MakeMinusTolerantSquare<true>(outerRad, cone.fOuterTolerance);
0395
0396     if (checkRminTreatment<coneTypeT>(cone)) {
0397       Float_t innerRad = GetRadiusOfConeAtPoint<Real_v, true>(cone, point.z());
0398       inside &= rsq > MakePlusTolerantSquare<true>(innerRad, cone.fInnerTolerance);
0399     }
0400     if (checkPhiTreatment<coneTypeT>(cone) && !vecCore::MaskEmpty(inside)) {
0401       Bool_t insector;
0402       PointInCyclicalSector<Real_v, coneTypeT, false, false>(cone, point.x(), point.y(), insector);
0403       inside &= insector;
0404     }
0405     done |= inside;
0406     if (vecCore::MaskFull(done)) return;
0407
0408     // Once it is checked that the point is inside or not, safety can be set to 0.
0409     // This will serve the case that the point is on the surface. So no need to check
0410     // that the point is really on surface.
0411     vecCore__MaskedAssignFunc(safety, !done, Float_t(0.));
0412
0413     // Now if the point is neither inside nor on surface, then it should be outside
0414     // and the safety should be set to some finite value, which is done by below logic
0415
0416     Float_t safeZ                = Abs(point.z()) - fDz;
0417     Float_t safeDistOuterSurface = -SafeDistanceToConicalSurface<Real_v, false>(cone, point);
0418
0419     Float_t safeDistInnerSurface(-kInfLength);
0420     if (checkRminTreatment<coneTypeT>(cone)) {
0421       safeDistInnerSurface = -SafeDistanceToConicalSurface<Real_v, true>(cone, point);
0422     }
0423
0424     vecCore__MaskedAssignFunc(safety, !done, Max(safeZ, Max(safeDistOuterSurface, safeDistInnerSurface)));
0425
0426     if (checkPhiTreatment<coneTypeT>(cone)) {
0427       Float_t safetyPhi = cone.fPhiWedge.SafetyToIn<Real_v>(point);
0428       vecCore__MaskedAssignFunc(safety, !done, Max(safetyPhi, safety));
0429     }
0430
0431     vecCore__MaskedAssignFunc(safety, vecCore::math::Abs(safety) < kTolerance, Float_t(0.));
0432   }
0433
0434   template <typename Real_v>
0435   VECGEOM_FORCE_INLINE
0436   VECCORE_ATT_HOST_DEVICE
0437   static void SafetyToOut(UnplacedStruct_t const &cone, Vector3D<Real_v> const &point, Real_v &safety)
0438   {
0439
0440     using namespace ConeUtilities;
0441     using namespace ConeTypes;
0442     safety = kInfLength;
0443     typedef typename vecCore::Mask_v<Real_v> Bool_t;
0444     typedef Real_v Float_t;
0445
0446     Bool_t done(false);
0447
0448     Float_t distz = Abs(point.z()) - cone.fDz;
0449     Float_t rsq   = point.Perp2();
0450
0451     //=== Next, check all dimensions of the cone --> for points outside z-range, return -1
0452     vecCore__MaskedAssignFunc(safety, !done, Float_t(-1.0));
0453
0454     // This logic to check if the point is outside is far better then
0455     // using GenericKernel and will improve performance.
0456     Bool_t outside = distz > Real_v(kConeTolerance);
0457
0458     Float_t outerRad  = GetRadiusOfConeAtPoint<Real_v, false>(cone, point.z());
0459     outside |= rsq > MakePlusTolerantSquare<true>(outerRad, cone.fOuterTolerance);
0460
0461     if (checkRminTreatment<coneTypeT>(cone)) {
0462       Float_t innerRad = GetRadiusOfConeAtPoint<Real_v, true>(cone, point.z());
0463       outside |= rsq < MakeMinusTolerantSquare<true>(innerRad, cone.fInnerTolerance);
0464     }
0465
0466     if (checkPhiTreatment<coneTypeT>(cone) && !vecCore::MaskEmpty(outside)) {
0467
0468       Bool_t insector;
0469       ConeUtilities::PointInCyclicalSector<Real_v, coneTypeT, false, false>(cone, point.x(), point.y(), insector);
0470       outside |= !insector;
0471     }
0472     done |= outside;
0473     if (vecCore::MaskFull(done)) return;
0474
0475     // Once it is checked that the point is inside or not, safety can be set to 0.
0476     // This will serve the case that the point is on the surface. So no need to check
0477     // that the point is really on surface.
0478     vecCore__MaskedAssignFunc(safety, !done, Float_t(0.));
0479
0480     // Now if the point is neither outside nor on surface, then it should be inside
0481     // and the safety should be set to some finite value, which is done by below logic
0482
0483     Precision fDz = cone.fDz;
0484     Float_t safeZ = fDz - Abs(point.z());
0485
0486     Float_t safeDistOuterSurface = SafeDistanceToConicalSurface<Real_v, false>(cone, point);
0487     Float_t safeDistInnerSurface(kInfLength);
0488     if (checkRminTreatment<coneTypeT>(cone)) {
0489       safeDistInnerSurface = SafeDistanceToConicalSurface<Real_v, true>(cone, point);
0490     }
0491
0492     vecCore__MaskedAssignFunc(safety, !done, Min(safeZ, Min(safeDistOuterSurface, safeDistInnerSurface)));
0493
0494     if (checkPhiTreatment<coneTypeT>(cone)) {
0495       Float_t safetyPhi = cone.fPhiWedge.SafetyToOut<Real_v>(point);
0496       vecCore__MaskedAssignFunc(safety, !done, Min(safetyPhi, safety));
0497     }
0498     vecCore__MaskedAssignFunc(safety, vecCore::math::Abs(safety) < kTolerance, Float_t(0.));
0499   }
0500
0501   template <typename Real_v, bool ForInnerSurface>
0502   VECGEOM_FORCE_INLINE
0503   VECCORE_ATT_HOST_DEVICE
0504   static Real_v SafeDistanceToConicalSurface(UnplacedStruct_t const &cone, Vector3D<Real_v> const &point)
0505   {
0506
0507     typedef Real_v Float_t;
0508     Float_t rho = point.Perp();
0509     if (ForInnerSurface) {
0510       Float_t pRMin = cone.fTanRMin * point.z() + (cone.fRmin1 + cone.fRmin2) * Float_t(0.5); // cone.fRminAv;
0511       return (rho - pRMin) * cone.fInvSecRMin;
0512     } else {
0513       Float_t pRMax = cone.fTanRMax * point.z() + (cone.fRmax1 + cone.fRmax2) * Float_t(0.5); // cone.fRmaxAv;
0514       return (pRMax - rho) * cone.fInvSecRMax;
0515     }
0516   }
0517 };
0518 } // namespace VECGEOM_IMPL_NAMESPACE
0519 } // namespace vecgeom
0520
0521 #endif