surfaces/surf/CylindricalImpl.h

0001 #ifndef VECGEOM_CYLINDRICAL_IMPL_H
0002 #define VECGEOM_CYLINDRICAL_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::kCylindrical, 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.Perp();
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   /// @param two_solutions whether there are two possible solutions
0038   /// @return Validity of the intersection
0039   bool Intersect(Vector3D<Real_t> const &point, Vector3D<Real_t> const &dir, bool left_side, Real_t &distance,
0040                  bool &two_solutions, Real_t &safety, Real_t radius)
0041   {
0042     QuadraticCoef<Real_t> coef;
0043     Real_t roots[2];
0044     int numroots      = 0;
0045     bool flip_exiting = left_side ^ (radius < 0);
0046     CylinderEq<Real_t>(point, dir, Abs(radius), coef);
0047     QuadraticSolver(coef, roots, numroots);
0048     two_solutions = (numroots == 2 && roots[0] > -vecgeom::kToleranceStrict<Real_t> &&
0049                      roots[1] > -vecgeom::kToleranceStrict<Real_t>);
0050     for (auto i = 0; i < numroots; ++i) {
0051       distance                = roots[i];
0052       Vector3D<Real_t> onsurf = point + distance * dir;
0053       Vector3D<Real_t> normal(onsurf[0], onsurf[1], 0);
0054       bool hit = flip_exiting ^ (dir.Dot(normal) < 0);
0055       // First solution giving a valid hit wins
0056       if (hit) {
0057         if (distance < -vecgeom::kToleranceStrict<Real_t> && distance < -Abs(radius)) {
0058           Real_t rho = point.Perp();
0059           safety     = Abs(radius) - rho;
0060         }
0061         return true;
0062       }
0063     }
0064     return false;
0065   }
0066
0067   VECGEOM_FORCE_INLINE
0068   VECCORE_ATT_HOST_DEVICE
0069   /// @brief Computes the isotropic safe distance to unplaced surfaces
0070   /// @param point Point in local surface coordinates
0071   /// @param left_side Flag specifying if the surface is intersected from the left-side that defines the normal
0072   /// @param distance Computed isotropic safety
0073   /// @param compute_onsurf Instructs to compute the projection of the point on surface
0074   /// @param onsurf Projection of the point on surface
0075   /// @return Validity of the calculation
0076   bool Safety(Vector3D<Real_t> const &point, bool left_side, Real_t &distance, Vector3D<Real_t> &onsurf,
0077               Real_t radius) const
0078   {
0079     Real_t rho        = point.Perp();
0080     bool flip_exiting = left_side ^ (radius < 0);
0081     distance          = flip_exiting ? Abs(radius) - rho : rho - Abs(radius);
0082     onsurf            = point; // we only need the z of the projected point
0083     return distance > -vecgeom::kToleranceDist<Real_t>;
0084   }
0085 };
0086
0087 } // namespace vgbrep
0088
0089 #endif