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