Back to home page

EIC code displayed by LXR
 
 

    

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

0001 /*
0002  * GenTrapImplementation.h
0003  *
0004  *  Created on: Aug 2, 2014
0005  *      Author: swenzel
0006  *   Review/completion: Nov 4, 2015
0007  *      Author: mgheata
0008  */
0009 
0010 #ifndef VECGEOM_VOLUMES_KERNEL_GENTRAPIMPLEMENTATION_H_
0011 #define VECGEOM_VOLUMES_KERNEL_GENTRAPIMPLEMENTATION_H_
0012 
0013 #include "VecGeom/base/Global.h"
0014 
0015 #include "VecGeom/volumes/kernel/GenericKernels.h"
0016 #include "VecGeom/volumes/kernel/BoxImplementation.h"
0017 #include "VecGeom/volumes/UnplacedGenTrap.h"
0018 
0019 #include <iostream>
0020 
0021 namespace vecgeom {
0022 
0023 VECGEOM_DEVICE_FORWARD_DECLARE(struct GenTrapImplementation;);
0024 VECGEOM_DEVICE_DECLARE_CONV(struct, GenTrapImplementation);
0025 
0026 inline namespace VECGEOM_IMPL_NAMESPACE {
0027 
0028 class PlacedGenTrap;
0029 class UnplacedGenTrap;
0030 
0031 template <typename T>
0032 struct GenTrapStruct;
0033 
0034 struct GenTrapImplementation {
0035 
0036   using Vertex_t         = Vector3D<Precision>;
0037   using PlacedShape_t    = PlacedGenTrap;
0038   using UnplacedStruct_t = GenTrapStruct<Precision>;
0039   using UnplacedVolume_t = UnplacedGenTrap;
0040 
0041   VECCORE_ATT_HOST_DEVICE
0042   static void PrintType()
0043   {
0044     // printf("SpecializedGenTrap<%i, %i>", transCodeT, rotCodeT);
0045   }
0046 
0047   template <typename Stream>
0048   static void PrintType(Stream &s, int transCodeT = translation::kGeneric, int rotCodeT = rotation::kGeneric)
0049   {
0050     s << "SpecializedGenTrap<" << transCodeT << "," << rotCodeT << ">";
0051   }
0052 
0053   template <typename Stream>
0054   static void PrintImplementationType(Stream & /*s*/)
0055   {
0056     // s << "GenTrapImplementation<" << transCodeT << "," << rotCodeT << ">";
0057   }
0058 
0059   template <typename Stream>
0060   static void PrintUnplacedType(Stream & /*s*/)
0061   {
0062     // s << "UnplacedGenTrap";
0063   }
0064 
0065   template <typename Real_v, typename Bool_v>
0066   VECGEOM_FORCE_INLINE
0067   VECCORE_ATT_HOST_DEVICE
0068   static void Contains(UnplacedStruct_t const &unplaced, Vector3D<Real_v> const &point, Bool_v &inside);
0069 
0070   template <typename Real_v, typename Bool_v>
0071   VECGEOM_FORCE_INLINE
0072   VECCORE_ATT_HOST_DEVICE
0073   static void ContainsGeneric(UnplacedStruct_t const &unplaced, Vector3D<Real_v> const &point, Bool_v &inside);
0074 
0075   template <typename Real_v, typename Inside_t>
0076   VECGEOM_FORCE_INLINE
0077   VECCORE_ATT_HOST_DEVICE
0078   static void Inside(UnplacedStruct_t const &unplaced, Vector3D<Real_v> const &point, Inside_t &inside);
0079 
0080   template <typename Real_v, typename Inside_t>
0081   VECGEOM_FORCE_INLINE
0082   VECCORE_ATT_HOST_DEVICE
0083   static void InsideGeneric(UnplacedStruct_t const &unplaced, Vector3D<Real_v> const &point, Inside_t &inside);
0084 
0085   template <typename Real_v, bool ForInside>
0086   VECGEOM_FORCE_INLINE
0087   VECCORE_ATT_HOST_DEVICE
0088   static void GenericKernelForContainsAndInside(UnplacedStruct_t const &unplaced, Vector3D<Real_v> const &point,
0089                                                 vecCore::Mask_v<Real_v> &completelyinside,
0090                                                 vecCore::Mask_v<Real_v> &completelyoutside);
0091 
0092   template <typename Real_v>
0093   VECGEOM_FORCE_INLINE
0094   VECCORE_ATT_HOST_DEVICE
0095   static void DistanceToIn(UnplacedStruct_t const &unplaced, Vector3D<Real_v> const &point,
0096                            Vector3D<Real_v> const &direction, Real_v const &stepMax, Real_v &distance);
0097 
0098   template <typename Real_v>
0099   VECGEOM_FORCE_INLINE
0100   VECCORE_ATT_HOST_DEVICE
0101   static void DistanceToInGeneric(UnplacedStruct_t const &unplaced, Vector3D<Real_v> const &point,
0102                                   Vector3D<Real_v> const &direction, Real_v const &stepMax, Real_v &distance);
0103 
0104   template <typename Real_v>
0105   VECGEOM_FORCE_INLINE
0106   VECCORE_ATT_HOST_DEVICE
0107   static void DistanceToOut(UnplacedStruct_t const &unplaced, Vector3D<Real_v> const &point,
0108                             Vector3D<Real_v> const &direction, Real_v const &stepMax, Real_v &distance);
0109 
0110   template <typename Real_v>
0111   VECGEOM_FORCE_INLINE
0112   VECCORE_ATT_HOST_DEVICE
0113   static void DistanceToOutGeneric(UnplacedStruct_t const &unplaced, Vector3D<Real_v> const &point,
0114                                    Vector3D<Real_v> const &direction, Real_v const &stepMax, Real_v &distance);
0115 
0116   template <typename Real_v>
0117   VECGEOM_FORCE_INLINE
0118   VECCORE_ATT_HOST_DEVICE
0119   static void SafetyToIn(UnplacedStruct_t const &unplaced, Vector3D<Real_v> const &point, Real_v &safety);
0120 
0121   template <typename Real_v>
0122   VECGEOM_FORCE_INLINE
0123   VECCORE_ATT_HOST_DEVICE
0124   static void SafetyToInGeneric(UnplacedStruct_t const &unplaced, Vector3D<Real_v> const &point, Real_v &safety);
0125 
0126   template <typename Real_v>
0127   VECGEOM_FORCE_INLINE
0128   VECCORE_ATT_HOST_DEVICE
0129   static void SafetyToOut(UnplacedStruct_t const &unplaced, Vector3D<Real_v> const &point, Real_v &safety);
0130 
0131   template <typename Real_v>
0132   VECGEOM_FORCE_INLINE
0133   VECCORE_ATT_HOST_DEVICE
0134   static void SafetyToOutGeneric(UnplacedStruct_t const &unplaced, Vector3D<Real_v> const &point, Real_v &safety);
0135 
0136   template <typename Real_v, typename Bool_v>
0137   VECGEOM_FORCE_INLINE
0138   VECCORE_ATT_HOST_DEVICE
0139   static void NormalKernel(UnplacedStruct_t const &unplaced, Vector3D<Real_v> const &point, Vector3D<Real_v> &normal,
0140                            Bool_v &valid);
0141 
0142   template <class Real_v>
0143   VECCORE_ATT_HOST_DEVICE
0144   static void GetClosestEdge(Vector3D<Real_v> const &point, Real_v vertexX[4], Real_v vertexY[4], Real_v &iseg,
0145                              Real_v &fraction);
0146 
0147   template <typename Real_v>
0148   VECGEOM_FORCE_INLINE
0149   VECCORE_ATT_HOST_DEVICE
0150   static vecCore::Mask_v<Real_v> IsInTopOrBottomPolygon(UnplacedStruct_t const &unplaced, Real_v const &pointx,
0151                                                         Real_v const &pointy, vecCore::Mask_v<Real_v> top);
0152 }; // End struct GenTrapImplementation
0153 
0154 //********************************
0155 //**** implementations start here
0156 //********************************/
0157 
0158 template <typename Real_v>
0159 struct FillPlaneDataHelper {
0160   VECCORE_ATT_HOST_DEVICE
0161   VECGEOM_FORCE_INLINE
0162   static void FillPlaneData(GenTrapStruct<Precision> const &unplaced, Real_v &cornerx, Real_v &cornery, Real_v &deltax,
0163                             Real_v &deltay, vecCore::Mask_v<Real_v> const &top, int edgeindex)
0164   {
0165 
0166     // no vectorized data lookup for SIMD
0167     // need to fill the SIMD types individually
0168 
0169     for (size_t i = 0; i < vecCore::VectorSize<Real_v>(); ++i) {
0170       int index  = edgeindex + top[i] * 4;
0171       deltax[i]  = unplaced.fDeltaX[index];
0172       deltay[i]  = unplaced.fDeltaY[index];
0173       cornerx[i] = unplaced.fVerticesX[index];
0174       cornery[i] = unplaced.fVerticesY[index];
0175     }
0176   }
0177 };
0178 
0179 //______________________________________________________________________________
0180 /** @brief A partial template specialization for nonSIMD cases (scalar, cuda, ... ) */
0181 template <>
0182 struct FillPlaneDataHelper<Precision> {
0183   VECCORE_ATT_HOST_DEVICE
0184   VECGEOM_FORCE_INLINE
0185   static void FillPlaneData(GenTrapStruct<Precision> const &unplaced, Precision &cornerx, Precision &cornery,
0186                             Precision &deltax, Precision &deltay, bool const &top, int edgeindex)
0187   {
0188     int index = edgeindex + top * 4;
0189     deltax    = unplaced.fDeltaX[index];
0190     deltay    = unplaced.fDeltaY[index];
0191     cornerx   = unplaced.fVerticesX[index];
0192     cornery   = unplaced.fVerticesY[index];
0193   }
0194 };
0195 
0196 //______________________________________________________________________________
0197 template <typename Real_v>
0198 VECGEOM_FORCE_INLINE
0199 VECCORE_ATT_HOST_DEVICE
0200 vecCore::Mask_v<Real_v> GenTrapImplementation::IsInTopOrBottomPolygon(UnplacedStruct_t const &unplaced,
0201                                                                       Real_v const &pointx, Real_v const &pointy,
0202                                                                       vecCore::Mask_v<Real_v> top)
0203 {
0204   // optimized "inside" check for top or bottom z-surfaces
0205   // this is a bit tricky if different tracks check different planes
0206   // ( for example in case of Backend = Vc when top is mixed )
0207   // ( this is because vector data lookup is tricky )
0208 
0209   // stripped down version of the Contains kernel ( not yet shared with that kernel )
0210   // std::cerr << "IsInTopOrBottom: pointx: " << pointx << "  pointy: " << pointy << "  top: " << top << "\n";
0211 
0212   using Bool_v             = vecCore::Mask_v<Real_v>;
0213   Bool_v completelyoutside = Bool_v(false);
0214   Bool_v degenerate        = Bool_v(true);
0215   for (int i = 0; i < 4; ++i) {
0216     Real_v deltaX;
0217     Real_v deltaY;
0218     Real_v cornerX;
0219     Real_v cornerY;
0220 
0221     // thats the only place where scalar and vector code diverge
0222     // IsSIMD misses...replaced with early_returns
0223     FillPlaneDataHelper<Real_v>::FillPlaneData(unplaced, cornerX, cornerY, deltaX, deltaY, top, i);
0224 
0225     // std::cerr << i << " CORNERS " << cornerX << " " << cornerY << " " << deltaX << " " << deltaY << "\n";
0226 
0227     Real_v cross = (pointx - cornerX) * deltaY;
0228     cross -= (pointy - cornerY) * deltaX;
0229     degenerate =
0230         degenerate && (deltaX < Real_v(MakePlusTolerant<true>(0.))) && (deltaY < Real_v(MakePlusTolerant<true>(0.)));
0231     completelyoutside = completelyoutside || (cross < Real_v(MakeMinusTolerantCrossProduct<true,Real_v>(0., Abs(deltaX) + Abs(deltaY))));
0232     // if (vecCore::EarlyReturnAllowed()) {
0233     if (vecCore::MaskFull(completelyoutside)) {
0234       return Bool_v(false); 
0235     }
0236     // }
0237   }
0238   completelyoutside = completelyoutside || degenerate;
0239   return (!completelyoutside);
0240 }
0241 
0242 //______________________________________________________________________________
0243 template <typename Real_v, bool ForInside>
0244 VECGEOM_FORCE_INLINE
0245 VECCORE_ATT_HOST_DEVICE
0246 void GenTrapImplementation::GenericKernelForContainsAndInside(UnplacedStruct_t const &unplaced,
0247                                                               Vector3D<Real_v> const &point,
0248                                                               vecCore::Mask_v<Real_v> &completelyinside,
0249                                                               vecCore::Mask_v<Real_v> &completelyoutside)
0250 {
0251 
0252   using Bool_v                        = vecCore::Mask_v<Real_v>;
0253   VECGEOM_CONST Precision tolerancesq = 10000. * kTolerance * kTolerance;
0254   // Local point has to be translated in the bbox local frame.
0255   Vector3D<Real_v> halfsize(unplaced.fBBdimensions[0], unplaced.fBBdimensions[1], unplaced.fBBdimensions[2]);
0256   BoxImplementation::GenericKernelForContainsAndInside<Real_v, ForInside>(halfsize, point - unplaced.fBBorigin,
0257                                                                           completelyinside, completelyoutside);
0258   //  if (vecCore::EarlyReturnAllowed()) {
0259   if (vecCore::MaskFull(completelyoutside)) {
0260     return;
0261   }
0262   //  }
0263 
0264   // analyse z
0265   Real_v cf = unplaced.fHalfInverseDz * (unplaced.fDz - point.z());
0266   // analyse if x-y coordinates of point are within polygon at z-height
0267 
0268   //  loop over edges connecting points i with i+4
0269   Real_v vertexX[4];
0270   Real_v vertexY[4];
0271   // vectorizes for scalar backend
0272   for (int i = 0; i < 4; i++) {
0273     // calculate x-y positions of vertex i at this z-height
0274     vertexX[i] = unplaced.fVerticesX[i + 4] + cf * unplaced.fConnectingComponentsX[i];
0275     vertexY[i] = unplaced.fVerticesY[i + 4] + cf * unplaced.fConnectingComponentsY[i];
0276   }
0277 
0278   for (int i = 0; i < 4; i++) {
0279     // this is based on the following idea:
0280     // we decided for each edge whether the point is above or below the
0281     // 2d line defined by that edge
0282     // In fact, this calculation is part of the calculation of the distance
0283     // of point to that line which is a cross product. In this case it is
0284     // an embedded cross product of 2D vectors in 3D. The resulting vector always points
0285     // in z-direction whose z-magnitude is directly related to the distance.
0286     // see, e.g.,  http://geomalgorithms.com/a02-_lines.html
0287     if (unplaced.fDegenerated[i]) continue;
0288     int j         = (i + 1) % 4;
0289     Real_v DeltaX = vertexX[j] - vertexX[i];
0290     Real_v DeltaY = vertexY[j] - vertexY[i];
0291     Real_v cross  = (point.x() - vertexX[i]) * DeltaY - (point.y() - vertexY[i]) * DeltaX;
0292     if (ForInside) {
0293       Bool_v onsurf     = (cross * cross < tolerancesq * (DeltaX * DeltaX + DeltaY * DeltaY));
0294       completelyoutside = completelyoutside || (((cross < Real_v(MakeMinusTolerant<ForInside>(0.))) && (!onsurf)));
0295       completelyinside  = completelyinside && (cross > Real_v(MakePlusTolerant<ForInside>(0.))) && (!onsurf);
0296     } else {
0297       completelyoutside = completelyoutside || (cross < Real_v(MakeMinusTolerant<ForInside>(0.)));
0298     }
0299 
0300     //    if (vecCore::EarlyReturnAllowed()) {
0301     if (vecCore::MaskFull(completelyoutside)) {
0302       return;
0303     }
0304     //    }
0305   }
0306 }
0307 
0308 //______________________________________________________________________________
0309 template <typename Real_v, typename Bool_v>
0310 VECGEOM_FORCE_INLINE
0311 VECCORE_ATT_HOST_DEVICE
0312 void GenTrapImplementation::ContainsGeneric(UnplacedStruct_t const &unplaced, Vector3D<Real_v> const &point,
0313                                             Bool_v &inside)
0314 {
0315   // Generic implementation for contains
0316   Bool_v unused;
0317   Bool_v outside;
0318   GenericKernelForContainsAndInside<Real_v, false>(unplaced, point, unused, outside);
0319   inside = !outside;
0320 }
0321 
0322 //______________________________________________________________________________
0323 template <typename Real_v, typename Bool_v>
0324 VECGEOM_FORCE_INLINE
0325 VECCORE_ATT_HOST_DEVICE
0326 void GenTrapImplementation::Contains(UnplacedStruct_t const &unplaced, Vector3D<Real_v> const &point, Bool_v &inside)
0327 {
0328   GenTrapImplementation::ContainsGeneric<Real_v, Bool_v>(unplaced, point, inside);
0329 }
0330 
0331 //______________________________________________________________________________
0332 template <>
0333 VECGEOM_FORCE_INLINE
0334 VECCORE_ATT_HOST_DEVICE
0335 void GenTrapImplementation::Contains<Precision, bool>(UnplacedStruct_t const &unplaced,
0336                                                       Vector3D<Precision> const &point, bool &inside)
0337 {
0338 // Scalar specialization for Contains function
0339 #ifndef VECCORE_CUDA
0340   // The tessellated section helper is built only for the planar cases
0341   if (unplaced.fTslHelper) {
0342     // Check Z range
0343     inside = false;
0344     if (vecCore::math::Abs(point.z()) > unplaced.fDz) return;
0345     inside = unplaced.fTslHelper->Contains(point);
0346     return;
0347   }
0348 #endif
0349   GenTrapImplementation::ContainsGeneric<Precision, bool>(unplaced, point, inside);
0350 }
0351 
0352 //______________________________________________________________________________
0353 template <typename Real_v, typename Inside_t>
0354 VECGEOM_FORCE_INLINE
0355 VECCORE_ATT_HOST_DEVICE
0356 void GenTrapImplementation::InsideGeneric(UnplacedStruct_t const &unplaced, Vector3D<Real_v> const &point,
0357                                           Inside_t &inside)
0358 {
0359 
0360   using Bool_v       = vecCore::Mask_v<Real_v>;
0361   using InsideBool_v = vecCore::Mask_v<Inside_t>;
0362   Bool_v completelyinside;
0363   Bool_v completelyoutside;
0364   GenericKernelForContainsAndInside<Real_v, true>(unplaced, point, completelyinside, completelyoutside);
0365 
0366   inside = Inside_t(EInside::kSurface);
0367   vecCore::MaskedAssign(inside, (InsideBool_v)completelyoutside, Inside_t(EInside::kOutside));
0368   vecCore::MaskedAssign(inside, (InsideBool_v)completelyinside, Inside_t(EInside::kInside));
0369 }
0370 
0371 //______________________________________________________________________________
0372 template <typename Real_v, typename Inside_t>
0373 VECGEOM_FORCE_INLINE
0374 VECCORE_ATT_HOST_DEVICE
0375 void GenTrapImplementation::Inside(UnplacedStruct_t const &unplaced, Vector3D<Real_v> const &point, Inside_t &inside)
0376 {
0377   // Generic dispatcher to Inside
0378   GenTrapImplementation::InsideGeneric<Real_v, Inside_t>(unplaced, point, inside);
0379 }
0380 
0381 //______________________________________________________________________________
0382 template <>
0383 VECGEOM_FORCE_INLINE
0384 VECCORE_ATT_HOST_DEVICE
0385 void GenTrapImplementation::Inside<Precision, Inside_t>(UnplacedStruct_t const &unplaced,
0386                                                         Vector3D<Precision> const &point, Inside_t &inside)
0387 {
0388 // Scalar specialization for Inside function
0389 #ifndef VECCORE_CUDA
0390   if (unplaced.fTslHelper) {
0391     inside         = EInside::kOutside;
0392     Precision safZ = vecCore::math::Abs(point.z()) - unplaced.fDz;
0393     if (safZ > kTolerance) return;
0394     bool insideZ = safZ < -kTolerance;
0395     inside       = unplaced.fTslHelper->Inside(point);
0396     if (insideZ || inside == EInside::kOutside) return;
0397     inside = EInside::kSurface;
0398     return;
0399   }
0400 #endif
0401   GenTrapImplementation::InsideGeneric<Precision, Inside_t>(unplaced, point, inside);
0402 }
0403 
0404 //______________________________________________________________________________
0405 template <typename Real_v>
0406 VECGEOM_FORCE_INLINE
0407 VECCORE_ATT_HOST_DEVICE
0408 void GenTrapImplementation::DistanceToInGeneric(UnplacedStruct_t const &unplaced, Vector3D<Real_v> const &point,
0409                                                 Vector3D<Real_v> const &direction, Real_v const &stepMax,
0410                                                 Real_v &distance)
0411 {
0412 
0413   using Bool_v = vecCore::Mask_v<Real_v>;
0414 
0415   // do a quick boundary box check (Arb8, USolids) is doing this
0416   // unplaced.GetBBox();
0417   // actually this could also give us some indication which face is likely to be hit
0418 
0419   Real_v bbdistance = Real_v(kInfLength);
0420   BoxImplementation::DistanceToIn(BoxStruct<Precision>(unplaced.fBBdimensions), point - unplaced.fBBorigin, direction,
0421                                   stepMax, bbdistance);
0422 
0423   distance = InfinityLength<Real_v>();
0424 
0425   // do a check on bbdistance
0426   // if none of the tracks can hit even the bounding box; just return
0427   Bool_v done = bbdistance >= InfinityLength<Real_v>();
0428   if (vecCore::MaskFull(done)) return;
0429 
0430   // some particle could hit z
0431   Real_v zsafety = Abs(point.z()) - unplaced.fDz;
0432   Bool_v canhitz = zsafety > Real_v(MakeMinusTolerant<true>(0.));
0433   canhitz        = canhitz && (point.z() * direction.z() < 0); // coming towards the origin
0434   canhitz        = canhitz && (!done);
0435 
0436   if (!vecCore::MaskEmpty(canhitz)) {
0437     //   std::cerr << "can potentially hit\n";
0438     // calculate distance to z-plane ( see Box algorithm )
0439     // check if hit point is inside top or bottom polygon
0440     Real_v next = zsafety / NonZeroAbs(direction.z());
0441     // transport to z-height of planes
0442     Real_v coord1 = point.x() + next * direction.x();
0443     Real_v coord2 = point.y() + next * direction.y();
0444     Bool_v top    = direction.z() < 0;
0445     Bool_v hits   = IsInTopOrBottomPolygon<Real_v>(unplaced, coord1, coord2, top);
0446     hits          = hits && canhitz;
0447     vecCore::MaskedAssign(distance, hits, bbdistance);
0448     done = done || hits;
0449     if (vecCore::MaskFull(done)) return;
0450   }
0451 
0452   // now treat lateral surfaces
0453   Real_v disttoplanes = unplaced.fSurfaceShell.DistanceToIn<Real_v>(point, direction, done);
0454   vecCore__MaskedAssignFunc(distance, !done, Min(disttoplanes, distance));
0455 }
0456 
0457 //______________________________________________________________________________
0458 template <typename Real_v>
0459 VECGEOM_FORCE_INLINE
0460 VECCORE_ATT_HOST_DEVICE
0461 void GenTrapImplementation::DistanceToIn(UnplacedStruct_t const &unplaced, Vector3D<Real_v> const &point,
0462                                          Vector3D<Real_v> const &direction, Real_v const &stepMax, Real_v &distance)
0463 {
0464   GenTrapImplementation::DistanceToInGeneric<Real_v>(unplaced, point, direction, stepMax, distance);
0465 }
0466 
0467 //______________________________________________________________________________
0468 template <>
0469 VECGEOM_FORCE_INLINE
0470 VECCORE_ATT_HOST_DEVICE
0471 void GenTrapImplementation::DistanceToIn(UnplacedStruct_t const &unplaced, Vector3D<Precision> const &point,
0472                                          Vector3D<Precision> const &direction, Precision const &stepMax,
0473                                          Precision &distance)
0474 {
0475 // Scalar specialization of DistanceToIn
0476 #ifndef VECCORE_CUDA
0477   if (unplaced.fTslHelper) {
0478     Precision invdirz = 1. / NonZero(direction.z());
0479     distance          = unplaced.fTslHelper->DistanceToIn<false>(point, direction, invdirz, stepMax);
0480     return;
0481   }
0482 #endif
0483   GenTrapImplementation::DistanceToInGeneric<Precision>(unplaced, point, direction, stepMax, distance);
0484 }
0485 
0486 //______________________________________________________________________________
0487 template <typename Real_v>
0488 VECGEOM_FORCE_INLINE
0489 VECCORE_ATT_HOST_DEVICE
0490 void GenTrapImplementation::DistanceToOutGeneric(UnplacedStruct_t const &unplaced, Vector3D<Real_v> const &point,
0491                                                  Vector3D<Real_v> const &direction, Real_v const & /* stepMax */,
0492                                                  Real_v &distance)
0493 {
0494 
0495   using Bool_v = vecCore::Mask_v<Real_v>;
0496 
0497   // we should check here the compilation condition
0498   // that treatNormal=true can only happen when Backend=kScalar
0499   // TODO: do this with some nice template features
0500 
0501   Bool_v negDirMask = direction.z() < 0;
0502   Real_v sign(1.0);
0503   vecCore__MaskedAssignFunc(sign, negDirMask, Real_v(-1.));
0504   //    Real_v invDirZ = 1./NonZero(direction.z());
0505   // this construct costs one multiplication more
0506   Real_v distmin = (sign * unplaced.fDz - point.z()) / NonZero(direction.z());
0507 
0508   Real_v distplane = unplaced.fSurfaceShell.DistanceToOut<Real_v>(point, direction);
0509   distance         = Min(distmin, distplane);
0510 }
0511 
0512 template <typename Real_v>
0513 VECGEOM_FORCE_INLINE
0514 VECCORE_ATT_HOST_DEVICE
0515 void GenTrapImplementation::DistanceToOut(UnplacedStruct_t const &unplaced, Vector3D<Real_v> const &point,
0516                                           Vector3D<Real_v> const &direction, Real_v const &stepMax, Real_v &distance)
0517 {
0518   GenTrapImplementation::DistanceToOutGeneric<Real_v>(unplaced, point, direction, stepMax, distance);
0519 }
0520 
0521 //______________________________________________________________________________
0522 template <>
0523 VECGEOM_FORCE_INLINE
0524 VECCORE_ATT_HOST_DEVICE
0525 void GenTrapImplementation::DistanceToOut(UnplacedStruct_t const &unplaced, Vector3D<Precision> const &point,
0526                                           Vector3D<Precision> const &direction, Precision const &stepMax,
0527                                           Precision &distance)
0528 
0529 {
0530 #ifndef VECCORE_CUDA
0531   if (unplaced.fTslHelper) {
0532     distance = unplaced.fTslHelper->DistanceToOut(point, direction);
0533     return;
0534   }
0535 #endif
0536   GenTrapImplementation::DistanceToOutGeneric<Precision>(unplaced, point, direction, stepMax, distance);
0537 }
0538 
0539 //______________________________________________________________________________
0540 template <typename Real_v>
0541 VECGEOM_FORCE_INLINE
0542 VECCORE_ATT_HOST_DEVICE
0543 void GenTrapImplementation::SafetyToInGeneric(UnplacedStruct_t const &unplaced, Vector3D<Real_v> const &point,
0544                                               Real_v &safety)
0545 {
0546   // Generic implementation for SafetyToIn
0547   using Bool_v = vecCore::Mask_v<Real_v>;
0548 
0549   Bool_v inside;
0550   // Check if all points are outside bounding box
0551   BoxImplementation::Contains(BoxStruct<Precision>(unplaced.fBBdimensions), point - unplaced.fBBorigin, inside);
0552   if (vecCore::MaskEmpty(inside)) {
0553     // All points outside, so compute safety using the bounding box
0554     // This is not optimal if top and bottom faces are not on top of each other
0555     BoxImplementation::SafetyToIn(BoxStruct<Precision>(unplaced.fBBdimensions), point - unplaced.fBBorigin, safety);
0556     return;
0557   }
0558 
0559   // Do Z
0560   safety = Abs(point[2]) - unplaced.fDz;
0561   safety = unplaced.fSurfaceShell.SafetyToIn<Real_v>(point, safety);
0562 }
0563 
0564 //______________________________________________________________________________
0565 template <typename Real_v>
0566 VECGEOM_FORCE_INLINE
0567 VECCORE_ATT_HOST_DEVICE
0568 void GenTrapImplementation::SafetyToIn(UnplacedStruct_t const &unplaced, Vector3D<Real_v> const &point, Real_v &safety)
0569 {
0570   // Dispatcher for SafetyToIn
0571   GenTrapImplementation::SafetyToInGeneric<Real_v>(unplaced, point, safety);
0572 }
0573 
0574 //______________________________________________________________________________
0575 template <>
0576 VECGEOM_FORCE_INLINE
0577 VECCORE_ATT_HOST_DEVICE
0578 void GenTrapImplementation::SafetyToIn(UnplacedStruct_t const &unplaced, Vector3D<Precision> const &point,
0579                                        Precision &safety)
0580 {
0581 // Scalar specialization for SafetyToIn
0582 #ifndef VECCORE_CUDA
0583   if (unplaced.fTslHelper) {
0584     safety = unplaced.fTslHelper->SafetyToIn(point);
0585     return;
0586   }
0587 #endif
0588   GenTrapImplementation::SafetyToInGeneric<Precision>(unplaced, point, safety);
0589 }
0590 
0591 //______________________________________________________________________________
0592 template <typename Real_v>
0593 VECCORE_ATT_HOST_DEVICE
0594 void GenTrapImplementation::SafetyToOutGeneric(UnplacedStruct_t const &unplaced, Vector3D<Real_v> const &point,
0595                                                Real_v &safety)
0596 {
0597   // Generic implementation for SafetyToOut
0598   // Do Z
0599   safety = unplaced.fDz - Abs(point[2]);
0600   safety = unplaced.fSurfaceShell.SafetyToOut<Real_v>(point, safety);
0601 }
0602 
0603 //______________________________________________________________________________
0604 template <typename Real_v>
0605 VECCORE_ATT_HOST_DEVICE
0606 void GenTrapImplementation::SafetyToOut(UnplacedStruct_t const &unplaced, Vector3D<Real_v> const &point, Real_v &safety)
0607 {
0608   // Dispatcher for SafetyToOut
0609   GenTrapImplementation::SafetyToOutGeneric<Real_v>(unplaced, point, safety);
0610 }
0611 
0612 //______________________________________________________________________________
0613 template <>
0614 VECGEOM_FORCE_INLINE
0615 VECCORE_ATT_HOST_DEVICE
0616 void GenTrapImplementation::SafetyToOut(UnplacedStruct_t const &unplaced, Vector3D<Precision> const &point,
0617                                         Precision &safety)
0618 {
0619 // Scalar specialization for SafetyToOut
0620 #ifndef VECCORE_CUDA
0621   if (unplaced.fTslHelper) {
0622     safety = unplaced.fTslHelper->SafetyToOut(point);
0623     return;
0624   }
0625 #endif
0626   GenTrapImplementation::SafetyToOutGeneric<Precision>(unplaced, point, safety);
0627 }
0628 
0629 //______________________________________________________________________________
0630 template <typename Real_v, typename Bool_v>
0631 VECCORE_ATT_HOST_DEVICE
0632 void GenTrapImplementation::NormalKernel(UnplacedStruct_t const &unplaced, Vector3D<Real_v> const &point,
0633                                          Vector3D<Real_v> &normal, Bool_v &valid)
0634 {
0635 
0636   // Computes the normal on a surface and returns it as a unit vector
0637   //   In case a point is further than tolerance_normal from a surface, set validNormal=false
0638   //   Must return a valid vector. (even if the point is not on the surface.)
0639 
0640   using Index_v = vecCore::Index_v<Real_v>;
0641 
0642   valid = Bool_v(true);
0643   normal.Set(0., 0., 0.);
0644   // Do bottom and top faces
0645   Real_v safz = Abs(unplaced.fDz - Abs(point.z()));
0646   Bool_v onZ  = (safz < 10. * kTolerance);
0647   vecCore__MaskedAssignFunc(normal[2], onZ && (point.z() > Real_v(0.)), Real_v(1.));
0648   vecCore__MaskedAssignFunc(normal[2], onZ && (point.z() < Real_v(0.)), Real_v(-1.));
0649 
0650   //    if (vecCore::EarlyReturnAllowed()) {
0651   if (vecCore::MaskFull(onZ)) {
0652     return;
0653   }
0654   //    }
0655   //  Real_v done = onZ;
0656   // Get the closest edge (point should be on this edge within tolerance)
0657   Real_v cf = unplaced.fHalfInverseDz * (unplaced.fDz - point.z());
0658   Real_v vertexX[4];
0659   Real_v vertexY[4];
0660   for (int i = 0; i < 4; i++) {
0661     // calculate x-y positions of vertex i at this z-height
0662     vertexX[i] = unplaced.fVerticesX[i + 4] + cf * unplaced.fConnectingComponentsX[i];
0663     vertexY[i] = unplaced.fVerticesY[i + 4] + cf * unplaced.fConnectingComponentsY[i];
0664   }
0665   Real_v seg;
0666   Real_v frac;
0667   GetClosestEdge<Real_v>(point, vertexX, vertexY, seg, frac);
0668   vecCore__MaskedAssignFunc(frac, frac < Real_v(0.), Real_v(0.));
0669   Index_v iseg = (Index_v)seg;
0670   if (unplaced.IsPlanar()) {
0671     // Normals for the planar case are pre-computed
0672     Vertex_t const *normals = unplaced.fSurfaceShell.GetNormals();
0673     normal                  = normals[iseg];
0674     return;
0675   }
0676   Index_v jseg = (iseg + 1) % 4;
0677   Real_v x0    = vertexX[iseg];
0678   Real_v y0    = vertexY[iseg];
0679   Real_v x2    = vertexX[jseg];
0680   Real_v y2    = vertexY[jseg];
0681   x0 += frac * (x2 - x0);
0682   y0 += frac * (y2 - y0);
0683   Real_v x1 = unplaced.fVerticesX[iseg + 4];
0684   Real_v y1 = unplaced.fVerticesY[iseg + 4];
0685   x1 += frac * (unplaced.fVerticesX[jseg + 4] - x1);
0686   y1 += frac * (unplaced.fVerticesY[jseg + 4] - y1);
0687   Real_v ax = x1 - x0;
0688   Real_v ay = y1 - y0;
0689   Real_v az = unplaced.fDz - point.z();
0690   Real_v bx = x2 - x0;
0691   Real_v by = y2 - y0;
0692   Real_v bz = Real_v(0.);
0693   // Cross product of the vector given by the section segment (that contains the
0694   // point) at z=point[2] and the vector connecting the point projection to its
0695   // correspondent on the top edge.
0696   normal.Set(ay * bz - az * by, az * bx - ax * bz, ax * by - ay * bx);
0697   normal.Normalize();
0698 }
0699 
0700 //______________________________________________________________________________
0701 template <typename Real_v>
0702 VECCORE_ATT_HOST_DEVICE
0703 void GenTrapImplementation::GetClosestEdge(Vector3D<Real_v> const &point, Real_v vertexX[4], Real_v vertexY[4],
0704                                            Real_v &iseg, Real_v &fraction)
0705 {
0706   /// Get index of the edge of the quadrilater represented by vert closest to point.
0707   /// If [P1,P2] is the closest segment and P is the point, the function returns the fraction of the
0708   /// projection of (P1P) over (P1P2). If projection of P is not in range [P1,P2] return -1.
0709 
0710   using Bool_v = vecCore::Mask_v<Real_v>;
0711 
0712   iseg = Real_v(0.);
0713   //  Real_v p1X, p1Y, p2X, p2Y;
0714   Real_v lsq, dx, dy, dpx, dpy, u;
0715   fraction    = Real_v(-1.);
0716   Real_v safe = InfinityLength<Real_v>();
0717   Real_v ssq  = InfinityLength<Real_v>();
0718   for (int i = 0; i < 4; ++i) {
0719     int j = (i + 1) % 4;
0720     dx    = vertexX[j] - vertexX[i];
0721     dy    = vertexY[j] - vertexY[i];
0722     dpx   = point.x() - vertexX[i];
0723     dpy   = point.y() - vertexY[i];
0724     lsq   = dx * dx + dy * dy;
0725     // Current segment collapsed to a point
0726     Bool_v collapsed = lsq < kTolerance;
0727     if (!vecCore::MaskEmpty(collapsed)) {
0728       vecCore__MaskedAssignFunc(ssq, lsq < kTolerance, dpx * dpx + dpy * dpy);
0729       // Missing a masked assign allowing to perform multiple assignments...
0730       vecCore::MaskedAssign(iseg, ssq < safe, (Precision)i);
0731       vecCore__MaskedAssignFunc(fraction, ssq < safe, Real_v(-1.));
0732       vecCore::MaskedAssign(safe, ssq < safe, ssq);
0733       if (vecCore::MaskFull(collapsed)) continue;
0734     }
0735     // Projection fraction
0736     u = (dpx * dx + dpy * dy) / NonZero(lsq);
0737     vecCore__MaskedAssignFunc(dpx, u > 1 && !collapsed, point.x() - vertexX[j]);
0738     vecCore__MaskedAssignFunc(dpy, u > 1 && !collapsed, point.y() - vertexY[j]);
0739     vecCore__MaskedAssignFunc(dpx, u >= 0 && u <= 1 && !collapsed, dpx - u * dx);
0740     vecCore__MaskedAssignFunc(dpy, u >= 0 && u <= 1 && !collapsed, dpy - u * dy);
0741     vecCore__MaskedAssignFunc(u, (u > 1 || u < 0) && !collapsed, Real_v(-1.));
0742     ssq = dpx * dpx + dpy * dpy;
0743     vecCore::MaskedAssign(iseg, ssq < safe, (Precision)i);
0744     vecCore::MaskedAssign(fraction, ssq < safe, u);
0745     vecCore::MaskedAssign(safe, ssq < safe, ssq);
0746   }
0747 }
0748 
0749 //*****************************
0750 //**** Implementations end here
0751 //*****************************
0752 } // namespace VECGEOM_IMPL_NAMESPACE
0753 } // namespace vecgeom
0754 
0755 #endif /* GENTRAPIMPLEMENTATION_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!