Back to home page

EIC code displayed by LXR
 
 

    

File indexing completed on 2026-11-13 11:00:37 

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