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 LoopUnplacedVolumeImplHelper<
0084           PolyhedronImplementation<Polyhedron::EInnerRadii::kGeneric, Polyhedron::EPhiCutout::kGeneric>>,
0085       public AlignedBase {
0086
0087 private:
0088   PolyhedronStruct<Precision> fPoly; ///< Structure holding polyhedron data
0089
0090 public:
0091   UnplacedPolyhedron() : fPoly() { ComputeBBox(); }
0092   /// \param sideCount Number of sides along phi in each Z-segment.
0093   /// \param zPlaneCount Number of Z-planes to draw segments between. The number
0094   ///                    of segments will always be this number minus one.
0095   /// \param zPlanes Z-coordinates of each Z-plane to draw segments between.
0096   /// \param rMin Radius to the sides (not to the corners!) of the inner shell
0097   ///             for the corresponding Z-plane.
0098   /// \param rMin Radius to the sides (not to the corners!) of the outer shell
0099   ///             for the corresponding Z-plane.
0100   UnplacedPolyhedron(const int sideCount, const int zPlaneCount, Precision const zPlanes[], Precision const rMin[],
0101                      Precision const rMax[]);
0102
0103   /// \param phiStart Angle in phi of first corner. This will be one phi angle
0104   ///                 of the phi cutout, if any cutout is specified. Specified
0105   ///                 in radians.
0106   /// \param phiDelta Total angle in phi over which the sides of each segment
0107   ///                 will be drawn. When added to the starting angle, this will
0108   ///                 mark one of the angles of the phi cutout, if any cutout is
0109   ///                 specified.
0110   /// \param sideCount Number of sides along phi in each Z-segment.
0111   /// \param zPlaneCount Number of Z-planes to draw segments between. The number
0112   ///                    of segments will always be this number minus one.
0113   /// \param zPlanes Z-coordinates of each Z-plane to draw segments between.
0114   /// \param rMin Radius to the sides (not to the corners!) of the inner shell
0115   ///             for the corresponding Z-plane.
0116   /// \param rMin Radius to the sides (not to the corners!) of the outer shell
0117   ///             for the corresponding Z-plane.
0118   VECCORE_ATT_HOST_DEVICE
0119   UnplacedPolyhedron(Precision phiStart, Precision phiDelta, const int sideCount, const int zPlaneCount,
0120                      Precision const zPlanes[], Precision const rMin[], Precision const rMax[]);
0121
0122   /// Alternative constructor, required for integration with Geant4.
0123   /// This constructor mirrors one in UnplacedPolycone(), for which the r[],z[] idea makes more sense.
0124   /// Input must be such that r[i],z[i] arrays describe the outer,inner or inner,outer envelope of the
0125   /// polyhedron, after connecting all adjacent points, and closing the polygon by connecting last -> first points.
0126   /// Hence z[] array must be symmetrical: z[0..Nz] = z[2Nz, 2Nz-1, ..., Nz+1], where Nz = zPlaneCount.
0127   ///
0128   /// \param phiStart Angle in phi of first corner. This will be one phi angle of the phi cutout, if any
0129   ///                 cutout is specified. Specified in radians.
0130   /// \param phiDelta Total angle in phi over which the sides of each segment will be drawn. When added to the
0131   ///                 starting angle, this will mark one of the angles of the phi cutout, if a cutout is specified.
0132   /// \param sideCount Number of sides along phi in each Z-segment.
0133   /// \param verticesCount Number of vertices describing inner/outer shape in r/z coordinates. The number
0134   ///                    of Z planes will be half of this number.
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 for the corresponding Z-plane.
0137   /// \param rMax Radius to the sides (not to the corners!) of the outer shell for the corresponding Z-plane.
0138   UnplacedPolyhedron(Precision phiStart, Precision phiDelta, const int sideCount, const int verticesCount,
0139                      Precision const r[], Precision const z[]);
0140
0141   VECCORE_ATT_HOST_DEVICE
0142   virtual ~UnplacedPolyhedron() {}
0143
0144   VECCORE_ATT_HOST_DEVICE
0145   PolyhedronStruct<Precision> const &GetStruct() const { return fPoly; }
0146
0147   VECCORE_ATT_HOST_DEVICE
0148   VECGEOM_FORCE_INLINE
0149   int GetSideCount() const { return fPoly.fSideCount; }
0150
0151   VECCORE_ATT_HOST_DEVICE
0152   VECGEOM_FORCE_INLINE
0153   int GetZSegmentCount() const { return fPoly.fZSegments.size(); }
0154
0155   VECCORE_ATT_HOST_DEVICE
0156   VECGEOM_FORCE_INLINE
0157   bool HasInnerRadii() const { return fPoly.fHasInnerRadii; }
0158
0159   VECCORE_ATT_HOST_DEVICE
0160   VECGEOM_FORCE_INLINE
0161   bool HasPhiCutout() const { return fPoly.fHasPhiCutout; }
0162
0163   VECCORE_ATT_HOST_DEVICE
0164   VECGEOM_FORCE_INLINE
0165   bool HasLargePhiCutout() const { return fPoly.fHasLargePhiCutout; }
0166
0167   VECCORE_ATT_HOST_DEVICE
0168   VECGEOM_FORCE_INLINE
0169   ZSegment const &GetZSegment(int i) const { return fPoly.fZSegments[i]; }
0170
0171   VECCORE_ATT_HOST_DEVICE
0172   VECGEOM_FORCE_INLINE
0173   Array<ZSegment> const &GetZSegments() const { return fPoly.fZSegments; }
0174
0175   VECCORE_ATT_HOST_DEVICE
0176   VECGEOM_FORCE_INLINE
0177   Precision GetZPlane(int i) const { return fPoly.fZPlanes[i]; }
0178
0179   VECCORE_ATT_HOST_DEVICE
0180   VECGEOM_FORCE_INLINE
0181   Array<Precision> const &GetZPlanes() const { return fPoly.fZPlanes; }
0182
0183   VECCORE_ATT_HOST_DEVICE
0184   VECGEOM_FORCE_INLINE
0185   Array<Precision> const &GetRMin() const { return fPoly.fRMin; }
0186
0187   VECCORE_ATT_HOST_DEVICE
0188   VECGEOM_FORCE_INLINE
0189   Array<Precision> const &GetRMax() const { return fPoly.fRMax; }
0190
0191   VECCORE_ATT_HOST_DEVICE
0192   VECGEOM_FORCE_INLINE
0193   Vector3D<Precision> GetPhiSection(int i) const { return fPoly.fPhiSections[i]; }
0194
0195   VECCORE_ATT_HOST_DEVICE
0196   VECGEOM_FORCE_INLINE
0197   SOA3D<Precision> const &GetPhiSections() const { return fPoly.fPhiSections; }
0198
0199   VECCORE_ATT_HOST_DEVICE
0200   VECGEOM_FORCE_INLINE
0201   evolution::Wedge const &GetPhiWedge() const { return fPoly.fPhiWedge; }
0202
0203   VECCORE_ATT_HOST_DEVICE
0204   VECGEOM_FORCE_INLINE
0205   TubeStruct<Precision> const &GetBoundingTube() const { return fPoly.fBoundingTube; }
0206
0207   VECCORE_ATT_HOST_DEVICE
0208   VECGEOM_FORCE_INLINE
0209   Precision GetBoundingTubeOffset() const { return fPoly.fBoundingTubeOffset; }
0210
0211 #ifndef VECCORE_CUDA
0212   VECCORE_ATT_HOST_DEVICE
0213   bool Normal(Vector3D<Precision> const &point, Vector3D<Precision> &normal) const override;
0214 #endif
0215
0216   // calculate array of triangle spanned by points v1,v2,v3
0217   // TODO: this function has nothing to do with a Polyhedron. It should live somewhere else ( indeed: the Quadriteral
0218   // seems to have such a function, too )
0219   VECCORE_ATT_HOST_DEVICE
0220   Precision GetTriangleArea(Vector3D<Precision> const &v1, Vector3D<Precision> const &v2,
0221                             Vector3D<Precision> const &v3) const;
0222
0223   Precision Capacity() const override;
0224
0225   Precision SurfaceArea() const override;
0226 #ifndef VECCORE_CUDA
0227   Precision DistanceSquarePointToSegment(Vector3D<Precision> &v1, Vector3D<Precision> &v2,
0228                                          const Vector3D<Precision> &p) const;
0229   bool InsideTriangle(Vector3D<Precision> &v1, Vector3D<Precision> &v2, Vector3D<Precision> &v3,
0230                       const Vector3D<Precision> &p) const;
0231
0232   // returns a random point inside the triangle described by v1,v2,v3
0233   // TODO: this function has nothing to do with a Polyhedron. It should live somewhere else ( indeed: the Quadriteral
0234   // seems to have such a function, too )
0235   Vector3D<Precision> GetPointOnTriangle(Vector3D<Precision> const &v1, Vector3D<Precision> const &v2,
0236                                          Vector3D<Precision> const &v3) const;
0237
0238   void Extent(Vector3D<Precision> &aMin, Vector3D<Precision> &aMax) const override;
0239
0240   Vector3D<Precision> SamplePointOnSurface() const override;
0241
0242   std::string GetEntityType() const { return "Polyhedron"; }
0243 #endif // !VECCORE_CUDA
0244
0245   /// Not a stored value, and should not be called from performance critical code.
0246   /// \return The angle along phi where the first corner is placed, specified in radians.
0247   VECCORE_ATT_HOST_DEVICE
0248   VECGEOM_FORCE_INLINE
0249   Precision GetPhiStart() const { return fPoly.fPhiStart; }
0250
0251   /// Not a stored value, and should not be called from performance critical code.
0252   /// \return The angle along phi where the last corner is placed, specified in degrees.
0253   VECCORE_ATT_HOST_DEVICE
0254   VECGEOM_FORCE_INLINE
0255   Precision GetPhiEnd() const { return fPoly.fPhiStart + fPoly.fPhiDelta; }
0256
0257   /// Not a stored value, and should not be called from performance critical code.
0258   /// \return The difference in angle along phi between the last corner and the first corner.
0259   VECCORE_ATT_HOST_DEVICE
0260   VECGEOM_FORCE_INLINE
0261   Precision GetPhiDelta() const { return fPoly.fPhiDelta; }
0262
0263   // \return the number of quadrilaterals (including triangles) that this
0264   // polyhedra consists of; this should be all visible surfaces except the endcaps
0265   VECCORE_ATT_HOST_DEVICE
0266   int GetNQuadrilaterals() const;
0267
0268   // reconstructs fZPlanes, fRmin, fRMax from Quadrilaterals
0269   template <typename PushableContainer>
0270   void ReconstructSectionArrays(PushableContainer &zplanes, PushableContainer &rmin, PushableContainer &rmax) const
0271   {
0272     // iterate over sections;
0273     // pick one inner quadrilateral and one outer quadrilateral
0274     // reconstruct rmin, rmax and z from these
0275
0276     // TODO: this might not yet be correct when we have degenerate
0277     // z-plane values
0278
0279     AOS3D<Precision> const *innercorners;
0280     AOS3D<Precision> const *outercorners;
0281
0282     // lambda function to recalculate the radii
0283     auto getradius = [](Vector3D<Precision> const &a, Vector3D<Precision> const &b) {
0284       return std::sqrt(a.Perp2() - (a - b).Mag2() / 4.);
0285     };
0286
0287     Array<ZSegment>::const_iterator s;
0288     Array<ZSegment>::const_iterator end = fPoly.fZSegments.cend();
0289
0290     for (s = fPoly.fZSegments.cbegin(); s != end; ++s) {
0291       outercorners          = (*s).outer.GetCorners();
0292       Vector3D<Precision> a = outercorners[0][0];
0293       Vector3D<Precision> b = outercorners[1][0];
0294       rmax.push_back(getradius(a, b));
0295       zplanes.push_back(a.z());
0296
0297       if (fPoly.fHasInnerRadii) {
0298         innercorners = (*s).inner.GetCorners();
0299         a            = innercorners[0][0];
0300         b            = innercorners[1][0];
0301         rmin.push_back(getradius(a, b));
0302       } else {
0303         rmin.push_back(0.);
0304       }
0305     }
0306     // for last segment need to add addidional plane
0307
0308     Vector3D<Precision> a = outercorners[2][0];
0309     Vector3D<Precision> b = outercorners[3][0];
0310     rmax.push_back(getradius(a, b));
0311     zplanes.push_back(a.z());
0312
0313     if (fPoly.fHasInnerRadii) {
0314       a = innercorners[2][0];
0315       b = innercorners[3][0];
0316       rmin.push_back(getradius(a, b));
0317     } else {
0318       rmin.push_back(0.);
0319     }
0320   }
0321
0322   VECCORE_ATT_HOST_DEVICE
0323   void DetectConvexity();
0324
0325   VECCORE_ATT_HOST_DEVICE
0326   virtual void Print() const final;
0327
0328   VECCORE_ATT_HOST_DEVICE
0329   void PrintSegments() const;
0330
0331   virtual void Print(std::ostream &os) const final;
0332
0333 #ifndef VECCORE_CUDA
0334   virtual SolidMesh *CreateMesh3D(Transformation3D const &trans, size_t nSegments) const override;
0335 #endif
0336
0337   template <TranslationCode transCodeT, RotationCode rotCodeT>
0338   VECCORE_ATT_DEVICE
0339   static VPlacedVolume *Create(LogicalVolume const *const logical_volume, Transformation3D const *const transformation,
0340 #ifdef VECCORE_CUDA
0341                                const int id, const int copy_no, const int child_id,
0342 #endif
0343                                VPlacedVolume *const placement = NULL);
0344
0345   VECCORE_ATT_DEVICE
0346   virtual VPlacedVolume *SpecializedVolume(LogicalVolume const *const volume,
0347                                            Transformation3D const *const transformation,
0348                                            const TranslationCode trans_code, const RotationCode rot_code,
0349 #ifdef VECCORE_CUDA
0350                                            const int id, const int copy_no, const int child_id,
0351 #endif
0352                                            VPlacedVolume *const placement = NULL) const final;
0353   VECGEOM_FORCE_INLINE
0354   virtual int MemorySize() const final { return sizeof(*this); }
0355
0356 #ifdef VECGEOM_CUDA_INTERFACE
0357   virtual size_t DeviceSizeOf() const override { return DevicePtr<cuda::UnplacedPolyhedron>::SizeOf(); }
0358   virtual DevicePtr<cuda::VUnplacedVolume> CopyToGpu() const override;
0359   virtual DevicePtr<cuda::VUnplacedVolume> CopyToGpu(DevicePtr<cuda::VUnplacedVolume> const gpu_ptr) const override;
0360   static void CopyToGpu(std::vector<VUnplacedVolume const *> const & volumes, std::vector<DevicePtr<cuda::VUnplacedVolume>> const & devPtrs);
0361 #endif
0362
0363 private:
0364   // This method does the proper construction of planes and segments.
0365   // Used by multiple constructors.
0366   VECCORE_ATT_HOST_DEVICE
0367   void Initialize(Precision phiStart, Precision phiDelta, const int sideCount, const int zPlaneCount,
0368                   Precision const zPlanes[], Precision const rMin[], Precision const rMax[]);
0369
0370 }; // End class UnplacedPolyhedron
0371
0372 } // namespace VECGEOM_IMPL_NAMESPACE
0373
0374 } // namespace vecgeom
0375
0376 #endif // VECGEOM_VOLUMES_UNPLACEDPOLYHEDRON_H_