volumes/kernel/TorusImplementation2.h

0001 /// @file TorusImplementation2.h
0002
0003 #ifndef VECGEOM_VOLUMES_KERNEL_TORUSIMPLEMENTATION2_H_
0004 #define VECGEOM_VOLUMES_KERNEL_TORUSIMPLEMENTATION2_H_
0005
0006 #include "VecGeom/base/Global.h"
0007 #include "VecGeom/base/Transformation3D.h"
0008 #include "VecGeom/volumes/kernel/GenericKernels.h"
0009 #include "VecGeom/volumes/kernel/TubeImplementation.h"
0010 #include "VecGeom/volumes/TorusStruct2.h"
0011
0012 #include <cstdio>
0013 #include <VecCore/VecCore>
0014
0015 namespace vecgeom {
0016
0017 VECGEOM_DEVICE_FORWARD_DECLARE(struct TorusImplementation2;);
0018 VECGEOM_DEVICE_DECLARE_CONV(struct, TorusImplementation2);
0019
0020 inline namespace VECGEOM_IMPL_NAMESPACE {
0021
0022 //_____________________________________________________________________________
0023 template <typename T>
0024 VECCORE_ATT_HOST_DEVICE unsigned int SolveCubic(T a, T b, T c, T *x)
0025 {
0026   // Find real solutions of the cubic equation : x^3 + a*x^2 + b*x + c = 0
0027   // Input: a,b,c
0028   // Output: x[3] real solutions
0029   // Returns number of real solutions (1 or 3)
0030   const T ott        = 1. / 3.;
0031   const T sq3        = Sqrt(3.);
0032   const T inv6sq3    = 1. / (6. * sq3);
0033   unsigned int ireal = 1;
0034   T p                = b - a * a * ott;
0035   T q                = c - a * b * ott + 2. * a * a * a * ott * ott * ott;
0036   T delta            = 4 * p * p * p + 27. * q * q;
0037   T t, u;
0038
0039   if (delta >= 0) {
0040     delta = Sqrt(delta);
0041     t     = (-3 * q * sq3 + delta) * inv6sq3;
0042     u     = (3 * q * sq3 + delta) * inv6sq3;
0043     x[0]  = CopySign(T(1.), t) * Cbrt(Abs(t)) - CopySign(T(1.), u) * Cbrt(Abs(u)) - a * ott;
0044   } else {
0045     delta = Sqrt(-delta);
0046     t     = -0.5 * q;
0047     u     = delta * inv6sq3;
0048     x[0]  = 2. * Pow(t * t + u * u, T(0.5) * ott) * cos(ott * ATan2(u, t));
0049     x[0] -= a * ott;
0050   }
0051
0052   t     = x[0] * x[0] + a * x[0] + b;
0053   u     = a + x[0];
0054   delta = u * u - T(4.) * t;
0055   if (delta >= 0) {
0056     ireal = 3;
0057     delta = Sqrt(delta);
0058     x[1]  = T(0.5) * (-u - delta);
0059     x[2]  = T(0.5) * (-u + delta);
0060   }
0061
0062   return ireal;
0063 }
0064
0065 template <typename T, unsigned int i, unsigned int j>
0066 VECGEOM_FORCE_INLINE VECCORE_ATT_HOST_DEVICE void CmpAndSwap(T *array)
0067 {
0068   if (vecCore::MaskFull(array[i] > array[j])) {
0069     T c      = array[j];
0070     array[j] = array[i];
0071     array[i] = c;
0072   }
0073 }
0074
0075 // a special function to sort a 4 element array
0076 // sorting is done inplace and in increasing order
0077 // implementation comes from a sorting network
0078 template <typename T>
0079 VECGEOM_FORCE_INLINE VECCORE_ATT_HOST_DEVICE void Sort4(T *array)
0080 {
0081   CmpAndSwap<T, 0, 2>(array);
0082   CmpAndSwap<T, 1, 3>(array);
0083   CmpAndSwap<T, 0, 1>(array);
0084   CmpAndSwap<T, 2, 3>(array);
0085   CmpAndSwap<T, 1, 2>(array);
0086 }
0087
0088 // solve quartic taken from ROOT/TGeo and adapted
0089 //_____________________________________________________________________________
0090
0091 template <typename T>
0092 VECCORE_ATT_HOST_DEVICE int SolveQuartic(T a, T b, T c, T d, T *x)
0093 {
0094   // Find real solutions of the quartic equation : x^4 + a*x^3 + b*x^2 + c*x + d = 0
0095   // Input: a,b,c,d
0096   // Output: x[4] - real solutions
0097   // Returns number of real solutions (0 to 3)
0098   T e     = b - 3. * a * a / 8.;
0099   T f     = c + a * a * a / 8. - 0.5 * a * b;
0100   T g     = d - 3. * a * a * a * a / 256. + a * a * b / 16. - a * c / 4.;
0101   T xx[4] = {vecgeom::kInfLength, vecgeom::kInfLength, vecgeom::kInfLength, vecgeom::kInfLength};
0102   T delta;
0103   T h                = 0.;
0104   unsigned int ireal = 0;
0105
0106   // special case when f is zero
0107   if (Abs(f) < 1E-6) {
0108     delta = e * e - 4. * g;
0109     if (delta < 0.) return 0;
0110     delta = Sqrt(delta);
0111     h     = 0.5 * (-e - delta);
0112     if (h >= 0) {
0113       h          = Sqrt(h);
0114       x[ireal++] = -h - 0.25 * a;
0115       x[ireal++] = h - 0.25 * a;
0116     }
0117     h = 0.5 * (-e + delta);
0118     if (h >= 0) {
0119       h          = Sqrt(h);
0120       x[ireal++] = -h - 0.25 * a;
0121       x[ireal++] = h - 0.25 * a;
0122     }
0123     Sort4(x);
0124     return ireal;
0125   }
0126
0127   if (Abs(g) < 1E-6) {
0128     x[ireal++] = -0.25 * a;
0129     // this actually wants to solve a second order equation
0130     // we should specialize if it happens often
0131     unsigned int ncubicroots = SolveCubic<T>(0, e, f, xx);
0132     // this loop is not nice
0133     for (unsigned int i = 0; i < ncubicroots; i++)
0134       x[ireal++] = xx[i] - 0.25 * a;
0135     Sort4(x); // could be Sort3
0136     return ireal;
0137   }
0138
0139   ireal = SolveCubic<T>(2. * e, e * e - 4. * g, -f * f, xx);
0140   if (ireal == 1) {
0141     if (xx[0] <= 0) return 0;
0142     h = Sqrt(xx[0]);
0143   } else {
0144     // 3 real solutions of the cubic
0145     for (unsigned int i = 0; i < 3; i++) {
0146       h = xx[i];
0147       if (h >= 0) break;
0148     }
0149     if (h <= 0) return 0;
0150     h = Sqrt(h);
0151   }
0152   T j   = 0.5 * (e + h * h - f / h);
0153   ireal = 0;
0154   delta = h * h - 4. * j;
0155   if (delta >= 0) {
0156     delta      = Sqrt(delta);
0157     x[ireal++] = 0.5 * (-h - delta) - 0.25 * a;
0158     x[ireal++] = 0.5 * (-h + delta) - 0.25 * a;
0159   }
0160   delta = h * h - 4. * g / j;
0161   if (delta >= 0) {
0162     delta      = Sqrt(delta);
0163     x[ireal++] = 0.5 * (h - delta) - 0.25 * a;
0164     x[ireal++] = 0.5 * (h + delta) - 0.25 * a;
0165   }
0166   Sort4(x);
0167   return ireal;
0168 }
0169
0170 class PlacedTorus2;
0171 template <typename T>
0172 struct TorusStruct2;
0173 class UnplacedTorus2;
0174
0175 class SIMDUnplacedTorus2;
0176
0177 struct TorusImplementation2 {
0178   using PlacedShape_t    = PlacedTorus2;
0179   using UnplacedStruct_t = TorusStruct2<Precision>;
0180   using UnplacedVolume_t = UnplacedTorus2;
0181
0182   template <class Real_v>
0183   VECCORE_ATT_HOST_DEVICE static Real_v DistSqrToTorusR(Vector3D<Real_v> const &point,
0184                                                         Vector3D<Real_v> const &dir, Real_v dist)
0185   {
0186     Vector3D<Real_v> p = point + dir * dist;
0187     Real_v rxy         = p.Perp();
0188     return (rxy - 1.) * (rxy - 1.) + p.z() * p.z();
0189   }
0190
0191   template <typename Real_v>
0192   VECGEOM_FORCE_INLINE VECCORE_ATT_HOST_DEVICE static void DistanceToOut(UnplacedStruct_t const &torus,
0193                                                                          Vector3D<Real_v> const &point,
0194                                                                          Vector3D<Real_v> const &dir,
0195                                                                          Real_v const & /*stepMax*/, Real_v &distance)
0196   {
0197     using Inside_v = vecCore::Index_v<Real_v>;
0198     using Bool_v   = vecCore::Mask_v<Real_v>;
0199
0200     bool hasphi  = (torus.dphi() < kTwoPi);
0201     bool hasrmin = (torus.rmin() > 0);
0202
0203     //=== First, for points outside --> return infinity
0204     Bool_v done = Bool_v(false);
0205     distance    = kInfLength;
0206
0207     // very simple calculations -- only if can save some time
0208     Real_v distz = Abs(point.z()) - torus.rmax();
0209     done |= distz > kHalfTolerance;
0210
0211     // outside of bounding tube?
0212     Real_v rsq = point.x() * point.x() + point.y() * point.y();
0213     // Real_v rdotv = point.x()*dir.x() + point.y()*dir.y();
0214     Precision outerExclRadius = torus.rtor() + torus.rmax() + kHalfTolerance;
0215     done |= rsq > outerExclRadius * outerExclRadius;
0216     Precision innerExclRadius = torus.rtor() - torus.rmax() - kHalfTolerance;
0217     done |= rsq < innerExclRadius * innerExclRadius;
0218     vecCore__MaskedAssignFunc(distance, done, Real_v(-1.));
0219     if (vecCore::EarlyReturnAllowed() && vecCore::MaskFull(done)) return;
0220
0221     //=== Use InsideKernel() for a quick check, and if outside --> return -1
0222     // Bool_t inside=false, outside=false;
0223     // GenericKernelForContainsAndInside<Backend,true,true>(torus, point, inside, outside);
0224     // MaskedAssign( inside, -1.0, &distance );
0225     // done |= inside;
0226     Inside_v locus;
0227     TorusImplementation2::InsideKernel<Real_v, Inside_v>(torus, point, locus);
0228     vecCore__MaskedAssignFunc(distance, locus == EInside::kOutside, Real_v(-1.));
0229     done |= locus == EInside::kOutside;
0230     if (vecCore::EarlyReturnAllowed() && vecCore::MaskFull(done)) return;
0231
0232     Real_v dout = ToBoundary<Real_v, false>(torus, point/torus.rtor(), dir, torus.rmax()/torus.rtor(), true);
0233     // ToBoundary<Backend, false, true>(torus, point, dir, torus.rmax());
0234     Real_v din(kInfLength);
0235     if (hasrmin) {
0236       din = ToBoundary<Real_v, true>(torus, point/torus.rtor(), dir, torus.rmin()/torus.rtor(), true);
0237       // ToBoundary<Backend, true, true>(torus, point, dir, torus.rmin());
0238     }
0239     distance = Min(dout, din);
0240     distance *= torus.rtor();
0241     // std::cerr << "dout, din: " << dout << ", " << din << '\n';
0242     // std::cerr << "distance = Min(dout, din): " << distance << '\n';
0243
0244     if (hasphi) {
0245       Real_v distPhi1;
0246       Real_v distPhi2;
0247       // torus.GetWedge().DistanceToOut<Backend>(point, dir, distPhi1, distPhi2);
0248       torus.GetWedge().DistanceToOut<Real_v>(point, dir, distPhi1, distPhi2);
0249       Bool_v smallerphi = distPhi1 < distance;
0250       if (!vecCore::MaskEmpty(smallerphi)) {
0251         Vector3D<Real_v> intersectionPoint = point + dir * distPhi1;
0252         Bool_v insideDisk;
0253         // UnplacedContainsDisk<Backend>(torus, intersectionPoint, insideDisk);
0254         UnplacedContainsDisk<Real_v, Bool_v>(torus, intersectionPoint, insideDisk);
0255
0256         if (!vecCore::MaskEmpty(insideDisk)) // Inside Disk
0257         {
0258           Real_v diri = intersectionPoint.x() * torus.GetWedge().GetAlong1().x() +
0259                         intersectionPoint.y() * torus.GetWedge().GetAlong1().y();
0260           Bool_v rightside = (diri >= 0);
0261
0262           vecCore__MaskedAssignFunc(distance, rightside && smallerphi && insideDisk, distPhi1);
0263         }
0264       }
0265       smallerphi = distPhi2 < distance;
0266       if (!vecCore::MaskEmpty(smallerphi)) {
0267
0268         Vector3D<Real_v> intersectionPoint = point + dir * distPhi2;
0269         Bool_v insideDisk;
0270         // UnplacedContainsDisk<Backend>(torus, intersectionPoint, insideDisk);
0271         UnplacedContainsDisk<Real_v, Bool_v>(torus, intersectionPoint, insideDisk);
0272         if (!vecCore::MaskEmpty(insideDisk)) // Inside Disk
0273         {
0274           Real_v diri2 = intersectionPoint.x() * torus.GetWedge().GetAlong2().x() +
0275                          intersectionPoint.y() * torus.GetWedge().GetAlong2().y();
0276           Bool_v rightside = (diri2 >= Real_v(0));
0277           vecCore__MaskedAssignFunc(distance, rightside && (distPhi2 < distance) && smallerphi && insideDisk, distPhi2);
0278         }
0279       }
0280     }
0281
0282     vecCore__MaskedAssignFunc(distance, done || distance >= kInfLength, Real_v(-1.));
0283   }
0284
0285   template <typename Real_v, typename Bool_v, bool notForDisk>
0286   VECGEOM_FORCE_INLINE VECCORE_ATT_HOST_DEVICE static void ContainsKernel(UnplacedStruct_t const &torus,
0287                                                                           Vector3D<Real_v> const &point, Bool_v &inside)
0288   {
0289     Bool_v unused(false);
0290     Bool_v outside(false);
0291     TorusImplementation2::GenericKernelForContainsAndInside<Real_v, false, notForDisk>(torus, point, unused, outside);
0292     inside = !outside;
0293   }
0294
0295   template <typename Real_v, typename Bool_v>
0296   VECGEOM_FORCE_INLINE VECCORE_ATT_HOST_DEVICE static void UnplacedContainsDisk(UnplacedStruct_t const &torus,
0297                                                                                 Vector3D<Real_v> const &point,
0298                                                                                 Bool_v &inside)
0299   {
0300     ContainsKernel<Real_v, Bool_v, false>(torus, point, inside);
0301   }
0302
0303   template <typename Real_v, typename Inside_t>
0304   VECGEOM_FORCE_INLINE VECCORE_ATT_HOST_DEVICE static void InsideKernel(UnplacedStruct_t const &torus,
0305                                                                         Vector3D<Real_v> const &point, Inside_t &inside)
0306   {
0307
0308     using Bool_v = vecCore::Mask_v<Real_v>;
0309     //
0310     Bool_v completelyinside, completelyoutside;
0311     TorusImplementation2::GenericKernelForContainsAndInside<Real_v, true, true>(torus, point, completelyinside,
0312                                                                                 completelyoutside);
0313     inside = Inside_t(EInside::kSurface);
0314     vecCore::MaskedAssign(inside, completelyoutside, Inside_t(EInside::kOutside));
0315     vecCore::MaskedAssign(inside, completelyinside, Inside_t(EInside::kInside));
0316   }
0317
0318   template <typename Real_v, bool ForInside, bool notForDisk>
0319   VECGEOM_FORCE_INLINE VECCORE_ATT_HOST_DEVICE static void GenericKernelForContainsAndInside(
0320       UnplacedStruct_t const &torus, Vector3D<Real_v> const &point, typename vecCore::Mask_v<Real_v> &completelyinside,
0321       typename vecCore::Mask_v<Real_v> &completelyoutside)
0322
0323   {
0324     // using vecgeom::GenericKernels;
0325     // here we are explicitely unrolling the loop since  a for statement will likely be a penality
0326     // check if second call to Abs is compiled away
0327     // and it can anyway not be vectorized
0328     /* rmax */
0329     using Bool_v                = vecCore::Mask_v<Real_v>;
0330     VECGEOM_CONST Precision tol = 100. * vecgeom::kTolerance;
0331
0332     Real_v rxy   = Sqrt(point[0] * point[0] + point[1] * point[1]);
0333     Real_v radsq = (rxy - torus.rtor()) * (rxy - torus.rtor()) + point[2] * point[2];
0334
0335     if (ForInside) {
0336       completelyoutside = radsq > (tol * torus.rmax() + torus.rmax2()); // rmax
0337       completelyinside  = radsq < (-tol * torus.rmax() + torus.rmax2());
0338     } else {
0339       completelyoutside = radsq > torus.rmax2();
0340     }
0341
0342     if (vecCore::EarlyReturnAllowed()) {
0343       if (vecCore::MaskFull(completelyoutside)) {
0344         return;
0345       }
0346     }
0347     /* rmin */
0348     if (ForInside) {
0349       completelyoutside |= radsq < (-tol * torus.rmin() + torus.rmin2()); // rmin
0350       completelyinside &= radsq > (tol * torus.rmin() + torus.rmin2());
0351     } else {
0352       completelyoutside |= radsq < torus.rmin2();
0353     }
0354
0355     // NOT YET NEEDED WHEN NOT PHI TREATMENT
0356     if (vecCore::EarlyReturnAllowed()) {
0357       if (vecCore::MaskFull(completelyoutside)) {
0358         return;
0359       }
0360     }
0361
0362     /* phi */
0363     if ((torus.dphi() < kTwoPi) && (notForDisk)) {
0364       Bool_v completelyoutsidephi;
0365       Bool_v completelyinsidephi;
0366       torus.GetWedge().GenericKernelForContainsAndInside<Real_v, ForInside>(point, completelyinsidephi,
0367                                                                             completelyoutsidephi);
0368
0369       completelyoutside |= completelyoutsidephi;
0370       if (ForInside) completelyinside &= completelyinsidephi;
0371     }
0372   }
0373
0374   template <typename Real_v, bool ForRmin>
0375   VECGEOM_FORCE_INLINE VECCORE_ATT_HOST_DEVICE static Real_v ToBoundary(UnplacedStruct_t const &torus,
0376                                                                         Vector3D<Real_v> const &pt,
0377                                                                         Vector3D<Real_v> const &dir, Real_v radius,
0378                                                                         bool out)
0379   {
0380     // to be taken from ROOT
0381     // Returns distance to the surface or the torus from a point, along
0382     // a direction. Point is close enough to the boundary so that the distance
0383     // to the torus is decreasing while moving along the given direction.
0384
0385     // Compute coeficients of the quartic
0386     Real_v s                 = vecgeom::kInfLength;
0387     VECGEOM_CONST Real_v tol = 100. * vecgeom::kTolerance;
0388     Real_v r0sq              = pt[0] * pt[0] + pt[1] * pt[1] + pt[2] * pt[2];
0389     Real_v rdotn             = pt[0] * dir[0] + pt[1] * dir[1] + pt[2] * dir[2];
0390     Real_v rsumsq            = 1. + radius * radius;
0391     Real_v a                 = 4. * rdotn;
0392     Real_v b                 = 2. * (r0sq + 2. * rdotn * rdotn - rsumsq + 2. * dir[2] * dir[2]);
0393     Real_v c                 = 4. * (r0sq * rdotn - rsumsq * rdotn + 2. * pt[2] * dir[2]);
0394     Real_v d                 = r0sq * r0sq - 2. * r0sq * rsumsq + 4. * pt[2] * pt[2] +
0395                (1. - radius * radius) * (1. - radius * radius);
0396
0397     Real_v x[4] = {vecgeom::kInfLength, vecgeom::kInfLength, vecgeom::kInfLength, vecgeom::kInfLength};
0398     int nsol    = 0;
0399
0400     // special condition
0401     if (vecCore::MaskFull(Abs(dir[2]) < 1E-3 && Abs(pt[2]) < 0.1 * radius)) {
0402       Real_v r0        = 1. - Sqrt((radius - pt[2]) * (radius + pt[2]));
0403       Real_v invdirxy2 = 1. / (1 - dir.z() * dir.z());
0404       Real_v b0        = (pt[0] * dir[0] + pt[1] * dir[1]) * invdirxy2;
0405       Real_v c0        = (pt[0] * pt[0] + (pt[1] - r0) * (pt[1] + r0)) * invdirxy2;
0406       Real_v delta     = b0 * b0 - c0;
0407       if (vecCore::MaskFull(delta > 0)) {
0408         x[nsol] = -b0 - Sqrt(delta);
0409         if (vecCore::MaskFull(x[nsol] > -tol)) nsol++;
0410         x[nsol] = -b0 + Sqrt(delta);
0411         if (vecCore::MaskFull(x[nsol] > -tol)) nsol++;
0412       }
0413       r0    = 1. + Sqrt((radius - pt[2]) * (radius + pt[2]));
0414       c0    = (pt[0] * pt[0] + (pt[1] - r0) * (pt[1] + r0)) * invdirxy2;
0415       delta = b0 * b0 - c0;
0416       if (vecCore::MaskFull(delta > 0)) {
0417         x[nsol] = -b0 - Sqrt(delta);
0418         if (vecCore::MaskFull(x[nsol] > -tol)) nsol++;
0419         x[nsol] = -b0 + Sqrt(delta);
0420         if (vecCore::MaskFull(x[nsol] > -tol)) nsol++;
0421       }
0422       if (nsol) {
0423         Sort4(x);
0424       }
0425     } else { // generic case
0426       nsol = SolveQuartic(a, b, c, d, x);
0427     }
0428     if (!nsol) {
0429       return vecgeom::kInfLength;
0430     }
0431
0432     // look for first positive solution
0433     Real_v ndotd;
0434     bool inner = vecCore::MaskFull(Abs(radius - torus.rmin()/torus.rtor()) < vecgeom::kTolerance);
0435     for (int i = 0; i < nsol; i++) {
0436       if (vecCore::MaskFull(x[i] < -100)) continue;
0437
0438       Vector3D<Real_v> r0   = pt + x[i] * dir;
0439       Vector3D<Real_v> norm = r0;
0440       r0.z()                = 0.;
0441       r0.Normalize();
0442       //r0 *= torus.rtor();
0443       norm -= r0;
0444       // norm = pt
0445       // for (unsigned int ipt = 0; ipt < 3; ipt++)
0446       //   norm[ipt] = pt[ipt] + x[i] * dir[ipt] - r0[ipt];
0447       // ndotd = norm[0] * dir[0] + norm[1] * dir[1] + norm[2] * dir[2];
0448       ndotd = norm.Dot(dir);
0449       if (inner ^ out) {
0450         if (vecCore::MaskFull(ndotd < 0)) continue; // discard this solution
0451       } else {
0452         if (vecCore::MaskFull(ndotd > 0)) continue; // discard this solution
0453       }
0454
0455       // The crossing point should be in the phi wedge
0456       if (torus.dphi() < vecgeom::kTwoPi) {
0457         if (!vecCore::MaskFull(torus.GetWedge().ContainsWithBoundary<Real_v>(r0))) continue;
0458       }
0459
0460       s = x[i];
0461       // refine solution with Newton iterations
0462       Real_v eps   = vecgeom::kInfLength;
0463       Real_v delta = s * s * s * s + a * s * s * s + b * s * s + c * s + d;
0464       Real_v eps0  = -delta / (4. * s * s * s + 3. * a * s * s + 2. * b * s + c);
0465       int ntry     = 0;
0466       while (vecCore::MaskFull(Abs(eps) > vecgeom::kTolerance)) {
0467         if (vecCore::MaskFull(Abs(eps0) > 200)) break;
0468         s += eps0;
0469         if (vecCore::MaskFull(Abs(s + eps0) < vecgeom::kTolerance)) break;
0470         delta = s * s * s * s + a * s * s * s + b * s * s + c * s + d;
0471         eps   = -delta / (4. * s * s * s + 3. * a * s * s + 2. * b * s + c);
0472         if (vecCore::MaskFull(Abs(eps) >= Abs(eps0))) break;
0473         ntry++;
0474         // Avoid infinite recursion
0475         if (ntry > 100) break;
0476         eps0 = eps;
0477       }
0478       // discard this solution
0479       if (vecCore::MaskFull(s < -tol)) continue;
0480       return Max(Real_v(0.), s);
0481     }
0482     return vecgeom::kInfLength;
0483   }
0484
0485   template <typename Real_v>
0486   VECGEOM_FORCE_INLINE VECCORE_ATT_HOST_DEVICE static void SafetyToOut(UnplacedStruct_t const &torus,
0487                                                                        Vector3D<Real_v> const &point, Real_v &safety)
0488   {
0489     Real_v rxy = Sqrt(point[0] * point[0] + point[1] * point[1]);
0490     Real_v rad = Sqrt((rxy - torus.rtor()) * (rxy - torus.rtor()) + point[2] * point[2]);
0491     safety     = torus.rmax() - rad;
0492     if (torus.rmin()) {
0493       safety = Min(rad - torus.rmin(), torus.rmax() - rad);
0494     }
0495
0496     // TODO: extend implementation for phi sector case
0497     bool hasphi = (torus.dphi() < kTwoPi);
0498     if (hasphi) {
0499       Real_v safetyPhi = torus.GetWedge().SafetyToOut<Real_v>(point);
0500       safety           = Min(safetyPhi, safety);
0501     }
0502   }
0503
0504   template <typename Real_v>
0505   VECGEOM_FORCE_INLINE VECCORE_ATT_HOST_DEVICE static void Contains(UnplacedStruct_t const &torus,
0506                                                                     Vector3D<Real_v> const &point,
0507                                                                     typename vecCore::Mask_v<Real_v> &contains)
0508   {
0509     using Bool_v = vecCore::Mask_v<Real_v>;
0510     Bool_v unused, outside;
0511     TorusImplementation2::GenericKernelForContainsAndInside<Real_v, true, true>(torus, point, unused, outside);
0512     contains = !outside;
0513   }
0514
0515   template <typename Real_v, typename Inside_t>
0516   VECGEOM_FORCE_INLINE VECCORE_ATT_HOST_DEVICE static void Inside(UnplacedStruct_t const &torus,
0517                                                                   Vector3D<Real_v> const &point, Inside_t &inside)
0518   {
0519     TorusImplementation2::InsideKernel<Real_v, Inside_t>(torus, point, inside);
0520   }
0521
0522   template <typename Real_v>
0523   VECGEOM_FORCE_INLINE VECCORE_ATT_HOST_DEVICE static void DistanceToIn(UnplacedStruct_t const &torus,
0524                                                                         Vector3D<Real_v> const &point,
0525                                                                         Vector3D<Real_v> const &direction,
0526                                                                         Real_v const &stepMax, Real_v &distance)
0527   {
0528
0529     Vector3D<Real_v> localPoint     = point;
0530     Vector3D<Real_v> localDirection = direction;
0531
0532     using Bool_v   = vecCore::Mask_v<Real_v>;
0533     using Inside_v = vecCore::Index_v<Real_v>;
0534
0535     ////////First naive implementation
0536     distance = kInfLength;
0537
0538     // Check Bounding Cylinder first
0539     Bool_v inBounds;
0540     Bool_v done         = Bool_v(false);
0541     Inside_v inside     = Inside_v(EInside::kOutside);
0542     Real_v tubeDistance = kInfLength;
0543
0544 #ifndef VECGEOM_NO_SPECIALIZATION
0545     // call the tube functionality -- first of all we check whether we are inside
0546     // bounding volume
0547     TubeImplementation<TubeTypes::HollowTube>::Contains(torus.GetBoundingTube().GetStruct(), localPoint, inBounds);
0548
0549     // only need to do this check if all particles (in vector) are outside ( otherwise useless )
0550     TubeImplementation<TubeTypes::HollowTube>::DistanceToIn(torus.GetBoundingTube().GetStruct(), localPoint,
0551                                                             localDirection, stepMax, tubeDistance);
0552 #else
0553     // call the tube functionality -- first of all we check whether we are inside
0554     // bounding volume
0555     TubeImplementation<TubeTypes::UniversalTube>::Contains(torus.GetBoundingTube().GetStruct(), localPoint, inBounds);
0556
0557     // only need to do this check if all particles (in vector) are outside ( otherwise useless )
0558     // vecCore::Mask_v<Real_v> notInBounds { !inBounds };
0559     if (!inBounds) {
0560       TubeImplementation<TubeTypes::UniversalTube>::DistanceToIn(torus.GetBoundingTube().GetStruct(), localPoint,
0561                                                                  localDirection, stepMax, tubeDistance);
0562     } else {
0563       tubeDistance = 0.;
0564     }
0565
0566 #endif // VECGEOM_NO_SPECIALIZATION
0567     if (inBounds) {
0568       // Check points on the wrong side (inside torus)
0569       TorusImplementation2::InsideKernel<Real_v, Inside_v>(torus, point, inside);
0570       if (vecCore::MaskFull(inside == Inside_v(EInside::kInside))) {
0571         done     = Bool_v(true);
0572         distance = Real_v(-1.);
0573       }
0574     } else {
0575       done = Bool_v(vecCore::MaskFull(tubeDistance == kInfLength));
0576     }
0577
0578     if (vecCore::EarlyReturnAllowed()) {
0579       if (vecCore::MaskFull(done)) {
0580         return;
0581       }
0582     }
0583
0584     // Propagate the point to the bounding tube, as this will reduce the
0585     // coefficients of the quartic and improve precision of the solutions
0586     localPoint += tubeDistance * localDirection;
0587     localPoint /= torus.rtor();
0588     Bool_v hasphi = Bool_v(torus.dphi() < vecgeom::kTwoPi);
0589     if (vecCore::MaskFull(hasphi)) {
0590       Real_v d1, d2;
0591
0592       auto wedge = torus.GetWedge();
0593       // checking distance to phi wedges
0594       // NOTE: if the tube told me its hitting surface, this would be unnessecary
0595       wedge.DistanceToIn<Real_v>(localPoint, localDirection, d1, d2);
0596
0597       // check phi intersections if bounding tube intersection is due to phi in which case we are done
0598       if (vecCore::MaskFull(d1 != kInfLength)) {
0599         Real_v daxis = DistSqrToTorusR(localPoint, localDirection, d1);
0600         if (vecCore::MaskFull(daxis >= torus.rmin2()/torus.rtor()/torus.rtor() && daxis < torus.rmax2()/torus.rtor()/torus.rtor() && d1 > -kTolerance)) {
0601           distance = d1;
0602         }
0603       }
0604
0605       if (vecCore::MaskFull(d2 != kInfLength)) {
0606         Real_v daxis = DistSqrToTorusR(localPoint, localDirection, d2);
0607         if (vecCore::MaskFull(daxis >= torus.rmin2()/torus.rtor()/torus.rtor() && daxis < torus.rmax2()/torus.rtor()/torus.rtor() && d2 > -kTolerance)) {
0608           distance = Min(distance, d2);
0609         }
0610       }
0611     }
0612
0613     Real_v dd = ToBoundary<Real_v, false>(torus, localPoint, localDirection, torus.rmax()/torus.rtor(), false);
0614
0615     // in case of a phi opening we also need to check the Rmin surface
0616     if (torus.rmin() > 0.) {
0617       Real_v ddrmin = ToBoundary<Real_v, true>(torus, localPoint, localDirection, torus.rmin()/torus.rtor(), false);
0618       dd            = Min(dd, ddrmin);
0619     }
0620     distance = Min(distance, dd);
0621     distance *= torus.rtor();
0622     distance += tubeDistance;
0623     // This has to be added because distance can become > kInfLength due to
0624     // missing early returns in CUDA. This makes comparisons to kInfLength fail.
0625     if (vecCore::MaskFull(Abs(distance) > kInfLength)) distance = kInfLength;
0626
0627     return;
0628   }
0629
0630   template <typename Real_v>
0631   VECGEOM_FORCE_INLINE VECCORE_ATT_HOST_DEVICE static void SafetyToIn(UnplacedStruct_t const &torus,
0632                                                                       Vector3D<Real_v> const &point, Real_v &safety)
0633   {
0634
0635     Vector3D<Real_v> localPoint = point;
0636
0637     // implementation taken from TGeoTorus
0638     Real_v rxy = Sqrt(localPoint[0] * localPoint[0] + localPoint[1] * localPoint[1]);
0639     Real_v rad = Sqrt((rxy - torus.rtor()) * (rxy - torus.rtor()) + localPoint[2] * localPoint[2]);
0640     safety     = rad - torus.rmax();
0641     if (torus.rmin()) {
0642       safety = Max(torus.rmin() - rad, rad - torus.rmax());
0643     }
0644
0645     bool hasphi = (torus.dphi() < kTwoPi);
0646     if (hasphi && vecCore::MaskFull(rxy != 0.)) {
0647       Real_v safetyPhi = torus.GetWedge().SafetyToIn<Real_v>(localPoint);
0648       safety           = Max(safetyPhi, safety);
0649     }
0650   }
0651
0652 }; // end struct
0653 } // namespace VECGEOM_IMPL_NAMESPACE
0654 } // namespace vecgeom
0655
0656 #endif // VECGEOM_VOLUMES_KERNEL_TORUSIMPLEMENTATION2_H_