Back to home page

EIC code displayed by LXR
 
 

    

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

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(0.)), sd1);
0184     } else {
0185       tmpPt = point + sd1 * direction;
0186       vecCore::MaskedAssign(outerDist,
0187                             !done && sphere.fPhiWedge.Contains<Real_v>(tmpPt) &&
0188                                 sphere.fThetaCone.Contains<Real_v>(tmpPt) && (sd1 >= Real_v(0.)),
0189                             sd1);
0190     }
0191 
0192     if (sphere.fRmin) {
0193       c  = rad2 - sphere.fRmin * sphere.fRmin;
0194       d2 = pDotV3d * pDotV3d - c;
0195       // Note: Abs(d2) was introduced to avoid Sqrt(negative) in other lanes than the ones satisfying d2>=0.
0196       vecCore__MaskedAssignFunc(sd2, d2 >= Real_v(0.), (-pDotV3d + Sqrt(Abs(d2))));
0197 
0198       if (sphere.fFullSphere) {
0199         vecCore::MaskedAssign(innerDist, !done && (sd2 >= Real_v(0.)), sd2);
0200       } else {
0201         //   std::cout<<" ---- Called by InnerRad ---- " << std::endl;
0202         tmpPt = point + sd2 * direction;
0203         vecCore::MaskedAssign(innerDist,
0204                               !done && (sd2 >= Real_v(0.)) && sphere.fPhiWedge.Contains<Real_v>(tmpPt) &&
0205                                   sphere.fThetaCone.Contains<Real_v>(tmpPt),
0206                               sd2);
0207       }
0208     }
0209 
0210     distance = Min(outerDist, innerDist);
0211 
0212     if (!fullPhiSphere) {
0213       GetMinDistFromPhi<Real_v, true>(sphere, point, direction, done, distance);
0214     }
0215 
0216     Real_v distThetaMin(kInfLength);
0217 
0218     if (!fullThetaSphere) {
0219       Bool_v intsect1(false);
0220       Bool_v intsect2(false);
0221       Real_v distTheta1(kInfLength);
0222       Real_v distTheta2(kInfLength);
0223 
0224       sphere.fThetaCone.DistanceToIn<Real_v>(point, direction, distTheta1, distTheta2, intsect1,
0225                                              intsect2); //,cone1IntSecPt, cone2IntSecPt);
0226       Vector3D<Real_v> coneIntSecPt1 = point + distTheta1 * direction;
0227       Real_v distCone1               = coneIntSecPt1.Mag2();
0228 
0229       Vector3D<Real_v> coneIntSecPt2 = point + distTheta2 * direction;
0230       Real_v distCone2               = coneIntSecPt2.Mag2();
0231 
0232       Bool_v isValidCone1 =
0233           (distCone1 >= sphere.fRmin * sphere.fRmin && distCone1 <= sphere.fRmax * sphere.fRmax) && intsect1;
0234       Bool_v isValidCone2 =
0235           (distCone2 >= sphere.fRmin * sphere.fRmin && distCone2 <= sphere.fRmax * sphere.fRmax) && intsect2;
0236 
0237       if (!fullPhiSphere) {
0238         isValidCone1 &= sphere.fPhiWedge.Contains<Real_v>(coneIntSecPt1);
0239         isValidCone2 &= sphere.fPhiWedge.Contains<Real_v>(coneIntSecPt2);
0240       }
0241       vecCore::MaskedAssign(distThetaMin, (!done && isValidCone2 && !isValidCone1), distTheta2);
0242       vecCore::MaskedAssign(distThetaMin, (!done && isValidCone1 && !isValidCone2), distTheta1);
0243       vecCore__MaskedAssignFunc(distThetaMin, (!done && isValidCone1 && isValidCone2), Min(distTheta1, distTheta2));
0244     }
0245 
0246     distance = Min(distThetaMin, distance);
0247 
0248     Vector3D<Real_v> directDir = (Vector3D<Real_v>(0., 0., 0.) - point);
0249     Real_v newDist             = directDir.Mag();
0250     vecCore__MaskedAssignFunc(distance,
0251                               Bool_v(sphere.fSTheta > kHalfTolerance || sphere.eTheta < (kPi - kHalfTolerance)) &&
0252                                   (Abs(directDir.Unit().x() - direction.x()) < kHalfTolerance) &&
0253                                   (Abs(directDir.Unit().y() - direction.y()) < kHalfTolerance) &&
0254                                   (Abs(directDir.Unit().z() - direction.z()) < kHalfTolerance),
0255                               Min(distance, newDist));
0256   }
0257 
0258   template <typename Real_v>
0259   VECGEOM_FORCE_INLINE
0260   VECCORE_ATT_HOST_DEVICE
0261   static void DistanceToOut(UnplacedStruct_t const &sphere, Vector3D<Real_v> const &point,
0262                             Vector3D<Real_v> const &direction, Real_v const & /* stepMax */, Real_v &distance)
0263   {
0264 
0265     using Bool_v = typename vecCore::Mask_v<Real_v>;
0266 
0267     distance = kInfLength;
0268     Bool_v done(false);
0269 
0270     Real_v snxt(kInfLength);
0271 
0272     // Intersection point
0273     Vector3D<Real_v> intSecPt;
0274     Real_v d2(0.);
0275 
0276     Real_v pDotV3d = point.Dot(direction);
0277 
0278     Real_v rad2 = point.Mag2();
0279     Real_v c    = rad2 - sphere.fRmax * sphere.fRmax;
0280 
0281     Real_v sd1(kInfLength);
0282     Real_v sd2(kInfLength);
0283 
0284     Bool_v cond = SphereUtilities::IsCompletelyOutside<Real_v>(sphere, point);
0285     vecCore__MaskedAssignFunc(distance, cond, Real_v(-1.0));
0286 
0287     done |= cond;
0288     if (vecCore::MaskFull(done)) return;
0289 
0290     cond = SphereUtilities::IsPointOnSurfaceAndMovingOut<Real_v, true>(sphere, point, direction);
0291     vecCore__MaskedAssignFunc(distance, !done && cond, Real_v(0.0));
0292     done |= cond;
0293     if (vecCore::MaskFull(done)) return;
0294 
0295     // Note: Abs(d2) was introduced to avoid Sqrt(negative) in other lanes than the ones satisfying d2>=0.
0296     d2 = (pDotV3d * pDotV3d - c);
0297     vecCore__MaskedAssignFunc(sd1, (!done && (d2 >= Real_v(0.))), (-pDotV3d + Sqrt(Abs(d2))));
0298 
0299     if (sphere.fRmin) {
0300       c  = rad2 - sphere.fRmin * sphere.fRmin;
0301       d2 = (pDotV3d * pDotV3d - c);
0302       vecCore__MaskedAssignFunc(sd2, (!done && (d2 >= Real_v(0.)) && (pDotV3d < Real_v(0.))),
0303                                 (-pDotV3d - Sqrt(Abs(d2))));
0304     }
0305 
0306     snxt = Min(sd1, sd2);
0307 
0308     Bool_v condSemi = (Bool_v(sphere.fSTheta == 0. && sphere.eTheta == kPi / 2.) && direction.z() >= Real_v(0.)) ||
0309                       (Bool_v(sphere.fSTheta == kPi / 2. && sphere.eTheta == kPi) && direction.z() <= Real_v(0.));
0310     vecCore::MaskedAssign(distance, !done && condSemi, snxt);
0311     done |= condSemi;
0312     if (vecCore::MaskFull(done)) return;
0313 
0314     Real_v distThetaMin(kInfLength);
0315     Real_v distPhiMin(kInfLength);
0316 
0317     if (!sphere.fFullThetaSphere) {
0318       Bool_v intsect1(false);
0319       Bool_v intsect2(false);
0320       Real_v distTheta1(kInfLength);
0321       Real_v distTheta2(kInfLength);
0322       sphere.fThetaCone.DistanceToOut<Real_v>(point, direction, distTheta1, distTheta2, intsect1, intsect2);
0323       vecCore::MaskedAssign(distThetaMin, (intsect2 && !intsect1), distTheta2);
0324       vecCore::MaskedAssign(distThetaMin, (!intsect2 && intsect1), distTheta1);
0325       vecCore__MaskedAssignFunc(distThetaMin, (intsect2 && intsect1), Min(distTheta1, distTheta2));
0326     }
0327 
0328     distance = Min(distThetaMin, snxt);
0329 
0330     if (!sphere.fFullPhiSphere) {
0331       if (sphere.fDPhi <= kPi) {
0332         Real_v distPhi1;
0333         Real_v distPhi2;
0334         sphere.fPhiWedge.DistanceToOut<Real_v>(point, direction, distPhi1, distPhi2);
0335         distPhiMin = Min(distPhi1, distPhi2);
0336         distance   = Min(distPhiMin, distance);
0337       } else {
0338         GetMinDistFromPhi<Real_v, false>(sphere, point, direction, done, distance);
0339       }
0340     }
0341   }
0342 
0343   template <typename Real_v, bool DistToIn>
0344   VECCORE_ATT_HOST_DEVICE
0345   static void GetMinDistFromPhi(UnplacedStruct_t const &sphere, Vector3D<Real_v> const &localPoint,
0346                                 Vector3D<Real_v> const &localDir, typename vecCore::Mask_v<Real_v> &done,
0347                                 Real_v &distance)
0348   {
0349     using Bool_v = typename vecCore::Mask_v<Real_v>;
0350     Real_v distPhi1(kInfLength);
0351     Real_v distPhi2(kInfLength);
0352     Real_v dist(kInfLength);
0353 
0354     if (DistToIn)
0355       sphere.fPhiWedge.DistanceToIn<Real_v>(localPoint, localDir, distPhi1, distPhi2);
0356     else
0357       sphere.fPhiWedge.DistanceToOut<Real_v>(localPoint, localDir, distPhi1, distPhi2);
0358 
0359     Bool_v containsCond1(false), containsCond2(false);
0360     // Min Face
0361     dist                   = Min(distPhi1, distPhi2);
0362     Vector3D<Real_v> tmpPt = localPoint + dist * localDir;
0363     Real_v rad2            = tmpPt.Mag2();
0364 
0365     Bool_v tempCond(false);
0366     tempCond = ((dist == distPhi1) && sphere.fPhiWedge.IsOnSurfaceGeneric<Real_v, true>(tmpPt)) ||
0367                ((dist == distPhi2) && sphere.fPhiWedge.IsOnSurfaceGeneric<Real_v, false>(tmpPt));
0368 
0369     containsCond1 = tempCond && (rad2 > sphere.fRmin * sphere.fRmin) && (rad2 < sphere.fRmax * sphere.fRmax) &&
0370                     sphere.fThetaCone.Contains<Real_v>(tmpPt);
0371 
0372     vecCore__MaskedAssignFunc(distance, !done && containsCond1, Min(dist, distance));
0373 
0374     // Max Face
0375     dist  = Max(distPhi1, distPhi2);
0376     tmpPt = localPoint + dist * localDir;
0377 
0378     rad2     = tmpPt.Mag2();
0379     tempCond = Bool_v(false);
0380     tempCond = ((dist == distPhi1) && sphere.fPhiWedge.IsOnSurfaceGeneric<Real_v, true>(tmpPt)) ||
0381                ((dist == distPhi2) && sphere.fPhiWedge.IsOnSurfaceGeneric<Real_v, false>(tmpPt));
0382 
0383     containsCond2 = tempCond && (rad2 > sphere.fRmin * sphere.fRmin) && (rad2 < sphere.fRmax * sphere.fRmax) &&
0384                     sphere.fThetaCone.Contains<Real_v>(tmpPt);
0385     vecCore__MaskedAssignFunc(distance, ((!done) && (!containsCond1) && containsCond2), Min(dist, distance));
0386   }
0387 
0388   template <typename Real_v>
0389   VECGEOM_FORCE_INLINE
0390   VECCORE_ATT_HOST_DEVICE
0391   static void SafetyToIn(UnplacedStruct_t const &sphere, Vector3D<Real_v> const &point, Real_v &safety)
0392   {
0393     using Bool_v = vecCore::Mask_v<Real_v>;
0394     Bool_v done(false);
0395 
0396     // General Precalcs
0397     Real_v rad = point.Mag();
0398 
0399     Real_v safeRMin(0.);
0400     Real_v safeRMax(0.);
0401 
0402     Bool_v completelyinside(false), completelyoutside(false);
0403     GenericKernelForContainsAndInside<Real_v, Bool_v, true>(sphere, point, completelyinside, completelyoutside);
0404 
0405     vecCore__MaskedAssignFunc(safety, completelyinside, Real_v(-1.0));
0406     done |= completelyinside;
0407     if (vecCore::MaskFull(done)) return;
0408 
0409     Bool_v isOnSurface = !completelyinside && !completelyoutside;
0410     vecCore__MaskedAssignFunc(safety, !done && isOnSurface, Real_v(0.0));
0411     done |= isOnSurface;
0412     if (vecCore::MaskFull(done)) return;
0413 
0414     if (sphere.fRmin) {
0415       safeRMin = sphere.fRmin - rad;
0416       safeRMax = rad - sphere.fRmax;
0417       safety   = vecCore::Blend(!done && (safeRMin > safeRMax), safeRMin, safeRMax);
0418     } else {
0419       vecCore__MaskedAssignFunc(safety, !done, (rad - sphere.fRmax));
0420     }
0421     // Distance to r shells over
0422 
0423     // Distance to phi extent
0424     if (!sphere.fFullPhiSphere) {
0425       Real_v safetyPhi = sphere.fPhiWedge.SafetyToIn<Real_v>(point);
0426       vecCore__MaskedAssignFunc(safety, !done, Max(safetyPhi, safety));
0427     }
0428 
0429     // Distance to Theta extent
0430     if (!sphere.fFullThetaSphere) {
0431       Real_v safetyTheta = sphere.fThetaCone.SafetyToIn<Real_v>(point);
0432       vecCore__MaskedAssignFunc(safety, !done, Max(safetyTheta, safety));
0433     }
0434   }
0435 
0436   template <typename Real_v>
0437   VECGEOM_FORCE_INLINE
0438   VECCORE_ATT_HOST_DEVICE
0439   static void SafetyToOut(UnplacedStruct_t const &sphere, Vector3D<Real_v> const &point, Real_v &safety)
0440   {
0441 
0442     using Bool_v = vecCore::Mask_v<Real_v>;
0443     Real_v rad   = point.Mag();
0444 
0445     Bool_v done(false);
0446 
0447     Bool_v completelyinside(false), completelyoutside(false);
0448     GenericKernelForContainsAndInside<Real_v, Bool_v, true>(sphere, point, completelyinside, completelyoutside);
0449     vecCore__MaskedAssignFunc(safety, completelyoutside, Real_v(-1.0));
0450     done |= completelyoutside;
0451     if (vecCore::MaskFull(done)) return;
0452 
0453     Bool_v isOnSurface = !completelyinside && !completelyoutside;
0454     vecCore__MaskedAssignFunc(safety, !done && isOnSurface, Real_v(0.0));
0455     done |= isOnSurface;
0456     if (vecCore::MaskFull(done)) return;
0457 
0458     // Distance to r shells
0459     if (sphere.fRmin) {
0460       Real_v safeRMin = (rad - sphere.fRmin);
0461       Real_v safeRMax = (sphere.fRmax - rad);
0462       safety          = vecCore::Blend(!done && (safeRMin < safeRMax), safeRMin, safeRMax);
0463     } else {
0464       vecCore__MaskedAssignFunc(safety, !done, (sphere.fRmax - rad));
0465     }
0466 
0467     // Distance to phi extent
0468     if (!sphere.fFullPhiSphere) {
0469       Real_v safetyPhi = sphere.fPhiWedge.SafetyToOut<Real_v>(point);
0470       vecCore__MaskedAssignFunc(safety, !done, Min(safetyPhi, safety));
0471     }
0472 
0473     // Distance to Theta extent
0474     Real_v safeTheta(0.);
0475     if (!sphere.fFullThetaSphere) {
0476       safeTheta = sphere.fThetaCone.SafetyToOut<Real_v>(point);
0477       vecCore__MaskedAssignFunc(safety, !done, Min(safeTheta, safety));
0478     }
0479   }
0480 
0481   /* This function should be called from NormalKernel, only for the
0482    * cases when the point is not on the surface and one want to calculate
0483    * the SurfaceNormal.
0484    *
0485    * Algo : Find the boundary which is closest to the point,
0486    * and return the normal to that boundary.
0487    *
0488    */
0489   template <typename Real_v>
0490   VECCORE_ATT_HOST_DEVICE
0491   static Vector3D<Real_v> ApproxSurfaceNormalKernel(UnplacedStruct_t const &sphere, Vector3D<Real_v> const &point)
0492   {
0493     using vecCore::math::Min;
0494     Vector3D<Real_v> norm(0., 0., 0.);
0495     Real_v radius     = point.Mag();
0496     Real_v distRMax   = Abs(radius - sphere.fRmax);
0497     Real_v distRMin   = InfinityLength<Real_v>();
0498     Real_v distPhi1   = InfinityLength<Real_v>();
0499     Real_v distPhi2   = InfinityLength<Real_v>();
0500     Real_v distTheta1 = InfinityLength<Real_v>();
0501     Real_v distTheta2 = InfinityLength<Real_v>();
0502     Real_v distMin    = distRMax;
0503 
0504     if (sphere.fRmin > 0.) {
0505       distRMin = Abs(sphere.fRmin - radius);
0506       distMin  = Min(distRMin, distRMax);
0507     }
0508 
0509     if (!sphere.fFullPhiSphere) {
0510       distPhi1 = Abs(point.x() * sphere.fPhiWedge.GetNormal1().x() + point.y() * sphere.fPhiWedge.GetNormal1().y());
0511       distPhi2 = Abs(point.x() * sphere.fPhiWedge.GetNormal2().x() + point.y() * sphere.fPhiWedge.GetNormal2().y());
0512       distMin  = Min(distMin, distPhi1, distPhi2);
0513     }
0514 
0515     if (!sphere.fFullThetaSphere) {
0516       Real_v rho = point.Perp();
0517       distTheta1 = sphere.fThetaCone.DistanceToLine<Real_v>(sphere.fThetaCone.GetSlope1(), rho, point.z());
0518       distTheta2 = sphere.fThetaCone.DistanceToLine<Real_v>(sphere.fThetaCone.GetSlope2(), rho, point.z());
0519       distMin    = Min(distMin, distTheta1, distTheta2);
0520     }
0521 
0522     vecCore__MaskedAssignFunc(norm, distMin == distRMax, point.Unit());
0523     vecCore__MaskedAssignFunc(norm, distMin == distRMin, -point.Unit());
0524 
0525     Vector3D<Real_v> normal1 = sphere.fPhiWedge.GetNormal1();
0526     Vector3D<Real_v> normal2 = sphere.fPhiWedge.GetNormal2();
0527     vecCore__MaskedAssignFunc(norm, distMin == distPhi1, -normal1);
0528     vecCore__MaskedAssignFunc(norm, distMin == distPhi2, -normal2);
0529 
0530     vecCore__MaskedAssignFunc(norm, distMin == distTheta1, norm + sphere.fThetaCone.GetNormal1<Real_v>(point));
0531     vecCore__MaskedAssignFunc(norm, distMin == distTheta2, norm + sphere.fThetaCone.GetNormal2<Real_v>(point));
0532 
0533     return norm;
0534   }
0535 
0536   template <typename Real_v>
0537   VECGEOM_FORCE_INLINE
0538   VECCORE_ATT_HOST_DEVICE
0539   static Vector3D<Real_v> Normal(UnplacedStruct_t const &sphere, Vector3D<Real_v> const &point,
0540                                  typename vecCore::Mask_v<Real_v> &valid)
0541   {
0542     Vector3D<Real_v> normal(0., 0., 0.);
0543     normal.Set(1e-30);
0544 
0545     using Bool_v = vecCore::Mask_v<Real_v>;
0546 
0547     /* Assumption : This function assumes that the point is on the surface.
0548      *
0549      * Algorithm :
0550      * Detect all those surfaces on which the point is at, and count the
0551      * numOfSurfaces. if(numOfSurfaces == 1) then normal corresponds to the
0552      * normal for that particular case.
0553      *
0554      * if(numOfSurfaces > 1 ), then add the normals corresponds to different
0555      * cases, and finally normalize it and return.
0556      *
0557      * We need following function
0558      * IsPointOnInnerRadius()
0559      * IsPointOnOuterRadius()
0560      * IsPointOnStartPhi()
0561      * IsPointOnEndPhi()
0562      * IsPointOnStartTheta()
0563      * IsPointOnEndTheta()
0564      *
0565      * set valid=true if numOfSurface > 0
0566      *
0567      * if above mentioned assumption not followed , ie.
0568      * In case the given point is outside, then find the closest boundary,
0569      * the required normal will be the normal to that boundary.
0570      * This logic is implemented in "ApproxSurfaceNormalKernel" function
0571      */
0572 
0573     Bool_v isPointOutside(false);
0574 
0575     // May be required Later
0576     /*
0577     if (!ForDistanceToOut) {
0578       Bool_v unused(false);
0579       GenericKernelForContainsAndInside<Real_v, true>(sphere, point, unused, isPointOutside);
0580       vecCore__MaskedAssignFunc(unused || isPointOutside, ApproxSurfaceNormalKernel<Real_v>(sphere, point), &normal);
0581     }
0582     */
0583 
0584     Bool_v isPointInside(false);
0585     GenericKernelForContainsAndInside<Real_v, Bool_v, true>(sphere, point, isPointInside, isPointOutside);
0586     vecCore__MaskedAssignFunc(normal, isPointInside || isPointOutside,
0587                               ApproxSurfaceNormalKernel<Real_v>(sphere, point));
0588 
0589     valid = Bool_v(false);
0590 
0591     Real_v noSurfaces(0.);
0592     Bool_v isPointOnOuterRadius = SphereUtilities::IsPointOnOuterRadius<Real_v>(sphere, point);
0593 
0594     vecCore__MaskedAssignFunc(noSurfaces, isPointOnOuterRadius, noSurfaces + 1);
0595     vecCore__MaskedAssignFunc(normal, !isPointOutside && isPointOnOuterRadius, normal + (point.Unit()));
0596 
0597     if (sphere.fRmin) {
0598       Bool_v isPointOnInnerRadius = SphereUtilities::IsPointOnInnerRadius<Real_v>(sphere, point);
0599       vecCore__MaskedAssignFunc(noSurfaces, isPointOnInnerRadius, noSurfaces + 1);
0600       vecCore__MaskedAssignFunc(normal, !isPointOutside && isPointOnInnerRadius, normal - point.Unit());
0601     }
0602 
0603     if (!sphere.fFullPhiSphere) {
0604       Bool_v isPointOnStartPhi = SphereUtilities::IsPointOnStartPhi<Real_v>(sphere, point);
0605       Bool_v isPointOnEndPhi   = SphereUtilities::IsPointOnEndPhi<Real_v>(sphere, point);
0606       vecCore__MaskedAssignFunc(noSurfaces, isPointOnStartPhi, noSurfaces + 1);
0607       vecCore__MaskedAssignFunc(noSurfaces, isPointOnEndPhi, noSurfaces + 1);
0608       vecCore__MaskedAssignFunc(normal, !isPointOutside && isPointOnStartPhi, normal - sphere.fPhiWedge.GetNormal1());
0609       vecCore__MaskedAssignFunc(normal, !isPointOutside && isPointOnEndPhi, normal - sphere.fPhiWedge.GetNormal2());
0610     }
0611 
0612     if (!sphere.fFullThetaSphere) {
0613       Bool_v isPointOnStartTheta = SphereUtilities::IsPointOnStartTheta<Real_v>(sphere, point);
0614       Bool_v isPointOnEndTheta   = SphereUtilities::IsPointOnEndTheta<Real_v>(sphere, point);
0615 
0616       vecCore__MaskedAssignFunc(noSurfaces, isPointOnStartTheta, noSurfaces + 1);
0617       vecCore__MaskedAssignFunc(normal, !isPointOutside && isPointOnStartTheta,
0618                                 normal + sphere.fThetaCone.GetNormal1<Real_v>(point));
0619 
0620       vecCore__MaskedAssignFunc(noSurfaces, isPointOnEndTheta, noSurfaces + 1);
0621       vecCore__MaskedAssignFunc(normal, !isPointOutside && isPointOnEndTheta,
0622                                 normal + sphere.fThetaCone.GetNormal2<Real_v>(point));
0623 
0624       Vector3D<Real_v> tempNormal(0., 0., -1.);
0625       vecCore__MaskedAssignFunc(
0626           normal, !isPointOutside && isPointOnStartTheta && isPointOnEndTheta && (sphere.eTheta <= kPi / 2.),
0627           tempNormal);
0628       Vector3D<Real_v> tempNormal2(0., 0., 1.);
0629       vecCore__MaskedAssignFunc(
0630           normal, !isPointOutside && isPointOnStartTheta && isPointOnEndTheta && (sphere.fSTheta >= kPi / 2.),
0631           tempNormal2);
0632     }
0633 
0634     normal.Normalize();
0635 
0636     valid = (noSurfaces > Real_v(0.));
0637 
0638     return normal;
0639   }
0640 };
0641 } // namespace VECGEOM_IMPL_NAMESPACE
0642 } // namespace vecgeom
0643 
0644 #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!