surfaces/mask/RingMask.h

0001 #ifndef VECGEOM_SURFACE_RINGMASK_H
0002 #define VECGEOM_SURFACE_RINGMASK_H
0003
0004 #include <VecGeom/surfaces/base/CommonTypes.h>
0005
0006 namespace vgbrep {
0007
0008 /// @brief Ring masks on plane surfaces.
0009 /// @tparam Real_t Precision used
0010 template <typename Real_t>
0011 struct RingMask {
0012   using value_type = Real_t;
0013   Range<Real_t> rangeR;        ///< Radius limits in the form of [Rmin, Rmax].
0014   bool isFullCirc;             ///< Does the phi cut exist here?
0015   AngleVector<Real_t> vecSPhi; ///< Cartesian coordinates of vectors that represents the start of the phi-cut.
0016   AngleVector<Real_t> vecEPhi; ///< Cartesian coordinates of vectors that represents the end of the phi-cut.
0017
0018   RingMask() = default;
0019
0020   template <typename Real_i>
0021   RingMask(Real_i rmin, Real_i rmax, bool isFullCircle, Real_i sphi = Real_i{0}, Real_i ephi = Real_i{0})
0022       : rangeR(static_cast<Real_t>(rmin), static_cast<Real_t>(rmax)), isFullCirc(isFullCircle)
0023   {
0024     // If there is no Phi cut, no need to set phi vectors.
0025     if (isFullCirc) return;
0026     vecSPhi.Set(static_cast<Real_t>(vecgeom::Cos(sphi)), static_cast<Real_t>(vecgeom::Sin(sphi)));
0027     vecEPhi.Set(static_cast<Real_t>(vecgeom::Cos(ephi)), static_cast<Real_t>(vecgeom::Sin(ephi)));
0028   }
0029
0030   template <typename Real_i>
0031   RingMask(const RingMask<Real_i> &other)
0032       : rangeR(static_cast<Real_t>(other.rangeR[0]), static_cast<Real_t>(other.rangeR[1])), isFullCirc(other.isFullCirc)
0033   {
0034     // If there is no Phi cut, no need to set phi vectors.
0035     if (isFullCirc) return;
0036     vecSPhi.Set(static_cast<Real_t>(other.vecSPhi[0]), static_cast<Real_t>(other.vecSPhi[1]));
0037     vecEPhi.Set(static_cast<Real_t>(other.vecEPhi[0]), static_cast<Real_t>(other.vecEPhi[1]));
0038   }
0039
0040   /// @brief Fills extents in X and Y for the ring mask
0041   /// @param xmin Minimum of the extent in X
0042   /// @param xmax Maximum of the extent in X
0043   /// @param ymin Minimum of the extent in Y
0044   /// @param ymax Maximum of the extent in Y
0045   void GetXYextent(Real_t &xmin, Real_t &xmax, Real_t &ymin, Real_t &ymax) const
0046   {
0047     auto const &rmax = rangeR[1];
0048     xmin             = -rmax;
0049     xmax             = rmax;
0050     ymin             = -rmax;
0051     ymax             = rmax;
0052     // The axis vector has to be between Rmin and Rmax and cannot be unit vector anymore
0053     Vector2D<Real_t> axis{1, 0};
0054     if (!isFullCirc) {
0055       // Projections of points that delimit vertices of the phi-cut ring
0056       Real_t x1, x2, x3, x4, y1, y2, y3, y4;
0057
0058       auto const &rmin = rangeR[0];
0059
0060       x1 = rmax * axis.Dot(vecSPhi); //< (sphi, Rmax)_x
0061       x2 = rmax * axis.Dot(vecEPhi); //< (ephi, Rmax)_x
0062       x3 = rmin * axis.Dot(vecSPhi); //< (sphi, Rmin)_x
0063       x4 = rmin * axis.Dot(vecEPhi); //< (ephi, Rmin)_x
0064       axis.Set(0, 1);
0065       y1 = rmax * axis.Dot(vecSPhi); //< (sphi, Rmax)_x
0066       y2 = rmax * axis.Dot(vecEPhi); //< (ephi, Rmax)_x
0067       y3 = rmin * axis.Dot(vecSPhi); //< (sphi, Rmin)_x
0068       y4 = rmin * axis.Dot(vecEPhi); //< (ephi, Rmin)_x
0069
0070       xmax = vecgeom::Max(vecgeom::Max(x1, x2), vecgeom::Max(x3, x4));
0071       ymax = vecgeom::Max(vecgeom::Max(y1, y2), vecgeom::Max(y3, y4));
0072       xmin = vecgeom::Min(vecgeom::Min(x1, x2), vecgeom::Min(x3, x4));
0073       ymin = vecgeom::Min(vecgeom::Min(y1, y2), vecgeom::Min(y3, y4));
0074       // If the axes lie within the circle
0075       if (InsidePhi(1, 0)) xmax = rmax;
0076       if (InsidePhi(0, 1)) ymax = rmax;
0077       if (InsidePhi(-1, 0)) xmin = -rmax;
0078       if (InsidePhi(0, -1)) ymin = -rmax;
0079     }
0080   }
0081
0082   /// @brief Fills the 3D extent of the ring
0083   /// @param aMin Bottom extent corner
0084   /// @param aMax Top extent corner
0085   void Extent3D(Vector3D<Real_t> &aMin, Vector3D<Real_t> &aMax) const
0086   {
0087     Real_t xmin, xmax, ymin, ymax;
0088     GetXYextent(xmin, xmax, ymin, ymax);
0089     aMin.Set(xmin, ymin, Real_t(0));
0090     aMax.Set(xmax, ymax, Real_t(0));
0091   }
0092
0093   /// @brief Returns the window extent of the ring
0094   /// @param window Extent window to be filled
0095   void GetExtent(WindowMask<Real_t> &window) const
0096   {
0097     Real_t xmin, xmax, ymin, ymax;
0098     GetXYextent(xmin, xmax, ymin, ymax);
0099     window.rangeU.Set(xmin, xmax);
0100     window.rangeV.Set(ymin, ymax);
0101   }
0102
0103   VECCORE_ATT_HOST_DEVICE
0104   VECGEOM_FORCE_INLINE
0105   bool HasRmin() const { return rangeR[0] > vecgeom::kToleranceDist<Real_t>; }
0106
0107   VECCORE_ATT_HOST_DEVICE
0108   VECGEOM_FORCE_INLINE
0109   bool IsConvex() const { return !HasRmin() && vecSPhi.CrossZ(vecEPhi) > -vecgeom::kToleranceDist<Real_t>; }
0110
0111   /// @brief Check if local point is in the radius range
0112   /// @param local Point in local coordinates
0113   /// @return Point inside range
0114   VECCORE_ATT_HOST_DEVICE
0115   VECGEOM_FORCE_INLINE
0116   bool InsideR(Vector3D<Real_t> const &local) const
0117   {
0118     Real_t rsq = local[0] * local[0] + local[1] * local[1];
0119     // The point must be inside the ring:
0120     if ((rsq < rangeR[0] * rangeR[0] - 2 * vecgeom::kRelTolerance<Real_t>(rangeR[0])) ||
0121         (rsq > rangeR[1] * rangeR[1] + 2 * vecgeom::kRelTolerance<Real_t>(rangeR[1])))
0122       return false;
0123     return true;
0124   }
0125
0126   /// @brief Check if local point is in the phi range
0127   /// @param x Local x coordinate
0128   /// @param y Local y coordinate
0129   /// @return Point inside phi
0130   VECCORE_ATT_HOST_DEVICE
0131   VECGEOM_FORCE_INLINE
0132   bool InsidePhi(Real_t const &x, Real_t const &y) const
0133   {
0134     if (isFullCirc) return true;
0135     AngleVector<Real_t> localAngle{x, y};
0136     auto convex = vecSPhi.CrossZ(vecEPhi) > Real_t(0);
0137     auto in1    = vecSPhi.CrossZ(localAngle) > -vecgeom::kToleranceStrict<Real_t>;
0138     auto in2    = localAngle.CrossZ(vecEPhi) > -vecgeom::kToleranceStrict<Real_t>;
0139     return convex ? in1 && in2 : in1 || in2;
0140   }
0141
0142   /// @brief Checks if the point is within mask.
0143   /// @details For the phi part, use the cross-products of the point vector with the unit
0144   ///  phi vectors, representing the signed safeties with respect to these vectors. The
0145   ///  point must also satisfy the radius limits.
0146   /// @param local Local coordinates of the point
0147   /// @return Inside the mask or not
0148   VECCORE_ATT_HOST_DEVICE
0149   VECGEOM_FORCE_INLINE
0150   bool Inside(Vector3D<Real_t> const &local) const { return InsideR(local) && InsidePhi(local[0], local[1]); }
0151
0152   /// @brief Computes safe distance to the frame combining surface and frame safeties (under-estimate)
0153   /// @details Computes first the maximum signed distance to each edge on a single axis. Two versions
0154   /// @param local Projected point in local coordinates
0155   /// @param safetySurf Safety from non-projected point to the frame support surface.
0156   /// @return Safety distance to the rectangle mask.
0157   VECCORE_ATT_HOST_DEVICE
0158   Real_t Safety(Vector3D<Real_t> const &local, Real_t safetySurf, bool &valid) const
0159   {
0160     valid       = true;
0161     Real_t rho  = local.Perp(); // still a square root, but less registers/branching
0162     Real_t safR = HasRmin() ? vecCore::math::Max(rangeR[0] - rho, rho - rangeR[1]) : rho - rangeR[1];
0163     if (isFullCirc || InsidePhi(local[0], local[1]))
0164 #ifdef SURF_ACCURATE_SAFETY
0165       return vecCore::math::Sqrt(safetySurf * safetySurf + ((safR > 0) ? safR * safR : 0));
0166 #else
0167       return vecCore::math::Max(safetySurf, safR);
0168 #endif
0169     AngleVector<Real_t> localAngle{local[0], local[1]};
0170     auto safphi1    = localAngle.CrossZ(vecSPhi);
0171     auto safphi2    = -localAngle.CrossZ(vecEPhi);
0172     bool largerPi   = vecSPhi.CrossZ(vecEPhi) < 0;
0173     auto safPhi     = largerPi ? vecCore::math::Min(safphi1, safphi2) : vecCore::math::Max(safphi1, safphi2);
0174     auto safInPlane = vecCore::math::Max(safR, safPhi);
0175 #ifdef SURF_ACCURATE_SAFETY
0176     auto safety = safetySurf * safetySurf + ((safInPlane > 0) ? safInPlane * safInPlane : 0);
0177     return vecCore::math::Sqrt(safety);
0178 #else
0179     return vecCore::math::Max(safetySurf, safInPlane);
0180 #endif
0181   }
0182
0183   /// @brief Safe distance from a point assumed inside the ring
0184   /// @details Used on host only for frame checks
0185   /// @param local Projected point in local coordinates
0186   /// @return Safe distance
0187   Real_t SafetyInside(Vector3D<Real_t> const &local) const
0188   {
0189     Real_t rho  = local.Perp();
0190     Real_t safR = HasRmin() ? vecCore::math::Min(rangeR[1] - rho, rho - rangeR[0]) : rangeR[1] - rho;
0191     safR        = vecCore::math::Max(safR, Real_t(0));
0192     if (isFullCirc) return safR;
0193     AngleVector<Real_t> localAngle{local[0], local[1]};
0194     auto saf1     = (vecSPhi.Dot(localAngle) > Real_t(0)) ? vecCore::math::Max(vecSPhi.CrossZ(localAngle), Real_t(0))
0195                                                           : vecgeom::InfinityLength<Real_t>();
0196     auto saf2     = (vecEPhi.Dot(localAngle) > Real_t(0)) ? vecCore::math::Max(localAngle.CrossZ(vecEPhi), Real_t(0))
0197                                                           : vecgeom::InfinityLength<Real_t>();
0198     Real_t safPhi = vecCore::math::Min(saf1, saf2);
0199     return vecCore::math::Min(safR, safPhi);
0200   }
0201
0202   /// @brief Safe distance from a point assumed outside the ring
0203   /// @details Used on host only for frame checks
0204   /// @param local Projected point in local coordinates
0205   /// @return Safe distance
0206   Real_t SafetyOutside(Vector3D<Real_t> const &local) const
0207   {
0208     Real_t rho  = local.Perp();
0209     Real_t safR = HasRmin() ? vecCore::math::Max(rho - rangeR[1], rangeR[0] - rho) : rho - rangeR[1];
0210     safR        = vecCore::math::Max(safR, Real_t(0));
0211     if (isFullCirc) return safR;
0212     AngleVector<Real_t> localAngle{local[0], local[1]};
0213     auto saf1     = (vecSPhi.Dot(localAngle) > Real_t(0)) ? vecCore::math::Max(localAngle.CrossZ(vecSPhi), Real_t(0))
0214                                                           : vecgeom::InfinityLength<Real_t>();
0215     auto saf2     = (vecEPhi.Dot(localAngle) > Real_t(0)) ? vecCore::math::Max(vecEPhi.CrossZ(localAngle), Real_t(0))
0216                                                           : vecgeom::InfinityLength<Real_t>();
0217     Real_t safPhi = vecCore::math::Min(saf1, saf2);
0218     if (safPhi == vecgeom::InfinityLength<Real_t>()) safPhi = Real_t(0);
0219     return vecCore::math::Max(safR, safPhi);
0220   }
0221 };
0222
0223 } // namespace vgbrep
0224
0225 #endif