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   VECCORE_ATT_HOST_DEVICE
0051   static void PrintType()
0052   {
0053     //  printf("SpecializedHype<%i, %i>", transCodeT, rotCodeT);
0054   }
0055
0056   template <typename Stream>
0057   static void PrintType(Stream &st, int transCodeT = translation::kGeneric, int rotCodeT = rotation::kGeneric)
0058   {
0059     st << "SpecializedHype<" << transCodeT << "," << rotCodeT << ">";
0060   }
0061
0062   template <typename Stream>
0063   static void PrintImplementationType(Stream &st)
0064   {
0065     (void)st;
0066     // st << "HypeImplementation<" << transCodeT << "," << rotCodeT << ">";
0067   }
0068
0069   template <typename Stream>
0070   static void PrintUnplacedType(Stream &st)
0071   {
0072     (void)st;
0073     // TODO: this is wrong
0074     // st << "UnplacedHype";
0075   }
0076
0077   template <typename Real_v, typename Bool_v>
0078   VECGEOM_FORCE_INLINE
0079   VECCORE_ATT_HOST_DEVICE
0080   static void Contains(UnplacedStruct_t const &hype, Vector3D<Real_v> const &point, Bool_v &inside)
0081   {
0082     Bool_v unused(false), outside(false);
0083     GenericKernelForContainsAndInside<Real_v, Bool_v, false>(hype, point, unused, outside);
0084     inside = !outside;
0085   }
0086
0087   template <typename Real_v, typename Inside_v>
0088   VECGEOM_FORCE_INLINE
0089   VECCORE_ATT_HOST_DEVICE
0090   static void Inside(UnplacedStruct_t const &hype, Vector3D<Real_v> const &point, Inside_v &inside)
0091   {
0092     using Bool_v       = vecCore::Mask_v<Real_v>;
0093     using InsideBool_v = vecCore::Mask_v<Inside_v>;
0094     Bool_v completelyinside(false), completelyoutside(false);
0095     GenericKernelForContainsAndInside<Real_v, Bool_v, true>(hype, point, completelyinside, completelyoutside);
0096     inside = EInside::kSurface;
0097     vecCore::MaskedAssign(inside, (InsideBool_v)completelyoutside, Inside_v(EInside::kOutside));
0098     vecCore::MaskedAssign(inside, (InsideBool_v)completelyinside, Inside_v(EInside::kInside));
0099   }
0100
0101   template <typename Real_v, typename Bool_v, bool ForInside>
0102   VECGEOM_FORCE_INLINE
0103   VECCORE_ATT_HOST_DEVICE
0104   static void GenericKernelForContainsAndInside(UnplacedStruct_t const &hype, Vector3D<Real_v> const &point,
0105                                                 Bool_v &completelyinside, Bool_v &completelyoutside)
0106   {
0107     using namespace ::vecgeom::HypeTypes;
0108     Real_v r2    = point.Perp2();
0109     Real_v oRad2 = (hype.fRmax2 + hype.fTOut2 * point.z() * point.z());
0110     Real_v iRad2(0.);
0111
0112     completelyoutside = (Abs(point.z()) > (hype.fDz + hype.zToleranceLevel));
0113     if (vecCore::MaskFull(completelyoutside)) return;
0114     completelyoutside |= (r2 > oRad2 + hype.outerRadToleranceLevel);
0115     if (vecCore::MaskFull(completelyoutside)) return;
0116     if (checkInnerSurfaceTreatment<hypeTypeT>(hype)) {
0117       iRad2 = (hype.fRmin2 + hype.fTIn2 * point.z() * point.z());
0118       completelyoutside |= (r2 < (iRad2 - hype.innerRadToleranceLevel));
0119     }
0120     if (vecCore::MaskFull(completelyoutside)) return;
0121
0122     if (ForInside) {
0123       completelyinside =
0124           (Abs(point.z()) < (hype.fDz - hype.zToleranceLevel)) && (r2 < oRad2 - hype.outerRadToleranceLevel);
0125
0126       if (checkInnerSurfaceTreatment<hypeTypeT>(hype)) completelyinside &= (r2 > (iRad2 + hype.innerRadToleranceLevel));
0127     }
0128   }
0129
0130   template <typename Real_v>
0131   VECGEOM_FORCE_INLINE
0132   VECCORE_ATT_HOST_DEVICE
0133   static void DistanceToIn(UnplacedStruct_t const &hype, Vector3D<Real_v> const &point,
0134                            Vector3D<Real_v> const &direction, Real_v const & /*stepMax*/, Real_v &distance)
0135   {
0136     using namespace ::vecgeom::HypeTypes;
0137     using Bool_v = vecCore::Mask_v<Real_v>;
0138     Real_v absZ  = Abs(point.z());
0139     distance     = kInfLength;
0140     Real_v zDist(kInfLength), dist(kInfLength);
0141     Real_v r = point.Perp2();
0142
0143     Bool_v done(false);
0144     Bool_v cond(false);
0145
0146     Bool_v surfaceCond = HypeUtilities::IsPointOnSurfaceAndMovingInside<Real_v, hypeTypeT>(hype, point, direction);
0147     vecCore__MaskedAssignFunc(distance, !done && surfaceCond, Real_v(0.0));
0148     done |= surfaceCond;
0149     if (vecCore::MaskFull(done)) return;
0150
0151     cond = HypeUtilities::IsCompletelyInside<Real_v, hypeTypeT>(hype, point);
0152     vecCore__MaskedAssignFunc(distance, !done && cond, Real_v(-1.0));
0153     done |= cond;
0154     if (vecCore::MaskFull(done)) return;
0155
0156     // checking whether point hits the Z Surface of hyperboloid
0157     Bool_v hittingZPlane =
0158         HypeUtilities::GetPointOfIntersectionWithZPlane<Real_v, hypeTypeT, true>(hype, point, direction, zDist);
0159     Bool_v isPointAboveOrBelowHypeAndGoingInside = (absZ > hype.fDz) && (point.z() * direction.z() < Real_v(0.));
0160     cond                                         = isPointAboveOrBelowHypeAndGoingInside && hittingZPlane;
0161     vecCore::MaskedAssign(distance, !done && cond, zDist);
0162     done |= cond;
0163     if (vecCore::MaskFull(done)) return;
0164
0165     // Moving the point to Z Surface
0166     Vector3D<Real_v> newPt = point + zDist * direction;
0167     Real_v rp2             = newPt.Perp2();
0168
0169     Bool_v hittingOuterSurfaceFromOutsideZRange =
0170         isPointAboveOrBelowHypeAndGoingInside && !hittingZPlane && (rp2 >= hype.fEndOuterRadius2);
0171     Bool_v hittingOuterSurfaceFromWithinZRange = ((r > ((hype.fRmax2 + hype.fTOut2 * absZ * absZ) + kHalfTolerance))) &&
0172                                                  (absZ >= Real_v(0.)) && (absZ <= hype.fDz);
0173
0174     cond = (hittingOuterSurfaceFromOutsideZRange || hittingOuterSurfaceFromWithinZRange ||
0175             (HypeUtilities::IsPointOnOuterSurfaceAndMovingOutside<Real_v>(hype, point, direction))) &&
0176            HypeHelpers<Real_v, true, false>::GetPointOfIntersectionWithHyperbolicSurface(hype, point, direction, dist);
0177     vecCore::MaskedAssign(distance, !done && cond, dist);
0178
0179     if (checkInnerSurfaceTreatment<hypeTypeT>(hype)) {
0180       done |= cond;
0181       if (vecCore::MaskFull(done)) return;
0182       Bool_v hittingInnerSurfaceFromOutsideZRange =
0183           isPointAboveOrBelowHypeAndGoingInside && !hittingZPlane && (rp2 <= hype.fEndInnerRadius2);
0184       Bool_v hittingInnerSurfaceFromWithinZRange = (r < ((hype.fRmin2 + hype.fTIn2 * absZ * absZ) - kHalfTolerance)) &&
0185                                                    (absZ >= Real_v(0.)) && (absZ <= hype.fDz);
0186
0187       // If it hits inner hyperbolic surface then distance will be the distance to inner hyperbolic surface
0188       // Or if the point is on the inner Hyperbolic surface but going out then the distance will be the distance to
0189       // opposite inner hyperbolic surface.
0190       cond = (hittingInnerSurfaceFromOutsideZRange || hittingInnerSurfaceFromWithinZRange ||
0191               (HypeUtilities::IsPointOnInnerSurfaceAndMovingOutside<Real_v, hypeTypeT>(hype, point, direction))) &&
0192              HypeHelpers<Real_v, true, true>::GetPointOfIntersectionWithHyperbolicSurface(hype, point, direction, dist);
0193       vecCore::MaskedAssign(distance, !done && cond, dist);
0194     }
0195   }
0196
0197   template <typename Real_v>
0198   VECGEOM_FORCE_INLINE
0199   VECCORE_ATT_HOST_DEVICE
0200   static void DistanceToOut(UnplacedStruct_t const &hype, Vector3D<Real_v> const &point,
0201                             Vector3D<Real_v> const &direction, Real_v const & /* stepMax */, Real_v &distance)
0202   {
0203     using namespace ::vecgeom::HypeTypes;
0204     using Bool_v = typename vecCore::Mask_v<Real_v>;
0205     distance     = kInfLength;
0206     Real_v zDist(kInfLength), dist(kInfLength);
0207
0208     Bool_v done(false);
0209
0210     Bool_v cond = HypeUtilities::IsPointOnSurfaceAndMovingOutside<Real_v, hypeTypeT>(hype, point, direction);
0211     vecCore__MaskedAssignFunc(distance, cond, Real_v(0.));
0212     done |= cond;
0213     if (vecCore::MaskFull(done)) return;
0214
0215     cond = HypeUtilities::IsCompletelyOutside<Real_v, hypeTypeT>(hype, point);
0216     vecCore__MaskedAssignFunc(distance, !done && cond, Real_v(-1.));
0217     done |= cond;
0218     if (vecCore::MaskFull(done)) return;
0219
0220     HypeUtilities::GetPointOfIntersectionWithZPlane<Real_v, hypeTypeT, false>(hype, point, direction, zDist);
0221     vecCore__MaskedAssignFunc(zDist, zDist < Real_v(0.), InfinityLength<Real_v>());
0222
0223     HypeHelpers<Real_v, false, false>::GetPointOfIntersectionWithHyperbolicSurface(hype, point, direction, dist);
0224     vecCore__MaskedAssignFunc(dist, dist < Real_v(0.), InfinityLength<Real_v>());
0225     vecCore__MaskedAssignFunc(distance, !done, Min(zDist, dist));
0226
0227     if (checkInnerSurfaceTreatment<hypeTypeT>(hype)) {
0228       HypeHelpers<Real_v, false, true>::GetPointOfIntersectionWithHyperbolicSurface(hype, point, direction, dist);
0229       vecCore__MaskedAssignFunc(dist, dist < Real_v(0.), InfinityLength<Real_v>());
0230       vecCore__MaskedAssignFunc(distance, !done, Min(distance, dist));
0231     }
0232   }
0233
0234   template <class Real_v>
0235   VECCORE_ATT_HOST_DEVICE
0236   static void SafetyToIn(UnplacedStruct_t const &hype, Vector3D<Real_v> const &point, Real_v &safety)
0237   {
0238
0239     using Bool_v = typename vecCore::Mask_v<Real_v>;
0240
0241     Real_v absZ = Abs(point.z());
0242     Real_v r2   = point.Perp2();
0243     Real_v r    = Sqrt(r2);
0244
0245     Bool_v done(false);
0246
0247     // New Simple Algo
0248     safety = 0.;
0249     // If point is inside then safety should be -1.
0250     Bool_v compIn(false), compOut(false);
0251     GenericKernelForContainsAndInside<Real_v, Bool_v, true>(hype, point, compIn, compOut);
0252     done = (!compIn && !compOut);
0253     if (vecCore::MaskFull(done)) return;
0254
0255     vecCore__MaskedAssignFunc(safety, compIn, Real_v(-1.0));
0256     done |= compIn;
0257     if (vecCore::MaskFull(done)) return;
0258
0259     Bool_v cond(false);
0260     Real_v sigz = absZ - hype.fDz;
0261     cond        = (sigz > kHalfTolerance) && (r < hype.fEndOuterRadius) && (r > hype.fEndInnerRadius);
0262     vecCore::MaskedAssign(safety, !done && cond, sigz);
0263     done |= cond;
0264     if (vecCore::MaskFull(done)) return;
0265
0266     cond = (sigz > kHalfTolerance) && (r > hype.fEndOuterRadius);
0267     vecCore__MaskedAssignFunc(safety, !done && cond,
0268                               Sqrt((r - hype.fEndOuterRadius) * (r - hype.fEndOuterRadius) + (sigz) * (sigz)));
0269     done |= cond;
0270     if (vecCore::MaskFull(done)) return;
0271
0272     cond = (sigz > kHalfTolerance) && (r < hype.fEndInnerRadius);
0273     vecCore__MaskedAssignFunc(safety, !done && cond,
0274                               Sqrt((r - hype.fEndInnerRadius) * (r - hype.fEndInnerRadius) + (sigz) * (sigz)));
0275     done |= cond;
0276     if (vecCore::MaskFull(done)) return;
0277
0278     cond =
0279         (r2 > ((hype.fRmax2 + hype.fTOut2 * absZ * absZ) + kHalfTolerance)) && (absZ > Real_v(0.)) && (absZ < hype.fDz);
0280     vecCore__MaskedAssignFunc(safety, !done && cond,
0281                               HypeUtilities::ApproxDistOutside<Real_v>(r, absZ, hype.fRmax, hype.fTOut));
0282     done |= cond;
0283     if (vecCore::MaskFull(done)) return;
0284
0285     vecCore__MaskedAssignFunc(safety,
0286                               !done && (r2 < ((hype.fRmin2 + hype.fTIn2 * absZ * absZ) - kHalfTolerance)) &&
0287                                   (absZ > Real_v(0.)) && (absZ < hype.fDz),
0288                               HypeUtilities::ApproxDistInside<Real_v>(r, absZ, hype.fRmin, hype.fTIn2));
0289   }
0290
0291   template <class Real_v>
0292   VECCORE_ATT_HOST_DEVICE
0293   static void SafetyToOut(UnplacedStruct_t const &hype, Vector3D<Real_v> const &point, Real_v &safety)
0294   {
0295     using namespace ::vecgeom::HypeTypes;
0296     using Bool_v = typename vecCore::Mask_v<Real_v>;
0297     safety       = 0.;
0298     Real_v r     = Sqrt(point.x() * point.x() + point.y() * point.y());
0299     Real_v absZ  = Abs(point.z());
0300     Bool_v inside(false), outside(false);
0301     Bool_v done(false);
0302
0303     Real_v distZ(0.), distInner(0.), distOuter(0.);
0304     safety = 0.;
0305     GenericKernelForContainsAndInside<Real_v, Bool_v, true>(hype, point, inside, outside);
0306     done = (!inside && !outside);
0307     if (vecCore::MaskFull(done)) return;
0308
0309     vecCore__MaskedAssignFunc(safety, outside, Real_v(-1.0));
0310     done |= outside;
0311     if (vecCore::MaskFull(done)) return;
0312
0313     vecCore__MaskedAssignFunc(distZ, !done && inside, Abs(Abs(point.z()) - hype.fDz));
0314     if (checkInnerSurfaceTreatment<hypeTypeT>(hype) && hype.fStIn) {
0315       vecCore__MaskedAssignFunc(distInner, !done && inside,
0316                                 HypeUtilities::ApproxDistOutside<Real_v>(r, absZ, hype.fRmin, hype.fTIn));
0317     }
0318
0319     if (checkInnerSurfaceTreatment<hypeTypeT>(hype) && !hype.fStIn) {
0320       vecCore__MaskedAssignFunc(distInner, !done && inside, (r - hype.fRmin));
0321     }
0322
0323     if (!checkInnerSurfaceTreatment<hypeTypeT>(hype) && !hype.fStIn) {
0324       vecCore__MaskedAssignFunc(distInner, !done && inside, InfinityLength<Real_v>());
0325     }
0326
0327     vecCore__MaskedAssignFunc(distOuter, !done && inside,
0328                               HypeUtilities::ApproxDistInside<Real_v>(r, absZ, hype.fRmax, hype.fTOut2));
0329     safety = Min(distInner, distOuter);
0330     safety = Min(safety, distZ);
0331   }
0332 };
0333 } // namespace VECGEOM_IMPL_NAMESPACE
0334 } // namespace vecgeom
0335
0336 #endif // VECGEOM_VOLUMES_KERNEL_HYPEIMPLEMENTATION_H_