Back to home page

EIC code displayed by LXR
 
 

    

File indexing completed on 2026-09-23 09:30:02 

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