surfaces/mask/ZPhiMask.h

0001 #ifndef VECGEOM_SURFACE_ZPHIMASK_H
0002 #define VECGEOM_SURFACE_ZPHIMASK_H
0003
0004 #include <VecGeom/management/Logger.h>
0005 #include <VecGeom/surfaces/base/CommonTypes.h>
0006
0007 namespace vgbrep {
0008
0009 /// @brief Mask in cylindrical coordinates, used on cylindrical surfaces.
0010 /// @tparam Real_t Precision type
0011 template <typename Real_t>
0012 struct ZPhiMask {
0013   //// ZPhi format:
0014   //// rangeZ             -> extent along z axis
0015   //// isFullCirc         -> Does the phi cut exist here?
0016   //// vecSPhi, vecEPhi   -> Cartesian coordinates of vectors that delimit phi-cut
0017
0018   using value_type = Real_t;
0019   Range<Real_t> rangeZ;        ///< Limits on the z-axis.
0020   Real_t r0;                   ///< Radius at z = 0
0021   Real_t invcalf;              ///< Inverse of the cosine of the surface angle with respect to z axis
0022   AngleVector<Real_t> vecSPhi; ///< Cartesian coordinates of vectors that represents the start of the phi-cut.
0023   AngleVector<Real_t> vecEPhi; ///< Cartesian coordinates of vectors that represents the end of the phi-cut.
0024   bool isFullCirc;             ///< Does the phi cut exist here?
0025
0026   ZPhiMask() = default;
0027   template <typename Real_i>
0028   ZPhiMask(Real_i zmin, Real_i zmax, bool isFullCircle, Real_i rbottom, Real_i rtop, Real_i sphi = Real_i{0},
0029            Real_i ephi = Real_i{0})
0030       : rangeZ(static_cast<Real_t>(zmin), static_cast<Real_t>(zmax)), isFullCirc(isFullCircle)
0031   {
0032     r0       = static_cast<Real_t>(0.5 * (rbottom + rtop));
0033     Real_t t = static_cast<Real_t>((rtop - rbottom) / (zmax - zmin));
0034     invcalf  = (t == Real_t(0)) ? Real_t(1) : vecCore::math::Sqrt(Real_t(1) + t * t);
0035     // If there is no Phi cut, we needn't worry about phi vectors.
0036     if (isFullCirc) return;
0037     vecSPhi.Set(static_cast<Real_t>(vecgeom::Cos(sphi)), static_cast<Real_t>(vecgeom::Sin(sphi)));
0038     vecEPhi.Set(static_cast<Real_t>(vecgeom::Cos(ephi)), static_cast<Real_t>(vecgeom::Sin(ephi)));
0039   }
0040   template <typename Real_i>
0041   ZPhiMask(const ZPhiMask<Real_i> &other)
0042       : r0(static_cast<Real_t>(other.r0)), invcalf(static_cast<Real_t>(other.invcalf)), isFullCirc(other.isFullCirc)
0043   {
0044     rangeZ  = Range<Real_t>(other.rangeZ);
0045     vecSPhi = AngleVector<Real_t>(other.vecSPhi);
0046     vecEPhi = AngleVector<Real_t>(other.vecEPhi);
0047   }
0048
0049   /// @brief Fills extents in X and Y for the ZPhi mask
0050   /// @param xmin Minimum of the extent in X
0051   /// @param xmax Maximum of the extent in X
0052   /// @param ymin Minimum of the extent in Y
0053   /// @param ymax Maximum of the extent in Y
0054   void GetXYextent(Real_t &xmin, Real_t &xmax, Real_t &ymin, Real_t &ymax) const
0055   {
0056     auto const rmax = vecCore::math::Max(Radius(rangeZ[0]), Radius(rangeZ[1]));
0057     xmin            = -rmax;
0058     xmax            = rmax;
0059     ymin            = -rmax;
0060     ymax            = rmax;
0061     // The axis vector has to be between Rmin and Rmax and cannot be unit vector anymore
0062     Vector2D<Real_t> axis{1, 0};
0063     if (!isFullCirc) {
0064       // Projections of points that delimit vertices of the phi-cut ring
0065       Real_t x1, x2, x3, x4, y1, y2, y3, y4;
0066
0067       auto const rmin = vecCore::math::Min(Radius(rangeZ[0]), Radius(rangeZ[1]));
0068
0069       x1 = rmax * axis.Dot(vecSPhi); //< (sphi, Rmax)_x
0070       x2 = rmax * axis.Dot(vecEPhi); //< (ephi, Rmax)_x
0071       x3 = rmin * axis.Dot(vecSPhi); //< (sphi, Rmin)_x
0072       x4 = rmin * axis.Dot(vecEPhi); //< (ephi, Rmin)_x
0073       axis.Set(0, 1);
0074       y1 = rmax * axis.Dot(vecSPhi); //< (sphi, Rmax)_x
0075       y2 = rmax * axis.Dot(vecEPhi); //< (ephi, Rmax)_x
0076       y3 = rmin * axis.Dot(vecSPhi); //< (sphi, Rmin)_x
0077       y4 = rmin * axis.Dot(vecEPhi); //< (ephi, Rmin)_x
0078
0079       xmax = vecgeom::Max(vecgeom::Max(x1, x2), vecgeom::Max(x3, x4));
0080       ymax = vecgeom::Max(vecgeom::Max(y1, y2), vecgeom::Max(y3, y4));
0081       xmin = vecgeom::Min(vecgeom::Min(x1, x2), vecgeom::Min(x3, x4));
0082       ymin = vecgeom::Min(vecgeom::Min(y1, y2), vecgeom::Min(y3, y4));
0083       // If the axes lie within the circle
0084       if (InsidePhi(1, 0)) xmax = rmax;
0085       if (InsidePhi(0, 1)) ymax = rmax;
0086       if (InsidePhi(-1, 0)) xmin = -rmax;
0087       if (InsidePhi(0, -1)) ymin = -rmax;
0088     }
0089   }
0090
0091   /// @brief Fills the 3D extent of the mask
0092   /// @param aMin Bottom extent corner
0093   /// @param aMax Top extent corner
0094   void Extent3D(Vector3D<Real_t> &aMin, Vector3D<Real_t> &aMax) const
0095   {
0096     Real_t xmin, xmax, ymin, ymax;
0097     GetXYextent(xmin, xmax, ymin, ymax);
0098     aMin.Set(xmin, ymin, rangeZ[0]);
0099     aMax.Set(xmax, ymax, rangeZ[1]);
0100   }
0101
0102   /// @brief Get radius of the mask dependent on z
0103   /// @param z Z value
0104   /// @return Radius of the mask
0105   Real_t Radius(Real_t z) const
0106   {
0107     Real_t t = vecCore::math::Sqrt((invcalf - Real_t(1)) * (invcalf + Real_t(1)));
0108     return r0 + t * z;
0109   }
0110
0111   /// @brief Check if local point is in the phi range
0112   /// @param x Local x coordinate
0113   /// @param y Local y coordinate
0114   /// @return Point inside phi
0115   VECCORE_ATT_HOST_DEVICE
0116   VECGEOM_FORCE_INLINE
0117   bool InsidePhi(Real_t const &x, Real_t const &y) const
0118   {
0119     if (isFullCirc) return true;
0120     AngleVector<Real_t> localAngle{x, y};
0121     auto convex = vecSPhi.CrossZ(vecEPhi) > Real_t(0);
0122     auto in1    = vecSPhi.CrossZ(localAngle) > -vecgeom::kToleranceStrict<Real_t>;
0123     auto in2    = localAngle.CrossZ(vecEPhi) > -vecgeom::kToleranceStrict<Real_t>;
0124     return convex ? in1 && in2 : in1 || in2;
0125   }
0126
0127   /// @brief Checks if the point is within mask.
0128   /// @details For the phi part, use the cross-products of the point vector with the unit
0129   ///  phi vectors, representing the signed safeties with respect to these vectors. The
0130   ///  point must also satisfy the Z limits.
0131   /// @param local Local coordinates of the point
0132   /// @return Inside the mask or not
0133   VECCORE_ATT_HOST_DEVICE
0134   bool Inside(Vector3D<Real_t> const &local) const
0135   {
0136     // The point must be inside z-span:
0137     if (local[2] < vecgeom::MakeMinusTolerant<true, Real_t>(rangeZ[0]) ||
0138         local[2] > vecgeom::MakePlusTolerant<true, Real_t>(rangeZ[1]))
0139       return false;
0140     return InsidePhi(local[0], local[1]);
0141   }
0142
0143   /// @brief Transform a ZPhi mask from a local reference defined by trans to the parent reference
0144   /// @param trans Transformation of the ZPhi mask with respect to the parent reference
0145   /// @return Transformed mask
0146   ZPhiMask<Real_t> InverseTransform(vecgeom::Transformation2DMP<Real_t> const &trans) const
0147   {
0148     ZPhiMask<Real_t> frame(*this);
0149     // Convert rangeZ
0150     Vector3D<Real_t> local;
0151     local           = trans.InverseTransform(Vector3D<Real_t>{0, 0, rangeZ[0]});
0152     frame.rangeZ[0] = local[2];
0153     local           = trans.InverseTransform(Vector3D<Real_t>{0, 0, rangeZ[1]});
0154     frame.rangeZ[1] = local[2];
0155
0156     // Convert radius at 0
0157     frame.r0 = Radius(-trans.Translation()[2]);
0158
0159     if (!isFullCirc) {
0160       // Convert phi range
0161       local = trans.InverseTransformDirection(Vector3D<Real_t>{vecSPhi[0], vecSPhi[1], 0});
0162       frame.vecSPhi.Set(local[0], local[1]);
0163       local = trans.InverseTransformDirection(Vector3D<Real_t>{vecEPhi[0], vecEPhi[1], 0});
0164       frame.vecEPhi.Set(local[0], local[1]);
0165     }
0166     return frame;
0167   }
0168
0169   /// @brief Combine this mask (1,2) with another (3,4) and get the resulting extent.
0170   /// @param other Another ZPhi frame mask
0171   bool CombineWith(ZPhiMask<Real_t> const &other)
0172   {
0173     // Check compatibility of frames
0174     if (vecCore::math::Abs(invcalf - other.invcalf) > vecgeom::kToleranceDist<Real_t>) return false;
0175     if (vecCore::math::Abs(Radius(0) - other.Radius(0)) > vecgeom::kToleranceDist<Real_t>) return false;
0176
0177     rangeZ[0] = vecCore::math::Min(rangeZ[0], other.rangeZ[0]);
0178     rangeZ[1] = vecCore::math::Max(rangeZ[1], other.rangeZ[1]);
0179     isFullCirc |= other.isFullCirc;
0180     if (isFullCirc) return true;
0181     // intersect the phi ranges
0182     // Matching start-start (1==3)
0183     if (ApproxEqualVector2(vecSPhi, other.vecSPhi)) {
0184       if (!InsidePhi(other.vecEPhi[0], other.vecEPhi[1])) vecEPhi = other.vecEPhi;
0185       return true;
0186     }
0187     // Matching end-end (2==4)
0188     if (ApproxEqualVector2(vecEPhi, other.vecEPhi)) {
0189       if (!InsidePhi(other.vecSPhi[0], other.vecSPhi[1])) vecSPhi = other.vecSPhi;
0190       return true;
0191     }
0192     // Matching end-start ranges (2==3)
0193     if (ApproxEqualVector2(vecEPhi, other.vecSPhi)) {
0194       if (InsidePhi(other.vecEPhi[0], other.vecEPhi[1])) {
0195         vecEPhi    = vecSPhi;
0196         isFullCirc = true;
0197       } else {
0198         vecEPhi = other.vecEPhi;
0199       }
0200       return true;
0201     }
0202     // Matching start-end ranges (1==4)
0203     if (ApproxEqualVector2(vecSPhi, other.vecEPhi)) {
0204       if (InsidePhi(other.vecSPhi[0], other.vecSPhi[1])) {
0205         vecEPhi    = vecSPhi;
0206         isFullCirc = true;
0207       } else {
0208         vecSPhi = other.vecSPhi;
0209       }
0210       return true;
0211     }
0212     // non-matching ends
0213     bool in1 = other.InsidePhi(vecSPhi[0], vecSPhi[1]);       // 1 inside (3,4)
0214     bool in2 = other.InsidePhi(vecEPhi[0], vecEPhi[1]);       // 2 inside (3,4)
0215     bool in3 = InsidePhi(other.vecSPhi[0], other.vecSPhi[1]); // 3 inside (1,2)
0216     bool in4 = InsidePhi(other.vecEPhi[0], other.vecEPhi[1]); // 4 inside (1,2)
0217
0218     if (!(in1 || in2 || in3 || in4)) {
0219       if ((vecSPhi - other.vecEPhi).Mag2() > (other.vecSPhi - vecEPhi).Mag2())
0220         vecEPhi = other.vecEPhi;
0221       else
0222         vecSPhi = other.vecSPhi;
0223       return true;
0224     }
0225
0226     if (in1 && in2 && in3 && in4) {
0227       vecEPhi    = vecSPhi;
0228       isFullCirc = true;
0229       return true;
0230     }
0231
0232     if (!(in1 || in2)) return true;
0233
0234     if (!(in3 || in4)) {
0235       vecSPhi = other.vecSPhi;
0236       vecEPhi = other.vecEPhi;
0237       return true;
0238     }
0239
0240     if (!in1 && in2) {
0241       vecEPhi = other.vecEPhi;
0242       return true;
0243     }
0244
0245     if (in1 && !in2) {
0246       vecSPhi = other.vecSPhi;
0247       return true;
0248     }
0249     VECGEOM_LOG(critical) << "CombineWith: wrong logic";
0250     return false;
0251   }
0252
0253   /// @brief Computes safe distance to the frame combining surface and frame safeties
0254   /// @details The safety to the cylindrical surface comes as safetySurf, it is positive since only
0255   ///  points outside the cylindrical shell are considered.
0256   /// @param local Projected point in local coordinates
0257   /// @param safetySurf Safety from non-projected point to the frame support surface.
0258   /// @return Safety distance to the mask.
0259   VECCORE_ATT_HOST_DEVICE
0260   Real_t Safety(Vector3D<Real_t> const &local, Real_t safetySurf, bool &valid) const
0261   {
0262     valid = true;
0263     if (!InsidePhi(local[0], local[1])) {
0264       // If the point is not in the phi range, there are other surfaces closer than this one
0265       // so this frame should not be part of the minimization process.
0266       valid = false;
0267       return 0.;
0268     }
0269     auto safetyZ = vecCore::math::Max(local[2] - rangeZ[1], rangeZ[0] - local[2]);
0270     return vecCore::math::Max(safetySurf, safetyZ);
0271   }
0272 };
0273
0274 } // namespace vgbrep
0275
0276 #endif