Back to home page

EIC code displayed by LXR
 
 

    

File indexing completed on 2025-09-17 09:18:24 

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