surfaces/surf/EllipticalImpl.h

0001 #ifndef VECGEOM_ELLIPTICAL_IMPL_H
0002 #define VECGEOM_ELLIPTICAL_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::kElliptical, Real_t> {
0011   EllipData<Real_t> const *fEllipData{nullptr};
0012
0013   VECGEOM_FORCE_INLINE
0014   VECCORE_ATT_HOST_DEVICE
0015   SurfaceHelper(EllipData<Real_t> const &ellipdata) { fEllipData = &ellipdata; }
0016
0017   VECGEOM_FORCE_INLINE
0018   VECCORE_ATT_HOST_DEVICE
0019   /// @brief Inside half-space function
0020   /// @param point Point in local surface coordinates
0021   /// @param tol tolerance for determining if a point is inside or not
0022   /// @return True if the point is behind the normal within kTolerance (surface is included)
0023   bool Inside(Vector3D<Real_t> const &point, bool flip, Real_t tol = vecgeom::kToleranceStrict<Real_t>)
0024   {
0025     int flipsign = !flip ? 1 : -1;
0026     return (point.x() * point.x() / fEllipData->fDDx + point.y() * point.y() / fEllipData->fDDy <
0027             static_cast<Real_t>(1) + 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)
0040   {
0041
0042     distance = vecgeom::InfinityLength<Real_t>();
0043
0044     Vector3D<Real_t> pcur(point);
0045
0046     // MISSING: optimsation to move point closer
0047     // Real_t offset(0.);
0048     // // Move point closer, if required
0049     // Real_v Rfar2(1024. * ellipticaltube.fRsph * ellipticaltube.fRsph); // 1024 = 32 * 32
0050     // vecCore__MaskedAssignFunc(pcur, ((pcur.Mag2() > Rfar2) && (direction.Dot(point) < Real_v(0.))),
0051     //                           pcur + (offset = pcur.Mag() - Real_v(2.) * ellipticaltube.fRsph) * direction);
0052
0053     // Scale elliptical tube to cylinder
0054     Real_t px = pcur.x() * fEllipData->fSx;
0055     Real_t py = pcur.y() * fEllipData->fSy;
0056     Real_t vx = dir.x() * fEllipData->fSx;
0057     Real_t vy = dir.y() * fEllipData->fSy;
0058
0059     // Find intersection with lateral surface, solve equation: A t^2 + 2B t + C = 0
0060     Real_t rr = px * px + py * py;
0061     Real_t A  = vx * vx + vy * vy;
0062     Real_t B  = px * vx + py * vy;
0063     Real_t C  = rr - fEllipData->R * fEllipData->R;
0064
0065     QuadraticCoef<Real_t> coef;
0066     coef.phalf = B / vecgeom::NonZero(A);
0067     coef.q     = C / vecgeom::NonZero(A);
0068
0069     Real_t roots[2];
0070     int numroots;
0071     QuadraticSolver(coef, roots, numroots);
0072     two_solutions = (numroots == 2 && roots[0] > -vecgeom::kToleranceStrict<Real_t> &&
0073                      roots[1] > -vecgeom::kToleranceStrict<Real_t>);
0074     for (auto i = 0; i < numroots; ++i) {
0075       distance                = roots[i];
0076       Vector3D<Real_t> onsurf = point + distance * dir;
0077       Vector3D<Real_t> normal(onsurf[0] * fEllipData->fDDy, onsurf[1] * fEllipData->fDDx, 0);
0078       bool hit = left_side ^ (dir.Dot(normal) < 0);
0079       // First solution giving a valid hit wins
0080       if (hit) {
0081         if (distance < -vecgeom::kToleranceStrict<Real_t>) {
0082           Real_t x   = point.x() * fEllipData->fSx;
0083           Real_t y   = point.y() * fEllipData->fSy;
0084           Real_t rho = vecCore::math::Sqrt(x * x + y * y);
0085           safety     = fEllipData->R - rho;
0086         }
0087         return true;
0088       }
0089     }
0090     return false;
0091   }
0092
0093   VECGEOM_FORCE_INLINE
0094   VECCORE_ATT_HOST_DEVICE
0095   /// @brief Computes the isotropic safe distance to unplaced surfaces
0096   /// @param point Point in local surface coordinates
0097   /// @param left_side Flag specifying if the surface is intersected from the left-side that defines the normal
0098   /// @param distance Computed isotropic safety
0099   /// @param compute_onsurf Instructs to compute the projection of the point on surface
0100   /// @param onsurf Projection of the point on surface
0101   /// @return Validity of the calculation
0102   bool Safety(Vector3D<Real_t> const &point, bool left_side, Real_t &distance, Vector3D<Real_t> &onsurf) const
0103   {
0104     Real_t x   = point.x() * fEllipData->fSx;
0105     Real_t y   = point.y() * fEllipData->fSy;
0106     Real_t rho = vecCore::math::Sqrt(x * x + y * y);
0107     distance   = left_side ? fEllipData->R - rho : rho - fEllipData->R;
0108     onsurf     = point; // we only need the z of the projected point
0109     return distance > -vecgeom::kToleranceDist<Real_t>;
0110   }
0111 };
0112
0113 } // namespace vgbrep
0114
0115 #endif