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