Back to home page

EIC code displayed by LXR
 
 

    

File indexing completed on 2025-01-18 10:14:02 

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 /// \brief This file implements the algorithms for Orb
0006 /// \file volumes/kernel/orbImplementation.h
0007 /// \author Raman Sehgal
0008 
0009 /// History notes:
0010 /// 2014 - 2015: original development (abstracted kernels); Raman Sehgal
0011 /// July 2016: revision + moving to new backend structure (Raman Sehgal)
0012 
0013 #ifndef VECGEOM_VOLUMES_KERNEL_ORBIMPLEMENTATION_H_
0014 #define VECGEOM_VOLUMES_KERNEL_ORBIMPLEMENTATION_H_
0015 
0016 #include "VecGeom/base/Vector3D.h"
0017 #include "VecGeom/volumes/OrbStruct.h"
0018 #include "VecGeom/volumes/kernel/GenericKernels.h"
0019 #include <VecCore/VecCore>
0020 
0021 #include <cstdio>
0022 
0023 namespace vecgeom {
0024 
0025 VECGEOM_DEVICE_FORWARD_DECLARE(struct OrbImplementation;);
0026 VECGEOM_DEVICE_DECLARE_CONV(struct, OrbImplementation);
0027 
0028 inline namespace VECGEOM_IMPL_NAMESPACE {
0029 
0030 class PlacedOrb;
0031 template <typename T>
0032 struct OrbStruct;
0033 class UnplacedOrb;
0034 
0035 struct OrbImplementation {
0036 
0037   using PlacedShape_t    = PlacedOrb;
0038   using UnplacedStruct_t = OrbStruct<Precision>;
0039   using UnplacedVolume_t = UnplacedOrb;
0040 
0041   VECCORE_ATT_HOST_DEVICE
0042   static void PrintType()
0043   {
0044     //  printf("SpecializedOrb<%i, %i>", transCodeT, rotCodeT);
0045   }
0046 
0047   template <typename Stream>
0048   static void PrintType(Stream &st, int transCodeT = translation::kGeneric, int rotCodeT = rotation::kGeneric)
0049   {
0050     st << "SpecializedOrb<" << transCodeT << "," << rotCodeT << ">";
0051   }
0052 
0053   template <typename Stream>
0054   static void PrintImplementationType(Stream &st)
0055   {
0056     (void)st;
0057     // st << "OrbImplementation<" << transCodeT << "," << rotCodeT << ">";
0058   }
0059 
0060   template <typename Stream>
0061   static void PrintUnplacedType(Stream &st)
0062   {
0063     (void)st;
0064     // TODO: this is wrong
0065     // st << "UnplacedOrb";
0066   }
0067 
0068   template <typename Real_v, typename Bool_v>
0069   VECGEOM_FORCE_INLINE
0070   VECCORE_ATT_HOST_DEVICE
0071   static void Contains(UnplacedStruct_t const &orb, Vector3D<Real_v> const &point, Bool_v &inside)
0072   {
0073     Bool_v unused, outside;
0074     GenericKernelForContainsAndInside<Real_v, Bool_v, false>(orb, point, unused, outside);
0075     inside = !outside;
0076   }
0077 
0078   // BIG QUESTION: DO WE WANT TO GIVE ALL 3 TEMPLATE PARAMETERS
0079   // -- OR -- DO WE WANT TO DEDUCE Bool_v, Index_t from Real_v???
0080   template <typename Real_v, typename Inside_t>
0081   VECGEOM_FORCE_INLINE
0082   VECCORE_ATT_HOST_DEVICE
0083   static void Inside(UnplacedStruct_t const &orb, Vector3D<Real_v> const &point, Inside_t &inside)
0084   {
0085 
0086     using Bool_v       = vecCore::Mask_v<Real_v>;
0087     using InsideBool_v = vecCore::Mask_v<Inside_t>;
0088     Bool_v completelyinside, completelyoutside;
0089     GenericKernelForContainsAndInside<Real_v, Bool_v, true>(orb, point, completelyinside, completelyoutside);
0090     inside = EInside::kSurface;
0091     vecCore::MaskedAssign(inside, (InsideBool_v)completelyoutside, Inside_t(EInside::kOutside));
0092     vecCore::MaskedAssign(inside, (InsideBool_v)completelyinside, Inside_t(EInside::kInside));
0093   }
0094 
0095   template <typename Real_v, typename Bool_v, bool ForInside>
0096   VECGEOM_FORCE_INLINE
0097   VECCORE_ATT_HOST_DEVICE
0098   static void GenericKernelForContainsAndInside(UnplacedStruct_t const &orb, Vector3D<Real_v> const &localPoint,
0099                                                 Bool_v &completelyinside, Bool_v &completelyoutside)
0100   {
0101     Precision fR = orb.fR;
0102     Real_v rad2  = localPoint.Mag2();
0103     Real_v tolR  = fR - Real_v(kTolerance);
0104     if (ForInside) completelyinside = (rad2 <= tolR * tolR);
0105     tolR              = fR + Real_v(kTolerance);
0106     completelyoutside = (rad2 >= tolR * tolR);
0107     return;
0108   }
0109 
0110   template <typename Real_v>
0111   VECGEOM_FORCE_INLINE
0112   VECCORE_ATT_HOST_DEVICE
0113   static void DistanceToIn(UnplacedStruct_t const &orb, Vector3D<Real_v> const &point,
0114                            Vector3D<Real_v> const &direction, Real_v const & /*stepMax*/, Real_v &distance)
0115   {
0116     using Bool_v         = vecCore::Mask_v<Real_v>;
0117     distance             = kInfLength;
0118     Real_v rad           = point.Mag();
0119     Bool_v isPointInside = (rad < Real_v(orb.fR - kTolerance));
0120     vecCore__MaskedAssignFunc(distance, isPointInside, Real_v(-1.));
0121     Bool_v done = isPointInside;
0122     if (vecCore::MaskFull(done)) return;
0123 
0124     Real_v pDotV3D          = point.Dot(direction);
0125     Bool_v isPointOnSurface = (rad >= Real_v(orb.fR - kTolerance)) && (rad <= Real_v(orb.fR + kTolerance));
0126     Bool_v cond             = (isPointOnSurface && (pDotV3D < Real_v(0.)));
0127     vecCore__MaskedAssignFunc(distance, !done && cond, Real_v(0.));
0128     done |= cond;
0129     if (vecCore::MaskFull(done)) return;
0130     Real_v dist(kInfLength);
0131     vecCore::MaskedAssign(
0132         distance, !done && DetectIntersectionAndCalculateDistance<Real_v, true>(orb, point, direction, dist), dist);
0133   }
0134 
0135   template <typename Real_v>
0136   VECGEOM_FORCE_INLINE
0137   VECCORE_ATT_HOST_DEVICE
0138   static void DistanceToOut(UnplacedStruct_t const &orb, Vector3D<Real_v> const &point,
0139                             Vector3D<Real_v> const &direction, Real_v const & /* stepMax */, Real_v &distance)
0140   {
0141     using Bool_v = vecCore::Mask_v<Real_v>;
0142 
0143     distance = kInfLength;
0144 
0145     Real_v rad            = point.Mag();
0146     Bool_v isPointOutside = (rad > Real_v(orb.fR + kTolerance));
0147     vecCore__MaskedAssignFunc(distance, isPointOutside, Real_v(-1.));
0148     Bool_v done = isPointOutside;
0149     if (vecCore::MaskFull(done)) return;
0150 
0151     Real_v pDotV3D          = point.Dot(direction);
0152     Bool_v isPointOnSurface = (rad >= Real_v(orb.fR - kTolerance)) && (rad <= Real_v(orb.fR + kTolerance));
0153     Bool_v cond             = (isPointOnSurface && (pDotV3D > Real_v(0.)));
0154     vecCore__MaskedAssignFunc(distance, !done && cond, Real_v(0.));
0155     done |= cond;
0156     if (vecCore::MaskFull(done)) return;
0157     Real_v dist(kInfLength);
0158     vecCore::MaskedAssign(
0159         distance, !done && DetectIntersectionAndCalculateDistance<Real_v, false>(orb, point, direction, dist), dist);
0160 
0161     return;
0162   }
0163 
0164   template <typename Real_v>
0165   VECGEOM_FORCE_INLINE
0166   VECCORE_ATT_HOST_DEVICE
0167   static void SafetyToIn(UnplacedStruct_t const &orb, Vector3D<Real_v> const &point, Real_v &safety)
0168   {
0169     using Bool_v         = vecCore::Mask_v<Real_v>;
0170     Real_v rad           = point.Mag();
0171     safety               = rad - Real_v(orb.fR);
0172     Bool_v isPointInside = (rad < Real_v(orb.fR - kTolerance));
0173     vecCore__MaskedAssignFunc(safety, isPointInside, Real_v(-1.));
0174     if (vecCore::MaskFull(isPointInside)) return;
0175 
0176     Bool_v isPointOnSurface = (rad > Real_v(orb.fR - kTolerance)) && (rad < Real_v(orb.fR + kTolerance));
0177     vecCore__MaskedAssignFunc(safety, isPointOnSurface, Real_v(0.));
0178   }
0179 
0180   template <typename Real_v>
0181   VECGEOM_FORCE_INLINE
0182   VECCORE_ATT_HOST_DEVICE
0183   static void SafetyToOut(UnplacedStruct_t const &orb, Vector3D<Real_v> const &point, Real_v &safety)
0184   {
0185     using Bool_v = vecCore::Mask_v<Real_v>;
0186 
0187     Real_v rad = point.Mag();
0188     safety     = Real_v(orb.fR) - rad;
0189 
0190     Bool_v isPointOutside = (rad > Real_v(orb.fR + kTolerance));
0191     vecCore__MaskedAssignFunc(safety, isPointOutside, Real_v(-1.));
0192     if (vecCore::MaskFull(isPointOutside)) return;
0193 
0194     Bool_v isPointOnSurface = (rad > Real_v(orb.fR - kTolerance)) && (rad < Real_v(orb.fR + kTolerance));
0195     vecCore__MaskedAssignFunc(safety, isPointOnSurface, Real_v(0.));
0196   }
0197 
0198   template <typename Real_v, bool ForDistanceToIn>
0199   VECGEOM_FORCE_INLINE
0200   VECCORE_ATT_HOST_DEVICE
0201   static typename vecCore::Mask_v<Real_v> DetectIntersectionAndCalculateDistance(UnplacedStruct_t const &orb,
0202                                                                                  Vector3D<Real_v> const &point,
0203                                                                                  Vector3D<Real_v> const &direction,
0204                                                                                  Real_v &distance)
0205   {
0206 
0207     using Bool_v   = vecCore::Mask_v<Real_v>;
0208     Real_v rad2    = point.Mag2();
0209     Real_v pDotV3D = point.Dot(direction);
0210     Precision fR   = orb.fR;
0211     Real_v c       = rad2 - fR * fR;
0212     Real_v d2      = (pDotV3D * pDotV3D - c);
0213 
0214     if (ForDistanceToIn) {
0215       Bool_v cond = ((d2 >= Real_v(0.)) && (pDotV3D <= Real_v(0.)));
0216       vecCore__MaskedAssignFunc(distance, cond, (-pDotV3D - Sqrt(vecCore::math::Abs(d2))));
0217       return cond;
0218     } else {
0219       vecCore__MaskedAssignFunc(distance, (d2 >= Real_v(0.)), (-pDotV3D + Sqrt(vecCore::math::Abs(d2))));
0220       return (d2 >= Real_v(0.));
0221     }
0222   }
0223 
0224   template <typename Real_v>
0225   VECGEOM_FORCE_INLINE
0226   VECCORE_ATT_HOST_DEVICE
0227   static Vector3D<Real_v> NormalKernel(UnplacedStruct_t const &orb, Vector3D<Real_v> const &point,
0228                                        typename vecCore::Mask_v<Real_v> &valid)
0229   {
0230     Real_v rad2             = point.Mag2();
0231     Real_v invRadius        = Real_v(1.) / Sqrt(rad2);
0232     Vector3D<Real_v> normal = point * invRadius;
0233 
0234     Real_v tolRMaxO = orb.fR + kTolerance;
0235     Real_v tolRMaxI = orb.fR - kTolerance;
0236 
0237     // Check radial surface
0238     valid = ((rad2 <= tolRMaxO * tolRMaxO) && (rad2 >= tolRMaxI * tolRMaxI)); // means we are on surface
0239     return normal;
0240   }
0241 };
0242 } // namespace VECGEOM_IMPL_NAMESPACE
0243 } // namespace vecgeom
0244 
0245 #endif // VECGEOM_VOLUMES_KERNEL_orbIMPLEMENTATION_H_
[ Source navigation ] [ Diff markup ] [ Identifier search ] [ general search ]

This page was automatically generated by the 2.3.7 LXR engine.

The LXR team
Valid CSS 2.1! Valid HTML 4.01!