VecGeom/volumes/Quadrilaterals.h

0001 /// \file Quadrilaterals.h
0002 /// \author Johannes de Fine Licht (johannes.definelicht@cern.ch)
0003
0004 #ifndef VECGEOM_VOLUMES_QUADRILATERALS_H_
0005 #define VECGEOM_VOLUMES_QUADRILATERALS_H_
0006
0007 #include "VecGeom/base/Global.h"
0008
0009 #include "VecGeom/base/Array.h"
0010 #include "VecGeom/base/AOS3D.h"
0011 #include "VecGeom/base/SOA3D.h"
0012 #include "VecGeom/base/Vector3D.h"
0013 #include "VecGeom/base/RNG.h"
0014 #include "VecGeom/volumes/kernel/GenericKernels.h"
0015 #include "VecGeom/volumes/Planes.h"
0016
0017 // Switches on/off explicit vectorization of algorithms using Vc
0018 // #define VECGEOM_QUADRILATERAL_ACCELERATION --> now done in CMakeFile
0019
0020 namespace vecgeom {
0021
0022 VECGEOM_DEVICE_FORWARD_DECLARE(class Quadrilaterals;);
0023 VECGEOM_DEVICE_DECLARE_CONV(class, Quadrilaterals);
0024
0025 inline namespace VECGEOM_IMPL_NAMESPACE {
0026
0027 class Quadrilaterals {
0028
0029 private:
0030   Planes fPlanes;               ///< The planes in which the quadrilaterals lie.
0031   Planes fSideVectors[4];       ///< Vectors pointing from a side constructed from two
0032                                 ///  corners to the origin, equivalent to
0033                                 ///  normal x (c1 - c0)
0034                                 ///  Used to check if an intersection is in bounds.
0035   AOS3D<Precision> fCorners[4]; ///< Four corners of the quadrilaterals. Used
0036                                 ///  for bounds checking.
0037
0038 public:
0039   typedef Planes Sides_t[4];
0040   typedef AOS3D<Precision> Corners_t[4];
0041
0042   Quadrilaterals() = default;
0043
0044   VECCORE_ATT_HOST_DEVICE
0045   Quadrilaterals(int size, bool convex = true);
0046
0047   /// @brief Construct aligned data using specialized allocator.
0048   /// @details To allocate the full object, call this constructor with placement new.
0049   /// @param size Size of the internal arrays
0050   /// @param a Aligned allocator, pre-initialized to fit the content
0051   /// @param convex Convexity of the quadrilaterals
0052   VECCORE_ATT_HOST_DEVICE
0053   Quadrilaterals(int size, AlignedAllocator &a, bool convex = true);
0054
0055   VECCORE_ATT_HOST_DEVICE
0056   ~Quadrilaterals();
0057
0058   VECCORE_ATT_HOST_DEVICE
0059   Quadrilaterals(Quadrilaterals const &other);
0060
0061   VECCORE_ATT_HOST_DEVICE
0062   Quadrilaterals &operator=(Quadrilaterals const &other);
0063
0064   // returns the number of quadrilaterals ( planes ) stored in this container
0065   VECCORE_ATT_HOST_DEVICE
0066   VECGEOM_FORCE_INLINE
0067   int size() const;
0068
0069   VECCORE_ATT_HOST_DEVICE
0070   VECGEOM_FORCE_INLINE
0071   static size_t aligned_sizeof_data(const size_t initSize);
0072
0073   VECCORE_ATT_HOST_DEVICE
0074   VECGEOM_FORCE_INLINE
0075   Planes const &GetPlanes() const;
0076
0077   VECCORE_ATT_HOST_DEVICE
0078   VECGEOM_FORCE_INLINE
0079   SOA3D<Precision> const &GetNormals() const;
0080
0081   VECCORE_ATT_HOST_DEVICE
0082   VECGEOM_FORCE_INLINE
0083   Vector3D<Precision> GetNormal(int i) const;
0084
0085   VECCORE_ATT_HOST_DEVICE
0086   VECGEOM_FORCE_INLINE
0087   Array<Precision> const &GetDistances() const;
0088
0089   VECCORE_ATT_HOST_DEVICE
0090   VECGEOM_FORCE_INLINE
0091   Precision GetDistance(int i) const;
0092
0093   VECCORE_ATT_HOST_DEVICE
0094   VECGEOM_FORCE_INLINE
0095   Sides_t const &GetSideVectors() const;
0096
0097   VECCORE_ATT_HOST_DEVICE
0098   VECGEOM_FORCE_INLINE
0099   Corners_t const &GetCorners() const;
0100
0101   VECCORE_ATT_HOST_DEVICE
0102   VECGEOM_FORCE_INLINE
0103   Precision GetTriangleArea(int index, int iCorner1, int iCorner2) const;
0104
0105   VECCORE_ATT_HOST_DEVICE
0106   VECGEOM_FORCE_INLINE
0107   Precision GetQuadrilateralArea(int index) const;
0108
0109   inline Vector3D<Precision> GetPointOnTriangle(int index, int iCorner0, int iCorner1, int iCorner2) const;
0110
0111   inline Vector3D<Precision> GetPointOnFace(int index) const;
0112
0113   VECCORE_ATT_HOST_DEVICE
0114   VECGEOM_FORCE_INLINE
0115   bool RayHitsQuadrilateral(int index, Vector3D<Precision> const &intersection) const
0116   {
0117     bool valid = true;
0118     for (int j = 0; j < 4; ++j) {
0119       valid &=
0120           intersection.Dot(fSideVectors[j].GetNormal(index)) + fSideVectors[j].GetDistances()[index] >= -kTolerance;
0121       if (vecCore::MaskEmpty(valid)) break;
0122     }
0123     return valid;
0124   }
0125
0126   /// Sets the corners of a pre-existing quadrilateral.
0127   /// \param corner0 First corner in counterclockwise order.
0128   /// \param corner1 Second corner in counterclockwise order.
0129   /// \param corner2 Third corner in counterclockwise order.
0130   /// \param corner3 Fourth corner in counterclockwise order.
0131   VECCORE_ATT_HOST_DEVICE
0132   void Set(int index, Vector3D<Precision> const &corner0, Vector3D<Precision> const &corner1,
0133            Vector3D<Precision> const &corner2, Vector3D<Precision> const &corner3);
0134
0135   /// Flips the sign of the normal and distance of the specified quadrilateral.
0136   VECCORE_ATT_HOST_DEVICE
0137   void FlipSign(int index);
0138
0139   template <typename Real_v>
0140   VECGEOM_FORCE_INLINE VECCORE_ATT_HOST_DEVICE vecCore::Mask_v<Real_v> Contains(Vector3D<Real_v> const &point) const;
0141
0142   template <typename Real_v, typename Inside_v>
0143   VECGEOM_FORCE_INLINE VECCORE_ATT_HOST_DEVICE Inside_v Inside(Vector3D<Real_v> const &point) const;
0144
0145   template <typename Real_v, typename Inside_v>
0146   VECGEOM_FORCE_INLINE VECCORE_ATT_HOST_DEVICE Inside_v Inside(Vector3D<Real_v> const &point, int i) const;
0147
0148   template <typename Real_v, bool behindPlanesT>
0149   VECGEOM_FORCE_INLINE VECCORE_ATT_HOST_DEVICE Real_v DistanceToIn(Vector3D<Real_v> const &point,
0150                                                                    Vector3D<Real_v> const &direction) const;
0151
0152   template <typename Real_v>
0153   VECGEOM_FORCE_INLINE VECCORE_ATT_HOST_DEVICE Real_v DistanceToOut(Vector3D<Real_v> const &point,
0154                                                                     Vector3D<Real_v> const &direction, Precision zMin,
0155                                                                     Precision zMax) const;
0156
0157   template <typename Real_v>
0158   VECGEOM_FORCE_INLINE VECCORE_ATT_HOST_DEVICE Real_v DistanceToOut(Vector3D<Real_v> const &point,
0159                                                                     Vector3D<Real_v> const &direction) const;
0160
0161   /// \param index Quadrilateral to compute distance to.
0162   VECGEOM_FORCE_INLINE
0163   VECCORE_ATT_HOST_DEVICE
0164   Precision ScalarDistanceSquared(int index, Vector3D<Precision> const &point) const;
0165
0166   VECCORE_ATT_HOST_DEVICE
0167   void Print() const;
0168
0169 }; // end of class declaration
0170
0171 VECCORE_ATT_HOST_DEVICE
0172 VECGEOM_FORCE_INLINE
0173 int Quadrilaterals::size() const { return fPlanes.size(); }
0174
0175 VECCORE_ATT_HOST_DEVICE
0176 VECGEOM_FORCE_INLINE
0177 size_t Quadrilaterals::aligned_sizeof_data(const size_t initSize)
0178 {
0179   return (5 * Planes::aligned_sizeof_data(initSize) + 4 * AOS3D<Precision>::aligned_sizeof_data(initSize));
0180 }
0181
0182 VECCORE_ATT_HOST_DEVICE VECGEOM_FORCE_INLINE Planes const &Quadrilaterals::GetPlanes() const { return fPlanes; }
0183
0184 VECCORE_ATT_HOST_DEVICE
0185 SOA3D<Precision> const &Quadrilaterals::GetNormals() const { return fPlanes.GetNormals(); }
0186
0187 VECCORE_ATT_HOST_DEVICE
0188 Vector3D<Precision> Quadrilaterals::GetNormal(int i) const { return fPlanes.GetNormal(i); }
0189
0190 VECCORE_ATT_HOST_DEVICE
0191 Array<Precision> const &Quadrilaterals::GetDistances() const { return fPlanes.GetDistances(); }
0192
0193 VECCORE_ATT_HOST_DEVICE
0194 Precision Quadrilaterals::GetDistance(int i) const { return fPlanes.GetDistance(i); }
0195
0196 VECCORE_ATT_HOST_DEVICE
0197 Quadrilaterals::Sides_t const &Quadrilaterals::GetSideVectors() const { return fSideVectors; }
0198
0199 VECCORE_ATT_HOST_DEVICE
0200 Quadrilaterals::Corners_t const &Quadrilaterals::GetCorners() const { return fCorners; }
0201
0202 VECCORE_ATT_HOST_DEVICE
0203 Precision Quadrilaterals::GetTriangleArea(int index, int iCorner1, int iCorner2) const
0204 {
0205   Precision fArea          = 0.;
0206   Vector3D<Precision> vec1 = fCorners[iCorner1][index] - fCorners[0][index];
0207   Vector3D<Precision> vec2 = fCorners[iCorner2][index] - fCorners[0][index];
0208
0209   fArea = 0.5 * (vec1.Cross(vec2)).Mag();
0210   return fArea;
0211 }
0212
0213 VECCORE_ATT_HOST_DEVICE
0214 Precision Quadrilaterals::GetQuadrilateralArea(int index) const
0215 {
0216   Precision fArea = 0.;
0217
0218   fArea = GetTriangleArea(index, 1, 2) + GetTriangleArea(index, 2, 3);
0219   return fArea;
0220 }
0221
0222 Vector3D<Precision> Quadrilaterals::GetPointOnTriangle(int index, int iCorner0, int iCorner1, int iCorner2) const
0223 {
0224   Precision r1 = RNG::Instance().uniform(0.0, 1.0);
0225   Precision r2 = RNG::Instance().uniform(0.0, 1.0);
0226   if (r1 + r2 > 1.) {
0227     r1 = 1. - r1;
0228     r2 = 1. - r2;
0229   }
0230   Vector3D<Precision> vec1 = fCorners[iCorner1][index] - fCorners[iCorner0][index];
0231   Vector3D<Precision> vec2 = fCorners[iCorner2][index] - fCorners[iCorner0][index];
0232   return fCorners[iCorner0][index] + r1 * vec1 + r2 * vec2;
0233 }
0234
0235 Vector3D<Precision> Quadrilaterals::GetPointOnFace(int index) const
0236 {
0237
0238   // Avoid degenerated surfaces
0239   int nvert = 0;
0240   int iCorners[4];
0241   for (int i = 0; i < 4; ++i) {
0242     if ((fCorners[(i + 1) % 4][index] - fCorners[i % 4][index]).Mag2() > kTolerance) {
0243       iCorners[nvert++] = i;
0244     }
0245   }
0246   if (nvert == 1) return fCorners[0][index];
0247   if (nvert == 2) return GetPointOnTriangle(index, iCorners[0], iCorners[0], iCorners[1]);
0248   if (nvert == 3) return GetPointOnTriangle(index, iCorners[0], iCorners[1], iCorners[2]);
0249   Precision choice = RNG::Instance().uniform(0, 1.0);
0250   if (choice < 0.5) {
0251     return GetPointOnTriangle(index, 0, 1, 2);
0252   } else {
0253     return GetPointOnTriangle(index, 0, 2, 3);
0254   }
0255 }
0256
0257 template <typename Real_v>
0258 VECCORE_ATT_HOST_DEVICE vecCore::Mask_v<Real_v> Quadrilaterals::Contains(Vector3D<Real_v> const &point) const
0259 {
0260   return fPlanes.Contains<Real_v>(point);
0261 }
0262
0263 template <typename Real_v, typename Inside_v>
0264 VECCORE_ATT_HOST_DEVICE Inside_v Quadrilaterals::Inside(Vector3D<Real_v> const &point) const
0265 {
0266   return fPlanes.Inside<Real_v, Inside_v>(point);
0267 }
0268
0269 template <typename Real_v, typename Inside_v>
0270 VECCORE_ATT_HOST_DEVICE Inside_v Quadrilaterals::Inside(Vector3D<Real_v> const &point, int i) const
0271 {
0272   return fPlanes.Inside<Real_v, Inside_v>(point, i);
0273 }
0274
0275 namespace {
0276
0277 template <class Real_v>
0278 struct AcceleratedDistanceToIn {
0279   template <bool behindPlanesT>
0280   VECGEOM_FORCE_INLINE VECCORE_ATT_HOST_DEVICE static void VectorLoop(
0281       int & /*i*/, const int /*n*/, Planes const & /*planes*/, Planes const (& /*sideVectors*/)[4],
0282       Vector3D<Real_v> const & /*point*/, Vector3D<Real_v> const & /*direction*/, Real_v & /*distance*/)
0283   {
0284     // Do nothing if not scalar backend
0285     return;
0286   }
0287 };
0288
0289 #if defined(VECGEOM_VC) && defined(VECGEOM_QUADRILATERAL_ACCELERATION)
0290 template <>
0291 struct AcceleratedDistanceToIn<Precision> {
0292
0293   template <bool behindPlanesT>
0294   VECGEOM_FORCE_INLINE VECCORE_ATT_HOST_DEVICE static void VectorLoop(int &i, const int n, Planes const &planes,
0295                                                                       Planes const (&sideVectors)[4],
0296                                                                       Vector3D<Precision> const &point,
0297                                                                       Vector3D<Precision> const &direction,
0298                                                                       Precision &distance)
0299   {
0300
0301     // Explicitly vectorize over quadrilaterals using Vc
0302     for (; i <= n - kVectorSize; i += kVectorSize) {
0303       Vector3D<VcPrecision> plane(VcPrecision(planes.GetNormals().x() + i), VcPrecision(planes.GetNormals().y() + i),
0304                                   VcPrecision(planes.GetNormals().z() + i));
0305       VcPrecision dPlane(&planes.GetDistances()[0] + i);
0306       VcPrecision distanceTest = plane.Dot(point) + dPlane;
0307
0308       // Check if the point is in front of/behind the plane according to the template parameter
0309       VcBool valid = Flip<behindPlanesT>::FlipSign(distanceTest) > -kTolerance;
0310       if (vecCore::MaskEmpty(valid)) continue;
0311
0312       VcPrecision directionProjection = plane.Dot(direction);
0313       valid &= Flip<!behindPlanesT>::FlipSign(directionProjection) > 0;
0314       if (vecCore::MaskEmpty(valid)) continue;
0315       VcPrecision tiny = Vc::copysign(VcPrecision(1E-20), directionProjection);
0316       distanceTest /= -(directionProjection + tiny);
0317       Vector3D<VcPrecision> intersection = Vector3D<VcPrecision>(direction) * distanceTest + point;
0318
0319       for (int j = 0; j < 4; ++j) {
0320         Vector3D<VcPrecision> sideVector(VcPrecision(sideVectors[j].GetNormals().x() + i),
0321                                          VcPrecision(sideVectors[j].GetNormals().y() + i),
0322                                          VcPrecision(sideVectors[j].GetNormals().z() + i));
0323         VcPrecision dSide(&sideVectors[j].GetDistances()[i]);
0324         valid &= sideVector.Dot(intersection) + dSide >= -kTolerance;
0325         // Where is your god now
0326         if (vecCore::MaskEmpty(valid)) goto distanceToInVcContinueOuter;
0327       }
0328       // If a hit is found, the algorithm can return, since only one side can
0329       // be hit for a convex set of quadrilaterals
0330       distanceTest(!valid) = InfinityLength<Precision>();
0331       distance             = Max(distanceTest.min(), Precision(0.));
0332       i                    = n;
0333       return;
0334     // Continue label of outer loop
0335     distanceToInVcContinueOuter:;
0336     }
0337     return;
0338   }
0339 };
0340 #endif
0341
0342 } // End anonymous namespace
0343
0344 template <typename Real_v, bool behindPlanesT>
0345 VECCORE_ATT_HOST_DEVICE Real_v Quadrilaterals::DistanceToIn(Vector3D<Real_v> const &point,
0346                                                             Vector3D<Real_v> const &direction) const
0347 {
0348
0349   // Looks for the shortest distance to one of the quadrilaterals.
0350   // The algorithm projects the position and direction onto the plane of the
0351   // quadrilateral, determines the intersection point, then checks if this point
0352   // is within the bounds of the quadrilateral. There are many opportunities to
0353   // perform early returns along the way, and the speed of this algorithm relies
0354   // heavily on this property.
0355   //
0356   // The code below is optimized for the Polyhedron, and will return as soon as
0357   // a valid intersection is found, since only one intersection will ever occur
0358   // per Z-segment in Polyhedron case. If used in other contexts, a template
0359   // parameter would have to be added to make a distinction.
0360
0361   using Bool_v = vecCore::Mask_v<Real_v>;
0362
0363   Real_v bestDistance = InfinityLength<Real_v>();
0364
0365   int i       = 0;
0366   const int n = size();
0367   AcceleratedDistanceToIn<Real_v>::template VectorLoop<behindPlanesT>(i, n, fPlanes, fSideVectors, point, direction,
0368                                                                       bestDistance);
0369
0370   // TODO: IN CASE QUADRILATERALS ARE PERPENDICULAR TO Z WE COULD SAVE MANY DIVISIONS
0371   for (; i < n; ++i) {
0372     Vector3D<Precision> normal = fPlanes.GetNormal(i);
0373     Real_v distance            = point.Dot(normal) + fPlanes.GetDistance(i);
0374     // Check if the point is in front of/behind the plane according to the
0375     // template parameter
0376     Bool_v valid = Flip<behindPlanesT>::FlipSign(distance) > -kTolerance;
0377     if (vecCore::MaskEmpty(valid)) continue;
0378     Real_v directionProjection = direction.Dot(normal);
0379     valid &= Flip<!behindPlanesT>::FlipSign(directionProjection) > kTolerance;
0380     if (vecCore::MaskEmpty(valid)) continue;
0381     distance /= -(directionProjection + CopySign(Real_v(1E-20), directionProjection));
0382     Vector3D<Real_v> intersection = point + direction * distance;
0383     for (int j = 0; j < 4; ++j) {
0384       valid &= intersection.Dot(fSideVectors[j].GetNormal(i)) + fSideVectors[j].GetDistances()[i] >= -kTolerance;
0385       if (vecCore::MaskEmpty(valid)) break;
0386     }
0387     vecCore::MaskedAssign(bestDistance, valid, distance);
0388     // If all hits are found, the algorithm can return, since only one side can
0389     // be hit for a convex set of quadrilaterals
0390     if (vecCore::MaskFull(bestDistance < InfinityLength<Real_v>())) break;
0391   }
0392
0393   return Max(Real_v(0.), bestDistance);
0394 }
0395
0396 namespace {
0397
0398 template <typename Real_v>
0399 VECGEOM_FORCE_INLINE VECCORE_ATT_HOST_DEVICE void AcceleratedDistanceToOut(
0400     int & /*i*/, const int /*n*/, Planes const & /*planes*/, Planes const (& /*sideVectors*/)[4],
0401     const Precision /*zMin*/, const Precision /*zMax*/, Vector3D<Real_v> const & /*point*/,
0402     Vector3D<Real_v> const & /*direction*/, Real_v & /*distance*/)
0403 {
0404   // Do nothing if the backend is not scalar
0405   return;
0406 }
0407
0408 #if defined(VECGEOM_VC) && defined(VECGEOM_QUADRILATERAL_ACCELERATION)
0409 template <>
0410 VECGEOM_FORCE_INLINE VECCORE_ATT_HOST_DEVICE void AcceleratedDistanceToOut<Precision>(
0411     int &i, const int n, Planes const &planes, Planes const (&sideVectors)[4], const Precision zMin,
0412     const Precision zMax, Vector3D<Precision> const &point, Vector3D<Precision> const &direction, Precision &distance)
0413 {
0414
0415   // Explicitly vectorize over quadrilaterals using Vc
0416   for (; i <= n - kVectorSize; i += kVectorSize) {
0417     Vector3D<VcPrecision> plane(VcPrecision(planes.GetNormals().x() + i), VcPrecision(planes.GetNormals().y() + i),
0418                                 VcPrecision(planes.GetNormals().z() + i));
0419     VcPrecision dPlane(&planes.GetDistances()[0] + i);
0420     VcPrecision distanceTest = plane.Dot(point) + dPlane;
0421     // Check if the point is behind the plane
0422     VcBool valid = distanceTest < kTolerance;
0423     if (vecCore::MaskEmpty(valid)) continue;
0424     VcPrecision directionProjection = plane.Dot(direction);
0425     // Because the point is behind the plane, the direction must be along the
0426     // normal
0427     valid &= directionProjection > 0;
0428     if (vecCore::MaskEmpty(valid)) continue;
0429     distanceTest /= -NonZero(directionProjection);
0430     valid &= distanceTest < distance;
0431     if (vecCore::MaskEmpty(valid)) continue;
0432
0433     if (zMin == zMax) { // need a careful treatment in case of degenerate Z planes
0434       // in this case need proper hit detection
0435       Vector3D<VcPrecision> intersection = Vector3D<VcPrecision>(direction) * distanceTest + point;
0436       for (int j = 0; j < 4; ++j) {
0437         Vector3D<VcPrecision> sideVector(VcPrecision(sideVectors[j].GetNormals().x() + i),
0438                                          VcPrecision(sideVectors[j].GetNormals().y() + i),
0439                                          VcPrecision(sideVectors[j].GetNormals().z() + i));
0440         VcPrecision dSide(&sideVectors[j].GetDistances()[i]);
0441         valid &= sideVector.Dot(intersection) + dSide >= -kTolerance;
0442         // Where is your god now
0443         if (vecCore::MaskEmpty(valid)) goto distanceToOutVcContinueOuter;
0444       }
0445     } else {
0446       VcPrecision zProjection = distanceTest * direction[2] + point[2];
0447       valid &= zProjection >= zMin && zProjection < zMax;
0448     }
0449   distanceToOutVcContinueOuter:
0450     if (vecCore::MaskEmpty(valid)) continue;
0451     distanceTest(!valid) = InfinityLength<Precision>();
0452     distance             = distanceTest.min();
0453   }
0454   distance = Max(Precision(0.), distance);
0455   return;
0456 }
0457 #endif
0458
0459 } // End anonymous namespace
0460
0461 template <typename Real_v>
0462 VECCORE_ATT_HOST_DEVICE Real_v Quadrilaterals::DistanceToOut(Vector3D<Real_v> const &point,
0463                                                              Vector3D<Real_v> const &direction, Precision zMin,
0464                                                              Precision zMax) const
0465 {
0466
0467   // The below computes the distance to the quadrilaterals similar to
0468   // DistanceToIn, but is optimized for Polyhedron, and as such can assume that
0469   // the quadrilaterals form a convex shell, and that the shortest distance to
0470   // one of the quadrilaterals will indeed be an intersection. The exception to
0471   // this is if the point leaves the Z-bounds specified in the input parameters.
0472   // If used for another purpose than Polyhedron, DistanceToIn should be used if
0473   // the set of quadrilaterals is not convex.
0474
0475   using Bool_v = vecCore::Mask_v<Real_v>;
0476
0477   Real_v bestDistance = InfinityLength<Real_v>();
0478
0479   int i       = 0;
0480   const int n = size();
0481   AcceleratedDistanceToOut<Real_v>(i, n, fPlanes, fSideVectors, zMin, zMax, point, direction, bestDistance);
0482
0483   for (; i < n; ++i) {
0484     Vector3D<Precision> normal = fPlanes.GetNormal(i);
0485     Real_v distanceTest        = point.Dot(normal) + fPlanes.GetDistance(i);
0486     // Check if the point is behind the plane
0487     Bool_v valid = distanceTest < kTolerance;
0488     if (vecCore::MaskEmpty(valid)) continue;
0489     Real_v directionProjection = direction.Dot(normal);
0490     // Because the point is behind the plane, the direction must be along the
0491     // normal
0492     valid &= directionProjection > kTolerance;
0493     if (vecCore::MaskEmpty(valid)) continue;
0494     distanceTest /= -directionProjection;
0495     valid &= distanceTest < bestDistance && distanceTest > -kTolerance / directionProjection;
0496     if (vecCore::MaskEmpty(valid)) continue;
0497
0498     // this is a tricky test when zMin == zMax ( degenerate planes )
0499     if (zMin == zMax) {
0500       // in this case need proper hit detection
0501       // valid &= zProjection >= zMin-1E-10 && zProjection <= zMax+1E-10;
0502       Vector3D<Real_v> intersection = point + distanceTest * direction;
0503
0504       valid = RayHitsQuadrilateral(i, intersection);
0505
0506     } else {
0507       Real_v zProjection = point[2] + distanceTest * direction[2];
0508       valid &= (zProjection >= zMin - kTolerance) && (zProjection < zMax + kTolerance);
0509     }
0510     if (vecCore::MaskEmpty(valid)) continue;
0511     vecCore::MaskedAssign(bestDistance, valid, distanceTest);
0512   }
0513
0514   if (bestDistance > -kTolerance) bestDistance = Max(bestDistance, Precision(0.));
0515   return bestDistance;
0516 }
0517
0518 template <typename Real_v>
0519 VECCORE_ATT_HOST_DEVICE Real_v Quadrilaterals::DistanceToOut(Vector3D<Real_v> const &point,
0520                                                              Vector3D<Real_v> const &direction) const
0521 {
0522   return DistanceToOut<Real_v>(point, direction, -InfinityLength<Real_v>(), InfinityLength<Real_v>());
0523 }
0524
0525 VECCORE_ATT_HOST_DEVICE
0526 Precision Quadrilaterals::ScalarDistanceSquared(int i, Vector3D<Precision> const &point) const
0527 {
0528
0529   // This function is used by the safety algorithms to return the exact distance
0530   // to the quadrilateral specified.
0531   // The algorithm has three stages, trying first to return the shortest
0532   // distance to the plane, then to the closest line segment, then to the
0533   // closest corner.
0534   VECGEOM_ASSERT(i < size());
0535
0536   Vector3D<Precision> planeNormal = fPlanes.GetNormal(i);
0537   Precision distance              = point.Dot(planeNormal) + fPlanes.GetDistance(i);
0538   // Find the projection of the point on the quadrilateral "i". There was
0539   // a bug below by adding a distance along the plane normal, while the
0540   // correct version should subtract.
0541   Vector3D<Precision> intersection = point - distance * planeNormal;
0542
0543   bool withinBound[4];
0544   for (int j = 0; j < 4; ++j) {
0545     // TODO: check if this autovectorizes. Otherwise it should be explicitly
0546     //       vectorized.
0547     withinBound[j] = intersection[0] * fSideVectors[j].GetNormals().x(i) +
0548                          intersection[1] * fSideVectors[j].GetNormals().y(i) +
0549                          intersection[2] * fSideVectors[j].GetNormals().z(i) + fSideVectors[j].GetDistances()[i] >=
0550                      0;
0551   }
0552   if (withinBound[0] && withinBound[1] && withinBound[2] && withinBound[3]) {
0553     return distance * distance;
0554   }
0555
0556   // If the closest point is not on the plane itself, it must either be the
0557   // distance to the closest line segment or to the closest corner.
0558   // Since it is already known whether the point is to the left or right of
0559   // each line, only one side and its corners have to be checked.
0560   //
0561   //           above
0562   //   corner3_______corner0
0563   //         |       |
0564   //   left  |       |  right
0565   //         |_______|
0566   //   corner2       corner1
0567   //           below
0568   //
0569   // The assumption above that only one side has to be checked is not always
0570   // true. If withinBound is false for 2 connected segments making an angle > 90
0571   // and if the point projection on each of these segments is outside bounds,
0572   // it can happen that the closest point is not the same depending on the
0573   // checked segment. Something like below:
0574   //
0575   //    Point  x
0576   //
0577   //     corner3 ___corner0
0578   //            |   -
0579   //            |     -
0580   //            |_______- corner1
0581   //
0582   // If the "above" segment is checked, corner3 will be selected closest, which
0583   // is correct, but if the "right" segment gets checked, corner0 will be
0584   // wrongly selected.
0585
0586   Precision distancesq = InfinityLength<Precision>();
0587   for (int j = 0; j < 4; ++j) {
0588     if (!withinBound[j]) {
0589       distance = DistanceToLineSegmentSquared1<Precision>(fCorners[j][i], fCorners[(j + 1) % 4][i], point);
0590       if (distance < distancesq) distancesq = distance;
0591     }
0592   }
0593   return distancesq;
0594 }
0595
0596 std::ostream &operator<<(std::ostream &os, Quadrilaterals const &quads);
0597
0598 } // namespace VECGEOM_IMPL_NAMESPACE
0599
0600 } // namespace vecgeom
0601
0602 #endif // VECGEOM_VOLUMES_QUADRILATERALS_H_