volumes/kernel/EllipticalTubeImplementation.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 EllipticalTube
0006 /// @file volumes/kernel/EllipticalTubeImplementation.h
0007 /// @author Raman Sehgal, Evgueni Tcherniaev
0008
0009 #ifndef VECGEOM_VOLUMES_KERNEL_ELLIPTICALTUBEIMPLEMENTATION_H_
0010 #define VECGEOM_VOLUMES_KERNEL_ELLIPTICALTUBEIMPLEMENTATION_H_
0011
0012 #include "VecGeom/base/Vector3D.h"
0013 #include "VecGeom/volumes/EllipticalTubeStruct.h"
0014 #include "VecGeom/volumes/kernel/GenericKernels.h"
0015 #include <VecCore/VecCore>
0016
0017 #include <cstdio>
0018
0019 namespace vecgeom {
0020
0021 VECGEOM_DEVICE_FORWARD_DECLARE(struct EllipticalTubeImplementation;);
0022 VECGEOM_DEVICE_DECLARE_CONV(struct, EllipticalTubeImplementation);
0023
0024 inline namespace VECGEOM_IMPL_NAMESPACE {
0025
0026 class PlacedEllipticalTube;
0027 template <typename T>
0028 struct EllipticalTubeStruct;
0029 class UnplacedEllipticalTube;
0030
0031 struct EllipticalTubeImplementation {
0032
0033   using PlacedShape_t    = PlacedEllipticalTube;
0034   using UnplacedStruct_t = EllipticalTubeStruct<Precision>;
0035   using UnplacedVolume_t = UnplacedEllipticalTube;
0036
0037   template <typename Real_v, typename Bool_v>
0038   VECGEOM_FORCE_INLINE VECCORE_ATT_HOST_DEVICE static void Contains(UnplacedStruct_t const &ellipticaltube,
0039                                                                     Vector3D<Real_v> const &point, Bool_v &inside)
0040   {
0041     Bool_v unused(false), outside(false);
0042     GenericKernelForContainsAndInside<Real_v, Bool_v, false>(ellipticaltube, point, unused, outside);
0043     inside = !outside;
0044   }
0045
0046   // BIG QUESTION: DO WE WANT TO GIVE ALL 3 TEMPLATE PARAMETERS
0047   // -- OR -- DO WE WANT TO DEDUCE Bool_v, Index_t from Real_v???
0048   template <typename Real_v, typename Inside_t>
0049   VECGEOM_FORCE_INLINE VECCORE_ATT_HOST_DEVICE static void Inside(UnplacedStruct_t const &ellipticaltube,
0050                                                                   Vector3D<Real_v> const &point, Inside_t &inside)
0051   {
0052
0053     using Bool_v       = vecCore::Mask_v<Real_v>;
0054     using InsideBool_v = vecCore::Mask_v<Inside_t>;
0055     Bool_v completelyinside, completelyoutside;
0056     GenericKernelForContainsAndInside<Real_v, Bool_v, true>(ellipticaltube, point, completelyinside, completelyoutside);
0057     inside = EInside::kSurface;
0058     vecCore::MaskedAssign(inside, (InsideBool_v)completelyoutside, Inside_t(EInside::kOutside));
0059     vecCore::MaskedAssign(inside, (InsideBool_v)completelyinside, Inside_t(EInside::kInside));
0060   }
0061
0062   template <typename Real_v, typename Bool_v, bool ForInside>
0063   VECGEOM_FORCE_INLINE VECCORE_ATT_HOST_DEVICE static void GenericKernelForContainsAndInside(
0064       UnplacedStruct_t const &ellipticaltube, Vector3D<Real_v> const &point, Bool_v &completelyinside,
0065       Bool_v &completelyoutside)
0066   {
0067     /* TODO : Logic to check where the point is inside or not.
0068     **
0069     ** if ForInside is false then it will only check if the point is outside,
0070     ** and is used by Contains function
0071     **
0072     ** if ForInside is true then it will check whether the point is inside or outside,
0073     ** and if neither inside nor outside then it is on the surface.
0074     ** and is used by Inside function
0075     */
0076     Real_v x      = point.x() * ellipticaltube.fSx;
0077     Real_v y      = point.y() * ellipticaltube.fSy;
0078     Real_v distR  = ellipticaltube.fQ1 * (x * x + y * y) - ellipticaltube.fQ2;
0079     Real_v distZ  = vecCore::math::Abs(point.z()) - ellipticaltube.fDz;
0080     Real_v safety = vecCore::math::Max(distR, distZ);
0081
0082     completelyoutside = safety > kHalfTolerance;
0083     if (ForInside) completelyinside = safety <= -kHalfTolerance;
0084     return;
0085   }
0086
0087   template <typename Real_v>
0088   VECGEOM_FORCE_INLINE VECCORE_ATT_HOST_DEVICE static void DistanceToIn(UnplacedStruct_t const &ellipticaltube,
0089                                                                         Vector3D<Real_v> const &point,
0090                                                                         Vector3D<Real_v> const &direction,
0091                                                                         Real_v const & /*stepMax*/, Real_v &distance)
0092   {
0093     /* TODO :  Logic to calculate Distance from outside point to the EllipticalTube surface */
0094     using Bool_v = vecCore::Mask_v<Real_v>;
0095     distance     = kInfLength;
0096     Real_v offset(0.);
0097     Vector3D<Real_v> pcur(point);
0098
0099     // Move point closer, if required
0100     Real_v Rfar2(1024. * ellipticaltube.fRsph * ellipticaltube.fRsph); // 1024 = 32 * 32
0101     vecCore__MaskedAssignFunc(pcur, ((pcur.Mag2() > Rfar2) && (direction.Dot(point) < Real_v(0.))),
0102                               pcur + (offset = pcur.Mag() - Real_v(2.) * ellipticaltube.fRsph) * direction);
0103
0104     // Scale elliptical tube to cylinder
0105     Real_v px = pcur.x() * ellipticaltube.fSx;
0106     Real_v py = pcur.y() * ellipticaltube.fSy;
0107     Real_v pz = pcur.z();
0108     Real_v vx = direction.x() * ellipticaltube.fSx;
0109     Real_v vy = direction.y() * ellipticaltube.fSy;
0110     Real_v vz = direction.z();
0111
0112     // Find intersection with Z planes
0113     Real_v invz  = Real_v(-1.) / NonZero(vz);
0114     Real_v dz    = vecCore::math::CopySign(Real_v(ellipticaltube.fDz), invz);
0115     Real_v tzmin = (pz - dz) * invz;
0116     Real_v tzmax = (pz + dz) * invz;
0117
0118     // Find intersection with lateral surface, solve equation: A t^2 + 2B t + C = 0
0119     Real_v rr = px * px + py * py;
0120     Real_v A  = vx * vx + vy * vy;
0121     Real_v B  = px * vx + py * vy;
0122     Real_v C  = rr - ellipticaltube.fR * ellipticaltube.fR;
0123     Real_v D  = B * B - A * C;
0124
0125     // Check if point leaving shape
0126     Real_v distZ       = vecCore::math::Abs(pz) - ellipticaltube.fDz;
0127     Real_v distR       = ellipticaltube.fQ1 * rr - ellipticaltube.fQ2;
0128     Bool_v parallelToZ = (A < kEpsilon || vecCore::math::Abs(vz) >= Real_v(1.));
0129     Bool_v leaving     = (distZ >= -kHalfTolerance && pz * vz >= Real_v(0.)) ||
0130                      (distR >= -kHalfTolerance && (B >= Real_v(0.) || parallelToZ));
0131
0132     // Two special cases where D <= 0:
0133     //   1) trajectory parallel to Z axis (A = 0, B = 0, C - any, D = 0)
0134     //   2) touch (D = 0) or no intersection (D < 0) with lateral surface
0135     vecCore__MaskedAssignFunc(distance, !leaving && parallelToZ, tzmin + offset);   // 1)
0136     Bool_v done = (leaving || parallelToZ || D <= A * A * ellipticaltube.fScratch); // 2)
0137
0138     // if (D <= A * A * ellipticaltube.fScratch) std::cerr << "=== SCRATCH D = " << D << std::endl;
0139
0140     // Find roots of the quadratic
0141     Real_v tmp(0.), t1(0.), t2(0.);
0142     vecCore__MaskedAssignFunc(tmp, !done, -B - vecCore::math::CopySign(vecCore::math::Sqrt(D), B));
0143     vecCore__MaskedAssignFunc(t1, !done, tmp / A);
0144     vecCore__MaskedAssignFunc(t2, !done, C / tmp);
0145     Real_v trmin = vecCore::math::Min(t1, t2);
0146     Real_v trmax = vecCore::math::Max(t1, t2);
0147
0148     // Set distance
0149     // No special check for inside points, for inside points distance will be negative
0150     Real_v tin  = vecCore::math::Max(tzmin, trmin);
0151     Real_v tout = vecCore::math::Min(tzmax, trmax);
0152     vecCore__MaskedAssignFunc(distance, !done && (tout - tin) >= kHalfTolerance, tin + offset);
0153   }
0154
0155   template <typename Real_v>
0156   VECGEOM_FORCE_INLINE VECCORE_ATT_HOST_DEVICE static void DistanceToOut(UnplacedStruct_t const &ellipticaltube,
0157                                                                          Vector3D<Real_v> const &point,
0158                                                                          Vector3D<Real_v> const &direction,
0159                                                                          Real_v const & /* stepMax */, Real_v &distance)
0160   {
0161     /* TODO :  Logic to calculate Distance from inside point to the EllipticalTube surface */
0162     using Bool_v = vecCore::Mask_v<Real_v>;
0163
0164     // Scale elliptical tube to cylinder
0165     Real_v px = point.x() * ellipticaltube.fSx;
0166     Real_v py = point.y() * ellipticaltube.fSy;
0167     Real_v pz = point.z();
0168     Real_v vx = direction.x() * ellipticaltube.fSx;
0169     Real_v vy = direction.y() * ellipticaltube.fSy;
0170     Real_v vz = direction.z();
0171
0172     // Check if point is outside ("wrong side")
0173     Real_v rr      = px * px + py * py;
0174     Real_v distR   = ellipticaltube.fQ1 * rr - ellipticaltube.fQ2;
0175     Real_v distZ   = vecCore::math::Abs(pz) - ellipticaltube.fDz;
0176     Bool_v outside = vecCore::math::Max(distR, distZ) > kHalfTolerance;
0177     distance       = Real_v(0.);
0178     vecCore__MaskedAssignFunc(distance, outside, Real_v(-1.));
0179
0180     // Find intersection with Z planes
0181     Real_v tzmax = kMaximum;
0182     vecCore__MaskedAssignFunc(tzmax, vz != Real_v(0.),
0183                               (vecCore::math::CopySign(Real_v(ellipticaltube.fDz), vz) - pz) / vz);
0184
0185     // Find intersection with lateral surface, solve equation: A t^2 + 2B t + C = 0
0186     Real_v A = vx * vx + vy * vy;
0187     Real_v B = px * vx + py * vy;
0188     Real_v C = rr - ellipticaltube.fR * ellipticaltube.fR;
0189     Real_v D = B * B - A * C;
0190
0191     // Two cases where D <= 0:
0192     //   1) trajectory parallel to Z axis (A = 0, B = 0, C - any, D = 0)
0193     //   2) touch (D = 0) or no intersection (D < 0) with lateral surface
0194     Bool_v parallelToZ = (A < kEpsilon || vecCore::math::Abs(vz) >= Real_v(1.));
0195     vecCore__MaskedAssignFunc(distance, (!outside && parallelToZ), tzmax); // 1)
0196     Bool_v done = (outside || parallelToZ || D <= Real_v(0.));             // 2)
0197     // Bool_v done = (outside || parallelToZ || D < A * A * ellipticaltube.fScratch); // alternative 2)
0198
0199     // Set distance
0200     vecCore__MaskedAssignFunc(distance, !done && B >= Real_v(0.),
0201                               vecCore::math::Min(tzmax, -C / (vecCore::math::Sqrt(D) + B)));
0202     vecCore__MaskedAssignFunc(distance, !done && B < Real_v(0.),
0203                               vecCore::math::Min(tzmax, (vecCore::math::Sqrt(D) - B) / A));
0204   }
0205
0206   template <typename Real_v>
0207   VECGEOM_FORCE_INLINE VECCORE_ATT_HOST_DEVICE static void SafetyToIn(UnplacedStruct_t const &ellipticaltube,
0208                                                                       Vector3D<Real_v> const &point, Real_v &safety)
0209   {
0210     /* TODO :  Logic to calculate Safety from outside point to the EllipticalTube surface */
0211     Real_v x     = point.x() * ellipticaltube.fSx;
0212     Real_v y     = point.y() * ellipticaltube.fSy;
0213     Real_v distR = vecCore::math::Sqrt(x * x + y * y) - ellipticaltube.fR;
0214     Real_v distZ = vecCore::math::Abs(point.z()) - ellipticaltube.fDz;
0215
0216     safety = vecCore::math::Max(distR, distZ);
0217     vecCore::MaskedAssign(safety, vecCore::math::Abs(safety) <= kHalfTolerance, Real_v(0.));
0218   }
0219
0220   template <typename Real_v>
0221   VECGEOM_FORCE_INLINE VECCORE_ATT_HOST_DEVICE static void SafetyToOut(UnplacedStruct_t const &ellipticaltube,
0222                                                                        Vector3D<Real_v> const &point, Real_v &safety)
0223   {
0224     /* TODO :  Logic to calculate Safety from inside point to the EllipticalTube surface */
0225     Real_v x     = point.x() * ellipticaltube.fSx;
0226     Real_v y     = point.y() * ellipticaltube.fSy;
0227     Real_v distR = ellipticaltube.fR - vecCore::math::Sqrt(x * x + y * y);
0228     Real_v distZ = ellipticaltube.fDz - vecCore::math::Abs(point.z());
0229
0230     safety = vecCore::math::Min(distR, distZ);
0231     vecCore::MaskedAssign(safety, vecCore::math::Abs(safety) <= kHalfTolerance, Real_v(0.));
0232   }
0233
0234   template <typename Real_v>
0235   VECGEOM_FORCE_INLINE VECCORE_ATT_HOST_DEVICE static Vector3D<Real_v> NormalKernel(
0236       UnplacedStruct_t const &ellipticaltube, Vector3D<Real_v> const &point, typename vecCore::Mask_v<Real_v> &valid)
0237   {
0238     // Computes the normal on a surface and returns it as a unit vector
0239     //   In case if the point is further than kHalfTolerance from the surface, set valid=false
0240     //   Must return a valid vector (even if the point is not on the surface)
0241     //
0242     //   On an edge provide an average normal of the corresponding base and lateral surface
0243     Vector3D<Real_v> normal(0.);
0244     valid = true;
0245
0246     Real_v x     = point.x() * ellipticaltube.fSx;
0247     Real_v y     = point.y() * ellipticaltube.fSy;
0248     Real_v distR = ellipticaltube.fQ1 * (x * x + y * y) - ellipticaltube.fQ2;
0249     vecCore__MaskedAssignFunc(
0250         normal, vecCore::math::Abs(distR) <= kHalfTolerance,
0251         Vector3D<Real_v>(point.x() * ellipticaltube.fDDy, point.y() * ellipticaltube.fDDx, 0.).Unit());
0252
0253     Real_v distZ = vecCore::math::Abs(point.z()) - ellipticaltube.fDz;
0254     vecCore__MaskedAssignFunc(normal[2], vecCore::math::Abs(distZ) <= kHalfTolerance, vecCore::math::Sign(point[2]));
0255     vecCore__MaskedAssignFunc(normal, normal.Mag2() > 1., normal.Unit());
0256
0257     vecCore::Mask_v<Real_v> done = normal.Mag2() > Real_v(0.);
0258     if (vecCore::MaskFull(done)) return normal;
0259
0260     // Point is not on the surface - normally, this should never be
0261     // Return normal to the nearest surface
0262     vecCore__MaskedAssignFunc(valid, !done, false);
0263     vecCore__MaskedAssignFunc(normal[2], !done, vecCore::math::Sign(point[2]));
0264     vecCore__MaskedAssignFunc(distR, !done, vecCore::math::Sqrt(x * x + y * y) - ellipticaltube.fR);
0265     vecCore__MaskedAssignFunc(
0266         normal, !done && distR > distZ && (x * x + y * y) > Real_v(0.),
0267         Vector3D<Real_v>(point.x() * ellipticaltube.fDDy, point.y() * ellipticaltube.fDDx, 0.).Unit());
0268     return normal;
0269   }
0270 };
0271 } // namespace VECGEOM_IMPL_NAMESPACE
0272 } // namespace vecgeom
0273
0274 #endif // VECGEOM_VOLUMES_KERNEL_ELLIPTICALTUBEIMPLEMENTATION_H_