Back to home page

EIC code displayed by LXR
 
 

    

File indexing completed on 2025-01-30 10:26:24 

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 
0254       using Bool_v       = vecCore::Mask_v<Real_v>;
0255       using InsideBool_v = vecCore::Mask_v<Inside_t>;
0256       Bool_v completelyinside, completelyoutside;
0257       GenericKernelForContainsAndInside<Real_v, true>(point, completelyinside, completelyoutside);
0258       Inside_t inside(EInside::kSurface);
0259       vecCore::MaskedAssign(inside, (InsideBool_v)completelyoutside, Inside_t(EInside::kOutside));
0260       vecCore::MaskedAssign(inside, (InsideBool_v)completelyinside, Inside_t(EInside::kInside));
0261       return inside;
0262     }
0263   */
0264   template <typename Real_v, typename Inside_t>
0265   VECCORE_ATT_HOST_DEVICE
0266   Inside_t Inside(Vector3D<Real_v> const &point) const
0267   {
0268     using Bool_v       = vecCore::Mask_v<Real_v>;
0269     using InsideBool_v = vecCore::Mask_v<Inside_t>;
0270     Bool_v completelyinside, completelyoutside;
0271     GenericKernelForContainsAndInside<Real_v, true>(point, completelyinside, completelyoutside);
0272     Inside_t inside(EInside::kSurface);
0273     vecCore::MaskedAssign(inside, (InsideBool_v)completelyoutside, Inside_t(EInside::kOutside));
0274     vecCore::MaskedAssign(inside, (InsideBool_v)completelyinside, Inside_t(EInside::kInside));
0275     return inside;
0276   }
0277 
0278   /**
0279    * estimate of the smallest distance to the ThetaCone boundary when
0280    * the point is located outside the ThetaCone
0281    */
0282   template <typename Real_v>
0283   VECCORE_ATT_HOST_DEVICE
0284   Real_v SafetyToIn(Vector3D<Real_v> const &point) const
0285   {
0286 
0287     using Bool_v = vecCore::Mask_v<Real_v>;
0288 
0289     Real_v safeTheta(0.);
0290     Real_v pointRad = Sqrt(point.x() * point.x() + point.y() * point.y());
0291     Real_v sfTh1    = DistanceToLine<Real_v>(slope1, pointRad, point.z());
0292     Real_v sfTh2    = DistanceToLine<Real_v>(slope2, pointRad, point.z());
0293 
0294     safeTheta   = Min(sfTh1, sfTh2);
0295     Bool_v done = Contains<Real_v>(point);
0296     vecCore__MaskedAssignFunc(safeTheta, done, Real_v(0.));
0297     if (vecCore::MaskFull(done)) return safeTheta;
0298 
0299     // Case 1 : Both cones are in Positive Z direction
0300     if (fSTheta < kPi / 2 + halfAngTolerance) {
0301       if (fETheta < kPi / 2 + halfAngTolerance) {
0302         if (fSTheta < fETheta) {
0303           vecCore::MaskedAssign(safeTheta, (!done && point.z() < Real_v(0.)), sfTh2);
0304         }
0305       }
0306 
0307       // Case 2 : First Cone is in Positive Z direction and Second is in Negative Z direction
0308       if (fETheta > kPi / 2 + halfAngTolerance) {
0309         if (fSTheta < fETheta) {
0310           vecCore::MaskedAssign(safeTheta, (!done && point.z() > Real_v(0.)), sfTh1);
0311           vecCore::MaskedAssign(safeTheta, (!done && point.z() < Real_v(0.)), sfTh2);
0312         }
0313       }
0314     }
0315 
0316     // Case 3 : Both cones are in Negative Z direction
0317     if (fETheta > kPi / 2 + halfAngTolerance) {
0318       if (fSTheta > kPi / 2 + halfAngTolerance) {
0319         if (fSTheta < fETheta) {
0320           vecCore::MaskedAssign(safeTheta, (!done && point.z() > Real_v(0.)), sfTh1);
0321         }
0322       }
0323     }
0324 
0325     return safeTheta;
0326   }
0327 
0328   /**
0329    * estimate of the smallest distance to the ThetaCone boundary when
0330    * the point is located inside the ThetaCone ( within the defining phi angle )
0331    */
0332   template <typename Real_v>
0333   VECCORE_ATT_HOST_DEVICE
0334   Real_v SafetyToOut(Vector3D<Real_v> const &point) const
0335   {
0336 
0337     Real_v pointRad    = Sqrt(point.x() * point.x() + point.y() * point.y());
0338     Real_v bisectorRad = Abs(point.z() * tanBisector);
0339 
0340     vecCore::Mask<Real_v> condition(false);
0341     Real_v sfTh1 = DistanceToLine<Real_v>(slope1, pointRad, point.z());
0342     Real_v sfTh2 = DistanceToLine<Real_v>(slope2, pointRad, point.z());
0343 
0344     // Case 1 : Both cones are in Positive Z direction
0345     if (fSTheta < kPi / 2 + halfAngTolerance) {
0346       if (fETheta < kPi / 2 + halfAngTolerance) {
0347         if (fSTheta < fETheta) {
0348           condition = (pointRad < bisectorRad) && (fSTheta != Real_v(0.));
0349         }
0350       }
0351 
0352       // Case 2 : First Cone is in Positive Z direction and Second is in Negative Z direction
0353       if (fETheta > kPi / 2 + halfAngTolerance) {
0354         if (fSTheta < fETheta) {
0355           condition = sfTh1 < sfTh2;
0356         }
0357       }
0358     }
0359 
0360     // Case 3 : Both cones are in Negative Z direction
0361     if (fETheta > kPi / 2 + halfAngTolerance) {
0362       if (fSTheta > kPi / 2 + halfAngTolerance) {
0363         if (fSTheta < fETheta) {
0364           condition = !((pointRad < bisectorRad) && (fETheta != Real_v(kPi)));
0365         }
0366       }
0367     }
0368 
0369     return vecCore::Blend(condition, sfTh1, sfTh2);
0370   }
0371 
0372   template <typename Real_v>
0373   VECCORE_ATT_HOST_DEVICE
0374   Real_v DistanceToLine(Precision const &slope, Real_v const &x, Real_v const &y) const
0375   {
0376 
0377     Real_v dist = (y - slope * x) / Sqrt(Real_v(1.) + slope * slope);
0378     return Abs(dist);
0379   }
0380 
0381   /**
0382    * estimate of the distance to the ThetaCone boundary with given direction
0383    */
0384   template <typename Real_v>
0385   VECCORE_ATT_HOST_DEVICE
0386   void DistanceToIn(Vector3D<Real_v> const &point, Vector3D<Real_v> const &dir, Real_v &distThetaCone1,
0387                     Real_v &distThetaCone2, typename vecCore::Mask_v<Real_v> &intsect1,
0388                     typename vecCore::Mask_v<Real_v> &intsect2) const
0389   {
0390 
0391     {
0392       using Bool_v = vecCore::Mask_v<Real_v>;
0393 
0394       Bool_v done(false);
0395       Bool_v fal(false);
0396 
0397       Real_v firstRoot(kInfLength), secondRoot(kInfLength);
0398 
0399       Real_v pDotV2d = point.x() * dir.x() + point.y() * dir.y();
0400       Real_v rho2    = point.x() * point.x() + point.y() * point.y();
0401       Real_v dirRho2 = dir.Perp2();
0402 
0403       Real_v b = pDotV2d - point.z() * dir.z() * tanSTheta2;
0404       // Real_v a = dir.x() * dir.x() + dir.y() * dir.y() - dir.z() * dir.z() * tanSTheta2;
0405       Real_v a    = dirRho2 - dir.z() * dir.z() * tanSTheta2;
0406       Real_v c    = rho2 - point.z() * point.z() * tanSTheta2;
0407       Real_v d2   = b * b - a * c;
0408       Real_v aInv = Real_v(1.) / NonZero(a);
0409 
0410       vecCore__MaskedAssignFunc(firstRoot, (d2 > Real_v(0.)), (-b + Sqrt(Abs(d2))) * aInv);
0411       done |= (Abs(firstRoot) < Real_v(3.) * kTolerance);
0412       vecCore__MaskedAssignFunc(firstRoot, ((Abs(firstRoot) < Real_v(3.) * kTolerance)), Real_v(0.));
0413       vecCore__MaskedAssignFunc(firstRoot, (!done && (firstRoot < Real_v(0.))), InfinityLength<Real_v>());
0414 
0415       Real_v b2 = pDotV2d - point.z() * dir.z() * tanETheta2;
0416       // Real_v a2 = dir.x() * dir.x() + dir.y() * dir.y() - dir.z() * dir.z() * tanETheta2;
0417       Real_v a2    = dirRho2 - dir.z() * dir.z() * tanETheta2;
0418       Real_v c2    = rho2 - point.z() * point.z() * tanETheta2;
0419       Real_v d22   = b2 * b2 - a2 * c2;
0420       Real_v a2Inv = Real_v(1.) / NonZero(a2);
0421 
0422       vecCore__MaskedAssignFunc(secondRoot, (d22 > Real_v(0.)), (-b2 - Sqrt(Abs(d22))) * a2Inv);
0423       vecCore__MaskedAssignFunc(secondRoot, (!done && (Abs(secondRoot) < Real_v(3.) * kTolerance)), Real_v(0.));
0424       done |= (Abs(secondRoot) < Real_v(3.) * kTolerance);
0425       vecCore__MaskedAssignFunc(secondRoot, !done && (secondRoot < Real_v(0.)), InfinityLength<Real_v>());
0426 
0427       if (fSTheta < kHalfPi + halfAngTolerance) {
0428         if (fETheta < kHalfPi + halfAngTolerance) {
0429           if (fSTheta < fETheta) {
0430             distThetaCone1          = firstRoot;
0431             distThetaCone2          = secondRoot;
0432             Real_v zOfIntSecPtCone1 = (point.z() + distThetaCone1 * dir.z());
0433             Real_v zOfIntSecPtCone2 = (point.z() + distThetaCone2 * dir.z());
0434 
0435             intsect1 = ((d2 > Real_v(0.)) && (zOfIntSecPtCone1 > Real_v(0.)));
0436             intsect2 = ((d22 > Real_v(0.)) && (zOfIntSecPtCone2 > Real_v(0.)));
0437           }
0438         }
0439 
0440         if (fETheta >= kHalfPi - halfAngTolerance && fETheta <= kHalfPi + halfAngTolerance) {
0441           vecCore__MaskedAssignFunc(distThetaCone2, (dir.z() > Real_v(0.)), -point.z() / dir.z());
0442           Real_v zOfIntSecPtCone2 = (point.z() + distThetaCone2 * dir.z());
0443           intsect2                = ((distThetaCone2 != kInfLength) && (Abs(zOfIntSecPtCone2) < halfAngTolerance));
0444         }
0445 
0446         if (fETheta > kHalfPi + halfAngTolerance) {
0447           if (fSTheta < fETheta) {
0448             distThetaCone1 = firstRoot;
0449             vecCore__MaskedAssignFunc(secondRoot, (d22 > Real_v(0.)), (-b2 + Sqrt(Abs(d22))) * a2Inv);
0450 
0451             done = fal;
0452             done |= (Abs(secondRoot) < Real_v(3.) * kTolerance);
0453             vecCore__MaskedAssignFunc(secondRoot, ((Abs(secondRoot) < Real_v(3.) * kTolerance)), Real_v(0.));
0454             vecCore__MaskedAssignFunc(secondRoot, !done && (secondRoot < Real_v(0.)), InfinityLength<Real_v>());
0455             distThetaCone2 = secondRoot;
0456 
0457             Real_v zOfIntSecPtCone1 = (point.z() + distThetaCone1 * dir.z());
0458             Real_v zOfIntSecPtCone2 = (point.z() + distThetaCone2 * dir.z());
0459 
0460             intsect1 = ((d2 > 0) && (distThetaCone1 != kInfLength) && (zOfIntSecPtCone1 > Real_v(0.)));
0461             intsect2 = ((d22 > 0) && (distThetaCone2 != kInfLength) && (zOfIntSecPtCone2 < Real_v(0.)));
0462           }
0463         }
0464       }
0465 
0466       if (fSTheta >= kHalfPi - halfAngTolerance) {
0467         if (fETheta > kHalfPi + halfAngTolerance) {
0468           if (fSTheta < fETheta) {
0469             vecCore__MaskedAssignFunc(firstRoot, (d2 > Real_v(0.)), (-b - Sqrt(Abs(d2))) * aInv);
0470             done = fal;
0471             done |= (Abs(firstRoot) < Real_v(3.) * kTolerance);
0472             vecCore__MaskedAssignFunc(firstRoot, ((Abs(firstRoot) < Real_v(3.) * kTolerance)), Real_v(0.));
0473             vecCore__MaskedAssignFunc(firstRoot, !done && (firstRoot < Real_v(0.)), InfinityLength<Real_v>());
0474             distThetaCone1 = firstRoot;
0475 
0476             vecCore__MaskedAssignFunc(secondRoot, (d22 > Real_v(0.)), (-b2 + Sqrt(Abs(d22))) * a2Inv);
0477             done = fal;
0478             done |= (Abs(secondRoot) < Real_v(3.) * kTolerance);
0479             vecCore__MaskedAssignFunc(secondRoot, ((Abs(secondRoot) < Real_v(3.) * kTolerance)), Real_v(0.));
0480             vecCore__MaskedAssignFunc(secondRoot, !done && (secondRoot < Real_v(0.)), InfinityLength<Real_v>());
0481             distThetaCone2 = secondRoot;
0482 
0483             Real_v zOfIntSecPtCone1 = (point.z() + distThetaCone1 * dir.z());
0484             Real_v zOfIntSecPtCone2 = (point.z() + distThetaCone2 * dir.z());
0485 
0486             intsect1 = ((d2 > 0) && (distThetaCone1 != kInfLength) && (zOfIntSecPtCone1 < Real_v(0.)));
0487             intsect2 = ((d22 > 0) && (distThetaCone2 != kInfLength) && (zOfIntSecPtCone2 < Real_v(0.)));
0488           }
0489         }
0490       }
0491 
0492       if (fSTheta >= kHalfPi - halfAngTolerance && fSTheta <= kHalfPi + halfAngTolerance) {
0493         vecCore__MaskedAssignFunc(distThetaCone1, (dir.z() < Real_v(0.)), -point.z() / dir.z());
0494         Real_v zOfIntSecPtCone1 = (point.z() + distThetaCone1 * dir.z());
0495         intsect1                = ((distThetaCone1 != kInfLength) && (Abs(zOfIntSecPtCone1) < halfAngTolerance));
0496       }
0497     }
0498   }
0499 
0500   template <typename Real_v>
0501   VECCORE_ATT_HOST_DEVICE
0502   void DistanceToOut(Vector3D<Real_v> const &point, Vector3D<Real_v> const &dir, Real_v &distThetaCone1,
0503                      Real_v &distThetaCone2, typename vecCore::Mask_v<Real_v> &intsect1,
0504                      typename vecCore::Mask_v<Real_v> &intsect2) const
0505   {
0506 
0507     using Bool_v = vecCore::Mask_v<Real_v>;
0508 
0509     Real_v inf(kInfLength);
0510 
0511     // Real_v firstRoot(kInfLength), secondRoot(kInfLength);
0512     Real_v firstRoot(kInfLength), secondRoot(kInfLength);
0513 
0514     Real_v pDotV2d = point.x() * dir.x() + point.y() * dir.y();
0515     Real_v rho2    = point.x() * point.x() + point.y() * point.y();
0516 
0517     Real_v b  = pDotV2d - point.z() * dir.z() * tanSTheta2;
0518     Real_v a  = dir.x() * dir.x() + dir.y() * dir.y() - dir.z() * dir.z() * tanSTheta2;
0519     Real_v c  = rho2 - point.z() * point.z() * tanSTheta2;
0520     Real_v d2 = b * b - a * c;
0521     vecCore__MaskedAssignFunc(d2, d2 < Real_v(0.) && Abs(d2) < kHalfTolerance, Real_v(0.));
0522     vecCore__MaskedAssignFunc(firstRoot, (d2 >= Real_v(0.)) && b >= Real_v(0.) && a != Real_v(0.),
0523                               ((-b - Sqrt(Abs(d2))) / NonZero(a)));
0524     vecCore__MaskedAssignFunc(firstRoot, (d2 >= Real_v(0.)) && b < Real_v(0.), ((c) / NonZero(-b + Sqrt(Abs(d2)))));
0525     vecCore__MaskedAssignFunc(firstRoot, firstRoot < Real_v(0.), InfinityLength<Real_v>());
0526 
0527     Real_v b2 = point.x() * dir.x() + point.y() * dir.y() - point.z() * dir.z() * tanETheta2;
0528     Real_v a2 = dir.x() * dir.x() + dir.y() * dir.y() - dir.z() * dir.z() * tanETheta2;
0529 
0530     Real_v c2  = point.x() * point.x() + point.y() * point.y() - point.z() * point.z() * tanETheta2;
0531     Real_v d22 = (b2 * b2) - (a2 * c2);
0532     vecCore__MaskedAssignFunc(d22, d22 < Real_v(0.) && Abs(d22) < kHalfTolerance, Real_v(0.));
0533 
0534     vecCore__MaskedAssignFunc(secondRoot, (d22 >= Real_v(0.)) && b2 >= Real_v(0.),
0535                               ((c2) / NonZero(-b2 - Sqrt(Abs(d22)))));
0536     vecCore__MaskedAssignFunc(secondRoot, (d22 >= Real_v(0.)) && b2 < Real_v(0.) && a2 != Real_v(0.),
0537                               ((-b2 + Sqrt(Abs(d22))) / NonZero(a2)));
0538 
0539     vecCore__MaskedAssignFunc(secondRoot, secondRoot < Real_v(0.) && Abs(secondRoot) > kTolerance,
0540                               InfinityLength<Real_v>());
0541     vecCore__MaskedAssignFunc(secondRoot, Abs(secondRoot) < kTolerance, Real_v(0.));
0542 
0543     if (fSTheta < kPi / 2 + halfAngTolerance) {
0544       if (fETheta < kPi / 2 + halfAngTolerance) {
0545         if (fSTheta < fETheta) {
0546           distThetaCone1          = firstRoot;
0547           distThetaCone2          = secondRoot;
0548           Real_v zOfIntSecPtCone1 = (point.z() + distThetaCone1 * dir.z());
0549           Real_v zOfIntSecPtCone2 = (point.z() + distThetaCone2 * dir.z());
0550 
0551           intsect1 = ((d2 > 0) && (distThetaCone1 != kInfLength) && ((zOfIntSecPtCone1) > -kHalfTolerance));
0552           intsect2 = ((d22 > 0) && (distThetaCone2 != kInfLength) && ((zOfIntSecPtCone2) > -kHalfTolerance));
0553 
0554           Real_v dirRho2 = dir.x() * dir.x() + dir.y() * dir.y();
0555           Real_v zs(kInfLength);
0556           if (fSTheta) zs = dirRho2 / tanSTheta;
0557           Real_v ze(kInfLength);
0558           if (fETheta) ze = dirRho2 / tanETheta;
0559           Bool_v cond = (point.x() == Real_v(0.) && point.y() == Real_v(0.) && point.z() == Real_v(0.) &&
0560                          dir.z() < zs && dir.z() < ze);
0561           vecCore__MaskedAssignFunc(distThetaCone1, cond, Real_v(0.));
0562           vecCore__MaskedAssignFunc(distThetaCone2, cond, Real_v(0.));
0563           intsect1 |= cond;
0564           intsect2 |= cond;
0565         }
0566       }
0567 
0568       // Second cone is the XY plane, compute distance to plane
0569       if (fETheta >= kPi / 2 - halfAngTolerance && fETheta <= kPi / 2 + halfAngTolerance) {
0570         distThetaCone1 = firstRoot;
0571         distThetaCone2 = inf;
0572         vecCore__MaskedAssignFunc(distThetaCone2, (dir.z() < Real_v(0.)), Real_v(-1.) * point.z() / dir.z());
0573         Real_v zOfIntSecPtCone2 = (point.z() + distThetaCone2 * dir.z());
0574         Real_v zOfIntSecPtCone1 = (point.z() + distThetaCone1 * dir.z());
0575         intsect2 =
0576             ((distThetaCone2 != kInfLength) && (Abs(zOfIntSecPtCone2) < kHalfTolerance) && !(dir.z() == Real_v(0.)));
0577         intsect1 = ((d2 >= 0) && (distThetaCone1 != kInfLength) && (Abs(zOfIntSecPtCone1) < kHalfTolerance) &&
0578                     !(dir.z() == Real_v(0.)));
0579       }
0580 
0581       if (fETheta > kPi / 2 + halfAngTolerance) {
0582         if (fSTheta < fETheta) {
0583           distThetaCone1 = firstRoot;
0584           vecCore__MaskedAssignFunc(secondRoot, (d22 >= Real_v(0.)) && b2 > Real_v(0.) && a2 != Real_v(0.),
0585                                     ((-b2 - Sqrt(Abs(d22))) / NonZero(a2)));
0586           vecCore__MaskedAssignFunc(secondRoot, (d22 >= Real_v(0.)) && b2 <= Real_v(0.),
0587                                     ((c2) / NonZero(-b2 + Sqrt(Abs(d22)))));
0588           vecCore__MaskedAssignFunc(secondRoot, secondRoot < Real_v(0.), InfinityLength<Real_v>());
0589           distThetaCone2          = secondRoot;
0590           Real_v zOfIntSecPtCone1 = (point.z() + distThetaCone1 * dir.z());
0591           Real_v zOfIntSecPtCone2 = (point.z() + distThetaCone2 * dir.z());
0592           intsect1 = ((d2 >= 0) && (distThetaCone1 != kInfLength) && ((zOfIntSecPtCone1) > -kHalfTolerance));
0593           intsect2 = ((d22 >= 0) && (distThetaCone2 != kInfLength) && ((zOfIntSecPtCone2) < kHalfTolerance));
0594         }
0595       }
0596     }
0597 
0598     if (fETheta > kPi / 2 + halfAngTolerance) {
0599       if (fSTheta < fETheta) {
0600         secondRoot = kInfLength;
0601         vecCore__MaskedAssignFunc(secondRoot, (d22 >= Real_v(0.)) && b2 > Real_v(0.) && a2 != Real_v(0.),
0602                                   ((-b2 - Sqrt(Abs(d22))) / NonZero(a2)));
0603         vecCore__MaskedAssignFunc(secondRoot, (d22 >= Real_v(0.)) && b2 <= Real_v(0.),
0604                                   ((c2) / NonZero(-b2 + Sqrt(Abs(d22)))));
0605         distThetaCone2 = secondRoot;
0606 
0607         if (fSTheta > kPi / 2 + halfAngTolerance) {
0608           vecCore__MaskedAssignFunc(firstRoot, (d2 >= Real_v(0.)) && b > Real_v(0.),
0609                                     ((c) / NonZero(-b - Sqrt(Abs(d2)))));
0610           vecCore__MaskedAssignFunc(firstRoot, (d2 >= Real_v(0.)) && b <= Real_v(0.) && a != Real_v(0.),
0611                                     ((-b + Sqrt(Abs(d2))) / NonZero(a)));
0612           vecCore__MaskedAssignFunc(firstRoot, firstRoot < Real_v(0.), InfinityLength<Real_v>());
0613           distThetaCone1          = firstRoot;
0614           Real_v zOfIntSecPtCone1 = (point.z() + distThetaCone1 * dir.z());
0615           intsect1 = ((d2 > 0) && (distThetaCone1 != kInfLength) && ((zOfIntSecPtCone1) < kHalfTolerance));
0616           Real_v zOfIntSecPtCone2 = (point.z() + distThetaCone2 * dir.z());
0617           intsect2 = ((d22 > 0) && (distThetaCone2 != kInfLength) && ((zOfIntSecPtCone2) < kHalfTolerance));
0618 
0619           Real_v dirRho2 = dir.x() * dir.x() + dir.y() * dir.y();
0620           Real_v zs(-kInfLength);
0621           if (tanSTheta) zs = -dirRho2 / tanSTheta;
0622           Real_v ze(-kInfLength);
0623           if (tanETheta) ze = -dirRho2 / tanETheta;
0624           Bool_v cond = (point.x() == Real_v(0.) && point.y() == Real_v(0.) && point.z() == Real_v(0.) &&
0625                          dir.z() > zs && dir.z() > ze);
0626           vecCore__MaskedAssignFunc(distThetaCone1, cond, Real_v(0.));
0627           vecCore__MaskedAssignFunc(distThetaCone2, cond, Real_v(0.));
0628           // intsect1 |= (cond && tr);
0629           // intsect2 |= (cond && tr);
0630           intsect1 |= cond;
0631           intsect2 |= cond;
0632         }
0633       }
0634     }
0635 
0636     // First cone is the XY plane, compute distance to plane
0637     if (fSTheta >= kPi / 2 - halfAngTolerance && fSTheta <= kPi / 2 + halfAngTolerance) {
0638       distThetaCone2 = secondRoot;
0639       distThetaCone1 = kInfLength;
0640       vecCore__MaskedAssignFunc(distThetaCone1, (dir.z() > Real_v(0.)), Real_v(-1.) * point.z() / NonZero(dir.z()));
0641       Real_v zOfIntSecPtCone1 = (point.z() + distThetaCone1 * dir.z());
0642 
0643       Real_v zOfIntSecPtCone2 = (point.z() + distThetaCone2 * dir.z());
0644 
0645       intsect1 =
0646           ((distThetaCone1 != kInfLength) && (Abs(zOfIntSecPtCone1) < kHalfTolerance) && (dir.z() != Real_v(0.)));
0647       intsect2 = ((d22 >= 0) && (distThetaCone2 != kInfLength) && (Abs(zOfIntSecPtCone2) < kHalfTolerance) &&
0648                   (dir.z() != Real_v(0.)));
0649     }
0650   }
0651 
0652   // This could be useful in case somebody just want to check whether point is completely inside ThetaRange
0653   template <typename Real_v>
0654   VECCORE_ATT_HOST_DEVICE
0655   typename vecCore::Mask_v<Real_v> IsCompletelyInside(Vector3D<Real_v> const &localPoint) const
0656   {
0657 
0658     using Bool_v       = vecCore::Mask_v<Real_v>;
0659     Real_v rho         = Sqrt(localPoint.Mag2() - (localPoint.z() * localPoint.z()));
0660     Real_v cone1Radius = Abs(localPoint.z() * tanSTheta);
0661     Real_v cone2Radius = Abs(localPoint.z() * tanETheta);
0662     Bool_v isPointOnZAxis =
0663         localPoint.z() != Real_v(0.) && localPoint.x() == Real_v(0.) && localPoint.y() == Real_v(0.);
0664 
0665     Bool_v isPointOnXYPlane =
0666         localPoint.z() == Real_v(0.) && (localPoint.x() != Real_v(0.) || localPoint.y() != Real_v(0.));
0667 
0668     Real_v startTheta(fSTheta), endTheta(fETheta);
0669 
0670     Bool_v completelyinside = (isPointOnZAxis && ((startTheta == Real_v(0.) && endTheta == kPi) ||
0671                                                   (localPoint.z() > Real_v(0.) && startTheta == Real_v(0.)) ||
0672                                                   (localPoint.z() < Real_v(0.) && endTheta == kPi)));
0673 
0674     completelyinside |=
0675         (!completelyinside && (isPointOnXYPlane && (startTheta < Real_v(kHalfPi) && endTheta > Real_v(kHalfPi) &&
0676                                                     (Real_v(kHalfPi) - startTheta) > kAngTolerance &&
0677                                                     (endTheta - Real_v(kHalfPi)) > kTolerance)));
0678 
0679     if (fSTheta < kHalfPi + halfAngTolerance) {
0680       if (fETheta < kHalfPi + halfAngTolerance) {
0681         if (fSTheta < fETheta) {
0682           Real_v tolAngMin = cone1Radius + Real_v(2 * kAngTolerance * 10.);
0683           Real_v tolAngMax = cone2Radius - Real_v(2 * kAngTolerance * 10.);
0684 
0685           completelyinside |=
0686               (!completelyinside &&
0687                (((rho <= tolAngMax) && (rho >= tolAngMin) && (localPoint.z() > Real_v(0.)) &&
0688                  Bool_v(fSTheta != Real_v(0.))) ||
0689                 ((rho <= tolAngMax) && Bool_v(fSTheta == Real_v(0.)) && (localPoint.z() > Real_v(0.)))));
0690         }
0691       }
0692 
0693       if (fETheta > kHalfPi + halfAngTolerance) {
0694         if (fSTheta < fETheta) {
0695           Real_v tolAngMin = cone1Radius + Real_v(2 * kAngTolerance * 10.);
0696           Real_v tolAngMax = cone2Radius + Real_v(2 * kAngTolerance * 10.);
0697 
0698           completelyinside |= (!completelyinside && (((rho >= tolAngMin) && (localPoint.z() > Real_v(0.))) ||
0699                                                      ((rho >= tolAngMax) && (localPoint.z() < Real_v(0.)))));
0700         }
0701       }
0702 
0703       if (fETheta >= kHalfPi - halfAngTolerance && fETheta <= kHalfPi + halfAngTolerance) {
0704 
0705         completelyinside &= !(Abs(localPoint.z()) < halfAngTolerance);
0706       }
0707     }
0708 
0709     if (fETheta > kHalfPi + halfAngTolerance) {
0710       if (fSTheta >= kHalfPi - halfAngTolerance && fSTheta <= kHalfPi + halfAngTolerance) {
0711 
0712         completelyinside &= !(Abs(localPoint.z()) < halfAngTolerance);
0713       }
0714 
0715       if (fSTheta > kHalfPi + halfAngTolerance) {
0716         if (fSTheta < fETheta) {
0717           Real_v tolAngMin = cone1Radius - Real_v(2 * kAngTolerance * 10.);
0718           Real_v tolAngMax = cone2Radius + Real_v(2 * kAngTolerance * 10.);
0719 
0720           completelyinside |=
0721               (!completelyinside &&
0722                (((rho <= tolAngMin) && (rho >= tolAngMax) && (localPoint.z() < Real_v(0.)) && Bool_v(fETheta != kPi)) ||
0723                 ((rho <= tolAngMin) && (localPoint.z() < Real_v(0.)) && Bool_v(fETheta == kPi))));
0724         }
0725       }
0726     }
0727     return completelyinside;
0728   }
0729 
0730   // This could be useful in case somebody just want to check whether point is completely outside ThetaRange
0731   template <typename Real_v>
0732   VECCORE_ATT_HOST_DEVICE
0733   typename vecCore::Mask_v<Real_v> IsCompletelyOutside(Vector3D<Real_v> const &localPoint) const
0734   {
0735 
0736     using Bool_v       = vecCore::Mask_v<Real_v>;
0737     Real_v diff        = localPoint.Perp2();
0738     Real_v rho         = Sqrt(diff);
0739     Real_v cone1Radius = Abs(localPoint.z() * tanSTheta);
0740     Real_v cone2Radius = Abs(localPoint.z() * tanETheta);
0741     Bool_v isPointOnZAxis =
0742         localPoint.z() != Real_v(0.) && localPoint.x() == Real_v(0.) && localPoint.y() == Real_v(0.);
0743 
0744     Bool_v isPointOnXYPlane =
0745         localPoint.z() == Real_v(0.) && (localPoint.x() != Real_v(0.) || localPoint.y() != Real_v(0.));
0746 
0747     // Real_v startTheta(fSTheta), endTheta(fETheta);
0748 
0749     Bool_v completelyoutside = (isPointOnZAxis && ((Bool_v(fSTheta != Real_v(0.)) && Bool_v(fETheta != kPi)) ||
0750                                                    (localPoint.z() > Real_v(0.) && Bool_v(fSTheta != Real_v(0.))) ||
0751                                                    (localPoint.z() < Real_v(0.) && Bool_v(fETheta != kPi))));
0752 
0753     completelyoutside |=
0754         (!completelyoutside &&
0755          (isPointOnXYPlane && Bool_v(((fSTheta < kHalfPi) && (fETheta < kHalfPi) &&
0756                                       ((kHalfPi - fSTheta) > kAngTolerance) && ((kHalfPi - fETheta) > kTolerance)) ||
0757                                      ((fSTheta > kHalfPi && fETheta > kHalfPi) &&
0758                                       ((fSTheta - kHalfPi) > kAngTolerance) && ((fETheta - kHalfPi) > kTolerance)))));
0759 
0760     if (fSTheta < kHalfPi + halfAngTolerance) {
0761       if (fETheta < kHalfPi + halfAngTolerance) {
0762         if (fSTheta < fETheta) {
0763 
0764           Real_v tolAngMin2 = cone1Radius - Real_v(2 * kAngTolerance * 10.);
0765           Real_v tolAngMax2 = cone2Radius + Real_v(2 * kAngTolerance * 10.);
0766 
0767           completelyoutside |=
0768               (!completelyoutside && ((rho < tolAngMin2) || (rho > tolAngMax2) || (localPoint.z() < Real_v(0.))));
0769         }
0770       }
0771 
0772       if (fETheta > kHalfPi + halfAngTolerance) {
0773         if (fSTheta < fETheta) {
0774           Real_v tolAngMin2 = cone1Radius - Real_v(2 * kAngTolerance * 10.);
0775           Real_v tolAngMax2 = cone2Radius - Real_v(2 * kAngTolerance * 10.);
0776           completelyoutside |= (!completelyoutside && (((rho < tolAngMin2) && (localPoint.z() > Real_v(0.))) ||
0777                                                        ((rho < tolAngMax2) && (localPoint.z() < Real_v(0.)))));
0778         }
0779       }
0780 
0781       if (fETheta >= kHalfPi - halfAngTolerance && fETheta <= kHalfPi + halfAngTolerance) {
0782 
0783         // completelyinside &= !(Abs(localPoint.z()) < halfAngTolerance);
0784         completelyoutside &= !(Abs(localPoint.z()) < halfAngTolerance);
0785       }
0786     }
0787 
0788     if (fETheta > kHalfPi + halfAngTolerance) {
0789       if (fSTheta >= kHalfPi - halfAngTolerance && fSTheta <= kHalfPi + halfAngTolerance) {
0790 
0791         // completelyinside &= !(Abs(localPoint.z()) < halfAngTolerance);
0792         completelyoutside &= !(Abs(localPoint.z()) < halfAngTolerance);
0793       }
0794 
0795       if (fSTheta > kHalfPi + halfAngTolerance) {
0796         if (fSTheta < fETheta) {
0797 
0798           Real_v tolAngMin2 = cone1Radius + Real_v(2 * kAngTolerance * 10.);
0799           Real_v tolAngMax2 = cone2Radius - Real_v(2 * kAngTolerance * 10.);
0800           completelyoutside |=
0801               (!completelyoutside && ((rho < tolAngMax2) || (rho > tolAngMin2) || (localPoint.z() > Real_v(0.))));
0802         }
0803       }
0804     }
0805 
0806     return completelyoutside;
0807   }
0808 
0809   template <typename Real_v, bool ForInside>
0810   VECCORE_ATT_HOST_DEVICE
0811   void GenericKernelForContainsAndInside(Vector3D<Real_v> const &localPoint,
0812                                          typename vecCore::Mask_v<Real_v> &completelyinside,
0813                                          typename vecCore::Mask_v<Real_v> &completelyoutside) const
0814   {
0815     if (ForInside) completelyinside = IsCompletelyInside<Real_v>(localPoint);
0816 
0817     completelyoutside = IsCompletelyOutside<Real_v>(localPoint);
0818   }
0819 
0820 }; // end of class ThetaCone
0821 } // namespace VECGEOM_IMPL_NAMESPACE
0822 } // namespace vecgeom
0823 
0824 #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!