Back to home page

EIC code displayed by LXR
 
 

    

File indexing completed on 2025-09-18 09:36:45 

0001 /*
0002  * ThetaCone.h
0003  *
0004  *      Author: Raman Sehgal
0005  */
0006 
0007 #ifndef VECGEOM_VOLUMES_THETACONE_H_
0008 #define VECGEOM_VOLUMES_THETACONE_H_
0009 
0010 #include "VecGeom/base/Global.h"
0011 #include "VecGeom/volumes/kernel/GenericKernels.h"
0012 #include <VecCore/VecCore>
0013 #define kHalfPi 0.5 * kPi
0014 namespace vecgeom {
0015 inline namespace VECGEOM_IMPL_NAMESPACE {
0016 
0017 /**
0018  * A class representing a ThetaCone (basically a double cone) which is represented by an angle theta ( 0 < theta < Pi).
0019  *It
0020  *
0021  * The ThetaCone has an "startTheta" and "endTheta" angle. For an angle = 90 degree, the ThetaCone is essentially
0022  * XY plane with circular boundary. Usually the ThetaCone is used to cut out "theta" sections along z-direction.
0023  *
0024  *
0025  * Note: This class is meant as an auxiliary class so it is a bit outside the ordinary volume
0026  * hierarchy.
0027  *
0028  *      \ ++++ /
0029  *       \    /
0030  *        \  /
0031  *         \/
0032  *         /\
0033  *        /  \
0034  *       /    \
0035  *      / ++++ \
0036  *
0037  *DistanceToIn and DistanceToOut provides distances with the First and Second ThetaCone in "distThetaCone1" and
0038  *"distThetaCone2" reference variables.
0039  *Reference bool variable "intsect1" and "intsect2" is used to detect the real intersection cone, i.e. whether the point
0040  *really intersects with a ThetaCone or not.
0041  */
0042 class ThetaCone {
0043 
0044 private:
0045   Precision fSTheta; // starting angle
0046   Precision fDTheta; // delta angle
0047   Precision kAngTolerance;
0048   Precision halfAngTolerance;
0049   Precision fETheta; // ending angle
0050   Precision tanSTheta;
0051   Precision tanETheta;
0052   Precision tanBisector;
0053   Precision slope1, slope2;
0054   Precision tanSTheta2;
0055   Precision tanETheta2;
0056 
0057 public:
0058   VECCORE_ATT_HOST_DEVICE
0059   ThetaCone() {}
0060 
0061   VECCORE_ATT_HOST_DEVICE
0062   ThetaCone(Precision sTheta, Precision dTheta) : fSTheta(sTheta), fDTheta(dTheta), kAngTolerance(kTolerance)
0063   {
0064 
0065     // initialize angles
0066     fETheta               = fSTheta + fDTheta;
0067     halfAngTolerance      = (0.5 * kAngTolerance);
0068     Precision tempfSTheta = fSTheta;
0069     Precision tempfETheta = fETheta;
0070 
0071     if (fSTheta > kPi / 2) tempfSTheta = kPi - fSTheta;
0072     if (fETheta > kPi / 2) tempfETheta = kPi - fETheta;
0073 
0074     tanSTheta   = tan(tempfSTheta);
0075     tanSTheta2  = tanSTheta * tanSTheta;
0076     tanETheta   = tan(tempfETheta);
0077     tanETheta2  = tanETheta * tanETheta;
0078     tanBisector = tan(tempfSTheta + (fDTheta / 2));
0079     if (fSTheta > kPi / 2 && fETheta > kPi / 2) tanBisector = tan(tempfSTheta - (fDTheta / 2));
0080     slope1 = tan(kPi / 2 - fSTheta);
0081     slope2 = tan(kPi / 2 - fETheta);
0082   }
0083 
0084   VECCORE_ATT_HOST_DEVICE
0085   ~ThetaCone() {}
0086 
0087   VECCORE_ATT_HOST_DEVICE
0088   Precision GetSlope1() const { return slope1; }
0089 
0090   VECCORE_ATT_HOST_DEVICE
0091   Precision GetSlope2() const { return slope2; }
0092 
0093   VECCORE_ATT_HOST_DEVICE
0094   Precision GetTanSTheta2() const { return tanSTheta2; }
0095 
0096   VECCORE_ATT_HOST_DEVICE
0097   Precision GetTanETheta2() const { return tanETheta2; }
0098 
0099   /* Function to calculate normal at a point to the Cone formed at
0100    * by StartTheta.
0101    *
0102    * @inputs : Vector3D : Point at which normal needs to be calculated
0103    *
0104    * @output : Vector3D : calculated normal at the input point.
0105    */
0106   template <typename Real_v>
0107   VECCORE_ATT_HOST_DEVICE
0108   Vector3D<Real_v> GetNormal1(Vector3D<Real_v> const &point) const
0109   {
0110 
0111     Vector3D<Real_v> normal(2. * point.x(), 2. * point.y(), -2. * tanSTheta2 * point.z());
0112 
0113     if (fSTheta <= kPi / 2.)
0114       return -normal;
0115     else
0116       return normal;
0117   }
0118 
0119   /* Function to calculate normal at a point to the Cone formed at
0120    *  by EndTheta.
0121    *
0122    * @inputs : Vector3D : Point at which normal needs to be calculated
0123    *
0124    * @output : Vector3D : calculated normal at the input point.
0125    */
0126 
0127   template <typename Real_v>
0128   VECCORE_ATT_HOST_DEVICE
0129   Vector3D<Real_v> GetNormal2(Vector3D<Real_v> const &point) const
0130   {
0131 
0132     Vector3D<Real_v> normal(2 * point.x(), 2 * point.y(), -2 * tanETheta2 * point.z());
0133 
0134     if (fETheta <= kPi / 2.)
0135       return normal;
0136     else
0137       return -normal;
0138   }
0139 
0140   /* Function Name : GetNormal<Real_v, ForStartTheta>()
0141    *
0142    * The function is the templatized version GetNormal1() and GetNormal2() function with one more template
0143    * parameter and will return the normal depending upon the boolean template parameter "ForStartTheta"
0144    * which if passed as true, will return normal to the StartingTheta of ThetaCone,
0145    * if passed as false, will return normal to the EndingTheta of ThetaCone
0146    *
0147    * from user point of view the same work can be done by calling GetNormal1() and GetNormal2()
0148    * functions, but this implementation will be used by "IsPointOnSurfaceAndMovingOut()" function
0149    */
0150   template <typename Real_v, bool ForStartTheta>
0151   VECCORE_ATT_HOST_DEVICE
0152   Vector3D<Real_v> GetNormal(Vector3D<Real_v> const &point) const
0153   {
0154 
0155     if (ForStartTheta) {
0156 
0157       Vector3D<Real_v> normal(point.x(), point.y(), -tanSTheta2 * point.z());
0158 
0159       if (fSTheta <= kHalfPi)
0160         return -normal;
0161       else
0162         return normal;
0163 
0164     } else {
0165 
0166       Vector3D<Real_v> normal(point.x(), point.y(), -tanETheta2 * point.z());
0167 
0168       if (fETheta <= kHalfPi)
0169         return normal;
0170       else
0171         return -normal;
0172     }
0173   }
0174 
0175   /* Function Name :  IsOnSurfaceGeneric<Real_v, ForStartTheta>()
0176    *
0177    * This version of IsOnSurfaceGeneric is having one more template parameter of type boolean,
0178    * which if passed as true, will check whether the point is on StartingTheta Surface of ThetaCone,
0179    * and if passed as false, will check whether the point is on EndingTheta Surface of ThetaCone
0180    *
0181    * this implementation will be used by "IsPointOnSurfaceAndMovingOut()" function.
0182    */
0183   template <typename Real_v, bool ForStartTheta>
0184   VECCORE_ATT_HOST_DEVICE
0185   typename vecCore::Mask_v<Real_v> IsOnSurfaceGeneric(Vector3D<Real_v> const &point) const
0186   {
0187     Real_v rhs(0.);
0188     if (ForStartTheta) {
0189       rhs = Abs(tanSTheta * point.z());
0190     } else {
0191       rhs = Abs(tanETheta * point.z());
0192     }
0193     Real_v rho2 = point.Perp2();
0194     return rho2 >= MakeMinusTolerantSquare<true>(rhs) &&
0195            rho2 <= MakePlusTolerantSquare<true>(rhs);
0196   }
0197 
0198   /* Function Name : IsPointOnSurfaceAndMovingOut<Real_v, ForStartTheta, MovingOut>
0199    *
0200    * This function is written to check if the point is on surface and is moving inside or outside.
0201    * This will basically be used by "DistanceToInKernel()" and "DistanceToOutKernel()" of the shapes,
0202    * which uses ThetaCone.
0203    *
0204    * It contains two extra template boolean parameters "ForStartTheta" and "MovingOut",
0205    * So call like "IsPointOnSurfaceAndMovingOut<Real_v,true,true>" will check whether the points is on
0206    * the StartingTheta Surface of Theta and moving outside.
0207    *
0208    * So overall can be called in following four combinations
0209    * 1) "IsPointOnSurfaceAndMovingOut<Real_v,true,true>" : Point on StartingTheta surface of ThetaCone and moving OUT
0210    * 2) "IsPointOnSurfaceAndMovingOut<Real_v,true,false>" : Point on StartingTheta surface of ThetaCone and moving IN
0211    * 3) "IsPointOnSurfaceAndMovingOut<Real_v,false,true>" : Point on EndingTheta surface of ThetaCone and moving OUT
0212    * 4) "IsPointOnSurfaceAndMovingOut<Real_v,false,false>" : Point on EndingTheta surface of ThetaCone and moving IN
0213    *
0214    * Very useful for DistanceToIn and DistanceToOut.
0215    */
0216   template <typename Real_v, bool ForStartTheta, bool MovingOut>
0217   VECCORE_ATT_HOST_DEVICE
0218   typename vecCore::Mask_v<Real_v> IsPointOnSurfaceAndMovingOut(Vector3D<Real_v> const &point,
0219                                                                 Vector3D<Real_v> const &dir) const
0220   {
0221 
0222     if (MovingOut) {
0223       return IsOnSurfaceGeneric<Real_v, ForStartTheta>(point) &&
0224              (dir.Dot(GetNormal<Real_v, ForStartTheta>(point)) > Real_v(0.));
0225     } else {
0226       return IsOnSurfaceGeneric<Real_v, ForStartTheta>(point) &&
0227              (dir.Dot(GetNormal<Real_v, ForStartTheta>(point)) < Real_v(0.));
0228     }
0229   }
0230 
0231   template <typename Real_v>
0232   VECCORE_ATT_HOST_DEVICE
0233   typename vecCore::Mask_v<Real_v> Contains(Vector3D<Real_v> const &point) const
0234   {
0235 
0236     using Bool_v = vecCore::Mask_v<Real_v>;
0237     Bool_v unused(false);
0238     Bool_v outside(false);
0239     GenericKernelForContainsAndInside<Real_v, false>(point, unused, outside);
0240     return !outside;
0241   }
0242 
0243   template <typename Real_v>
0244   VECCORE_ATT_HOST_DEVICE
0245   typename vecCore::Mask_v<Real_v> ContainsWithBoundary(Vector3D<Real_v> const & /*point*/) const
0246   {
0247   }
0248 
0249   template <typename Real_v, typename Inside_t>
0250   VECCORE_ATT_HOST_DEVICE
0251   Inside_t Inside(Vector3D<Real_v> const &point) const
0252   {
0253     using Bool_v       = vecCore::Mask_v<Real_v>;
0254     using InsideBool_v = vecCore::Mask_v<Inside_t>;
0255     Bool_v completelyinside, completelyoutside;
0256     GenericKernelForContainsAndInside<Real_v, true>(point, completelyinside, completelyoutside);
0257     Inside_t inside(EInside::kSurface);
0258     vecCore::MaskedAssign(inside, (InsideBool_v)completelyoutside, Inside_t(EInside::kOutside));
0259     vecCore::MaskedAssign(inside, (InsideBool_v)completelyinside, Inside_t(EInside::kInside));
0260     return inside;
0261   }
0262 
0263   /**
0264    * estimate of the smallest distance to the ThetaCone boundary when
0265    * the point is located outside the ThetaCone
0266    */
0267   template <typename Real_v>
0268   VECCORE_ATT_HOST_DEVICE
0269   Real_v SafetyToIn(Vector3D<Real_v> const &point) const
0270   {
0271 
0272     using Bool_v = vecCore::Mask_v<Real_v>;
0273 
0274     Real_v safeTheta(0.);
0275     Real_v pointRad = Sqrt(point.x() * point.x() + point.y() * point.y());
0276     Real_v sfTh1    = DistanceToLine<Real_v>(slope1, pointRad, point.z());
0277     Real_v sfTh2    = DistanceToLine<Real_v>(slope2, pointRad, point.z());
0278 
0279     safeTheta   = Min(sfTh1, sfTh2);
0280     Bool_v done = Contains<Real_v>(point);
0281     vecCore__MaskedAssignFunc(safeTheta, done, Real_v(0.));
0282     if (vecCore::MaskFull(done)) return safeTheta;
0283 
0284     // Case 1 : Both cones are in Positive Z direction
0285     if (fSTheta < kPi / 2 + halfAngTolerance) {
0286       if (fETheta < kPi / 2 + halfAngTolerance) {
0287         if (fSTheta < fETheta) {
0288           vecCore::MaskedAssign(safeTheta, (!done && point.z() < Real_v(0.)), sfTh2);
0289         }
0290       }
0291 
0292       // Case 2 : First Cone is in Positive Z direction and Second is in Negative Z direction
0293       if (fETheta > kPi / 2 + halfAngTolerance) {
0294         if (fSTheta < fETheta) {
0295           vecCore::MaskedAssign(safeTheta, (!done && point.z() > Real_v(0.)), sfTh1);
0296           vecCore::MaskedAssign(safeTheta, (!done && point.z() < Real_v(0.)), sfTh2);
0297         }
0298       }
0299     }
0300 
0301     // Case 3 : Both cones are in Negative Z direction
0302     if (fETheta > kPi / 2 + halfAngTolerance) {
0303       if (fSTheta > kPi / 2 + halfAngTolerance) {
0304         if (fSTheta < fETheta) {
0305           vecCore::MaskedAssign(safeTheta, (!done && point.z() > Real_v(0.)), sfTh1);
0306         }
0307       }
0308     }
0309 
0310     return safeTheta;
0311   }
0312 
0313   /**
0314    * estimate of the smallest distance to the ThetaCone boundary when
0315    * the point is located inside the ThetaCone ( within the defining phi angle )
0316    */
0317   template <typename Real_v>
0318   VECCORE_ATT_HOST_DEVICE
0319   Real_v SafetyToOut(Vector3D<Real_v> const &point) const
0320   {
0321 
0322     Real_v pointRad    = Sqrt(point.x() * point.x() + point.y() * point.y());
0323     Real_v bisectorRad = Abs(point.z() * tanBisector);
0324 
0325     vecCore::Mask<Real_v> condition(false);
0326     Real_v sfTh1 = DistanceToLine<Real_v>(slope1, pointRad, point.z());
0327     Real_v sfTh2 = DistanceToLine<Real_v>(slope2, pointRad, point.z());
0328 
0329     // Case 1 : Both cones are in Positive Z direction
0330     if (fSTheta < kPi / 2 + halfAngTolerance) {
0331       if (fETheta < kPi / 2 + halfAngTolerance) {
0332         if (fSTheta < fETheta) {
0333           condition = (pointRad < bisectorRad) && (fSTheta != Real_v(0.));
0334         }
0335       }
0336 
0337       // Case 2 : First Cone is in Positive Z direction and Second is in Negative Z direction
0338       if (fETheta > kPi / 2 + halfAngTolerance) {
0339         if (fSTheta < fETheta) {
0340           condition = sfTh1 < sfTh2;
0341         }
0342       }
0343     }
0344 
0345     // Case 3 : Both cones are in Negative Z direction
0346     if (fETheta > kPi / 2 + halfAngTolerance) {
0347       if (fSTheta > kPi / 2 + halfAngTolerance) {
0348         if (fSTheta < fETheta) {
0349           condition = !((pointRad < bisectorRad) && (fETheta != Real_v(kPi)));
0350         }
0351       }
0352     }
0353 
0354     return vecCore::Blend(condition, sfTh1, sfTh2);
0355   }
0356 
0357   template <typename Real_v>
0358   VECCORE_ATT_HOST_DEVICE
0359   Real_v DistanceToLine(Precision const &slope, Real_v const &x, Real_v const &y) const
0360   {
0361 
0362     Real_v dist = (y - slope * x) / Sqrt(Real_v(1.) + slope * slope);
0363     return Abs(dist);
0364   }
0365 
0366   /**
0367    * estimate of the distance to the ThetaCone boundary with given direction
0368    */
0369   template <typename Real_v>
0370   VECCORE_ATT_HOST_DEVICE
0371   void DistanceToIn(Vector3D<Real_v> const &point, Vector3D<Real_v> const &dir, Real_v &distThetaCone1,
0372                     Real_v &distThetaCone2, typename vecCore::Mask_v<Real_v> &intsect1,
0373                     typename vecCore::Mask_v<Real_v> &intsect2) const
0374   {
0375 
0376     {
0377       using Bool_v = vecCore::Mask_v<Real_v>;
0378 
0379       Bool_v done(false);
0380       Bool_v fal(false);
0381 
0382       Real_v firstRoot(kInfLength), secondRoot(kInfLength);
0383 
0384       Real_v pDotV2d = point.x() * dir.x() + point.y() * dir.y();
0385       Real_v rho2    = point.x() * point.x() + point.y() * point.y();
0386       Real_v dirRho2 = dir.Perp2();
0387 
0388       Real_v b = pDotV2d - point.z() * dir.z() * tanSTheta2;
0389       // Real_v a = dir.x() * dir.x() + dir.y() * dir.y() - dir.z() * dir.z() * tanSTheta2;
0390       Real_v a    = dirRho2 - dir.z() * dir.z() * tanSTheta2;
0391       Real_v c    = rho2 - point.z() * point.z() * tanSTheta2;
0392       Real_v d2   = b * b - a * c;
0393       Real_v aInv = Real_v(1.) / NonZero(a);
0394 
0395       vecCore__MaskedAssignFunc(firstRoot, (d2 > Real_v(0.)), (-b + Sqrt(Abs(d2))) * aInv);
0396       done |= (Abs(firstRoot) < Real_v(3.) * kTolerance);
0397       vecCore__MaskedAssignFunc(firstRoot, ((Abs(firstRoot) < Real_v(3.) * kTolerance)), Real_v(0.));
0398       vecCore__MaskedAssignFunc(firstRoot, (!done && (firstRoot < Real_v(0.))), InfinityLength<Real_v>());
0399 
0400       Real_v b2 = pDotV2d - point.z() * dir.z() * tanETheta2;
0401       // Real_v a2 = dir.x() * dir.x() + dir.y() * dir.y() - dir.z() * dir.z() * tanETheta2;
0402       Real_v a2    = dirRho2 - dir.z() * dir.z() * tanETheta2;
0403       Real_v c2    = rho2 - point.z() * point.z() * tanETheta2;
0404       Real_v d22   = b2 * b2 - a2 * c2;
0405       Real_v a2Inv = Real_v(1.) / NonZero(a2);
0406 
0407       vecCore__MaskedAssignFunc(secondRoot, (d22 > Real_v(0.)), (-b2 - Sqrt(Abs(d22))) * a2Inv);
0408       vecCore__MaskedAssignFunc(secondRoot, (!done && (Abs(secondRoot) < Real_v(3.) * kTolerance)), Real_v(0.));
0409       done |= (Abs(secondRoot) < Real_v(3.) * kTolerance);
0410       vecCore__MaskedAssignFunc(secondRoot, !done && (secondRoot < Real_v(0.)), InfinityLength<Real_v>());
0411 
0412       if (fSTheta < kHalfPi + halfAngTolerance) {
0413         if (fETheta < kHalfPi + halfAngTolerance) {
0414           if (fSTheta < fETheta) {
0415             distThetaCone1          = firstRoot;
0416             distThetaCone2          = secondRoot;
0417             Real_v zOfIntSecPtCone1 = (point.z() + distThetaCone1 * dir.z());
0418             Real_v zOfIntSecPtCone2 = (point.z() + distThetaCone2 * dir.z());
0419 
0420             intsect1 = ((d2 > Real_v(0.)) && (zOfIntSecPtCone1 > Real_v(0.)));
0421             intsect2 = ((d22 > Real_v(0.)) && (zOfIntSecPtCone2 > Real_v(0.)));
0422           }
0423         }
0424 
0425         if (fETheta >= kHalfPi - halfAngTolerance && fETheta <= kHalfPi + halfAngTolerance) {
0426           vecCore__MaskedAssignFunc(distThetaCone2, (dir.z() > Real_v(0.)), -point.z() / dir.z());
0427           Real_v zOfIntSecPtCone2 = (point.z() + distThetaCone2 * dir.z());
0428           intsect2                = ((distThetaCone2 != kInfLength) && (Abs(zOfIntSecPtCone2) < halfAngTolerance));
0429         }
0430 
0431         if (fETheta > kHalfPi + halfAngTolerance) {
0432           if (fSTheta < fETheta) {
0433             distThetaCone1 = firstRoot;
0434             vecCore__MaskedAssignFunc(secondRoot, (d22 > Real_v(0.)), (-b2 + Sqrt(Abs(d22))) * a2Inv);
0435 
0436             done = fal;
0437             done |= (Abs(secondRoot) < Real_v(3.) * kTolerance);
0438             vecCore__MaskedAssignFunc(secondRoot, ((Abs(secondRoot) < Real_v(3.) * kTolerance)), Real_v(0.));
0439             vecCore__MaskedAssignFunc(secondRoot, !done && (secondRoot < Real_v(0.)), InfinityLength<Real_v>());
0440             distThetaCone2 = secondRoot;
0441 
0442             Real_v zOfIntSecPtCone1 = (point.z() + distThetaCone1 * dir.z());
0443             Real_v zOfIntSecPtCone2 = (point.z() + distThetaCone2 * dir.z());
0444 
0445             intsect1 = ((d2 > 0) && (distThetaCone1 != kInfLength) && (zOfIntSecPtCone1 > Real_v(0.)));
0446             intsect2 = ((d22 > 0) && (distThetaCone2 != kInfLength) && (zOfIntSecPtCone2 < Real_v(0.)));
0447           }
0448         }
0449       }
0450 
0451       if (fSTheta >= kHalfPi - halfAngTolerance) {
0452         if (fETheta > kHalfPi + halfAngTolerance) {
0453           if (fSTheta < fETheta) {
0454             vecCore__MaskedAssignFunc(firstRoot, (d2 > Real_v(0.)), (-b - Sqrt(Abs(d2))) * aInv);
0455             done = fal;
0456             done |= (Abs(firstRoot) < Real_v(3.) * kTolerance);
0457             vecCore__MaskedAssignFunc(firstRoot, ((Abs(firstRoot) < Real_v(3.) * kTolerance)), Real_v(0.));
0458             vecCore__MaskedAssignFunc(firstRoot, !done && (firstRoot < Real_v(0.)), InfinityLength<Real_v>());
0459             distThetaCone1 = firstRoot;
0460 
0461             vecCore__MaskedAssignFunc(secondRoot, (d22 > Real_v(0.)), (-b2 + Sqrt(Abs(d22))) * a2Inv);
0462             done = fal;
0463             done |= (Abs(secondRoot) < Real_v(3.) * kTolerance);
0464             vecCore__MaskedAssignFunc(secondRoot, ((Abs(secondRoot) < Real_v(3.) * kTolerance)), Real_v(0.));
0465             vecCore__MaskedAssignFunc(secondRoot, !done && (secondRoot < Real_v(0.)), InfinityLength<Real_v>());
0466             distThetaCone2 = secondRoot;
0467 
0468             Real_v zOfIntSecPtCone1 = (point.z() + distThetaCone1 * dir.z());
0469             Real_v zOfIntSecPtCone2 = (point.z() + distThetaCone2 * dir.z());
0470 
0471             intsect1 = ((d2 > 0) && (distThetaCone1 != kInfLength) && (zOfIntSecPtCone1 < Real_v(0.)));
0472             intsect2 = ((d22 > 0) && (distThetaCone2 != kInfLength) && (zOfIntSecPtCone2 < Real_v(0.)));
0473           }
0474         }
0475       }
0476 
0477       if (fSTheta >= kHalfPi - halfAngTolerance && fSTheta <= kHalfPi + halfAngTolerance) {
0478         vecCore__MaskedAssignFunc(distThetaCone1, (dir.z() < Real_v(0.)), -point.z() / dir.z());
0479         Real_v zOfIntSecPtCone1 = (point.z() + distThetaCone1 * dir.z());
0480         intsect1                = ((distThetaCone1 != kInfLength) && (Abs(zOfIntSecPtCone1) < halfAngTolerance));
0481       }
0482     }
0483   }
0484 
0485   template <typename Real_v>
0486   VECCORE_ATT_HOST_DEVICE
0487   void DistanceToOut(Vector3D<Real_v> const &point, Vector3D<Real_v> const &dir, Real_v &distThetaCone1,
0488                      Real_v &distThetaCone2, typename vecCore::Mask_v<Real_v> &intsect1,
0489                      typename vecCore::Mask_v<Real_v> &intsect2) const
0490   {
0491 
0492     using Bool_v = vecCore::Mask_v<Real_v>;
0493 
0494     Real_v inf(kInfLength);
0495 
0496     // Real_v firstRoot(kInfLength), secondRoot(kInfLength);
0497     Real_v firstRoot(kInfLength), secondRoot(kInfLength);
0498 
0499     Real_v pDotV2d = point.x() * dir.x() + point.y() * dir.y();
0500     Real_v rho2    = point.x() * point.x() + point.y() * point.y();
0501 
0502     Real_v b  = pDotV2d - point.z() * dir.z() * tanSTheta2;
0503     Real_v a  = dir.x() * dir.x() + dir.y() * dir.y() - dir.z() * dir.z() * tanSTheta2;
0504     Real_v c  = rho2 - point.z() * point.z() * tanSTheta2;
0505     Real_v d2 = b * b - a * c;
0506     vecCore__MaskedAssignFunc(d2, d2 < Real_v(0.) && Abs(d2) < kHalfTolerance, Real_v(0.));
0507     vecCore__MaskedAssignFunc(firstRoot, (d2 >= Real_v(0.)) && b >= Real_v(0.) && a != Real_v(0.),
0508                               ((-b - Sqrt(Abs(d2))) / NonZero(a)));
0509     vecCore__MaskedAssignFunc(firstRoot, (d2 >= Real_v(0.)) && b < Real_v(0.), ((c) / NonZero(-b + Sqrt(Abs(d2)))));
0510     vecCore__MaskedAssignFunc(firstRoot, firstRoot < Real_v(0.), InfinityLength<Real_v>());
0511 
0512     Real_v b2 = point.x() * dir.x() + point.y() * dir.y() - point.z() * dir.z() * tanETheta2;
0513     Real_v a2 = dir.x() * dir.x() + dir.y() * dir.y() - dir.z() * dir.z() * tanETheta2;
0514 
0515     Real_v c2  = point.x() * point.x() + point.y() * point.y() - point.z() * point.z() * tanETheta2;
0516     Real_v d22 = (b2 * b2) - (a2 * c2);
0517     vecCore__MaskedAssignFunc(d22, d22 < Real_v(0.) && Abs(d22) < kHalfTolerance, Real_v(0.));
0518 
0519     vecCore__MaskedAssignFunc(secondRoot, (d22 >= Real_v(0.)) && b2 >= Real_v(0.),
0520                               ((c2) / NonZero(-b2 - Sqrt(Abs(d22)))));
0521     vecCore__MaskedAssignFunc(secondRoot, (d22 >= Real_v(0.)) && b2 < Real_v(0.) && a2 != Real_v(0.),
0522                               ((-b2 + Sqrt(Abs(d22))) / NonZero(a2)));
0523 
0524     vecCore__MaskedAssignFunc(secondRoot, secondRoot < Real_v(0.) && Abs(secondRoot) > kTolerance,
0525                               InfinityLength<Real_v>());
0526     vecCore__MaskedAssignFunc(secondRoot, Abs(secondRoot) < kTolerance, Real_v(0.));
0527 
0528     if (fSTheta < kPi / 2 + halfAngTolerance) {
0529       if (fETheta < kPi / 2 + halfAngTolerance) {
0530         if (fSTheta < fETheta) {
0531           distThetaCone1          = firstRoot;
0532           distThetaCone2          = secondRoot;
0533           Real_v zOfIntSecPtCone1 = (point.z() + distThetaCone1 * dir.z());
0534           Real_v zOfIntSecPtCone2 = (point.z() + distThetaCone2 * dir.z());
0535 
0536           intsect1 = ((d2 > 0) && (distThetaCone1 != kInfLength) && ((zOfIntSecPtCone1) > -kHalfTolerance));
0537           intsect2 = ((d22 > 0) && (distThetaCone2 != kInfLength) && ((zOfIntSecPtCone2) > -kHalfTolerance));
0538 
0539           Real_v dirRho2 = dir.x() * dir.x() + dir.y() * dir.y();
0540           Real_v zs(kInfLength);
0541           if (fSTheta) zs = dirRho2 / tanSTheta;
0542           Real_v ze(kInfLength);
0543           if (fETheta) ze = dirRho2 / tanETheta;
0544           Bool_v cond = (point.x() == Real_v(0.) && point.y() == Real_v(0.) && point.z() == Real_v(0.) &&
0545                          dir.z() < zs && dir.z() < ze);
0546           vecCore__MaskedAssignFunc(distThetaCone1, cond, Real_v(0.));
0547           vecCore__MaskedAssignFunc(distThetaCone2, cond, Real_v(0.));
0548           intsect1 |= cond;
0549           intsect2 |= cond;
0550         }
0551       }
0552 
0553       // Second cone is the XY plane, compute distance to plane
0554       if (fETheta >= kPi / 2 - halfAngTolerance && fETheta <= kPi / 2 + halfAngTolerance) {
0555         distThetaCone1 = firstRoot;
0556         distThetaCone2 = inf;
0557         vecCore__MaskedAssignFunc(distThetaCone2, (dir.z() < Real_v(0.)), Real_v(-1.) * point.z() / dir.z());
0558         Real_v zOfIntSecPtCone2 = (point.z() + distThetaCone2 * dir.z());
0559         Real_v zOfIntSecPtCone1 = (point.z() + distThetaCone1 * dir.z());
0560         intsect2 =
0561             ((distThetaCone2 != kInfLength) && (Abs(zOfIntSecPtCone2) < kHalfTolerance) && !(dir.z() == Real_v(0.)));
0562         intsect1 = ((d2 >= 0) && (distThetaCone1 != kInfLength) && (Abs(zOfIntSecPtCone1) < kHalfTolerance) &&
0563                     !(dir.z() == Real_v(0.)));
0564       }
0565 
0566       if (fETheta > kPi / 2 + halfAngTolerance) {
0567         if (fSTheta < fETheta) {
0568           distThetaCone1 = firstRoot;
0569           vecCore__MaskedAssignFunc(secondRoot, (d22 >= Real_v(0.)) && b2 > Real_v(0.) && a2 != Real_v(0.),
0570                                     ((-b2 - Sqrt(Abs(d22))) / NonZero(a2)));
0571           vecCore__MaskedAssignFunc(secondRoot, (d22 >= Real_v(0.)) && b2 <= Real_v(0.),
0572                                     ((c2) / NonZero(-b2 + Sqrt(Abs(d22)))));
0573           vecCore__MaskedAssignFunc(secondRoot, secondRoot < Real_v(0.), InfinityLength<Real_v>());
0574           distThetaCone2          = secondRoot;
0575           Real_v zOfIntSecPtCone1 = (point.z() + distThetaCone1 * dir.z());
0576           Real_v zOfIntSecPtCone2 = (point.z() + distThetaCone2 * dir.z());
0577           intsect1 = ((d2 >= 0) && (distThetaCone1 != kInfLength) && ((zOfIntSecPtCone1) > -kHalfTolerance));
0578           intsect2 = ((d22 >= 0) && (distThetaCone2 != kInfLength) && ((zOfIntSecPtCone2) < kHalfTolerance));
0579         }
0580       }
0581     }
0582 
0583     if (fETheta > kPi / 2 + halfAngTolerance) {
0584       if (fSTheta < fETheta) {
0585         secondRoot = kInfLength;
0586         vecCore__MaskedAssignFunc(secondRoot, (d22 >= Real_v(0.)) && b2 > Real_v(0.) && a2 != Real_v(0.),
0587                                   ((-b2 - Sqrt(Abs(d22))) / NonZero(a2)));
0588         vecCore__MaskedAssignFunc(secondRoot, (d22 >= Real_v(0.)) && b2 <= Real_v(0.),
0589                                   ((c2) / NonZero(-b2 + Sqrt(Abs(d22)))));
0590         distThetaCone2 = secondRoot;
0591 
0592         if (fSTheta > kPi / 2 + halfAngTolerance) {
0593           vecCore__MaskedAssignFunc(firstRoot, (d2 >= Real_v(0.)) && b > Real_v(0.),
0594                                     ((c) / NonZero(-b - Sqrt(Abs(d2)))));
0595           vecCore__MaskedAssignFunc(firstRoot, (d2 >= Real_v(0.)) && b <= Real_v(0.) && a != Real_v(0.),
0596                                     ((-b + Sqrt(Abs(d2))) / NonZero(a)));
0597           vecCore__MaskedAssignFunc(firstRoot, firstRoot < Real_v(0.), InfinityLength<Real_v>());
0598           distThetaCone1          = firstRoot;
0599           Real_v zOfIntSecPtCone1 = (point.z() + distThetaCone1 * dir.z());
0600           intsect1 = ((d2 > 0) && (distThetaCone1 != kInfLength) && ((zOfIntSecPtCone1) < kHalfTolerance));
0601           Real_v zOfIntSecPtCone2 = (point.z() + distThetaCone2 * dir.z());
0602           intsect2 = ((d22 > 0) && (distThetaCone2 != kInfLength) && ((zOfIntSecPtCone2) < kHalfTolerance));
0603 
0604           Real_v dirRho2 = dir.x() * dir.x() + dir.y() * dir.y();
0605           Real_v zs(-kInfLength);
0606           if (tanSTheta) zs = -dirRho2 / tanSTheta;
0607           Real_v ze(-kInfLength);
0608           if (tanETheta) ze = -dirRho2 / tanETheta;
0609           Bool_v cond = (point.x() == Real_v(0.) && point.y() == Real_v(0.) && point.z() == Real_v(0.) &&
0610                          dir.z() > zs && dir.z() > ze);
0611           vecCore__MaskedAssignFunc(distThetaCone1, cond, Real_v(0.));
0612           vecCore__MaskedAssignFunc(distThetaCone2, cond, Real_v(0.));
0613           // intsect1 |= (cond && tr);
0614           // intsect2 |= (cond && tr);
0615           intsect1 |= cond;
0616           intsect2 |= cond;
0617         }
0618       }
0619     }
0620 
0621     // First cone is the XY plane, compute distance to plane
0622     if (fSTheta >= kPi / 2 - halfAngTolerance && fSTheta <= kPi / 2 + halfAngTolerance) {
0623       distThetaCone2 = secondRoot;
0624       distThetaCone1 = kInfLength;
0625       vecCore__MaskedAssignFunc(distThetaCone1, (dir.z() > Real_v(0.)), Real_v(-1.) * point.z() / NonZero(dir.z()));
0626       Real_v zOfIntSecPtCone1 = (point.z() + distThetaCone1 * dir.z());
0627 
0628       Real_v zOfIntSecPtCone2 = (point.z() + distThetaCone2 * dir.z());
0629 
0630       intsect1 =
0631           ((distThetaCone1 != kInfLength) && (Abs(zOfIntSecPtCone1) < kHalfTolerance) && (dir.z() != Real_v(0.)));
0632       intsect2 = ((d22 >= 0) && (distThetaCone2 != kInfLength) && (Abs(zOfIntSecPtCone2) < kHalfTolerance) &&
0633                   (dir.z() != Real_v(0.)));
0634     }
0635   }
0636 
0637   // This could be useful in case somebody just want to check whether point is completely inside ThetaRange
0638   template <typename Real_v>
0639   VECCORE_ATT_HOST_DEVICE
0640   typename vecCore::Mask_v<Real_v> IsCompletelyInside(Vector3D<Real_v> const &localPoint) const
0641   {
0642 
0643     using Bool_v       = vecCore::Mask_v<Real_v>;
0644     Real_v rho         = Sqrt(localPoint.Mag2() - (localPoint.z() * localPoint.z()));
0645     Real_v cone1Radius = Abs(localPoint.z() * tanSTheta);
0646     Real_v cone2Radius = Abs(localPoint.z() * tanETheta);
0647     Bool_v isPointOnZAxis =
0648         localPoint.z() != Real_v(0.) && localPoint.x() == Real_v(0.) && localPoint.y() == Real_v(0.);
0649 
0650     Bool_v isPointOnXYPlane =
0651         localPoint.z() == Real_v(0.) && (localPoint.x() != Real_v(0.) || localPoint.y() != Real_v(0.));
0652 
0653     Real_v startTheta(fSTheta), endTheta(fETheta);
0654 
0655     Bool_v completelyinside = (isPointOnZAxis && ((startTheta == Real_v(0.) && endTheta == kPi) ||
0656                                                   (localPoint.z() > Real_v(0.) && startTheta == Real_v(0.)) ||
0657                                                   (localPoint.z() < Real_v(0.) && endTheta == kPi)));
0658 
0659     completelyinside |=
0660         (!completelyinside && (isPointOnXYPlane && (startTheta < Real_v(kHalfPi) && endTheta > Real_v(kHalfPi) &&
0661                                                     (Real_v(kHalfPi) - startTheta) > kAngTolerance &&
0662                                                     (endTheta - Real_v(kHalfPi)) > kTolerance)));
0663 
0664     if (fSTheta < kHalfPi + halfAngTolerance) {
0665       if (fETheta < kHalfPi + halfAngTolerance) {
0666         if (fSTheta < fETheta) {
0667           Real_v tolAngMin = cone1Radius + Real_v(2 * kAngTolerance * 10.);
0668           Real_v tolAngMax = cone2Radius - Real_v(2 * kAngTolerance * 10.);
0669 
0670           completelyinside |=
0671               (!completelyinside &&
0672                (((rho <= tolAngMax) && (rho >= tolAngMin) && (localPoint.z() > Real_v(0.)) &&
0673                  Bool_v(fSTheta != Real_v(0.))) ||
0674                 ((rho <= tolAngMax) && Bool_v(fSTheta == Real_v(0.)) && (localPoint.z() > Real_v(0.)))));
0675         }
0676       }
0677 
0678       if (fETheta > kHalfPi + halfAngTolerance) {
0679         if (fSTheta < fETheta) {
0680           Real_v tolAngMin = cone1Radius + Real_v(2 * kAngTolerance * 10.);
0681           Real_v tolAngMax = cone2Radius + Real_v(2 * kAngTolerance * 10.);
0682 
0683           completelyinside |= (!completelyinside && (((rho >= tolAngMin) && (localPoint.z() > Real_v(0.))) ||
0684                                                      ((rho >= tolAngMax) && (localPoint.z() < Real_v(0.)))));
0685         }
0686       }
0687 
0688       if (fETheta >= kHalfPi - halfAngTolerance && fETheta <= kHalfPi + halfAngTolerance) {
0689 
0690         completelyinside &= !(Abs(localPoint.z()) < halfAngTolerance);
0691       }
0692     }
0693 
0694     if (fETheta > kHalfPi + halfAngTolerance) {
0695       if (fSTheta >= kHalfPi - halfAngTolerance && fSTheta <= kHalfPi + halfAngTolerance) {
0696 
0697         completelyinside &= !(Abs(localPoint.z()) < halfAngTolerance);
0698       }
0699 
0700       if (fSTheta > kHalfPi + halfAngTolerance) {
0701         if (fSTheta < fETheta) {
0702           Real_v tolAngMin = cone1Radius - Real_v(2 * kAngTolerance * 10.);
0703           Real_v tolAngMax = cone2Radius + Real_v(2 * kAngTolerance * 10.);
0704 
0705           completelyinside |=
0706               (!completelyinside &&
0707                (((rho <= tolAngMin) && (rho >= tolAngMax) && (localPoint.z() < Real_v(0.)) && Bool_v(fETheta != kPi)) ||
0708                 ((rho <= tolAngMin) && (localPoint.z() < Real_v(0.)) && Bool_v(fETheta == kPi))));
0709         }
0710       }
0711     }
0712     return completelyinside;
0713   }
0714 
0715   // This could be useful in case somebody just want to check whether point is completely outside ThetaRange
0716   template <typename Real_v>
0717   VECCORE_ATT_HOST_DEVICE
0718   typename vecCore::Mask_v<Real_v> IsCompletelyOutside(Vector3D<Real_v> const &localPoint) const
0719   {
0720 
0721     using Bool_v       = vecCore::Mask_v<Real_v>;
0722     Real_v diff        = localPoint.Perp2();
0723     Real_v rho         = Sqrt(diff);
0724     Real_v cone1Radius = Abs(localPoint.z() * tanSTheta);
0725     Real_v cone2Radius = Abs(localPoint.z() * tanETheta);
0726     Bool_v isPointOnZAxis =
0727         localPoint.z() != Real_v(0.) && localPoint.x() == Real_v(0.) && localPoint.y() == Real_v(0.);
0728 
0729     Bool_v isPointOnXYPlane =
0730         localPoint.z() == Real_v(0.) && (localPoint.x() != Real_v(0.) || localPoint.y() != Real_v(0.));
0731 
0732     // Real_v startTheta(fSTheta), endTheta(fETheta);
0733 
0734     Bool_v completelyoutside = (isPointOnZAxis && ((Bool_v(fSTheta != Real_v(0.)) && Bool_v(fETheta != kPi)) ||
0735                                                    (localPoint.z() > Real_v(0.) && Bool_v(fSTheta != Real_v(0.))) ||
0736                                                    (localPoint.z() < Real_v(0.) && Bool_v(fETheta != kPi))));
0737 
0738     completelyoutside |=
0739         (!completelyoutside &&
0740          (isPointOnXYPlane && Bool_v(((fSTheta < kHalfPi) && (fETheta < kHalfPi) &&
0741                                       ((kHalfPi - fSTheta) > kAngTolerance) && ((kHalfPi - fETheta) > kTolerance)) ||
0742                                      ((fSTheta > kHalfPi && fETheta > kHalfPi) &&
0743                                       ((fSTheta - kHalfPi) > kAngTolerance) && ((fETheta - kHalfPi) > kTolerance)))));
0744 
0745     if (fSTheta < kHalfPi + halfAngTolerance) {
0746       if (fETheta < kHalfPi + halfAngTolerance) {
0747         if (fSTheta < fETheta) {
0748 
0749           Real_v tolAngMin2 = cone1Radius - Real_v(2 * kAngTolerance * 10.);
0750           Real_v tolAngMax2 = cone2Radius + Real_v(2 * kAngTolerance * 10.);
0751 
0752           completelyoutside |=
0753               (!completelyoutside && ((rho < tolAngMin2) || (rho > tolAngMax2) || (localPoint.z() < Real_v(0.))));
0754         }
0755       }
0756 
0757       if (fETheta > kHalfPi + halfAngTolerance) {
0758         if (fSTheta < fETheta) {
0759           Real_v tolAngMin2 = cone1Radius - Real_v(2 * kAngTolerance * 10.);
0760           Real_v tolAngMax2 = cone2Radius - Real_v(2 * kAngTolerance * 10.);
0761           completelyoutside |= (!completelyoutside && (((rho < tolAngMin2) && (localPoint.z() > Real_v(0.))) ||
0762                                                        ((rho < tolAngMax2) && (localPoint.z() < Real_v(0.)))));
0763         }
0764       }
0765 
0766       if (fETheta >= kHalfPi - halfAngTolerance && fETheta <= kHalfPi + halfAngTolerance) {
0767 
0768         // completelyinside &= !(Abs(localPoint.z()) < halfAngTolerance);
0769         completelyoutside &= !(Abs(localPoint.z()) < halfAngTolerance);
0770       }
0771     }
0772 
0773     if (fETheta > kHalfPi + halfAngTolerance) {
0774       if (fSTheta >= kHalfPi - halfAngTolerance && fSTheta <= kHalfPi + halfAngTolerance) {
0775 
0776         // completelyinside &= !(Abs(localPoint.z()) < halfAngTolerance);
0777         completelyoutside &= !(Abs(localPoint.z()) < halfAngTolerance);
0778       }
0779 
0780       if (fSTheta > kHalfPi + halfAngTolerance) {
0781         if (fSTheta < fETheta) {
0782 
0783           Real_v tolAngMin2 = cone1Radius + Real_v(2 * kAngTolerance * 10.);
0784           Real_v tolAngMax2 = cone2Radius - Real_v(2 * kAngTolerance * 10.);
0785           completelyoutside |=
0786               (!completelyoutside && ((rho < tolAngMax2) || (rho > tolAngMin2) || (localPoint.z() > Real_v(0.))));
0787         }
0788       }
0789     }
0790 
0791     return completelyoutside;
0792   }
0793 
0794   template <typename Real_v, bool ForInside>
0795   VECCORE_ATT_HOST_DEVICE
0796   void GenericKernelForContainsAndInside(Vector3D<Real_v> const &localPoint,
0797                                          typename vecCore::Mask_v<Real_v> &completelyinside,
0798                                          typename vecCore::Mask_v<Real_v> &completelyoutside) const
0799   {
0800     if (ForInside) completelyinside = IsCompletelyInside<Real_v>(localPoint);
0801 
0802     completelyoutside = IsCompletelyOutside<Real_v>(localPoint);
0803   }
0804 
0805 }; // end of class ThetaCone
0806 } // namespace VECGEOM_IMPL_NAMESPACE
0807 } // namespace vecgeom
0808 
0809 #endif /* VECGEOM_VOLUMES_THETACONE_H_ */
[ 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!