VecGeom/volumes/UnplacedPolyhedron.h

0001 /// \file UnplacedPolyhedron.h
0002 /// \author Johannes de Fine Licht (johannes.definelicht@cern.ch)
0003
0004 #ifndef VECGEOM_VOLUMES_UNPLACEDPOLYHEDRON_H_
0005 #define VECGEOM_VOLUMES_UNPLACEDPOLYHEDRON_H_
0006
0007 #include "VecGeom/base/Cuda.h"
0008 #include "VecGeom/base/Global.h"
0009 #include "VecGeom/base/AlignedBase.h"
0010 #include "VecGeom/volumes/UnplacedVolume.h"
0011 #include "VecGeom/volumes/PolyhedronStruct.h"
0012 #include "VecGeom/volumes/kernel/PolyhedronImplementation.h"
0013 #include "VecGeom/volumes/UnplacedVolumeImplHelper.h"
0014
0015 namespace vecgeom {
0016
0017 VECGEOM_DEVICE_FORWARD_DECLARE(class UnplacedPolyhedron;);
0018 VECGEOM_DEVICE_DECLARE_CONV(class, UnplacedPolyhedron);
0019
0020 template <typename Stream>
0021 Stream &operator<<(Stream &st, EInnerRadii a)
0022 {
0023   if (a == EInnerRadii::kFalse) st << "EInnerRadii::kFalse";
0024   if (a == EInnerRadii::kGeneric) st << "EInnerRadii::kGeneric";
0025   if (a == EInnerRadii::kTrue) st << "EInnerRadii::kTrue";
0026   return st;
0027 }
0028
0029 template <typename Stream>
0030 Stream &operator<<(Stream &st, EPhiCutout a)
0031 {
0032   if (a == EPhiCutout::kFalse) st << "EPhiCutout::kFalse";
0033   if (a == EPhiCutout::kGeneric) st << "EPhiCutout::kGeneric";
0034   if (a == EPhiCutout::kTrue) st << "EPhiCutout::kTrue";
0035   if (a == EPhiCutout::kLarge) st << "EPhiCutout::kLarge";
0036   return st;
0037 }
0038
0039 inline namespace VECGEOM_IMPL_NAMESPACE {
0040
0041 /// \class UnplacedPolyhedron
0042 /// \brief A series of regular n-sided segments along the Z-axis with varying
0043 ///        radii and mutual distance in Z.
0044 ///
0045 ///
0046 /// ---- Cross section of single Z segment ----
0047 ///
0048 /// R/Phi--->    -o- Z
0049 /// |        ________________
0050 /// v       /        ^      .\,
0051 ///        /    rMax |     .  \,
0052 ///       /          |    . <------ fPhiSections[1]
0053 ///      /       ____|___.      \,
0054 ///     /       /    ^   \       \,
0055 ///    /       /     |rMin\       \,
0056 ///   /       /      |     \_______\ phiStart/fPhiSections[0]
0057 ///   \       \                ^
0058 ///    \       \               |
0059 ///     \       \________      |
0060 ///      \           ^   \<---fZSegments.phi
0061 ///      fZSegments.inner \,
0062 ///        \               \,
0063 ///         \_______________\,
0064 ///           ^              phiStart+phiDelta/fPhiSections[n-1]
0065 /// zSegment.outer
0066 ///
0067 ///
0068 /// ---- Segments along Z ----
0069 ///
0070 ///                          fZPlanes[size-1]
0071 /// fRMax[1]_____fRMax[2] __       |
0072 ///       /|     |\     /|  \___   v
0073 ///      / |     | \___/ |  |   |\.
0074 ///     |  |     | |   | |  |   | \.
0075 ///     |  |     | |   | |  |   |  |
0076 ///     |  |     | |___| |  |   | /
0077 ///      \ |     | /   \ |  |___|/    ^ R/Phi
0078 ///     ^ \|_____|/     \|__/         |
0079 ///     |                             |     Z
0080 ///     fZPlanes[0]/fRMax[0]           ----->
0081
0082 class UnplacedPolyhedron
0083     : public UnplacedVolumeImplHelper<
0084           PolyhedronImplementation<Polyhedron::EInnerRadii::kGeneric, Polyhedron::EPhiCutout::kGeneric>>,
0085       public AlignedBase {
0086
0087 private:
0088   PolyhedronStruct<Precision> *fPoly{nullptr}; ///< Structure holding polyhedron data
0089   bool fAllocated{false};
0090   const void *fBuffer{nullptr};
0091
0092 public:
0093   UnplacedPolyhedron() = default;
0094   /// \param sideCount Number of sides along phi in each Z-segment.
0095   /// \param zPlaneCount Number of Z-planes to draw segments between. The number
0096   ///                    of segments will always be this number minus one.
0097   /// \param zPlanes Z-coordinates of each Z-plane to draw segments between.
0098   /// \param rMin Radius to the sides (not to the corners!) of the inner shell
0099   ///             for the corresponding Z-plane.
0100   /// \param rMin Radius to the sides (not to the corners!) of the outer shell
0101   ///             for the corresponding Z-plane.
0102   UnplacedPolyhedron(const int sideCount, const int zPlaneCount, Precision const zPlanes[], Precision const rMin[],
0103                      Precision const rMax[]);
0104
0105   /// \param phiStart Angle in phi of first corner. This will be one phi angle
0106   ///                 of the phi cutout, if any cutout is specified. Specified
0107   ///                 in radians.
0108   /// \param phiDelta Total angle in phi over which the sides of each segment
0109   ///                 will be drawn. When added to the starting angle, this will
0110   ///                 mark one of the angles of the phi cutout, if any cutout is
0111   ///                 specified.
0112   /// \param sideCount Number of sides along phi in each Z-segment.
0113   /// \param zPlaneCount Number of Z-planes to draw segments between. The number
0114   ///                    of segments will always be this number minus one.
0115   /// \param zPlanes Z-coordinates of each Z-plane to draw segments between.
0116   /// \param rMin Radius to the sides (not to the corners!) of the inner shell
0117   ///             for the corresponding Z-plane.
0118   /// \param rMax Radius to the sides (not to the corners!) of the outer shell
0119   ///             for the corresponding Z-plane.
0120   VECCORE_ATT_HOST_DEVICE
0121   UnplacedPolyhedron(Precision phiStart, Precision phiDelta, const int sideCount, const int zPlaneCount,
0122                      Precision const zPlanes[], Precision const rMin[], Precision const rMax[]);
0123
0124   /// @brief CPU constructor, constructing the PolyhedronStruct in a preallocated buffer
0125   /// \param phiStart Angle in phi of first corner. This will be one phi angle
0126   ///                 of the phi cutout, if any cutout is specified. Specified
0127   ///                 in radians.
0128   /// \param phiDelta Total angle in phi over which the sides of each segment
0129   ///                 will be drawn. When added to the starting angle, this will
0130   ///                 mark one of the angles of the phi cutout, if any cutout is
0131   ///                 specified.
0132   /// \param sideCount Number of sides along phi in each Z-segment.
0133   /// \param zPlaneCount Number of Z-planes to draw segments between. The number
0134   ///                    of segments will always be this number minus one.
0135   /// \param zPlanes Z-coordinates of each Z-plane to draw segments between.
0136   /// \param rMin Radius to the sides (not to the corners!) of the inner shell
0137   ///             for the corresponding Z-plane.
0138   /// \param rMax Radius to the sides (not to the corners!) of the outer shell
0139   ///             for the corresponding Z-plane.
0140   /// @param buff_size Size of the buffer to hold the PolyhedronStruct
0141   /// @param buffer Buffer on GPU already pre-allocated
0142   VECCORE_ATT_HOST_DEVICE
0143   UnplacedPolyhedron(Precision phiStart, Precision phiDelta, const int sideCount, const int zPlaneCount,
0144                      Precision const zPlanes[], Precision const rMin[], Precision const rMax[], size_t const buff_size,
0145                      void const *buffer);
0146
0147   /// Alternative constructor, required for integration with Geant4.
0148   /// This constructor mirrors one in UnplacedPolycone(), for which the r[],z[] idea makes more sense.
0149   /// Input must be such that r[i],z[i] arrays describe the outer,inner or inner,outer envelope of the
0150   /// polyhedron, after connecting all adjacent points, and closing the polygon by connecting last -> first points.
0151   /// Hence z[] array must be symmetrical: z[0..Nz] = z[2Nz, 2Nz-1, ..., Nz+1], where Nz = zPlaneCount.
0152   ///
0153   /// \param phiStart Angle in phi of first corner. This will be one phi angle of the phi cutout, if any
0154   ///                 cutout is specified. Specified in radians.
0155   /// \param phiDelta Total angle in phi over which the sides of each segment will be drawn. When added to the
0156   ///                 starting angle, this will mark one of the angles of the phi cutout, if a cutout is specified.
0157   /// \param sideCount Number of sides along phi in each Z-segment.
0158   /// \param verticesCount Number of vertices describing inner/outer shape in r/z coordinates. The number
0159   ///                    of Z planes will be half of this number.
0160   /// \param zPlanes Z-coordinates of each Z-plane to draw segments between.
0161   /// \param rMin Radius to the sides (not to the corners!) of the inner shell for the corresponding Z-plane.
0162   /// \param rMax Radius to the sides (not to the corners!) of the outer shell for the corresponding Z-plane.
0163   UnplacedPolyhedron(Precision phiStart, Precision phiDelta, const int sideCount, const int verticesCount,
0164                      Precision const r[], Precision const z[]);
0165
0166   VECCORE_ATT_HOST_DEVICE
0167   virtual ~UnplacedPolyhedron()
0168   {
0169     if (fAllocated) delete[] (char *)fBuffer;
0170   }
0171
0172   VECCORE_ATT_HOST_DEVICE
0173   VECGEOM_FORCE_INLINE
0174   virtual ESolidType GetType() const override { return ESolidType::polyhedron; }
0175
0176   VECCORE_ATT_HOST_DEVICE
0177   PolyhedronStruct<Precision> const &GetStruct() const { return *fPoly; }
0178
0179   VECCORE_ATT_HOST_DEVICE
0180   VECGEOM_FORCE_INLINE
0181   int GetSideCount() const { return fPoly->fSideCount; }
0182
0183   VECCORE_ATT_HOST_DEVICE
0184   VECGEOM_FORCE_INLINE
0185   int GetZSegmentCount() const { return fPoly->fZSegments.size(); }
0186
0187   VECCORE_ATT_HOST_DEVICE
0188   VECGEOM_FORCE_INLINE
0189   bool HasInnerRadii() const { return fPoly->fHasInnerRadii; }
0190
0191   VECCORE_ATT_HOST_DEVICE
0192   VECGEOM_FORCE_INLINE
0193   bool HasPhiCutout() const { return fPoly->fHasPhiCutout; }
0194
0195   VECCORE_ATT_HOST_DEVICE
0196   VECGEOM_FORCE_INLINE
0197   bool HasLargePhiCutout() const { return fPoly->fHasLargePhiCutout; }
0198
0199   VECCORE_ATT_HOST_DEVICE
0200   VECGEOM_FORCE_INLINE
0201   ZSegment const &GetZSegment(int i) const { return fPoly->fZSegments[i]; }
0202
0203   VECCORE_ATT_HOST_DEVICE
0204   VECGEOM_FORCE_INLINE
0205   Array<ZSegment> const &GetZSegments() const { return fPoly->fZSegments; }
0206
0207   VECCORE_ATT_HOST_DEVICE
0208   VECGEOM_FORCE_INLINE
0209   Precision GetZPlane(int i) const { return fPoly->fZPlanes[i]; }
0210
0211   VECCORE_ATT_HOST_DEVICE
0212   VECGEOM_FORCE_INLINE
0213   Array<Precision> const &GetZPlanes() const { return fPoly->fZPlanes; }
0214
0215   VECCORE_ATT_HOST_DEVICE
0216   VECGEOM_FORCE_INLINE
0217   Array<Precision> const &GetRMin() const { return fPoly->fRMin; }
0218
0219   VECCORE_ATT_HOST_DEVICE
0220   VECGEOM_FORCE_INLINE
0221   Array<Precision> const &GetRMax() const { return fPoly->fRMax; }
0222
0223   VECCORE_ATT_HOST_DEVICE
0224   VECGEOM_FORCE_INLINE
0225   Vector3D<Precision> GetPhiSection(int i) const { return fPoly->fPhiSections[i]; }
0226
0227   VECCORE_ATT_HOST_DEVICE
0228   VECGEOM_FORCE_INLINE
0229   SOA3D<Precision> const &GetPhiSections() const { return fPoly->fPhiSections; }
0230
0231   VECCORE_ATT_HOST_DEVICE
0232   VECGEOM_FORCE_INLINE
0233   evolution::Wedge const &GetPhiWedge() const { return fPoly->fPhiWedge; }
0234
0235   VECCORE_ATT_HOST_DEVICE
0236   VECGEOM_FORCE_INLINE
0237   TubeStruct<Precision> const &GetBoundingTube() const { return fPoly->fBoundingTube; }
0238
0239   VECCORE_ATT_HOST_DEVICE
0240   VECGEOM_FORCE_INLINE
0241   Precision GetBoundingTubeOffset() const { return fPoly->fBoundingTubeOffset; }
0242
0243   VECCORE_ATT_HOST_DEVICE
0244   bool Normal(Vector3D<Precision> const &point, Vector3D<Precision> &normal) const override;
0245
0246   // calculate array of triangle spanned by points v1,v2,v3
0247   // TODO: this function has nothing to do with a Polyhedron. It should live somewhere else ( indeed: the Quadriteral
0248   // seems to have such a function, too )
0249   VECCORE_ATT_HOST_DEVICE
0250   Precision GetTriangleArea(Vector3D<Precision> const &v1, Vector3D<Precision> const &v2,
0251                             Vector3D<Precision> const &v3) const;
0252
0253   Precision Capacity() const override;
0254
0255   Precision SurfaceArea() const override;
0256
0257   VECCORE_ATT_HOST_DEVICE
0258   void Extent(Vector3D<Precision> &aMin, Vector3D<Precision> &aMax) const override;
0259
0260 #ifndef VECCORE_CUDA
0261   Precision DistanceSquarePointToSegment(Vector3D<Precision> &v1, Vector3D<Precision> &v2,
0262                                          const Vector3D<Precision> &p) const;
0263   bool InsideTriangle(Vector3D<Precision> &v1, Vector3D<Precision> &v2, Vector3D<Precision> &v3,
0264                       const Vector3D<Precision> &p) const;
0265
0266   // returns a random point inside the triangle described by v1,v2,v3
0267   // TODO: this function has nothing to do with a Polyhedron. It should live somewhere else ( indeed: the Quadriteral
0268   // seems to have such a function, too )
0269   Vector3D<Precision> GetPointOnTriangle(Vector3D<Precision> const &v1, Vector3D<Precision> const &v2,
0270                                          Vector3D<Precision> const &v3) const;
0271
0272   Vector3D<Precision> SamplePointOnSurface() const override;
0273
0274   std::string GetEntityType() const { return "Polyhedron"; }
0275 #endif // !VECCORE_CUDA
0276
0277   /// Not a stored value, and should not be called from performance critical code.
0278   /// \return The angle along phi where the first corner is placed, specified in radians.
0279   VECCORE_ATT_HOST_DEVICE
0280   VECGEOM_FORCE_INLINE
0281   Precision GetPhiStart() const { return fPoly->fPhiStart; }
0282
0283   /// Not a stored value, and should not be called from performance critical code.
0284   /// \return The angle along phi where the last corner is placed, specified in degrees.
0285   VECCORE_ATT_HOST_DEVICE
0286   VECGEOM_FORCE_INLINE
0287   Precision GetPhiEnd() const { return fPoly->fPhiStart + fPoly->fPhiDelta; }
0288
0289   /// Not a stored value, and should not be called from performance critical code.
0290   /// \return The difference in angle along phi between the last corner and the first corner.
0291   VECCORE_ATT_HOST_DEVICE
0292   VECGEOM_FORCE_INLINE
0293   Precision GetPhiDelta() const { return fPoly->fPhiDelta; }
0294
0295   // \return the number of quadrilaterals (including triangles) that this
0296   // polyhedra consists of; this should be all visible surfaces except the endcaps
0297   VECCORE_ATT_HOST_DEVICE
0298   int GetNQuadrilaterals() const;
0299
0300   // reconstructs fZPlanes, fRmin, fRMax from Quadrilaterals
0301   template <typename PushableContainer>
0302   void ReconstructSectionArrays(PushableContainer &zplanes, PushableContainer &rmin, PushableContainer &rmax) const
0303   {
0304     // iterate over sections;
0305     // pick one inner quadrilateral and one outer quadrilateral
0306     // reconstruct rmin, rmax and z from these
0307
0308     // TODO: this might not yet be correct when we have degenerate
0309     // z-plane values
0310
0311     AOS3D<Precision> const *innercorners;
0312     AOS3D<Precision> const *outercorners;
0313
0314     // lambda function to recalculate the radii
0315     auto getradius = [](Vector3D<Precision> const &a, Vector3D<Precision> const &b) {
0316       return std::sqrt(a.Perp2() - (a - b).Mag2() / 4.);
0317     };
0318
0319     Array<ZSegment>::const_iterator s;
0320     Array<ZSegment>::const_iterator end = fPoly->fZSegments.cend();
0321
0322     for (s = fPoly->fZSegments.cbegin(); s != end; ++s) {
0323       outercorners          = (*s).outer.GetCorners();
0324       Vector3D<Precision> a = outercorners[0][0];
0325       Vector3D<Precision> b = outercorners[1][0];
0326       rmax.push_back(getradius(a, b));
0327       zplanes.push_back(a.z());
0328
0329       if (fPoly->fHasInnerRadii) {
0330         innercorners = (*s).inner.GetCorners();
0331         a            = innercorners[0][0];
0332         b            = innercorners[1][0];
0333         rmin.push_back(getradius(a, b));
0334       } else {
0335         rmin.push_back(0.);
0336       }
0337     }
0338     // for last segment need to add addidional plane
0339
0340     Vector3D<Precision> a = outercorners[2][0];
0341     Vector3D<Precision> b = outercorners[3][0];
0342     rmax.push_back(getradius(a, b));
0343     zplanes.push_back(a.z());
0344
0345     if (fPoly->fHasInnerRadii) {
0346       a = innercorners[2][0];
0347       b = innercorners[3][0];
0348       rmin.push_back(getradius(a, b));
0349     } else {
0350       rmin.push_back(0.);
0351     }
0352   }
0353
0354   VECCORE_ATT_HOST_DEVICE
0355   void DetectConvexity();
0356
0357   VECCORE_ATT_HOST_DEVICE
0358   virtual void Print() const final;
0359
0360   VECCORE_ATT_HOST_DEVICE
0361   void PrintSegments() const;
0362
0363   virtual void Print(std::ostream &os) const final;
0364
0365 #ifndef VECCORE_CUDA
0366   virtual SolidMesh *CreateMesh3D(Transformation3D const &trans, size_t nSegments) const override;
0367 #endif
0368
0369   VECCORE_ATT_DEVICE
0370   static VPlacedVolume *Create(LogicalVolume const *const logical_volume, Transformation3D const *const transformation,
0371 #ifdef VECCORE_CUDA
0372                                const int id, const int copy_no, const int child_id,
0373 #endif
0374                                VPlacedVolume *const placement = NULL);
0375
0376   VECCORE_ATT_DEVICE
0377   virtual VPlacedVolume *SpecializedVolume(LogicalVolume const *const volume,
0378                                            Transformation3D const *const transformation,
0379 #ifdef VECCORE_CUDA
0380                                            const int id, const int copy_no, const int child_id,
0381 #endif
0382                                            VPlacedVolume *const placement = NULL) const final;
0383   VECGEOM_FORCE_INLINE
0384   virtual int MemorySize() const final { return sizeof(*this); }
0385
0386 #ifdef VECGEOM_CUDA_INTERFACE
0387   virtual size_t DeviceSizeOf() const override { return DevicePtr<cuda::UnplacedPolyhedron>::SizeOf(); }
0388   virtual DevicePtr<cuda::VUnplacedVolume> CopyToGpu() const override;
0389   virtual DevicePtr<cuda::VUnplacedVolume> CopyToGpu(DevicePtr<cuda::VUnplacedVolume> const gpu_ptr) const override;
0390   static void CopyToGpu(std::vector<VUnplacedVolume const *> const &volumes,
0391                         std::vector<DevicePtr<cuda::VUnplacedVolume>> const &devPtrs);
0392 #endif
0393
0394 private:
0395   // This method does the proper construction of planes and segments.
0396   // Used by multiple constructors.
0397   VECCORE_ATT_HOST_DEVICE
0398   void Initialize(Precision phiStart, Precision phiDelta, const int sideCount, const int zPlaneCount,
0399                   Precision const zPlanes[], Precision const rMin[], Precision const rMax[]);
0400
0401 }; // End class UnplacedPolyhedron
0402
0403 } // namespace VECGEOM_IMPL_NAMESPACE
0404
0405 } // namespace vecgeom
0406
0407 #endif // VECGEOM_VOLUMES_UNPLACEDPOLYHEDRON_H_