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