EIC code displayed by LXR master/​include/​VecGeom/​surfaces/​base/​CpuTypes.h	Source navigation Diff markup Identifier search General search
	Version: master test Revert   Change

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

 #ifndef VECGEOM_SURFACE_CPUTYPES_H
 #define VECGEOM_SURFACE_CPUTYPES_H
 
 #include <VecGeom/surfaces/base/CommonTypes.h>
 #include <VecGeom/surfaces/Model.h>
 #include <set>
 
 namespace vgbrep {
 
 template <typename T>
 VECGEOM_FORCE_INLINE char const *to_cstring(T type)
 {
   return nullptr;
 }
 
 template <>
 VECGEOM_FORCE_INLINE char const *to_cstring<SurfaceType>(SurfaceType type)
 {
   static const char *const data[] = {"no_surf",   "planar", "cylindrical", "conical",
                                      "spherical", "torus",  "elliptical",  "arb4"};
   VECGEOM_ASSERT(size_t(type) * sizeof(const char *) < sizeof(data));
   return data[static_cast<int>(type)];
 }
 
 template <>
 VECGEOM_FORCE_INLINE char const *to_cstring<bool>(bool type)
 {
   if (type) return "true";
   return "false";
 }
 
 template <>
 VECGEOM_FORCE_INLINE char const *to_cstring<FrameType>(FrameType type)
 {
   static const char *const data[] = {"no_frame", "rangeZ", "ring", "z_phi", "rangeSph", "window", "triangle", "quad"};
   VECGEOM_ASSERT(size_t(type) * sizeof(const char *) < sizeof(data));
   return data[static_cast<int>(type)];
 }
 
 template <>
 VECGEOM_FORCE_INLINE char const *to_cstring<AxisType>(AxisType type)
 {
   static const char *const data[] = {"x-axis", "y-axis", "z-axis", "r-axis", "phi-axis", "xy-grid", "no-axis"};
   VECGEOM_ASSERT(size_t(type) * sizeof(const char *) < sizeof(data));
   return data[static_cast<int>(type)];
 }
 
 using LogicExpressionCPU = std::vector<logic_int>;
 
 // Placeholder (on host) for all surfaces belonging to a volume. An array of those will be indexed
 // by the logical volume id. Also an intermediate helper for building portals.
 // Note: the local surfaces defined by solids will have local references that will be changed by
 // the flattening process, depending on the scene on which the parent volume will be flattened
 struct VolumeShellCPU {
   std::vector<int> fSurfaces;             ///< Local surface id's for this volume
   std::vector<int> fExitingSurfaces;      ///< Local surface id's for this volume, excluding virtual ones
   std::vector<int> fEnteringSurfaces;     ///< Local surface id's for all surfaces of daughters, excluding virtual ones
   std::vector<int> fEnteringSurfacesPvol; ///< Global PVol ids for the daughter surfaces of this volume
   std::vector<int> fEnteringSurfacesPvolTrans; ///< Array of ids to daughter placed volume transformations
   std::vector<int> fEnteringSurfacesLvolIds;   ///< Array of ids to daughter logical volumes
   std::vector<int> fDaughterPvolIds;           ///< Global PV Ids of the daughter PVs of this Volume
   std::vector<int> fDaughterPvolTrans;         ///< Transformations of the daughter PVs of this Volume
 
   LogicExpressionCPU fLogic; ///< Logic expression for the solid
   bool fSimplified{false};   ///< The logic was simplified
   int fBVH{0};
 };
 
 /// @brief Helper for dividing a side in equal slices along one axis, keeping frame candidates in each slice
 /// @details The range given by the side extent is divided in equal slices along an axis. Each
 ///          slice intersects a number of frames. The slice is found besed on the crossing point, then the full
 ///          frame loops are reduced to the list of candidates in that slice.
 struct SideDivisionCPU {
   using VecInt_t = std::vector<int>;
   double fStartU{0.};                ///< Division start on primary division axis
   double fStepU{0.};                 ///< Division step on primary division axis
   double fStartV{0.};                ///< Division start on secondary axis (if any)
   double fStepV{0.};                 ///< Division step on secondary axis (if any)
   unsigned short fNslices{0};        ///< Total number of slices
   unsigned short fNslicesU{0};       ///< Number of slices on primary division axis
   unsigned short fNslicesV{1};       ///< Number of slices on secondary division axis (if no axis this must be 1)
   AxisType fAxis{AxisType::kNoAxis}; ///< Division axis
   std::vector<VecInt_t> fSlices;     ///< Array of slices
 
   // Methods
   SideDivisionCPU() = default;
   SideDivisionCPU(AxisType axis, double coord_min, double coord_max, int ncand)
       : fStartU(coord_min), fStepU((coord_max - coord_min) / ncand), fNslices(ncand), fAxis(axis)
   {
     fSlices.insert(fSlices.end(), ncand, {});
   }
 
   SideDivisionCPU(AxisType axis, double coord_minU, double coord_maxU, double coord_minV, double coord_maxV, int ncand)
       : fAxis(axis)
   {
     fNslicesU = 1 + std::sqrt(static_cast<double>(ncand)); // rounding happens
     fNslicesV = ncand / fNslicesU;                         // rounding happens
     fNslices  = fNslicesU * fNslicesV;
     fStartU   = coord_minU;
     fStepU    = (coord_maxU - coord_minU) / fNslicesU;
     fStartV   = coord_minV;
     fStepV    = (coord_maxV - coord_minV) / fNslicesV;
     fSlices.insert(fSlices.end(), fNslices, {});
   }
 
   void AddCandidate(int icand, double coord_min, double coord_max)
   {
     int istart = (coord_min - fStartU - vecgeom::kToleranceDist<double>) / fStepU;
     int iend   = (coord_max - fStartU + vecgeom::kToleranceDist<double>) / fStepU;
     for (int i = istart; i <= iend && i < fNslices; ++i) {
       if (i < 0) continue;
       fSlices[i].push_back(icand);
     }
   }
 
   /// @brief Get the indices of the slices intersected by a bounding box
   /// @param amin lower bbox corner
   /// @param amax upper bbox corner
   /// @param umin lower first index
   /// @param umax upper first index
   /// @param vmin lower second index
   /// @param vmax upper second index
   /// @return number of crossed slices
   int GetNslices(Vector3D<double> const &amin, Vector3D<double> const &amax, int &umin, int &umax, int &vmin,
                  int &vmax) const
   {
     umin = umax = vmin = vmax = 0;
     switch (fAxis) {
     case AxisType::kXY: {
       umin = std::max(0, int((amin[0] - fStartU) / fStepU));
       if (umin >= fNslicesU) return 0;
       umax = std::min(fNslicesU - 1, int((amax[0] - fStartU) / fStepU));
       if (umax < 0) return 0;
       vmin = std::max(0, int((amin[1] - fStartV) / fStepV));
       if (vmin >= fNslicesV) return 0;
       vmax = std::min(fNslicesV - 1, int((amax[1] - fStartV) / fStepV));
       if (vmax < 0) return 0;
       return (umax - umin + 1) * (vmax - vmin + 1);
     }
     case AxisType::kR: {
       auto rmin = amin[0];
       auto rmax = amax[0];
       umin      = std::max(0, int((rmin - fStartU) / fStepU));
       if (umin >= fNslices) return 0;
       umax = std::min(fNslices - 1, int((rmax - fStartU) / fStepU));
       if (umax < 0) return 0;
       return (umax - umin + 1);
     }
     case AxisType::kPhi: {
       VECGEOM_LOG(error) << "GetNdlices not supported for phi division";
       return 0;
     }
     default:
       umin = std::max(0, int((amin[fAxis] - fStartU) / fStepU));
       if (umin >= fNslices) return 0;
       umax = std::min(fNslices - 1, int((amax[fAxis] - fStartU) / fStepU));
       if (umax < 0) return 0;
       return (umax - umin + 1);
     }
   }
 
   void AddCandidate(int icand, double umin, double umax, double vmin, double vmax)
   {
     int istartU = (umin - fStartU - vecgeom::kToleranceDist<double>) / fStepU;
     int iendU   = (umax - fStartU + vecgeom::kToleranceDist<double>) / fStepU;
     int istartV = (vmin - fStartV - vecgeom::kToleranceDist<double>) / fStepV;
     int iendV   = (vmax - fStartV + vecgeom::kToleranceDist<double>) / fStepV;
     for (int i = istartU; i <= iendU && i < fNslicesU; ++i) {
       for (int j = istartV; j <= iendV && j < fNslicesV; ++j) {
         if (i < 0 || j < 0) continue;
         fSlices[i * fNslicesV + j].push_back(icand);
       }
     }
   }
 
   const VecInt_t *GetSlice(int u, int v) const { return &fSlices[u * fNslicesV + v]; }
 
   template <typename Real_t>
   void CopyTo(SideDivision<Real_t> &div, SliceCand *slices, int *candidates) const
   {
     div.fStartU   = fStartU;
     div.fStepU    = fStepU;
     div.fStartV   = fStartV;
     div.fStepV    = fStepV;
     div.fNslices  = fNslices;
     div.fNslicesU = fNslicesU;
     div.fNslicesV = fNslicesV;
     div.fAxis     = fAxis;
     div.fSlices   = slices;
     auto cand     = candidates;
     for (size_t i = 0; i < fSlices.size(); ++i) {
       slices[i].fNcand      = fSlices[i].size();
       slices[i].fCandidates = cand;
       for (int j = 0; j < slices[i].fNcand; ++j)
         cand[j] = fSlices[i][j];
       cand += slices[i].fNcand;
     }
   }
 
   template <typename Real_t>
   size_t GetSizeObject() const
   {
     return sizeof(SideDivision<Real_t>);
   }
 
   template <typename Real_t>
   size_t GetSizeSlices() const
   {
     return fNslices * sizeof(SliceCand);
   }
 
   size_t GetNcandidates() const
   {
     size_t size = 0;
     for (auto i = 0; i < fNslices; ++i)
       size += fSlices[i].size();
     return size;
   }
 
   size_t GetSizeCandidates() const { return GetNcandidates() * sizeof(int); }
 
   template <typename Real_t>
   size_t GetSize() const
   {
     size_t size = GetSizeObject<Real_t>() + GetSizeSlices<Real_t>() + GetSizeCandidates();
     return size;
   }
 
   void Print() const
   {
     if (fAxis == AxisType::kXY) {
       std::cout << to_cstring(fAxis) << " division: fNslices = " << fNslices << ", fNslicesU = " << fNslicesU
                 << ", fNslicesV = " << fNslicesV << ", fStartU = " << fStartU << ", fStepU = " << fStepU
                 << ", fStartV = " << fStartV << ", fStepV = " << fStepV << ", efficiency = " << 100 * Efficiency()
                 << " %\n";
     } else {
       std::cout << to_cstring(fAxis) << " division: fNslices = " << fNslices << ", fStart = " << fStartU
                 << ", fStep = " << fStepU << ", efficiency = " << 100 * Efficiency() << " %\n";
     }
     for (auto i = 0; i < fNslices; ++i)
       std::cout << fSlices[i].size() << " ";
     std::cout << "\n";
   }
 
   double Efficiency() const
   {
     // VECGEOM_ASSERT(fSlices[0].size() * fSlices[fNslices - 1].size() > 0); // except kXY
     double average = 0;
     for (auto i = 0; i < fNslices; ++i)
       average += fSlices[i].size();
     average /= fNslices;
     double eff = (fNslices - average) / (average * (fNslices - 1));
     return eff;
   }
 };
 
 // Surface data used only on CPU during the conversion process
 template <typename Real_t>
 struct CPUsurfData {
   using VecInt_t       = std::vector<int>;
   using VecChar_t      = std::vector<char>;
   using MultimapInt_t  = std::multimap<long, int>;
   using SurfData_t     = SurfData<Real_t>;
   using EllipData_t    = EllipData<Real_t>;
   using TorusData_t    = TorusData<Real_t>;
   using Arb4Data_t     = Arb4Data<Real_t>;
   using WindowMask_t   = WindowMask<Real_t>;
   using RingMask_t     = RingMask<Real_t>;
   using ZPhiMask_t     = ZPhiMask<Real_t>;
   using TriangleMask_t = TriangleMask<Real_t>;
   using QuadMask_t     = QuadrilateralMask<Real_t>;
 
   std::vector<WindowMask_t> fWindowMasks; ///< rectangular masks
   std::vector<RingMask_t> fRingMasks;     ///< ring masks
   std::vector<ZPhiMask_t> fZPhiMasks;     ///< cylindrical masks
   std::vector<TriangleMask_t> fTriangleMasks;
   std::vector<QuadMask_t> fQuadMasks; ///< quadrilateral masks
 
   std::vector<EllipData_t> fEllipData;                                         ///< data for elliptical surfaces
   std::vector<TorusData_t> fTorusData;                                         ///< data for torus surfaces
   std::vector<Arb4Data_t> fArb4Data;                                           ///< data for Arb4 surfaces
   std::vector<TransformationMP<Real_t>> fPVolTrans;                            ///< Transformations to placed volumes
   std::vector<FramedSurface<Real_t, TransformationMP<Real_t>>> fLocalSurfaces; ///< local surfaces per logical volume
   std::vector<FramedSurface<Real_t, TransformationMP<Real_t>>> fFramedSurf;    ///< global surfaces
   std::vector<CommonSurface<Real_t>> fCommonSurfaces;                          ///< common surfaces
   std::vector<VolumeShellCPU> fShells;                                         ///< vector of local volume surfaces
   std::vector<VolumeShellCPU> fSceneShells;                                    ///< vector of scene volume surfaces
 
   VecInt_t fSceneStartIndex;            ///< Start indices for data indexed by state id (per scene)
   VecInt_t fSceneTouchables;            ///< Number of touchables (per scene)
   std::vector<MultimapInt_t> fSurfHash; ///< maps rotation hash index to a list of common surface id's (per scene)
 
   std::vector<VecInt_t> fCandidatesEntering; ///< list of entering candidates: scene0...,scene1...
   std::vector<VecInt_t> fCandidatesExiting;  ///< list of exiting candidates: scene0...,scene1...
   std::vector<VecInt_t> fFrameIndEntering; ///< list of start frame indices for entering candidates: scene0...,scene1...
   std::vector<VecInt_t> fFrameIndExiting;  ///< list of start frame indices for exiting candidates: scene0...,scene1...
   std::vector<VecChar_t> fSidesEntering;   ///< list of relevant sides for entering candidates: scene0...,scene1...
   std::vector<VecChar_t> fSidesExiting;    ///< list of relevant sides for exiting candidates: scene0...,scene1...
 
   std::vector<SideDivisionCPU> fSideDivisions; ///< list of divisions for all sides
 
 private:
   CPUsurfData() = default;
 
 public:
   static VECGEOM_FORCE_INLINE CPUsurfData<Real_t> &Instance()
   {
     static CPUsurfData<Real_t> gCPUsurfdata;
     return gCPUsurfdata;
   }
 
   void Clear()
   {
     // Dispose of surface data and shrink the container
     std::vector<WindowMask_t>().swap(fWindowMasks);
     std::vector<RingMask_t>().swap(fRingMasks);
     std::vector<ZPhiMask_t>().swap(fZPhiMasks);
     std::vector<TriangleMask_t>().swap(fTriangleMasks);
     std::vector<QuadMask_t>().swap(fQuadMasks);
     std::vector<EllipData_t>().swap(fEllipData);
     std::vector<TorusData_t>().swap(fTorusData);
     std::vector<TransformationMP<Real_t>>().swap(fPVolTrans);
     std::vector<FramedSurface<Real_t, TransformationMP<Real_t>>>().swap(fLocalSurfaces);
     std::vector<FramedSurface<Real_t, TransformationMP<Real_t>>>().swap(fFramedSurf);
     std::vector<CommonSurface<Real_t>>().swap(fCommonSurfaces);
     std::vector<VolumeShellCPU>().swap(fShells);
     std::vector<VolumeShellCPU>().swap(fSceneShells);
     VecInt_t().swap(fSceneStartIndex);
     VecInt_t().swap(fSceneTouchables);
     std::vector<MultimapInt_t>().swap(fSurfHash);
     std::vector<VecInt_t>().swap(fCandidatesEntering);
     std::vector<VecInt_t>().swap(fCandidatesExiting);
     std::vector<VecInt_t>().swap(fFrameIndEntering);
     std::vector<VecInt_t>().swap(fFrameIndExiting);
     std::vector<VecChar_t>().swap(fSidesEntering);
     std::vector<VecChar_t>().swap(fSidesExiting);
     std::vector<SideDivisionCPU>().swap(fSideDivisions);
   }
 
   VecInt_t &GetCandidatesEntering(int scene_id, int state_id)
   {
     return fCandidatesEntering[fSceneStartIndex[scene_id] + state_id];
   }
 
   VecInt_t &GetCandidatesExiting(int scene_id, int state_id)
   {
     return fCandidatesExiting[fSceneStartIndex[scene_id] + state_id];
   }
 
   VecInt_t &GetFrameIndEntering(int scene_id, int state_id)
   {
     return fFrameIndEntering[fSceneStartIndex[scene_id] + state_id];
   }
 
   VecInt_t &GetFrameIndExiting(int scene_id, int state_id)
   {
     return fFrameIndExiting[fSceneStartIndex[scene_id] + state_id];
   }
 
   VecChar_t &GetSidesEntering(int scene_id, int state_id)
   {
     return fSidesEntering[fSceneStartIndex[scene_id] + state_id];
   }
 
   VecChar_t &GetSidesExiting(int scene_id, int state_id)
   {
     return fSidesExiting[fSceneStartIndex[scene_id] + state_id];
   }
 
   WindowMask_t const &GetWindowMask(int id) const { return fWindowMasks[id]; }
   RingMask_t const &GetRingMask(int id) const { return fRingMasks[id]; }
   ZPhiMask_t const &GetZPhiMask(int id) const { return fZPhiMasks[id]; }
   TriangleMask_t const &GetTriangleMask(int id) const { return fTriangleMasks[id]; }
   QuadMask_t const &GetQuadMask(int id) const { return fQuadMasks[id]; }
 
   EllipData_t const &GetEllipData(int id) const { return fEllipData[id]; }
   TorusData_t const &GetTorusData(int id) const { return fTorusData[id]; }
   Arb4Data_t const &GetArb4Data(int id) const { return fArb4Data[id]; }
 
   /// @brief Trampoline function to to the frame checker
   /// @param f1 Parent framed surface
   /// @param f2 Child framed surface
   /// @return Parent frame classification with respect to child
   FrameIntersect CheckFrames(FramedSurface<Real_t, TransformationMP<Real_t>> const &f1,
                              FramedSurface<Real_t, TransformationMP<Real_t>> const &f2,
                              TransformationMP<Real_t> *tscene = nullptr)
   {
     auto log_not_supported = [&]() {
       VECGEOM_LOG(error) << "Embedding check " << to_cstring(f1.fFrame.type) << " - " << to_cstring(f2.fFrame.type)
                          << " not supported";
     };
     TransformationMP<Real_t> t1 = f1.fTrans;
     TransformationMP<Real_t> t2 = f2.fTrans;
     if (tscene) t2 *= (*tscene);
     TransformationMP<Real_t> trans = t2 * t1.Inverse();
     trans.SetProperties();
 
     switch (f1.fFrame.type) {
     case FrameType::kRing: {
       RingMask_t const mask1 = GetRingMask(f1.fFrame.id);
       switch (f2.fFrame.type) {
       case FrameType::kRing: {
         RingMask_t const mask2 = GetRingMask(f2.fFrame.id);
         return FrameChecker<Real_t, RingMask_t, RingMask_t>::CheckFrames(mask1, mask2, trans);
       }
       case FrameType::kWindow: {
         WindowMask_t const mask2 = GetWindowMask(f2.fFrame.id);
         return FrameChecker<Real_t, RingMask_t, WindowMask_t>::CheckFrames(mask1, mask2, trans);
       }
       case FrameType::kTriangle: {
         TriangleMask_t const mask2 = GetTriangleMask(f2.fFrame.id);
         return FrameChecker<Real_t, RingMask_t, TriangleMask_t>::CheckFrames(mask1, mask2, trans);
       }
       case FrameType::kQuadrilateral: {
         QuadMask_t const mask2 = GetQuadMask(f2.fFrame.id);
         return FrameChecker<Real_t, RingMask_t, QuadMask_t>::CheckFrames(mask1, mask2, trans);
       }
       default:
         log_not_supported();
       };
       break;
     }
     case FrameType::kZPhi: {
       ZPhiMask_t const mask1 = GetZPhiMask(f1.fFrame.id);
       switch (f2.fFrame.type) {
       case FrameType::kZPhi: {
         ZPhiMask_t const mask2 = GetZPhiMask(f2.fFrame.id);
         return FrameChecker<Real_t, ZPhiMask_t, ZPhiMask_t>::CheckFrames(mask1, mask2, trans);
       }
       default:
         log_not_supported();
       };
       break;
     }
     case FrameType::kWindow: {
       WindowMask_t const mask1 = GetWindowMask(f1.fFrame.id);
       switch (f2.fFrame.type) {
       case FrameType::kRing: {
         RingMask_t const mask2 = GetRingMask(f2.fFrame.id);
         return FrameChecker<Real_t, WindowMask_t, RingMask_t>::CheckFrames(mask1, mask2, trans);
       }
       case FrameType::kWindow: {
         WindowMask_t const mask2 = GetWindowMask(f2.fFrame.id);
         return FrameChecker<Real_t, WindowMask_t, WindowMask_t>::CheckFrames(mask1, mask2, trans);
       }
       case FrameType::kTriangle: {
         TriangleMask_t const mask2 = GetTriangleMask(f2.fFrame.id);
         return FrameChecker<Real_t, WindowMask_t, TriangleMask_t>::CheckFrames(mask1, mask2, trans);
       }
       case FrameType::kQuadrilateral: {
         QuadMask_t const mask2 = GetQuadMask(f2.fFrame.id);
         return FrameChecker<Real_t, WindowMask_t, QuadMask_t>::CheckFrames(mask1, mask2, trans);
       }
       default:
         log_not_supported();
       };
       break;
     }
     case FrameType::kTriangle: {
       TriangleMask_t const mask1 = GetTriangleMask(f1.fFrame.id);
       switch (f2.fFrame.type) {
       case FrameType::kRing: {
         RingMask_t const mask2 = GetRingMask(f2.fFrame.id);
         return FrameChecker<Real_t, TriangleMask_t, RingMask_t>::CheckFrames(mask1, mask2, trans);
       }
       case FrameType::kWindow: {
         WindowMask_t const mask2 = GetWindowMask(f2.fFrame.id);
         return FrameChecker<Real_t, TriangleMask_t, WindowMask_t>::CheckFrames(mask1, mask2, trans);
       }
       case FrameType::kTriangle: {
         TriangleMask_t const mask2 = GetTriangleMask(f2.fFrame.id);
         return FrameChecker<Real_t, TriangleMask_t, TriangleMask_t>::CheckFrames(mask1, mask2, trans);
       }
       case FrameType::kQuadrilateral: {
         QuadMask_t const mask2 = GetQuadMask(f2.fFrame.id);
         return FrameChecker<Real_t, TriangleMask_t, QuadMask_t>::CheckFrames(mask1, mask2, trans);
       }
       default:
         log_not_supported();
       };
       break;
     }
     case FrameType::kQuadrilateral: {
       QuadMask_t const mask1 = GetQuadMask(f1.fFrame.id);
       switch (f2.fFrame.type) {
       case FrameType::kRing: {
         RingMask_t const mask2 = GetRingMask(f2.fFrame.id);
         return FrameChecker<Real_t, QuadMask_t, RingMask_t>::CheckFrames(mask1, mask2, trans);
       }
       case FrameType::kWindow: {
         WindowMask_t const mask2 = GetWindowMask(f2.fFrame.id);
         return FrameChecker<Real_t, QuadMask_t, WindowMask_t>::CheckFrames(mask1, mask2, trans);
       }
       case FrameType::kTriangle: {
         TriangleMask_t const mask2 = GetTriangleMask(f2.fFrame.id);
         return FrameChecker<Real_t, QuadMask_t, TriangleMask_t>::CheckFrames(mask1, mask2, trans);
       }
       case FrameType::kQuadrilateral: {
         QuadMask_t const mask2 = GetQuadMask(f2.fFrame.id);
         return FrameChecker<Real_t, QuadMask_t, QuadMask_t>::CheckFrames(mask1, mask2, trans);
       }
       default:
         log_not_supported();
       };
       break;
     }
     default:
       log_not_supported();
     };
     return FrameIntersect::kNoIntersect;
   }
 };
 
 } // namespace vgbrep
 
 #endif

This page was automatically generated by the 2.3.7 LXR engine. The LXR team