Back to home page

EIC code displayed by LXR
 
 

    

File indexing completed on 2026-05-09 08:59:36 

0001 
0002 /// @file SphereImplementation.h
0003 /// @author Raman Sehgal (raman.sehgal@cern.ch)
0004 
0005 #ifndef VECGEOM_VOLUMES_KERNEL_SPHEREIMPLEMENTATION_H_
0006 #define VECGEOM_VOLUMES_KERNEL_SPHEREIMPLEMENTATION_H_
0007 
0008 #include "VecGeom/base/Vector3D.h"
0009 #include "VecGeom/volumes/SphereStruct.h"
0010 #include "VecGeom/volumes/kernel/GenericKernels.h"
0011 #include <VecCore/VecCore>
0012 #include "VecGeom/volumes/kernel/OrbImplementation.h"
0013 #include "VecGeom/volumes/SphereUtilities.h"
0014 
0015 #include <cstdio>
0016 
0017 namespace vecgeom {
0018 
0019 VECGEOM_DEVICE_FORWARD_DECLARE(struct SphereImplementation;);
0020 VECGEOM_DEVICE_DECLARE_CONV(struct, SphereImplementation);
0021 
0022 inline namespace VECGEOM_IMPL_NAMESPACE {
0023 
0024 class PlacedSphere;
0025 template <typename T>
0026 struct SphereStruct;
0027 class UnplacedSphere;
0028 
0029 struct SphereImplementation {
0030 
0031   using PlacedShape_t    = PlacedSphere;
0032   using UnplacedStruct_t = SphereStruct<Precision>;
0033   using UnplacedVolume_t = UnplacedSphere;
0034 
0035   template <typename Real_v, typename Bool_v>
0036   VECGEOM_FORCE_INLINE VECCORE_ATT_HOST_DEVICE static void Contains(UnplacedStruct_t const &sphere,
0037                                                                     Vector3D<Real_v> const &point, Bool_v &inside)
0038   {
0039     Bool_v unused(false), outside(false);
0040     GenericKernelForContainsAndInside<Real_v, Bool_v, false>(sphere, point, unused, outside);
0041     inside = !outside;
0042   }
0043 
0044   // BIG QUESTION: DO WE WANT TO GIVE ALL 3 TEMPLATE PARAMETERS
0045   // -- OR -- DO WE WANT TO DEDUCE Bool_v, Index_t from Real_v???
0046   template <typename Real_v, typename Inside_t>
0047   VECGEOM_FORCE_INLINE VECCORE_ATT_HOST_DEVICE static void Inside(UnplacedStruct_t const &sphere,
0048                                                                   Vector3D<Real_v> const &point, Inside_t &inside)
0049   {
0050 
0051     using Bool_v       = vecCore::Mask_v<Real_v>;
0052     using InsideBool_v = vecCore::Mask_v<Inside_t>;
0053     Bool_v completelyinside, completelyoutside;
0054     GenericKernelForContainsAndInside<Real_v, Bool_v, true>(sphere, point, completelyinside, completelyoutside);
0055     inside = EInside::kSurface;
0056     vecCore::MaskedAssign(inside, (InsideBool_v)completelyoutside, Inside_t(EInside::kOutside));
0057     vecCore::MaskedAssign(inside, (InsideBool_v)completelyinside, Inside_t(EInside::kInside));
0058   }
0059 
0060   template <typename Real_v, typename Bool_v, bool ForInside>
0061   VECGEOM_FORCE_INLINE VECCORE_ATT_HOST_DEVICE static void GenericKernelForContainsAndInside(
0062       UnplacedStruct_t const &sphere, Vector3D<Real_v> const &localPoint, Bool_v &completelyinside,
0063       Bool_v &completelyoutside)
0064   {
0065     Real_v rad2 = localPoint.Mag2();
0066 
0067     // Check radial surfaces
0068     // Radial check for GenericKernel Start
0069     if (sphere.fRmin)
0070       completelyinside =
0071           rad2 <= MakeMinusTolerantSquare<true>(sphere.fRmax) && rad2 >= MakePlusTolerantSquare<true>(sphere.fRmin);
0072     else
0073       completelyinside = rad2 <= MakeMinusTolerantSquare<true>(sphere.fRmax);
0074 
0075     if (sphere.fRmin)
0076       completelyoutside =
0077           rad2 >= MakePlusTolerantSquare<true>(sphere.fRmax) || rad2 <= MakeMinusTolerantSquare<true>(sphere.fRmin);
0078     else
0079       completelyoutside = rad2 >= MakePlusTolerantSquare<true>(sphere.fRmax);
0080 
0081     // Phi boundaries  : Do not check if it has no phi boundary!
0082     if (!sphere.fFullPhiSphere) {
0083 
0084       Bool_v completelyoutsidephi(false);
0085       Bool_v completelyinsidephi(false);
0086       sphere.fPhiWedge.GenericKernelForContainsAndInside<Real_v, ForInside>(localPoint, completelyinsidephi,
0087                                                                             completelyoutsidephi);
0088       completelyoutside |= completelyoutsidephi;
0089 
0090       if (ForInside) completelyinside &= completelyinsidephi;
0091     }
0092     // Phi Check for GenericKernel Over
0093 
0094     // Theta bondaries
0095     if (!sphere.fFullThetaSphere) {
0096 
0097       Bool_v completelyoutsidetheta(false);
0098       Bool_v completelyinsidetheta(false);
0099       sphere.fThetaCone.GenericKernelForContainsAndInside<Real_v, ForInside>(localPoint, completelyinsidetheta,
0100                                                                              completelyoutsidetheta);
0101       completelyoutside |= completelyoutsidetheta;
0102 
0103       if (ForInside) completelyinside &= completelyinsidetheta;
0104     }
0105     return;
0106   }
0107 
0108   template <typename Real_v>
0109   VECGEOM_FORCE_INLINE VECCORE_ATT_HOST_DEVICE static void DistanceToIn(UnplacedStruct_t const &sphere,
0110                                                                         Vector3D<Real_v> const &point,
0111                                                                         Vector3D<Real_v> const &direction,
0112                                                                         Real_v const & /* stepMax */, Real_v &distance)
0113   {
0114     using Bool_v = vecCore::Mask_v<Real_v>;
0115     distance     = kInfLength;
0116 
0117     Bool_v done(false);
0118 
0119     bool fullPhiSphere   = sphere.fFullPhiSphere;
0120     bool fullThetaSphere = sphere.fFullThetaSphere;
0121 
0122     Vector3D<Real_v> tmpPt;
0123     // General Precalcs
0124     Real_v rad2    = point.Mag2();
0125     Real_v pDotV3d = point.Dot(direction);
0126 
0127     Real_v c = rad2 - sphere.fRmax * sphere.fRmax;
0128 
0129     Bool_v cond = SphereUtilities::IsCompletelyInside<Real_v>(sphere, point);
0130     vecCore__MaskedAssignFunc(distance, cond, Real_v(-1.0));
0131     done |= cond;
0132     if (vecCore::MaskFull(done)) return;
0133 
0134     cond = SphereUtilities::IsPointOnSurfaceAndMovingOut<Real_v, false>(sphere, point, direction);
0135     vecCore__MaskedAssignFunc(distance, !done && cond, Real_v(0.0));
0136     done |= cond;
0137     if (vecCore::MaskFull(done)) return;
0138 
0139     Real_v sd1(kInfLength);
0140     Real_v sd2(kInfLength);
0141     Real_v d2 = (pDotV3d * pDotV3d - c);
0142     cond      = (d2 < Real_v(0.) || ((c > Real_v(0.)) && (pDotV3d > Real_v(0.))));
0143     done |= cond;
0144     if (vecCore::MaskFull(done)) return; // Returning in case of no intersection with outer shell
0145 
0146     // Note: Abs(d2) was introduced to avoid Sqrt(negative) in other lanes than the ones satisfying d2>=0.
0147     vecCore__MaskedAssignFunc(sd1, d2 >= Real_v(0.), (-pDotV3d - Sqrt(Abs(d2))));
0148 
0149     Real_v outerDist(kInfLength);
0150     Real_v innerDist(kInfLength);
0151 
0152     if (sphere.fFullSphere) {
0153       vecCore::MaskedAssign(outerDist, !done && (sd1 >= Real_v(-kTolerance)), sd1);
0154     } else {
0155       tmpPt = point + sd1 * direction;
0156       vecCore::MaskedAssign(outerDist,
0157                             !done && (sd1 >= Real_v(-kTolerance) && sd1 < Real_v(kInfLength)) &&
0158                                 sphere.fPhiWedge.Inside<Real_v, Inside_t>(tmpPt) != EInside::kOutside &&
0159                                 sphere.fThetaCone.Inside<Real_v, Inside_t>(tmpPt) != EInside::kOutside,
0160                             sd1);
0161     }
0162 
0163     if (sphere.fRmin) {
0164       c  = rad2 - sphere.fRmin * sphere.fRmin;
0165       d2 = pDotV3d * pDotV3d - c;
0166       // Note: Abs(d2) was introduced to avoid Sqrt(negative) in other lanes than the ones satisfying d2>=0.
0167       vecCore__MaskedAssignFunc(sd2, d2 >= Real_v(0.), (-pDotV3d + Sqrt(Abs(d2))));
0168 
0169       if (sphere.fFullSphere) {
0170         vecCore::MaskedAssign(innerDist, !done && (sd2 >= Real_v(-kTolerance)), sd2);
0171       } else {
0172         //   std::cerr<<" ---- Called by InnerRad ---- " << std::endl;
0173         tmpPt = point + sd2 * direction;
0174         vecCore::MaskedAssign(innerDist,
0175                               !done && (sd2 >= Real_v(-kTolerance) && sd2 < Real_v(kInfLength)) &&
0176                                   sphere.fPhiWedge.Inside<Real_v, Inside_t>(tmpPt) != EInside::kOutside &&
0177                                   sphere.fThetaCone.Inside<Real_v, Inside_t>(tmpPt) != EInside::kOutside,
0178                               sd2);
0179       }
0180     }
0181 
0182     distance = Min(outerDist, innerDist);
0183 
0184     if (!fullPhiSphere) {
0185       GetMinDistFromPhi<Real_v, true>(sphere, point, direction, done, distance);
0186     }
0187 
0188     Real_v distThetaMin(kInfLength);
0189 
0190     if (!fullThetaSphere) {
0191       Bool_v intsect1(false);
0192       Bool_v intsect2(false);
0193       Real_v distTheta1(kInfLength), distCone1(kInfLength);
0194       Real_v distTheta2(kInfLength), distCone2(kInfLength);
0195 
0196       sphere.fThetaCone.DistanceToIn<Real_v>(point, direction, distTheta1, distTheta2, intsect1,
0197                                              intsect2); //,cone1IntSecPt, cone2IntSecPt);
0198       Vector3D<Real_v> coneIntSecPt1 = point + distTheta1 * direction;
0199       if (vecCore::MaskFull(intsect1)) distCone1 = coneIntSecPt1.Mag2(); // avoid FPE due to kInfLength * kInfLength
0200 
0201       Vector3D<Real_v> coneIntSecPt2 = point + distTheta2 * direction;
0202       if (vecCore::MaskFull(intsect2)) distCone2 = coneIntSecPt2.Mag2(); // avoid FPE due to kInfLength * kInfLength
0203 
0204       Bool_v isValidCone1 =
0205           (distCone1 >= sphere.fRmin * sphere.fRmin && distCone1 <= sphere.fRmax * sphere.fRmax) && intsect1;
0206       Bool_v isValidCone2 =
0207           (distCone2 >= sphere.fRmin * sphere.fRmin && distCone2 <= sphere.fRmax * sphere.fRmax) && intsect2;
0208 
0209       if (!fullPhiSphere) {
0210         isValidCone1 &= sphere.fPhiWedge.Contains<Real_v>(coneIntSecPt1);
0211         isValidCone2 &= sphere.fPhiWedge.Contains<Real_v>(coneIntSecPt2);
0212       }
0213       vecCore::MaskedAssign(distThetaMin, (!done && isValidCone2 && !isValidCone1), distTheta2);
0214       vecCore::MaskedAssign(distThetaMin, (!done && isValidCone1 && !isValidCone2), distTheta1);
0215       vecCore__MaskedAssignFunc(distThetaMin, (!done && isValidCone1 && isValidCone2), Min(distTheta1, distTheta2));
0216     }
0217 
0218     distance = Min(distThetaMin, distance);
0219   }
0220 
0221   template <typename Real_v>
0222   VECGEOM_FORCE_INLINE VECCORE_ATT_HOST_DEVICE static void DistanceToOut(UnplacedStruct_t const &sphere,
0223                                                                          Vector3D<Real_v> const &point,
0224                                                                          Vector3D<Real_v> const &direction,
0225                                                                          Real_v const & /* stepMax */, Real_v &distance)
0226   {
0227 
0228     using Bool_v = typename vecCore::Mask_v<Real_v>;
0229 
0230     distance = kInfLength;
0231     Bool_v done(false);
0232 
0233     Real_v snxt(kInfLength);
0234 
0235     // Intersection point
0236     Vector3D<Real_v> intSecPt;
0237     Real_v d2(0.);
0238 
0239     Real_v pDotV3d = point.Dot(direction);
0240 
0241     Real_v rad2 = point.Mag2();
0242     Real_v c    = rad2 - sphere.fRmax * sphere.fRmax;
0243 
0244     Real_v sd1(kInfLength);
0245     Real_v sd2(kInfLength);
0246 
0247     Bool_v cond = SphereUtilities::IsCompletelyOutside<Real_v>(sphere, point);
0248     vecCore__MaskedAssignFunc(distance, cond, Real_v(-1.0));
0249 
0250     done |= cond;
0251     if (vecCore::MaskFull(done)) return;
0252 
0253     cond = SphereUtilities::IsPointOnSurfaceAndMovingOut<Real_v, true>(sphere, point, direction);
0254     vecCore__MaskedAssignFunc(distance, !done && cond, Real_v(0.0));
0255     done |= cond;
0256     if (vecCore::MaskFull(done)) return;
0257 
0258     // Note: Abs(d2) was introduced to avoid Sqrt(negative) in other lanes than the ones satisfying d2>=0.
0259     d2 = (pDotV3d * pDotV3d - c);
0260     vecCore__MaskedAssignFunc(sd1, (!done && (d2 >= Real_v(0.))), (-pDotV3d + Sqrt(Abs(d2))));
0261 
0262     if (sphere.fRmin) {
0263       c  = rad2 - sphere.fRmin * sphere.fRmin;
0264       d2 = (pDotV3d * pDotV3d - c);
0265       vecCore__MaskedAssignFunc(sd2, (!done && (d2 >= Real_v(0.)) && (pDotV3d < Real_v(0.))),
0266                                 (-pDotV3d - Sqrt(Abs(d2))));
0267     }
0268 
0269     snxt = Min(sd1, sd2);
0270 
0271     Bool_v condSemi = (Bool_v(sphere.fSTheta == 0. && sphere.eTheta == kPi / 2.) && direction.z() >= Real_v(0.)) ||
0272                       (Bool_v(sphere.fSTheta == kPi / 2. && sphere.eTheta == kPi) && direction.z() <= Real_v(0.));
0273     vecCore::MaskedAssign(distance, !done && condSemi, snxt);
0274     done |= condSemi;
0275     if (vecCore::MaskFull(done)) return;
0276 
0277     Real_v distThetaMin(kInfLength);
0278     Real_v distPhiMin(kInfLength);
0279 
0280     if (!sphere.fFullThetaSphere) {
0281       Bool_v intsect1(false);
0282       Bool_v intsect2(false);
0283       Real_v distTheta1(kInfLength);
0284       Real_v distTheta2(kInfLength);
0285       sphere.fThetaCone.DistanceToOut<Real_v>(point, direction, distTheta1, distTheta2, intsect1, intsect2);
0286       vecCore::MaskedAssign(distThetaMin, (intsect2 && !intsect1), distTheta2);
0287       vecCore::MaskedAssign(distThetaMin, (!intsect2 && intsect1), distTheta1);
0288       vecCore__MaskedAssignFunc(distThetaMin, (intsect2 && intsect1), Min(distTheta1, distTheta2));
0289     }
0290 
0291     distance = Min(distThetaMin, snxt);
0292 
0293     if (!sphere.fFullPhiSphere) {
0294       if (sphere.fDPhi <= kPi) {
0295         Real_v distPhi1;
0296         Real_v distPhi2;
0297         sphere.fPhiWedge.DistanceToOut<Real_v>(point, direction, distPhi1, distPhi2);
0298         distPhiMin = Min(distPhi1, distPhi2);
0299         distance   = Min(distPhiMin, distance);
0300       } else {
0301         GetMinDistFromPhi<Real_v, false>(sphere, point, direction, done, distance);
0302       }
0303     }
0304   }
0305 
0306   template <typename Real_v, bool DistToIn>
0307   VECCORE_ATT_HOST_DEVICE static void GetMinDistFromPhi(UnplacedStruct_t const &sphere,
0308                                                         Vector3D<Real_v> const &localPoint,
0309                                                         Vector3D<Real_v> const &localDir,
0310                                                         typename vecCore::Mask_v<Real_v> &done, Real_v &distance)
0311   {
0312     using Bool_v = typename vecCore::Mask_v<Real_v>;
0313     Real_v distPhi1(kInfLength);
0314     Real_v distPhi2(kInfLength);
0315     Real_v dist(kInfLength);
0316 
0317     if (DistToIn)
0318       sphere.fPhiWedge.DistanceToIn<Real_v>(localPoint, localDir, distPhi1, distPhi2);
0319     else
0320       sphere.fPhiWedge.DistanceToOut<Real_v>(localPoint, localDir, distPhi1, distPhi2);
0321 
0322     Bool_v containsCond1(false), containsCond2(false);
0323     Bool_v tempCond(false);
0324     // Min Face
0325     dist = Min(distPhi1, distPhi2);
0326     if (dist < kInfLength) {
0327       Vector3D<Real_v> tmpPt = localPoint + dist * localDir;
0328       Real_v rad2            = tmpPt.Mag2();
0329 
0330       tempCond = ((dist == distPhi1) && sphere.fPhiWedge.IsOnSurfaceGeneric<Real_v, true>(tmpPt)) ||
0331                  ((dist == distPhi2) && sphere.fPhiWedge.IsOnSurfaceGeneric<Real_v, false>(tmpPt));
0332 
0333       containsCond1 = tempCond && (rad2 > sphere.fRmin * sphere.fRmin) && (rad2 < sphere.fRmax * sphere.fRmax) &&
0334                       sphere.fThetaCone.Contains<Real_v>(tmpPt);
0335     }
0336 
0337     vecCore__MaskedAssignFunc(distance, !done && containsCond1, Min(dist, distance));
0338 
0339     // Max Face
0340     dist = Max(distPhi1, distPhi2);
0341     if (dist < kInfLength) {
0342       Vector3D<Real_v> tmpPt = localPoint + dist * localDir;
0343 
0344       Real_v rad2 = tmpPt.Mag2();
0345       tempCond    = ((dist == distPhi1) && sphere.fPhiWedge.IsOnSurfaceGeneric<Real_v, true>(tmpPt)) ||
0346                  ((dist == distPhi2) && sphere.fPhiWedge.IsOnSurfaceGeneric<Real_v, false>(tmpPt));
0347 
0348       containsCond2 = tempCond && (rad2 > sphere.fRmin * sphere.fRmin) && (rad2 < sphere.fRmax * sphere.fRmax) &&
0349                       sphere.fThetaCone.Contains<Real_v>(tmpPt);
0350     }
0351     vecCore__MaskedAssignFunc(distance, ((!done) && (!containsCond1) && containsCond2), Min(dist, distance));
0352   }
0353 
0354   template <typename Real_v>
0355   VECGEOM_FORCE_INLINE VECCORE_ATT_HOST_DEVICE static void SafetyToIn(UnplacedStruct_t const &sphere,
0356                                                                       Vector3D<Real_v> const &point, Real_v &safety)
0357   {
0358     using Bool_v = vecCore::Mask_v<Real_v>;
0359     Bool_v done(false);
0360 
0361     // General Precalcs
0362     Real_v rad = point.Mag();
0363 
0364     Real_v safeRMin(0.);
0365     Real_v safeRMax(0.);
0366 
0367     Bool_v completelyinside(false), completelyoutside(false);
0368     GenericKernelForContainsAndInside<Real_v, Bool_v, true>(sphere, point, completelyinside, completelyoutside);
0369 
0370     vecCore__MaskedAssignFunc(safety, completelyinside, Real_v(-1.0));
0371     done |= completelyinside;
0372     if (vecCore::MaskFull(done)) return;
0373 
0374     Bool_v isOnSurface = !completelyinside && !completelyoutside;
0375     vecCore__MaskedAssignFunc(safety, !done && isOnSurface, Real_v(0.0));
0376     done |= isOnSurface;
0377     if (vecCore::MaskFull(done)) return;
0378 
0379     if (sphere.fRmin) {
0380       safeRMin = sphere.fRmin - rad;
0381       safeRMax = rad - sphere.fRmax;
0382       safety   = vecCore::Blend(!done && (safeRMin > safeRMax), safeRMin, safeRMax);
0383     } else {
0384       vecCore__MaskedAssignFunc(safety, !done, (rad - sphere.fRmax));
0385     }
0386     // Distance to r shells over
0387 
0388     // Distance to phi extent
0389     if (!sphere.fFullPhiSphere) {
0390       Real_v safetyPhi = sphere.fPhiWedge.SafetyToIn<Real_v>(point);
0391       vecCore__MaskedAssignFunc(safety, !done, Max(safetyPhi, safety));
0392     }
0393 
0394     // Distance to Theta extent
0395     if (!sphere.fFullThetaSphere) {
0396       Real_v safetyTheta = sphere.fThetaCone.SafetyToIn<Real_v>(point);
0397       vecCore__MaskedAssignFunc(safety, !done, Max(safetyTheta, safety));
0398     }
0399   }
0400 
0401   template <typename Real_v>
0402   VECGEOM_FORCE_INLINE VECCORE_ATT_HOST_DEVICE static void SafetyToOut(UnplacedStruct_t const &sphere,
0403                                                                        Vector3D<Real_v> const &point, Real_v &safety)
0404   {
0405 
0406     using Bool_v = vecCore::Mask_v<Real_v>;
0407     Real_v rad   = point.Mag();
0408 
0409     Bool_v done(false);
0410 
0411     Bool_v completelyinside(false), completelyoutside(false);
0412     GenericKernelForContainsAndInside<Real_v, Bool_v, true>(sphere, point, completelyinside, completelyoutside);
0413     vecCore__MaskedAssignFunc(safety, completelyoutside, Real_v(-1.0));
0414     done |= completelyoutside;
0415     if (vecCore::MaskFull(done)) return;
0416 
0417     Bool_v isOnSurface = !completelyinside && !completelyoutside;
0418     vecCore__MaskedAssignFunc(safety, !done && isOnSurface, Real_v(0.0));
0419     done |= isOnSurface;
0420     if (vecCore::MaskFull(done)) return;
0421 
0422     // Distance to r shells
0423     if (sphere.fRmin) {
0424       Real_v safeRMin = (rad - sphere.fRmin);
0425       Real_v safeRMax = (sphere.fRmax - rad);
0426       safety          = vecCore::Blend(!done && (safeRMin < safeRMax), safeRMin, safeRMax);
0427     } else {
0428       vecCore__MaskedAssignFunc(safety, !done, (sphere.fRmax - rad));
0429     }
0430 
0431     // Distance to phi extent
0432     if (!sphere.fFullPhiSphere) {
0433       Real_v safetyPhi = sphere.fPhiWedge.SafetyToOut<Real_v>(point);
0434       vecCore__MaskedAssignFunc(safety, !done, Min(safetyPhi, safety));
0435     }
0436 
0437     // Distance to Theta extent
0438     Real_v safeTheta(0.);
0439     if (!sphere.fFullThetaSphere) {
0440       safeTheta = sphere.fThetaCone.SafetyToOut<Real_v>(point);
0441       vecCore__MaskedAssignFunc(safety, !done, Min(safeTheta, safety));
0442     }
0443   }
0444 
0445   /* This function should be called from NormalKernel, only for the
0446    * cases when the point is not on the surface and one want to calculate
0447    * the SurfaceNormal.
0448    *
0449    * Algo : Find the boundary which is closest to the point,
0450    * and return the normal to that boundary.
0451    *
0452    */
0453   template <typename Real_v>
0454   VECCORE_ATT_HOST_DEVICE static Vector3D<Real_v> ApproxSurfaceNormalKernel(UnplacedStruct_t const &sphere,
0455                                                                             Vector3D<Real_v> const &point)
0456   {
0457     using vecCore::math::Min;
0458     Vector3D<Real_v> norm(0., 0., 0.);
0459     Real_v radius     = point.Mag();
0460     Real_v distRMax   = Abs(radius - sphere.fRmax);
0461     Real_v distRMin   = InfinityLength<Real_v>();
0462     Real_v distPhi1   = InfinityLength<Real_v>();
0463     Real_v distPhi2   = InfinityLength<Real_v>();
0464     Real_v distTheta1 = InfinityLength<Real_v>();
0465     Real_v distTheta2 = InfinityLength<Real_v>();
0466     Real_v distMin    = distRMax;
0467 
0468     if (sphere.fRmin > 0.) {
0469       distRMin = Abs(sphere.fRmin - radius);
0470       distMin  = Min(distRMin, distRMax);
0471     }
0472 
0473     if (!sphere.fFullPhiSphere) {
0474       distPhi1 = Abs(point.x() * sphere.fPhiWedge.GetNormal1().x() + point.y() * sphere.fPhiWedge.GetNormal1().y());
0475       distPhi2 = Abs(point.x() * sphere.fPhiWedge.GetNormal2().x() + point.y() * sphere.fPhiWedge.GetNormal2().y());
0476       distMin  = Min(distMin, distPhi1, distPhi2);
0477     }
0478 
0479     if (!sphere.fFullThetaSphere) {
0480       Real_v rho = point.Perp();
0481       distTheta1 = sphere.fThetaCone.DistanceToLine<Real_v>(sphere.fThetaCone.GetSlope1(), rho, point.z());
0482       distTheta2 = sphere.fThetaCone.DistanceToLine<Real_v>(sphere.fThetaCone.GetSlope2(), rho, point.z());
0483       distMin    = Min(distMin, distTheta1, distTheta2);
0484     }
0485 
0486     vecCore__MaskedAssignFunc(norm, distMin == distRMax, point.Unit());
0487     vecCore__MaskedAssignFunc(norm, distMin == distRMin, -point.Unit());
0488 
0489     Vector3D<Real_v> normal1 = sphere.fPhiWedge.GetNormal1();
0490     Vector3D<Real_v> normal2 = sphere.fPhiWedge.GetNormal2();
0491     vecCore__MaskedAssignFunc(norm, distMin == distPhi1, -normal1);
0492     vecCore__MaskedAssignFunc(norm, distMin == distPhi2, -normal2);
0493 
0494     vecCore__MaskedAssignFunc(norm, distMin == distTheta1, norm + sphere.fThetaCone.GetNormal1<Real_v>(point));
0495     vecCore__MaskedAssignFunc(norm, distMin == distTheta2, norm + sphere.fThetaCone.GetNormal2<Real_v>(point));
0496 
0497     return norm;
0498   }
0499 
0500   template <typename Real_v>
0501   VECGEOM_FORCE_INLINE VECCORE_ATT_HOST_DEVICE static Vector3D<Real_v> Normal(UnplacedStruct_t const &sphere,
0502                                                                               Vector3D<Real_v> const &point,
0503                                                                               typename vecCore::Mask_v<Real_v> &valid)
0504   {
0505     Vector3D<Real_v> normal(0., 0., 0.);
0506     normal.Set(1e-30);
0507 
0508     using Bool_v = vecCore::Mask_v<Real_v>;
0509 
0510     /* Assumption : This function assumes that the point is on the surface.
0511      *
0512      * Algorithm :
0513      * Detect all those surfaces on which the point is at, and count the
0514      * numOfSurfaces. if(numOfSurfaces == 1) then normal corresponds to the
0515      * normal for that particular case.
0516      *
0517      * if(numOfSurfaces > 1 ), then add the normals corresponds to different
0518      * cases, and finally normalize it and return.
0519      *
0520      * We need following function
0521      * IsPointOnInnerRadius()
0522      * IsPointOnOuterRadius()
0523      * IsPointOnStartPhi()
0524      * IsPointOnEndPhi()
0525      * IsPointOnStartTheta()
0526      * IsPointOnEndTheta()
0527      *
0528      * set valid=true if numOfSurface > 0
0529      *
0530      * if above mentioned assumption not followed , ie.
0531      * In case the given point is outside, then find the closest boundary,
0532      * the required normal will be the normal to that boundary.
0533      * This logic is implemented in "ApproxSurfaceNormalKernel" function
0534      */
0535 
0536     Bool_v isPointOutside(false);
0537 
0538     // May be required Later
0539     /*
0540     if (!ForDistanceToOut) {
0541       Bool_v unused(false);
0542       GenericKernelForContainsAndInside<Real_v, true>(sphere, point, unused, isPointOutside);
0543       vecCore__MaskedAssignFunc(unused || isPointOutside, ApproxSurfaceNormalKernel<Real_v>(sphere, point), &normal);
0544     }
0545     */
0546 
0547     Bool_v isPointInside(false);
0548     GenericKernelForContainsAndInside<Real_v, Bool_v, true>(sphere, point, isPointInside, isPointOutside);
0549     vecCore__MaskedAssignFunc(normal, isPointInside || isPointOutside,
0550                               ApproxSurfaceNormalKernel<Real_v>(sphere, point));
0551 
0552     valid = Bool_v(false);
0553 
0554     Real_v noSurfaces(0.);
0555     Bool_v isPointOnOuterRadius = SphereUtilities::IsPointOnOuterRadius<Real_v>(sphere, point);
0556 
0557     vecCore__MaskedAssignFunc(noSurfaces, isPointOnOuterRadius, noSurfaces + 1);
0558     vecCore__MaskedAssignFunc(normal, !isPointOutside && isPointOnOuterRadius, normal + (point.Unit()));
0559 
0560     if (sphere.fRmin) {
0561       Bool_v isPointOnInnerRadius = SphereUtilities::IsPointOnInnerRadius<Real_v>(sphere, point);
0562       vecCore__MaskedAssignFunc(noSurfaces, isPointOnInnerRadius, noSurfaces + 1);
0563       vecCore__MaskedAssignFunc(normal, !isPointOutside && isPointOnInnerRadius, normal - point.Unit());
0564     }
0565 
0566     if (!sphere.fFullPhiSphere) {
0567       Bool_v isPointOnStartPhi = SphereUtilities::IsPointOnStartPhi<Real_v>(sphere, point);
0568       Bool_v isPointOnEndPhi   = SphereUtilities::IsPointOnEndPhi<Real_v>(sphere, point);
0569       vecCore__MaskedAssignFunc(noSurfaces, isPointOnStartPhi, noSurfaces + 1);
0570       vecCore__MaskedAssignFunc(noSurfaces, isPointOnEndPhi, noSurfaces + 1);
0571       vecCore__MaskedAssignFunc(normal, !isPointOutside && isPointOnStartPhi, normal - sphere.fPhiWedge.GetNormal1());
0572       vecCore__MaskedAssignFunc(normal, !isPointOutside && isPointOnEndPhi, normal - sphere.fPhiWedge.GetNormal2());
0573     }
0574 
0575     if (!sphere.fFullThetaSphere) {
0576       Bool_v isPointOnStartTheta = SphereUtilities::IsPointOnStartTheta<Real_v>(sphere, point);
0577       Bool_v isPointOnEndTheta   = SphereUtilities::IsPointOnEndTheta<Real_v>(sphere, point);
0578 
0579       vecCore__MaskedAssignFunc(noSurfaces, isPointOnStartTheta, noSurfaces + 1);
0580       vecCore__MaskedAssignFunc(normal, !isPointOutside && isPointOnStartTheta,
0581                                 normal + sphere.fThetaCone.GetNormal1<Real_v>(point));
0582 
0583       vecCore__MaskedAssignFunc(noSurfaces, isPointOnEndTheta, noSurfaces + 1);
0584       vecCore__MaskedAssignFunc(normal, !isPointOutside && isPointOnEndTheta,
0585                                 normal + sphere.fThetaCone.GetNormal2<Real_v>(point));
0586 
0587       Vector3D<Real_v> tempNormal(0., 0., -1.);
0588       vecCore__MaskedAssignFunc(
0589           normal, !isPointOutside && isPointOnStartTheta && isPointOnEndTheta && (sphere.eTheta <= kPi / 2.),
0590           tempNormal);
0591       Vector3D<Real_v> tempNormal2(0., 0., 1.);
0592       vecCore__MaskedAssignFunc(
0593           normal, !isPointOutside && isPointOnStartTheta && isPointOnEndTheta && (sphere.fSTheta >= kPi / 2.),
0594           tempNormal2);
0595     }
0596 
0597     normal.Normalize();
0598 
0599     valid = (noSurfaces > Real_v(0.));
0600 
0601     return normal;
0602   }
0603 };
0604 } // namespace VECGEOM_IMPL_NAMESPACE
0605 } // namespace vecgeom
0606 
0607 #endif // VECGEOM_VOLUMES_KERNEL_sphereIMPLEMENTATION_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!