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 SIMDUnplacedVolumeImplHelper<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
0316   void GetParametersList(int aNumber, T *aArray) const;
0317
0318   VECCORE_ATT_HOST_DEVICE
0319   UnplacedTrapezoid *Clone() const;
0320
0321 public:
0322   VECCORE_ATT_HOST_DEVICE
0323   virtual void Print() const override;
0324
0325   virtual void Print(std::ostream &os) const override;
0326
0327 #ifndef VECCORE_CUDA
0328   virtual SolidMesh *CreateMesh3D(Transformation3D const &trans, size_t nSegments) const override;
0329 #endif
0330
0331 #ifdef VECGEOM_CUDA_INTERFACE
0332   virtual size_t DeviceSizeOf() const override { return DevicePtr<cuda::UnplacedTrapezoid>::SizeOf(); }
0333
0334   virtual DevicePtr<cuda::VUnplacedVolume> CopyToGpu() const override;
0335
0336   virtual DevicePtr<cuda::VUnplacedVolume> CopyToGpu(DevicePtr<cuda::VUnplacedVolume> const gpu_ptr) const override;
0337 #endif
0338
0339 #ifndef VECCORE_CUDA
0340   // this is the function called from the VolumeFactory, it may be specific to the trapezoid
0341   template <TranslationCode trans_code, RotationCode rot_code>
0342   static VPlacedVolume *Create(LogicalVolume const *const logical_volume, Transformation3D const *const transformation,
0343                                VPlacedVolume *const placement = NULL);
0344
0345   VPlacedVolume *SpecializedVolume(LogicalVolume const *const volume, Transformation3D const *const transformation,
0346                                    const TranslationCode trans_code, const RotationCode rot_code,
0347                                    VPlacedVolume *const placement) const override;
0348
0349 #else
0350
0351   template <TranslationCode trans_code, RotationCode rot_code>
0352   VECCORE_ATT_DEVICE
0353   static VPlacedVolume *Create(LogicalVolume const *const logical_volume, Transformation3D const *const transformation,
0354                                const int id, const int copy_no, const int child_id,
0355                                VPlacedVolume *const placement = NULL);
0356
0357   VECCORE_ATT_DEVICE VPlacedVolume *SpecializedVolume(LogicalVolume const *const volume,
0358                                                       Transformation3D const *const transformation,
0359                                                       const TranslationCode trans_code, const RotationCode rot_code,
0360                                                       const int id, const int copy_no, const int child_id,
0361                                                       VPlacedVolume *const placement) const override;
0362 #endif
0363
0364   // Note: use of ATan() makes this one slow -- to be avoided
0365   VECCORE_ATT_HOST_DEVICE
0366   VECGEOM_FORCE_INLINE
0367   Precision GetAlpha1() const { return vecCore::math::ATan(fTrap.fTanAlpha1); }
0368
0369   // Note: use of Atan() makes this one slow -- to be avoided
0370   VECCORE_ATT_HOST_DEVICE
0371   VECGEOM_FORCE_INLINE
0372   Precision GetAlpha2() const { return vecCore::math::ATan(fTrap.fTanAlpha2); }
0373
0374   // The next functions force upon the user insider knowledge about how the side planes should be used
0375   VECCORE_ATT_HOST_DEVICE
0376   TrapezoidStruct<Precision> const &GetStruct() const { return fTrap; }
0377
0378   // #ifndef VECGEOM_PLANESHELL_DISABLE
0379   //   VECCORE_ATT_HOST_DEVICE
0380   //   VECGEOM_FORCE_INLINE
0381   //   PlaneShell<4,Precision> const *GetPlanes() const { return fTrap.GetPlanes(); }
0382
0383   // #else
0384   //   using TrapSidePlane = TrapezoidStruct<Precision>::TrapSidePlane;
0385   //   VECCORE_ATT_HOST_DEVICE
0386   //   TrapSidePlane const *GetPlanes() const { return fTrap.GetPlanes(); }
0387   // #endif
0388
0389   /// \brief Calculate trapezoid parameters when user provides the 8 corners
0390   VECCORE_ATT_HOST_DEVICE
0391   void fromCornersToParameters(TrapCorners const pt);
0392
0393 private:
0394   /// \brief Calculate the 8 corner points using pre-stored parameters, then use corners to build planes
0395   VECCORE_ATT_HOST_DEVICE
0396   void FromParametersToCorners(TrapCorners pt) const;
0397
0398   // \brief Determine corner points using intersections of the pre-calculated planes
0399   VECCORE_ATT_HOST_DEVICE
0400   void fromPlanesToCorners(TrapCorners pt) const;
0401
0402   /// \brief Construct the four side planes from input corner points
0403   VECCORE_ATT_HOST_DEVICE
0404   bool MakePlanes(TrapCorners const corners);
0405
0406   /// \brief Construct the four side planes by converting stored parameters into TrapCorners object
0407   VECCORE_ATT_HOST_DEVICE
0408   bool MakePlanes();
0409
0410 /// \brief Construct the four side planes from input corner points
0411 #ifndef VECGEOM_PLANESHELL_DISABLE
0412   VECCORE_ATT_HOST_DEVICE
0413   bool MakeAPlane(Vector3D<Precision> const &p1, Vector3D<Precision> const &p2, Vector3D<Precision> const &p3,
0414                   Vector3D<Precision> const &p4, unsigned int iplane);
0415 #else
0416   VECCORE_ATT_HOST_DEVICE
0417   bool MakeAPlane(Vector3D<Precision> const &p1, Vector3D<Precision> const &p2, Vector3D<Precision> const &p3,
0418                   Vector3D<Precision> const &p4, TrapezoidStruct<Precision>::TrapSidePlane &plane);
0419 #endif
0420
0421 public:
0422 #ifndef VECCORE_CUDA
0423 #ifdef VECGEOM_ROOT
0424   TGeoShape const *ConvertToRoot(char const *label) const;
0425 #endif
0426
0427 #ifdef VECGEOM_GEANT4
0428   G4VSolid const *ConvertToGeant4(char const *label) const;
0429 #endif
0430 #endif
0431 };
0432
0433 // Adding specialized factory.
0434 template <>
0435 struct Maker<UnplacedTrapezoid> {
0436   template <typename... ArgTypes>
0437   static UnplacedTrapezoid *MakeInstance(const Precision dz, const Precision theta, const Precision phi,
0438                                          const Precision dy1, const Precision dx1, const Precision dx2,
0439                                          const Precision Alpha1, const Precision dy2, const Precision dx3,
0440                                          const Precision dx4, const Precision Alpha2);
0441
0442   template <typename... ArgTypes>
0443   static UnplacedTrapezoid *MakeInstance(TrapCorners const pt);
0444 };
0445
0446 // Helper function to be used by all the Factories of Trapezoid
0447 #ifndef VECGEOM_NO_SPECIALIZATION
0448 UnplacedTrapezoid *GetSpecialized(const Precision dz, const Precision theta, const Precision phi, const Precision dy1,
0449                                   const Precision dx1, const Precision dx2, const Precision Alpha1, const Precision dy2,
0450                                   const Precision dx3, const Precision dx4, const Precision Alpha2);
0451 #endif
0452
0453 using GenericUnplacedTrapezoid = UnplacedTrapezoid;
0454
0455 } // namespace VECGEOM_IMPL_NAMESPACE
0456 } // namespace vecgeom
0457
0458 #endif