surfaces/conv/ConeConverter.h

0001 #ifndef VECGEOM_SURFACE_CONECONVERTER_H_
0002 #define VECGEOM_SURFACE_CONECONVERTER_H_
0003
0004 #include <VecGeom/surfaces/conv/Builder.h>
0005 #include <VecGeom/surfaces/Model.h>
0006
0007 #include <VecGeom/volumes/Cone.h>
0008
0009 namespace vgbrep {
0010 namespace conv {
0011
0012 /// @brief Converter for Cone
0013 /// @tparam Real_t Precision type
0014 /// @param cone Cone solid to be converted
0015 /// @param logical_id Id of the logical volume
0016 /// @return Conversion success
0017 template <typename Real_t>
0018 bool CreateConeSurfaces(vecgeom::UnplacedCone const &cone, int logical_id, bool intersection = false)
0019 {
0020   using RingMask_t = RingMask<Real_t>;
0021   using ZPhiMask_t = ZPhiMask<Real_t>;
0022   using Vector3D   = vecgeom::Vector3D<vecgeom::Precision>;
0023
0024   LogicExpressionCPU logic; // top & bottom & [rmin] & rmax & (dphi < 180) ? sphi * ephi : sphi | ephi
0025   auto rmin1 = cone.GetRmin1();
0026   auto rmax1 = cone.GetRmax1();
0027   auto rmin2 = cone.GetRmin2();
0028   auto rmax2 = cone.GetRmax2();
0029   auto dz    = cone.GetDz();
0030   auto sphi  = cone.GetSPhi();
0031   auto dphi  = cone.GetDPhi();
0032   auto ephi  = dphi + sphi;
0033
0034   VECGEOM_ASSERT(dphi > vecgeom::kTolerance);
0035   VECGEOM_ASSERT(rmax1 - rmin1 > -vecgeom::kTolerance);
0036   VECGEOM_ASSERT(rmax2 - rmin2 > -vecgeom::kTolerance);
0037
0038   // if rmax == rmin then for safety it is set to rmax = rmin + ConeTolerance already in the solid model
0039   // For using mixed precision, the compiled tolerance must be replaced by the mixed tolerance
0040   if (rmax1 - rmin1 == vecgeom::kConeTolerance) rmax1 = rmin1 + vecgeom::kToleranceCone<Real_t>;
0041   if (rmax2 - rmin2 == vecgeom::kConeTolerance) rmax2 = rmin2 + vecgeom::kToleranceCone<Real_t>;
0042
0043   bool fullCirc  = ApproxEqual(dphi, vecgeom::kTwoPi);
0044   bool smallerPi = dphi < (vecgeom::kPi - vecgeom::kTolerance);
0045
0046   int isurf;
0047   vecgeom::Precision surfdata[2];
0048
0049   // We need angles in degrees for transformations
0050   auto sphid = vecgeom::kRadToDeg * sphi;
0051   auto ephid = vecgeom::kRadToDeg * ephi;
0052   vecgeom::Transformation3DMP<Real_t> identity;
0053
0054   // surface at +Dz
0055   isurf = builder::CreateLocalSurface<Real_t>(
0056       builder::CreateUnplacedSurface<Real_t>(SurfaceType::kPlanar),
0057       builder::CreateFrame<Real_t>(FrameType::kRing, RingMask_t{rmin2, rmax2, fullCirc, sphi, ephi}),
0058       vecgeom::Transformation3DMP<Precision>(0., 0., dz, 0., 0., 0.));
0059   builder::AddSurfaceToShell<Real_t>(logical_id, isurf);
0060   if (intersection) builder::GetSurface<Real_t>(isurf).fSkipConvexity = true;
0061   logic.push_back(isurf);
0062   // surface at -Dz
0063   isurf = builder::CreateLocalSurface<Real_t>(
0064       builder::CreateUnplacedSurface<Real_t>(SurfaceType::kPlanar),
0065       builder::CreateFrame<Real_t>(FrameType::kRing, RingMask_t{rmin1, rmax1, fullCirc, sphi, ephi}),
0066       vecgeom::Transformation3DMP<Precision>(0., 0., -dz, 0., 180., -sphid - ephid));
0067   builder::AddSurfaceToShell<Real_t>(logical_id, isurf);
0068   if (intersection) builder::GetSurface<Real_t>(isurf).fSkipConvexity = true;
0069   logic.push_back(land);
0070   logic.push_back(isurf);
0071   // inner cone
0072   if (rmin2 > vecgeom::kTolerance || rmin1 > vecgeom::kTolerance) {
0073     surfdata[0] = 0.5 * (rmin1 + rmin2);
0074     surfdata[1] = 0.5 * (rmin2 - rmin1) / dz;
0075     auto stype  = std::abs(surfdata[1]) > vecgeom::kTolerance ? SurfaceType::kConical : SurfaceType::kCylindrical;
0076     isurf       = builder::CreateLocalSurface<Real_t>(
0077         builder::CreateUnplacedSurface<Real_t>(stype, surfdata, /*flipped=*/true),
0078         builder::CreateFrame<Real_t>(FrameType::kZPhi, ZPhiMask_t{-dz, dz, fullCirc, rmin1, rmin2, sphi, ephi}),
0079         /*identity transformation*/ identity);
0080     builder::AddSurfaceToShell<Real_t>(logical_id, isurf);
0081     if (intersection) builder::GetSurface<Real_t>(isurf).fSkipConvexity = true;
0082     logic.push_back(land);
0083     logic.push_back(isurf);
0084   }
0085   // outer cone
0086   surfdata[0] = 0.5 * (rmax1 + rmax2);
0087   surfdata[1] = 0.5 * (rmax2 - rmax1) / dz;
0088   auto stype  = std::abs(surfdata[1]) > vecgeom::kTolerance ? SurfaceType::kConical : SurfaceType::kCylindrical;
0089   isurf       = builder::CreateLocalSurface<Real_t>(
0090       builder::CreateUnplacedSurface<Real_t>(stype, surfdata),
0091       builder::CreateFrame<Real_t>(FrameType::kZPhi, ZPhiMask_t{-dz, dz, fullCirc, rmax1, rmax2, sphi, ephi}),
0092       /*identity transformation*/ identity);
0093   builder::AddSurfaceToShell<Real_t>(logical_id, isurf);
0094   if (intersection) builder::GetSurface<Real_t>(isurf).fSkipConvexity = true;
0095   logic.push_back(land);
0096   logic.push_back(isurf);
0097
0098   if (ApproxEqual(dphi, vecgeom::kTwoPi)) {
0099     builder::AddLogicToShell<Real_t>(logical_id, logic);
0100     return true;
0101   }
0102
0103   vecgeom::Transformation3D transformation;
0104   std::vector<Vector3D> vert; // Stores 3D coordinates of the four corners that create a quadrilateral.
0105   auto csphi = std::cos(sphi);
0106   auto ssphi = std::sin(sphi);
0107   auto cephi = std::cos(ephi);
0108   auto sephi = std::sin(ephi);
0109
0110   // corners of the sphi and ephi faces
0111   std::vector<Vector3D> corners = {{rmin1 * csphi, rmin1 * ssphi, -dz}, {rmax1 * csphi, rmax1 * ssphi, -dz},
0112                                    {rmax2 * csphi, rmax2 * ssphi, dz},  {rmin2 * csphi, rmin2 * ssphi, dz},
0113                                    {rmax1 * cephi, rmax1 * sephi, -dz}, {rmin1 * cephi, rmin1 * sephi, -dz},
0114                                    {rmin2 * cephi, rmin2 * sephi, dz},  {rmax2 * cephi, rmax2 * sephi, dz}};
0115
0116   // plane cap at sphi
0117   vert  = {corners[0], corners[1], corners[2], corners[3]};
0118   isurf = builder::CreateLocalSurfaceFromVertices<Real_t>(vert, logical_id);
0119   if (intersection) builder::GetSurface<Real_t>(isurf).fSkipConvexity = true;
0120   if (isurf >= 0) {
0121     if (!smallerPi) builder::GetSurface<Real_t>(isurf).fEmbedding = false;
0122     logic.push_back(land);
0123     logic.push_back(lplus); // '('
0124     logic.push_back(isurf);
0125   }
0126
0127   // plane cap at sphi+dphi
0128   vert  = {corners[4], corners[5], corners[6], corners[7]};
0129   isurf = builder::CreateLocalSurfaceFromVertices<Real_t>(vert, logical_id);
0130   if (intersection) builder::GetSurface<Real_t>(isurf).fSkipConvexity = true;
0131   if (isurf >= 0) {
0132     if (!smallerPi) builder::GetSurface<Real_t>(isurf).fEmbedding = false;
0133     logic.push_back(smallerPi ? land : lor);
0134     logic.push_back(isurf);
0135     logic.push_back(lminus); // ')'
0136   }
0137
0138   builder::AddLogicToShell<Real_t>(logical_id, logic);
0139   return true;
0140 }
0141
0142 } // namespace conv
0143 } // namespace vgbrep
0144 #endif