Back to home page

EIC code displayed by LXR
 
 

    

File indexing completed on 2026-04-17 08:35:36 

0001 #ifndef VECGEOM_TORUS_IMPL_H
0002 #define VECGEOM_TORUS_IMPL_H
0003 
0004 #include <VecGeom/surfaces/surf/SurfaceHelper.h>
0005 #include <VecGeom/surfaces/base/Equations.h>
0006 
0007 namespace vgbrep {
0008 
0009 template <typename Real_t>
0010 struct SurfaceHelper<SurfaceType::kTorus, Real_t> {
0011   TorusData<Real_t> const *fTorusData{nullptr};
0012 
0013   VECGEOM_FORCE_INLINE
0014   VECCORE_ATT_HOST_DEVICE
0015   SurfaceHelper(TorusData<Real_t> const &torusdata) { fTorusData = &torusdata; }
0016 
0017   /// @brief Fills the 3D extent of the Arb4
0018   /// @param aMin Bottom extent corner
0019   /// @param aMax Top extent corner
0020   void Extent3D(Vector3D<Real_t> &aMin, Vector3D<Real_t> &aMax) const
0021   {
0022     auto rtor  = fTorusData->Radius();
0023     auto rtube = fTorusData->RadiusTube();
0024     aMin.Set(-rtor - rtube, -rtor - rtube, -rtube);
0025     aMax.Set(rtor + rtube, rtor + rtube, rtube);
0026   }
0027 
0028   VECGEOM_FORCE_INLINE
0029   VECCORE_ATT_HOST_DEVICE
0030   /// @brief Inside half-space function
0031   /// @param point Point in local surface coordinates
0032   /// @param tol tolerance for determining if a point is inside or not
0033   /// @return True if the point is behind the normal within kTolerance (surface is included)
0034   bool Inside(Vector3D<Real_t> const &point, bool flip = false, Real_t tol = vecgeom::kToleranceStrict<Real_t>)
0035   {
0036     if (!fTorusData->InsidePhi(point, flip)) return false;
0037 
0038     int flipsign = (flip ^ fTorusData->IsFlipped()) ? -1 : 1;
0039     bool flipped = fTorusData->IsFlipped() ? 1 : 0;
0040     Real_t rho   = Sqrt(point.x() * point.x() + point.y() * point.y());
0041     Real_t rTor  = fTorusData->Radius();
0042     Real_t rTube = fTorusData->RadiusTube();
0043 
0044     bool check1 = (rho - (rTor - rTube)) > -flipsign * tol;
0045     bool check2 = (rho - (rTor + rTube)) < flipsign * tol;
0046     bool check3 = Sqrt(point.z() * point.z() + (rho - rTor) * (rho - rTor)) - rTube < flipsign * tol;
0047 
0048     if (!flipped) {
0049       return check1 && check2 && check3;
0050     } else {
0051       return !check1 || !check2 || !check3;
0052     }
0053   }
0054 
0055   VECGEOM_FORCE_INLINE
0056   VECCORE_ATT_HOST_DEVICE
0057   /// @brief Find signed distance to next intersection from local point.
0058   /// @param point Point in local surface coordinates
0059   /// @param dir Direction in the local surface coordinates
0060   /// @param left_side Flag specifying if the surface is intersected from the left-side that defines the normal
0061   /// @param distance Computed distance to surface
0062   /// @return Validity of the intersection
0063   bool Intersect(Vector3D<Real_t> const &point, Vector3D<Real_t> const &dir, bool left_side, Real_t &distance,
0064                  bool &two_solutions, Real_t &safety)
0065   {
0066 
0067     // RESCALING, from now on RadTube is normalized to TorusRadius and TorusRadius = 1!
0068     Vector3D<Real_t> localpoint = point / fTorusData->Radius();
0069     Real_t RadTube_R0           = fTorusData->RadiusTube() / fTorusData->Radius();
0070 
0071     Real_t tubeDistance = 0;
0072 
0073     // if point is outside bounding tube of the torus, propagate to the bounding tube.
0074     if (!(SurfaceHelper<SurfaceType::kCylindrical, Real_t>().Inside(localpoint, false, fTorusData->InnerBCRadius()) &&
0075           SurfaceHelper<SurfaceType::kCylindrical, Real_t>().Inside(localpoint, false, fTorusData->OuterBCRadius())) ||
0076         (Abs(localpoint[2]) > Abs(RadTube_R0 + vecgeom::kTolerance))) {
0077 
0078       Real_t tmp   = vecgeom::InfinityLength<Real_t>();
0079       tubeDistance = -1; // ensure it returns if no hit
0080       // check upper plane
0081       // emulate transformation of upper surface
0082       localpoint[2] = point[2] / fTorusData->Radius() - RadTube_R0;
0083       Vector3D<Real_t> tmppoint;
0084       Real_t tmp_rho;
0085       if (SurfaceHelper<SurfaceType::kPlanar, Real_t>().Intersect(localpoint, dir, false, tmp, two_solutions, safety)) {
0086         tmppoint = point / fTorusData->Radius() + tmp * dir;
0087         tmp_rho  = Sqrt(tmppoint[0] * tmppoint[0] + tmppoint[1] * tmppoint[1]);
0088         if ((tmp > -vecgeom::kTolerance) && (tmp_rho - (1. - RadTube_R0) > -vecgeom::kTolerance) &&
0089             (tmp_rho - (1. + RadTube_R0) < vecgeom::kTolerance)) {
0090           tubeDistance = tmp;
0091         }
0092       }
0093 
0094       // check lower plane
0095       // emulate transformation of lower surface
0096       Vector3D<Real_t> localdir = {dir[0], dir[1], -dir[2]};
0097       localpoint[2]             = -point[2] / fTorusData->Radius() - RadTube_R0;
0098       if (SurfaceHelper<SurfaceType::kPlanar, Real_t>().Intersect(localpoint, localdir, false, tmp, two_solutions,
0099                                                                   safety)) {
0100         tmppoint = point / fTorusData->Radius() + tmp * dir;
0101         tmp_rho  = Sqrt(tmppoint[0] * tmppoint[0] + tmppoint[1] * tmppoint[1]);
0102         if ((tmp > -vecgeom::kTolerance) && (tmp_rho - (1. - RadTube_R0) > -vecgeom::kTolerance) &&
0103             (tmp_rho - (1. + RadTube_R0) < vecgeom::kTolerance) && (tmp > -vecgeom::kTolerance)) {
0104           if (tubeDistance > -vecgeom::kTolerance) {
0105             tubeDistance = Min(tubeDistance, tmp);
0106           } else {
0107             tubeDistance = tmp;
0108           }
0109         }
0110       }
0111 
0112       // check outer cylinder
0113       localpoint = point / fTorusData->Radius();
0114       if (SurfaceHelper<SurfaceType::kCylindrical, Real_t>().Intersect(localpoint, dir, false, tmp, two_solutions,
0115                                                                        safety, fTorusData->OuterBCRadius())) {
0116         tmppoint = point / fTorusData->Radius() + tmp * dir;
0117         if (tmp > -vecgeom::kTolerance && (Abs(tmppoint[2]) < Abs(RadTube_R0 + vecgeom::kTolerance))) {
0118           if (tubeDistance > -vecgeom::kTolerance) {
0119             tubeDistance = Min(tubeDistance, tmp);
0120           } else {
0121             tubeDistance = tmp;
0122           }
0123         }
0124       }
0125 
0126       // check inner cylinder
0127       localpoint = point / fTorusData->Radius();
0128       if (SurfaceHelper<SurfaceType::kCylindrical, Real_t>().Intersect(localpoint, dir, false, tmp, two_solutions,
0129                                                                        safety, fTorusData->InnerBCRadius())) {
0130         tmppoint = point / fTorusData->Radius() + tmp * dir;
0131         if (tmp > -vecgeom::kTolerance && (Abs(tmppoint[2]) < Abs(RadTube_R0 + vecgeom::kTolerance))) {
0132           if (tubeDistance > -vecgeom::kTolerance) {
0133             tubeDistance = Min(tubeDistance, tmp);
0134           } else {
0135             tubeDistance = tmp;
0136           }
0137         }
0138       }
0139 
0140       // bounding tube not hit, cannot hit torus either, return
0141       if (tubeDistance < -vecgeom::kTolerance) return false;
0142 
0143       // transform to hit point on bounding tube
0144       localpoint = point / fTorusData->Radius() + tubeDistance * dir;
0145     }
0146 
0147     QuarticCoef<Real_t> coef;
0148     Real_t roots[4]          = {vecgeom::InfinityLength<Real_t>(), vecgeom::InfinityLength<Real_t>(),
0149                                 vecgeom::InfinityLength<Real_t>(), vecgeom::InfinityLength<Real_t>()};
0150     int numroots             = 0;
0151     Real_t s                 = vecgeom::InfinityLength<Real_t>();
0152     VECGEOM_CONST Real_t tol = 100. * vecgeom::kTolerance;
0153 
0154     bool flip_exiting = left_side ^ fTorusData->IsFlipped();
0155     TorusEq<Real_t>(localpoint, dir, 1, RadTube_R0, coef);
0156 
0157     // special condition
0158     if (Abs(dir[2]) < 1E-3 && Abs(localpoint[2]) < 0.1 * RadTube_R0) {
0159       Real_t r0        = 1 - Sqrt((RadTube_R0 - localpoint[2]) * (RadTube_R0 + localpoint[2]));
0160       Real_t invdirxy2 = 1. / (1 - dir.z() * dir.z());
0161       Real_t b0        = (localpoint[0] * dir[0] + localpoint[1] * dir[1]) * invdirxy2;
0162       Real_t c0        = (localpoint[0] * localpoint[0] + (localpoint[1] - r0) * (localpoint[1] + r0)) * invdirxy2;
0163       Real_t delta     = b0 * b0 - c0;
0164       if (delta > 0) {
0165         roots[numroots] = -b0 - Sqrt(delta);
0166         if (roots[numroots] > -tol) numroots++;
0167         roots[numroots] = -b0 + Sqrt(delta);
0168         if (roots[numroots] > -tol) numroots++;
0169       }
0170       r0    = 1 + Sqrt((RadTube_R0 - localpoint[2]) * (RadTube_R0 + localpoint[2]));
0171       c0    = (localpoint[0] * localpoint[0] + (localpoint[1] - r0) * (localpoint[1] + r0)) * invdirxy2;
0172       delta = b0 * b0 - c0;
0173       if (delta > 0) {
0174         roots[numroots] = -b0 - Sqrt(delta);
0175         if (roots[numroots] > -tol) numroots++;
0176         roots[numroots] = -b0 + Sqrt(delta);
0177         if (roots[numroots] > -tol) numroots++;
0178       }
0179       if (numroots) {
0180         Sort4(roots);
0181       }
0182     } else {
0183       QuarticSolver(coef, roots, numroots);
0184     }
0185 
0186     // iterate over possible solutions to find the first valid one
0187     for (auto i = 0; i < numroots; ++i) {
0188       distance = roots[i];
0189 
0190       // discard obviously incorrect solutions.
0191       // Before Newton refinement is applied a macroscopic number of -10 is used
0192       if (distance < -100) continue;
0193 
0194       Vector3D<Real_t> onsurf = localpoint + distance * dir;
0195 
0196       // apply phi cut
0197       if (!fTorusData->InsidePhi(onsurf, false)) continue;
0198 
0199       // calculate normal
0200       // note: using the normal before the Newton refinement could be dangerous but has proven safe so far.
0201       Vector3D<Real_t> normal = onsurf;
0202       onsurf.z()              = 0.;
0203       onsurf.Normalize();
0204       // onsurf *= fTorusData->Radius();
0205       normal -= onsurf;
0206 
0207       bool hit = flip_exiting ^ (dir.Dot(normal) < 0);
0208       // First solution giving a valid hit wins
0209       if (hit) {
0210 
0211         // refine solution with Newton iterations
0212         s            = roots[i];
0213         Real_t eps   = vecgeom::InfinityLength<Real_t>();
0214         Real_t delta = s * s * s * s + coef.a * s * s * s + coef.b * s * s + coef.c * s + coef.d;
0215         Real_t eps0  = -delta / (4. * s * s * s + 3. * coef.a * s * s + 2. * coef.b * s + coef.c);
0216         int ntry     = 0;
0217         while (Abs(eps) > vecgeom::kTolerance) {
0218           if (Abs(eps0) > 200) break;
0219           s += eps0;
0220           if (Abs(s + eps0) < vecgeom::kTolerance) break;
0221           delta = s * s * s * s + coef.a * s * s * s + coef.b * s * s + coef.c * s + coef.d;
0222           eps   = -delta / (4. * s * s * s + 3. * coef.a * s * s + 2. * coef.b * s + coef.c);
0223           if (Abs(eps) >= Abs(eps0)) break;
0224           ntry++;
0225           // Avoid infinite recursion
0226           if (ntry > 100) break;
0227           eps0 = eps;
0228         }
0229         // discard this solution
0230         if (s < -tol) continue;
0231         // use more accurate solution
0232         distance = Max(Real_t(0.), s);
0233 
0234         distance += tubeDistance;         // add distance to bounding tube (0 if inside)
0235         distance *= fTorusData->Radius(); // denormalize
0236 
0237         return true;
0238       } else {
0239         continue;
0240       }
0241     }
0242     return false;
0243   }
0244 
0245   VECGEOM_FORCE_INLINE
0246   VECCORE_ATT_HOST_DEVICE
0247   /// @brief Computes the isotropic safe distance to unplaced surfaces
0248   /// @param point Point in local surface coordinates
0249   /// @param left_side Flag specifying if the surface is intersected from the left-side that defines the normal
0250   /// @param distance Computed isotropic safety
0251   /// @param onsurf Projection of the point on surface
0252   /// @return Validity of the calculation
0253   bool Safety(Vector3D<Real_t> const &point, bool left_side, Real_t &distance, Vector3D<Real_t> &onsurf) const
0254   {
0255     if (!fTorusData->InsidePhi(point, false)) {
0256       // If the point is not in the phi range, there are other surfaces closer than this one
0257       distance = vecgeom::InfinityLength<Real_t>();
0258       return true;
0259     }
0260     Real_t rho        = point.Perp();
0261     Real_t rTor       = fTorusData->Radius();
0262     Real_t rTube      = fTorusData->RadiusTube();
0263     auto safety       = Sqrt(point.z() * point.z() + (rho - rTor) * (rho - rTor)) - rTube;
0264     bool flip_exiting = left_side ^ fTorusData->IsFlipped();
0265     distance          = flip_exiting ? -safety : safety;
0266 
0267     return distance > -vecgeom::kToleranceDist<Real_t>;
0268   }
0269 };
0270 
0271 } // namespace vgbrep
0272 
0273 #endif
[ Source navigation ] [ Diff markup ] [ Identifier search ] [ general search ]

This page was automatically generated by the 2.3.7 LXR engine.

The LXR team
Valid CSS 2.1! Valid HTML 4.01!