Back to home page

EIC code displayed by LXR
 
 

    

File indexing completed on 2026-09-14 09:26:27 

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