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