VecGeom/volumes/ThetaCone.h

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