volumes/kernel/HypeImplementation.h

0001 //===-- kernel/HypeImplementation.h - Instruction class definition -------*- C++ -*-===//
0002 //
0003 //                     GeantV - VecGeom
0004 //
0005 // This file is distributed under the LGPL
0006 // License. See LICENSE.TXT for details.
0007 //
0008 //===----------------------------------------------------------------------===//
0009 ///
0010 /// \file
0011 /// \author Marilena Bandieramonte (marilena.bandieramonte@cern.ch)
0012 /// \brief This file implements the Hype shape
0013 ///
0014
0015 #ifndef VECGEOM_VOLUMES_KERNEL_HYPEIMPLEMENTATION_H_
0016 #define VECGEOM_VOLUMES_KERNEL_HYPEIMPLEMENTATION_H_
0017
0018 #include "VecGeom/base/Vector3D.h"
0019 #include "VecGeom/volumes/HypeStruct.h"
0020 #include "VecGeom/volumes/kernel/GenericKernels.h"
0021 #include <VecCore/VecCore>
0022 #include "VecGeom/volumes/HypeUtilities.h"
0023 #include "VecGeom/volumes/kernel/shapetypes/HypeTypes.h"
0024
0025 // different SafetyToIn implementations
0026 // #define ACCURATE_BB
0027 #define ACCURATE_BC
0028
0029 namespace vecgeom {
0030
0031 VECGEOM_DEVICE_DECLARE_CONV_TEMPLATE(struct, HypeImplementation, typename);
0032
0033 inline namespace VECGEOM_IMPL_NAMESPACE {
0034
0035 template <typename T>
0036 struct HypeStruct;
0037
0038 template <typename T>
0039 class SPlacedHype;
0040 template <typename T>
0041 class SUnplacedHype;
0042
0043 template <typename hypeTypeT>
0044 struct HypeImplementation {
0045
0046   using UnplacedStruct_t = HypeStruct<Precision>;
0047   using UnplacedVolume_t = SUnplacedHype<hypeTypeT>;
0048   using PlacedShape_t    = SPlacedHype<UnplacedVolume_t>;
0049
0050   template <typename Real_v, typename Bool_v>
0051   VECGEOM_FORCE_INLINE VECCORE_ATT_HOST_DEVICE static void Contains(UnplacedStruct_t const &hype,
0052                                                                     Vector3D<Real_v> const &point, Bool_v &inside)
0053   {
0054     Bool_v unused(false), outside(false);
0055     GenericKernelForContainsAndInside<Real_v, Bool_v, false>(hype, point, unused, outside);
0056     inside = !outside;
0057   }
0058
0059   template <typename Real_v, typename Inside_v>
0060   VECGEOM_FORCE_INLINE VECCORE_ATT_HOST_DEVICE static void Inside(UnplacedStruct_t const &hype,
0061                                                                   Vector3D<Real_v> const &point, Inside_v &inside)
0062   {
0063     using Bool_v       = vecCore::Mask_v<Real_v>;
0064     using InsideBool_v = vecCore::Mask_v<Inside_v>;
0065     Bool_v completelyinside(false), completelyoutside(false);
0066     GenericKernelForContainsAndInside<Real_v, Bool_v, true>(hype, point, completelyinside, completelyoutside);
0067     inside = EInside::kSurface;
0068     vecCore::MaskedAssign(inside, (InsideBool_v)completelyoutside, Inside_v(EInside::kOutside));
0069     vecCore::MaskedAssign(inside, (InsideBool_v)completelyinside, Inside_v(EInside::kInside));
0070   }
0071
0072   template <typename Real_v, typename Bool_v, bool ForInside>
0073   VECGEOM_FORCE_INLINE VECCORE_ATT_HOST_DEVICE static void GenericKernelForContainsAndInside(
0074       UnplacedStruct_t const &hype, Vector3D<Real_v> const &point, Bool_v &completelyinside, Bool_v &completelyoutside)
0075   {
0076     using namespace ::vecgeom::HypeTypes;
0077     Real_v r2    = point.Perp2();
0078     Real_v oRad2 = (hype.fRmax2 + hype.fTOut2 * point.z() * point.z());
0079     Real_v iRad2(0.);
0080
0081     completelyoutside = (Abs(point.z()) > (hype.fDz + hype.zToleranceLevel));
0082     if (vecCore::MaskFull(completelyoutside)) return;
0083     completelyoutside |= (r2 > oRad2 + hype.outerRadToleranceLevel);
0084     if (vecCore::MaskFull(completelyoutside)) return;
0085     if (checkInnerSurfaceTreatment<hypeTypeT>(hype)) {
0086       iRad2 = (hype.fRmin2 + hype.fTIn2 * point.z() * point.z());
0087       completelyoutside |= (r2 < (iRad2 - hype.innerRadToleranceLevel));
0088     }
0089     if (vecCore::MaskFull(completelyoutside)) return;
0090
0091     if (ForInside) {
0092       completelyinside =
0093           (Abs(point.z()) < (hype.fDz - hype.zToleranceLevel)) && (r2 < oRad2 - hype.outerRadToleranceLevel);
0094
0095       if (checkInnerSurfaceTreatment<hypeTypeT>(hype)) completelyinside &= (r2 > (iRad2 + hype.innerRadToleranceLevel));
0096     }
0097   }
0098
0099   template <typename Real_v>
0100   VECGEOM_FORCE_INLINE VECCORE_ATT_HOST_DEVICE static void DistanceToIn(UnplacedStruct_t const &hype,
0101                                                                         Vector3D<Real_v> const &point,
0102                                                                         Vector3D<Real_v> const &direction,
0103                                                                         Real_v const & /*stepMax*/, Real_v &distance)
0104   {
0105     using namespace ::vecgeom::HypeTypes;
0106     using Bool_v = vecCore::Mask_v<Real_v>;
0107     Real_v absZ  = Abs(point.z());
0108     distance     = kInfLength;
0109     Real_v zDist(kInfLength), dist(kInfLength);
0110     Real_v r = point.Perp2();
0111
0112     Bool_v done(false);
0113     Bool_v cond(false);
0114
0115     Bool_v surfaceCond = HypeUtilities::IsPointOnSurfaceAndMovingInside<Real_v, hypeTypeT>(hype, point, direction);
0116     vecCore__MaskedAssignFunc(distance, !done && surfaceCond, Real_v(0.0));
0117     done |= surfaceCond;
0118     if (vecCore::MaskFull(done)) return;
0119
0120     cond = HypeUtilities::IsCompletelyInside<Real_v, hypeTypeT>(hype, point);
0121     vecCore__MaskedAssignFunc(distance, !done && cond, Real_v(-1.0));
0122     done |= cond;
0123     if (vecCore::MaskFull(done)) return;
0124
0125     // checking whether point hits the Z Surface of hyperboloid
0126     Bool_v hittingZPlane =
0127         HypeUtilities::GetPointOfIntersectionWithZPlane<Real_v, hypeTypeT, true>(hype, point, direction, zDist);
0128     Bool_v isPointAboveOrBelowHypeAndGoingInside = (absZ > hype.fDz) && (point.z() * direction.z() < Real_v(0.));
0129     cond                                         = isPointAboveOrBelowHypeAndGoingInside && hittingZPlane;
0130     vecCore::MaskedAssign(distance, !done && cond, zDist);
0131     done |= cond;
0132     if (vecCore::MaskFull(done)) return;
0133
0134     // Moving the point to Z Surface
0135     Vector3D<Real_v> newPt = point + zDist * direction;
0136     Real_v rp2             = newPt.Perp2();
0137
0138     Bool_v hittingOuterSurfaceFromOutsideZRange =
0139         isPointAboveOrBelowHypeAndGoingInside && !hittingZPlane && (rp2 >= hype.fEndOuterRadius2);
0140     Bool_v hittingOuterSurfaceFromWithinZRange = ((r > ((hype.fRmax2 + hype.fTOut2 * absZ * absZ) + kHalfTolerance))) &&
0141                                                  (absZ >= Real_v(0.)) && (absZ <= hype.fDz);
0142
0143     cond = (hittingOuterSurfaceFromOutsideZRange || hittingOuterSurfaceFromWithinZRange ||
0144             (HypeUtilities::IsPointOnOuterSurfaceAndMovingOutside<Real_v>(hype, point, direction))) &&
0145            HypeHelpers<Real_v, true, false>::GetPointOfIntersectionWithHyperbolicSurface(hype, point, direction, dist);
0146     vecCore::MaskedAssign(distance, !done && cond, dist);
0147
0148     if (checkInnerSurfaceTreatment<hypeTypeT>(hype)) {
0149       done |= cond;
0150       if (vecCore::MaskFull(done)) return;
0151       Bool_v hittingInnerSurfaceFromOutsideZRange =
0152           isPointAboveOrBelowHypeAndGoingInside && !hittingZPlane && (rp2 <= hype.fEndInnerRadius2);
0153       Bool_v hittingInnerSurfaceFromWithinZRange = (r < ((hype.fRmin2 + hype.fTIn2 * absZ * absZ) - kHalfTolerance)) &&
0154                                                    (absZ >= Real_v(0.)) && (absZ <= hype.fDz);
0155
0156       // If it hits inner hyperbolic surface then distance will be the distance to inner hyperbolic surface
0157       // Or if the point is on the inner Hyperbolic surface but going out then the distance will be the distance to
0158       // opposite inner hyperbolic surface.
0159       cond = (hittingInnerSurfaceFromOutsideZRange || hittingInnerSurfaceFromWithinZRange ||
0160               (HypeUtilities::IsPointOnInnerSurfaceAndMovingOutside<Real_v, hypeTypeT>(hype, point, direction))) &&
0161              HypeHelpers<Real_v, true, true>::GetPointOfIntersectionWithHyperbolicSurface(hype, point, direction, dist);
0162       vecCore::MaskedAssign(distance, !done && cond, dist);
0163     }
0164   }
0165
0166   template <typename Real_v>
0167   VECGEOM_FORCE_INLINE VECCORE_ATT_HOST_DEVICE static void DistanceToOut(UnplacedStruct_t const &hype,
0168                                                                          Vector3D<Real_v> const &point,
0169                                                                          Vector3D<Real_v> const &direction,
0170                                                                          Real_v const & /* stepMax */, Real_v &distance)
0171   {
0172     using namespace ::vecgeom::HypeTypes;
0173     using Bool_v = typename vecCore::Mask_v<Real_v>;
0174     distance     = kInfLength;
0175     Real_v zDist(kInfLength), dist(kInfLength);
0176
0177     Bool_v done(false);
0178
0179     Bool_v cond = HypeUtilities::IsPointOnSurfaceAndMovingOutside<Real_v, hypeTypeT>(hype, point, direction);
0180     vecCore__MaskedAssignFunc(distance, cond, Real_v(0.));
0181     done |= cond;
0182     if (vecCore::MaskFull(done)) return;
0183
0184     cond = HypeUtilities::IsCompletelyOutside<Real_v, hypeTypeT>(hype, point);
0185     vecCore__MaskedAssignFunc(distance, !done && cond, Real_v(-1.));
0186     done |= cond;
0187     if (vecCore::MaskFull(done)) return;
0188
0189     HypeUtilities::GetPointOfIntersectionWithZPlane<Real_v, hypeTypeT, false>(hype, point, direction, zDist);
0190     vecCore__MaskedAssignFunc(zDist, zDist < Real_v(0.), InfinityLength<Real_v>());
0191
0192     HypeHelpers<Real_v, false, false>::GetPointOfIntersectionWithHyperbolicSurface(hype, point, direction, dist);
0193     vecCore__MaskedAssignFunc(dist, dist < Real_v(0.), InfinityLength<Real_v>());
0194     vecCore__MaskedAssignFunc(distance, !done, Min(zDist, dist));
0195
0196     if (checkInnerSurfaceTreatment<hypeTypeT>(hype)) {
0197       HypeHelpers<Real_v, false, true>::GetPointOfIntersectionWithHyperbolicSurface(hype, point, direction, dist);
0198       vecCore__MaskedAssignFunc(dist, dist < Real_v(0.), InfinityLength<Real_v>());
0199       vecCore__MaskedAssignFunc(distance, !done, Min(distance, dist));
0200     }
0201   }
0202
0203   template <class Real_v>
0204   VECCORE_ATT_HOST_DEVICE static void SafetyToIn(UnplacedStruct_t const &hype, Vector3D<Real_v> const &point,
0205                                                  Real_v &safety)
0206   {
0207
0208     using Bool_v = typename vecCore::Mask_v<Real_v>;
0209
0210     Real_v absZ = Abs(point.z());
0211     Real_v r2   = point.Perp2();
0212     Real_v r    = Sqrt(r2);
0213
0214     Bool_v done(false);
0215
0216     // New Simple Algo
0217     safety = 0.;
0218     // If point is inside then safety should be -1.
0219     Bool_v compIn(false), compOut(false);
0220     GenericKernelForContainsAndInside<Real_v, Bool_v, true>(hype, point, compIn, compOut);
0221     done = (!compIn && !compOut);
0222     if (vecCore::MaskFull(done)) return;
0223
0224     vecCore__MaskedAssignFunc(safety, compIn, Real_v(-1.0));
0225     done |= compIn;
0226     if (vecCore::MaskFull(done)) return;
0227
0228     Bool_v cond(false);
0229     Real_v sigz = absZ - hype.fDz;
0230     cond        = (sigz > kHalfTolerance) && (r < hype.fEndOuterRadius) && (r > hype.fEndInnerRadius);
0231     vecCore::MaskedAssign(safety, !done && cond, sigz);
0232     done |= cond;
0233     if (vecCore::MaskFull(done)) return;
0234
0235     cond = (sigz > kHalfTolerance) && (r > hype.fEndOuterRadius);
0236     vecCore__MaskedAssignFunc(safety, !done && cond,
0237                               Sqrt((r - hype.fEndOuterRadius) * (r - hype.fEndOuterRadius) + (sigz) * (sigz)));
0238     done |= cond;
0239     if (vecCore::MaskFull(done)) return;
0240
0241     cond = (sigz > kHalfTolerance) && (r < hype.fEndInnerRadius);
0242     vecCore__MaskedAssignFunc(safety, !done && cond,
0243                               Sqrt((r - hype.fEndInnerRadius) * (r - hype.fEndInnerRadius) + (sigz) * (sigz)));
0244     done |= cond;
0245     if (vecCore::MaskFull(done)) return;
0246
0247     cond =
0248         (r2 > ((hype.fRmax2 + hype.fTOut2 * absZ * absZ) + kHalfTolerance)) && (absZ > Real_v(0.)) && (absZ < hype.fDz);
0249     vecCore__MaskedAssignFunc(safety, !done && cond,
0250                               HypeUtilities::ApproxDistOutside<Real_v>(r, absZ, hype.fRmax, hype.fTOut));
0251     done |= cond;
0252     if (vecCore::MaskFull(done)) return;
0253
0254     vecCore__MaskedAssignFunc(safety,
0255                               !done && (r2 < ((hype.fRmin2 + hype.fTIn2 * absZ * absZ) - kHalfTolerance)) &&
0256                                   (absZ > Real_v(0.)) && (absZ < hype.fDz),
0257                               HypeUtilities::ApproxDistInside<Real_v>(r, absZ, hype.fRmin, hype.fTIn2));
0258   }
0259
0260   template <class Real_v>
0261   VECCORE_ATT_HOST_DEVICE static void SafetyToOut(UnplacedStruct_t const &hype, Vector3D<Real_v> const &point,
0262                                                   Real_v &safety)
0263   {
0264     using namespace ::vecgeom::HypeTypes;
0265     using Bool_v = typename vecCore::Mask_v<Real_v>;
0266     safety       = 0.;
0267     Real_v r     = Sqrt(point.x() * point.x() + point.y() * point.y());
0268     Real_v absZ  = Abs(point.z());
0269     Bool_v inside(false), outside(false);
0270     Bool_v done(false);
0271
0272     Real_v distZ(0.), distInner(0.), distOuter(0.);
0273     safety = 0.;
0274     GenericKernelForContainsAndInside<Real_v, Bool_v, true>(hype, point, inside, outside);
0275     done = (!inside && !outside);
0276     if (vecCore::MaskFull(done)) return;
0277
0278     vecCore__MaskedAssignFunc(safety, outside, Real_v(-1.0));
0279     done |= outside;
0280     if (vecCore::MaskFull(done)) return;
0281
0282     vecCore__MaskedAssignFunc(distZ, !done && inside, Abs(Abs(point.z()) - hype.fDz));
0283     if (checkInnerSurfaceTreatment<hypeTypeT>(hype) && hype.fStIn) {
0284       vecCore__MaskedAssignFunc(distInner, !done && inside,
0285                                 HypeUtilities::ApproxDistOutside<Real_v>(r, absZ, hype.fRmin, hype.fTIn));
0286     }
0287
0288     if (checkInnerSurfaceTreatment<hypeTypeT>(hype) && !hype.fStIn) {
0289       vecCore__MaskedAssignFunc(distInner, !done && inside, (r - hype.fRmin));
0290     }
0291
0292     if (!checkInnerSurfaceTreatment<hypeTypeT>(hype) && !hype.fStIn) {
0293       vecCore__MaskedAssignFunc(distInner, !done && inside, InfinityLength<Real_v>());
0294     }
0295
0296     vecCore__MaskedAssignFunc(distOuter, !done && inside,
0297                               HypeUtilities::ApproxDistInside<Real_v>(r, absZ, hype.fRmax, hype.fTOut2));
0298     safety = Min(distInner, distOuter);
0299     safety = Min(safety, distZ);
0300   }
0301 };
0302 } // namespace VECGEOM_IMPL_NAMESPACE
0303 } // namespace vecgeom
0304
0305 #endif // VECGEOM_VOLUMES_KERNEL_HYPEIMPLEMENTATION_H_