surfaces/surf/ConicalImpl.h

0001 #ifndef VECGEOM_CONICAL_IMPL_H
0002 #define VECGEOM_CONICAL_IMPL_H
0003
0004 #include <VecGeom/surfaces/surf/SurfaceHelper.h>
0005 #include <VecGeom/surfaces/base/Equations.h>
0006
0007 namespace vgbrep {
0008
0009 template <typename Real_t>
0010 struct SurfaceHelper<SurfaceType::kConical, Real_t> {
0011
0012   VECGEOM_FORCE_INLINE
0013   VECCORE_ATT_HOST_DEVICE
0014   SurfaceHelper() {}
0015
0016   VECGEOM_FORCE_INLINE
0017   VECCORE_ATT_HOST_DEVICE
0018   /// @brief Inside half-space function
0019   /// @param point Point in local surface coordinates
0020   /// @param tol tolerance for determining if a point is inside or not
0021   /// @return True if the point is behind the normal within kTolerance (surface is included)
0022   bool Inside(Vector3D<Real_t> const &point, bool flip, Real_t radius, Real_t slope,
0023               Real_t tol = vecgeom::kToleranceStrict<Real_t>)
0024   {
0025     int flipsign      = (radius < 0) ? -1 : 1;
0026     int bool_flipsign = !flip ? 1 : -1;
0027     Real_t coneR      = Abs(radius) + point.z() * slope;
0028     Real_t rho        = point.Perp();
0029     return flipsign * (rho - coneR) < bool_flipsign * tol;
0030   }
0031
0032   VECGEOM_FORCE_INLINE
0033   VECCORE_ATT_HOST_DEVICE
0034   /// @brief Find signed distance to next intersection from local point.
0035   /// @param point Point in local surface coordinates
0036   /// @param dir Direction in the local surface coordinates
0037   /// @param left_side Flag specifying if the surface is intersected from the left-side that defines the normal
0038   /// @param distance Computed distance to surface
0039   /// @return Validity of the intersection
0040   bool Intersect(Vector3D<Real_t> const &point, Vector3D<Real_t> const &dir, bool left_side, Real_t &distance,
0041                  bool &two_solutions, Real_t &safety, Real_t radius, Real_t slope)
0042   {
0043     QuadraticCoef<Real_t> coef;
0044     Real_t roots[2];
0045     int numroots      = 0;
0046     bool flip_exiting = left_side ^ (radius < 0);
0047     ConeEq<Real_t>(point, dir, Abs(radius), slope, coef);
0048     QuadraticSolver(coef, roots, numroots);
0049     two_solutions = (numroots == 2 && roots[0] > -vecgeom::kToleranceStrict<Real_t> &&
0050                      roots[1] > -vecgeom::kToleranceStrict<Real_t>);
0051     for (auto i = 0; i < numroots; ++i) {
0052       distance                = roots[i];
0053       Vector3D<Real_t> onsurf = point + distance * dir;
0054       // Exclude solutions beyond the tip of the cone. What if the tip is included? TODO
0055       if (Abs(radius) + onsurf[2] * slope < 0) continue;
0056       Vector3D<Real_t> normal(onsurf[0], onsurf[1], -std::sqrt(onsurf[0] * onsurf[0] + onsurf[1] * onsurf[1]) * slope);
0057       bool hit = flip_exiting ^ (dir.Dot(normal) < 0);
0058       // First solution giving a valid hit wins
0059       if (hit) {
0060         if (distance < -vecgeom::kToleranceStrict<Real_t>) {
0061           Real_t rho     = point.Perp();
0062           auto distanceR = Abs(radius) + point[2] * slope - rho;
0063           Real_t calf    = static_cast<Real_t>(1) / std::sqrt(static_cast<Real_t>(1) + slope * slope);
0064           safety         = distanceR * calf;
0065         }
0066         return true;
0067       }
0068     }
0069     return false;
0070   }
0071
0072   VECGEOM_FORCE_INLINE
0073   VECCORE_ATT_HOST_DEVICE
0074   /// @brief Computes the isotropic safe distance to unplaced surfaces
0075   /// @param point Point in local surface coordinates
0076   /// @param left_side Flag specifying if the surface is intersected from the left-side that defines the normal
0077   /// @param distance Computed isotropic safety
0078   /// @param onsurf Projection of the point on surface
0079   /// @return Validity of the calculation
0080   bool Safety(Vector3D<Real_t> const &point, bool left_side, Real_t &distance, Vector3D<Real_t> &onsurf, Real_t radius,
0081               Real_t slope) const
0082   {
0083     Real_t coneR      = Abs(radius) + point[2] * slope;
0084     Real_t rho        = point.Perp();
0085     bool flip_exiting = left_side ^ (radius < 0);
0086     auto distanceR    = flip_exiting ? coneR - rho : rho - coneR;
0087     Real_t calf       = static_cast<Real_t>(1) / std::sqrt(static_cast<Real_t>(1) + slope * slope);
0088     distance          = distanceR * calf;
0089     // We only use for the ZPhi frame safety the z of the point propagated on the cone surface
0090     onsurf = point;
0091     onsurf[2] += distance * slope * calf;
0092     return distance > -vecgeom::kToleranceDist<Real_t>;
0093   }
0094 };
0095
0096 } // namespace vgbrep
0097
0098 #endif