volumes/kernel/ParaboloidImplementation.h

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