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