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