VecGeom/surfaces/Model.h

0001 #ifndef VECGEOM_SURFACE_MODEL_H_
0002 #define VECGEOM_SURFACE_MODEL_H_
0003
0004 #include <VecGeom/navigation/NavigationState.h>
0005 #include <VecGeom/surfaces/base/CommonTypes.h>
0006 #include <VecGeom/surfaces/base/Equations.h>
0007 #include <VecGeom/surfaces/surf/SurfaceImpl.h>
0008 #include <VecGeom/surfaces/mask/FrameMasks.h>
0009 #include <VecGeom/surfaces/cuda/DeviceStorage.h>
0010
0011 namespace vgbrep {
0012
0013 ///< Forward declaration of surface data structure used in navigation
0014 template <typename Real_t>
0015 struct SurfData;
0016
0017 struct Frame;
0018 using Extent = Frame;
0019
0020 /// @brief Unplaced half-space surface type.
0021 /// @details Unplaced surfaces are infinite half-spaces, having a normal side convention:
0022 ///   - All unplaced planes are (xOy), having the normal oriented on positive z
0023 ///   - Unplaced cylinders, cones and tori have the z axis as axis of symmetry. Normals pointing outwards.
0024 ///   - Unplaced spheres have the origin as center, normal pointing outwards.
0025 /// The type does not store the surface data, but only an id to an external storage.
0026 template <typename Real_t>
0027 struct UnplacedSurface {
0028   SurfaceType type{SurfaceType::kPlanar}; ///< surface type
0029   int id{-1};                             ///< surface id
0030   Real_t fRadius{0.};
0031   Real_t fSlope{0.};
0032
0033   VECCORE_ATT_HOST_DEVICE
0034   Real_t Radius() const { return std::abs(Real_t(fRadius)); }
0035   VECCORE_ATT_HOST_DEVICE
0036   Real_t RadiusZ(Real_t z) const { return std::abs(fRadius) + z * fSlope; }
0037   VECCORE_ATT_HOST_DEVICE
0038   bool IsFlipped() const { return fRadius < 0; }
0039
0040   UnplacedSurface() = default;
0041
0042   UnplacedSurface(SurfaceType stype, int sid = -1)
0043   {
0044     type    = stype;
0045     id      = sid;
0046     fRadius = static_cast<Real_t>(0);
0047     fSlope  = static_cast<Real_t>(0);
0048   }
0049
0050   template <typename Real_i>
0051   UnplacedSurface(SurfaceType stype, int sid = -1, Real_i radius = 0, Real_i slope = 0, bool flip = false)
0052   {
0053     type    = stype;
0054     id      = sid;
0055     fRadius = flip ? static_cast<Real_t>(-radius) : static_cast<Real_t>(radius);
0056     fSlope  = static_cast<Real_t>(slope);
0057   }
0058
0059   template <typename Real_i>
0060   UnplacedSurface(const UnplacedSurface<Real_i> &other)
0061       : type(other.type), id(other.id), fRadius(static_cast<Real_t>(other.fRadius)),
0062         fSlope(static_cast<Real_t>(other.fSlope))
0063   {
0064   }
0065
0066   /// @brief A local point is inside if behind the normal within tolerance
0067   /// @tparam Real_t Floating-point precision type
0068   /// @param point Point in the local surface coordinates
0069   /// @param surfdata data container with the surface data
0070   /// @param flip flipping the tolerance for inside (needed for negated booleans)
0071   /// @return Inside half-space
0072   template <typename DataContainer>
0073   VECCORE_ATT_HOST_DEVICE bool Inside(Vector3D<Real_t> const &point, DataContainer const &data, bool flip,
0074                                       Real_t tol = vecgeom::kToleranceStrict<Real_t>) const
0075   {
0076     switch (type) {
0077     case SurfaceType::kPlanar:
0078       return SurfaceHelper<SurfaceType::kPlanar, Real_t>().Inside(point, flip, tol);
0079     case SurfaceType::kCylindrical:
0080       return SurfaceHelper<SurfaceType::kCylindrical, Real_t>().Inside(point, flip, fRadius, tol);
0081     case SurfaceType::kConical:
0082       return SurfaceHelper<SurfaceType::kConical, Real_t>().Inside(point, flip, fRadius, fSlope, tol);
0083     case SurfaceType::kElliptical:
0084       return SurfaceHelper<SurfaceType::kElliptical, Real_t>(data.GetEllipData(id)).Inside(point, flip, tol);
0085     case SurfaceType::kSpherical:
0086       return SurfaceHelper<SurfaceType::kSpherical, Real_t>().Inside(point, flip, fRadius, tol);
0087     case SurfaceType::kTorus:
0088       return SurfaceHelper<SurfaceType::kTorus, Real_t>(data.GetTorusData(id)).Inside(point, flip, tol);
0089     case SurfaceType::kArb4:
0090       return SurfaceHelper<SurfaceType::kArb4, Real_t>(data.GetArb4Data(id)).Inside(point, flip, tol);
0091     default:
0092       return false;
0093     };
0094   }
0095
0096   /// @brief Find signed distance to next intersection from local point.
0097   /// @tparam Real_t Floating-point precision type
0098   /// @param point Point in the local surface coordinates
0099   /// @param dir Direction in the local surface coordinates
0100   /// @param left_side Flag specifying if the surface is intersected from the left-side that defines the normal
0101   /// @param surfdata Surface data storage.
0102   /// @param distance Computed distance to surface
0103   /// @return Validity of the intersection
0104   VECCORE_ATT_HOST_DEVICE bool Intersect(Vector3D<Real_t> const &point, Vector3D<Real_t> const &dir, bool left_side,
0105                                          SurfData<Real_t> const &surfdata, Real_t &distance, bool &two_solutions,
0106                                          Real_t &safety) const
0107   {
0108     switch (type) {
0109     case SurfaceType::kPlanar:
0110       return SurfaceHelper<SurfaceType::kPlanar, Real_t>().Intersect(point, dir, left_side, distance, two_solutions,
0111                                                                      safety);
0112     case SurfaceType::kCylindrical:
0113       return SurfaceHelper<SurfaceType::kCylindrical, Real_t>().Intersect(point, dir, left_side, distance,
0114                                                                           two_solutions, safety, fRadius);
0115     case SurfaceType::kConical:
0116       return SurfaceHelper<SurfaceType::kConical, Real_t>().Intersect(point, dir, left_side, distance, two_solutions,
0117                                                                       safety, fRadius, fSlope);
0118     case SurfaceType::kElliptical:
0119       return SurfaceHelper<SurfaceType::kElliptical, Real_t>(surfdata.GetEllipData(id))
0120           .Intersect(point, dir, left_side, distance, two_solutions, safety);
0121     case SurfaceType::kSpherical:
0122       return SurfaceHelper<SurfaceType::kSpherical, Real_t>().Intersect(point, dir, left_side, distance, two_solutions,
0123                                                                         safety, fRadius);
0124     case SurfaceType::kTorus:
0125       return SurfaceHelper<SurfaceType::kTorus, Real_t>(surfdata.GetTorusData(id))
0126           .Intersect(point, dir, left_side, distance, two_solutions, safety);
0127     case SurfaceType::kArb4:
0128       return SurfaceHelper<SurfaceType::kArb4, Real_t>(surfdata.GetArb4Data(id))
0129           .Intersect(point, dir, left_side, distance, two_solutions, safety);
0130     default:
0131       return false;
0132     };
0133   }
0134
0135   /// @brief Computes the isotropic safe distance to unplaced surfaces
0136   /// @tparam Real_t Precision type for parameters
0137   /// @param point Point in local surface coordinates
0138   /// @param left_side Flag specifying if the surface is intersected from the left-side that defines the normal
0139   /// @param surfdata Surface data storage
0140   /// @param distance Computed isotropic safety
0141   /// @param onsurf Projection of the point on surface
0142   /// @return
0143   VECCORE_ATT_HOST_DEVICE bool Safety(Vector3D<Real_t> const &point, bool left_side, SurfData<Real_t> const &surfdata,
0144                                       Real_t &distance, Vector3D<Real_t> &onsurf) const
0145   {
0146     switch (type) {
0147     case SurfaceType::kPlanar:
0148       return SurfaceHelper<SurfaceType::kPlanar, Real_t>().Safety(point, left_side, distance, onsurf);
0149     case SurfaceType::kCylindrical:
0150       return SurfaceHelper<SurfaceType::kCylindrical, Real_t>().Safety(point, left_side, distance, onsurf, fRadius);
0151     case SurfaceType::kConical:
0152       return SurfaceHelper<SurfaceType::kConical, Real_t>().Safety(point, left_side, distance, onsurf, fRadius, fSlope);
0153     case SurfaceType::kElliptical:
0154       return SurfaceHelper<SurfaceType::kElliptical, Real_t>(surfdata.GetEllipData(id))
0155           .Safety(point, left_side, distance, onsurf);
0156     case SurfaceType::kSpherical:
0157       return SurfaceHelper<SurfaceType::kSpherical, Real_t>().Safety(point, left_side, distance, onsurf, fRadius);
0158     case SurfaceType::kTorus:
0159       return SurfaceHelper<SurfaceType::kTorus, Real_t>(surfdata.GetTorusData(id))
0160           .Safety(point, left_side, distance, onsurf);
0161     case SurfaceType::kArb4:
0162       return SurfaceHelper<SurfaceType::kArb4, Real_t>(surfdata.GetArb4Data(id))
0163           .Safety(point, left_side, distance, onsurf);
0164     default:
0165       return false;
0166     };
0167   }
0168 };
0169
0170 /// @brief A frame delimiting the real solid surface on an infinite half-space
0171 struct Frame {
0172   FrameType type{FrameType::kWindow}; ///< frame type
0173   int id{-1};                         ///< frame mask id
0174
0175   Frame() = default;
0176   Frame(FrameType mtype, int mid = -1) : type(mtype), id(mid) {}
0177
0178   // A function to check if local point is within the Frame's mask.
0179   template <typename Real_t>
0180   VECCORE_ATT_HOST_DEVICE bool Inside(Vector3D<Real_t> const &local, SurfData<Real_t> const &surfdata) const
0181   {
0182     switch (type) {
0183     case FrameType::kRing:
0184       return surfdata.GetRingMask(id).Inside(local);
0185     case FrameType::kZPhi:
0186       return surfdata.GetZPhiMask(id).Inside(local);
0187     case FrameType::kWindow:
0188       return surfdata.GetWindowMask(id).Inside(local);
0189     // TODO: Support these
0190     case FrameType::kTriangle:
0191       return surfdata.GetTriangleMask(id).Inside(local);
0192     case FrameType::kQuadrilateral:
0193       return surfdata.GetQuadMask(id).Inside(local);
0194     case FrameType::kRangeZ:
0195       /*return (local[2] > vecgeom::MakeMinusTolerant<true>(u[0]) &&
0196               local[2] < vecgeom::MakePlusTolerant<true>(u[1]));*/
0197     case FrameType::kRangeSph:
0198       /*return (rsq > vecgeom::MakeMinusTolerantSquare<true>(u[0]) &&
0199               rsq < vecgeom::MakePlusTolerantSquare<true>(u[1]));*/
0200     default:
0201       // unhandled
0202       return false;
0203     };
0204     return false;
0205   }
0206
0207   /// @brief A function dispatcher to compute the safety for the frame.
0208   /// @tparam Real_t
0209   /// @param local Point on surface in local coordinates
0210   /// @param safetySurf Safety to the support surface
0211   /// @param surfdata Surface data storage
0212   /// @return Safety to the framed surface.
0213   template <typename Real_t>
0214   VECCORE_ATT_HOST_DEVICE Real_t Safety(Vector3D<Real_t> const &local, Real_t safetySurf,
0215                                         SurfData<Real_t> const &surfdata, bool &valid) const
0216   {
0217     switch (type) {
0218     case FrameType::kRing:
0219       return surfdata.GetRingMask(id).Safety(local, safetySurf, valid);
0220     case FrameType::kZPhi:
0221       return surfdata.GetZPhiMask(id).Safety(local, safetySurf, valid);
0222     case FrameType::kWindow:
0223       return surfdata.GetWindowMask(id).Safety(local, safetySurf, valid);
0224     case FrameType::kTriangle:
0225       return surfdata.GetTriangleMask(id).Safety(local, safetySurf, valid);
0226     case FrameType::kQuadrilateral:
0227       return surfdata.GetQuadMask(id).Safety(local, safetySurf, valid);
0228     case FrameType::kRangeZ:
0229     case FrameType::kRangeSph:
0230     default:
0231       // unhandled
0232       valid = false;
0233     };
0234     return Real_t(0);
0235   }
0236
0237   template <typename Real_t, typename DataContainer>
0238   VECCORE_ATT_HOST_DEVICE void Extent3D(Vector3D<Real_t> &aMin, Vector3D<Real_t> &aMax, DataContainer const &data) const
0239   {
0240     switch (type) {
0241     case FrameType::kRing:
0242       // printf("Frame type: Ring. ");
0243       data.GetRingMask(id).Extent3D(aMin, aMax);
0244       break;
0245     case FrameType::kZPhi:
0246       // printf("Frame type: ZPhi. ");
0247       data.GetZPhiMask(id).Extent3D(aMin, aMax);
0248       break;
0249     case FrameType::kWindow:
0250       // printf("Frame type: Window. ");
0251       data.GetWindowMask(id).Extent3D(aMin, aMax);
0252       break;
0253     case FrameType::kTriangle:
0254       // printf("Frame type: Triangle. ");
0255       data.GetTriangleMask(id).Extent3D(aMin, aMax);
0256       break;
0257     case FrameType::kQuadrilateral:
0258       // printf("Frame type: Quadrilateral. ");
0259       data.GetQuadMask(id).Extent3D(aMin, aMax);
0260       break;
0261     case FrameType::kRangeZ:;
0262     case FrameType::kRangeSph:;
0263     default:
0264       // unhandled
0265       return;
0266     };
0267   }
0268 };
0269
0270 // This holds the transformation of the surface
0271 // with respect to the frame of the ancestor volume onto which this surface is flattened.
0272 //
0273 // Example of a 4-level hierarchy flattened on 2 levels:
0274 // A, B, C, D, E, F = logical volumes
0275 // t1, t2, t3, t4, t5 = local volume transformations
0276 //
0277 //            A                     A + C(t2) + E(t2*t4)   (one scene given by volume A)
0278 //      (t1) / \ (t2)                      |
0279 //          B   C                     (t1) |
0280 //    (t3) /     \ (t4)   -->              |
0281 //        D       E                 B + D(t3) + F(t3*t5)   (another scene given by volume B)
0282 //   (t5) |
0283 //        F
0284 //
0285 // In the above example, say F has a surface S positioned with a local ransformation (ts).
0286 // The global corresponding to S will have a total transformation (t1 * t3 * t5 * ts).
0287 // However, it the tree is flattened on two scenes as above (one level is A and the other level is B)
0288 // then the local transformation of S will be just (t3 * t5 * ts). Its global transformation
0289 // will be just (t1). In this approach, the transformation for a global surface is always obtained
0290 // by multiplying the full scene transformation with the local surface transformation.
0291 //
0292 // The advantage of this approach is that it gives full flexibility for chosing the flattened
0293 // volumes, and a given local surface can be referenced by multiple portals (less memory)
0294
0295 /// @brief A placed surface on a scene having a frame and a navigation state associated to a touchable
0296 template <typename Real_t, typename Transformation_t>
0297 struct FramedSurface {
0298   using NavState_t = vecgeom::NavigationState::Value_t;
0299   UnplacedSurface<Real_t> fSurface; ///< Surface identifier
0300   Frame fFrame;                     ///< Frame
0301   Transformation_t fTrans;  ///< 3D Transformation of the surface in the compacted sub-hierarchy top volume frame or 2D
0302                             ///< transformation of the frame with respect to its common surface
0303   int fTraversal{-2};       ///< Pre-computed traversal surface on the other side
0304                             ///<   -2       = unknown
0305                             ///<   -1       = no frame on the other side
0306                             ///<    i       = frame index on the other side
0307   int fParent{-1};          ///< Index of the first parent frame on the common surface
0308   int fLogicId{0};          ///< Logic flag for surface:
0309                             ///<   0        = non-Bool
0310                             ///<   positive = true logic surface
0311                             ///<   negative = negated logic surface
0312   int fSceneCS{0};          ///< The frame may belong to a daughter scene common surface
0313   int fSceneCSind{0};       ///< Index of the corresponding frame on the scene CS
0314   unsigned fSurfIndex{0};   ///< Surface index in the volume shell (can be optimized by compacting with fLogicId)
0315   NavIndex_t fState{0};     ///< sub-path navigation state id in the parent scene
0316   bool fNeverCheck{false};  ///< The frame should never be checked
0317   bool fEmbedded{false};    ///< The surface is embedded in the parent surface if any
0318   bool fEmbedding{true};    ///< The frame always embeds daughter state frames if on the same CS
0319   bool fOverlapping{false}; ///< The frame is overlapping another frame and requires a relocation after crossing
0320   bool fVirtualParent{false}; ///< The parent frame is a virtual surface of a boolean
0321   bool fSkipConvexity{false}; ///< whether the convexity check for booleans can be skipped (only the case for end caps
0322                               ///< of elliptical tubes)
0323
0324   FramedSurface() = default;
0325   FramedSurface(UnplacedSurface<Real_t> const &unplaced, Frame const &frame, Transformation_t trans,
0326                 NavIndex_t index = 0, const bool never_check = 0)
0327       : fSurface(unplaced), fFrame(frame), fTrans(trans), fState(index), fNeverCheck(never_check)
0328   {
0329   }
0330
0331   template <typename Real_i>
0332   FramedSurface(const FramedSurface<Real_i, TransformationMP<Real_i>> &other)
0333       : fSurface(other.fSurface), fFrame(other.fFrame), fTrans(other.fTrans), fTraversal(other.fTraversal),
0334         fParent(other.fParent), fLogicId(other.fLogicId), fSceneCS(other.fSceneCS), fSceneCSind(other.fSceneCSind),
0335         fSurfIndex(other.fSurfIndex), fState(other.fState), fNeverCheck(other.fNeverCheck), fEmbedded(other.fEmbedded),
0336         fEmbedding(other.fEmbedding), fOverlapping(other.fOverlapping), fVirtualParent(other.fVirtualParent),
0337         fSkipConvexity(other.fSkipConvexity)
0338   {
0339   }
0340
0341   /// Sorting by decreasing state depth and increasing state index
0342   bool operator<(FramedSurface const &other) const
0343   {
0344     using vecgeom::NavigationState;
0345     auto level1 = NavigationState::GetLevelImpl(fState);
0346     auto level2 = NavigationState::GetLevelImpl(other.fState);
0347     if (level1 > level2) return true;
0348     if (level1 < level2) return false;
0349     if (!fEmbedding && other.fEmbedding) return true;
0350     if (fEmbedding && !other.fEmbedding) return false;
0351     if (fState < other.fState) return true;
0352     if (fState > other.fState) return false;
0353     // If states are identical, sort by fSurfIndex
0354     if (fSurfIndex == other.fSurfIndex) VECGEOM_LOG(critical) << "### Found frames with same state and surface index";
0355     if (fSurfIndex < other.fSurfIndex) return true;
0356     return false;
0357   }
0358
0359   template <typename DataContainer>
0360   VECCORE_ATT_HOST_DEVICE void Extent3D(Vector3D<Real_t> &aMin, Vector3D<Real_t> &aMax, DataContainer const &data) const
0361   {
0362     aMin.Set(0, 0, 0);
0363     aMax.Set(0, 0, 0);
0364     if (fFrame.type == FrameType::kNoFrame) {
0365       VECGEOM_LOG(critical) << "Extent3D should not be called for kNoFrame frames";
0366       return;
0367     }
0368     if (fNeverCheck) {
0369       // Special cases where the bounding box is computed based on the UnplacedSurface
0370       if (fSurface.type == SurfaceType::kTorus)
0371         SurfaceHelper<SurfaceType::kTorus, Real_t>(data.GetTorusData(fSurface.id)).Extent3D(aMin, aMax);
0372       else if (fSurface.type == SurfaceType::kArb4)
0373         SurfaceHelper<SurfaceType::kArb4, Real_t>(data.GetArb4Data(fSurface.id)).Extent3D(aMin, aMax);
0374     } else
0375       fFrame.Extent3D(aMin, aMax, data);
0376   }
0377
0378   /// @brief Get the parent state index for the framed surface
0379   /// @param parent_state Parent state variable filled with the return value
0380   /// @return True if there is a parent, false if this is a top scene surface
0381   VECCORE_ATT_HOST_DEVICE
0382   VECGEOM_FORCE_INLINE
0383   bool GetParentState(NavIndex_t &parent_state) const
0384   {
0385     using vecgeom::NavigationState;
0386     parent_state = fState;
0387     NavigationState::PopImpl(parent_state);
0388     return (fState > 0);
0389   }
0390
0391   /// @brief Get logical volume id for this frame
0392   /// @return Volume id.
0393   VECCORE_ATT_HOST_DEVICE
0394   VECGEOM_FORCE_INLINE
0395   int VolumeId() const
0396   {
0397     // We may need to cache the logical volume id in the surface directly
0398     return vecgeom::NavigationState::GetLogicalIdImpl(fState);
0399   }
0400
0401   VECCORE_ATT_HOST_DEVICE
0402   VECCORE_FORCE_NOINLINE
0403   void PrintState() const { vecgeom::NavigationState::PrintTopImpl(fState); }
0404
0405   /// Transform point and direction to the local frame
0406   void Transform(Vector3D<Real_t> const &point, Vector3D<Real_t> const &dir, Vector3D<Real_t> &localpoint,
0407                  Vector3D<Real_t> &localdir) const
0408   {
0409     localpoint = fTrans.Transform(point);
0410     localdir   = fTrans.TransformDirection(dir);
0411   }
0412
0413   ///< Check if the propagated point on surface is within the frame
0414   VECCORE_ATT_HOST_DEVICE bool InsideFrame(Vector3D<Real_t> const &point, SurfData<Real_t> const &surfdata) const
0415   {
0416     Vector3D<Real_t> localpoint(point);
0417     // For single-frame surfaces, fTrans is zero, so it may be worth testing this.
0418     if (!fTrans.IsIdentity()) localpoint = fTrans.Transform(point);
0419     return fFrame.Inside(localpoint, surfdata);
0420   }
0421
0422   /// @brief Calculate the shortest distance (or an underestimate) from a point on the surface
0423   ///  to the surface frame.
0424   /// @tparam Real_t Precision type for parameters
0425   /// @param point Point on the surface
0426   /// @param safetySurf Safety to the surface
0427   /// @param surfdata Surface data storage
0428   /// @return Combined safety surface+frame
0429   VECCORE_ATT_HOST_DEVICE Real_t SafetyFrame(Vector3D<Real_t> const &point, Real_t safetySurf,
0430                                              SurfData<Real_t> const &surfdata, bool &valid) const
0431   {
0432     Vector3D<Real_t> localpoint(point);
0433     // For single-frame surfaces, fTrans is zero, so it may be worth testing this.
0434     if (!fTrans.IsIdentity()) localpoint = fTrans.Transform(point);
0435     return fFrame.Safety(localpoint, safetySurf, surfdata, valid);
0436   }
0437
0438   /// @brief Calculate the shortest distance (or an underestimate) from a point on the surface
0439   ///  to the surface frame.
0440   /// @tparam Real_t Precision type for parameters
0441   /// @param point Point on the surface, in the surface reference frame
0442   /// @param safetySurf Safety to the surface
0443   /// @param surfdata Surface data storage
0444   /// @return Combined safety surface+frame
0445   VECCORE_ATT_HOST_DEVICE Real_t LocalSafetyFrame(Vector3D<Real_t> const &point, Real_t safetySurf,
0446                                                   SurfData<Real_t> const &surfdata, bool &valid) const
0447   {
0448     return fFrame.Safety(point, safetySurf, surfdata, valid);
0449   }
0450 };
0451
0452 /// @brief A side represents all common placed surfaces
0453 struct Side {
0454   int fNumParents{0};      ///< number of different parent volumes contributing to this side
0455   int fNsurf{0};           ///< Number of placed surfaces on this side
0456   int fDivision{-1};       ///< Division helper for the side
0457   int *fSurfaces{nullptr}; ///< [fNsurf] Array of placed surfaces on this side
0458
0459   Side() = default;
0460
0461   // Add existing placed surface to this side
0462   int AddSurface(int isurf)
0463   {
0464     // Re-allocate policy to keep memory footprint low
0465     // Sides have to be relocated for GPU in contiguous memory
0466     int *surfaces = new int[fNsurf + 1];
0467     for (auto i = 0; i < fNsurf; ++i)
0468       surfaces[i] = fSurfaces[i];
0469     surfaces[fNsurf++] = isurf;
0470     delete[] fSurfaces;
0471     fSurfaces = surfaces;
0472     return fNsurf - 1;
0473   }
0474
0475   VECCORE_ATT_HOST_DEVICE VECGEOM_FORCE_INLINE int GetTopSurfaceIndex() const
0476   {
0477     return (fNsurf > 0) ? GetSurfaceIndex(fNsurf - 1) : -1;
0478   }
0479   VECCORE_ATT_HOST_DEVICE VECGEOM_FORCE_INLINE int GetSurfaceIndex(int isurf) const { return fSurfaces[isurf]; }
0480   VECCORE_ATT_HOST_DEVICE VECGEOM_FORCE_INLINE bool HasChildren() const { return fNsurf > fNumParents; }
0481
0482   template <typename Real_t>
0483   VECCORE_ATT_HOST_DEVICE VECGEOM_FORCE_INLINE FramedSurface<Real_t, vecgeom::Transformation2DMP<Real_t>> &GetSurface(
0484       int index, SurfData<Real_t> const &surfdata) const
0485   {
0486     return surfdata.fFramedSurf[fSurfaces[index]];
0487   }
0488
0489   template <typename Real_t>
0490   VECCORE_ATT_HOST_DEVICE VECGEOM_FORCE_INLINE FramedSurface<Real_t, vecgeom::Transformation2DMP<Real_t>> const &
0491   TopSurface(SurfData<Real_t> const &surfdata) const
0492   {
0493     return surfdata.fFramedSurf[fSurfaces[fNsurf - 1]];
0494   }
0495
0496   size_t size() const { return sizeof(Side) + fNsurf * sizeof(int); }
0497
0498   void CopyTo(char *buffer)
0499   {
0500     // to be implemented
0501   }
0502 };
0503
0504 template <typename Real_t, typename Real_i>
0505 struct SideIterator {
0506   SliceCand const *fSlice{nullptr}; ///< Division slice to be iterated
0507   int fNsurf{0};                    ///< number of surfaces to iterate
0508   int fCrt{0};                      ///< Current index
0509
0510   VECCORE_ATT_HOST_DEVICE
0511   VECGEOM_FORCE_INLINE
0512   SideIterator(const Side &side, Vector3D<Real_i> const &onsurf, SurfData<Real_t> const &surfdata, int increment = 1,
0513                int istart = 0)
0514       : fNsurf(side.fNsurf), fCrt(istart)
0515   {
0516     if (side.fDivision >= 0) {
0517       auto const &div = surfdata.fSideDivisions[side.fDivision];
0518       auto ind        = div.GetSliceIndex(onsurf);
0519       if (ind >= 0) {
0520         fSlice    = &div.fSlices[ind];
0521         fNsurf    = fSlice->fNcand;
0522         int inc01 = (increment + 1) >> 1; // [-1/1 -> 0/1]
0523         int start = (1 - inc01) * (fSlice->fNcand - 1);
0524         int end   = inc01 * (fSlice->fNcand - 1);
0525         for (fCrt = start; fCrt * increment < end; fCrt += increment)
0526           if (increment * fSlice->fCandidates[fCrt] >= increment * istart) break;
0527       }
0528     }
0529   }
0530
0531   VECCORE_ATT_HOST_DEVICE
0532   VECGEOM_FORCE_INLINE
0533   bool Done() const { return (fCrt < 0) || (fCrt >= fNsurf); }
0534
0535   VECCORE_ATT_HOST_DEVICE
0536   VECGEOM_FORCE_INLINE
0537   void SetNextIndex(int ind)
0538   {
0539     if (fSlice) {
0540       int i;
0541       for (i = fCrt + 1; i < fSlice->fNcand; ++i)
0542         if (fSlice->fCandidates[i] >= ind) break;
0543       fCrt = i - 1;
0544     } else {
0545       fCrt = ind - 1;
0546     }
0547   }
0548
0549   VECCORE_ATT_HOST_DEVICE VECGEOM_FORCE_INLINE int operator()() { return fSlice ? fSlice->fCandidates[fCrt] : fCrt; }
0550
0551   VECCORE_ATT_HOST_DEVICE
0552   VECGEOM_FORCE_INLINE
0553   SideIterator &operator++()
0554   {
0555     fCrt++;
0556     return *this;
0557   }
0558
0559   VECCORE_ATT_HOST_DEVICE
0560   VECGEOM_FORCE_INLINE
0561   SideIterator &operator--()
0562   {
0563     fCrt--;
0564     return *this;
0565   }
0566 };
0567
0568 /// @brief A common surface made of two sides, having a global transformation.
0569 template <typename Real_t>
0570 struct CommonSurface {
0571   SurfaceType fType{SurfaceType::kPlanar}; ///< Type of surface
0572   int fSceneId{0};                         ///< Scene id. if negative, it is a top scene id
0573   TransformationMP<Real_t> fTrans;         ///< Transformation of the first left frame
0574   NavIndex_t fDefaultState{0};             ///< The default state for this surface (deepest mother)
0575   Side fLeftSide;                          ///< Left-side (behind normal)
0576   Side fRightSide;                         ///< Right-side (alongside normal)
0577   bool fFlipped;                           ///< whether the common surface is flipped due to boolean negation
0578
0579   CommonSurface() = default;
0580
0581   CommonSurface(SurfaceType type, int global_surf, bool flipped) : fType(type)
0582   {
0583     // Add by default the first surface to the left side
0584     fLeftSide.AddSurface(global_surf);
0585     fFlipped = flipped;
0586   };
0587
0588   template <typename Real_i>
0589   CommonSurface(const CommonSurface<Real_i> &other)
0590       : fType(other.fType), fSceneId(other.fSceneId), fTrans(other.fTrans), fDefaultState(other.fDefaultState),
0591         fLeftSide(other.fLeftSide), fRightSide(other.fRightSide), fFlipped(other.fFlipped)
0592   {
0593   }
0594
0595   VECCORE_ATT_HOST_DEVICE
0596   VECGEOM_FORCE_INLINE
0597   int GetSceneId() const { return std::abs(fSceneId); }
0598
0599   VECCORE_ATT_HOST_DEVICE
0600   VECGEOM_FORCE_INLINE
0601   bool IsSceneSurface() const { return (fSceneId < 0); }
0602
0603   ///< Get the normal to the surface from a point on surface
0604   void GetNormal(Vector3D<Real_t> const &point, Vector3D<Real_t> &normal, SurfData<Real_t> const &surfdata) const
0605   {
0606     Vector3D<Real_t> localnorm;
0607     // point to local frame
0608     auto const &trans      = fTrans;
0609     auto localpoint        = trans.Transform(point);
0610     auto const &framedsurf = fLeftSide.GetSurface(0, surfdata);
0611     framedsurf.fSurface.GetNormal(localpoint, localnorm, surfdata);
0612     trans.InverseTransformDirection(localnorm, normal);
0613   }
0614 };
0615
0616 /// @brief A volume shell holding indices for all placed surfaces belonging to a volume.
0617 struct VolumeShell {
0618   LogicExpression fLogic;                   ///< Logic expression for local surfaces
0619   int fBVH{0};                              ///< The BVH index for this logical volume
0620   int fNsurf{0};                            ///< Number of local surfaces
0621   int fNExitingSurfaces{0};                 ///< Number of framed exiting surfaces
0622   int fNEnteringSurfaces{0};                ///< Number of framed entering surfaces
0623   int fNDaughterPvols{0};                   ///< Number of daughter volumes
0624   int *fSurfaces{nullptr};                  ///< Local surface id's
0625   int *fExitingSurfaces{nullptr};           ///< Indices of real surfaces in the fSurfaces array
0626   int *fEnteringSurfaces{nullptr};          ///< List of framed entering surfaces
0627   int *fEnteringSurfacesPvol{nullptr};      ///< Pvol id per framed surface
0628   int *fEnteringSurfacesPvolTrans{nullptr}; ///< Array of ids to daughter placed volume transformations
0629   int *fEnteringSurfacesLvolIds{nullptr};   ///< Array of ids to daughter logical volumes
0630   int *fDaughterPvolIds{nullptr};           ///< Global PV Ids of the daughter PVs of this Volume
0631   int *fDaughterPvolTrans{nullptr};         ///< Transformations of the daughter PVs of this Volume
0632 };
0633 } // namespace vgbrep
0634
0635 #endif