volumes/kernel/EllipticalConeImplementation.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 EllipticalCone
0006 /// @file volumes/kernel/EllipticalConeImplementation.h
0007 /// @author Raman Sehgal, Evgueni Tcherniaev
0008
0009 #ifndef VECGEOM_VOLUMES_KERNEL_ELLIPTICALCONEIMPLEMENTATION_H_
0010 #define VECGEOM_VOLUMES_KERNEL_ELLIPTICALCONEIMPLEMENTATION_H_
0011
0012 #include "VecGeom/base/Vector3D.h"
0013 #include "VecGeom/volumes/EllipticalConeStruct.h"
0014 #include "VecGeom/volumes/kernel/GenericKernels.h"
0015 #include <VecCore/VecCore>
0016
0017 #include <cstdio>
0018 #include <iomanip>
0019
0020 namespace vecgeom {
0021
0022 VECGEOM_DEVICE_FORWARD_DECLARE(struct EllipticalConeImplementation;);
0023 VECGEOM_DEVICE_DECLARE_CONV(struct, EllipticalConeImplementation);
0024
0025 inline namespace VECGEOM_IMPL_NAMESPACE {
0026
0027 class PlacedEllipticalCone;
0028 template <typename T>
0029 struct EllipticalConeStruct;
0030 class UnplacedEllipticalCone;
0031
0032 struct EllipticalConeImplementation {
0033
0034   using PlacedShape_t    = PlacedEllipticalCone;
0035   using UnplacedStruct_t = EllipticalConeStruct<Precision>;
0036   using UnplacedVolume_t = UnplacedEllipticalCone;
0037
0038   VECCORE_ATT_HOST_DEVICE
0039   static void PrintType()
0040   {
0041     //  printf("SpecializedEllipticalCone<%i, %i>", transCodeT, rotCodeT);
0042   }
0043
0044   template <typename Stream>
0045   static void PrintType(Stream &st, int transCodeT = translation::kGeneric, int rotCodeT = rotation::kGeneric)
0046   {
0047     st << "SpecializedEllipticalCone<" << transCodeT << "," << rotCodeT << ">";
0048   }
0049
0050   template <typename Stream>
0051   static void PrintImplementationType(Stream &st)
0052   {
0053     (void)st;
0054     // st << "EllipticalConeImplementation<" << transCodeT << "," << rotCodeT << ">";
0055   }
0056
0057   template <typename Stream>
0058   static void PrintUnplacedType(Stream &st)
0059   {
0060     (void)st;
0061     // TODO: this is wrong
0062     // st << "UnplacedEllipticalCone";
0063   }
0064
0065   template <typename Real_v, typename Bool_v>
0066   VECGEOM_FORCE_INLINE
0067   VECCORE_ATT_HOST_DEVICE
0068   static void Contains(UnplacedStruct_t const &ellipticalcone, Vector3D<Real_v> const &point, Bool_v &inside)
0069   {
0070     Bool_v unused, outside;
0071     GenericKernelForContainsAndInside<Real_v, Bool_v, false>(ellipticalcone, point, unused, outside);
0072     inside = !outside;
0073   }
0074
0075   // BIG QUESTION: DO WE WANT TO GIVE ALL 3 TEMPLATE PARAMETERS
0076   // -- OR -- DO WE WANT TO DEDUCE Bool_v, Index_t from Real_v???
0077   template <typename Real_v, typename Inside_t>
0078   VECGEOM_FORCE_INLINE
0079   VECCORE_ATT_HOST_DEVICE
0080   static void Inside(UnplacedStruct_t const &ellipticalcone, Vector3D<Real_v> const &point, Inside_t &inside)
0081   {
0082
0083     using Bool_v       = vecCore::Mask_v<Real_v>;
0084     using InsideBool_v = vecCore::Mask_v<Inside_t>;
0085     Bool_v completelyinside, completelyoutside;
0086     GenericKernelForContainsAndInside<Real_v, Bool_v, true>(ellipticalcone, point, completelyinside, completelyoutside);
0087     inside = EInside::kSurface;
0088     vecCore::MaskedAssign(inside, (InsideBool_v)completelyoutside, Inside_t(EInside::kOutside));
0089     vecCore::MaskedAssign(inside, (InsideBool_v)completelyinside, Inside_t(EInside::kInside));
0090   }
0091
0092   template <typename Real_v, typename Bool_v, bool ForInside>
0093   VECGEOM_FORCE_INLINE
0094   VECCORE_ATT_HOST_DEVICE
0095   static void GenericKernelForContainsAndInside(UnplacedStruct_t const &ellipticalcone, Vector3D<Real_v> const &point,
0096                                                 Bool_v &completelyinside, Bool_v &completelyoutside)
0097   {
0098     /* TODO : Logic to check where the point is inside or not.
0099     **
0100     ** if ForInside is false then it will only check if the point is outside,
0101     ** and is used by Contains function
0102     **
0103     ** if ForInside is true then it will check whether the point is inside or outside,
0104     ** and if neither inside nor outside then it is on the surface.
0105     ** and is used by Inside function
0106     */
0107     Real_v px     = point.x() * ellipticalcone.invDx;
0108     Real_v py     = point.y() * ellipticalcone.invDy;
0109     Real_v pz     = point.z();
0110     Real_v hp     = vecCore::math::Sqrt(px * px + py * py) + pz;
0111     Real_v ds     = (hp - ellipticalcone.fDz) * ellipticalcone.cosAxisMin;
0112     Real_v dz     = vecCore::math::Abs(pz) - ellipticalcone.fZCut;
0113     Real_v safety = vecCore::math::Max(ds, dz);
0114
0115     completelyoutside = safety > kHalfTolerance;
0116     if (ForInside) completelyinside = safety <= -kHalfTolerance;
0117     return;
0118   }
0119
0120   template <typename Real_v>
0121   VECGEOM_FORCE_INLINE
0122   VECCORE_ATT_HOST_DEVICE
0123   static void DistanceToIn(UnplacedStruct_t const &ellipticalcone, Vector3D<Real_v> const &point,
0124                            Vector3D<Real_v> const &direction, Real_v const & /*stepMax*/, Real_v &distance)
0125   {
0126     /* TODO :  Logic to calculate Distance from outside point to the EllipticalCone surface */
0127     using Bool_v       = vecCore::Mask_v<Real_v>;
0128     Real_v kTwoEpsilon = 2. * kEpsilon;
0129     distance           = Real_v(kInfLength);
0130     Real_v offset(0.);
0131     Vector3D<Real_v> p(point);
0132
0133     // Move point closer, if required
0134     Real_v Rfar2(1024. * ellipticalcone.fRsph * ellipticalcone.fRsph); // 1024 = 32 * 32
0135     vecCore__MaskedAssignFunc(offset, ((p.Mag2() > Rfar2) && (direction.Dot(p) < Real_v(0.))),
0136                               p.Mag() - Real_v(2.) * ellipticalcone.fRsph);
0137     p += offset * direction;
0138
0139     // Special cases to keep in mind:
0140     //   0) Point is on the surface and leaving the solid
0141     //   1) Trajectory is parallel to the surface (A = 0, single root at t = -C/2B)
0142     //   2) No intersection (D < 0) or touch (D < eps) with lateral surface
0143     //   3) Exception: when the trajectory traverses the apex (D < eps) and A < 0
0144     //      then always there is an intersection with the solid
0145
0146     // Set working variables, transform elliptican cone to cone
0147     Real_v px  = p.x() * ellipticalcone.invDx;
0148     Real_v py  = p.y() * ellipticalcone.invDy;
0149     Real_v pz  = p.z();
0150     Real_v pz0 = p.z() - ellipticalcone.fDz; // pz if apex would be in origin
0151     Real_v vx  = direction.x() * ellipticalcone.invDx;
0152     Real_v vy  = direction.y() * ellipticalcone.invDy;
0153     Real_v vz  = direction.z();
0154
0155     // Compute coefficients of the quadratic equation: A t^2 + 2B t + C = 0
0156     Real_v Ar = vx * vx + vy * vy;
0157     Real_v Br = px * vx + py * vy;
0158     Real_v Cr = px * px + py * py;
0159     // 1) Check if A = 0
0160     // If so, slightly modify vz to avoid degeneration of the quadratic equation
0161     // The magnitude of vz will be modified in a way that preserves correct behavior when 0) point is leaving the solid
0162     Real_v vzvz  = vz * vz;
0163     Bool_v tinyA = vecCore::math::Abs(Ar - vzvz) < kTwoEpsilon * vzvz;
0164     vecCore__MaskedAssignFunc(vz, tinyA, vz + vecCore::math::Abs(vz) * kTwoEpsilon);
0165
0166     Real_v Az = vz * vz;
0167     Real_v Bz = pz0 * vz;
0168     Real_v Cz = pz0 * pz0;
0169     Real_v A  = Ar - Az;
0170     Real_v B  = Br - Bz;
0171     Real_v B0 = Br - pz0 * direction.z(); // B calculated with original v.z()
0172     Real_v C  = Cr - Cz;
0173     Real_v D  = B * B - A * C;
0174
0175     // 0) Check if point is leaving the solid
0176     Real_v sfz = vecCore::math::Abs(pz) - ellipticalcone.fZCut;
0177     Real_v nz  = vecCore::math::Sqrt(Cr);
0178     Real_v sfr = (nz + pz0) * ellipticalcone.cosAxisMin;
0179     vecCore::MaskedAssign(nz, (vecCore::math::Abs(p.x()) + vecCore::math::Abs(p.y()) < Real_v(0.1) * kHalfTolerance),
0180                           Real_v(1.));       // point is on z-axis
0181     Real_v pzA = pz0 + ellipticalcone.dApex; // slightly shifted apex position for "flying away" check
0182     Bool_v done =
0183         (sfz >= -kHalfTolerance && pz * vz >= Real_v(0.)) || (sfr >= -kHalfTolerance && Br + nz * vz >= Real_v(0.)) ||
0184         (pz0 * ellipticalcone.cosAxisMin > -kHalfTolerance && (Cr - pzA * pzA) <= Real_v(0.) && A >= Real_v(0.));
0185
0186     // 2) Check if scratching (D < eps & A > 0) or no intersection (D < 0)
0187     // 3) if (D < eps & A < 0) then trajectory traverses the apex area - continue calculation
0188     vecCore__MaskedAssignFunc(D, (sfr <= Real_v(0.) && D < Real_v(0.)), Real_v(0.));
0189     done |= (D < Real_v(0.)) || ((D < kTwoEpsilon * B * B) && (A >= Real_v(0.)));
0190
0191     // Find intersection with Z planes
0192     Real_v invz  = Real_v(-1.) / NonZero(vz);
0193     Real_v dz    = vecCore::math::CopySign(Real_v(ellipticalcone.fZCut), invz);
0194     Real_v tzin  = (pz - dz) * invz;
0195     Real_v tzout = (pz + dz) * invz;
0196
0197     // Find roots of the quadratic equation
0198     Real_v tmp(0.), t1(0.), t2(0.);
0199     vecCore__MaskedAssignFunc(tmp, !done, -B - vecCore::math::CopySign(vecCore::math::Sqrt(D), B));
0200     vecCore__MaskedAssignFunc(t1, !done, tmp / A);
0201     vecCore__MaskedAssignFunc(t2, !done && tmp != 0, C / tmp);
0202     vecCore__MaskedAssignFunc(t2, !done && tinyA && B != Real_v(0.), -C / (Real_v(2.) * B0)); // A ~ 0, t = -C / 2B
0203     Real_v tmin = vecCore::math::Min(t1, t2);
0204     Real_v tmax = vecCore::math::Max(t1, t2);
0205
0206     // Set default - intersection with lower nappe (A > 0)
0207     Real_v trin  = tmin;
0208     Real_v trout = tmax;
0209     // Check if intersection with upper nappe only, return infinity
0210     done |= (A >= Real_v(0.) && pz0 + vz * tmin >= Real_v(0.));
0211
0212     // Check if intersection with both nappes (A < 0)
0213     vecCore__MaskedAssignFunc(trin, (!done && A < Real_v(0.)), Real_v(-kInfLength));
0214     vecCore__MaskedAssignFunc(trout, (!done && A < Real_v(0.)), Real_v(kInfLength));
0215     vecCore__MaskedAssignFunc(trin, (!done && A < Real_v(0.) && vz < Real_v(0.)), tmax);
0216     vecCore__MaskedAssignFunc(trout, (!done && A < Real_v(0.) && vz > Real_v(0.)), tmin);
0217
0218     // Set distance
0219     // No special check for inside points, distance for inside points will be negative
0220     Real_v tin  = vecCore::math::Max(tzin, trin);
0221     Real_v tout = vecCore::math::Min(tzout, trout);
0222     vecCore__MaskedAssignFunc(distance, !done && (tout - tin) >= kHalfTolerance, tin + offset);
0223   }
0224
0225   template <typename Real_v>
0226   VECGEOM_FORCE_INLINE
0227   VECCORE_ATT_HOST_DEVICE
0228   static void DistanceToOut(UnplacedStruct_t const &ellipticalcone, Vector3D<Real_v> const &point,
0229                             Vector3D<Real_v> const &direction, Real_v const & /* stepMax */, Real_v &distance)
0230   {
0231     /* TODO :  Logic to calculate Distance from inside point to the EllipticalCone surface */
0232     using Bool_v       = vecCore::Mask_v<Real_v>;
0233     Real_v kTwoEpsilon = 2. * kEpsilon;
0234     distance           = Real_v(0.);
0235
0236     // Special cases to keep in mind:
0237     //   0) Point is on the surface and leaving the solid
0238     //   1) Trajectory is parallel to the surface (A = 0, single root at t = -C/2B)
0239     //   2) No intersection (D < 0) or touch (D < eps) with lateral surface
0240     //   3) Exception: when the trajectory traverses the apex (D < eps) and A < 0
0241     //      then always there is an intersection with the solid
0242
0243     // Set working variables, transform elliptican cone to cone
0244     Real_v px  = point.x() * ellipticalcone.invDx;
0245     Real_v py  = point.y() * ellipticalcone.invDy;
0246     Real_v pz  = point.z();
0247     Real_v pz0 = pz - ellipticalcone.fDz; // pz if apex would be in origin
0248     Real_v hp  = vecCore::math::Sqrt(px * px + py * py) + pz;
0249     Real_v sfr = (hp - ellipticalcone.fDz) * ellipticalcone.cosAxisMin;
0250     Real_v sfz = vecCore::math::Abs(pz) - ellipticalcone.fZCut;
0251
0252     // Check if point is outside
0253     Bool_v outside = vecCore::math::Max(sfr, sfz) > kHalfTolerance;
0254     vecCore__MaskedAssignFunc(distance, outside, Real_v(-1.));
0255     Bool_v done = outside;
0256
0257     // Compute coefficients of the quadratic equation: A t^2 + 2B t + C = 0
0258     Real_v vx = direction.x() * ellipticalcone.invDx;
0259     Real_v vy = direction.y() * ellipticalcone.invDy;
0260     Real_v vz = direction.z();
0261     Real_v Ar = vx * vx + vy * vy;
0262     Real_v Br = px * vx + py * vy;
0263     Real_v Cr = px * px + py * py;
0264     // 1) Check if A = 0
0265     // If so, slightly modify vz to avoid degeneration of the quadratic equation
0266     // The magnitude of vz will be modified in a way that point is leaving the solid
0267     Bool_v tinyA = vecCore::math::Abs(Ar - vz * vz) < kTwoEpsilon * vz * vz;
0268     vecCore__MaskedAssignFunc(vz, tinyA, vz + vecCore::math::Abs(vz) * kTwoEpsilon);
0269
0270     Real_v Az = vz * vz;
0271     Real_v Bz = pz0 * vz;
0272     Real_v Cz = pz0 * pz0;
0273     Real_v A  = Ar - Az;
0274     Real_v B  = Br - Bz;
0275     Real_v B0 = Br - pz0 * direction.z(); // B calculated with original v.z()
0276     Real_v C  = Cr - Cz;
0277     Real_v D  = B * B - A * C;
0278     vecCore__MaskedAssignFunc(D, (sfr <= Real_v(0.) && D < Real_v(0.)), Real_v(0.));
0279
0280     // 2) Check if scratching (D < eps & A > 0) or no intersection (D < 0)
0281     // 3) if (D < eps & A < 0) then trajectory traverses the apex area - continue calculation
0282     done |= (D < Real_v(0.)) || (D < kTwoEpsilon * B * B && A >= Real_v(0.));
0283
0284     // Find intersection with Z planes
0285     Real_v tzout = kMaximum;
0286     vecCore__MaskedAssignFunc(tzout, vz != Real_v(0.),
0287                               (vecCore::math::CopySign(Real_v(ellipticalcone.fZCut), vz) - pz) / direction.z());
0288
0289     // Find roots of the quadratic equation
0290     Real_v tmp(0.), t1(0.), t2(0.);
0291     vecCore__MaskedAssignFunc(tmp, !done, -B - vecCore::math::CopySign(vecCore::math::Sqrt(D), B));
0292     vecCore__MaskedAssignFunc(t1, !done, tmp / A);
0293     vecCore__MaskedAssignFunc(t2, !done && tmp != Real_v(0.), C / tmp);
0294     vecCore__MaskedAssignFunc(t2, !done && tinyA && B0 != Real_v(0.), -C / (Real_v(2.) * B0)); // A ~ 0, t = -C / 2B
0295     Real_v tmin = vecCore::math::Min(t1, t2);
0296     Real_v tmax = vecCore::math::Max(t1, t2);
0297
0298     // Set default - intersection with lower nappe (A > 0)
0299     Real_v trout = tmax;
0300     // Check if intersection with upper nappe only or flying away, return 0
0301     done |= ((A >= Real_v(0.) && pz0 + vz * tmax >= Real_v(0.)) || (pz0 >= Real_v(0.) && vz >= Real_v(0.)));
0302
0303     // Check if intersection with both nappes (A < 0)
0304     vecCore__MaskedAssignFunc(trout, (!done && A < Real_v(0.)), Real_v(kInfLength));
0305     vecCore__MaskedAssignFunc(trout, (!done && A < Real_v(0.) && vz > Real_v(0.)), tmin);
0306
0307     // Set distance
0308     // No special check for inside points, distance for inside points will be negative
0309     vecCore__MaskedAssignFunc(distance, !done, vecCore::math::Min(tzout, trout));
0310   }
0311
0312   template <typename Real_v>
0313   VECGEOM_FORCE_INLINE
0314   VECCORE_ATT_HOST_DEVICE
0315   static void SafetyToIn(UnplacedStruct_t const &ellipticalcone, Vector3D<Real_v> const &point, Real_v &safety)
0316   {
0317     /* TODO :  Logic to calculate Safety from outside point to the EllipticalCone surface */
0318     Real_v px = point.x() * ellipticalcone.invDx;
0319     Real_v py = point.y() * ellipticalcone.invDy;
0320     Real_v pz = point.z();
0321     Real_v hp = vecCore::math::Sqrt(px * px + py * py) + pz;
0322     Real_v ds = (hp - ellipticalcone.fDz) * ellipticalcone.cosAxisMin;
0323     Real_v dz = vecCore::math::Abs(pz) - ellipticalcone.fZCut;
0324     safety    = vecCore::math::Max(ds, dz);
0325     vecCore::MaskedAssign(safety, vecCore::math::Abs(safety) <= kHalfTolerance, Real_v(0.));
0326   }
0327
0328   template <typename Real_v>
0329   VECGEOM_FORCE_INLINE
0330   VECCORE_ATT_HOST_DEVICE
0331   static void SafetyToOut(UnplacedStruct_t const &ellipticalcone, Vector3D<Real_v> const &point, Real_v &safety)
0332   {
0333     /* TODO :  Logic to calculate Safety from inside point to the EllipticalCone surface */
0334     Real_v px = point.x() * ellipticalcone.invDx;
0335     Real_v py = point.y() * ellipticalcone.invDy;
0336     Real_v pz = point.z();
0337     Real_v hp = vecCore::math::Sqrt(px * px + py * py) + pz;
0338     Real_v ds = (ellipticalcone.fDz - hp) * ellipticalcone.cosAxisMin;
0339     Real_v dz = ellipticalcone.fZCut - vecCore::math::Abs(pz);
0340     safety    = vecCore::math::Min(ds, dz);
0341     vecCore::MaskedAssign(safety, vecCore::math::Abs(safety) <= kHalfTolerance, Real_v(0.));
0342   }
0343
0344   template <typename Real_v>
0345   VECGEOM_FORCE_INLINE
0346   VECCORE_ATT_HOST_DEVICE
0347   static Vector3D<Real_v> NormalKernel(UnplacedStruct_t const &ellipticalcone, Vector3D<Real_v> const &point,
0348                                        typename vecCore::Mask_v<Real_v> &valid)
0349   {
0350     // Computes the normal on a surface and returns it as a unit vector
0351     //   In case if the point is further than kHalfTolerance from the surface, set valid=false
0352     //   Must return a valid vector (even if the point is not on the surface)
0353     //
0354     //   On an edge provide an average normal of the corresponding base and lateral surface
0355     Vector3D<Real_v> normal(0., 0., 0.);
0356     valid = true;
0357
0358     // Check z planes
0359     Real_v px = point.x();
0360     Real_v py = point.y();
0361     Real_v pz = point.z();
0362     Real_v dz = vecCore::math::Abs(pz) - ellipticalcone.fZCut;
0363     vecCore__MaskedAssignFunc(normal[2], vecCore::math::Abs(dz) <= kHalfTolerance, vecCore::math::Sign(pz));
0364
0365     // Check lateral surface
0366     Real_v nx = px * ellipticalcone.invDx * ellipticalcone.invDx;
0367     Real_v ny = py * ellipticalcone.invDy * ellipticalcone.invDy;
0368     Real_v nz = vecCore::math::Sqrt(px * nx + py * ny);
0369     vecCore__MaskedAssignFunc(nz, (nx * nx + ny * ny) == Real_v(0.), Real_v(1.)); // z-axis
0370     Vector3D<Real_v> nside(nx, ny, nz);
0371     Real_v ds = (nz + pz - ellipticalcone.fDz) * ellipticalcone.cosAxisMin;
0372     vecCore__MaskedAssignFunc(normal, vecCore::math::Abs(ds) <= kHalfTolerance, (normal + nside.Unit()).Unit());
0373
0374     // Check if done
0375     vecCore::Mask_v<Real_v> done = normal.Mag2() > Real_v(0.);
0376     if (vecCore::MaskFull(done)) return normal;
0377
0378     // Point is not on the surface - normally, this should never be
0379     // Return normal to the nearest surface
0380     vecCore__MaskedAssignFunc(valid, !done, false);
0381     vecCore__MaskedAssignFunc(normal[2], !done, vecCore::math::Sign(pz));
0382     vecCore__MaskedAssignFunc(normal, !done && ds > dz, nside.Unit());
0383     return normal;
0384   }
0385 };
0386 } // namespace VECGEOM_IMPL_NAMESPACE
0387 } // namespace vecgeom
0388
0389 #endif // VECGEOM_VOLUMES_KERNEL_ELLIPTICALCONEIMPLEMENTATION_H_