volumes/kernel/EllipsoidImplementation.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 Ellipsoid
0006 /// @file volumes/kernel/EllipsoidImplementation.h
0007 /// @author Evgueni Tcherniaev
0008
0009 #ifndef VECGEOM_VOLUMES_KERNEL_ELLIPSOIDIMPLEMENTATION_H_
0010 #define VECGEOM_VOLUMES_KERNEL_ELLIPSOIDIMPLEMENTATION_H_
0011
0012 #include "VecGeom/base/Vector3D.h"
0013 #include "VecGeom/volumes/EllipsoidStruct.h"
0014 #include "VecGeom/volumes/kernel/GenericKernels.h"
0015 #include <VecCore/VecCore>
0016
0017 #include <cstdio>
0018 #include <iomanip>
0019
0020 namespace vecgeom {
0021
0022 VECGEOM_DEVICE_FORWARD_DECLARE(struct EllipsoidImplementation;);
0023 VECGEOM_DEVICE_DECLARE_CONV(struct, EllipsoidImplementation);
0024
0025 inline namespace VECGEOM_IMPL_NAMESPACE {
0026
0027 class PlacedEllipsoid;
0028 template <typename T>
0029 struct EllipsoidStruct;
0030 class UnplacedEllipsoid;
0031
0032 struct EllipsoidImplementation {
0033
0034   using PlacedShape_t    = PlacedEllipsoid;
0035   using UnplacedStruct_t = EllipsoidStruct<Precision>;
0036   using UnplacedVolume_t = UnplacedEllipsoid;
0037
0038   template <typename Real_v, typename Bool_v>
0039   VECGEOM_FORCE_INLINE VECCORE_ATT_HOST_DEVICE static void Contains(UnplacedStruct_t const &ellipsoid,
0040                                                                     Vector3D<Real_v> const &point, Bool_v &inside)
0041   {
0042     Bool_v unused(false), outside(false);
0043     GenericKernelForContainsAndInside<Real_v, Bool_v, false>(ellipsoid, point, unused, outside);
0044     inside = !outside;
0045   }
0046
0047   // BIG QUESTION: DO WE WANT TO GIVE ALL 3 TEMPLATE PARAMETERS
0048   // -- OR -- DO WE WANT TO DEDUCE Bool_v, Index_t from Real_v???
0049   template <typename Real_v, typename Inside_t>
0050   VECGEOM_FORCE_INLINE VECCORE_ATT_HOST_DEVICE static void Inside(UnplacedStruct_t const &ellipsoid,
0051                                                                   Vector3D<Real_v> const &point, Inside_t &inside)
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>(ellipsoid, 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 &ellipsoid, 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() * ellipsoid.fSx;
0077     Real_v y      = point.y() * ellipsoid.fSy;
0078     Real_v z      = point.z() * ellipsoid.fSz;
0079     Real_v distZ  = vecCore::math::Abs(z - ellipsoid.fScZMidCut) - ellipsoid.fScZDimCut;
0080     Real_v distR  = ellipsoid.fQ1 * (x * x + y * y + z * z) - ellipsoid.fQ2;
0081     Real_v safety = vecCore::math::Max(distZ, distR);
0082
0083     completelyoutside = safety > kHalfTolerance;
0084     if (ForInside) completelyinside = safety <= -kHalfTolerance;
0085     return;
0086   }
0087
0088   template <typename Real_v>
0089   VECGEOM_FORCE_INLINE VECCORE_ATT_HOST_DEVICE static void DistanceToIn(UnplacedStruct_t const &ellipsoid,
0090                                                                         Vector3D<Real_v> const &point,
0091                                                                         Vector3D<Real_v> const &direction,
0092                                                                         Real_v const & /*stepMax*/, Real_v &distance)
0093   {
0094     /* TODO :  Logic to calculate Distance from outside point to the Ellipsoid surface */
0095     using Bool_v = vecCore::Mask_v<Real_v>;
0096     distance     = kInfLength;
0097     Real_v offset(0.);
0098     Vector3D<Real_v> pcur(point);
0099
0100     // Move point closer, if required
0101     Real_v Rfar2(1024. * ellipsoid.fRsph * ellipsoid.fRsph); // 1024 = 32 * 32
0102     vecCore__MaskedAssignFunc(pcur, ((pcur.Mag2() > Rfar2) && (direction.Dot(point) < Real_v(0.))),
0103                               pcur + (offset = pcur.Mag() - Real_v(2.) * ellipsoid.fRsph) * direction);
0104
0105     // Scale ellipsoid to sphere
0106     Real_v px = pcur.x() * ellipsoid.fSx;
0107     Real_v py = pcur.y() * ellipsoid.fSy;
0108     Real_v pz = pcur.z() * ellipsoid.fSz;
0109     Real_v vx = direction.x() * ellipsoid.fSx;
0110     Real_v vy = direction.y() * ellipsoid.fSy;
0111     Real_v vz = direction.z() * ellipsoid.fSz;
0112
0113     // Check if point is leaving the solid
0114     Real_v pzcut = pz - ellipsoid.fScZMidCut;
0115     Real_v dzcut = Real_v(ellipsoid.fScZDimCut);
0116     Real_v distZ = vecCore::math::Abs(pzcut) - dzcut;
0117
0118     Real_v rr    = px * px + py * py + pz * pz;
0119     Real_v vv    = vx * vx + vy * vy + vz * vz;
0120     Real_v pv    = px * vx + py * vy + pz * vz;
0121     Real_v distR = ellipsoid.fQ1 * rr - ellipsoid.fQ2;
0122     Bool_v leaving =
0123         (distZ >= -kHalfTolerance && pzcut * vz >= Real_v(0.)) || (distR >= -kHalfTolerance && pv >= Real_v(0.));
0124
0125     // Find intersection with Z planes
0126     Real_v invz  = Real_v(-1.) / NonZero(vz);
0127     Real_v dz    = vecCore::math::CopySign(dzcut, invz);
0128     Real_v tzmin = (pzcut - dz) * invz;
0129     Real_v tzmax = (pzcut + dz) * invz;
0130
0131     // Find intersection with sphere
0132     Real_v A     = vv;
0133     Real_v B     = pv;
0134     Real_v C     = (rr - ellipsoid.fR * ellipsoid.fR);
0135     Real_v D     = B * B - A * C;
0136     Real_v sqrtD = vecCore::math::Sqrt(vecCore::math::Abs(D));
0137     Real_v trmin = (-B - sqrtD) / A;
0138     Real_v trmax = (-B + sqrtD) / A;
0139
0140     // Set preliminary distances to in/out
0141     Real_v tmin = vecCore::math::Max(tzmin, trmin);
0142     Real_v tmax = vecCore::math::Min(tzmax, trmax);
0143
0144     // Check if no intersection
0145     Real_v EPS  = Real_v(2.) * rr * vv * kEpsilon;
0146     Bool_v done = leaving || (D <= EPS) || ((tmax - tmin) <= kHalfTolerance);
0147
0148     // Set distance
0149     vecCore__MaskedAssignFunc(distance, !done, tmin + offset);
0150   }
0151
0152   template <typename Real_v>
0153   VECGEOM_FORCE_INLINE VECCORE_ATT_HOST_DEVICE static void DistanceToOut(UnplacedStruct_t const &ellipsoid,
0154                                                                          Vector3D<Real_v> const &point,
0155                                                                          Vector3D<Real_v> const &direction,
0156                                                                          Real_v const & /* stepMax */, Real_v &distance)
0157   {
0158     /* TODO :  Logic to calculate Distance from inside point to the Ellipsoid surface */
0159     using Bool_v = vecCore::Mask_v<Real_v>;
0160
0161     // Scale ellipsoid to sphere
0162     Real_v px = point.x() * ellipsoid.fSx;
0163     Real_v py = point.y() * ellipsoid.fSy;
0164     Real_v pz = point.z() * ellipsoid.fSz;
0165     Real_v vx = direction.x() * ellipsoid.fSx;
0166     Real_v vy = direction.y() * ellipsoid.fSy;
0167     Real_v vz = direction.z() * ellipsoid.fSz;
0168
0169     // Check if point is outside ("wrong side")
0170     Real_v pzcut = pz - ellipsoid.fScZMidCut;
0171     Real_v dzcut = Real_v(ellipsoid.fScZDimCut);
0172     Real_v distZ = vecCore::math::Abs(pzcut) - dzcut;
0173
0174     Real_v rr      = px * px + py * py + pz * pz;
0175     Real_v vv      = vx * vx + vy * vy + vz * vz;
0176     Real_v pv      = px * vx + py * vy + pz * vz;
0177     Real_v distR   = ellipsoid.fQ1 * rr - ellipsoid.fQ2;
0178     Bool_v outside = vecCore::math::Max(distR, distZ) > kHalfTolerance;
0179
0180     distance = Real_v(0.);
0181     vecCore__MaskedAssignFunc(distance, outside, Real_v(-1.));
0182
0183     // Find intersection with Z planes
0184     Real_v tzmax = kMaximum;
0185     vecCore__MaskedAssignFunc(tzmax, vz != Real_v(0.), (vecCore::math::CopySign(Real_v(dzcut), vz) - pzcut) / vz);
0186
0187     // Find intersection with sphere
0188     Real_v B     = pv / vv;
0189     Real_v C     = (rr - ellipsoid.fR * ellipsoid.fR) / vv;
0190     Real_v D     = B * B - C;
0191     Real_v sqrtD = vecCore::math::Sqrt(vecCore::math::Abs(D));
0192     Real_v trmax = -B + sqrtD;
0193
0194     // Check if no intersection
0195     Real_v EPS  = Real_v(2.) * rr * vv * kEpsilon;
0196     Bool_v done = outside || (D <= EPS);
0197
0198     // Set distance
0199     vecCore__MaskedAssignFunc(distance, !done, vecCore::math::Min(tzmax, trmax));
0200   }
0201
0202   template <typename Real_v>
0203   VECGEOM_FORCE_INLINE VECCORE_ATT_HOST_DEVICE static void SafetyToIn(UnplacedStruct_t const &ellipsoid,
0204                                                                       Vector3D<Real_v> const &point, Real_v &safety)
0205   {
0206     /* TODO :  Logic to calculate Safety from outside point to the Ellipsoid surface */
0207     Real_v x = point.x() * ellipsoid.fSx;
0208     Real_v y = point.y() * ellipsoid.fSy;
0209     Real_v z = point.z() * ellipsoid.fSz;
0210     Real_v r = vecCore::math::Sqrt(x * x + y * y + z * z);
0211     // Set safety to zero if point is on surface
0212     Real_v safeZ = vecCore::math::Abs(z - ellipsoid.fScZMidCut) - ellipsoid.fScZDimCut;
0213     Real_v safeR = r - ellipsoid.fR;
0214     safety       = vecCore::math::Max(safeZ, safeR);
0215     vecCore::MaskedAssign(safety, vecCore::math::Abs(safety) <= kHalfTolerance, Real_v(0.));
0216     // Adjust safety using bounding box
0217     Real_v distZ  = vecCore::math::Max(point.z() - ellipsoid.fZTopCut, ellipsoid.fZBottomCut - point.z());
0218     Real_v distXY = vecCore::math::Max(vecCore::math::Abs(point.x()) - ellipsoid.fXmax,
0219                                        vecCore::math::Abs(point.y()) - ellipsoid.fYmax);
0220     vecCore__MaskedAssignFunc(safety, safety > Real_v(0.), vecCore::math::Max(safety, distZ, distXY));
0221   }
0222
0223   template <typename Real_v>
0224   VECGEOM_FORCE_INLINE VECCORE_ATT_HOST_DEVICE static void SafetyToOut(UnplacedStruct_t const &ellipsoid,
0225                                                                        Vector3D<Real_v> const &point, Real_v &safety)
0226   {
0227     /* TODO :  Logic to calculate Safety from inside point to the Ellipsoid surface */
0228     Real_v x = point.x() * ellipsoid.fSx;
0229     Real_v y = point.y() * ellipsoid.fSy;
0230     Real_v z = point.z() * ellipsoid.fSz;
0231     // Set safety to zero if point is on surface
0232     Real_v safeR = ellipsoid.fR - vecCore::math::Sqrt(x * x + y * y + z * z);
0233     Real_v safeZ = ellipsoid.fScZDimCut - vecCore::math::Abs(z - ellipsoid.fScZMidCut);
0234     safety       = vecCore::math::Min(safeZ, safeR);
0235     vecCore::MaskedAssign(safety, vecCore::math::Abs(safety) <= kHalfTolerance, Real_v(0.));
0236     // Adjust safety in z direction
0237     Real_v distZ = vecCore::math::Min(ellipsoid.fZTopCut - point.z(), point.z() - ellipsoid.fZBottomCut);
0238     vecCore__MaskedAssignFunc(safety, safety > Real_v(0.), vecCore::math::Min(safeR, distZ));
0239   }
0240
0241   template <typename Real_v>
0242   VECGEOM_FORCE_INLINE VECCORE_ATT_HOST_DEVICE static Vector3D<Real_v> NormalKernel(
0243       UnplacedStruct_t const &ellipsoid, Vector3D<Real_v> const &point, typename vecCore::Mask_v<Real_v> &valid)
0244   {
0245     // Computes the normal on a surface and returns it as a unit vector
0246     //   In case if the point is further than kHalfTolerance from the surface, set valid=false
0247     //   Must return a valid vector (even if the point is not on the surface)
0248     //
0249     //   On an edge provide an average normal of the corresponding base and lateral surface
0250     Vector3D<Real_v> normal(0.);
0251     valid = true;
0252
0253     Real_v px   = point.x();
0254     Real_v py   = point.y();
0255     Real_v pz   = point.z();
0256     Real_v A    = ellipsoid.fDx;
0257     Real_v B    = ellipsoid.fDy;
0258     Real_v C    = ellipsoid.fDz;
0259     Real_v x    = px * ellipsoid.fSx;
0260     Real_v y    = py * ellipsoid.fSy;
0261     Real_v z    = pz * ellipsoid.fSz;
0262     Real_v mag2 = x * x + y * y + z * z;
0263
0264     // Check lateral surface
0265     Real_v distR = ellipsoid.fQ1 * mag2 - ellipsoid.fQ2;
0266     vecCore__MaskedAssignFunc(normal, vecCore::math::Abs(distR) <= kHalfTolerance,
0267                               Vector3D<Real_v>(px / (A * A), py / (B * B), pz / (C * C)).Unit());
0268
0269     // Check z cuts
0270     Real_v distZ = vecCore::math::Abs(z - ellipsoid.fScZMidCut) - ellipsoid.fScZDimCut;
0271     vecCore__MaskedAssignFunc(normal[2], vecCore::math::Abs(distZ) <= kHalfTolerance,
0272                               normal[2] + vecCore::math::Sign(z - ellipsoid.fScZMidCut));
0273
0274     // Average normal, if required
0275     vecCore__MaskedAssignFunc(normal, normal.Mag2() > 1., normal.Unit());
0276     vecCore::Mask_v<Real_v> done = normal.Mag2() > Real_v(0.);
0277     if (vecCore::MaskFull(done)) return normal;
0278
0279     // Point is not on the surface - normally, this should never be
0280     // Return normal to the nearest surface
0281     vecCore__MaskedAssignFunc(valid, !done, false);
0282     vecCore__MaskedAssignFunc(normal[2], !done, vecCore::math::Sign(z - ellipsoid.fScZMidCut));
0283     vecCore__MaskedAssignFunc(distR, !done, vecCore::math::Sqrt(mag2) - ellipsoid.fR);
0284     vecCore__MaskedAssignFunc(normal, !done && distR > distZ && mag2 > Real_v(0.),
0285                               Vector3D<Real_v>(px / (A * A), py / (B * B), pz / (C * C)).Unit());
0286     return normal;
0287   }
0288 };
0289 } // namespace VECGEOM_IMPL_NAMESPACE
0290 } // namespace vecgeom
0291
0292 #endif // VECGEOM_VOLUMES_KERNEL_ELLIPSOIDIMPLEMENTATION_H_