surfaces/surf/SphericalImpl.h

0001 #ifndef VECGEOM_SPHERICAL_IMPL_H
0002 #define VECGEOM_SPHERICAL_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::kSpherical, 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 tol = vecgeom::kToleranceStrict<Real_t>)
0023   {
0024     int flipsign      = (radius < 0) ? -1 : 1;
0025     int bool_flipsign = !flip ? 1 : -1;
0026     Real_t rho        = point.Mag();
0027     return flipsign * (rho - Abs(radius)) < bool_flipsign * tol;
0028   }
0029
0030   VECGEOM_FORCE_INLINE
0031   VECCORE_ATT_HOST_DEVICE
0032   /// @brief Find signed distance to next intersection from local point.
0033   /// @param point Point in local surface coordinates
0034   /// @param dir Direction in the local surface coordinates
0035   /// @param left_side Flag specifying if the surface is intersected from the left-side that defines the normal
0036   /// @param distance Computed distance to surface
0037   /// @return Validity of the intersection
0038   bool Intersect(Vector3D<Real_t> const &point, Vector3D<Real_t> const &dir, bool left_side, Real_t &distance,
0039                  bool &two_solutions, Real_t &safety, Real_t radius)
0040   {
0041     QuadraticCoef<Real_t> coef;
0042     Real_t roots[2];
0043     int numroots      = 0;
0044     bool flip_exiting = left_side ^ (radius < 0);
0045     SphereEq<Real_t>(point, dir, radius, coef);
0046     QuadraticSolver(coef, roots, numroots);
0047     two_solutions = (numroots == 2 && roots[0] > -vecgeom::kToleranceStrict<Real_t> &&
0048                      roots[1] > -vecgeom::kToleranceStrict<Real_t>);
0049     for (auto i = 0; i < numroots; ++i) {
0050       distance                = roots[i];
0051       Vector3D<Real_t> onsurf = point + distance * dir;
0052       Vector3D<Real_t> normal(onsurf[0], onsurf[1], onsurf[2]);
0053       bool hit = flip_exiting ^ (dir.Dot(normal) < 0);
0054       // First solution giving a valid hit wins
0055       if (hit) {
0056         if (distance < -vecgeom::kToleranceStrict<Real_t> && distance < -radius) {
0057           Real_t rho = point.Mag();
0058           safety     = Abs(radius) - rho;
0059         }
0060         return true;
0061       }
0062     }
0063     return false;
0064   }
0065
0066   VECGEOM_FORCE_INLINE
0067   VECCORE_ATT_HOST_DEVICE
0068   /// @brief Computes the isotropic safe distance to unplaced surfaces
0069   /// @param point Point in local surface coordinates
0070   /// @param left_side Flag specifying if the surface is intersected from the left-side that defines the normal
0071   /// @param distance Computed isotropic safety
0072   /// @param onsurf Projection of the point on surface
0073   /// @return Validity of the calculation
0074   bool Safety(Vector3D<Real_t> const &point, bool left_side, Real_t &distance, Vector3D<Real_t> &onsurf,
0075               Real_t radius) const
0076   {
0077     Real_t rho = point.Mag();
0078     distance   = left_side ? Abs(radius) - rho : rho - Abs(radius);
0079     // the onsurf computation code is missing below
0080
0081     return distance > -vecgeom::kToleranceDist<Real_t>;
0082   }
0083 };
0084
0085 } // namespace vgbrep
0086
0087 #endif