VecGeom/volumes/Rectangles.h

0001 /// \file Rectangles.h
0002 /// \author Johannes de Fine Licht (johannes.definelicht@cern.ch)
0003
0004 #ifndef VECGEOM_VOLUMES_RECTANGLES_H_
0005 #define VECGEOM_VOLUMES_RECTANGLES_H_
0006
0007 #include "VecGeom/base/Global.h"
0008
0009 #include "VecGeom/base/AlignedBase.h"
0010 #include "VecGeom/base/SOA3D.h"
0011 #include "VecGeom/base/Vector3D.h"
0012 #include "VecGeom/volumes/Planes.h"
0013
0014 namespace vecgeom {
0015
0016 VECGEOM_DEVICE_FORWARD_DECLARE(class Rectangles;);
0017 VECGEOM_DEVICE_DECLARE_CONV(class, Rectangles);
0018
0019 inline namespace VECGEOM_IMPL_NAMESPACE {
0020
0021 /// \class Rectangles
0022 ///
0023 /// \brief Stores a number of rectangles in SOA form to allow vectorized
0024 ///        operations.
0025 ///
0026 /// To allow efficient computation, two corners, one normalized side and
0027 /// the plane equation in which the rectangle lies is stored.
0028 /// If the set of rectangles are assumed to be convex, the convex methods
0029 /// can be called for faster computation, falling back on the implementations
0030 /// for planes.
0031 class Rectangles : public AlignedBase {
0032
0033   //   fCorners[0]
0034   //             p0----------
0035   //       |      |         |
0036   // (Normalized) |         |
0037   //    fSides    |   -o-   |
0038   //       |      |         |
0039   //       v      |         |
0040   //             p1---------p2
0041   //                          fCorners[1]
0042
0043 private:
0044   Planes fPlanes;
0045   SOA3D<Precision> fSides;
0046   SOA3D<Precision> fCorners[2];
0047
0048 public:
0049   typedef SOA3D<Precision> Corners_t[2];
0050
0051   VECCORE_ATT_HOST_DEVICE
0052   Rectangles(int size);
0053
0054   VECCORE_ATT_HOST_DEVICE
0055   ~Rectangles();
0056
0057   VECCORE_ATT_HOST_DEVICE
0058   VECGEOM_FORCE_INLINE
0059   int size() const;
0060
0061   VECCORE_ATT_HOST_DEVICE
0062   VECGEOM_FORCE_INLINE
0063   Vector3D<Precision> GetNormal(int i) const;
0064
0065   VECCORE_ATT_HOST_DEVICE
0066   VECGEOM_FORCE_INLINE
0067   SOA3D<Precision> const &GetNormals() const;
0068
0069   VECCORE_ATT_HOST_DEVICE
0070   VECGEOM_FORCE_INLINE
0071   Precision GetDistance(int i) const;
0072
0073   VECCORE_ATT_HOST_DEVICE
0074   VECGEOM_FORCE_INLINE
0075   Array<Precision> const &GetDistances() const;
0076
0077   VECCORE_ATT_HOST_DEVICE
0078   inline Vector3D<Precision> GetCenter(int i) const;
0079
0080   VECCORE_ATT_HOST_DEVICE
0081   inline Vector3D<Precision> GetCorner(int i, int j) const;
0082
0083   VECCORE_ATT_HOST_DEVICE
0084   VECGEOM_FORCE_INLINE
0085   Corners_t const &GetCorners() const;
0086
0087   VECCORE_ATT_HOST_DEVICE
0088   inline Vector3D<Precision> GetSide(int i) const;
0089
0090   VECCORE_ATT_HOST_DEVICE
0091   VECGEOM_FORCE_INLINE
0092   SOA3D<Precision> const &GetSides() const;
0093
0094   VECCORE_ATT_HOST_DEVICE
0095   void Set(int index, Vector3D<Precision> const &p0, Vector3D<Precision> const &p1, Vector3D<Precision> const &p2);
0096
0097   VECCORE_ATT_HOST_DEVICE
0098   inline Precision Distance(Vector3D<Precision> const &point, Vector3D<Precision> const &direction) const;
0099 };
0100
0101 VECCORE_ATT_HOST_DEVICE
0102 int Rectangles::size() const
0103 {
0104   return fPlanes.size();
0105 }
0106
0107 VECCORE_ATT_HOST_DEVICE
0108 Vector3D<Precision> Rectangles::GetNormal(int i) const
0109 {
0110   return fPlanes.GetNormal(i);
0111 }
0112
0113 VECCORE_ATT_HOST_DEVICE
0114 SOA3D<Precision> const &Rectangles::GetNormals() const
0115 {
0116   return fPlanes.GetNormals();
0117 }
0118
0119 VECCORE_ATT_HOST_DEVICE
0120 Precision Rectangles::GetDistance(int i) const
0121 {
0122   return fPlanes.GetDistance(i);
0123 }
0124
0125 VECCORE_ATT_HOST_DEVICE
0126 Array<Precision> const &Rectangles::GetDistances() const
0127 {
0128   return fPlanes.GetDistances();
0129 }
0130
0131 VECCORE_ATT_HOST_DEVICE
0132 Vector3D<Precision> Rectangles::GetCenter(int i) const
0133 {
0134   return -GetDistance(i) * GetNormal(i);
0135 }
0136
0137 VECCORE_ATT_HOST_DEVICE
0138 Vector3D<Precision> Rectangles::GetCorner(int i, int j) const
0139 {
0140   return Vector3D<Precision>(fCorners[i][0][j], fCorners[i][1][j], fCorners[i][2][j]);
0141 }
0142
0143 VECCORE_ATT_HOST_DEVICE
0144 Rectangles::Corners_t const &Rectangles::GetCorners() const
0145 {
0146   return fCorners;
0147 }
0148
0149 VECCORE_ATT_HOST_DEVICE
0150 Vector3D<Precision> Rectangles::GetSide(int i) const
0151 {
0152   return Vector3D<Precision>(fSides[0][i], fSides[1][i], fSides[2][i]);
0153 }
0154
0155 VECCORE_ATT_HOST_DEVICE
0156 SOA3D<Precision> const &Rectangles::GetSides() const
0157 {
0158   return fSides;
0159 }
0160
0161 VECCORE_ATT_HOST_DEVICE
0162 void Rectangles::Set(int index, Vector3D<Precision> const &p0, Vector3D<Precision> const &p1,
0163                      Vector3D<Precision> const &p2)
0164 {
0165
0166   // Store corners and sides
0167   fCorners[0].set(index, p0);
0168   fCorners[1].set(index, p2);
0169   Vector3D<Precision> side = p1 - p0;
0170   side.Normalize();
0171   fSides.set(index, side);
0172
0173   // Compute plane equation to retrieve normal and distance to origin
0174   // ax + by + cz + d = 0
0175   Precision a, b, c, d;
0176   a = p0[1] * (p1[2] - p2[2]) + p1[1] * (p2[2] - p0[2]) + p2[1] * (p0[2] - p1[2]);
0177   b = p0[2] * (p1[0] - p2[0]) + p1[2] * (p2[0] - p0[0]) + p2[2] * (p0[0] - p1[0]);
0178   c = p0[0] * (p1[1] - p2[1]) + p1[0] * (p2[1] - p0[1]) + p2[0] * (p0[1] - p1[1]);
0179   d = -p0[0] * (p1[1] * p2[2] - p2[1] * p1[2]) - p1[0] * (p2[1] * p0[2] - p0[1] * p2[2]) -
0180       p2[0] * (p0[1] * p1[2] - p1[1] * p0[2]);
0181   Vector3D<Precision> normal(a, b, c);
0182   // Normalize the plane equation
0183   // (ax + by + cz + d) / sqrt(a^2 + b^2 + c^2) = 0 =>
0184   // n0*x + n1*x + n2*x + p = 0
0185   Precision inverseLength = 1. / normal.Length();
0186   normal *= inverseLength;
0187   d *= inverseLength;
0188   if (d >= 0) {
0189     // Ensure normal is pointing away from origin
0190     normal = -normal;
0191     d      = -d;
0192   }
0193
0194   fPlanes.Set(index, normal, d);
0195 }
0196
0197 VECCORE_ATT_HOST_DEVICE
0198 Precision Rectangles::Distance(Vector3D<Precision> const &point, Vector3D<Precision> const &direction) const
0199 {
0200   Precision bestDistance = kInfLength;
0201   for (int i = 0, iMax = size(); i < iMax; ++i) {
0202     Vector3D<Precision> normal       = GetNormal(i);
0203     Vector3D<Precision> side         = GetSide(i);
0204     Precision t                      = -(normal.Dot(point) + GetDistance(i)) / normal.Dot(direction);
0205     Vector3D<Precision> intersection = point + t * direction;
0206     Vector3D<Precision> fromP0       = intersection - fCorners[0][i];
0207     Vector3D<Precision> fromP2       = intersection - fCorners[1][i];
0208     if (t >= 0 && t < bestDistance && side.Dot(fromP0) >= 0 && (-side).Dot(fromP2) >= 0) {
0209       bestDistance = t;
0210     }
0211   }
0212   return bestDistance;
0213 }
0214
0215 std::ostream &operator<<(std::ostream &os, Rectangles const &rhs);
0216
0217 } // End inline impl namespace
0218
0219 } // End global namespace
0220
0221 #endif // VECGEOM_VOLUMES_RECTANGLES_H_