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   VECCORE_ATT_HOST_DEVICE
0039   static void PrintType()
0040   {
0041     //  printf("SpecializedEllipsoid<%i, %i>", transCodeT, rotCodeT);
0042   }
0043
0044   template <typename Stream>
0045   static void PrintType(Stream &st, int transCodeT = translation::kGeneric, int rotCodeT = rotation::kGeneric)
0046   {
0047     st << "SpecializedEllipsoid<" << transCodeT << "," << rotCodeT << ">";
0048   }
0049
0050   template <typename Stream>
0051   static void PrintImplementationType(Stream &st)
0052   {
0053     (void)st;
0054     // st << "EllipsoidImplementation<" << transCodeT << "," << rotCodeT << ">";
0055   }
0056
0057   template <typename Stream>
0058   static void PrintUnplacedType(Stream &st)
0059   {
0060     (void)st;
0061     // TODO: this is wrong
0062     // st << "UnplacedEllipsoid";
0063   }
0064
0065   template <typename Real_v, typename Bool_v>
0066   VECGEOM_FORCE_INLINE
0067   VECCORE_ATT_HOST_DEVICE
0068   static void Contains(UnplacedStruct_t const &ellipsoid, Vector3D<Real_v> const &point, Bool_v &inside)
0069   {
0070     Bool_v unused, outside;
0071     GenericKernelForContainsAndInside<Real_v, Bool_v, false>(ellipsoid, point, unused, outside);
0072     inside = !outside;
0073   }
0074
0075   // BIG QUESTION: DO WE WANT TO GIVE ALL 3 TEMPLATE PARAMETERS
0076   // -- OR -- DO WE WANT TO DEDUCE Bool_v, Index_t from Real_v???
0077   template <typename Real_v, typename Inside_t>
0078   VECGEOM_FORCE_INLINE
0079   VECCORE_ATT_HOST_DEVICE
0080   static void Inside(UnplacedStruct_t const &ellipsoid, Vector3D<Real_v> const &point, Inside_t &inside)
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>(ellipsoid, 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 &ellipsoid, 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() * ellipsoid.fSx;
0107     Real_v y      = point.y() * ellipsoid.fSy;
0108     Real_v z      = point.z() * ellipsoid.fSz;
0109     Real_v distZ  = vecCore::math::Abs(z - ellipsoid.fScZMidCut) - ellipsoid.fScZDimCut;
0110     Real_v distR  = ellipsoid.fQ1 * (x * x + y * y + z * z) - ellipsoid.fQ2;
0111     Real_v safety = vecCore::math::Max(distZ, distR);
0112
0113     completelyoutside = safety > kHalfTolerance;
0114     if (ForInside) completelyinside = safety <= -kHalfTolerance;
0115     return;
0116   }
0117
0118   template <typename Real_v>
0119   VECGEOM_FORCE_INLINE
0120   VECCORE_ATT_HOST_DEVICE
0121   static void DistanceToIn(UnplacedStruct_t const &ellipsoid, Vector3D<Real_v> const &point,
0122                            Vector3D<Real_v> const &direction, Real_v const & /*stepMax*/, Real_v &distance)
0123   {
0124     /* TODO :  Logic to calculate Distance from outside point to the Ellipsoid surface */
0125     using Bool_v = vecCore::Mask_v<Real_v>;
0126     distance     = kInfLength;
0127     Real_v offset(0.);
0128     Vector3D<Real_v> pcur(point);
0129
0130     // Move point closer, if required
0131     Real_v Rfar2(1024. * ellipsoid.fRsph * ellipsoid.fRsph); // 1024 = 32 * 32
0132     vecCore__MaskedAssignFunc(pcur, ((pcur.Mag2() > Rfar2) && (direction.Dot(point) < Real_v(0.))),
0133                               pcur + (offset = pcur.Mag() - Real_v(2.) * ellipsoid.fRsph) * direction);
0134
0135     // Scale ellipsoid to sphere
0136     Real_v px = pcur.x() * ellipsoid.fSx;
0137     Real_v py = pcur.y() * ellipsoid.fSy;
0138     Real_v pz = pcur.z() * ellipsoid.fSz;
0139     Real_v vx = direction.x() * ellipsoid.fSx;
0140     Real_v vy = direction.y() * ellipsoid.fSy;
0141     Real_v vz = direction.z() * ellipsoid.fSz;
0142
0143     // Check if point is leaving the solid
0144     Real_v pzcut = pz - ellipsoid.fScZMidCut;
0145     Real_v dzcut = Real_v(ellipsoid.fScZDimCut);
0146     Real_v distZ = vecCore::math::Abs(pzcut) - dzcut;
0147
0148     Real_v rr    = px * px + py * py + pz * pz;
0149     Real_v vv    = vx * vx + vy * vy + vz * vz;
0150     Real_v pv    = px * vx + py * vy + pz * vz;
0151     Real_v distR = ellipsoid.fQ1 * rr - ellipsoid.fQ2;
0152     Bool_v leaving =
0153         (distZ >= -kHalfTolerance && pzcut * vz >= Real_v(0.)) || (distR >= -kHalfTolerance && pv >= Real_v(0.));
0154
0155     // Find intersection with Z planes
0156     Real_v invz  = Real_v(-1.) / NonZero(vz);
0157     Real_v dz    = vecCore::math::CopySign(dzcut, invz);
0158     Real_v tzmin = (pzcut - dz) * invz;
0159     Real_v tzmax = (pzcut + dz) * invz;
0160
0161     // Find intersection with sphere
0162     Real_v A     = vv;
0163     Real_v B     = pv;
0164     Real_v C     = (rr - ellipsoid.fR * ellipsoid.fR);
0165     Real_v D     = B * B - A * C;
0166     Real_v sqrtD = vecCore::math::Sqrt(vecCore::math::Abs(D));
0167     Real_v trmin = (-B - sqrtD) / A;
0168     Real_v trmax = (-B + sqrtD) / A;
0169
0170     // Set preliminary distances to in/out
0171     Real_v tmin = vecCore::math::Max(tzmin, trmin);
0172     Real_v tmax = vecCore::math::Min(tzmax, trmax);
0173
0174     // Check if no intersection
0175     Real_v EPS  = Real_v(2.) * rr * vv * kEpsilon;
0176     Bool_v done = leaving || (D <= EPS) || ((tmax - tmin) <= kHalfTolerance);
0177
0178     // Set distance
0179     vecCore__MaskedAssignFunc(distance, !done, tmin + offset);
0180   }
0181
0182   template <typename Real_v>
0183   VECGEOM_FORCE_INLINE
0184   VECCORE_ATT_HOST_DEVICE
0185   static void DistanceToOut(UnplacedStruct_t const &ellipsoid, Vector3D<Real_v> const &point,
0186                             Vector3D<Real_v> const &direction, Real_v const & /* stepMax */, Real_v &distance)
0187   {
0188     /* TODO :  Logic to calculate Distance from inside point to the Ellipsoid surface */
0189     using Bool_v = vecCore::Mask_v<Real_v>;
0190
0191     // Scale ellipsoid to sphere
0192     Real_v px = point.x() * ellipsoid.fSx;
0193     Real_v py = point.y() * ellipsoid.fSy;
0194     Real_v pz = point.z() * ellipsoid.fSz;
0195     Real_v vx = direction.x() * ellipsoid.fSx;
0196     Real_v vy = direction.y() * ellipsoid.fSy;
0197     Real_v vz = direction.z() * ellipsoid.fSz;
0198
0199     // Check if point is outside ("wrong side")
0200     Real_v pzcut = pz - ellipsoid.fScZMidCut;
0201     Real_v dzcut = Real_v(ellipsoid.fScZDimCut);
0202     Real_v distZ = vecCore::math::Abs(pzcut) - dzcut;
0203
0204     Real_v rr      = px * px + py * py + pz * pz;
0205     Real_v vv      = vx * vx + vy * vy + vz * vz;
0206     Real_v pv      = px * vx + py * vy + pz * vz;
0207     Real_v distR   = ellipsoid.fQ1 * rr - ellipsoid.fQ2;
0208     Bool_v outside = vecCore::math::Max(distR, distZ) > kHalfTolerance;
0209
0210     distance = Real_v(0.);
0211     vecCore__MaskedAssignFunc(distance, outside, Real_v(-1.));
0212
0213     // Find intersection with Z planes
0214     Real_v tzmax = kMaximum;
0215     vecCore__MaskedAssignFunc(tzmax, vz != Real_v(0.), (vecCore::math::CopySign(Real_v(dzcut), vz) - pzcut) / vz);
0216
0217     // Find intersection with sphere
0218     Real_v B     = pv / vv;
0219     Real_v C     = (rr - ellipsoid.fR * ellipsoid.fR) / vv;
0220     Real_v D     = B * B - C;
0221     Real_v sqrtD = vecCore::math::Sqrt(vecCore::math::Abs(D));
0222     Real_v trmax = -B + sqrtD;
0223
0224     // Check if no intersection
0225     Real_v EPS  = Real_v(2.) * rr * vv * kEpsilon;
0226     Bool_v done = outside || (D <= EPS);
0227
0228     // Set distance
0229     vecCore__MaskedAssignFunc(distance, !done, vecCore::math::Min(tzmax, trmax));
0230   }
0231
0232   template <typename Real_v>
0233   VECGEOM_FORCE_INLINE
0234   VECCORE_ATT_HOST_DEVICE
0235   static void SafetyToIn(UnplacedStruct_t const &ellipsoid, Vector3D<Real_v> const &point, Real_v &safety)
0236   {
0237     /* TODO :  Logic to calculate Safety from outside point to the Ellipsoid surface */
0238     Real_v x = point.x() * ellipsoid.fSx;
0239     Real_v y = point.y() * ellipsoid.fSy;
0240     Real_v z = point.z() * ellipsoid.fSz;
0241     Real_v r = vecCore::math::Sqrt(x * x + y * y + z * z);
0242     // Set safety to zero if point is on surface
0243     Real_v safeZ = vecCore::math::Abs(z - ellipsoid.fScZMidCut) - ellipsoid.fScZDimCut;
0244     Real_v safeR = r - ellipsoid.fR;
0245     safety       = vecCore::math::Max(safeZ, safeR);
0246     vecCore::MaskedAssign(safety, vecCore::math::Abs(safety) <= kHalfTolerance, Real_v(0.));
0247     // Adjust safety using bounding box
0248     Real_v distZ  = vecCore::math::Max(point.z() - ellipsoid.fZTopCut, ellipsoid.fZBottomCut - point.z());
0249     Real_v distXY = vecCore::math::Max(vecCore::math::Abs(point.x()) - ellipsoid.fXmax,
0250                                        vecCore::math::Abs(point.y()) - ellipsoid.fYmax);
0251     vecCore__MaskedAssignFunc(safety, safety > Real_v(0.), vecCore::math::Max(safety, distZ, distXY));
0252   }
0253
0254   template <typename Real_v>
0255   VECGEOM_FORCE_INLINE
0256   VECCORE_ATT_HOST_DEVICE
0257   static void SafetyToOut(UnplacedStruct_t const &ellipsoid, Vector3D<Real_v> const &point, Real_v &safety)
0258   {
0259     /* TODO :  Logic to calculate Safety from inside point to the Ellipsoid surface */
0260     Real_v x = point.x() * ellipsoid.fSx;
0261     Real_v y = point.y() * ellipsoid.fSy;
0262     Real_v z = point.z() * ellipsoid.fSz;
0263     // Set safety to zero if point is on surface
0264     Real_v safeR = ellipsoid.fR - vecCore::math::Sqrt(x * x + y * y + z * z);
0265     Real_v safeZ = ellipsoid.fScZDimCut - vecCore::math::Abs(z - ellipsoid.fScZMidCut);
0266     safety       = vecCore::math::Min(safeZ, safeR);
0267     vecCore::MaskedAssign(safety, vecCore::math::Abs(safety) <= kHalfTolerance, Real_v(0.));
0268     // Adjust safety in z direction
0269     Real_v distZ = vecCore::math::Min(ellipsoid.fZTopCut - point.z(), point.z() - ellipsoid.fZBottomCut);
0270     vecCore__MaskedAssignFunc(safety, safety > Real_v(0.), vecCore::math::Min(safeR, distZ));
0271   }
0272
0273   template <typename Real_v>
0274   VECGEOM_FORCE_INLINE
0275   VECCORE_ATT_HOST_DEVICE
0276   static Vector3D<Real_v> NormalKernel(UnplacedStruct_t const &ellipsoid, Vector3D<Real_v> const &point,
0277                                        typename vecCore::Mask_v<Real_v> &valid)
0278   {
0279     // Computes the normal on a surface and returns it as a unit vector
0280     //   In case if the point is further than kHalfTolerance from the surface, set valid=false
0281     //   Must return a valid vector (even if the point is not on the surface)
0282     //
0283     //   On an edge provide an average normal of the corresponding base and lateral surface
0284     Vector3D<Real_v> normal(0.);
0285     valid = true;
0286
0287     Real_v px   = point.x();
0288     Real_v py   = point.y();
0289     Real_v pz   = point.z();
0290     Real_v A    = ellipsoid.fDx;
0291     Real_v B    = ellipsoid.fDy;
0292     Real_v C    = ellipsoid.fDz;
0293     Real_v x    = px * ellipsoid.fSx;
0294     Real_v y    = py * ellipsoid.fSy;
0295     Real_v z    = pz * ellipsoid.fSz;
0296     Real_v mag2 = x * x + y * y + z * z;
0297
0298     // Check lateral surface
0299     Real_v distR = ellipsoid.fQ1 * mag2 - ellipsoid.fQ2;
0300     vecCore__MaskedAssignFunc(normal, vecCore::math::Abs(distR) <= kHalfTolerance,
0301                               Vector3D<Real_v>(px / (A * A), py / (B * B), pz / (C * C)).Unit());
0302
0303     // Check z cuts
0304     Real_v distZ = vecCore::math::Abs(z - ellipsoid.fScZMidCut) - ellipsoid.fScZDimCut;
0305     vecCore__MaskedAssignFunc(normal[2], vecCore::math::Abs(distZ) <= kHalfTolerance,
0306                               normal[2] + vecCore::math::Sign(z - ellipsoid.fScZMidCut));
0307
0308     // Average normal, if required
0309     vecCore__MaskedAssignFunc(normal, normal.Mag2() > 1., normal.Unit());
0310     vecCore::Mask_v<Real_v> done = normal.Mag2() > Real_v(0.);
0311     if (vecCore::MaskFull(done)) return normal;
0312
0313     // Point is not on the surface - normally, this should never be
0314     // Return normal to the nearest surface
0315     vecCore__MaskedAssignFunc(valid, !done, false);
0316     vecCore__MaskedAssignFunc(normal[2], !done, vecCore::math::Sign(z - ellipsoid.fScZMidCut));
0317     vecCore__MaskedAssignFunc(distR, !done, vecCore::math::Sqrt(mag2) - ellipsoid.fR);
0318     vecCore__MaskedAssignFunc(normal, !done && distR > distZ && mag2 > Real_v(0.),
0319                               Vector3D<Real_v>(px / (A * A), py / (B * B), pz / (C * C)).Unit());
0320     return normal;
0321   }
0322 };
0323 } // namespace VECGEOM_IMPL_NAMESPACE
0324 } // namespace vecgeom
0325
0326 #endif // VECGEOM_VOLUMES_KERNEL_ELLIPSOIDIMPLEMENTATION_H_