Back to home page

EIC code displayed by LXR
 
 

    

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

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