surfaces/conv/CutTubeConverter.h

0001 #ifndef VECGEOM_SURFACE_CUTTUBECONVERTER_H_
0002 #define VECGEOM_SURFACE_CUTTUBECONVERTER_H_
0003
0004 #include <VecGeom/surfaces/conv/Builder.h>
0005 #include <VecGeom/surfaces/Model.h>
0006
0007 #include <VecGeom/volumes/CutTube.h>
0008
0009 namespace vgbrep {
0010 namespace conv {
0011
0012 // sign function
0013 template <typename Real_t>
0014 int sgn(Real_t val)
0015 {
0016   return (Real_t(0) < val) - (val < Real_t(0));
0017 }
0018
0019 /// @brief Converter for Tube
0020 /// @tparam Real_t Precision type
0021 /// @param tube Tube solid to be converted
0022 /// @param logical_id Id of the logical volume
0023 /// @return Conversion success
0024 template <typename Real_t>
0025 bool CreateTubeSurfaces(vecgeom::UnplacedCutTube const &tube, int logical_id, bool intersection = false)
0026 {
0027   using ZPhiMask_t   = ZPhiMask<Real_t>;
0028   using WindowMask_t = WindowMask<Real_t>;
0029   using Vector3      = vecgeom::Vector3D<vecgeom::Precision>;
0030
0031   LogicExpressionCPU logic; // top & bottom & [rmin] & rmax & (dphi < 180) ? sphi * ephi : sphi | ephi
0032   auto sphi  = tube.sphi();
0033   auto dphi  = tube.dphi();
0034   auto ephi  = tube.sphi() + tube.dphi();
0035   auto csphi = vecCore::math::Cos(sphi);
0036   auto ssphi = vecCore::math::Sin(sphi);
0037   auto cephi = vecCore::math::Cos(ephi);
0038   auto sephi = vecCore::math::Sin(ephi);
0039   auto rmin  = tube.rmin();
0040   auto rmax  = tube.rmax();
0041
0042   // get normal of the end caps
0043   auto bottom_normal = tube.BottomNormal();
0044   auto top_normal    = tube.TopNormal();
0045
0046   // get maximum extent
0047   Vector3 aMin, aMax;
0048   tube.Extent(aMin, aMax);
0049
0050   VECGEOM_ASSERT(dphi > vecgeom::kTolerance);
0051
0052   bool fullCirc  = ApproxEqual(dphi, vecgeom::kTwoPi);
0053   bool smallerPi = dphi < (vecgeom::kPi - vecgeom::kTolerance);
0054
0055   int isurf;
0056   vecgeom::Precision surfdata[2];
0057
0058   // We need angles in degrees for transformations
0059   auto thetad_top    = top_normal.Theta() * vecgeom::kRadToDeg;
0060   auto phid_top      = top_normal.Phi() * vecgeom::kRadToDeg;
0061   auto thetad_bottom = bottom_normal.Theta() * vecgeom::kRadToDeg;
0062   auto phid_bottom   = bottom_normal.Phi() * vecgeom::kRadToDeg;
0063   // assert for 0 <= theta top <= 90 and 90 <= theta bottom <= 180
0064   VECGEOM_ASSERT(top_normal[2] >= 0 && bottom_normal[2] <= 0);
0065
0066   auto &cpudata = CPUsurfData<Real_t>::Instance();
0067
0068   // surface at +dz
0069   // due to the cut the top and bottom caps are elliptical, which is handled by using a rectangle frame
0070   // that fully contains the ellipse and making the surface logical such that the full boolean expression
0071   // must be evaluated to decide whether it is a hit or not.
0072   // As a consequence, all surfaces of this volume must be logical surfaces.
0073   VECGEOM_ASSERT(std::abs(top_normal.z()) > vecgeom::kTolerance); // assert before division
0074   isurf = builder::CreateLocalSurface<Real_t>(
0075       builder::CreateUnplacedSurface<Real_t>(SurfaceType::kPlanar),
0076       builder::CreateFrame<Real_t>(FrameType::kWindow, WindowMask_t{tube.rmax(), tube.rmax() / top_normal.z()}),
0077       vecgeom::Transformation3DMP<Precision>(0., 0., tube.z(), phid_top - 90, -thetad_top, 0.));
0078   auto &surf = cpudata.fLocalSurfaces[isurf];
0079   builder::AddSurfaceToShell<Real_t>(logical_id, isurf);
0080   if (intersection) surf.fSkipConvexity = true;
0081   // Make the surface "logical"
0082   surf.fLogicId = isurf;
0083   logic.push_back(isurf);
0084   logic.push_back(land);
0085
0086   // surface at -dz
0087   VECGEOM_ASSERT(std::abs(std::cos(vecgeom::kPi - bottom_normal.Theta())) >
0088                  vecgeom::kTolerance); // assert before division
0089   isurf = builder::CreateLocalSurface<Real_t>(
0090       builder::CreateUnplacedSurface<Real_t>(SurfaceType::kPlanar),
0091       builder::CreateFrame<Real_t>(
0092           FrameType::kWindow,
0093           WindowMask_t{
0094               tube.rmax(),
0095               static_cast<vecgeom::Precision>(
0096                   tube.rmax() / cos(vecgeom::kPi - bottom_normal.Theta()))}), // static_cast since cos returns double
0097       vecgeom::Transformation3DMP<Precision>(0., 0., -tube.z(), phid_bottom - 90, -thetad_bottom, 0.));
0098   auto &surf2 = cpudata.fLocalSurfaces[isurf];
0099   builder::AddSurfaceToShell<Real_t>(logical_id, isurf);
0100   if (intersection) surf2.fSkipConvexity = true;
0101   // Make the surface "logical"
0102   surf2.fLogicId = isurf;
0103   logic.push_back(isurf);
0104
0105   vecgeom::Transformation3DMP<Real_t> identity;
0106
0107   // inner cylinder
0108   if (tube.rmin() > vecgeom::kTolerance) {
0109     surfdata[0] = tube.rmin();
0110     isurf       = builder::CreateLocalSurface<Real_t>(
0111         builder::CreateUnplacedSurface<Real_t>(SurfaceType::kCylindrical, surfdata, /*flipped=*/true),
0112         builder::CreateFrame<Real_t>(FrameType::kZPhi,
0113                                      ZPhiMask_t{aMin[2], aMax[2], fullCirc, tube.rmin(), tube.rmin(), sphi, ephi}),
0114         /*identity transformation*/ identity);
0115     auto &surf3 = cpudata.fLocalSurfaces[isurf];
0116     if (intersection) surf3.fSkipConvexity = true;
0117     builder::AddSurfaceToShell<Real_t>(logical_id, isurf);
0118     // Make the surface "logical"
0119     surf3.fLogicId = isurf;
0120     logic.push_back(land);
0121     logic.push_back(isurf);
0122   }
0123   // outer cylinder
0124   surfdata[0] = tube.rmax();
0125   isurf       = builder::CreateLocalSurface<Real_t>(
0126       builder::CreateUnplacedSurface<Real_t>(SurfaceType::kCylindrical, surfdata),
0127       builder::CreateFrame<Real_t>(FrameType::kZPhi,
0128                                    ZPhiMask_t{aMin[2], aMax[2], fullCirc, tube.rmax(), tube.rmax(), sphi, ephi}),
0129       /*identity transformation*/ identity);
0130   auto &surf4 = cpudata.fLocalSurfaces[isurf];
0131   if (intersection) surf4.fSkipConvexity = true;
0132   builder::AddSurfaceToShell<Real_t>(logical_id, isurf);
0133   // Make the surface "logical"
0134   surf4.fLogicId = isurf;
0135   logic.push_back(land);
0136   logic.push_back(isurf);
0137
0138   if (ApproxEqual(dphi, vecgeom::kTwoPi)) {
0139     builder::AddLogicToShell<Real_t>(logical_id, logic);
0140     return true;
0141   }
0142
0143   // Implementation of the plane caps at Sphi and Ephi by using a window with the size of the extent of the
0144   // quadrilaterals. Depending on the phi cut position, this window is a bit larger than the actual surface, but we
0145   // found this implementation to be faster than a quadrilateral with the correct size plane cap at Sphi
0146   auto zmin1 = std::min(-tube.z() - rmin * (bottom_normal.x() * csphi + bottom_normal.y() * ssphi) / bottom_normal.z(),
0147                         -tube.z() - rmax * (bottom_normal.x() * csphi + bottom_normal.y() * ssphi) / bottom_normal.z());
0148   auto zmax1 = std::max(tube.z() - rmin * (top_normal.x() * csphi + top_normal.y() * ssphi) / top_normal.z(),
0149                         tube.z() - rmax * (top_normal.x() * csphi + top_normal.y() * ssphi) / top_normal.z());
0150   auto zmin2 = std::min(-tube.z() - rmin * (bottom_normal.x() * cephi + bottom_normal.y() * sephi) / bottom_normal.z(),
0151                         -tube.z() - rmax * (bottom_normal.x() * cephi + bottom_normal.y() * sephi) / bottom_normal.z());
0152   auto zmax2 = std::max(tube.z() - rmin * (top_normal.x() * cephi + top_normal.y() * sephi) / top_normal.z(),
0153                         tube.z() - rmax * (top_normal.x() * cephi + top_normal.y() * sephi) / top_normal.z());
0154
0155   std::vector<Vector3> vert(4);
0156   vert[0].Set(rmin * csphi, rmin * ssphi, zmin1);
0157   vert[1].Set(rmax * csphi, rmax * ssphi, zmin1);
0158   vert[2].Set(rmax * csphi, rmax * ssphi, zmax1);
0159   vert[3].Set(rmin * csphi, rmin * ssphi, zmax1);
0160   isurf = builder::CreateLocalSurfaceFromVertices<Real_t>(vert, logical_id);
0161   if (intersection) builder::GetSurface<Real_t>(isurf).fSkipConvexity = true;
0162   VECGEOM_ASSERT(isurf >= 0);
0163   if (!smallerPi) builder::GetSurface<Real_t>(isurf).fEmbedding = false;
0164   // Make the surface "logical"
0165   cpudata.fLocalSurfaces[isurf].fLogicId = isurf;
0166   logic.push_back(land);
0167   logic.push_back(lplus); // '('
0168   logic.push_back(isurf);
0169
0170   // plane cap at Sphi+Dphi
0171   vert[0].Set(rmax * cephi, rmax * sephi, zmin2);
0172   vert[1].Set(rmin * cephi, rmin * sephi, zmin2);
0173   vert[2].Set(rmin * cephi, rmin * sephi, zmax2);
0174   vert[3].Set(rmax * cephi, rmax * sephi, zmax2);
0175   isurf = builder::CreateLocalSurfaceFromVertices<Real_t>(vert, logical_id);
0176   if (intersection) builder::GetSurface<Real_t>(isurf).fSkipConvexity = true;
0177   VECGEOM_ASSERT(isurf >= 0);
0178   if (!smallerPi) builder::GetSurface<Real_t>(isurf).fEmbedding = false;
0179   // Make the surface "logical"
0180   cpudata.fLocalSurfaces[isurf].fLogicId = isurf;
0181   logic.push_back(smallerPi ? land : lor);
0182   logic.push_back(isurf);
0183   logic.push_back(lminus); // ')'
0184   builder::AddLogicToShell<Real_t>(logical_id, logic);
0185   return true;
0186 }
0187
0188 } // namespace conv
0189 } // namespace vgbrep
0190 #endif