Back to home page

EIC code displayed by LXR
 
 

    

File indexing completed on 2026-05-17 08:25:52 

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