surfaces/surf/Arb4Impl.h

0001 #ifndef VECGEOM_ARB4_IMPL_H
0002 #define VECGEOM_ARB4_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::kArb4, Real_t> {
0011   Arb4Data<Real_t> const *fArb4Data{nullptr};
0012
0013   VECGEOM_FORCE_INLINE
0014   VECCORE_ATT_HOST_DEVICE
0015   SurfaceHelper(Arb4Data<Real_t> const &arbdata) { fArb4Data = &arbdata; }
0016
0017   /// @brief Fills the 3D extent of the Arb4
0018   /// @param aMin Bottom extent corner
0019   /// @param aMax Top extent corner
0020   void Extent3D(Vector3D<Real_t> &aMin, Vector3D<Real_t> &aMax) const
0021   {
0022     Real_t dz = fArb4Data->halfH;
0023     aMin.Set(fArb4Data->verticesX[0], fArb4Data->verticesY[0], -dz);
0024     aMax = aMin;
0025     for (auto i = 0; i < 4; ++i) {
0026       aMin = vecCore::math::Min(aMin, Vector3D<Real_t>(fArb4Data->verticesX[i], fArb4Data->verticesY[i], -dz));
0027       aMax = vecCore::math::Max(aMax, Vector3D<Real_t>(fArb4Data->verticesX[i], fArb4Data->verticesY[i], -dz));
0028       aMin = vecCore::math::Min(aMin, Vector3D<Real_t>(fArb4Data->verticesX[i], fArb4Data->verticesY[i], dz));
0029       aMax = vecCore::math::Max(aMax, Vector3D<Real_t>(fArb4Data->verticesX[i], fArb4Data->verticesY[i], dz));
0030     }
0031   }
0032
0033   VECGEOM_FORCE_INLINE
0034   VECCORE_ATT_HOST_DEVICE
0035   /// @brief Inside half-space function
0036   /// @param point Point in local surface coordinates
0037   /// @param flip flipping the tolerance for inside (needed for negated booleans)
0038   /// @param tol tolerance for determining if a point is inside or not
0039   /// @return True if the point is behind the normal within kTolerance (surface is included)
0040   bool Inside(Vector3D<Real_t> const &point, bool flip, Real_t tol = vecgeom::kToleranceStrict<Real_t>)
0041   {
0042     Real_t dz     = fArb4Data->halfH;
0043     Real_t dz_inv = fArb4Data->halfH_inv;
0044
0045     // distance along the top-bottom connecting line to reach the z value of the point.
0046     Real_t cf = 0.5 * dz_inv * (dz - point.z());
0047
0048     //  loop over edges connecting points i with i+4
0049     Real_t vertexX[2];
0050     Real_t vertexY[2];
0051
0052     for (int i = 0; i < 2; i++) {
0053       // calculate x-y positions of vertex i at this z-height
0054       vertexX[i] = fArb4Data->verticesX[i + 2] + cf * fArb4Data->connecting_compX[i];
0055       vertexY[i] = fArb4Data->verticesY[i + 2] + cf * fArb4Data->connecting_compY[i];
0056     }
0057
0058     Real_t dx = vertexX[1] - vertexX[0];
0059     Real_t dy = vertexY[1] - vertexY[0];
0060     Real_t px = point.x() - vertexX[0];
0061     Real_t py = point.y() - vertexY[0];
0062
0063     int sign = !flip ? -1 : 1;
0064     return (dx * py - dy * px) > sign * tol;
0065   }
0066
0067   VECGEOM_FORCE_INLINE
0068   VECCORE_ATT_HOST_DEVICE
0069   /// @brief Find signed distance to next intersection from local point.
0070   /// @param point Point in local surface coordinates
0071   /// @param dir Direction in the local surface coordinates
0072   /// @param left_side Flag specifying if the surface is intersected from the left-side that defines the normal
0073   /// @param distance Computed distance to surface
0074   /// @return Validity of the intersection
0075   bool Intersect(Vector3D<Real_t> const &point, Vector3D<Real_t> const &dir, bool left_side, Real_t &distance,
0076                  bool &two_solutions, Real_t &safety)
0077   {
0078
0079     using Vector3D = vecgeom::Vector3D<Real_t>;
0080
0081     Real_t dzp    = point[2] + fArb4Data->halfH;
0082     Real_t dz_inv = fArb4Data->halfH_inv;
0083
0084     // positions in x and y at the connection between top and bottom points as a function of z
0085     Real_t xs1 = fArb4Data->verticesX[1] + fArb4Data->ftx1 * dzp;
0086     Real_t ys1 = fArb4Data->verticesY[1] + fArb4Data->fty1 * dzp;
0087     Real_t xs2 = fArb4Data->verticesX[0] + fArb4Data->ftx2 * dzp;
0088     Real_t ys2 = fArb4Data->verticesY[0] + fArb4Data->fty2 * dzp;
0089     Real_t dxs = xs2 - xs1;
0090     Real_t dys = ys2 - ys1;
0091
0092     // calculate quadratic coefficients
0093     Real_t a = (fArb4Data->fDeltatx * dir[1] - fArb4Data->fDeltaty * dir[0] + fArb4Data->ft1crosst2 * dir[2]) * dir[2];
0094     Real_t b = dxs * dir[1] - dys * dir[0] +
0095                (fArb4Data->fDeltatx * point[1] - fArb4Data->fDeltaty * point[0] + fArb4Data->fty2 * xs1 -
0096                 fArb4Data->fty1 * xs2 + fArb4Data->ftx1 * ys2 - fArb4Data->ftx2 * ys1) *
0097                    dir[2];
0098     Real_t c = dxs * point[1] - dys * point[0] + xs1 * ys2 - xs2 * ys1;
0099
0100     QuadraticCoef<Real_t> coef;
0101     coef.phalf = Real_t(1.) * b / (Real_t(2.) * vecgeom::NonZero(a));
0102     coef.q     = c / vecgeom::NonZero(a);
0103
0104     Real_t roots[2];
0105     int numroots;
0106     QuadraticSolver(coef, roots, numroots);
0107
0108     // validate hits with normal and frame check
0109     for (auto i = 0; i < numroots; ++i) {
0110       if (roots[i] < -vecgeom::kTolerance) continue;
0111       distance        = roots[i];
0112       Vector3D onsurf = point + distance * dir;
0113       if (Abs(onsurf.z()) > fArb4Data->halfH + vecgeom::kTolerance) continue;
0114
0115       // calculate the normal
0116       // when calculating the normal, we also calculate the vertical ratio rz and the horizontal ratios r in x and y.
0117       // They can be immediately used as a framecheck, so that the frame does not need to be checked later again.
0118       Real_t r   = -1;
0119       Real_t rz  = 0.5 * dz_inv * (onsurf[2] + fArb4Data->halfH);
0120       Real_t num = (onsurf.x() - fArb4Data->verticesX[1]) - rz * (fArb4Data->verticesX[3] - fArb4Data->verticesX[1]);
0121       Real_t denom =
0122           (fArb4Data->verticesX[0] - fArb4Data->verticesX[1]) +
0123           rz * (fArb4Data->verticesX[2] - fArb4Data->verticesX[0] - fArb4Data->verticesX[3] + fArb4Data->verticesX[1]);
0124       if (Abs(denom) > 1e-6) {
0125         r = num / vecgeom::NonZero(denom);
0126         if (!((r >= -vecgeom::kTolerance) && (r <= 1. + vecgeom::kTolerance))) continue;
0127       }
0128       num = (onsurf.y() - fArb4Data->verticesY[1]) - rz * (fArb4Data->verticesY[3] - fArb4Data->verticesY[1]);
0129       denom =
0130           (fArb4Data->verticesY[0] - fArb4Data->verticesY[1]) +
0131           rz * (fArb4Data->verticesY[2] - fArb4Data->verticesY[0] - fArb4Data->verticesY[3] + fArb4Data->verticesY[1]);
0132       if (Abs(denom) > 1e-6) r = num / vecgeom::NonZero(denom);
0133
0134       // frame check?
0135       if (!((r >= -vecgeom::kTolerance) && (r <= 1. + vecgeom::kTolerance))) continue;
0136       if (!((rz >= -vecgeom::kTolerance) && (rz <= 1. + vecgeom::kTolerance))) continue;
0137
0138       // calculating normal
0139       Vector3D unorm = fArb4Data->fViCrossHi0 + rz * fArb4Data->fViCrossVj + r * fArb4Data->fHi1CrossHi0;
0140
0141       bool hit = left_side ^ (dir.Dot(unorm) < 0);
0142       // First solution giving a valid hit wins
0143       if (hit) {
0144         if (distance < -vecgeom::kTolerance) safety = Safety(point, left_side, distance, onsurf);
0145         return true;
0146       }
0147     }
0148     return false;
0149   }
0150
0151   VECGEOM_FORCE_INLINE
0152   VECCORE_ATT_HOST_DEVICE
0153   /// @brief Computes the isotropic safe distance to unplaced surfaces
0154   /// @param point Point in local surface coordinates
0155   /// @param left_side Flag specifying if the surface is intersected from the left-side that defines the normal
0156   /// @param distance Computed isotropic safety
0157   /// @param onsurf Projection of the point on surface
0158   /// @return Validity of the calculation
0159   bool Safety(Vector3D<Real_t> const &point, bool left_side, Real_t &distance, Vector3D<Real_t> &onsurf) const
0160   {
0161
0162     Real_t dzp = fArb4Data->halfH + point[2];
0163
0164     // safety in z
0165     distance = (fabs(point[2]) > fArb4Data->halfH + vecgeom::kTolerance) ? fabs(point[2]) - fArb4Data->halfH : 0.;
0166
0167 #ifdef SURF_ACCURATE_SAFETY
0168     Real_t v_safety         = vecgeom::InfinityLength<Real_t>();
0169     Real_t tmp              = 0.;
0170     Vector3D<Real_t> vplane = {fArb4Data->verticesX[0], fArb4Data->verticesY[0], -fArb4Data->halfH};
0171     tmp                     = (point - vplane).Dot(fArb4Data->normal0);
0172     if (tmp > 0) v_safety = tmp;
0173
0174     vplane = {fArb4Data->verticesX[1], fArb4Data->verticesY[1], -fArb4Data->halfH};
0175     tmp    = (point - vplane).Dot(fArb4Data->normal1);
0176     if (tmp > 0) v_safety = Min(v_safety, tmp);
0177
0178     vplane = {fArb4Data->verticesX[2], fArb4Data->verticesY[2], fArb4Data->halfH};
0179     tmp    = (point - vplane).Dot(fArb4Data->normal2);
0180     if (tmp > 0) v_safety = Min(v_safety, tmp);
0181
0182     vplane = {fArb4Data->verticesX[3], fArb4Data->verticesY[3], -fArb4Data->halfH};
0183     tmp    = (point - vplane).Dot(fArb4Data->normal3);
0184     if (tmp > 0) v_safety = Min(v_safety, tmp);
0185
0186     if (v_safety < vecgeom::InfinityLength<Real_t>()) {
0187       distance = Max(distance, v_safety);
0188     }
0189
0190 #endif
0191
0192     Real_t xs1 = fArb4Data->verticesX[1] + fArb4Data->ftx1 * dzp;
0193     Real_t ys1 = fArb4Data->verticesY[1] + fArb4Data->fty1 * dzp;
0194     Real_t xs2 = fArb4Data->verticesX[0] + fArb4Data->ftx2 * dzp;
0195     Real_t ys2 = fArb4Data->verticesY[0] + fArb4Data->fty2 * dzp;
0196     Real_t dxs = xs2 - xs1;
0197     Real_t dys = ys2 - ys1;
0198
0199     Real_t umin   = 0.;
0200     Real_t safety = vecgeom::InfinityLength<Real_t>();
0201
0202     Real_t dx1 = 0.;
0203     Real_t dx2 = 0.;
0204     Real_t dy1 = 0.;
0205     Real_t dy2 = 0.;
0206
0207     Real_t dpx = point[0] - xs1;
0208     Real_t dpy = point[1] - ys1;
0209     Real_t lsq = dxs * dxs + dys * dys;
0210     Real_t u   = (dpx * dxs + dpy * dys) / vecgeom::NonZero(lsq);
0211     if (u > 1 + vecgeom::kTolerance) {
0212       dpx = point[0] - xs2;
0213       dpy = point[1] - ys2;
0214     }
0215     if (u >= -vecgeom::kTolerance && u <= 1 + vecgeom::kTolerance) {
0216       dpx = dpx - u * dxs;
0217       dpy = dpy - u * dys;
0218     }
0219     Real_t ssq = dpx * dpx + dpy * dpy; // safety squared
0220     if (ssq < safety) {
0221       dx1    = fArb4Data->verticesX[0] - fArb4Data->verticesX[1];
0222       dx2    = fArb4Data->verticesX[2] - fArb4Data->verticesX[3];
0223       dy1    = fArb4Data->verticesY[0] - fArb4Data->verticesY[1];
0224       dy2    = fArb4Data->verticesY[2] - fArb4Data->verticesY[3];
0225       umin   = u;
0226       safety = ssq;
0227     }
0228     if (umin < -vecgeom::kTolerance || umin > 1 + vecgeom::kTolerance) umin = 0.;
0229     dxs = dx1 + umin * (dx2 - dx1);
0230     dys = dy1 + umin * (dy2 - dy1);
0231     // Denominator below always positive as fDz>0
0232     safety *= Real_t(1.) - Real_t(4.) * fArb4Data->halfH * fArb4Data->halfH /
0233                                (dxs * dxs + dys * dys + Real_t(4.) * fArb4Data->halfH * fArb4Data->halfH);
0234     distance = Max(distance, Sqrt(safety));
0235
0236     return distance > -vecgeom::kToleranceDist<Real_t>;
0237   }
0238 };
0239
0240 } // namespace vgbrep
0241
0242 #endif