volumes/kernel/TrdImplementation.h

0001 // This file is part of VecGeom and is distributed under the
0002 // conditions in the file LICENSE.txt in the top directory.
0003 // For the full list of authors see CONTRIBUTORS.txt and `git log`.
0004
0005 /// This file implements the algorithms for Trd
0006 /// @file volumes/kernel/TrdImplementation.h
0007 /// @author Georgios Bitzes
0008
0009 #ifndef VECGEOM_VOLUMES_KERNEL_TRDIMPLEMENTATION_H_
0010 #define VECGEOM_VOLUMES_KERNEL_TRDIMPLEMENTATION_H_
0011
0012 #include "VecGeom/base/Global.h"
0013 #include "VecGeom/volumes/kernel/GenericKernels.h"
0014 #include "VecGeom/volumes/TrdStruct.h"
0015 #include "VecGeom/volumes/kernel/shapetypes/TrdTypes.h"
0016 #include <stdlib.h>
0017 #include <cstdio>
0018
0019 namespace vecgeom {
0020
0021 VECGEOM_DEVICE_DECLARE_CONV_TEMPLATE(struct, TrdImplementation, typename);
0022
0023 inline namespace VECGEOM_IMPL_NAMESPACE {
0024
0025 namespace TrdUtilities {
0026
0027 /*
0028  * Checks whether a point (x, y) falls on the left or right half-plane
0029  * of a line. The line is defined by a (vx, vy) vector, extended to infinity.
0030  *
0031  * Of course this can only be used for lines that pass through (0, 0), but
0032  * you can supply transformed coordinates for the point to check for any line.
0033  *
0034  * This simply calculates the magnitude of the cross product of vectors (px, py)
0035  * and (vx, vy), which is defined as |x| * |v| * sin theta.
0036  *
0037  * If the cross product is positive, the point is clockwise of V, or the "right"
0038  * half-plane. If it's negative, the point is CCW and on the "left" half-plane.
0039  */
0040
0041 template <typename Real_v>
0042 VECGEOM_FORCE_INLINE
0043 VECCORE_ATT_HOST_DEVICE
0044 void PointLineOrientation(Real_v const &px, Real_v const &py, Precision const &vx, Precision const &vy,
0045                           Real_v &crossProduct)
0046 {
0047   crossProduct = vx * py - vy * px;
0048 }
0049
0050 /*
0051  * Check intersection of the trajectory of a particle with a segment
0052  * that's bound from -Ylimit to +Ylimit.j
0053  *
0054  * All points of the along-vector of a plane lie on
0055  * s * (alongX, alongY)
0056  * All points of the trajectory of the particle lie on
0057  * (x, y) + t * (vx, vy)
0058  * Thefore, it must hold that s * (alongX, alongY) == (x, y) + t * (vx, vy)
0059  * Solving by t we get t = (alongY*x - alongX*y) / (vy*alongX - vx*alongY)
0060  *
0061  * t gives the distance, but how to make sure hitpoint is inside the
0062  * segment and not just the infinite line defined by the segment?
0063  *
0064  * Check that |hity| <= Ylimit
0065  */
0066
0067 template <typename Real_v>
0068 VECGEOM_FORCE_INLINE
0069 VECCORE_ATT_HOST_DEVICE
0070 void PlaneTrajectoryIntersection(Real_v const &alongX, Real_v const &alongY, Real_v const &ylimit, Real_v const &posx,
0071                                  Real_v const &posy, Real_v const &dirx, Real_v const &diry, Real_v &dist,
0072                                  vecCore::Mask_v<Real_v> &ok)
0073 {
0074   dist = (alongY * posx - alongX * posy) / (diry * alongX - dirx * alongY);
0075
0076   Real_v hity = posy + dist * diry;
0077   ok          = vecCore::math::Abs(hity) <= ylimit && dist > 0;
0078 }
0079
0080 template <typename Real_v, bool forY, bool mirroredPoint, bool toInside>
0081 VECGEOM_FORCE_INLINE
0082 VECCORE_ATT_HOST_DEVICE
0083 void FaceTrajectoryIntersection(TrdStruct<Precision> const &trd, Vector3D<Real_v> const &pos,
0084                                 Vector3D<Real_v> const &dir, Real_v &dist, vecCore::Mask_v<Real_v> &ok)
0085 {
0086   Real_v alongV, posV, dirV, posK, dirK, fV, fK, halfKplus, v1, ndotv;
0087   //    fNormals[0].Set(-fCalfX, 0., fFx*fCalfX);
0088   //    fNormals[1].Set(fCalfX, 0., fFx*fCalfX);
0089   //    fNormals[2].Set(0., -fCalfY, fFy*fCalfY);
0090   //    fNormals[3].Set(0., fCalfY, fFy*fCalfY);
0091   if (forY) {
0092     alongV    = trd.fY2minusY1;
0093     v1        = trd.fDY1;
0094     posV      = pos.y();
0095     posK      = pos.x();
0096     dirV      = dir.y();
0097     dirK      = dir.x();
0098     fK        = trd.fFx;
0099     fV        = trd.fFy;
0100     halfKplus = trd.fHalfX1plusX2;
0101   } else {
0102     alongV    = trd.fX2minusX1;
0103     v1        = trd.fDX1;
0104     posV      = pos.x();
0105     posK      = pos.y();
0106     dirV      = dir.x();
0107     dirK      = dir.y();
0108     fK        = trd.fFy;
0109     fV        = trd.fFx;
0110     halfKplus = trd.fHalfY1plusY2;
0111   }
0112   if (mirroredPoint) {
0113     posV *= Real_v(-1.);
0114     dirV *= Real_v(-1.);
0115   }
0116
0117   ndotv = dirV + fV * dir.z();
0118   if (toInside)
0119     ok = ndotv < Real_v(0.);
0120   else
0121     ok = ndotv > Real_v(0.);
0122   if (vecCore::MaskEmpty(ok)) return;
0123   Real_v alongZ = Real_v(2.0) * trd.fDZ;
0124
0125   // distance from trajectory to face
0126   dist = (alongZ * (posV - v1) - alongV * (pos.z() + trd.fDZ)) / (dir.z() * alongV - dirV * alongZ + kTiny);
0127   ok &= dist > Real_v(MakeMinusTolerant<true>(0.));
0128   if (!vecCore::MaskEmpty(ok)) {
0129     // need to make sure z hit falls within bounds
0130     Real_v hitz = pos.z() + dist * dir.z();
0131     ok &= vecCore::math::Abs(hitz) <= trd.fDZ;
0132     // need to make sure hit on varying dimension falls within bounds
0133     Real_v hitk = posK + dist * dirK;
0134     Real_v dK   = halfKplus - fK * hitz; // calculate the width of the varying dimension at hitz
0135     ok &= vecCore::math::Abs(hitk) <= dK;
0136     vecCore::MaskedAssign(dist, ok & (vecCore::math::Abs(dist) < kHalfTolerance), Real_v(0.0));
0137   }
0138 }
0139
0140 template <typename Real_v, typename trdTypeT, bool inside>
0141 VECGEOM_FORCE_INLINE
0142 VECCORE_ATT_HOST_DEVICE
0143 void Safety(TrdStruct<Precision> const &trd, Vector3D<Real_v> const &pos, Real_v &dist)
0144 {
0145   using namespace TrdTypes;
0146   using Bool_v = vecCore::Mask_v<Real_v>;
0147
0148   Real_v safz = trd.fDZ - vecCore::math::Abs(pos.z());
0149   // std::cout << "safz: " << safz << std::endl;
0150   dist = safz;
0151
0152   Real_v distx = trd.fHalfX1plusX2 - trd.fFx * pos.z();
0153   Bool_v okx   = distx >= 0;
0154   Real_v safx  = (distx - vecCore::math::Abs(pos.x())) * trd.fCalfX;
0155   vecCore::MaskedAssign(dist, okx && safx < dist, safx);
0156   // std::cout << "safx: " << safx << std::endl;
0157
0158   if (checkVaryingY<trdTypeT>(trd)) {
0159     Real_v disty = trd.fHalfY1plusY2 - trd.fFy * pos.z();
0160     Bool_v oky   = disty >= 0;
0161     Real_v safy  = (disty - vecCore::math::Abs(pos.y())) * trd.fCalfY;
0162     vecCore::MaskedAssign(dist, oky && safy < dist, safy);
0163   } else {
0164     Real_v safy = trd.fDY1 - vecCore::math::Abs(pos.y());
0165     vecCore::MaskedAssign(dist, safy < dist, safy);
0166   }
0167   if (!inside) dist = -dist;
0168 }
0169
0170 template <typename Real_v, typename trdTypeT, bool surfaceT>
0171 VECGEOM_FORCE_INLINE
0172 VECCORE_ATT_HOST_DEVICE
0173 static void UnplacedInside(TrdStruct<Precision> const &trd, Vector3D<Real_v> const &point,
0174                            vecCore::Mask_v<Real_v> &completelyinside, vecCore::Mask_v<Real_v> &completelyoutside)
0175 {
0176
0177   using namespace TrdUtilities;
0178   using namespace TrdTypes;
0179
0180   Real_v pzPlusDz = point.z() + trd.fDZ;
0181
0182   // inside Z?
0183   completelyoutside = vecCore::math::Abs(point.z()) > MakePlusTolerant<surfaceT>(trd.fDZ);
0184   if (surfaceT) completelyinside = vecCore::math::Abs(point.z()) < MakeMinusTolerant<surfaceT>(trd.fDZ);
0185
0186   // inside X?
0187   Real_v cross;
0188   // Note: we cannot compare directly the cross product with the surface tolerance, but with
0189   // the tolerance multiplied by the length of the lateral segment connecting dx1 and dx2
0190   PointLineOrientation<Real_v>(vecCore::math::Abs(point.x()) - trd.fDX1, pzPlusDz, trd.fX2minusX1, 2.0 * trd.fDZ,
0191                                cross);
0192   if (surfaceT) {
0193     completelyoutside |= cross < -trd.fToleranceX;
0194     completelyinside &= cross > trd.fToleranceX;
0195   } else {
0196     completelyoutside |= cross < 0;
0197   }
0198
0199   // inside Y?
0200   if (HasVaryingY<trdTypeT>::value != TrdTypes::kNo) {
0201     // If Trd type is unknown don't bother with a runtime check, assume the general case
0202     PointLineOrientation<Real_v>(vecCore::math::Abs(point.y()) - trd.fDY1, pzPlusDz, trd.fY2minusY1, 2.0 * trd.fDZ,
0203                                  cross);
0204     if (surfaceT) {
0205       completelyoutside |= cross < -trd.fToleranceY;
0206       completelyinside &= cross > trd.fToleranceY;
0207     } else {
0208       completelyoutside |= cross < 0;
0209     }
0210   } else {
0211     completelyoutside |= vecCore::math::Abs(point.y()) > MakePlusTolerant<surfaceT>(trd.fDY1);
0212     if (surfaceT) completelyinside &= vecCore::math::Abs(point.y()) < MakeMinusTolerant<surfaceT>(trd.fDY1);
0213   }
0214 }
0215
0216 } // namespace TrdUtilities
0217
0218 template <typename T>
0219 class SPlacedTrd;
0220 template <typename T>
0221 class SUnplacedTrd;
0222
0223 template <typename T>
0224 struct TrdStruct;
0225
0226 template <typename trdTypeT>
0227 struct TrdImplementation {
0228
0229   using UnplacedStruct_t = TrdStruct<Precision>;
0230   using UnplacedVolume_t = SUnplacedTrd<trdTypeT>;
0231   using PlacedShape_t    = SPlacedTrd<UnplacedVolume_t>;
0232
0233   VECCORE_ATT_HOST_DEVICE
0234   static void PrintType() {}
0235
0236   template <typename Stream>
0237   static void PrintType(Stream &s, int transCodeT = translation::kGeneric, int rotCodeT = rotation::kGeneric)
0238   {
0239     s << "SpecializedTrd<" << transCodeT << "," << rotCodeT << ">";
0240   }
0241
0242   template <typename Stream>
0243   static void PrintImplementationType(Stream & /*s*/)
0244   {
0245   }
0246
0247   template <typename Stream>
0248   static void PrintUnplacedType(Stream & /*s*/)
0249   {
0250   }
0251
0252   template <typename Real_v, typename Bool_v>
0253   VECGEOM_FORCE_INLINE
0254   VECCORE_ATT_HOST_DEVICE
0255   static void UnplacedContains(UnplacedStruct_t const &trd, Vector3D<Real_v> const &point, Bool_v &inside)
0256   {
0257
0258     Bool_v unused;
0259     TrdUtilities::UnplacedInside<Real_v, trdTypeT, false>(trd, point, unused, inside);
0260     inside = !inside;
0261   }
0262
0263   template <typename Real_v, typename Bool_v>
0264   VECGEOM_FORCE_INLINE
0265   VECCORE_ATT_HOST_DEVICE
0266   static void Contains(UnplacedStruct_t const &trd, Vector3D<Real_v> const &point, Bool_v &inside)
0267   {
0268
0269     Bool_v unused;
0270     TrdUtilities::UnplacedInside<Real_v, trdTypeT, false>(trd, point, unused, inside);
0271     inside = !inside;
0272   }
0273
0274   template <typename Real_v, typename Inside_v>
0275   VECGEOM_FORCE_INLINE
0276   VECCORE_ATT_HOST_DEVICE
0277   static void Inside(UnplacedStruct_t const &trd, Vector3D<Real_v> const &point, Inside_v &inside)
0278   {
0279     // use double-based vector for result, as Bool_v is a mask for precision_v
0280     using Bool_v = vecCore::Mask_v<Real_v>;
0281     const Real_v in(EInside::kInside);
0282     const Real_v out(EInside::kOutside);
0283     Bool_v inmask  = Bool_v(false);
0284     Bool_v outmask = Bool_v(false);
0285     Real_v result(EInside::kSurface);
0286
0287     TrdUtilities::UnplacedInside<Real_v, trdTypeT, true>(trd, point, inmask, outmask);
0288
0289     vecCore::MaskedAssign(result, inmask, in);
0290     vecCore::MaskedAssign(result, outmask, out);
0291
0292     // Manual conversion from double to int here is necessary because int_v and
0293     // precision_v have different number of elements in SIMD vector, so Bool_v
0294     // (mask for precision_v) cannot be cast to mask for inside, which is a
0295     // different type and does not exist in the current backend system
0296     for (size_t i = 0; i < vecCore::VectorSize(result); i++)
0297       vecCore::Set(inside, i, vecCore::Get(result, i));
0298   }
0299
0300   template <typename Real_v>
0301   VECGEOM_FORCE_INLINE
0302   VECCORE_ATT_HOST_DEVICE
0303   static void DistanceToIn(UnplacedStruct_t const &trd, Vector3D<Real_v> const &point,
0304                            Vector3D<Real_v> const &direction, Real_v const & /*stepMax*/, Real_v &distance)
0305   {
0306
0307     using namespace TrdUtilities;
0308     using namespace TrdTypes;
0309     using Bool_v = vecCore::Mask_v<Real_v>;
0310
0311     Real_v hitx, hity;
0312     // Real_v hitz;
0313
0314     Vector3D<Real_v> pos_local;
0315     Vector3D<Real_v> dir_local;
0316     distance = InfinityLength<Real_v>();
0317
0318     // hit Z faces?
0319     Bool_v inz   = vecCore::math::Abs(point.z()) < Real_v(MakeMinusTolerant<true>(trd.fDZ));
0320     Real_v distx = trd.fHalfX1plusX2 - trd.fFx * point.z();
0321     Bool_v inx   = (distx - vecCore::math::Abs(point.x())) * trd.fCalfX > Real_v(MakePlusTolerant<true>(0.));
0322     Real_v disty;
0323     Bool_v iny;
0324     if (checkVaryingY<trdTypeT>(trd)) {
0325       disty = trd.fHalfY1plusY2 - trd.fFy * point.z();
0326       iny   = (disty - vecCore::math::Abs(point.y())) * trd.fCalfY > Real_v(MakePlusTolerant<true>(0.));
0327     } else {
0328       disty = vecCore::math::Abs(point.y()) - trd.fDY1;
0329       iny   = disty < Real_v(MakeMinusTolerant<true>(0.));
0330     }
0331     Bool_v inside = inx & iny & inz;
0332     vecCore__MaskedAssignFunc(distance, inside, Real_v(-1.));
0333     Bool_v done = inside;
0334     Bool_v okz  = point.z() * direction.z() < Real_v(0.);
0335     okz &= !inz;
0336     if (!vecCore::MaskEmpty(okz)) {
0337       Real_v distz = (vecCore::math::Abs(point.z()) - trd.fDZ) / vecCore::math::Abs(direction.z());
0338       // exclude case in which particle is going away
0339       hitx = vecCore::math::Abs(point.x() + distz * direction.x());
0340       hity = vecCore::math::Abs(point.y() + distz * direction.y());
0341
0342       // hitting top face?
0343       Bool_v okzt = point.z() > (trd.fDZ - kHalfTolerance) && hitx <= trd.fDX2 && hity <= trd.fDY2;
0344       // hitting bottom face?
0345       Bool_v okzb = point.z() < (-trd.fDZ + kHalfTolerance) && hitx <= trd.fDX1 && hity <= trd.fDY1;
0346
0347       okz &= (okzt | okzb);
0348       vecCore::MaskedAssign(distance, okz, distz);
0349     }
0350     done |= okz;
0351     if (vecCore::MaskFull(done)) {
0352       vecCore::MaskedAssign(distance, vecCore::math::Abs(distance) < kHalfTolerance, Real_v(0.0));
0353       return;
0354     }
0355
0356     // hitting X faces?
0357     Bool_v okx = Bool_v(false);
0358     if (!vecCore::MaskFull(inx)) {
0359
0360       FaceTrajectoryIntersection<Real_v, false, false, true>(trd, point, direction, distx, okx);
0361       vecCore::MaskedAssign(distance, okx, distx);
0362
0363       FaceTrajectoryIntersection<Real_v, false, true, true>(trd, point, direction, distx, okx);
0364       vecCore::MaskedAssign(distance, okx, distx);
0365     }
0366     done |= okx;
0367     if (vecCore::MaskFull(done)) {
0368       vecCore::MaskedAssign(distance, vecCore::math::Abs(distance) < kHalfTolerance, Real_v(0.0));
0369       return;
0370     }
0371
0372     // hitting Y faces?
0373     Bool_v oky;
0374     if (checkVaryingY<trdTypeT>(trd)) {
0375       if (!vecCore::MaskFull(iny)) {
0376         FaceTrajectoryIntersection<Real_v, true, false, true>(trd, point, direction, disty, oky);
0377         vecCore::MaskedAssign(distance, oky, disty);
0378
0379         FaceTrajectoryIntersection<Real_v, true, true, true>(trd, point, direction, disty, oky);
0380         vecCore::MaskedAssign(distance, oky, disty);
0381       }
0382     } else {
0383       if (!vecCore::MaskFull(iny)) {
0384         disty /= vecCore::math::Abs(direction.y());
0385         Real_v zhit = point.z() + disty * direction.z();
0386         Real_v xhit = point.x() + disty * direction.x();
0387         Real_v dx   = trd.fHalfX1plusX2 - trd.fFx * zhit;
0388         oky         = point.y() * direction.y() < 0 && disty > -kHalfTolerance && vecCore::math::Abs(xhit) < dx &&
0389               vecCore::math::Abs(zhit) < trd.fDZ;
0390         vecCore::MaskedAssign(distance, oky, disty);
0391       }
0392     }
0393     vecCore::MaskedAssign(distance, vecCore::math::Abs(distance) < kHalfTolerance, Real_v(0.0));
0394   }
0395
0396   template <typename Real_v>
0397   VECGEOM_FORCE_INLINE
0398   VECCORE_ATT_HOST_DEVICE
0399   static void DistanceToOut(UnplacedStruct_t const &trd, Vector3D<Real_v> const &point, Vector3D<Real_v> const &dir,
0400                             Real_v const & /*stepMax*/, Real_v &distance)
0401   {
0402
0403     using namespace TrdUtilities;
0404     using namespace TrdTypes;
0405     using Bool_v = vecCore::Mask_v<Real_v>;
0406
0407     Real_v hitx, hity;
0408     // Real_v hitz;
0409     distance = Real_v(0.0);
0410
0411     // hit top Z face?
0412     Real_v invdir = Real_v(1.) / vecCore::math::Abs(dir.z() + kTiny);
0413     Real_v safz   = trd.fDZ - vecCore::math::Abs(point.z());
0414     Bool_v out    = safz < Real_v(MakeMinusTolerant<true>(0.));
0415     Real_v distx  = trd.fHalfX1plusX2 - trd.fFx * point.z();
0416     out |= (distx - vecCore::math::Abs(point.x())) * trd.fCalfX < Real_v(MakeMinusTolerant<true>(0.));
0417     Real_v disty;
0418     if (checkVaryingY<trdTypeT>(trd)) {
0419       disty = trd.fHalfY1plusY2 - trd.fFy * point.z();
0420       out |= (disty - vecCore::math::Abs(point.y())) * trd.fCalfY < Real_v(MakeMinusTolerant<true>(0.));
0421     } else {
0422       disty = trd.fDY1 - vecCore::math::Abs(point.y());
0423       out |= disty < Real_v(MakeMinusTolerant<true>(0.));
0424     }
0425     if (/*vecCore::EarlyReturnAllowed() && */ vecCore::MaskFull(out)) {
0426       distance = Real_v(-1.);
0427       return;
0428     }
0429     Bool_v okzt = dir.z() > Real_v(0.);
0430     if (!vecCore::MaskEmpty(okzt)) {
0431       Real_v distz = (trd.fDZ - point.z()) * invdir;
0432       hitx         = vecCore::math::Abs(point.x() + distz * dir.x());
0433       hity         = vecCore::math::Abs(point.y() + distz * dir.y());
0434       okzt &= hitx <= trd.fDX2 && hity <= trd.fDY2;
0435       vecCore::MaskedAssign(distance, okzt, distz);
0436       if (vecCore::EarlyReturnAllowed() && vecCore::MaskFull(okzt)) {
0437         vecCore::MaskedAssign(distance, vecCore::math::Abs(distance) < kHalfTolerance, Real_v(0.0));
0438         return;
0439       }
0440     }
0441
0442     // hit bottom Z face?
0443     Bool_v okzb = dir.z() < Real_v(0.);
0444     if (!vecCore::MaskEmpty(okzb)) {
0445       Real_v distz = (point.z() + trd.fDZ) * invdir;
0446       hitx         = vecCore::math::Abs(point.x() + distz * dir.x());
0447       hity         = vecCore::math::Abs(point.y() + distz * dir.y());
0448       okzb &= hitx <= trd.fDX1 && hity <= trd.fDY1;
0449       vecCore::MaskedAssign(distance, okzb, distz);
0450       if (vecCore::EarlyReturnAllowed() && vecCore::MaskFull(okzb)) {
0451         vecCore::MaskedAssign(distance, vecCore::math::Abs(distance) < kHalfTolerance, Real_v(0.0));
0452         return;
0453       }
0454     }
0455
0456     // hitting X faces?
0457     Bool_v okx;
0458
0459     FaceTrajectoryIntersection<Real_v, false, false, false>(trd, point, dir, distx, okx);
0460
0461     vecCore::MaskedAssign(distance, okx, distx);
0462     if (vecCore::EarlyReturnAllowed() && vecCore::MaskFull(okx)) {
0463       vecCore::MaskedAssign(distance, vecCore::math::Abs(distance) < kHalfTolerance, Real_v(0.0));
0464       return;
0465     }
0466
0467     FaceTrajectoryIntersection<Real_v, false, true, false>(trd, point, dir, distx, okx);
0468     vecCore::MaskedAssign(distance, okx, distx);
0469     if (vecCore::EarlyReturnAllowed() && vecCore::MaskFull(okx)) {
0470       vecCore::MaskedAssign(distance, vecCore::math::Abs(distance) < kHalfTolerance, Real_v(0.0));
0471       return;
0472     }
0473
0474     // hitting Y faces?
0475     Bool_v oky;
0476
0477     if (checkVaryingY<trdTypeT>(trd)) {
0478       FaceTrajectoryIntersection<Real_v, true, false, false>(trd, point, dir, disty, oky);
0479       vecCore::MaskedAssign(distance, oky, disty);
0480       if (vecCore::EarlyReturnAllowed() && vecCore::MaskFull(oky)) {
0481         vecCore::MaskedAssign(distance, vecCore::math::Abs(distance) < kHalfTolerance, Real_v(0.0));
0482         return;
0483       }
0484
0485       FaceTrajectoryIntersection<Real_v, true, true, false>(trd, point, dir, disty, oky);
0486       vecCore::MaskedAssign(distance, oky, disty);
0487     } else {
0488       Real_v plane = trd.fDY1;
0489       vecCore__MaskedAssignFunc(plane, dir.y() < Real_v(0.), Real_v(-trd.fDY1));
0490       disty       = (plane - point.y()) / dir.y();
0491       Real_v zhit = point.z() + disty * dir.z();
0492       Real_v xhit = point.x() + disty * dir.x();
0493       Real_v dx   = trd.fHalfX1plusX2 - trd.fFx * zhit;
0494       oky         = vecCore::math::Abs(xhit) < dx && vecCore::math::Abs(zhit) < trd.fDZ;
0495       vecCore::MaskedAssign(distance, oky, disty);
0496     }
0497     vecCore__MaskedAssignFunc(distance, vecCore::math::Abs(distance) < kHalfTolerance, Real_v(0.0));
0498     vecCore__MaskedAssignFunc(distance, out, Real_v(-1.0));
0499   }
0500
0501   template <typename Real_v>
0502   VECGEOM_FORCE_INLINE
0503   VECCORE_ATT_HOST_DEVICE
0504   static void SafetyToIn(UnplacedStruct_t const &trd, Vector3D<Real_v> const &point, Real_v &safety)
0505   {
0506     using namespace TrdUtilities;
0507     Safety<Real_v, trdTypeT, false>(trd, point, safety);
0508   }
0509
0510   template <typename Real_v>
0511   VECGEOM_FORCE_INLINE
0512   VECCORE_ATT_HOST_DEVICE
0513   static void SafetyToOut(UnplacedStruct_t const &trd, Vector3D<Real_v> const &point, Real_v &safety)
0514   {
0515     using namespace TrdUtilities;
0516     Safety<Real_v, trdTypeT, true>(trd, point, safety);
0517   }
0518 };
0519 } // namespace VECGEOM_IMPL_NAMESPACE
0520 } // namespace vecgeom
0521
0522 #endif // VECGEOM_VOLUMES_KERNEL_TRDIMPLEMENTATION_H_