VecGeom/volumes/UnplacedTrapezoid.h

0001 ///===-- volumes/UnplacedTrapezoid.h ----------------------------------*- C++ -*-===//
0002 ///
0003 /// \file volumes/UnplacedTrapezoid.h
0004 /// \author Guilherme Lima (lima@fnal.gov)
0005 /// \brief This file contains the declaration of the UnplacedTrapezoid class
0006 ///
0007 /// _____________________________________________________________________________
0008 /// A trapezoid is the solid bounded by the following surfaces:
0009 /// - 2 XY-parallel quadrilaterals (trapezoids) cutting the Z axis at Z=-dz and Z=+dz
0010 /// - 4 additional *planar* quadrilaterals intersecting the faces at +/-dz at their edges
0011 ///
0012 /// The easiest way to think about the trapezoid is starting from a box with faces at +/-dz,
0013 /// and moving its eight corners in such a way as to define the two +/-dz trapezoids, without
0014 /// destroying the coplanarity of the other four faces.  Then the (x,y,z) coordinates of the
0015 /// eight corners can be used to build this trapezoid.
0016 ///
0017 /// If the four side faces are not coplanar, a generic trapezoid must be used (GenTrap).
0018 //===------------------------------------------------------------------------===//
0019 ///
0020 /// 140520 G. Lima   Created based on USolids algorithms and vectorized types
0021 /// 160722 G. Lima   Migration to new helpers and VecCore-based scheme
0022
0023 #ifndef VECGEOM_VOLUMES_UNPLACEDTRAPEZOID_H_
0024 #define VECGEOM_VOLUMES_UNPLACEDTRAPEZOID_H_
0025
0026 #include "VecGeom/base/Cuda.h"
0027 #include "VecGeom/base/Global.h"
0028 #include "VecGeom/base/AlignedBase.h"
0029 #include "VecGeom/base/Vector3D.h"
0030 #include "VecGeom/volumes/UnplacedVolumeImplHelper.h"
0031
0032 #include "VecGeom/volumes/TrapezoidStruct.h" // the pure Trapezoid struct
0033 #include "VecGeom/volumes/kernel/TrapezoidImplementation.h"
0034
0035 namespace vecgeom {
0036
0037 VECGEOM_DEVICE_FORWARD_DECLARE(class UnplacedTrapezoid;);
0038 VECGEOM_DEVICE_DECLARE_CONV(class, UnplacedTrapezoid);
0039
0040 inline namespace VECGEOM_IMPL_NAMESPACE {
0041
0042 typedef Vector3D<Precision> TrapCorners[8];
0043
0044 class UnplacedTrapezoid : public UnplacedVolumeImplHelper<TrapezoidImplementation>, public AlignedBase {
0045
0046 private:
0047   TrapezoidStruct<Precision> fTrap;
0048
0049   // variables to store cached values for Capacity and SurfaceArea
0050   // Precision fCubicVolume, fSurfaceArea;
0051
0052 public:
0053   // full constructor
0054   // Note: theta, phi are assumed to be in radians, for compatibility with Geant4
0055   VECCORE_ATT_HOST_DEVICE
0056   UnplacedTrapezoid(const Precision dz, const Precision theta, const Precision phi, const Precision dy1,
0057                     const Precision dx1, const Precision dx2, const Precision Alpha1, const Precision dy2,
0058                     const Precision dx3, const Precision dx4, const Precision Alpha2)
0059       : fTrap(dz, theta, phi, dy1, dx1, dx2, std::tan(Alpha1), dy2, dx3, dx4, std::tan(Alpha2))
0060   {
0061     fGlobalConvexity = true;
0062     MakePlanes();
0063     ComputeBBox();
0064   }
0065
0066   // default constructor
0067   VECCORE_ATT_HOST_DEVICE
0068   UnplacedTrapezoid() : fTrap(0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0.)
0069   {
0070     fGlobalConvexity = true;
0071     ComputeBBox();
0072   }
0073
0074   /// \brief Fast constructor: all parameters from one array
0075   VECCORE_ATT_HOST_DEVICE
0076   UnplacedTrapezoid(Precision const *params)
0077       : UnplacedTrapezoid(params[0], params[1], params[2], params[3], params[4], params[5], params[6], params[7],
0078                           params[8], params[9], params[10])
0079   {
0080     ComputeBBox();
0081   }
0082
0083   /// \brief Constructor based on 8 corner points
0084   // convention: p0(---); p1(+--); p2(-+-); p3(++-); p4(--+); p5(+-+); p6(-++); p7(+++)
0085   VECCORE_ATT_HOST_DEVICE
0086   UnplacedTrapezoid(TrapCorners const corners);
0087
0088   /// \brief Constructor for masquerading a box (test purposes)
0089   VECCORE_ATT_HOST_DEVICE
0090   UnplacedTrapezoid(Precision xbox, Precision ybox, Precision zbox);
0091
0092   /// \brief Constructor required by Geant4
0093   VECCORE_ATT_HOST_DEVICE
0094   // Constructor corresponding to Trd1
0095   UnplacedTrapezoid(Precision dx1, Precision dx2, Precision dy, Precision dz);
0096
0097   // Constructor corresponding to Trd2
0098   /// \brief Constructor for a Trd-like trapezoid
0099   UnplacedTrapezoid(Precision dx1, Precision dx2, Precision dy1, Precision dy2, Precision dz);
0100
0101   /// \brief Constructor for a Parallelepiped-like trapezoid (Note: still to be validated)
0102   UnplacedTrapezoid(Precision dx, Precision dy, Precision dz, Precision alpha, Precision theta, Precision phi);
0103
0104   /// \brief Accessors
0105   /// @{
0106   // VECCORE_ATT_HOST_DEVICE
0107   // TrapParameters const& GetParameters() const { return _params; }
0108
0109   VECCORE_ATT_HOST_DEVICE
0110   VECGEOM_FORCE_INLINE
0111   Precision dz() const { return fTrap.fDz; }
0112
0113   VECCORE_ATT_HOST_DEVICE
0114   VECGEOM_FORCE_INLINE
0115   Precision theta() const { return fTrap.fTheta; }
0116
0117   VECCORE_ATT_HOST_DEVICE
0118   VECGEOM_FORCE_INLINE
0119   Precision phi() const { return fTrap.fPhi; }
0120
0121   VECCORE_ATT_HOST_DEVICE
0122   VECGEOM_FORCE_INLINE
0123   Precision dy1() const { return fTrap.fDy1; }
0124
0125   VECCORE_ATT_HOST_DEVICE
0126   VECGEOM_FORCE_INLINE
0127   Precision dy2() const { return fTrap.fDy2; }
0128
0129   VECCORE_ATT_HOST_DEVICE
0130   VECGEOM_FORCE_INLINE
0131   Precision dx1() const { return fTrap.fDx1; }
0132
0133   VECCORE_ATT_HOST_DEVICE
0134   VECGEOM_FORCE_INLINE
0135   Precision dx2() const { return fTrap.fDx2; }
0136
0137   VECCORE_ATT_HOST_DEVICE
0138   VECGEOM_FORCE_INLINE
0139   Precision dx3() const { return fTrap.fDx3; }
0140
0141   VECCORE_ATT_HOST_DEVICE
0142   VECGEOM_FORCE_INLINE
0143   Precision dx4() const { return fTrap.fDx4; }
0144
0145   VECCORE_ATT_HOST_DEVICE
0146   VECGEOM_FORCE_INLINE
0147   Precision tanAlpha1() const { return fTrap.fTanAlpha1; }
0148
0149   VECCORE_ATT_HOST_DEVICE
0150   VECGEOM_FORCE_INLINE
0151   Precision tanAlpha2() const { return fTrap.fTanAlpha2; }
0152
0153   VECCORE_ATT_HOST_DEVICE
0154   VECGEOM_FORCE_INLINE
0155   Precision alpha1() const { return GetAlpha1(); } // note: slow, avoid using it
0156
0157   VECCORE_ATT_HOST_DEVICE
0158   VECGEOM_FORCE_INLINE
0159   Precision alpha2() const { return GetAlpha2(); } // note: slow, avoid using it
0160
0161   VECCORE_ATT_HOST_DEVICE
0162   VECGEOM_FORCE_INLINE
0163   Precision tanThetaCosPhi() const { return fTrap.fTthetaCphi; }
0164
0165   VECCORE_ATT_HOST_DEVICE
0166   VECGEOM_FORCE_INLINE
0167   Precision tanThetaSinPhi() const { return fTrap.fTthetaSphi; }
0168
0169   VECCORE_ATT_HOST_DEVICE
0170   VECGEOM_FORCE_INLINE
0171   Precision GetDz() const { return fTrap.fDz; }
0172
0173   VECCORE_ATT_HOST_DEVICE
0174   VECGEOM_FORCE_INLINE
0175   Precision GetTheta() const { return this->theta(); }
0176
0177   VECCORE_ATT_HOST_DEVICE
0178   VECGEOM_FORCE_INLINE
0179   Precision GetPhi() const { return this->phi(); }
0180
0181   VECCORE_ATT_HOST_DEVICE
0182   VECGEOM_FORCE_INLINE
0183   Precision GetDy1() const { return fTrap.fDy1; }
0184
0185   VECCORE_ATT_HOST_DEVICE
0186   VECGEOM_FORCE_INLINE
0187   Precision GetDx1() const { return fTrap.fDx1; }
0188
0189   VECCORE_ATT_HOST_DEVICE
0190   VECGEOM_FORCE_INLINE
0191   Precision GetDx2() const { return fTrap.fDx2; }
0192
0193   VECCORE_ATT_HOST_DEVICE
0194   VECGEOM_FORCE_INLINE
0195   Precision GetTanAlpha1() const { return fTrap.fTanAlpha1; }
0196
0197   VECCORE_ATT_HOST_DEVICE
0198   VECGEOM_FORCE_INLINE
0199   Precision GetDy2() const { return fTrap.fDy2; }
0200
0201   VECCORE_ATT_HOST_DEVICE
0202   VECGEOM_FORCE_INLINE
0203   Precision GetDx3() const { return fTrap.fDx3; }
0204
0205   VECCORE_ATT_HOST_DEVICE
0206   VECGEOM_FORCE_INLINE
0207   Precision GetDx4() const { return fTrap.fDx4; }
0208
0209   VECCORE_ATT_HOST_DEVICE
0210   VECGEOM_FORCE_INLINE
0211   Precision GetTanAlpha2() const { return fTrap.fTanAlpha2; }
0212
0213   VECCORE_ATT_HOST_DEVICE
0214   VECGEOM_FORCE_INLINE
0215   Precision GetTanThetaSinPhi() const { return fTrap.fTthetaSphi; }
0216
0217   VECCORE_ATT_HOST_DEVICE
0218   VECGEOM_FORCE_INLINE
0219   Precision GetTanThetaCosPhi() const { return fTrap.fTthetaCphi; }
0220   /// @}
0221
0222   VECCORE_ATT_HOST_DEVICE
0223   VECGEOM_FORCE_INLINE
0224   void SetDz(Precision val) { fTrap.fDz = val; }
0225
0226   VECCORE_ATT_HOST_DEVICE
0227   VECGEOM_FORCE_INLINE
0228   void SetTheta(Precision val)
0229   {
0230     fTrap.fTheta = val;
0231     fTrap.CalculateCached();
0232     this->MakePlanes();
0233   }
0234
0235   VECCORE_ATT_HOST_DEVICE
0236   VECGEOM_FORCE_INLINE
0237   void SetPhi(Precision val)
0238   {
0239     fTrap.fPhi = val;
0240     fTrap.CalculateCached();
0241     this->MakePlanes();
0242   }
0243
0244   VECCORE_ATT_HOST_DEVICE
0245   VECGEOM_FORCE_INLINE
0246   void SetDy1(Precision val) { fTrap.fDy1 = val; }
0247
0248   VECCORE_ATT_HOST_DEVICE
0249   VECGEOM_FORCE_INLINE
0250   void SetDy2(Precision val) { fTrap.fDy2 = val; }
0251
0252   VECCORE_ATT_HOST_DEVICE
0253   VECGEOM_FORCE_INLINE
0254   void SetDx1(Precision val) { fTrap.fDx1 = val; }
0255
0256   VECCORE_ATT_HOST_DEVICE
0257   VECGEOM_FORCE_INLINE
0258   void SetDx2(Precision val) { fTrap.fDx2 = val; }
0259
0260   VECCORE_ATT_HOST_DEVICE
0261   VECGEOM_FORCE_INLINE
0262   void SetDx3(Precision val) { fTrap.fDx3 = val; }
0263
0264   VECCORE_ATT_HOST_DEVICE
0265   VECGEOM_FORCE_INLINE
0266   void SetDx4(Precision val) { fTrap.fDx4 = val; }
0267
0268   VECCORE_ATT_HOST_DEVICE
0269   VECGEOM_FORCE_INLINE
0270   void SetTanAlpha1(Precision val) { fTrap.fTanAlpha1 = val; }
0271
0272   VECCORE_ATT_HOST_DEVICE
0273   VECGEOM_FORCE_INLINE
0274   void SetTanAlpha2(Precision val) { fTrap.fTanAlpha2 = val; }
0275
0276   VECCORE_ATT_HOST_DEVICE
0277   void Extent(Vector3D<Precision> &, Vector3D<Precision> &) const override;
0278
0279   // VECCORE_ATT_HOST_DEVICE
0280   // void CalcCapacity();
0281
0282   // VECCORE_ATT_HOST_DEVICE
0283   // void CalcSurfaceArea();
0284
0285   // VECCORE_ATT_HOST_DEVICE
0286   Precision Capacity() const override;
0287   // {
0288   //   Assert(!fOutdated);
0289   //   return fCubicVolume;
0290   // }
0291
0292   // VECCORE_ATT_HOST_DEVICE
0293   Precision SurfaceArea() const override;
0294   // {
0295   //   Assert(!fOutdated);
0296   //   return fSurfaceArea;
0297   // }
0298
0299   Vector3D<Precision> SamplePointOnSurface() const override;
0300
0301   Vector3D<Precision> GetPointOnPlane(Vector3D<Precision> const &p0, Vector3D<Precision> const &p1,
0302                                       Vector3D<Precision> const &p2, Vector3D<Precision> const &p3) const;
0303
0304   VECCORE_ATT_HOST_DEVICE
0305   virtual bool Normal(Vector3D<Precision> const &point, Vector3D<Precision> &normal) const override
0306   {
0307     bool valid = false;
0308     normal     = TrapezoidImplementation::NormalKernel(fTrap, point, valid);
0309     return valid;
0310   }
0311
0312   std::string GetEntityType() const { return "Trapezoid"; }
0313
0314   template <typename T>
0315   VECCORE_ATT_HOST_DEVICE void GetParametersList(int aNumber, T *aArray) const;
0316
0317   VECCORE_ATT_HOST_DEVICE
0318   UnplacedTrapezoid *Clone() const;
0319
0320 public:
0321   VECCORE_ATT_HOST_DEVICE
0322   virtual void Print() const override;
0323
0324   virtual void Print(std::ostream &os) const override;
0325
0326 #ifndef VECCORE_CUDA
0327   virtual SolidMesh *CreateMesh3D(Transformation3D const &trans, size_t nSegments) const override;
0328 #endif
0329
0330 #ifdef VECGEOM_CUDA_INTERFACE
0331   virtual size_t DeviceSizeOf() const override { return DevicePtr<cuda::UnplacedTrapezoid>::SizeOf(); }
0332
0333   virtual DevicePtr<cuda::VUnplacedVolume> CopyToGpu() const override;
0334
0335   virtual DevicePtr<cuda::VUnplacedVolume> CopyToGpu(DevicePtr<cuda::VUnplacedVolume> const gpu_ptr) const override;
0336 #endif
0337
0338 #ifndef VECCORE_CUDA
0339   // this is the function called from the VolumeFactory, it may be specific to the trapezoid
0340   static VPlacedVolume *Create(LogicalVolume const *const logical_volume, Transformation3D const *const transformation,
0341                                VPlacedVolume *const placement = NULL);
0342
0343   VPlacedVolume *SpecializedVolume(LogicalVolume const *const volume, Transformation3D const *const transformation,
0344                                    VPlacedVolume *const placement) const override;
0345
0346 #else
0347
0348   VECCORE_ATT_DEVICE
0349   static VPlacedVolume *Create(LogicalVolume const *const logical_volume, Transformation3D const *const transformation,
0350                                const int id, const int copy_no, const int child_id,
0351                                VPlacedVolume *const placement = NULL);
0352
0353   VECCORE_ATT_DEVICE VPlacedVolume *SpecializedVolume(LogicalVolume const *const volume,
0354                                                       Transformation3D const *const transformation, const int id,
0355                                                       const int copy_no, const int child_id,
0356                                                       VPlacedVolume *const placement) const override;
0357 #endif
0358
0359   // Note: use of ATan() makes this one slow -- to be avoided
0360   VECCORE_ATT_HOST_DEVICE
0361   VECGEOM_FORCE_INLINE
0362   Precision GetAlpha1() const { return vecCore::math::ATan(fTrap.fTanAlpha1); }
0363
0364   // Note: use of Atan() makes this one slow -- to be avoided
0365   VECCORE_ATT_HOST_DEVICE
0366   VECGEOM_FORCE_INLINE
0367   Precision GetAlpha2() const { return vecCore::math::ATan(fTrap.fTanAlpha2); }
0368
0369   VECCORE_ATT_HOST_DEVICE
0370   VECGEOM_FORCE_INLINE
0371   virtual ESolidType GetType() const override { return ESolidType::trapezoid; }
0372
0373   // The next functions force upon the user insider knowledge about how the side planes should be used
0374   VECCORE_ATT_HOST_DEVICE
0375   TrapezoidStruct<Precision> const &GetStruct() const { return fTrap; }
0376
0377   // #ifdef VECGEOM_PLANESHELL
0378   //   VECCORE_ATT_HOST_DEVICE
0379   //   VECGEOM_FORCE_INLINE
0380   //   PlaneShell<4,Precision> const *GetPlanes() const { return fTrap.GetPlanes(); }
0381
0382   // #else
0383   //   using TrapSidePlane = TrapezoidStruct<Precision>::TrapSidePlane;
0384   //   VECCORE_ATT_HOST_DEVICE
0385   //   TrapSidePlane const *GetPlanes() const { return fTrap.GetPlanes(); }
0386   // #endif
0387
0388   /// \brief Calculate trapezoid parameters when user provides the 8 corners
0389   VECCORE_ATT_HOST_DEVICE
0390   void fromCornersToParameters(TrapCorners const pt);
0391
0392 private:
0393   /// \brief Calculate the 8 corner points using pre-stored parameters, then use corners to build planes
0394   VECCORE_ATT_HOST_DEVICE
0395   void FromParametersToCorners(TrapCorners pt) const;
0396
0397   // \brief Determine corner points using intersections of the pre-calculated planes
0398   VECCORE_ATT_HOST_DEVICE
0399   void fromPlanesToCorners(TrapCorners pt) const;
0400
0401   /// \brief Construct the four side planes from input corner points
0402   VECCORE_ATT_HOST_DEVICE
0403   bool MakePlanes(TrapCorners const corners);
0404
0405   /// \brief Construct the four side planes by converting stored parameters into TrapCorners object
0406   VECCORE_ATT_HOST_DEVICE
0407   bool MakePlanes();
0408
0409 /// \brief Construct the four side planes from input corner points
0410 #ifdef VECGEOM_PLANESHELL
0411   VECCORE_ATT_HOST_DEVICE
0412   bool MakeAPlane(Vector3D<Precision> const &p1, Vector3D<Precision> const &p2, Vector3D<Precision> const &p3,
0413                   Vector3D<Precision> const &p4, unsigned int iplane);
0414 #else
0415   VECCORE_ATT_HOST_DEVICE
0416   bool MakeAPlane(Vector3D<Precision> const &p1, Vector3D<Precision> const &p2, Vector3D<Precision> const &p3,
0417                   Vector3D<Precision> const &p4, TrapezoidStruct<Precision>::TrapSidePlane &plane);
0418 #endif
0419
0420 public:
0421 #ifndef VECCORE_CUDA
0422 #ifdef VECGEOM_ROOT
0423   TGeoShape const *ConvertToRoot(char const *label) const;
0424 #endif
0425
0426 #ifdef VECGEOM_GEANT4
0427   G4VSolid const *ConvertToGeant4(char const *label) const;
0428 #endif
0429 #endif
0430 };
0431
0432 // Adding specialized factory.
0433 template <>
0434 struct Maker<UnplacedTrapezoid> {
0435   template <typename... ArgTypes>
0436   static UnplacedTrapezoid *MakeInstance(const Precision dz, const Precision theta, const Precision phi,
0437                                          const Precision dy1, const Precision dx1, const Precision dx2,
0438                                          const Precision Alpha1, const Precision dy2, const Precision dx3,
0439                                          const Precision dx4, const Precision Alpha2);
0440
0441   template <typename... ArgTypes>
0442   static UnplacedTrapezoid *MakeInstance(TrapCorners const pt);
0443 };
0444
0445 // Helper function to be used by all the Factories of Trapezoid
0446 #ifndef VECGEOM_NO_SPECIALIZATION
0447 UnplacedTrapezoid *GetSpecialized(const Precision dz, const Precision theta, const Precision phi, const Precision dy1,
0448                                   const Precision dx1, const Precision dx2, const Precision Alpha1, const Precision dy2,
0449                                   const Precision dx3, const Precision dx4, const Precision Alpha2);
0450 #endif
0451
0452 using GenericUnplacedTrapezoid = UnplacedTrapezoid;
0453
0454 } // namespace VECGEOM_IMPL_NAMESPACE
0455 } // namespace vecgeom
0456
0457 #endif