Back to home page

EIC code displayed by LXR
 
 

    

File indexing completed on 2025-03-13 09:29:36 

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