volumes/kernel/TrapezoidImplementation.h

0001 //===-- kernel/TrapezoidImplementation.h ----------------------------*- C++ -*-===//
0002 //===--------------------------------------------------------------------------===//
0003 ///
0004 /// \file   kernel/TrapezoidImplementation.h
0005 /// \author Guilherme Lima (lima@fnal.gov)
0006 /// \brief This file implements the algorithms for the trapezoid
0007 ///
0008 /// Implementation details: initially based on USolids algorithms and vectorized types.
0009 ///
0010 //===--------------------------------------------------------------------------===//
0011 ///
0012 /// 140520  G. Lima   Created from USolids' UTrap algorithms
0013 /// 160722  G. Lima   Revision + moving to new backend structure
0014
0015 #ifndef VECGEOM_VOLUMES_KERNEL_TRAPEZOIDIMPLEMENTATION_H_
0016 #define VECGEOM_VOLUMES_KERNEL_TRAPEZOIDIMPLEMENTATION_H_
0017
0018 #include "VecGeom/base/Vector3D.h"
0019 #include "VecGeom/volumes/TrapezoidStruct.h"
0020 #include "VecGeom/volumes/kernel/GenericKernels.h"
0021 #include <VecCore/VecCore>
0022
0023 #include <cstdio>
0024
0025 namespace vecgeom {
0026
0027 VECGEOM_DEVICE_FORWARD_DECLARE(struct TrapezoidImplementation;);
0028 VECGEOM_DEVICE_DECLARE_CONV(struct, TrapezoidImplementation);
0029
0030 inline namespace VECGEOM_IMPL_NAMESPACE {
0031
0032 class PlacedTrapezoid;
0033 class UnplacedTrapezoid;
0034
0035 struct TrapezoidImplementation {
0036
0037   using PlacedShape_t    = PlacedTrapezoid;
0038   using UnplacedStruct_t = TrapezoidStruct<Precision>;
0039   using UnplacedVolume_t = UnplacedTrapezoid;
0040 #ifndef VECGEOM_PLANESHELL
0041   using TrapSidePlane = TrapezoidStruct<Precision>::TrapSidePlane;
0042 #endif
0043
0044 #ifndef VECGEOM_PLANESHELL
0045   template <typename Real_v>
0046   VECGEOM_FORCE_INLINE VECCORE_ATT_HOST_DEVICE static void EvaluateTrack(UnplacedStruct_t const &unplaced,
0047                                                                          Vector3D<Real_v> const &point,
0048                                                                          Vector3D<Real_v> const &dir, Real_v *pdist,
0049                                                                          Real_v *proj, Real_v *vdist)
0050   {
0051     TrapSidePlane const *fPlanes = unplaced.GetPlanes();
0052     // loop over side planes - find pdist,proj for each side plane
0053     // auto-vectorizable part of loop
0054     for (unsigned int i = 0; i < 4; ++i) {
0055       // Note: normal vector is pointing outside the volume (convention), therefore
0056       // pdist>0 if point is outside  and  pdist<0 means inside
0057       pdist[i] = fPlanes[i].fA * point.x() + fPlanes[i].fB * point.y() + fPlanes[i].fC * point.z() + fPlanes[i].fD;
0058
0059       // proj is projection of dir over the normal vector of side plane, hence
0060       // proj > 0 if pointing ~same direction as normal and proj<0 if ~opposite to normal
0061       proj[i] = fPlanes[i].fA * dir.x() + fPlanes[i].fB * dir.y() + fPlanes[i].fC * dir.z();
0062
0063       vdist[i] = -pdist[i] / NonZero(proj[i]);
0064     }
0065   }
0066 #endif
0067
0068   template <typename Real_v, typename Bool_v>
0069   VECGEOM_FORCE_INLINE VECCORE_ATT_HOST_DEVICE static void Contains(UnplacedStruct_t const &unplaced,
0070                                                                     Vector3D<Real_v> const &point, Bool_v &inside)
0071   {
0072     Bool_v unused(false), outside(false);
0073     GenericKernelForContainsAndInside<Real_v, Bool_v, false>(unplaced, point, unused, outside);
0074     inside = !outside;
0075   }
0076
0077   // BIG QUESTION: DO WE WANT TO GIVE ALL 3 TEMPLATE PARAMETERS
0078   // -- OR -- DO WE WANT TO DEDUCE Bool_v, Index_t from Real_v???
0079   template <typename Real_v, typename Inside_t>
0080   VECGEOM_FORCE_INLINE VECCORE_ATT_HOST_DEVICE static void Inside(UnplacedStruct_t const &unplaced,
0081                                                                   Vector3D<Real_v> const &point, Inside_t &inside)
0082   {
0083     using Bool_v = vecCore::Mask_v<Real_v>;
0084     Bool_v completelyInside(false), completelyOutside(false);
0085     GenericKernelForContainsAndInside<Real_v, Bool_v, true>(unplaced, point, completelyInside, completelyOutside);
0086
0087     using InsideBool_v = vecCore::Mask_v<Inside_t>;
0088     inside             = Inside_t(EInside::kSurface);
0089     vecCore__MaskedAssignFunc(inside, (InsideBool_v)completelyOutside, Inside_t(EInside::kOutside));
0090     vecCore__MaskedAssignFunc(inside, (InsideBool_v)completelyInside, Inside_t(EInside::kInside));
0091
0092     return;
0093   }
0094
0095   template <typename Real_v, typename Bool_v, bool ForInside>
0096   VECGEOM_FORCE_INLINE VECCORE_ATT_HOST_DEVICE static void GenericKernelForContainsAndInside(
0097       UnplacedStruct_t const &unplaced, Vector3D<Real_v> const &point, Bool_v &completelyInside,
0098       Bool_v &completelyOutside)
0099   {
0100     // z-region
0101     completelyOutside = Abs(point[2]) > MakePlusTolerant<true>(unplaced.fDz);
0102     if (ForInside) {
0103       completelyInside = Abs(point[2]) < MakeMinusTolerant<true>(unplaced.fDz);
0104     }
0105
0106 #ifdef VECGEOM_PLANESHELL
0107     unplaced.GetPlanes()->GenericKernelForContainsAndInside<Real_v, true>(point, completelyInside, completelyOutside);
0108 #else
0109     // here for PLANESHELL=OFF (disabled)
0110     TrapSidePlane const *fPlanes = unplaced.GetPlanes();
0111     Real_v dist[4];
0112     for (unsigned int i = 0; i < 4; ++i) {
0113       dist[i] = fPlanes[i].fA * point.x() + fPlanes[i].fB * point.y() + fPlanes[i].fC * point.z() + fPlanes[i].fD;
0114     }
0115
0116     for (unsigned int i = 0; i < 4; ++i) {
0117       // is it outside of this side plane?
0118       completelyOutside = completelyOutside || dist[i] > Real_v(MakePlusTolerant<true>(0.));
0119       if (ForInside) {
0120         completelyInside = completelyInside && dist[i] < Real_v(MakeMinusTolerant<true>(0.));
0121       }
0122       // if (vecCore::EarlyReturnMaxLength(completelyOutside,1) && vecCore::MaskFull(completelyOutside)) return;
0123     }
0124 #endif
0125
0126     return;
0127   }
0128
0129   ////////////////////////////////////////////////////////////////////////////
0130   //
0131   // Calculate distance to shape from outside - return kInfLength if no
0132   // intersection.
0133   //
0134   // ALGORITHM: For each component (z-planes, side planes), calculate
0135   // pair of minimum (smin) and maximum (smax) intersection values for
0136   // which the particle is in the extent of the shape.  The point of
0137   // entrance (exit) is found by the largest smin (smallest smax).
0138   //
0139   //  If largest smin > smallest smax, the trajectory does not reach
0140   //  inside the shape.
0141   //
0142   template <typename Real_v>
0143   VECGEOM_FORCE_INLINE VECCORE_ATT_HOST_DEVICE static void DistanceToIn(UnplacedStruct_t const &unplaced,
0144                                                                         Vector3D<Real_v> const &point,
0145                                                                         Vector3D<Real_v> const &dir,
0146                                                                         Real_v const &stepMax, Real_v &distance)
0147   {
0148     (void)stepMax;
0149     using Bool_v = vecCore::Mask_v<Real_v>;
0150     distance     = kInfLength;
0151
0152     //
0153     // Step 1: find range of distances along dir between Z-planes (smin, smax)
0154     //
0155
0156     // step 1.a) input particle is moving away --> return infinity
0157     Real_v signZdir = Sign(dir.z());
0158     Real_v max      = signZdir * unplaced.fDz - point.z(); // z-dist to farthest z-plane
0159
0160     // done = done || (dir.z()>0.0 && max < MakePlusTolerant<true>(0.));  // check if moving away towards +z
0161     // done = done || (dir.z()<0.0 && max > MakeMinusTolerant<true>(0.)); // check if moving away towards -z
0162     Bool_v done(signZdir * max < Real_v(MakePlusTolerant<true>(0.0))); // if outside + moving away towards +/-z
0163
0164     // if all particles moving away, we're done
0165     if (vecCore::EarlyReturnMaxLength(done, 1) && vecCore::MaskFull(done)) return;
0166
0167     // Step 1.b) General case:
0168     //   smax,smin are range of distances within z-range, taking direction into account.
0169     //   smin<smax - smax is positive, but smin may be either positive or negative
0170     Real_v invdir = Real_v(1.0) / NonZero(dir.z()); // convert distances from z to dir
0171     Real_v smax   = max * invdir;
0172     Real_v smin   = -(signZdir * unplaced.fDz + point.z()) * invdir;
0173
0174     //
0175     // Step 2: find distances for intersections with side planes.
0176     //
0177
0178 #ifdef VECGEOM_PLANESHELL
0179     // If disttoplanes is such that smin < dist < smax, then distance=disttoplanes
0180     Real_v disttoplanes = unplaced.GetPlanes()->DistanceToIn(point, dir, smin, smax);
0181     vecCore::MaskedAssign(distance, !done, disttoplanes);
0182
0183 #else
0184
0185     // here for VECGEOM_PLANESHELL_DISABLE
0186
0187     // loop over side planes - find pdist,Comp for each side plane
0188     Real_v pdist[4], comp[4], vdist[4];
0189     // EvaluateTrack<Real_v>(unplaced, point, dir, pdist, comp, vdist);
0190
0191     // auto-vectorizable part of loop
0192     TrapSidePlane const *fPlanes = unplaced.GetPlanes();
0193     for (unsigned int i = 0; i < 4; ++i) {
0194       // Note: normal vector is pointing outside the volume (convention), therefore
0195       // pdist>0 if point is outside  and  pdist<0 means inside
0196       pdist[i] = fPlanes[i].fA * point.x() + fPlanes[i].fB * point.y() + fPlanes[i].fC * point.z() + fPlanes[i].fD;
0197
0198       // Comp is projection of dir over the normal vector of side plane, hence
0199       // Comp > 0 if pointing ~same direction as normal and Comp<0 if ~opposite to normal
0200       comp[i] = fPlanes[i].fA * dir.x() + fPlanes[i].fB * dir.y() + fPlanes[i].fC * dir.z();
0201
0202       vdist[i] = -pdist[i] / NonZero(comp[i]);
0203     }
0204
0205     // check special cases
0206     for (int i = 0; i < 4; ++i) {
0207       // points fully outside a plane and moving away or parallel to that plane
0208       done = done || (pdist[i] > Real_v(MakePlusTolerant<true>(0.)) && comp[i] >= Real_v(0.));
0209       // points at a plane surface and exiting
0210       done = done || (pdist[i] > Real_v(MakeMinusTolerant<true>(0.)) && comp[i] > Real_v(0.));
0211     }
0212     // if all particles moving away, we're done
0213     if (vecCore::EarlyReturnMaxLength(done, 1) && vecCore::MaskFull(done)) return;
0214
0215     // this part does not auto-vectorize
0216     for (unsigned int i = 0; i < 4; ++i) {
0217       // if outside and moving away, return infinity
0218       Bool_v posPoint = pdist[i] > Real_v(MakeMinusTolerant<true>(0.));
0219       Bool_v posDir   = comp[i] > 0;
0220
0221       // check if trajectory will intercept plane within current range (smin,smax), otherwise track misses shape
0222       Bool_v interceptFromInside  = (!posPoint && posDir);
0223       Bool_v interceptFromOutside = (posPoint && !posDir);
0224
0225       //.. If dist is such that smin < dist < smax, then adjust either smin or smax
0226       vecCore__MaskedAssignFunc(smax, interceptFromInside && vdist[i] < smax, vdist[i]);
0227       vecCore__MaskedAssignFunc(smin, interceptFromOutside && vdist[i] > smin, vdist[i]);
0228     }
0229
0230     vecCore::MaskedAssign(distance, !done && smin <= smax, smin);
0231     vecCore__MaskedAssignFunc(distance, distance < Real_v(MakeMinusTolerant<true>(0.0)), Real_v(-1.));
0232 #endif
0233   }
0234
0235   template <typename Real_v>
0236   VECGEOM_FORCE_INLINE VECCORE_ATT_HOST_DEVICE static void DistanceToOut(UnplacedStruct_t const &unplaced,
0237                                                                          Vector3D<Real_v> const &point,
0238                                                                          Vector3D<Real_v> const &dir,
0239                                                                          Real_v const &stepMax, Real_v &distance)
0240   {
0241     (void)stepMax;
0242     using Bool_v = vecCore::Mask_v<Real_v>;
0243
0244     // step 0: if point is outside any plane --> return -1, otherwise initialize at Infinity
0245     Bool_v outside = Abs(point.z()) > MakePlusTolerant<true>(unplaced.fDz);
0246     distance       = vecCore::Blend(outside, Real_v(-1.0), InfinityLength<Real_v>());
0247     Bool_v done(outside);
0248     if (vecCore::EarlyReturnMaxLength(done, 1) && vecCore::MaskFull(done)) return;
0249
0250     //
0251     // Step 1: find range of distances along dir between Z-planes (smin, smax)
0252     //
0253
0254     Real_v distz = (Sign(dir.z()) * unplaced.fDz - point.z()) / NonZero(dir.z());
0255     vecCore__MaskedAssignFunc(distance, !done && Abs(dir.z()) /** maxXY*/ > kTolerance, distz);
0256
0257     //
0258     // Step 2: find distances for intersections with side planes.
0259     //
0260
0261 #ifdef VECGEOM_PLANESHELL
0262     Real_v disttoplanes = unplaced.GetPlanes()->DistanceToOut(point, dir);
0263     vecCore::MaskedAssign(distance, disttoplanes < distance, disttoplanes);
0264
0265 #else
0266     //=== Here for VECGEOM_PLANESHELL_DISABLE
0267
0268     // loop over side planes - find pdist,Proj for each side plane
0269     Real_v pdist[4], proj[4], vdist[4];
0270     // Real_v dist1(distance);
0271     // EvaluateTrack<Real_v>(unplaced, point, dir, pdist, proj, vdist);
0272
0273     TrapSidePlane const *fPlanes = unplaced.GetPlanes();
0274     for (unsigned int i = 0; i < 4; ++i) {
0275       // Note: normal vector is pointing outside the volume (convention), therefore
0276       // pdist>0 if point is outside  and  pdist<0 means inside
0277       pdist[i] = fPlanes[i].fA * point.x() + fPlanes[i].fB * point.y() + fPlanes[i].fC * point.z() + fPlanes[i].fD;
0278
0279       // Proj is projection of dir over the normal vector of side plane, hence
0280       // Proj > 0 if pointing ~same direction as normal and Proj<0 if pointing ~opposite to normal
0281       proj[i] = fPlanes[i].fA * dir.x() + fPlanes[i].fB * dir.y() + fPlanes[i].fC * dir.z();
0282
0283       vdist[i] = -pdist[i] / NonZero(proj[i]);
0284     }
0285
0286     // early return if point is outside of plane
0287     // for (unsigned int i = 0; i < 4; ++i) {
0288     //   done = done || (pdist[i] > MakePlusTolerant<true>(0.));
0289     // }
0290     // vecCore::MaskedAssign(dist1, done, Real_v(-1.0));
0291     // if (vecCore::EarlyReturnMaxLength(done,1) && vecCore::MaskFull(done)) return;
0292
0293     // std::cerr<<"=== point="<< point <<", dir="<< dir <<", distance="<< distance <<"\n";
0294     for (unsigned int i = 0; i < 4; ++i) {
0295       // if track is pointing towards plane and vdist<distance, then distance=vdist
0296       // vecCore__MaskedAssignFunc(dist1, !done && proj[i] > 0.0 && vdist[i] < dist1, vdist[i]);
0297       vecCore__MaskedAssignFunc(distance, pdist[i] > MakePlusTolerant<true>(0.), Real_v(-1.0));
0298       vecCore__MaskedAssignFunc(distance, proj[i] * unplaced.fDz > kTolerance && -Sign(pdist[i]) * vdist[i] < distance,
0299                                 -Sign(pdist[i]) * vdist[i]);
0300       // std::cerr<<"i="<< i <<", pdist="<< pdist[i] <<", proj="<< proj[i] <<", vdist="<< vdist[i] <<" --> dist="<<
0301       // dist1 <<", "<< distance <<"\n";
0302     }
0303 #endif
0304   }
0305
0306   template <typename Real_v>
0307   VECGEOM_FORCE_INLINE VECCORE_ATT_HOST_DEVICE static void SafetyToIn(UnplacedStruct_t const &unplaced,
0308                                                                       Vector3D<Real_v> const &point, Real_v &safety)
0309   {
0310     safety = Abs(point.z()) - unplaced.fDz;
0311
0312 #ifdef VECGEOM_PLANESHELL
0313     // Get safety over side planes
0314     unplaced.GetPlanes()->SafetyToIn(point, safety);
0315 #else
0316     // Loop over side planes
0317     TrapSidePlane const *fPlanes = unplaced.GetPlanes();
0318     Real_v dist[4];
0319     for (int i = 0; i < 4; ++i) {
0320       dist[i] = fPlanes[i].fA * point.x() + fPlanes[i].fB * point.y() + fPlanes[i].fC * point.z() + fPlanes[i].fD;
0321     }
0322
0323     // for (int i = 0; i < 4; ++i) {
0324     //   vecCore::MaskedAssign(safety, dist[i] > safety, dist[i]);
0325     // }
0326     Real_v safmax = Max(Max(dist[0], dist[1]), Max(dist[2], dist[3]));
0327     vecCore::MaskedAssign(safety, safmax > safety, safmax);
0328 #endif
0329   }
0330
0331   template <typename Real_v>
0332   VECGEOM_FORCE_INLINE VECCORE_ATT_HOST_DEVICE static void SafetyToOut(UnplacedStruct_t const &unplaced,
0333                                                                        Vector3D<Real_v> const &point, Real_v &safety)
0334   {
0335     // If point is outside (wrong-side) --> safety to negative value
0336     safety = unplaced.fDz - Abs(point.z());
0337
0338     // If all test points are outside, we're done
0339     // if (vecCore::EarlyReturnMaxLength(safety,1)) {
0340     //   if (vecCore::MaskFull(safety < kHalfTolerance)) return;
0341     // }
0342
0343 #ifdef VECGEOM_PLANESHELL
0344     // Get safety over side planes
0345     unplaced.GetPlanes()->SafetyToOut(point, safety);
0346 #else
0347     // Loop over side planes
0348     TrapSidePlane const *fPlanes = unplaced.GetPlanes();
0349
0350     // auto-vectorizable loop
0351     Real_v dist[4];
0352     for (int i = 0; i < 4; ++i) {
0353       dist[i] = -(fPlanes[i].fA * point.x() + fPlanes[i].fB * point.y() + fPlanes[i].fC * point.z() + fPlanes[i].fD);
0354     }
0355
0356     // unvectorizable loop
0357     // for (int i = 0; i < 4; ++i) {
0358     //   vecCore::MaskedAssign(safety, dist[i] < safety, dist[i]);
0359     // }
0360
0361     Real_v safmin = Min(Min(dist[0], dist[1]), Min(dist[2], dist[3]));
0362     vecCore::MaskedAssign(safety, safmin < safety, safmin);
0363 #endif
0364   }
0365
0366   template <typename Real_v>
0367   VECGEOM_FORCE_INLINE VECCORE_ATT_HOST_DEVICE static Vector3D<Real_v> NormalKernel(
0368       UnplacedStruct_t const &unplaced, Vector3D<Real_v> const &point, typename vecCore::Mask_v<Real_v> &valid)
0369   {
0370
0371     VECGEOM_CONST Precision delta = 1000. * kTolerance;
0372     Vector3D<Real_v> normal, cornerNormal;
0373     Real_v safety = InfinityLength<Real_v>();
0374     bool edge     = false;
0375
0376 #ifdef VECGEOM_PLANESHELL
0377     // Get normal from side planes -- PlaneShell case
0378     safety = unplaced.GetPlanes()->NormalKernel(point, normal, edge);
0379
0380 #else
0381     // Loop over side planes
0382     Vector3D<Real_v> cornerNormal;
0383     unsigned char surfaces = 0;
0384     // vectorizable loop
0385     TrapSidePlane const *fPlanes = unplaced.GetPlanes();
0386     Real_v dist[4];
0387     for (int i = 0; i < 4; ++i) {
0388       dist[i] = Abs(fPlanes[i].fA * point.x() + fPlanes[i].fB * point.y() + fPlanes[i].fC * point.z() + fPlanes[i].fD);
0389       // If closest update normal
0390       if (dist[i] < safety) {
0391         normal = unplaced.normals[i];
0392         safety = dist[i];
0393       }
0394       // If on surface add to separate vector
0395       if (dist[i] < kTolerance) {
0396         surfaces++;
0397         cornerNormal += unplaced.normals[i];
0398       }
0399     }
0400
0401     if (surfaces > 1) {
0402       // The point is on the edge - do not normalize the vector
0403       normal = cornerNormal;
0404       edge   = true;
0405     }
0406 #endif
0407
0408     // check if normal is valid w.r.t. z-planes, and define normals based on safety (see above)
0409     Real_v safz = vecCore::math::Abs(vecCore::math::Abs(point[2]) - unplaced.fDz);
0410     if (edge && safz < kTolerance) {
0411       // The point is on a corner
0412       normal += Vector3D<Real_v>(0., 0., Sign(point.z()));
0413     } else {
0414       if (safz < safety) {
0415         normal.Set(0., 0., Sign(point.z()));
0416         safety = safz;
0417       }
0418     }
0419     valid = Abs(safety) <= delta;
0420     // returned vector must be normalized
0421     normal.Normalize();
0422     return normal;
0423   }
0424 };
0425
0426 } // namespace VECGEOM_IMPL_NAMESPACE
0427 } // namespace vecgeom
0428
0429 #endif // VECGEOM_VOLUMES_KERNEL_TRAPEZOIDIMPLEMENTATION_H_