VecGeom/volumes/Wedge.h

0001 /*
0002  * Wedge.h
0003  *
0004  *  Created on: 09.10.2014
0005  *      Author: swenzel
0006  */
0007
0008 #ifndef VECGEOM_VOLUMES_WEDGE_H_
0009 #define VECGEOM_VOLUMES_WEDGE_H_
0010
0011 #include "VecGeom/base/Global.h"
0012 #include "VecGeom/volumes/kernel/GenericKernels.h"
0013
0014 namespace vecgeom {
0015 inline namespace VECGEOM_IMPL_NAMESPACE {
0016
0017 /**
0018  * A class representing a wedge which is represented by an angle. It
0019  * can be used to divide 3D spaces or to clip wedges from solids.
0020  * The wedge has an "inner" and "outer" side. For an angle = 180 degree, the wedge is essentially
0021  * an ordinary halfspace. Usually the wedge is used to cut out "phi" sections along z-direction.
0022  *
0023  * Idea: should have Unplaced and PlacedWegdes, should have specializations
0024  * for "PhiWegde" and which are used in symmetric
0025  * shapes such as tubes or spheres.
0026  *
0027  * Note: This class is meant as an auxiliary class so it is a bit outside the ordinary volume
0028  * hierarchy.
0029  *
0030  *       / +++++++++++
0031  *      / ++++++++++++
0032  *     / +++++++++++++
0033  *    / +++++ INSIDE +
0034  *   / +++++++++++++++
0035  *  / fDPhi +++++++++
0036  * x------------------ ( this is at angle fSPhi )
0037  *
0038  *     OUTSIDE
0039  *
0040  */
0041 class Wedge {
0042
0043 private:
0044   Precision fSPhi{0.};               // starting angle
0045   Precision fDPhi{0.};               // delta angle representing/defining the wedge
0046   Vector3D<Precision> fAlongVector1; // vector along the first plane
0047   Vector3D<Precision> fAlongVector2; // vector aling the second plane
0048
0049   Vector3D<Precision> fNormalVector1; // normal vector for first plane
0050   // convention is that it points inwards
0051
0052   Vector3D<Precision> fNormalVector2; // normal vector for second plane
0053                                       // convention is that it points inwards
0054
0055 public:
0056   Wedge() = default;
0057
0058   VECCORE_ATT_HOST_DEVICE
0059   Wedge(Precision angle, Precision zeroangle = 0);
0060
0061   VECCORE_ATT_HOST_DEVICE
0062   ~Wedge() {}
0063
0064   VECCORE_ATT_HOST_DEVICE
0065   void SetStartPhi(Precision const &arg) { fSPhi = arg; }
0066
0067   VECCORE_ATT_HOST_DEVICE
0068   void SetDeltaPhi(Precision const &arg) { fDPhi = arg; }
0069
0070   VECCORE_ATT_HOST_DEVICE
0071   void Set(Precision const &dphi, Precision const &sphi)
0072   {
0073     SetStartPhi(sphi);
0074     SetDeltaPhi(dphi);
0075   }
0076
0077   VECCORE_ATT_HOST_DEVICE
0078   Vector3D<Precision> GetAlong1() const { return fAlongVector1; }
0079
0080   VECCORE_ATT_HOST_DEVICE
0081   Vector3D<Precision> GetAlong2() const { return fAlongVector2; }
0082
0083   VECCORE_ATT_HOST_DEVICE
0084   Vector3D<Precision> GetNormal1() const { return fNormalVector1; }
0085
0086   VECCORE_ATT_HOST_DEVICE
0087   Vector3D<Precision> GetNormal2() const { return fNormalVector2; }
0088
0089   /* Function Name : GetNormal<ForStartPhi>()
0090    *
0091    * The function is the templatized version GetNormal1() and GetNormal2() function and will
0092    * return the normal depending upon the boolean template parameter "ForStartPhi"
0093    * which if passed as true, will return normal to the StartingPhi of Wedge,
0094    * if passed as false, will return normal to the EndingPhi of Wedge
0095    *
0096    * from user point of view the same work can be done by calling GetNormal1() and GetNormal2()
0097    * functions, but this implementation will be used by "IsPointOnSurfaceAndMovingOut()" function
0098    */
0099   template <bool ForStartPhi>
0100   VECGEOM_FORCE_INLINE VECCORE_ATT_HOST_DEVICE Vector3D<Precision> GetNormal() const;
0101
0102   // very important:
0103   template <typename Backend>
0104   VECCORE_ATT_HOST_DEVICE typename Backend::bool_v Contains(Vector3D<typename Backend::precision_v> const &point) const;
0105
0106   // GL note: for tubes, use of TubeImpl::PointInCyclicalSector outperformed next two methods in vector mode
0107   template <typename Backend>
0108   VECCORE_ATT_HOST_DEVICE typename Backend::bool_v ContainsWithBoundary(
0109       Vector3D<typename Backend::precision_v> const &point) const;
0110
0111   template <typename Backend>
0112   VECCORE_ATT_HOST_DEVICE typename Backend::bool_v ContainsWithoutBoundary(
0113       Vector3D<typename Backend::precision_v> const &point) const;
0114
0115   template <typename Backend>
0116   VECCORE_ATT_HOST_DEVICE typename Backend::inside_v Inside(Vector3D<typename Backend::precision_v> const &point) const;
0117
0118   // static function determining if input points are on a plane surface which is part of a wedge
0119   // ( given by along and normal )
0120   template <typename Backend>
0121   VECCORE_ATT_HOST_DEVICE static typename Backend::bool_v IsOnSurfaceGeneric(
0122       Vector3D<Precision> const &alongVector, Vector3D<Precision> const &normalVector,
0123       Vector3D<typename Backend::precision_v> const &point);
0124
0125   /* Function Name :  IsOnSurfaceGeneric<Backend, ForStartPhi>()
0126    *
0127    * This version of IsOnSurfaceGeneric is having one more template parameter of type boolean,
0128    * which if passed as true, will check whether the point is on StartingPhi Surface of Wedge,
0129    * and if passed as false, will check whether the point is on EndingPhi Surface of Wedge
0130    *
0131    * this implementation will be used by "IsPointOnSurfaceAndMovingOut()" function.
0132    */
0133   template <typename Backend, bool ForStartPhi>
0134   VECGEOM_FORCE_INLINE VECCORE_ATT_HOST_DEVICE typename Backend::bool_v IsOnSurfaceGeneric(
0135       Vector3D<typename Backend::precision_v> const &point) const;
0136
0137   /* Function Name : IsPointOnSurfaceAndMovingOut<Backend, ForStartPhi, MovingOut>
0138    *
0139    * This function is written to check if the point is on surface and is moving inside or outside.
0140    * This will basically be used by "DistanceToInKernel()" and "DistanceToOutKernel()" of the shapes,
0141    * which uses wedge.
0142    *
0143    * It contains two extra template boolean parameters "ForStartPhi" and "MovingOut",
0144    * So call like "IsPointOnSurfaceAndMovingOut<Backend,true,true>" will check whether the points is on
0145    * the StartingPhi Surface of wedge and moving outside.
0146    *
0147    * So overall can be called in following four combinations
0148    * 1) "IsPointOnSurfaceAndMovingOut<Backend,true,true>" : Point on StartingPhi surface of wedge and moving OUT
0149    * 2) "IsPointOnSurfaceAndMovingOut<Backend,true,false>" : Point on StartingPhi surface of wedge and moving IN
0150    * 3) "IsPointOnSurfaceAndMovingOut<Backend,false,true>" : Point on EndingPhi surface of wedge and moving OUT
0151    * 2) "IsPointOnSurfaceAndMovingOut<Backend,false,false>" : Point on EndingPhi surface of wedge and moving IN
0152    *
0153    * Very useful for DistanceToIn and DistanceToOut.
0154    */
0155   template <typename Backend, bool ForStartPhi, bool MovingOut>
0156   VECGEOM_FORCE_INLINE VECCORE_ATT_HOST_DEVICE typename Backend::bool_v IsPointOnSurfaceAndMovingOut(
0157       Vector3D<typename Backend::precision_v> const &point, Vector3D<typename Backend::precision_v> const &dir) const;
0158
0159   VECCORE_ATT_HOST_DEVICE
0160   bool IsOnSurface1(Vector3D<Precision> const &point) const
0161   {
0162     return Wedge::IsOnSurfaceGeneric<kScalar>(fAlongVector1, fNormalVector1, point);
0163   }
0164
0165   VECCORE_ATT_HOST_DEVICE
0166   bool IsOnSurface2(Vector3D<Precision> const &point) const
0167   {
0168     return Wedge::IsOnSurfaceGeneric<kScalar>(fAlongVector2, fNormalVector2, point);
0169   }
0170
0171   /**
0172    * estimate of the smallest distance to the Wedge boundary when
0173    * the point is located outside the Wedge
0174    */
0175   template <typename Backend>
0176   VECCORE_ATT_HOST_DEVICE typename Backend::precision_v SafetyToIn(
0177       Vector3D<typename Backend::precision_v> const &point) const;
0178
0179   /**
0180    * estimate of the smallest distance to the Wedge boundary when
0181    * the point is located inside the Wedge ( within the defining phi angle )
0182    */
0183   template <typename Backend>
0184   VECCORE_ATT_HOST_DEVICE typename Backend::precision_v SafetyToOut(
0185       Vector3D<typename Backend::precision_v> const &point) const;
0186
0187   /**
0188    * estimate of the distance to the Wedge boundary with given direction
0189    */
0190   template <typename Backend>
0191   VECCORE_ATT_HOST_DEVICE void DistanceToIn(Vector3D<typename Backend::precision_v> const &point,
0192                                             Vector3D<typename Backend::precision_v> const &dir,
0193                                             typename Backend::precision_v &distWedge1,
0194                                             typename Backend::precision_v &distWedge2) const;
0195
0196   template <typename Backend>
0197   VECCORE_ATT_HOST_DEVICE void DistanceToOut(Vector3D<typename Backend::precision_v> const &point,
0198                                              Vector3D<typename Backend::precision_v> const &dir,
0199                                              typename Backend::precision_v &distWedge1,
0200                                              typename Backend::precision_v &distWedge2) const;
0201
0202   // this could be useful to be public such that other shapes can directly
0203   // use completelyinside + completelyoutside
0204
0205   template <typename Backend, bool ForInside>
0206   VECCORE_ATT_HOST_DEVICE void GenericKernelForContainsAndInside(
0207       Vector3D<typename Backend::precision_v> const &localPoint, typename Backend::bool_v &completelyinside,
0208       typename Backend::bool_v &completelyoutside) const;
0209
0210 }; // end of class Wedge
0211
0212 template <bool ForStartPhi>
0213 VECGEOM_FORCE_INLINE VECCORE_ATT_HOST_DEVICE Vector3D<Precision> Wedge::GetNormal() const
0214 {
0215   if (ForStartPhi)
0216     return fNormalVector1;
0217   else
0218     return fNormalVector2;
0219 }
0220
0221 template <typename Backend, bool ForStartPhi, bool MovingOut>
0222 VECGEOM_FORCE_INLINE VECCORE_ATT_HOST_DEVICE typename Backend::bool_v Wedge::IsPointOnSurfaceAndMovingOut(
0223     Vector3D<typename Backend::precision_v> const &point, Vector3D<typename Backend::precision_v> const &dir) const
0224 {
0225
0226   if (MovingOut)
0227     return IsOnSurfaceGeneric<Backend, ForStartPhi>(point) &&
0228            (dir.Dot(-GetNormal<ForStartPhi>()) > 0.005 * kHalfTolerance);
0229   else
0230     return IsOnSurfaceGeneric<Backend, ForStartPhi>(point) &&
0231            (dir.Dot(-GetNormal<ForStartPhi>()) < 0.005 * kHalfTolerance);
0232 }
0233
0234 template <typename Backend, bool ForStartPhi>
0235 VECGEOM_FORCE_INLINE VECCORE_ATT_HOST_DEVICE typename Backend::bool_v Wedge::IsOnSurfaceGeneric(
0236     Vector3D<typename Backend::precision_v> const &point) const
0237 {
0238
0239   if (ForStartPhi)
0240     return IsOnSurfaceGeneric<Backend>(fAlongVector1, fNormalVector1, point);
0241   else
0242     return IsOnSurfaceGeneric<Backend>(fAlongVector2, fNormalVector2, point);
0243 }
0244
0245 template <typename Backend>
0246 VECCORE_ATT_HOST_DEVICE typename Backend::inside_v Wedge::Inside(
0247     Vector3D<typename Backend::precision_v> const &point) const
0248 {
0249
0250   typedef typename Backend::bool_v Bool_t;
0251   Bool_t completelyinside, completelyoutside;
0252   GenericKernelForContainsAndInside<Backend, true>(point, completelyinside, completelyoutside);
0253   typename Backend::inside_v inside = EInside::kSurface;
0254   vecCore::MaskedAssign(inside, completelyoutside, EInside::kOutside);
0255   vecCore::MaskedAssign(inside, completelyinside, EInside::kInside);
0256   return inside;
0257 }
0258
0259 template <typename Backend>
0260 VECCORE_ATT_HOST_DEVICE typename Backend::bool_v Wedge::ContainsWithBoundary(
0261     Vector3D<typename Backend::precision_v> const &point) const
0262 {
0263   typedef typename Backend::bool_v Bool_t;
0264   Bool_t completelyinside, completelyoutside;
0265   GenericKernelForContainsAndInside<Backend, true>(point, completelyinside, completelyoutside);
0266   return !completelyoutside;
0267 }
0268
0269 template <typename Backend>
0270 VECCORE_ATT_HOST_DEVICE typename Backend::bool_v Wedge::ContainsWithoutBoundary(
0271     Vector3D<typename Backend::precision_v> const &point) const
0272 {
0273   typedef typename Backend::bool_v Bool_t;
0274   Bool_t completelyinside, completelyoutside;
0275   GenericKernelForContainsAndInside<Backend, true>(point, completelyinside, completelyoutside);
0276   return completelyinside;
0277 }
0278
0279 template <typename Backend>
0280 VECCORE_ATT_HOST_DEVICE typename Backend::bool_v Wedge::Contains(
0281     Vector3D<typename Backend::precision_v> const &point) const
0282 {
0283   typedef typename Backend::bool_v Bool_t;
0284   Bool_t unused;
0285   Bool_t outside;
0286   GenericKernelForContainsAndInside<Backend, false>(point, unused, outside);
0287   return !outside;
0288 }
0289
0290 // Implementation follows
0291 template <typename Backend, bool ForInside>
0292 VECCORE_ATT_HOST_DEVICE void Wedge::GenericKernelForContainsAndInside(
0293     Vector3D<typename Backend::precision_v> const &localPoint, typename Backend::bool_v &completelyinside,
0294     typename Backend::bool_v &completelyoutside) const
0295 {
0296   typedef typename Backend::precision_v Real_v;
0297
0298   // this part of the code assumes some symmetry knowledge and is currently only
0299   // correct for a PhiWedge assumed to be aligned along the z-axis.
0300   Real_v x      = localPoint.x();
0301   Real_v y      = localPoint.y();
0302   Real_v startx = fAlongVector1.x();
0303   Real_v starty = fAlongVector1.y();
0304   Real_v endx   = fAlongVector2.x();
0305   Real_v endy   = fAlongVector2.y();
0306
0307   Real_v startCheck = (-x * starty + y * startx);
0308   Real_v endCheck   = (-endx * y + endy * x);
0309
0310   completelyoutside = startCheck < 0.;
0311   if (fDPhi < kPi)
0312     completelyoutside |= endCheck < 0.;
0313   else
0314     completelyoutside &= endCheck < 0.;
0315   if (ForInside) {
0316     // TODO: see if the compiler optimizes across these function calls since
0317     // a couple of multiplications inside IsOnSurfaceGeneric are already done previously
0318     typename Backend::bool_v onSurface =
0319         Wedge::IsOnSurfaceGeneric<Backend>(fAlongVector1, fNormalVector1, localPoint) ||
0320         Wedge::IsOnSurfaceGeneric<Backend>(fAlongVector2, fNormalVector2, localPoint);
0321     completelyoutside &= !onSurface;
0322     completelyinside = !onSurface && !completelyoutside;
0323   }
0324 }
0325
0326 template <typename Backend>
0327 VECCORE_ATT_HOST_DEVICE typename Backend::bool_v Wedge::IsOnSurfaceGeneric(
0328     Vector3D<Precision> const &alongVector, Vector3D<Precision> const &normalVector,
0329     Vector3D<typename Backend::precision_v> const &point)
0330 {
0331   // on right side of half plane ??
0332   typedef typename Backend::bool_v Bool_v;
0333   Bool_v condition1 = alongVector.x() * point.x() + alongVector.y() * point.y() >= 0.;
0334   if (vecCore::MaskEmpty(condition1)) return Bool_v(false);
0335   // within the right distance to the plane ??
0336   Bool_v condition2 = Abs(normalVector.x() * point.x() + normalVector.y() * point.y()) < kTolerance;
0337   return condition1 && condition2;
0338 }
0339
0340 template <typename Backend>
0341 VECCORE_ATT_HOST_DEVICE typename Backend::precision_v Wedge::SafetyToOut(
0342     Vector3D<typename Backend::precision_v> const &point) const
0343 {
0344   typedef typename Backend::precision_v Float_t;
0345   // algorithm: calculate projections to both planes
0346   // return minimum / maximum depending on fAngle < PI or not
0347
0348   // assuming that we have z wedge and the planes pass through the origin
0349   Float_t dist1 = point.x() * fNormalVector1.x() + point.y() * fNormalVector1.y();
0350   Float_t dist2 = point.x() * fNormalVector2.x() + point.y() * fNormalVector2.y();
0351
0352   // std::cerr << "d1 " << dist1<<"  "<<point << "\n";
0353   // std::cerr << "d2 " << dist2<<"  "<<point << "\n";
0354
0355   if (fDPhi < kPi) {
0356     return Min(dist1, dist2);
0357   } else {
0358     // Float_t disttocorner = Sqrt(point.x()*point.x() + point.y()*point.y());
0359     // return Max(dist1,Max(dist2,disttocorner));
0360     return Max(dist1, dist2);
0361   }
0362 }
0363
0364 template <typename Backend>
0365 VECCORE_ATT_HOST_DEVICE typename Backend::precision_v Wedge::SafetyToIn(
0366     Vector3D<typename Backend::precision_v> const &point) const
0367 {
0368   typedef typename Backend::precision_v Float_t;
0369   // algorithm: calculate projections to both planes
0370   // return maximum / minimum depending on fAngle < PI or not
0371   // assuming that we have z wedge and the planes pass through the origin
0372
0373   // actually we
0374
0375   Float_t dist1 = point.x() * fNormalVector1.x() + point.y() * fNormalVector1.y();
0376   Float_t dist2 = point.x() * fNormalVector2.x() + point.y() * fNormalVector2.y();
0377
0378   // std::cerr << "d1 " << dist1<<"  "<<point << "\n";
0379   // std::cerr << "d2 " << dist2<<"  "<<point << "\n";
0380
0381   if (fDPhi < kPi) {
0382     // Float_t disttocorner = Sqrt(point.x()*point.x() + point.y()*point.y());
0383     // commented out DistanceToCorner in order to not have a differences with Geant4 and Root
0384     return Max(-1 * dist1, -1 * dist2);
0385   } else {
0386     return Min(-1 * dist1, -1 * dist2);
0387   }
0388 }
0389
0390 template <class Backend>
0391 VECCORE_ATT_HOST_DEVICE void Wedge::DistanceToIn(Vector3D<typename Backend::precision_v> const &point,
0392                                                  Vector3D<typename Backend::precision_v> const &dir,
0393                                                  typename Backend::precision_v &distWedge1,
0394                                                  typename Backend::precision_v &distWedge2) const
0395 {
0396   typedef typename Backend::precision_v Float_t;
0397   typedef typename Backend::bool_v Bool_t;
0398   // algorithm::first calculate projections of direction to both planes,
0399   // then calculate real distance along given direction,
0400   // distance can be negative
0401
0402   distWedge1 = kInfLength;
0403   distWedge2 = kInfLength;
0404
0405   Float_t comp1 = dir.x() * fNormalVector1.x() + dir.y() * fNormalVector1.y();
0406   Float_t comp2 = dir.x() * fNormalVector2.x() + dir.y() * fNormalVector2.y();
0407
0408   Bool_t cmp1 = comp1 > 0.;
0409   if (!vecCore::MaskEmpty(cmp1)) {
0410     Float_t tmp = -(point.x() * fNormalVector1.x() + point.y() * fNormalVector1.y()) / comp1;
0411     vecCore::MaskedAssign(distWedge1, cmp1 && tmp > -kTolerance, tmp);
0412   }
0413   Bool_t cmp2 = comp2 > 0.;
0414   if (!vecCore::MaskEmpty(cmp2)) {
0415     Float_t tmp = -(point.x() * fNormalVector2.x() + point.y() * fNormalVector2.y()) / comp2;
0416     vecCore::MaskedAssign(distWedge2, cmp2 && tmp > -kTolerance, tmp);
0417   }
0418
0419   // std::cerr << "c1 " << comp1 <<" d1="<<distWedge1<<" p="<<point<< "\n";
0420   // std::cerr << "c2 " << comp2 <<" d2="<<distWedge2<< "\n";
0421 }
0422
0423 template <class Backend>
0424 VECCORE_ATT_HOST_DEVICE void Wedge::DistanceToOut(Vector3D<typename Backend::precision_v> const &point,
0425                                                   Vector3D<typename Backend::precision_v> const &dir,
0426                                                   typename Backend::precision_v &distWedge1,
0427                                                   typename Backend::precision_v &distWedge2) const
0428 {
0429
0430   typedef typename Backend::precision_v Float_t;
0431   typedef typename Backend::bool_v Bool_t;
0432
0433   // algorithm::first calculate projections of direction to both planes,
0434   // then calculate real distance along given direction,
0435   // distance can be negative
0436
0437   Float_t comp1 = dir.x() * fNormalVector1.x() + dir.y() * fNormalVector1.y();
0438   Float_t comp2 = dir.x() * fNormalVector2.x() + dir.y() * fNormalVector2.y();
0439
0440   // std::cerr << "c1 " << comp1 << "\n";
0441   // std::cerr << "c2 " << comp2 << "\n";
0442   distWedge1 = kInfLength;
0443   distWedge2 = kInfLength;
0444
0445   Bool_t cmp1 = comp1 < 0.;
0446   if (!vecCore::MaskEmpty(cmp1)) {
0447     Float_t tmp = -(point.x() * fNormalVector1.x() + point.y() * fNormalVector1.y()) / comp1;
0448     vecCore::MaskedAssign(distWedge1, cmp1 && tmp > -kTolerance, tmp);
0449   }
0450
0451   Bool_t cmp2 = comp2 < 0.;
0452   if (!vecCore::MaskEmpty(cmp2)) {
0453     Float_t tmp = -(point.x() * fNormalVector2.x() + point.y() * fNormalVector2.y()) / comp2;
0454     vecCore::MaskedAssign(distWedge2, cmp2 && tmp > -kTolerance, tmp);
0455   }
0456
0457   // std::cerr << "c1 " << comp1 <<" d1="<<distWedge1<<" "<<point<< "\n";
0458   // std::cerr << "c2 " << comp2 <<" d2=" <<distWedge2<<" "<<point<<"\n";
0459 }
0460 } // namespace VECGEOM_IMPL_NAMESPACE
0461 } // namespace vecgeom
0462
0463 #endif /* VECGEOM_VOLUMES_WEDGE_H_ */