Back to home page

EIC code displayed by LXR
 
 

    

File indexing completed on 2026-09-27 09:30:13 

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 /// This file implements the algorithms for Paraboloid shape
0006 /// @file volumes/kernel/ParaboloidImplementation.h
0007 /// @author Marilena Bandieramonte
0008 ///
0009 /// A paraboloid is the solid bounded by the following surfaces:
0010 /// - 2 planes parallel with XY cutting the Z axis at z = -dz and z = +dz
0011 /// - the surface of revolution of a parabola described by: z = a * (x^2 + y^2) + b
0012 ///
0013 /// The parameters a and b are automatically computed from:
0014 /// - rlo - radius of the circle of intersection between the
0015 /// parabolic surface and the plane z = -dz
0016 /// - rhi - the radius of the circle of intersection between the
0017 /// parabolic surface and the plane z = +dz
0018 /// - dz = a * rhi^2 + b and  -dz = a * rlo^2 + b, where rhi > rlo, both >= 0
0019 /// - a = 2 * dz * dd and b = -dz * (rlo^2 + rhi^2) * dd, where dd = 1 / (rhi^2 - rlo^2)
0020 
0021 #ifndef VECGEOM_VOLUMES_KERNEL_PARABOLOIDIMPLEMENTATION_H_
0022 #define VECGEOM_VOLUMES_KERNEL_PARABOLOIDIMPLEMENTATION_H_
0023 
0024 #include "VecGeom/base/Vector3D.h"
0025 #include "VecGeom/volumes/ParaboloidStruct.h"
0026 #include "VecGeom/volumes/kernel/GenericKernels.h"
0027 #include <VecCore/VecCore>
0028 
0029 #include <cstdio>
0030 
0031 namespace vecgeom {
0032 
0033 VECGEOM_DEVICE_FORWARD_DECLARE(struct ParaboloidImplementation;);
0034 VECGEOM_DEVICE_DECLARE_CONV(struct, ParaboloidImplementation);
0035 
0036 inline namespace VECGEOM_IMPL_NAMESPACE {
0037 
0038 class PlacedParaboloid;
0039 template <typename T>
0040 struct ParaboloidStruct;
0041 class UnplacedParaboloid;
0042 
0043 struct ParaboloidImplementation {
0044 
0045   using PlacedShape_t    = PlacedParaboloid;
0046   using UnplacedStruct_t = ParaboloidStruct<Precision>;
0047   using UnplacedVolume_t = UnplacedParaboloid;
0048 
0049   template <typename Real_v, typename Bool_v>
0050   VECGEOM_FORCE_INLINE VECCORE_ATT_HOST_DEVICE static void Contains(UnplacedStruct_t const &paraboloid,
0051                                                                     Vector3D<Real_v> const &point, Bool_v &inside)
0052   {
0053     Bool_v unused(false), outside(false);
0054     GenericKernelForContainsAndInside<Real_v, Bool_v, false>(paraboloid, point, unused, outside);
0055     inside = !outside;
0056   }
0057 
0058   // BIG QUESTION: DO WE WANT TO GIVE ALL 3 TEMPLATE PARAMETERS
0059   // -- OR -- DO WE WANT TO DEDUCE Bool_v, Index_t from Real_v???
0060   template <typename Real_v, typename Inside_t>
0061   VECGEOM_FORCE_INLINE VECCORE_ATT_HOST_DEVICE static void Inside(UnplacedStruct_t const &paraboloid,
0062                                                                   Vector3D<Real_v> const &point, Inside_t &inside)
0063   {
0064 
0065     using Bool_v       = vecCore::Mask_v<Real_v>;
0066     using InsideBool_v = vecCore::Mask_v<Inside_t>;
0067     Bool_v completelyinside, completelyoutside;
0068     GenericKernelForContainsAndInside<Real_v, Bool_v, true>(paraboloid, point, completelyinside, completelyoutside);
0069     inside = EInside::kSurface;
0070     vecCore::MaskedAssign(inside, (InsideBool_v)completelyoutside, Inside_t(EInside::kOutside));
0071     vecCore::MaskedAssign(inside, (InsideBool_v)completelyinside, Inside_t(EInside::kInside));
0072   }
0073 
0074   template <typename Real_v, typename Bool_v, bool ForInside>
0075   VECGEOM_FORCE_INLINE VECCORE_ATT_HOST_DEVICE static void GenericKernelForContainsAndInside(
0076       UnplacedStruct_t const &paraboloid, Vector3D<Real_v> const &point, Bool_v &completelyinside,
0077       Bool_v &completelyoutside)
0078   {
0079     // using Bool_v = vecCore::Mask_v<Real_v>;
0080     completelyinside  = Bool_v(false);
0081     completelyoutside = Bool_v(false);
0082 
0083     Real_v rho2       = point.Perp2();
0084     Real_v paraRho2   = paraboloid.fK1 * point.z() + paraboloid.fK2;
0085     Real_v diff       = rho2 - paraRho2;
0086     Real_v absZ       = Abs(point.z());
0087     completelyoutside = (absZ > Real_v(paraboloid.fDz + kTolerance)) || (diff > kTolerance);
0088     if (vecCore::MaskFull(completelyoutside)) return;
0089     if (ForInside) completelyinside = (absZ < Real_v(paraboloid.fDz - kTolerance)) && (diff < -kTolerance);
0090   }
0091 
0092   template <typename Real_v, bool ForTopZPlane>
0093   VECCORE_ATT_HOST_DEVICE static vecCore::Mask_v<Real_v> IsOnZPlane(UnplacedStruct_t const &paraboloid,
0094                                                                     Vector3D<Real_v> const &point)
0095   {
0096     Real_v rho2 = point.Perp2();
0097     if (ForTopZPlane) {
0098       return Abs(point.z() - paraboloid.fDz) < kTolerance && rho2 < (paraboloid.fRhi2 + kHalfTolerance);
0099     } else {
0100       return Abs(point.z() + paraboloid.fDz) < kTolerance && rho2 < (paraboloid.fRlo2 + kHalfTolerance);
0101     }
0102   }
0103 
0104   template <typename Real_v>
0105   VECCORE_ATT_HOST_DEVICE static vecCore::Mask_v<Real_v> IsOnParabolicSurface(UnplacedStruct_t const &paraboloid,
0106                                                                               Vector3D<Real_v> const &point)
0107   {
0108 
0109     using Bool_v = vecCore::Mask_v<Real_v>;
0110 
0111     Real_v value              = paraboloid.fA * point.Perp2() + paraboloid.fB - point.z();
0112     Bool_v onParabolicSurface = value > -kTolerance && value < kTolerance;
0113     return onParabolicSurface;
0114   }
0115 
0116   template <typename Real_v>
0117   VECGEOM_FORCE_INLINE VECCORE_ATT_HOST_DEVICE static void DistanceToIn(UnplacedStruct_t const &paraboloid,
0118                                                                         Vector3D<Real_v> const &point,
0119                                                                         Vector3D<Real_v> const &direction,
0120                                                                         Real_v const & /* stepMax */, Real_v &distance)
0121   {
0122 
0123     using Bool_v = vecCore::Mask_v<Real_v>;
0124 
0125     Bool_v done(false);
0126     distance = InfinityLength<Real_v>();
0127     Real_v offset(0.);
0128     Vector3D<Real_v> p(point);
0129 
0130     // Move point closer, if required
0131     Precision Rsph = 1.5 * vecCore::math::Max(paraboloid.fDx, paraboloid.fDz);
0132     Real_v Rfar2(1024. * Rsph * Rsph); // 1024 = 32 * 32
0133     vecCore__MaskedAssignFunc(offset, ((p.Mag2() > Rfar2) && (direction.Dot(p) < Real_v(0.))),
0134                               p.Mag() - Real_v(2.) * Rsph);
0135     p += offset * direction;
0136 
0137     Real_v absZ   = Abs(p.z());
0138     Real_v rho2   = p.Perp2(); // p.x()*p.x()+p.y()*p.y();
0139     Bool_v checkZ = p.z() * direction.z() >= Real_v(0.);
0140 
0141     // check if the point is distancing in Z
0142     Bool_v isDistancingInZ = (absZ > paraboloid.fDz && checkZ);
0143     done |= isDistancingInZ;
0144     if (vecCore::MaskFull(done)) return;
0145 
0146     Real_v paraRho2 = paraboloid.fK1 * p.z() + paraboloid.fK2;
0147     Real_v diff     = rho2 - paraRho2;
0148 
0149     vecCore__MaskedAssignFunc(distance, !done, Real_v(-1.));
0150     Bool_v insideZ                              = absZ < Real_v(paraboloid.fDz - kTolerance);
0151     Bool_v insideParabolicSurfaceOuterTolerance = (diff < -kTolerance);
0152     done |= !done && (insideZ && insideParabolicSurfaceOuterTolerance);
0153     if (vecCore::MaskFull(done)) return;
0154 
0155     Bool_v isOnZPlaneAndMovingInside = (IsOnZPlane<Real_v, true>(paraboloid, point) && direction.z() < Real_v(0.)) ||
0156                                        (IsOnZPlane<Real_v, false>(paraboloid, point) && direction.z() > Real_v(0.));
0157     vecCore__MaskedAssignFunc(distance, !done && isOnZPlaneAndMovingInside, Real_v(0.));
0158     done |= isOnZPlaneAndMovingInside;
0159     if (vecCore::MaskFull(done)) return;
0160 
0161     Vector3D<Real_v> normal(p.x(), p.y(), Real_v(-paraboloid.fK1 * Real_v(0.5)));
0162     Bool_v isOnParabolicSurfaceAndMovingInside =
0163         diff > -kTolerance && diff < kTolerance && direction.Dot(normal) < Real_v(0.);
0164     vecCore__MaskedAssignFunc(distance, !done && isOnParabolicSurfaceAndMovingInside, Real_v(0.));
0165     done |= isOnParabolicSurfaceAndMovingInside;
0166     if (vecCore::MaskFull(done)) return;
0167 
0168     vecCore__MaskedAssignFunc(distance, !done, InfinityLength<Real_v>());
0169 
0170     /* Intersection tests with Z planes are not required if the point is within Z Range
0171      * In this case it will either intsect with parabolic surface or not intersect at all.
0172      */
0173     if (!vecCore::MaskFull(absZ < paraboloid.fDz)) {
0174       Real_v distZ(InfinityLength<Real_v>());                            // = (absZ - paraboloid.fDz) / absDirZ;
0175       Bool_v bottomPlane = p.z() < -paraboloid.fDz && direction.z() > 0; //(true);
0176       Bool_v topPlane    = p.z() > paraboloid.fDz && direction.z() < 0;
0177       vecCore__MaskedAssignFunc(distZ, topPlane, (paraboloid.fDz - p.z()) / NonZero(direction.z()));
0178       vecCore__MaskedAssignFunc(distZ, bottomPlane, (-paraboloid.fDz - p.z()) / NonZero(direction.z()));
0179       Real_v xHit    = p.x() + distZ * direction.x();
0180       Real_v yHit    = p.y() + distZ * direction.y();
0181       Real_v rhoHit2 = xHit * xHit + yHit * yHit;
0182 
0183       vecCore::MaskedAssign(distance, !done && topPlane && rhoHit2 <= paraboloid.fRhi2, distZ + offset);
0184       done |= topPlane && rhoHit2 < paraboloid.fRhi2;
0185       if (vecCore::MaskFull(done)) return;
0186 
0187       vecCore::MaskedAssign(distance, !done && bottomPlane && rhoHit2 <= paraboloid.fRlo2, distZ + offset);
0188       done |= (bottomPlane && rhoHit2 <= paraboloid.fRlo2); // || (topPlane && rhoHit2 <= paraboloid.fRhi2);
0189       if (vecCore::MaskFull(done)) return;
0190     }
0191 
0192     /* Intersection tests with Parabolic surface are not required if the point is above
0193      * top Z plane Radius of point is less the Rhi. In this case depending upon the
0194      * direction it will either intersect with top Z plane or not intersect at all
0195      */
0196     if (!vecCore::MaskFull(p.z() > paraboloid.fDz && rho2 < paraboloid.fRhi2)) {
0197       // Quadratic Solver for Parabolic surface
0198       Real_v dirRho2 = direction.Perp2();
0199       Real_v pDotV2D = p.x() * direction.x() + p.y() * direction.y();
0200       Real_v a       = paraboloid.fA * dirRho2;
0201       Real_v b       = Real_v(0.5) * direction.z() - paraboloid.fA * pDotV2D;
0202       Real_v c       = (paraboloid.fB + paraboloid.fA * p.Perp2() - p.z());
0203       Real_v d2      = b * b - a * c;
0204       done |= d2 < Real_v(0.);
0205       if (vecCore::MaskFull(done)) return;
0206 
0207       Real_v distParab = InfinityLength<Real_v>();
0208       vecCore__MaskedAssignFunc(distParab, !done && (b <= Real_v(0.)), (b - Sqrt(d2)) / NonZero(a));
0209       vecCore__MaskedAssignFunc(distParab, !done && (b > Real_v(0.)), (c / NonZero(b + Sqrt(d2))));
0210       Real_v zHit = p.z() + distParab * direction.z();
0211       vecCore::MaskedAssign(distance, Abs(zHit) <= paraboloid.fDz && distParab > Real_v(0.), distParab + offset);
0212     }
0213   }
0214 
0215   template <typename Real_v>
0216   VECGEOM_FORCE_INLINE VECCORE_ATT_HOST_DEVICE static void DistanceToOut(UnplacedStruct_t const &paraboloid,
0217                                                                          Vector3D<Real_v> const &point,
0218                                                                          Vector3D<Real_v> const &direction,
0219                                                                          Real_v const & /* stepMax */, Real_v &distance)
0220   {
0221 
0222     using Bool_v = vecCore::Mask_v<Real_v>;
0223 
0224     // setting distance to -1. for wrong side points
0225     distance = -1.;
0226     Bool_v done(false);
0227 
0228     // Outside Z range
0229     Bool_v outsideZ = Abs(point.z()) > paraboloid.fDz + kTolerance;
0230     done |= outsideZ;
0231     if (vecCore::MaskFull(done)) return;
0232 
0233     // Outside Parabolic surface
0234     Real_v rho2                                  = point.Perp2();
0235     Real_v paraRho2                              = paraboloid.fK1 * point.z() + paraboloid.fK2;
0236     Real_v value                                 = rho2 - paraRho2;
0237     Bool_v outsideParabolicSurfaceOuterTolerance = (value > kHalfTolerance);
0238     done |= outsideParabolicSurfaceOuterTolerance;
0239     if (vecCore::MaskFull(done)) return;
0240 
0241     // On Z Plane and moving outside;
0242     Bool_v isOnZPlaneAndMovingOutside = (IsOnZPlane<Real_v, true>(paraboloid, point) && direction.z() > Real_v(0.)) ||
0243                                         (IsOnZPlane<Real_v, false>(paraboloid, point) && direction.z() < Real_v(0.));
0244     vecCore__MaskedAssignFunc(distance, !done && isOnZPlaneAndMovingOutside, Real_v(0.));
0245     done |= isOnZPlaneAndMovingOutside;
0246     if (vecCore::MaskFull(done)) return;
0247 
0248     // On Parabolic Surface and moving outside
0249     Vector3D<Real_v> normal(point.x(), point.y(), Real_v(-paraboloid.fK1 * Real_v(0.5)));
0250     Bool_v isOnParabolicSurfaceAndMovingInside =
0251         value > -kTolerance && value < kTolerance && direction.Dot(normal) > Real_v(0.);
0252     vecCore__MaskedAssignFunc(distance, !done && isOnParabolicSurfaceAndMovingInside, Real_v(0.));
0253     done |= isOnParabolicSurfaceAndMovingInside;
0254     if (vecCore::MaskFull(done)) return;
0255 
0256     vecCore__MaskedAssignFunc(distance, !done, InfinityLength<Real_v>());
0257 
0258     Real_v distZ   = InfinityLength<Real_v>();
0259     Real_v dirZinv = Real_v(1.) / NonZero(direction.z());
0260 
0261     Bool_v dir_mask = direction.z() < 0;
0262     vecCore__MaskedAssignFunc(distZ, dir_mask, -(paraboloid.fDz + point.z()) * dirZinv);
0263     vecCore__MaskedAssignFunc(distZ, !dir_mask, (paraboloid.fDz - point.z()) * dirZinv);
0264 
0265     Real_v dirRho2 = direction.Perp2();
0266     Real_v pDotV2D = point.x() * direction.x() + point.y() * direction.y();
0267     Real_v a       = Real_v(paraboloid.fA * dirRho2);
0268     Real_v b       = Real_v(0.5) * direction.z() - Real_v(paraboloid.fA) * pDotV2D;
0269     Real_v c       = paraboloid.fB + paraboloid.fA * point.Perp2() - point.z();
0270     Real_v d2      = b * b - a * c;
0271 
0272     Real_v distParab = InfinityLength<Real_v>();
0273     vecCore__MaskedAssignFunc(distParab, d2 >= Real_v(0.) && (b > Real_v(0.)),
0274                               (b + Sqrt(d2)) * (Real_v(1.) / NonZero(a)));
0275     vecCore__MaskedAssignFunc(distParab, d2 >= Real_v(0.) && (b <= Real_v(0.)), (c / NonZero(b - Sqrt(d2))));
0276     distance = Min(distParab, distZ);
0277   }
0278 
0279   template <typename Real_v>
0280   VECGEOM_FORCE_INLINE VECCORE_ATT_HOST_DEVICE static void SafetyToIn(UnplacedStruct_t const &paraboloid,
0281                                                                       Vector3D<Real_v> const &point, Real_v &safety)
0282   {
0283 
0284     using Bool_v = vecCore::Mask_v<Real_v>;
0285     Real_v absZ  = Abs(point.z());
0286     Real_v safeZ = absZ - paraboloid.fDz;
0287 
0288     safety = -1.;
0289     Bool_v done(false);
0290     Bool_v insideZ = absZ < paraboloid.fDz - kTolerance;
0291 
0292     Real_v rho2                                 = point.Perp2();
0293     Real_v value                                = paraboloid.fA * rho2 + paraboloid.fB - point.z();
0294     Bool_v insideParabolicSurfaceOuterTolerance = (value < -kHalfTolerance);
0295     done |= (insideZ && insideParabolicSurfaceOuterTolerance);
0296     if (vecCore::MaskFull(done)) return;
0297 
0298     Bool_v onZPlane =
0299         Abs(Abs(point.z()) - paraboloid.fDz) < kTolerance &&
0300         (rho2 < Real_v(paraboloid.fRhi2 + kHalfTolerance) || rho2 < Real_v(paraboloid.fRlo2 + kHalfTolerance));
0301     vecCore__MaskedAssignFunc(safety, onZPlane, Real_v(0.));
0302     done |= onZPlane;
0303     if (vecCore::MaskFull(done)) return;
0304 
0305     Bool_v onParabolicSurface = value > -kTolerance && value < kTolerance;
0306     vecCore__MaskedAssignFunc(safety, !done && onParabolicSurface, Real_v(0.));
0307     done |= onParabolicSurface;
0308     if (vecCore::MaskFull(done)) return;
0309 
0310     vecCore__MaskedAssignFunc(safety, !done, InfinityLength<Real_v>());
0311 
0312     Real_v r0sq = (point.z() - paraboloid.fB) * paraboloid.fInvA;
0313 
0314     safety = safeZ;
0315 
0316     Bool_v underParaboloid = (r0sq < 0);
0317     done |= underParaboloid;
0318     if (vecCore::MaskFull(done)) return;
0319 
0320     Real_v safeR = InfinityLength<Real_v>();
0321     Real_v ro2   = point.x() * point.x() + point.y() * point.y();
0322     Real_v dr    = Sqrt(ro2) - Sqrt(r0sq);
0323 
0324     Bool_v drCloseToZero = (dr < Real_v(1.E-8));
0325     done |= drCloseToZero;
0326     if (vecCore::MaskFull(done)) return;
0327 
0328     // then go for the tangent
0329     Real_v talf = Real_v(-2.) * paraboloid.fA * Sqrt(r0sq);
0330     Real_v salf = talf / Sqrt(Real_v(1.) + talf * talf);
0331     safeR       = Abs(dr * salf);
0332 
0333     Real_v max_safety = Max(safeR, safeZ);
0334     vecCore::MaskedAssign(safety, !done, max_safety);
0335   }
0336 
0337   template <typename Real_v>
0338   VECGEOM_FORCE_INLINE VECCORE_ATT_HOST_DEVICE static void SafetyToOut(UnplacedStruct_t const &paraboloid,
0339                                                                        Vector3D<Real_v> const &point, Real_v &safety)
0340   {
0341 
0342     using Bool_v = vecCore::Mask_v<Real_v>;
0343 
0344     Real_v absZ = Abs(point.z());
0345     Real_v safZ = (paraboloid.fDz - absZ);
0346 
0347     safety = -1.;
0348     Bool_v done(false);
0349     Bool_v outsideZ = absZ > Real_v(paraboloid.fDz + kTolerance);
0350     done |= outsideZ;
0351     if (vecCore::MaskFull(done)) return;
0352 
0353     Real_v rho2                                  = point.Perp2();
0354     Real_v value                                 = paraboloid.fA * rho2 + paraboloid.fB - point.z();
0355     Bool_v outsideParabolicSurfaceOuterTolerance = (value > kHalfTolerance);
0356     done |= outsideParabolicSurfaceOuterTolerance;
0357     if (vecCore::MaskFull(done)) return;
0358 
0359     Bool_v onZPlane =
0360         Abs(Abs(point.z()) - paraboloid.fDz) < kTolerance &&
0361         (rho2 < Real_v(paraboloid.fRhi2 + kHalfTolerance) || rho2 < Real_v(paraboloid.fRlo2 + kHalfTolerance));
0362     vecCore__MaskedAssignFunc(safety, onZPlane, Real_v(0.));
0363     done |= onZPlane;
0364     if (vecCore::MaskFull(done)) return;
0365 
0366     Bool_v onSurface = value > -kTolerance && value < kTolerance;
0367     vecCore__MaskedAssignFunc(safety, !done && onSurface, Real_v(0.));
0368     done |= onSurface;
0369     if (vecCore::MaskFull(done)) return;
0370     Real_v r0sq = (point.z() - paraboloid.fB) * paraboloid.fInvA;
0371 
0372     safety = 0.;
0373 
0374     Bool_v closeToParaboloid = (r0sq < 0);
0375     done |= closeToParaboloid;
0376     if (vecCore::MaskFull(done)) return;
0377 
0378     Real_v safR = InfinityLength<Real_v>();
0379     Real_v ro2  = point.x() * point.x() + point.y() * point.y();
0380     Real_v z0   = paraboloid.fA * ro2 + paraboloid.fB;
0381     Real_v dr   = Sqrt(ro2) - Sqrt(r0sq); // avoid square root of a negative number
0382 
0383     Bool_v drCloseToZero = (dr > Real_v(-1.E-8));
0384     done |= drCloseToZero;
0385     if (vecCore::MaskFull(done)) return;
0386 
0387     Real_v dz = Abs(point.z() - z0);
0388     safR      = -dr * dz / Sqrt(dr * dr + dz * dz);
0389 
0390     Real_v min_safety = Min(safR, safZ);
0391     vecCore::MaskedAssign(safety, !done, min_safety);
0392   }
0393 
0394   template <typename Real_v>
0395   VECGEOM_FORCE_INLINE VECCORE_ATT_HOST_DEVICE static Vector3D<Real_v> NormalKernel(
0396       UnplacedStruct_t const &paraboloid, Vector3D<Real_v> const &point, typename vecCore::Mask_v<Real_v> &valid)
0397   {
0398     using Bool_v = vecCore::Mask_v<Real_v>;
0399 
0400     // used to store the normal that needs to be returned
0401     Vector3D<Real_v> normal(0., 0., 0.);
0402     Real_v nsurf(0.); // used to store the number of surfaces on which the point lie.
0403     // in case of paraboloid it can maximum go upto 2
0404 
0405     // The interface is in fact only scalar, so do the correct treatment for points on the axis of symmetry
0406     Vector3D<Real_v> normParabolic(0., 0., vecCore::math::Sign(-paraboloid.fA));
0407     Real_v r = point.Perp();
0408     if (r > kTolerance) {
0409       Real_v talf = -2 * paraboloid.fA * r;
0410       Real_v calf = 1. / Sqrt(1. + talf * talf);
0411       Real_v salf = talf * calf;
0412       Vector3D<Real_v> normParabolic((salf * point.x() / NonZero(point.Perp())),
0413                                      (salf * point.y() / NonZero(point.Perp())), calf);
0414       normParabolic.Normalize();
0415     }
0416 
0417     // Logic for Valid Normal i.e. when point is on the surface
0418     Bool_v isOnZPlane = IsOnZPlane<Real_v, true>(paraboloid, point) || IsOnZPlane<Real_v, false>(paraboloid, point);
0419     Bool_v isOnParabolicSurface = IsOnParabolicSurface<Real_v>(paraboloid, point);
0420 
0421     vecCore__MaskedAssignFunc(nsurf, isOnZPlane, nsurf + 1);
0422     vecCore__MaskedAssignFunc(normal[2], (IsOnZPlane<Real_v, true>(paraboloid, point)), Real_v(1.));
0423     vecCore__MaskedAssignFunc(normal[2], (IsOnZPlane<Real_v, false>(paraboloid, point)), Real_v(-1.));
0424 
0425     vecCore__MaskedAssignFunc(nsurf, isOnParabolicSurface, nsurf + 1);
0426     vecCore__MaskedAssignFunc(normal[0], isOnParabolicSurface, normal[0] - normParabolic[0]);
0427     vecCore__MaskedAssignFunc(normal[1], isOnParabolicSurface, normal[1] - normParabolic[1]);
0428     vecCore__MaskedAssignFunc(normal[2], isOnParabolicSurface, normal[2] - normParabolic[2]);
0429 
0430     valid = Bool_v(true);
0431     valid &= (nsurf > 0);
0432 
0433     if (vecCore::MaskFull(valid)) return normal.Normalized();
0434 
0435     // This block is used to calculate the Approximate normal
0436     Vector3D<Real_v> norm(0., 0., 0.); // used to store approximate normal
0437     vecCore__MaskedAssignFunc(norm[2], point.z() > Real_v(0.), Real_v(1.));
0438     vecCore__MaskedAssignFunc(norm[2], point.z() < Real_v(0.), Real_v(-1.));
0439 
0440     Real_v safz = paraboloid.fDz - Abs(point.z());
0441     Real_v safr = Abs(r - Sqrt((point.z() - paraboloid.fB) * paraboloid.fInvA));
0442     vecCore__MaskedAssignFunc(norm[0], safz >= Real_v(0.) && safr < safz, normParabolic.x());
0443     vecCore__MaskedAssignFunc(norm[1], safz >= Real_v(0.) && safr < safz, normParabolic.y());
0444     vecCore__MaskedAssignFunc(norm[2], safz >= Real_v(0.) && safr < safz, normParabolic.z());
0445 
0446     // If Valid is not set, that means the point is NOT on the surface,
0447     // So in that case we have to rely on Approximate normal
0448     vecCore__MaskedAssignFunc(normal[0], !valid, norm.x());
0449     vecCore__MaskedAssignFunc(normal[1], !valid, norm.y());
0450     vecCore__MaskedAssignFunc(normal[2], !valid, norm.z());
0451 
0452     return normal.Normalized();
0453   }
0454 };
0455 } // namespace VECGEOM_IMPL_NAMESPACE
0456 } // namespace vecgeom
0457 
0458 #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!