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