Back to home page

EIC code displayed by LXR
 
 

    

Warning, file /include/VecGeom/volumes/kernel/ConeImplementation.h was not indexed or was modified since last indexation (in which case cross-reference links may be missing, inaccurate or erroneous).

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

This page was automatically generated by the 2.3.7 LXR engine.

The LXR team
Valid CSS 2.1! Valid HTML 4.01!