VecGeom/volumes/SphereUtilities.h

0001 #ifndef VECGEOM_VOLUMES_SPHEREUTILITIES_H_
0002 #define VECGEOM_VOLUMES_SPHEREUTILITIES_H_
0003
0004 #include "VecGeom/base/Global.h"
0005
0006 #ifndef VECCORE_CUDA
0007 #include "VecGeom/base/RNG.h"
0008 #endif
0009
0010 #include "VecGeom/base/Global.h"
0011 #include "VecGeom/volumes/Wedge_Evolution.h"
0012 #include "VecGeom/base/Vector3D.h"
0013 #include "VecGeom/volumes/SphereStruct.h"
0014 #include "VecGeom/volumes/kernel/GenericKernels.h"
0015 #include <cstdio>
0016
0017 namespace vecgeom {
0018 inline namespace VECGEOM_IMPL_NAMESPACE {
0019
0020 class UnplacedSphere;
0021 template <typename T>
0022 struct SphereStruct;
0023
0024 template <typename T>
0025 VECGEOM_FORCE_INLINE
0026 VECCORE_ATT_HOST_DEVICE
0027 T sqr(T x)
0028 {
0029   return x * x;
0030 }
0031
0032 #ifndef VECCORE_CUDA
0033 // Generate radius in annular ring according to uniform area
0034 template <typename T>
0035 inline T GetRadiusInRing(T rmin, T rmax)
0036 {
0037   if (rmin == rmax) return rmin;
0038
0039   T rng(RNG::Instance().uniform(0.0, 1.0));
0040
0041   if (rmin <= T(0.0)) return rmax * Sqrt(rng);
0042
0043   T rmin2 = rmin * rmin;
0044   T rmax2 = rmax * rmax;
0045
0046   return Sqrt(rng * (rmax2 - rmin2) + rmin2);
0047 }
0048 #endif
0049
0050 namespace SphereUtilities {
0051 using UnplacedStruct_t = SphereStruct<Precision>;
0052
0053 template <class Real_v>
0054 VECCORE_ATT_HOST_DEVICE
0055 typename vecCore::Mask_v<Real_v> IsPointOnInnerRadius(UnplacedStruct_t const &unplaced, Vector3D<Real_v> const &point)
0056 {
0057   auto mag2 = point.Mag2();
0058   return mag2 <= MakePlusTolerantSquare<true>(unplaced.fRmin) &&
0059          mag2 >= MakeMinusTolerantSquare<true>(unplaced.fRmin);
0060 }
0061
0062 template <class Real_v>
0063 VECCORE_ATT_HOST_DEVICE
0064 typename vecCore::Mask_v<Real_v> IsPointOnOuterRadius(UnplacedStruct_t const &unplaced, Vector3D<Real_v> const &point)
0065 {
0066   auto mag2 = point.Mag2();
0067   return mag2 <= MakePlusTolerantSquare<true>(unplaced.fRmax) &&
0068          mag2 >= MakeMinusTolerantSquare<true>(unplaced.fRmax);
0069 }
0070
0071 template <class Real_v>
0072 VECGEOM_FORCE_INLINE
0073 VECCORE_ATT_HOST_DEVICE
0074 typename vecCore::Mask_v<Real_v> IsPointOnStartPhi(UnplacedStruct_t const &unplaced, Vector3D<Real_v> const &point)
0075 {
0076
0077   return unplaced.fPhiWedge.IsOnSurfaceGeneric<Real_v>(unplaced.fPhiWedge.GetAlong1(), unplaced.fPhiWedge.GetNormal1(),
0078                                                        point);
0079 }
0080
0081 template <class Real_v>
0082 VECGEOM_FORCE_INLINE
0083 VECCORE_ATT_HOST_DEVICE
0084 typename vecCore::Mask_v<Real_v> IsPointOnEndPhi(UnplacedStruct_t const &unplaced, Vector3D<Real_v> const &point)
0085 {
0086
0087   return unplaced.fPhiWedge.IsOnSurfaceGeneric<Real_v>(unplaced.fPhiWedge.GetAlong2(), unplaced.fPhiWedge.GetNormal2(),
0088                                                        point);
0089 }
0090
0091 template <class Real_v>
0092 VECGEOM_FORCE_INLINE
0093 VECCORE_ATT_HOST_DEVICE
0094 typename vecCore::Mask_v<Real_v> IsPointOnStartTheta(UnplacedStruct_t const &unplaced, Vector3D<Real_v> const &point)
0095 {
0096
0097   return unplaced.fThetaCone.IsOnSurfaceGeneric<Real_v, true>(point);
0098 }
0099
0100 template <class Real_v>
0101 VECGEOM_FORCE_INLINE
0102 VECCORE_ATT_HOST_DEVICE
0103 typename vecCore::Mask_v<Real_v> IsPointOnEndTheta(UnplacedStruct_t const &unplaced, Vector3D<Real_v> const &point)
0104 {
0105
0106   return unplaced.fThetaCone.IsOnSurfaceGeneric<Real_v, false>(point);
0107 }
0108
0109 template <class Real_v>
0110 VECGEOM_FORCE_INLINE
0111 VECCORE_ATT_HOST_DEVICE
0112 typename vecCore::Mask_v<Real_v> IsCompletelyOutside(UnplacedStruct_t const &unplaced,
0113                                                      Vector3D<Real_v> const &localPoint)
0114 {
0115
0116   using Bool_v              = vecCore::Mask_v<Real_v>;
0117   Real_v rad                = localPoint.Mag();
0118   Bool_v outsideRadiusRange = (rad > (unplaced.fRmax + kTolerance)) || (rad < (unplaced.fRmin - kTolerance));
0119   Bool_v outsidePhiRange(false), insidePhiRange(false);
0120   unplaced.fPhiWedge.GenericKernelForContainsAndInside<Real_v, true>(localPoint, insidePhiRange, outsidePhiRange);
0121   Bool_v outsideThetaRange = unplaced.fThetaCone.IsCompletelyOutside<Real_v>(localPoint);
0122   Bool_v completelyoutside = outsideRadiusRange || outsidePhiRange || outsideThetaRange;
0123   return completelyoutside;
0124 }
0125
0126 template <class Real_v>
0127 VECGEOM_FORCE_INLINE
0128 VECCORE_ATT_HOST_DEVICE
0129 typename vecCore::Mask_v<Real_v> IsCompletelyInside(UnplacedStruct_t const &unplaced,
0130                                                     Vector3D<Real_v> const &localPoint)
0131 {
0132   using Bool_v             = vecCore::Mask_v<Real_v>;
0133   Real_v rad               = localPoint.Mag();
0134   Bool_v insideRadiusRange = (rad < (unplaced.fRmax - kTolerance)) && (rad > (unplaced.fRmin + kTolerance));
0135   Bool_v outsidePhiRange(false), insidePhiRange(false);
0136   unplaced.fPhiWedge.GenericKernelForContainsAndInside<Real_v, true>(localPoint, insidePhiRange, outsidePhiRange);
0137   Bool_v insideThetaRange = unplaced.fThetaCone.IsCompletelyInside<Real_v>(localPoint);
0138   Bool_v completelyinside = insideRadiusRange && insidePhiRange && insideThetaRange;
0139   return completelyinside;
0140 }
0141
0142 template <class Real_v, bool ForInnerRadius, bool MovingOut>
0143 VECCORE_ATT_HOST_DEVICE
0144 typename vecCore::Mask_v<Real_v> IsPointOnRadialSurfaceAndMovingOut(UnplacedStruct_t const &unplaced,
0145                                                                     Vector3D<Real_v> const &point,
0146                                                                     Vector3D<Real_v> const &dir)
0147 {
0148   // Rays from rmax+tolerance or rmin-tolerance can be moving out even if not going fully "backward"
0149   if (MovingOut) {
0150     if (ForInnerRadius) {
0151       return IsPointOnInnerRadius<Real_v>(unplaced, point) && (dir.Dot(-point) > Real_v(kSqrtTolerance));
0152     } else {
0153       return IsPointOnOuterRadius<Real_v>(unplaced, point) && (dir.Dot(point) > Real_v(kSqrtTolerance));
0154     }
0155   } else {
0156     if (ForInnerRadius) {
0157       return IsPointOnInnerRadius<Real_v>(unplaced, point) && (dir.Dot(-point) < Real_v(-kSqrtTolerance));
0158     } else
0159       return IsPointOnOuterRadius<Real_v>(unplaced, point) && (dir.Dot(point) < Real_v(-kSqrtTolerance));
0160   }
0161 }
0162
0163 template <class Real_v, bool MovingOut>
0164 VECGEOM_FORCE_INLINE
0165 VECCORE_ATT_HOST_DEVICE
0166 typename vecCore::Mask_v<Real_v> IsPointOnSurfaceAndMovingOut(UnplacedStruct_t const &unplaced,
0167                                                               Vector3D<Real_v> const &point,
0168                                                               Vector3D<Real_v> const &dir)
0169 {
0170
0171   using Bool_v = vecCore::Mask_v<Real_v>;
0172
0173   Bool_v tempOuterRad = IsPointOnRadialSurfaceAndMovingOut<Real_v, false, MovingOut>(unplaced, point, dir);
0174   Bool_v tempInnerRad(false), tempStartPhi(false), tempEndPhi(false), tempStartTheta(false), tempEndTheta(false);
0175   if (unplaced.fRmin) tempInnerRad = IsPointOnRadialSurfaceAndMovingOut<Real_v, true, MovingOut>(unplaced, point, dir);
0176   if (unplaced.fDPhi < (kTwoPi - kHalfTolerance)) {
0177     tempStartPhi = unplaced.fPhiWedge.IsPointOnSurfaceAndMovingOut<Real_v, true, MovingOut>(point, dir);
0178     tempEndPhi   = unplaced.fPhiWedge.IsPointOnSurfaceAndMovingOut<Real_v, false, MovingOut>(point, dir);
0179   }
0180   if (unplaced.fDTheta < (kPi - kHalfTolerance)) {
0181     tempStartTheta = unplaced.fThetaCone.IsPointOnSurfaceAndMovingOut<Real_v, true, MovingOut>(point, dir);
0182     tempEndTheta   = unplaced.fThetaCone.IsPointOnSurfaceAndMovingOut<Real_v, false, MovingOut>(point, dir);
0183   }
0184
0185   Bool_v isPointOnSurfaceAndMovingOut =
0186       ((tempOuterRad || tempInnerRad) && unplaced.fPhiWedge.Contains<Real_v>(point) &&
0187        unplaced.fThetaCone.Contains<Real_v>(point)) ||
0188       ((tempStartPhi || tempEndPhi) && (point.Mag2() >= unplaced.fRmin * unplaced.fRmin) &&
0189        (point.Mag2() <= unplaced.fRmax * unplaced.fRmax) && unplaced.fThetaCone.Contains<Real_v>(point)) ||
0190       ((tempStartTheta || tempEndTheta) && (point.Mag2() >= unplaced.fRmin * unplaced.fRmin) &&
0191        (point.Mag2() <= unplaced.fRmax * unplaced.fRmax) && unplaced.fPhiWedge.Contains<Real_v>(point));
0192
0193   return isPointOnSurfaceAndMovingOut;
0194 }
0195
0196 } // namespace SphereUtilities
0197 } // namespace VECGEOM_IMPL_NAMESPACE
0198 } // namespace vecgeom
0199
0200 #endif // VECGEOM_VOLUMES_SPHEREUTILITIES_H_