VecGeom/volumes/ConeUtilities.h

0001 /*
0002  * ConeUtilities.h
0003  *
0004  *  Created on: June 01, 2017
0005  *      Author: Raman Sehgal
0006  */
0007 #ifndef VECGEOM_CONEUTILITIES_H_
0008 #define VECGEOM_CONEUTILITIES_H_
0009
0010 #include "VecGeom/base/Global.h"
0011 #include "VecGeom/volumes/Wedge_Evolution.h"
0012 #include "VecGeom/base/Vector3D.h"
0013 #include "VecGeom/volumes/ConeStruct.h"
0014 #include "VecGeom/volumes/kernel/GenericKernels.h"
0015 #include "VecGeom/volumes/kernel/shapetypes/ConeTypes.h"
0016 #include "VecGeom/volumes/kernel/TubeImplementation.h"
0017 #include <cstdio>
0018
0019 namespace vecgeom {
0020
0021 inline namespace VECGEOM_IMPL_NAMESPACE {
0022
0023 class UnplacedCone;
0024 template <typename T>
0025 struct ConeStruct;
0026 using UnplacedStruct_t = ConeStruct<Precision>;
0027
0028 namespace ConeUtilities {
0029
0030 /**
0031  * Returns whether a point is inside a cylindrical sector, as defined
0032  * by the two vectors that go along the endpoints of the sector
0033  *
0034  * The same could be achieved using atan2 to calculate the angle formed
0035  * by the point, the origin and the X-axes, but this is a lot faster,
0036  * using only multiplications and comparisons
0037  *
0038  * (-x*starty + y*startx) >= 0: calculates whether going from the start vector to the point
0039  * we are traveling in the CCW direction (taking the shortest direction, of course)
0040  *
0041  * (-endx*y + endy*x) >= 0: calculates whether going from the point to the end vector
0042  * we are traveling in the CCW direction (taking the shortest direction, of course)
0043  *
0044  * For a sector smaller than pi, we need that BOTH of them hold true - if going from start, to the
0045  * point, and then to the end we are travelling in CCW, it's obvious the point is inside the
0046  * cylindrical sector.
0047  *
0048  * For a sector bigger than pi, only one of the conditions needs to be true. This is less obvious why.
0049  * Since the sector angle is greater than pi, it can be that one of the two vectors might be
0050  * farther than pi away from the point. In that case, the shortest direction will be CW, so even
0051  * if the point is inside, only one of the two conditions need to hold.
0052  *
0053  * If going from start to point is CCW, then certainly the point is inside as the sector
0054  * is larger than pi.
0055  *
0056  * If going from point to end is CCW, again, the point is certainly inside.
0057  *
0058  * This function is a frankensteinian creature that can determine which of the two cases (smaller vs
0059  * larger than pi) to use either at compile time (if it has enough information, saving an ifVolumeType
0060  * statement) or at runtime.
0061  **/
0062
0063 #if (1)
0064 template <typename Real_v, typename ShapeType, bool onSurfaceT, bool includeSurface = true>
0065 VECGEOM_FORCE_INLINE
0066 VECCORE_ATT_HOST_DEVICE
0067 static void PointInCyclicalSector(UnplacedStruct_t const &volume, Real_v const &x, Real_v const &y,
0068                                   typename vecCore::Mask_v<Real_v> &ret)
0069 {
0070
0071   using namespace ::vecgeom::ConeTypes;
0072   // assert(SectorType<ShapeType>::value != kNoAngle && "ShapeType without a
0073   // sector passed to PointInCyclicalSector");
0074
0075   // typedef Real_v Real_v;
0076   // using vecgeom::ConeTypes::SectorType;
0077   // using vecgeom::ConeTypes::EAngleType;
0078
0079   Real_v startx = volume.fAlongPhi1x; // GetAlongPhi1X();
0080   Real_v starty = volume.fAlongPhi1y; // GetAlongPhi1Y();
0081
0082   Real_v endx = volume.fAlongPhi2x; // GetAlongPhi2X();
0083   Real_v endy = volume.fAlongPhi2y; // GetAlongPhi2Y();
0084
0085   bool smallerthanpi;
0086
0087   if (SectorType<ShapeType>::value == kUnknownAngle)
0088     smallerthanpi = volume.fDPhi <= M_PI;
0089   else
0090     smallerthanpi = SectorType<ShapeType>::value == kOnePi || SectorType<ShapeType>::value == kSmallerThanPi;
0091
0092   Real_v startCheck = (-x * starty) + (y * startx);
0093   Real_v endCheck   = (-endx * y) + (endy * x);
0094
0095   if (onSurfaceT) {
0096     // in this case, includeSurface is irrelevant
0097     ret = (Abs(startCheck) <= kConeTolerance) || (Abs(endCheck) <= kConeTolerance);
0098   } else {
0099     if (smallerthanpi) {
0100       if (includeSurface)
0101         ret = (startCheck >= -kConeTolerance) & (endCheck >= -kConeTolerance);
0102       else
0103         ret = (startCheck >= kConeTolerance) & (endCheck >= kConeTolerance);
0104     } else {
0105       if (includeSurface)
0106         ret = (startCheck >= -kConeTolerance) || (endCheck >= -kConeTolerance);
0107       else
0108         ret = (startCheck >= kConeTolerance) || (endCheck >= kConeTolerance);
0109     }
0110   }
0111 }
0112
0113 #endif
0114
0115 #if (1)
0116 template <typename Real_v, bool ForInnerRadius>
0117 VECGEOM_FORCE_INLINE
0118 VECCORE_ATT_HOST_DEVICE
0119 static Real_v GetRadiusOfConeAtPoint(UnplacedStruct_t const &cone, Real_v const pointZ)
0120 {
0121
0122   if (ForInnerRadius) {
0123     if (cone.fRmin1 == cone.fRmin2) {
0124       return cone.fRmin1;
0125     } else {
0126       return cone.fInnerSlope * pointZ + cone.fInnerOffset;
0127     }
0128
0129   } else {
0130     if (cone.fOriginalRmax1 == cone.fOriginalRmax2) {
0131       return cone.fOriginalRmax1;
0132     } else {
0133       return cone.fOuterSlope * pointZ + cone.fOuterOffset;
0134     }
0135   }
0136 }
0137
0138 #endif
0139
0140 /*
0141  * Check intersection of the trajectory with a phi-plane
0142  * All points of the along-vector of a phi plane lie on
0143  * s * (alongX, alongY)
0144  * All points of the trajectory of the particle lie on
0145  * (x, y) + t * (vx, vy)
0146  * Thefore, it must hold that s * (alongX, alongY) == (x, y) + t * (vx, vy)
0147  * Solving by t we get t = (alongY*x - alongX*y) / (vy*alongX - vx*alongY)
0148  * s = (x + t*vx) / alongX = (newx) / alongX
0149  *
0150  * If we have two non colinear phi-planes, need to make sure
0151  * point falls on its positive direction <=> dot product between
0152  * along vector and hit-point is positive <=> hitx*alongX + hity*alongY > 0
0153  */
0154
0155 template <typename Real_v, typename ConeType, bool PositiveDirectionOfPhiVector, bool insectorCheck>
0156 VECGEOM_FORCE_INLINE
0157 VECCORE_ATT_HOST_DEVICE
0158 static void PhiPlaneTrajectoryIntersection(Precision alongX, Precision alongY, Precision normalX, Precision normalY,
0159                                            UnplacedStruct_t const &cone, Vector3D<Real_v> const &pos,
0160                                            Vector3D<Real_v> const &dir, Real_v &dist,
0161                                            typename vecCore::Mask_v<Real_v> &ok)
0162 {
0163   const Real_v zero(0.0);
0164   dist = kInfLength;
0165
0166   // approaching phi plane from the right side?
0167   // this depends whether we use it for DistanceToIn or DistanceToOut
0168   // Note: wedge normals poing towards the wedge inside, by convention!
0169   if (insectorCheck)
0170     ok = ((dir.x() * normalX) + (dir.y() * normalY) > zero); // DistToIn  -- require tracks entering volume
0171   else
0172     ok = ((dir.x() * normalX) + (dir.y() * normalY) < zero); // DistToOut -- require tracks leaving volume
0173
0174   // if( /*Backend::early_returns &&*/ vecCore::MaskEmpty(ok) ) return;
0175
0176   Real_v dirDotXY = (dir.y() * alongX) - (dir.x() * alongY);
0177   vecCore__MaskedAssignFunc(dist, dirDotXY != 0, ((alongY * pos.x()) - (alongX * pos.y())) / NonZero(dirDotXY));
0178   ok &= dist > -kConeTolerance;
0179   // if( /*Backend::early_returns &&*/ vecCore::MaskEmpty(ok) ) return;
0180
0181   if (insectorCheck) {
0182     Real_v hitx          = pos.x() + dist * dir.x();
0183     Real_v hity          = pos.y() + dist * dir.y();
0184     Real_v hitz          = pos.z() + dist * dir.z();
0185     Real_v r2            = (hitx * hitx) + (hity * hity);
0186     Real_v innerRadIrTol = GetRadiusOfConeAtPoint<Real_v, true>(cone, hitz) + kTolerance;
0187     Real_v outerRadIrTol = GetRadiusOfConeAtPoint<Real_v, false>(cone, hitz) - kTolerance;
0188
0189     ok &= Abs(hitz) <= cone.fTolIz && (r2 >= innerRadIrTol * innerRadIrTol) && (r2 <= outerRadIrTol * outerRadIrTol);
0190
0191     // GL: tested with this if(PosDirPhiVec) around if(insector), so
0192     // if(insector){} requires PosDirPhiVec==true to run
0193     //  --> shapeTester still finishes OK (no mismatches) (some cycles saved...)
0194     if (PositiveDirectionOfPhiVector) {
0195       ok = ok && ((hitx * alongX) + (hity * alongY)) > zero;
0196     }
0197   } else {
0198     if (PositiveDirectionOfPhiVector) {
0199       Real_v hitx = pos.x() + dist * dir.x();
0200       Real_v hity = pos.y() + dist * dir.y();
0201       ok          = ok && ((hitx * alongX) + (hity * alongY)) >= zero;
0202     }
0203   }
0204 }
0205
0206 template <typename Real_v, bool ForInnerSurface>
0207 VECCORE_ATT_HOST_DEVICE
0208 static Vector3D<Real_v> GetNormal(UnplacedStruct_t const &cone, Vector3D<Real_v> const &point)
0209 {
0210
0211   // typedef Real_v Real_v;
0212   Real_v rho = point.Perp();
0213   Vector3D<Real_v> norm(0., 0., 0.);
0214
0215   if (ForInnerSurface) {
0216     // Handling inner conical surface
0217     Precision rmin1 = cone.fRmin1;
0218     Precision rmin2 = cone.fRmin2;
0219     if ((rmin1 == rmin2) && (rmin1 != 0.)) {
0220       // cone act like tube
0221       norm.Set(-point.x(), -point.y(), 0.);
0222     } else {
0223       Precision secRMin = cone.fSecRMin;
0224       norm.Set(-point.x(), -point.y(), cone.fZNormInner * (rho * secRMin));
0225     }
0226   } else {
0227     Precision rmax1 = cone.fRmax1;
0228     Precision rmax2 = cone.fRmax2;
0229     if ((rmax1 == rmax2) && (rmax1 != 0.)) {
0230       // cone act like tube
0231       norm.Set(point.x(), point.y(), 0.);
0232     } else {
0233       Precision secRMax = cone.fSecRMax;
0234       norm.Set(point.x(), point.y(), cone.fZNormOuter * (rho * secRMax));
0235     }
0236   }
0237   return norm;
0238 }
0239
0240 template <typename Real_v, bool ForInnerSurface>
0241 VECGEOM_FORCE_INLINE
0242 VECCORE_ATT_HOST_DEVICE
0243 static typename vecCore::Mask_v<Real_v> IsOnConicalSurface(UnplacedStruct_t const &cone, Vector3D<Real_v> const &point)
0244 {
0245
0246   using namespace ConeUtilities;
0247   using namespace ConeTypes;
0248   const Real_v rho       = point.Perp2();
0249   const Real_v coneRad   = GetRadiusOfConeAtPoint<Real_v, ForInnerSurface>(cone, point.z());
0250   const Real_v coneRad2  = coneRad * coneRad;
0251   const Real_v tolerance = (ForInnerSurface) ? cone.fInnerTolerance : cone.fOuterTolerance;
0252   return (rho >= (coneRad2 - tolerance * coneRad)) && (rho <= (coneRad2 + tolerance * coneRad)) &&
0253          (Abs(point.z()) < (cone.fDz + kConeTolerance));
0254 }
0255
0256 // precondition: point is on cone surface - as returned from IsOnConicalSurface()
0257 template <typename Real_v, bool ForInnerSurface>
0258 VECGEOM_FORCE_INLINE
0259 VECCORE_ATT_HOST_DEVICE
0260 static typename vecCore::Mask_v<Real_v> IsMovingOutsideConicalSurface(UnplacedStruct_t const &cone,
0261                                                                       Vector3D<Real_v> const &point,
0262                                                                       Vector3D<Real_v> const &direction)
0263 {
0264   return direction.Dot(GetNormal<Real_v, ForInnerSurface>(cone, point)) >= Real_v(0.);
0265 }
0266
0267 // precondition: point is on cone surface - as returned from IsOnConicalSurface()
0268 template <typename Real_v, bool ForInnerSurface>
0269 VECGEOM_FORCE_INLINE
0270 VECCORE_ATT_HOST_DEVICE
0271 static typename vecCore::Mask_v<Real_v> IsMovingInsideConicalSurface(UnplacedStruct_t const &cone,
0272                                                                      Vector3D<Real_v> const &point,
0273                                                                      Vector3D<Real_v> const &direction)
0274 {
0275   return direction.Dot(GetNormal<Real_v, ForInnerSurface>(cone, point)) <= Real_v(0.);
0276 }
0277
0278 template <typename Real_v>
0279 VECGEOM_FORCE_INLINE
0280 VECCORE_ATT_HOST_DEVICE
0281 static typename vecCore::Mask_v<Real_v> IsOnStartPhi(UnplacedStruct_t const &cone, Vector3D<Real_v> const &point)
0282 {
0283   //  class evolution::Wedge;
0284   return cone.fPhiWedge.IsOnSurfaceGeneric(cone.fPhiWedge.GetAlong1(), cone.fPhiWedge.GetNormal1(), point);
0285 }
0286
0287 template <typename Real_v>
0288 VECGEOM_FORCE_INLINE
0289 VECCORE_ATT_HOST_DEVICE
0290 static typename vecCore::Mask_v<Real_v> IsOnEndPhi(UnplacedStruct_t const &cone, Vector3D<Real_v> const &point)
0291 {
0292
0293   return cone.fPhiWedge.IsOnSurfaceGeneric(cone.fPhiWedge.GetAlong2(), cone.fPhiWedge.GetNormal2(), point);
0294 }
0295
0296 template <typename Real_v, bool ForTopPlane>
0297 VECCORE_ATT_HOST_DEVICE
0298 static typename vecCore::Mask_v<Real_v> IsOnZPlaneAndMovingInside(UnplacedStruct_t const &cone,
0299                                                                   Vector3D<Real_v> const &point,
0300                                                                   Vector3D<Real_v> const &direction)
0301 {
0302
0303   Real_v rho    = point.Perp2();
0304   Precision fDz = cone.fDz;
0305
0306   if (ForTopPlane) {
0307     return (rho > (cone.fSqRmin2 - kConeTolerance)) && (rho < (cone.fSqRmax2 + kConeTolerance)) &&
0308            (point.z() < (fDz + kConeTolerance)) && (point.z() > (fDz - kConeTolerance)) &&
0309            (direction.z() < Real_v(0.));
0310   } else {
0311     return (rho > (cone.fSqRmin1 - kConeTolerance)) && (rho < (cone.fSqRmax1 + kConeTolerance)) &&
0312            (point.z() < (-fDz + kConeTolerance)) && (point.z() > (-fDz - kConeTolerance)) &&
0313            (direction.z() > Real_v(0.));
0314   }
0315 }
0316
0317 template <typename Real_v, bool ForTopPlane>
0318 VECCORE_ATT_HOST_DEVICE
0319 static typename vecCore::Mask_v<Real_v> IsOnZPlaneAndMovingOutside(UnplacedStruct_t const &cone,
0320                                                                    Vector3D<Real_v> const &point,
0321                                                                    Vector3D<Real_v> const &direction)
0322 {
0323
0324   Real_v rho    = point.Perp2();
0325   Precision fDz = cone.fDz;
0326
0327   if (ForTopPlane) {
0328     return (rho > (cone.fSqRmin2 - kConeTolerance)) && (rho < (cone.fSqRmax2 + kConeTolerance)) &&
0329            (point.z() < (fDz + kConeTolerance)) && (point.z() > (fDz - kConeTolerance)) &&
0330            (direction.z() > Real_v(0.));
0331   } else {
0332     return (rho > (cone.fSqRmin1 - kConeTolerance)) && (rho < (cone.fSqRmax1 + kConeTolerance)) &&
0333            (point.z() < (-fDz + kConeTolerance)) && (point.z() > (-fDz - kConeTolerance)) &&
0334            (direction.z() < Real_v(0.));
0335   }
0336 }
0337
0338 } // namespace ConeUtilities
0339
0340 /* This class is introduced to allow Partial Specialization of selected functions,
0341 ** and will be very much useful when running Cone and Polycone in Scalar mode
0342 */
0343 template <class Real_v, class coneTypeT>
0344 class ConeHelpers {
0345
0346 public:
0347   ConeHelpers() {}
0348   ~ConeHelpers() {}
0349   template <bool ForDistToIn, bool ForInnerSurface>
0350   VECCORE_ATT_HOST_DEVICE
0351   static typename vecCore::Mask_v<Real_v> DetectIntersectionAndCalculateDistanceToConicalSurface(
0352       UnplacedStruct_t const &cone, Vector3D<Real_v> const &point, Vector3D<Real_v> const &direction, Real_v &distance)
0353   {
0354
0355     using namespace ConeUtilities;
0356     using namespace ConeTypes;
0357     typedef typename vecCore::Mask_v<Real_v> Bool_t;
0358     const Real_v zero(0.0);
0359
0360     distance                = kInfLength;
0361     Bool_t onConicalSurface = IsOnConicalSurface<Real_v, ForInnerSurface>(cone, point);
0362     Vector3D<Real_v> normal = ConeUtilities::GetNormal<Real_v, ForInnerSurface>(cone, point);
0363     Bool_t tangentToSurface = vecCore::math::Abs(direction.Dot(normal)) == zero;
0364     Bool_t done             = onConicalSurface && tangentToSurface;
0365     if (vecCore::MaskFull(done)) return Bool_t(false);
0366
0367     // if precond is all false, save some CPU
0368     const Bool_t precond = !done && onConicalSurface;
0369     if (!vecCore::MaskEmpty(precond)) {
0370       if (ForDistToIn) {
0371         Bool_t isOnSurfaceAndMovingInside = precond
0372           & ConeUtilities::IsMovingInsideConicalSurface<Real_v, ForInnerSurface>(cone, point, direction);
0373
0374         if (!checkPhiTreatment<coneTypeT>(cone)) {
0375           vecCore__MaskedAssignFunc(distance, isOnSurfaceAndMovingInside, zero);
0376           done |= isOnSurfaceAndMovingInside;
0377           if (vecCore::MaskFull(done)) return done;
0378         } else {
0379           Bool_t insector(false);
0380           ConeUtilities::PointInCyclicalSector<Real_v, coneTypeT, false, true>(cone, point.x(), point.y(), insector);
0381           vecCore__MaskedAssignFunc(distance, insector && isOnSurfaceAndMovingInside, zero);
0382           done |= (insector && isOnSurfaceAndMovingInside);
0383           if (vecCore::MaskFull(done)) return done;
0384         }
0385
0386       } else {
0387     Bool_t isOnSurfaceAndMovingOutside = precond
0388           & ConeUtilities::IsMovingOutsideConicalSurface<Real_v, ForInnerSurface>(cone, point, direction);
0389
0390         if (!checkPhiTreatment<coneTypeT>(cone)) {
0391           vecCore__MaskedAssignFunc(distance, isOnSurfaceAndMovingOutside, zero);
0392           done |= isOnSurfaceAndMovingOutside;
0393           if (vecCore::MaskFull(done)) return done;
0394         } else {
0395           Bool_t insector(false);
0396           ConeUtilities::PointInCyclicalSector<Real_v, coneTypeT, false, true>(cone, point.x(), point.y(), insector);
0397           vecCore__MaskedAssignFunc(distance, insector && isOnSurfaceAndMovingOutside, zero);
0398           done |= (insector && isOnSurfaceAndMovingOutside);
0399           if (vecCore::MaskFull(done)) return done;
0400         }
0401       }
0402     }
0403
0404     Real_v pDotV2D = point.x() * direction.x() + point.y() * direction.y();
0405
0406     Real_v a(0.), b(0.), c(0.);
0407     Bool_t ok(false);
0408     Precision fDz = cone.fDz;
0409     if (ForInnerSurface) {
0410
0411       Precision rmin1 = cone.fRmin1;
0412       Precision rmin2 = cone.fRmin2;
0413       if (rmin1 == rmin2) {
0414         b = pDotV2D;
0415         a = direction.Perp2();
0416         c = point.Perp2() - rmin2 * rmin2;
0417       } else {
0418
0419         Precision t = cone.fTanInnerApexAngle;
0420         Real_v newPz(0.);
0421         if (cone.fRmin2 > cone.fRmin1)
0422           newPz = (point.z() + fDz + cone.fInnerConeApex) * t;
0423         else
0424           newPz = (point.z() - fDz - cone.fInnerConeApex) * t;
0425
0426         Real_v dirT = direction.z() * t;
0427         a           = (direction.x() * direction.x()) + (direction.y() * direction.y()) - dirT * dirT;
0428
0429         b = pDotV2D - (newPz * dirT);
0430         c = point.Perp2() - (newPz * newPz);
0431       }
0432
0433       Real_v b2 = b * b;
0434       Real_v ac = a * c;
0435       if (vecCore::MaskFull(b2 < ac)) return Bool_t(false);
0436       Real_v d2 = b2 - ac;
0437
0438       Real_v delta = Sqrt(vecCore::math::Abs(d2));
0439       if (ForDistToIn) {
0440         vecCore__MaskedAssignFunc(distance, !done && d2 >= zero && (b >= zero), (c / NonZero(-b - delta)));
0441         vecCore__MaskedAssignFunc(distance, !done && d2 >= zero && (b < zero), (-b + delta) / NonZero(a));
0442       } else {
0443         vecCore__MaskedAssignFunc(distance, !done && d2 >= zero && (b >= zero), (-b - delta) / NonZero(a));
0444         vecCore__MaskedAssignFunc(distance, !done && d2 >= zero && (b < zero), (c / NonZero(-b + delta)));
0445       }
0446
0447       if (vecCore::MaskFull(distance < zero)) return Bool_t(false);
0448       Real_v newZ = point.z() + (direction.z() * distance);
0449       ok          = (Abs(newZ) < fDz);
0450
0451     } else {
0452
0453       // if (rmax1 == rmax2) {
0454       if (cone.fOriginalRmax1 == cone.fOriginalRmax2) {
0455         b = pDotV2D;
0456         a = direction.Perp2();
0457         c = point.Perp2() - cone.fOriginalRmax2 * cone.fOriginalRmax2;
0458       } else {
0459
0460         Precision t = cone.fTanOuterApexAngle;
0461         Real_v newPz(0.);
0462         // if (cone.fRmax2 > cone.fRmax1)
0463         if (cone.fOriginalRmax2 > cone.fOriginalRmax1)
0464           newPz = (point.z() + fDz + cone.fOuterConeApex) * t;
0465         else
0466           newPz = (point.z() - fDz - cone.fOuterConeApex) * t;
0467         Real_v dirT = direction.z() * t;
0468         a           = direction.x() * direction.x() + direction.y() * direction.y() - dirT * dirT;
0469         b           = pDotV2D - (newPz * dirT);
0470         c           = point.Perp2() - (newPz * newPz);
0471       }
0472       Real_v b2 = b * b;
0473       Real_v ac = a * c;
0474       if (vecCore::MaskFull(b2 < ac)) return Bool_t(false);
0475       Real_v d2    = b2 - ac;
0476       Real_v delta = Sqrt(vecCore::math::Abs(d2));
0477
0478       if (ForDistToIn) {
0479         vecCore__MaskedAssignFunc(distance, !done && d2 >= zero && (b > zero), (-b - delta) / NonZero(a));
0480         vecCore__MaskedAssignFunc(distance, !done && d2 >= zero && (b < zero), (c / NonZero(-b + delta)));
0481       } else {
0482         vecCore__MaskedAssignFunc(distance, !done && d2 >= zero && (b < zero), (-b + delta) / NonZero(a));
0483         vecCore__MaskedAssignFunc(distance, !done && d2 >= zero && (b >= zero), (c / NonZero(-b - delta)));
0484         ok = distance > zero;
0485       }
0486
0487       if (vecCore::MaskFull(distance < zero)) return Bool_t(false);
0488       if (ForDistToIn) {
0489         Real_v newZ = point.z() + (direction.z() * distance);
0490         ok          = (Abs(newZ) < cone.fDz + kHalfTolerance);
0491       }
0492     }
0493     vecCore__MaskedAssignFunc(distance, distance < zero, Real_v(kInfLength));
0494
0495     if (checkPhiTreatment<coneTypeT>(cone)) {
0496       Real_v hitx(0), hity(0), hitz(0);
0497       Bool_t insector(false); // = Backend::kFalse;
0498       vecCore__MaskedAssignFunc(hitx, distance < kInfLength, point.x() + distance * direction.x());
0499       vecCore__MaskedAssignFunc(hity, distance < kInfLength, point.y() + distance * direction.y());
0500       vecCore__MaskedAssignFunc(hitz, distance < kInfLength, point.z() + distance * direction.z());
0501
0502       ConeUtilities::PointInCyclicalSector<Real_v, coneTypeT, false, true>(cone, hitx, hity, insector);
0503       ok &= ((insector) && (distance < kInfLength));
0504     }
0505     return ok;
0506   }
0507
0508   template <bool ForInside>
0509   VECGEOM_FORCE_INLINE
0510   VECCORE_ATT_HOST_DEVICE
0511   static void GenericKernelForContainsAndInside(UnplacedStruct_t const &cone, Vector3D<Real_v> const &point,
0512                                                 typename vecCore::Mask_v<Real_v> &completelyinside,
0513                                                 typename vecCore::Mask_v<Real_v> &completelyoutside)
0514   {
0515
0516     typedef typename vecCore::Mask_v<Real_v> Bool_t;
0517
0518     // very fast check on z-height
0519     Real_v absz       = Abs(point[2]);
0520     completelyoutside = absz > MakePlusTolerant<ForInside>(cone.fDz, kConeTolerance);
0521     if (ForInside) {
0522       completelyinside = absz < MakeMinusTolerant<ForInside>(cone.fDz, kConeTolerance);
0523     }
0524     if (vecCore::MaskFull(completelyoutside)) {
0525       return;
0526     }
0527
0528     // check on RMAX
0529     Real_v rmax(0.);
0530     Real_v r2 = point.x() * point.x() + point.y() * point.y();
0531     // calculate cone radius at the z-height of position
0532     if (cone.fOriginalRmax1 == cone.fOriginalRmax2)
0533       rmax = Real_v(cone.fOriginalRmax1);
0534     else
0535       rmax = cone.fOuterSlope * point.z() + cone.fOuterOffset;
0536
0537     completelyoutside |= r2 > MakePlusTolerantSquare<ForInside>(rmax, cone.fOuterTolerance);
0538     if (ForInside) {
0539       completelyinside &= r2 < MakeMinusTolerantSquare<ForInside>(rmax, cone.fOuterTolerance);
0540     }
0541     if (vecCore::MaskFull(completelyoutside)) {
0542       return;
0543     }
0544
0545     // check on RMIN
0546     if (ConeTypes::checkRminTreatment<coneTypeT>(cone)) {
0547       Real_v rmin  = cone.fInnerSlope * point.z() + cone.fInnerOffset;
0548
0549       completelyoutside |= r2 <= MakeMinusTolerantSquare<ForInside>(rmin, cone.fInnerTolerance);
0550       if (ForInside) {
0551         completelyinside &= r2 > MakePlusTolerantSquare<ForInside>(rmin, cone.fInnerTolerance);
0552       }
0553       if (vecCore::MaskFull(completelyoutside)) {
0554         return;
0555       }
0556     }
0557
0558     if (ConeTypes::checkPhiTreatment<coneTypeT>(cone)) {
0559       Bool_t completelyoutsidephi;
0560       Bool_t completelyinsidephi;
0561       cone.fPhiWedge.GenericKernelForContainsAndInside<Real_v, ForInside>(point, completelyinsidephi,
0562                                                                           completelyoutsidephi);
0563       completelyoutside |= completelyoutsidephi;
0564       if (ForInside) completelyinside &= completelyinsidephi;
0565     }
0566   }
0567
0568   template <typename Inside_v>
0569   VECGEOM_FORCE_INLINE
0570   VECCORE_ATT_HOST_DEVICE
0571   static void Inside(UnplacedStruct_t const &cone, Vector3D<Real_v> const &point, Inside_v &inside)
0572   {
0573
0574     using Bool_v       = vecCore::Mask_v<Real_v>;
0575     using InsideBool_v = vecCore::Mask_v<Inside_v>;
0576     Bool_v completelyinside(false), completelyoutside(false);
0577     GenericKernelForContainsAndInside<true>(cone, point, completelyinside, completelyoutside);
0578     inside = EInside::kSurface;
0579     vecCore::MaskedAssign(inside, (InsideBool_v)completelyoutside, Inside_v(EInside::kOutside));
0580     vecCore::MaskedAssign(inside, (InsideBool_v)completelyinside, Inside_v(EInside::kInside));
0581   }
0582 };
0583
0584 template <class coneTypeT>
0585 class ConeHelpers<Precision, coneTypeT> {
0586
0587 public:
0588   ConeHelpers() {}
0589   ~ConeHelpers() {}
0590
0591   template <bool ForDistToIn, bool ForInnerSurface>
0592   VECCORE_ATT_HOST_DEVICE
0593   static bool DetectIntersectionAndCalculateDistanceToConicalSurface(UnplacedStruct_t const &cone,
0594                                                                      Vector3D<Precision> const &point,
0595                                                                      Vector3D<Precision> const &direction,
0596                                                                      Precision &distance)
0597   {
0598
0599     using namespace ConeUtilities;
0600     using namespace ConeTypes;
0601     distance              = kInfLength;
0602     bool onConicalSurface = IsOnConicalSurface<Precision, ForInnerSurface>(cone, point);
0603     Vector3D<Precision> normal = ConeUtilities::GetNormal<Precision, ForInnerSurface>(cone, point);
0604     bool tangentToSurface = vecCore::math::Abs(direction.Dot(normal)) == 0.;
0605     if (onConicalSurface && tangentToSurface) {
0606       return false;
0607     }
0608
0609     if (onConicalSurface) {
0610       if (ForDistToIn) {
0611     bool isMovingInside = IsMovingInsideConicalSurface<Precision, ForInnerSurface>(cone, point, direction);
0612
0613     if (!checkPhiTreatment<coneTypeT>(cone)) {
0614       if (isMovingInside) { // && onConicalSurface
0615         distance = 0.;
0616         return true;
0617       }
0618     } else {
0619       bool insector(false);
0620       ConeUtilities::PointInCyclicalSector<Precision, coneTypeT, false, true>(cone, point.x(), point.y(), insector);
0621       if (insector && isMovingInside) { // && onConicalSurface
0622         distance = 0.;
0623         return true;
0624       }
0625     }
0626       }
0627
0628       else {  // !ForDistToIn
0629     bool isMovingOutside = IsMovingOutsideConicalSurface<Precision, ForInnerSurface>(cone, point, direction);
0630
0631     if (!checkPhiTreatment<coneTypeT>(cone)) {
0632       if (isMovingOutside) { // && onConicalSurface
0633         distance = 0.;
0634         return true;
0635       }
0636     } else {
0637       bool insector(false);
0638       ConeUtilities::PointInCyclicalSector<Precision, coneTypeT, false, true>(cone, point.x(), point.y(), insector);
0639
0640       if (insector && isMovingOutside) {  // && onConicalSurface
0641         distance = 0.;
0642         return true;
0643       }
0644     }
0645       }
0646     }
0647
0648     bool ok(false);
0649     Precision pDotV2D = point.x() * direction.x() + point.y() * direction.y();
0650
0651     Precision a(kInfLength), b(kInfLength), c(kInfLength);
0652     if (ForInnerSurface) {
0653
0654       if (cone.fRmin1 == cone.fRmin2) {
0655         b = pDotV2D;
0656         a = direction.Perp2();
0657         c = point.Perp2() - cone.fRmin2 * cone.fRmin2;
0658       } else {
0659
0660         Precision newPz(0.);
0661         if (cone.fRmin2 > cone.fRmin1)
0662           newPz = (point.z() + cone.fDz + cone.fInnerConeApex) * cone.fTanInnerApexAngle;
0663         else
0664           newPz = (point.z() - cone.fDz - cone.fInnerConeApex) * cone.fTanInnerApexAngle;
0665
0666         Precision dirT = direction.z() * cone.fTanInnerApexAngle;
0667         a              = (direction.x() * direction.x()) + (direction.y() * direction.y()) - dirT * dirT;
0668
0669         b = pDotV2D - (newPz * dirT);
0670         c = point.Perp2() - (newPz * newPz);
0671       }
0672
0673       Precision b2 = b * b;
0674       Precision ac = a * c;
0675       if (b2 < ac) return false;
0676
0677       Precision d2 = b2 - ac;
0678
0679       Precision delta = Sqrt(d2);
0680       if (ForDistToIn) {
0681         if (b >= 0.) {
0682           distance = (c / NonZero(-b - delta));
0683         } else {
0684           distance = (-b + delta) / NonZero(a);
0685         }
0686       } else {
0687         if (b == 0. && delta == 0.) return false;
0688         if (b >= 0.) {
0689           distance = (-b - delta) / NonZero(a);
0690         } else {
0691           distance = (c / NonZero(-b + delta));
0692         }
0693       }
0694
0695       if (distance < 0.) return false;
0696       Precision newZ = point.z() + (direction.z() * distance);
0697       ok             = (Abs(newZ) < cone.fDz);
0698
0699     } else {
0700
0701       /*if (cone.fRmax1 == cone.fRmax2) {*/
0702       if (cone.fOriginalRmax1 == cone.fOriginalRmax2) {
0703
0704         a = direction.Perp2();
0705         b = pDotV2D;
0706         c = (point.Perp2() - cone.fOriginalRmax2 * cone.fOriginalRmax2);
0707       } else {
0708
0709         Precision newPz(0.);
0710         // if (cone.fRmax2 > cone.fRmax1)
0711         if (cone.fOriginalRmax2 > cone.fOriginalRmax1)
0712           newPz = (point.z() + cone.fDz + cone.fOuterConeApex) * cone.fTanOuterApexAngle;
0713         else
0714           newPz = (point.z() - cone.fDz - cone.fOuterConeApex) * cone.fTanOuterApexAngle;
0715         Precision dirT = direction.z() * cone.fTanOuterApexAngle;
0716         a              = direction.x() * direction.x() + direction.y() * direction.y() - dirT * dirT;
0717         b              = (pDotV2D - (newPz * dirT));
0718         c              = (point.Perp2() - (newPz * newPz));
0719       }
0720       Precision b2 = b * b;
0721       Precision ac = a * c;
0722       Precision d2 = b2 - ac;
0723       if (d2 < 0) return false;
0724       Precision delta = Sqrt(d2);
0725
0726       if (ForDistToIn) {
0727         if (b == 0. && delta == 0.) return false;
0728         if (b > 0.) {
0729           distance = (-b - delta) / NonZero(a); // BE ATTENTIVE, not covers the condition for b==0.
0730         } else {
0731           distance = (c / NonZero(-b + delta));
0732         }
0733         Precision newZ = point.z() + (direction.z() * distance);
0734         ok             = (Abs(newZ) < cone.fDz + kHalfTolerance);
0735       } else {
0736         if (b < 0.) {
0737           distance = (-b + delta) / NonZero(a);
0738         } else {
0739           distance = (c / NonZero(-b - delta));
0740         }
0741         ok = distance > 0.;
0742       }
0743
0744       if (distance < 0.) return false;
0745     }
0746     /*   if (distance < 0.) {
0747          distance = kInfLength;
0748        }
0749    */
0750     if (checkPhiTreatment<coneTypeT>(cone)) {
0751       Precision hitx(0), hity(0);
0752       bool insector(false);
0753       if (distance < kInfLength) {
0754         hitx = point.x() + distance * direction.x();
0755         hity = point.y() + distance * direction.y();
0756       }
0757
0758       ConeUtilities::PointInCyclicalSector<Precision, coneTypeT, false, true>(cone, hitx, hity, insector);
0759       ok &= ((insector) && (distance < kInfLength));
0760     }
0761     return ok;
0762   }
0763
0764   template <bool ForInside>
0765   VECGEOM_FORCE_INLINE
0766   VECCORE_ATT_HOST_DEVICE
0767   static void GenericKernelForContainsAndInside(UnplacedStruct_t const &cone, Vector3D<Precision> const &point,
0768                                                 bool &completelyinside, bool &completelyoutside)
0769   {
0770
0771     // very fast check on z-height
0772     Precision absz    = Abs(point[2]);
0773     completelyoutside = absz > MakePlusTolerant<ForInside>(cone.fDz, kConeTolerance);
0774     if (ForInside) {
0775       completelyinside = absz < MakeMinusTolerant<ForInside>(cone.fDz, kConeTolerance);
0776     }
0777     if (completelyoutside) return;
0778
0779     // check on RMAX
0780     Precision r2 = point.x() * point.x() + point.y() * point.y();
0781     // calculate cone radius at the z-height of position
0782     Precision rmax = 0.;
0783     if (cone.fOriginalRmax1 == cone.fOriginalRmax2)
0784       rmax = cone.fOriginalRmax1;
0785     else
0786       rmax = cone.fOuterSlope * point.z() + cone.fOuterOffset;
0787
0788     completelyoutside |= r2 > MakePlusTolerantSquare<ForInside>(rmax, cone.fOuterTolerance);
0789     if (ForInside) {
0790       completelyinside &= r2 < MakeMinusTolerantSquare<ForInside>(rmax, cone.fOuterTolerance);
0791     }
0792     if (completelyoutside) return;
0793
0794     // check on RMIN
0795     if (ConeTypes::checkRminTreatment<coneTypeT>(cone)) {
0796       Precision rmin  = cone.fInnerSlope * point.z() + cone.fInnerOffset;
0797
0798       completelyoutside |= r2 <= MakeMinusTolerantSquare<ForInside>(rmin, cone.fInnerTolerance);
0799       if (ForInside) {
0800         completelyinside &= r2 > MakePlusTolerantSquare<ForInside>(rmin, cone.fInnerTolerance);
0801       }
0802       if (completelyoutside) return;
0803     }
0804
0805     if (ConeTypes::checkPhiTreatment<coneTypeT>(cone)) {
0806       bool completelyoutsidephi(false);
0807       bool completelyinsidephi(false);
0808       cone.fPhiWedge.GenericKernelForContainsAndInside<Precision, ForInside>(point, completelyinsidephi,
0809                                                                              completelyoutsidephi);
0810       completelyoutside |= completelyoutsidephi;
0811       if (ForInside) completelyinside &= completelyinsidephi;
0812     }
0813   }
0814
0815   template <typename Inside_v>
0816   VECGEOM_FORCE_INLINE
0817   VECCORE_ATT_HOST_DEVICE
0818   static void Inside(UnplacedStruct_t const &cone, Vector3D<Precision> const &point, Inside_v &inside)
0819   {
0820     bool completelyinside(false), completelyoutside(false);
0821     GenericKernelForContainsAndInside<true>(cone, point, completelyinside, completelyoutside);
0822     inside = EInside::kSurface;
0823     if (completelyoutside) inside = EInside::kOutside;
0824     if (completelyinside) inside = EInside::kInside;
0825   }
0826 };
0827 } // namespace VECGEOM_IMPL_NAMESPACE
0828 } // namespace vecgeom
0829
0830 #endif