volumes/kernel/CutTubeImplementation.h

0001 /*
0002  * CutTubeImplementation.h
0003  *
0004  *  Created on: 03.11.2016
0005  *      Author: mgheata
0006  */
0007
0008 #ifndef VECGEOM_VOLUMES_KERNEL_CUTTUBEIMPLEMENTATION_H_
0009 #define VECGEOM_VOLUMES_KERNEL_CUTTUBEIMPLEMENTATION_H_
0010
0011 #include <cstdio>
0012
0013 #include "VecGeom/base/Vector3D.h"
0014 #include "VecGeom/volumes/kernel/GenericKernels.h"
0015 #include "VecGeom/volumes/CutTubeStruct.h"
0016 #include "TubeImplementation.h"
0017
0018 namespace vecgeom {
0019
0020 VECGEOM_DEVICE_FORWARD_DECLARE(struct CutTubeImplementation;);
0021 VECGEOM_DEVICE_DECLARE_CONV(struct, CutTubeImplementation);
0022
0023 inline namespace VECGEOM_IMPL_NAMESPACE {
0024
0025 class PlacedCutTube;
0026 class UnplacedCutTube;
0027
0028 template <typename T>
0029 struct CutTubeStruct;
0030
0031 struct CutTubeImplementation {
0032
0033   using PlacedShape_t    = PlacedCutTube;
0034   using UnplacedStruct_t = CutTubeStruct<Precision>;
0035   using UnplacedVolume_t = UnplacedCutTube;
0036
0037   template <typename Real_v, typename Bool_v>
0038   VECGEOM_FORCE_INLINE VECCORE_ATT_HOST_DEVICE static void Contains(UnplacedStruct_t const &unplaced,
0039                                                                     Vector3D<Real_v> const &point, Bool_v &inside);
0040
0041   template <typename Real_v, typename Inside_t>
0042   VECGEOM_FORCE_INLINE VECCORE_ATT_HOST_DEVICE static void Inside(UnplacedStruct_t const &unplaced,
0043                                                                   Vector3D<Real_v> const &point, Inside_t &inside);
0044
0045   template <typename Real_v>
0046   VECGEOM_FORCE_INLINE VECCORE_ATT_HOST_DEVICE static void DistanceToIn(UnplacedStruct_t const &unplaced,
0047                                                                         Vector3D<Real_v> const &point,
0048                                                                         Vector3D<Real_v> const &direction,
0049                                                                         Real_v const &stepMax, Real_v &distance);
0050
0051   template <typename Real_v>
0052   VECGEOM_FORCE_INLINE VECCORE_ATT_HOST_DEVICE static void DistanceToInKernel(UnplacedStruct_t const &unplaced,
0053                                                                               Vector3D<Real_v> const &point,
0054                                                                               Vector3D<Real_v> const &direction,
0055                                                                               Real_v const &stepMax, Real_v &distance);
0056
0057   template <typename Real_v>
0058   VECGEOM_FORCE_INLINE VECCORE_ATT_HOST_DEVICE static void DistanceToOut(UnplacedStruct_t const &unplaced,
0059                                                                          Vector3D<Real_v> const &point,
0060                                                                          Vector3D<Real_v> const &direction,
0061                                                                          Real_v const &stepMax, Real_v &distance);
0062
0063   template <typename Real_v>
0064   VECGEOM_FORCE_INLINE VECCORE_ATT_HOST_DEVICE static void SafetyToIn(UnplacedStruct_t const &unplaced,
0065                                                                       Vector3D<Real_v> const &point, Real_v &safety);
0066
0067   template <typename Real_v>
0068   VECGEOM_FORCE_INLINE VECCORE_ATT_HOST_DEVICE static void SafetyToOut(UnplacedStruct_t const &unplaced,
0069                                                                        Vector3D<Real_v> const &point, Real_v &safety);
0070
0071   template <typename Real_v, typename Bool_v>
0072   VECGEOM_FORCE_INLINE VECCORE_ATT_HOST_DEVICE static void NormalKernel(UnplacedStruct_t const &unplaced,
0073                                                                         Vector3D<Real_v> const &point,
0074                                                                         Vector3D<Real_v> &normal, Bool_v &valid);
0075 }; // End struct CutTubeImplementation
0076
0077 //********************************
0078 //**** implementations start here
0079 //********************************/
0080
0081 //______________________________________________________________________________
0082 template <typename Real_v, typename Bool_v>
0083 VECGEOM_FORCE_INLINE VECCORE_ATT_HOST_DEVICE void CutTubeImplementation::Contains(UnplacedStruct_t const &unplaced,
0084                                                                                   Vector3D<Real_v> const &point,
0085                                                                                   Bool_v &contains)
0086 {
0087   contains                = Bool_v(false);
0088   Bool_v inside_cutplanes = Bool_v(false);
0089   // Check the cut planes first
0090   unplaced.GetCutPlanes().Contains<Real_v, Bool_v>(point, inside_cutplanes);
0091   if (vecCore::EarlyReturnAllowed()) {
0092     if (vecCore::MaskEmpty(inside_cutplanes)) return;
0093   }
0094   // Check the tube
0095   TubeImplementation<TubeTypes::UniversalTube>::Contains<Real_v>(unplaced.GetTubeStruct(), point, contains);
0096   contains &= inside_cutplanes;
0097 }
0098
0099 //______________________________________________________________________________
0100 template <typename Real_v, typename Inside_t>
0101 VECGEOM_FORCE_INLINE VECCORE_ATT_HOST_DEVICE void CutTubeImplementation::Inside(UnplacedStruct_t const &unplaced,
0102                                                                                 Vector3D<Real_v> const &point,
0103                                                                                 Inside_t &inside)
0104 {
0105   inside                    = Inside_t(EInside::kOutside);
0106   Inside_t inside_cutplanes = Inside_t(EInside::kOutside);
0107   // Check the cut planes first
0108   unplaced.GetCutPlanes().Inside<Real_v, Inside_t>(point, inside_cutplanes);
0109   if (vecCore::EarlyReturnAllowed()) {
0110     if (vecCore::MaskFull(inside_cutplanes == inside)) return;
0111   }
0112   // Check the tube
0113   TubeImplementation<TubeTypes::UniversalTube>::Inside<Real_v, Inside_t>(unplaced.GetTubeStruct(), point, inside);
0114   vecCore::MaskedAssign(inside,
0115                         (inside_cutplanes == EInside::kOutside) ||
0116                             (inside_cutplanes == EInside::kSurface && inside != EInside::kOutside),
0117                         inside_cutplanes);
0118 }
0119
0120 template <typename Real_v>
0121 VECGEOM_FORCE_INLINE VECCORE_ATT_HOST_DEVICE void CutTubeImplementation::DistanceToIn(UnplacedStruct_t const &unplaced,
0122                                                                                       Vector3D<Real_v> const &pointT,
0123                                                                                       Vector3D<Real_v> const &dir,
0124                                                                                       Real_v const &stepMax,
0125                                                                                       Real_v &distance)
0126 {
0127
0128   Vector3D<Real_v> point = pointT;
0129   Real_v ptDist          = point.Mag();
0130   Real_v distToMove(0.);
0131   using Bool_v    = vecCore::Mask_v<Real_v>;
0132   Precision order = 100.;
0133   // Bool_v cond = (ptDist > order*unplaced.fTubeStruct.fMaxVal);
0134   Bool_v cond = (ptDist > order * unplaced.fMaxVal);
0135   /* if the point is at a distance (DIST) of more than 100 times of the maximum dimension
0136    * (of the shape) from the origin of shape, then before calculating distance, first
0137    * manually move the point with distance ( distToMove = DIST-100.*maxDim) along the
0138    * direction, and then calculate DistanceToIn of new moved point using DistanceToInKernel,
0139    *
0140    * The final distance will be (distToMove + DistanceToIn),
0141    */
0142   // vecCore__MaskedAssignFunc(distToMove,cond,(ptDist-Real_v(order*unplaced.fTubeStruct.fMaxVal)));
0143   vecCore__MaskedAssignFunc(distToMove, cond, (ptDist - Real_v(order * unplaced.fMaxVal)));
0144   vecCore__MaskedAssignFunc(point, cond, point + distToMove * dir);
0145   DistanceToInKernel<Real_v>(unplaced, point, dir, stepMax, distance);
0146   distance += distToMove;
0147 }
0148 //______________________________________________________________________________
0149 template <typename Real_v>
0150 VECGEOM_FORCE_INLINE VECCORE_ATT_HOST_DEVICE void CutTubeImplementation::DistanceToInKernel(
0151     UnplacedStruct_t const &unplaced, Vector3D<Real_v> const &point, Vector3D<Real_v> const &direction,
0152     Real_v const &stepMax, Real_v &distance)
0153 {
0154
0155 #define USE_CONV_WRONG_SIDE 1
0156 #define USE_CONV_FROM_BOUNDARY 1
0157   // Compute distance to cut planes
0158   using Bool_v = vecCore::Mask_v<Real_v>;
0159   distance     = InfinityLength<Real_v>();
0160
0161 #if USE_CONV_WRONG_SIDE == 1
0162   using Inside_v = vecCore::Index_v<Real_v>;
0163   // Check the cut planes first
0164   Inside_v inside_cutplanes;
0165   unplaced.GetCutPlanes().Inside<Real_v, Inside_v>(point, inside_cutplanes);
0166   Inside_v instart;
0167   // Check the tube
0168   TubeImplementation<TubeTypes::UniversalTube>::Inside<Real_v, Inside_v>(unplaced.GetTubeStruct(), point, instart);
0169   vecCore::MaskedAssign(instart,
0170                         (inside_cutplanes == EInside::kOutside) ||
0171                             (inside_cutplanes == EInside::kSurface && instart != EInside::kOutside),
0172                         inside_cutplanes);
0173   // Points already inside have to return negative distance
0174   vecCore__MaskedAssignFunc(distance, instart == EInside::kInside, Real_v(-1.));
0175 #endif
0176
0177   Bool_v inside;
0178   Real_v dplanes, safplanes;
0179   // Compute distance to cut planes
0180   unplaced.GetCutPlanes().DistanceToIn<Real_v>(point, direction, dplanes);
0181   // We need the safety to cut out points already inside the cut planes.
0182   // Unfortunately there is no easy way to get the right normal in vector mode to avoid safety calculation
0183   unplaced.GetCutPlanes().SafetyToIn<Real_v>(point, safplanes);
0184
0185   // Mark tracks hitting the planes
0186   Bool_v hitplanes = vecCore::math::Abs(dplanes) < stepMax && safplanes > Real_v(-kTolerance);
0187
0188 #if USE_CONV_WRONG_SIDE == 1
0189   if (vecCore::EarlyReturnAllowed()) {
0190     if (vecCore::MaskFull((inside_cutplanes != EInside::kInside) && !hitplanes)) // No particles are hitting
0191       return;
0192   }
0193 #endif
0194
0195   // Propagate with dplanes only the particles that are hitting
0196   Vector3D<Real_v> propagated = point;
0197   vecCore__MaskedAssignFunc(dplanes, !hitplanes, Real_v(0.));
0198   if (vecCore::EarlyReturnAllowed()) {
0199     if (!vecCore::MaskEmpty(hitplanes)) {
0200       propagated += dplanes * direction;
0201       // Hitting the planes does not guarantee that the hit point is between them
0202       unplaced.GetCutPlanes().Inside<Real_v, Inside_v>(propagated, inside_cutplanes);
0203     }
0204   } else {
0205     // In CUDA we have to propagate in all cases
0206     propagated += dplanes * direction;
0207     // Hitting the planes does not guarantee that the hit point is between them
0208     unplaced.GetCutPlanes().Inside<Real_v, Inside_v>(propagated, inside_cutplanes);
0209   }
0210   Bool_v done = inside_cutplanes == EInside::kOutside;
0211 #if USE_CONV_WRONG_SIDE == 1
0212   done |= instart == EInside::kInside;
0213 #endif
0214   // Tracks that cannot get between the cut planes cannot hit
0215   if (vecCore::EarlyReturnAllowed()) {
0216     if (vecCore::MaskFull(done)) // No particles are hitting
0217       return;
0218   }
0219
0220   // All propagated points not yet marked as done should now lie between the cut planes
0221   // Check now which of the propagated points already entering the solid
0222   TubeImplementation<TubeTypes::UniversalTube>::Inside<Real_v>(unplaced.GetTubeStruct(), propagated, instart);
0223   inside = hitplanes && (instart != EInside::kOutside);
0224   vecCore::MaskedAssign(distance, inside && !done, dplanes);
0225   done |= inside;
0226
0227   if (vecCore::EarlyReturnAllowed()) {
0228     if (vecCore::MaskFull(done)) { // Some particles hit top/bottom
0229       // The line below is needed for the convention
0230       vecCore__MaskedAssignFunc(distance, vecCore::math::Abs(distance) < Real_v(kTolerance), Real_v(0.));
0231       return;
0232     }
0233   }
0234
0235   // The limit distance for tube crossing cannot exceed the distance to
0236   // exiting the cut planes
0237   Real_v dexit = InfinityLength<Real_v>();
0238   unplaced.GetCutPlanes().DistanceToOut<Real_v>(propagated, direction, dexit);
0239
0240   // Compute distance to tube
0241   Real_v dtube = InfinityLength<Real_v>();
0242   TubeImplementation<TubeTypes::UniversalTube>::DistanceToInKernel<Real_v>(unplaced.GetTubeStruct(), propagated,
0243                                                                            direction, stepMax, dtube);
0244   // A.G Propagation to cut planes can put the point inside the tube, so DistanceToIn may return -1
0245   // In such case we need to set dtube to 0, otherwise we may get wrong negative answers
0246   vecCore__MaskedAssignFunc(dtube, dtube < Real_v(0.), Real_v(0.));
0247   vecCore__MaskedAssignFunc(dtube, dexit < dtube, InfinityLength<Real_v>());
0248   vecCore__MaskedAssignFunc(distance, !done && (dtube + dplanes) < stepMax, dtube + dplanes);
0249 // The line below is needed for the convention
0250 #if USE_CONV_FROM_BOUNDARY == 1
0251 //  vecCore__MaskedAssignFunc(distance, vecCore::math::Abs(distance) < Real_v(kTolerance), Real_v(0.));
0252 #endif
0253 }
0254
0255 //______________________________________________________________________________
0256 template <typename Real_v>
0257 VECGEOM_FORCE_INLINE VECCORE_ATT_HOST_DEVICE void CutTubeImplementation::DistanceToOut(
0258     UnplacedStruct_t const &unplaced, Vector3D<Real_v> const &point, Vector3D<Real_v> const &direction,
0259     Real_v const &stepMax, Real_v &distance)
0260 {
0261   // Compute distance to cut planes
0262   distance = InfinityLength<Real_v>();
0263   unplaced.GetCutPlanes().DistanceToOut<Real_v>(point, direction, distance);
0264
0265   // Compute distance to tube
0266   Real_v dtube = InfinityLength<Real_v>();
0267   TubeImplementation<TubeTypes::UniversalTube>::DistanceToOut<Real_v>(unplaced.GetTubeStruct(), point, direction,
0268                                                                       stepMax, dtube);
0269   vecCore::MaskedAssign(distance, dtube < distance, dtube);
0270   // The line below is needed to avoid din=dout=0 when starting from a boundary
0271   // A.G. There should be a dir.dot.norm check to avoid the condition above. Rounding to tolerance
0272   //      prevents G4UnionSolid to exit correctly from the boundary crossing loop
0273   // vecCore__MaskedAssignFunc(distance, distance >= Real_v(0.) && distance < Real_v(kTolerance), Real_v(kTolerance));
0274 }
0275
0276 //______________________________________________________________________________
0277 template <typename Real_v>
0278 VECGEOM_FORCE_INLINE VECCORE_ATT_HOST_DEVICE void CutTubeImplementation::SafetyToIn(UnplacedStruct_t const &unplaced,
0279                                                                                     Vector3D<Real_v> const &point,
0280                                                                                     Real_v &safety)
0281 {
0282   // Compute safety to cut planes. These will contain the sign, i.e. if on
0283   // the wrong side they will be negative
0284   unplaced.GetCutPlanes().SafetyToIn<Real_v>(point, safety);
0285   Real_v saftube;
0286   // Compute safety to tube
0287   TubeImplementation<TubeTypes::UniversalTube>::SafetyToIn<Real_v>(unplaced.GetTubeStruct(), point, saftube);
0288   // The safety is the maximum of the 2 values
0289   vecCore::MaskedAssign(safety, saftube > safety, saftube);
0290   // The line below is needed for the rounding convention
0291   //  vecCore__MaskedAssignFunc(safety, vecCore::math::Abs(safety) < Real_v(kTolerance), Real_v(0.));
0292 }
0293
0294 //______________________________________________________________________________
0295 template <typename Real_v>
0296 VECCORE_ATT_HOST_DEVICE void CutTubeImplementation::SafetyToOut(UnplacedStruct_t const &unplaced,
0297                                                                 Vector3D<Real_v> const &point, Real_v &safety)
0298 {
0299   // Compute safety to cut planes. These will contain the sign, i.e. if on
0300   // the wrong side they will be negative
0301   unplaced.GetCutPlanes().SafetyToOut<Real_v>(point, safety);
0302   Real_v saftube;
0303   // Compute safety to tube
0304   TubeImplementation<TubeTypes::UniversalTube>::SafetyToOut<Real_v>(unplaced.GetTubeStruct(), point, saftube);
0305   // The safety is the minimum of the 2 values
0306   vecCore::MaskedAssign(safety, saftube < safety, saftube);
0307   // The line below is needed for the rounding convention
0308   //  vecCore__MaskedAssignFunc(safety, vecCore::math::Abs(safety) < Real_v(kTolerance), Real_v(0.));
0309 }
0310
0311 //______________________________________________________________________________
0312 template <typename Real_v, typename Bool_v>
0313 VECCORE_ATT_HOST_DEVICE void CutTubeImplementation::NormalKernel(UnplacedStruct_t const &unplaced,
0314                                                                  Vector3D<Real_v> const &point,
0315                                                                  Vector3D<Real_v> &normal, Bool_v &valid)
0316 {
0317   // Compute safety to cut planes
0318   valid = Bool_v(true);
0319   Real_v safcut;
0320   unplaced.GetCutPlanes().SafetyToOut<Real_v>(point, safcut);
0321   // Compute safety to tube
0322   Real_v saftube;
0323   TubeImplementation<TubeTypes::UniversalTube>::SafetyToOut<Real_v>(unplaced.GetTubeStruct(), point, saftube);
0324   Bool_v istube = vecCore::math::Abs(saftube) < vecCore::math::Abs(safcut);
0325   // The statement below works only as long as Real_v is scalar
0326   if (istube) {
0327     TubeImplementation<TubeTypes::UniversalTube>::NormalKernel<Real_v, Bool_v>(unplaced.GetTubeStruct(), point, normal,
0328                                                                                valid);
0329     return;
0330   }
0331   // Select the correct cut plane
0332   if (point.z() < 0)
0333     normal = unplaced.fCutPlanes.GetNormal(0);
0334   else
0335     normal = unplaced.fCutPlanes.GetNormal(1);
0336 }
0337
0338 //*****************************
0339 //**** Implementations end here
0340 //*****************************
0341 } // namespace VECGEOM_IMPL_NAMESPACE
0342 } // namespace vecgeom
0343
0344 #endif /* VECGEOM_VOLUMES_KERNEL_CUTTUBEIMPLEMENTATION_H_ */