Add the CADsupport module: CAD geometries as exact TGeo solids

Detectors/Base/CMakeLists.txt

@@ -37,7 +37,6 @@ o2_add_library(DetectorsBase
                        src/GlobalParams.cxx
                        src/O2Tessellated.cxx
                        src/TGeoGeometryUtils.cxx
-                       src/CADGeometryUtils.cxx
                PUBLIC_LINK_LIBRARIES FairRoot::Base
                                      O2::CommonUtils
                                      O2::DetectorsCommonDataFormats
@@ -90,6 +89,13 @@ o2_add_test(
   PUBLIC_LINK_LIBRARIES O2::DetectorsBase
   LABELS detectorsbase)
 
+o2_add_test(
+  O2Tessellated
+  SOURCES test/testO2Tessellated.cxx
+  COMPONENT_NAME DetectorsBase
+  PUBLIC_LINK_LIBRARIES O2::DetectorsBase
+  LABELS detectorsbase)
+
 if(BUILD_SIMULATION)
   if (NOT APPLE)
     o2_add_test(

Detectors/Base/include/DetectorsBase/O2Tessellated.h

@@ -88,6 +88,10 @@ class O2Tessellated : public TGeoBBox
   const TBuffer3D& GetBuffer3D(int reqSections, Bool_t localFrame) const override;
   void GetMeshNumbers(int& nvert, int& nsegs, int& npols) const override;
   int GetNmeshVertices() const override { return fNvert; }
+
+  /// Fill \a array with \a npoints points on this solid's boundary: every vertex, then deterministic R2 samples on facet interiors.
+  Bool_t GetPointsOnSegments(Int_t npoints, Double_t* array) const override;
+
   void InspectShape() const override {}
   TBuffer3D* MakeBuffer3D() const override;
   void Print(Option_t* option = "") const override;

Detectors/Base/src/O2Tessellated.cxx

@@ -484,6 +484,67 @@ void O2Tessellated::GetMeshNumbers(int& nvert, int& nsegs, int& npols) const
   npols = GetNfacets();
 }
 
+////////////////////////////////////////////////////////////////////////////////
+/// Fill array with npoints points on the solid's boundary. See the header.
+
+Bool_t O2Tessellated::GetPointsOnSegments(Int_t npoints, Double_t* array) const
+{
+  if (array == nullptr || npoints <= 0 || fVertices.empty()) {
+    return kFALSE;
+  }
+  const int vertexCount = static_cast<int>(fVertices.size());
+  if (npoints < vertexCount) {
+    // Hand the caller back to SetPoints(), which gives it every vertex -- more points than asked
+    // for, all of them exactly on the shape.
+    return kFALSE;
+  }
+  for (int vertexIndex = 0; vertexIndex < vertexCount; ++vertexIndex) {
+    fVertices[vertexIndex].CopyTo(&array[3 * vertexIndex]);
+  }
+
+  const int extraCount = npoints - vertexCount;
+  const int facetCount = static_cast<int>(fFacets.size());
+  if (extraCount == 0) {
+    return kTRUE;
+  }
+  if (facetCount == 0) {
+    for (int extraIndex = 0; extraIndex < extraCount; ++extraIndex) {
+      fVertices[extraIndex % vertexCount].CopyTo(&array[3 * (vertexCount + extraIndex)]);
+    }
+    return kTRUE;
+  }
+
+  // The same deterministic R2 low-discrepancy pair O2BVHSurfaceSolid::GetPointsOnSegments uses:
+  // what a shape hands out must depend on the shape and on nothing else.
+  constexpr double kAlpha1 = 0.7548776662466927;
+  constexpr double kAlpha2 = 0.5698402909980532;
+  for (int extraIndex = 0; extraIndex < extraCount; ++extraIndex) {
+    const int facetIndex =
+      static_cast<int>((static_cast<long long>(extraIndex) * facetCount) / extraCount) % facetCount;
+    const TGeoFacet& facet = fFacets[facetIndex];
+    const int facetVertices = facet.GetNvert();
+    double first = std::fmod(0.5 + kAlpha1 * (extraIndex + 1), 1.);
+    double second = std::fmod(0.5 + kAlpha2 * (extraIndex + 1), 1.);
+    if (first + second > 1.) {
+      first = 1. - first;
+      second = 1. - second;
+    }
+    // A quad facet is two triangles sharing vertex 0; pick one by the parity of the sample index
+    // so both halves are covered.
+    const int cornerB = (facetVertices > 3 && (extraIndex & 1)) ? 2 : 1;
+    const int cornerC = (facetVertices > 3 && (extraIndex & 1)) ? 3 : ((facetVertices > 2) ? 2 : 1);
+    const Vertex_t& vertexA = fVertices[facet[0]];
+    const Vertex_t& vertexB = fVertices[facet[cornerB]];
+    const Vertex_t& vertexC = fVertices[facet[cornerC]];
+    const double weightA = 1. - first - second;
+    double* slot = &array[3 * (vertexCount + extraIndex)];
+    slot[0] = weightA * vertexA.x() + first * vertexB.x() + second * vertexC.x();
+    slot[1] = weightA * vertexA.y() + first * vertexB.y() + second * vertexC.y();
+    slot[2] = weightA * vertexA.z() + first * vertexB.z() + second * vertexC.z();
+  }
+  return kTRUE;
+}
+
 ////////////////////////////////////////////////////////////////////////////////
 /// Creates a TBuffer3D describing *this* shape.
 /// Coordinates are in local reference frame.
@@ -901,6 +962,27 @@ inline Vec3f<T> triangleNormal(const Vec3f<T>& a, const Vec3f<T>& b, const Vec3f
   return normalize(cross(e1, e2));
 }
 
+/// Outward pad of every BVH leaf box, so a facet lies strictly inside the box that stands for it.
+constexpr float kFacetBoxPad = 0.001f;
+
+/// Lowering the ray bound cannot drop a nearer facet while |origin| + |box| + distance stays below
+/// this: the float rounding of ray, box and traversal then stays well inside kFacetBoxPad.
+constexpr double kMaxPruneScale = kFacetBoxPad * (1 << 24) / 8.;
+
+/// The largest hit distance that may be used as a ray bound for this origin and root box.
+template <typename BBox>
+double pruneLimit(const BBox& bbox, const double* point)
+{
+  double origin = 0.;
+  double box = 0.;
+  for (int index = 0; index < 3; ++index) {
+    origin = std::max(origin, std::abs(point[index]));
+    box = std::max({box, std::abs(static_cast<double>(bbox.min[index])),
+                    std::abs(static_cast<double>(bbox.max[index]))});
+  }
+  return kMaxPruneScale - origin - box;
+}
+
 } // end anonymous namespace
 
 ////////////////////////////////////////////////////////////////////////////////
@@ -960,6 +1042,10 @@ Double_t O2Tessellated::DistFromOutside(const Double_t* point, const Double_t* d
 
   static constexpr bool use_robust_traversal = true;
 
+  // the ray object is ours and mutable: bvh2 re-reads tmax at every box test, so lowering it on a
+  // hit prunes the rest of the traversal
+  const double prune_limit = pruneLimit(topnode_bbox, point);
+
   Vertex_t dir_v{dir[0], dir[1], dir[2]};
   // Traverse the BVH and apply concrete object intersection in BVH leafs
   bvh::v2::GrowingStack<Bvh::Index> stack;
@@ -979,6 +1065,9 @@ Double_t O2Tessellated::DistFromOutside(const Double_t* point, const Double_t* d
 
       if (thisdist < local_step) {
         local_step = thisdist;
+        if (local_step <= prune_limit) {
+          ray.tmax = truncate_roundup(local_step);
+        }
       }
     }
     return false; // go on after this
@@ -1023,6 +1112,10 @@ Double_t O2Tessellated::DistFromInside(const Double_t* point, const Double_t* di
 
   static constexpr bool use_robust_traversal = true;
 
+  // as in DistFromOutside: lowering the ray's own tmax on a hit prunes the rest of the traversal
+  const auto rootbox = mybvh->get_root().get_bbox();
+  const double prune_limit = pruneLimit(rootbox, point);
+
   Vertex_t dir_v{dir[0], dir[1], dir[2]};
   // Traverse the BVH and apply concrete object intersection in BVH leafs
   bvh::v2::GrowingStack<Bvh::Index> stack;
@@ -1045,6 +1138,9 @@ Double_t O2Tessellated::DistFromInside(const Double_t* point, const Double_t* di
         rayTriangle(Vertex_t{point[0], point[1], point[2]}, dir_v, v0, v1, v2, 0.);
       if (t < local_step) {
         local_step = t;
+        if (local_step <= prune_limit) {
+          ray.tmax = truncate_roundup(local_step);
+        }
       }
     }
     return false; // go on after this
@@ -1095,12 +1191,12 @@ void O2Tessellated::BuildBVH()
     const auto& v2 = fVertices[facet[1]];
     const auto& v3 = fVertices[facet[2]];
     BBox bbox;
-    bbox.min[0] = std::min(std::min(v1[0], v2[0]), v3[0]) - 0.001f;
-    bbox.min[1] = std::min(std::min(v1[1], v2[1]), v3[1]) - 0.001f;
-    bbox.min[2] = std::min(std::min(v1[2], v2[2]), v3[2]) - 0.001f;
-    bbox.max[0] = std::max(std::max(v1[0], v2[0]), v3[0]) + 0.001f;
-    bbox.max[1] = std::max(std::max(v1[1], v2[1]), v3[1]) + 0.001f;
-    bbox.max[2] = std::max(std::max(v1[2], v2[2]), v3[2]) + 0.001f;
+    bbox.min[0] = std::min(std::min(v1[0], v2[0]), v3[0]) - kFacetBoxPad;
+    bbox.min[1] = std::min(std::min(v1[1], v2[1]), v3[1]) - kFacetBoxPad;
+    bbox.min[2] = std::min(std::min(v1[2], v2[2]), v3[2]) - kFacetBoxPad;
+    bbox.max[0] = std::max(std::max(v1[0], v2[0]), v3[0]) + kFacetBoxPad;
+    bbox.max[1] = std::max(std::max(v1[1], v2[1]), v3[1]) + kFacetBoxPad;
+    bbox.max[2] = std::max(std::max(v1[2], v2[2]), v3[2]) + kFacetBoxPad;
     return bbox;
   };
 

Detectors/Base/src/TGeoGeometryUtils.cxx

@@ -136,7 +136,8 @@ TGeoTessellated* MakeTessellated(const TBuffer3D& buf)
 }
 } // end anonymous namespace
 
-///< Transform any (primitive) TGeoShape to a TGeoTessellated
+///< Transform any (primitive) TGeoShape to a TGeoTessellated.
+/// Display and export only: TGeoTessellated does not navigate (it is tracked as its bounding box); use O2Tessellated for transport.
 TGeoTessellated* TGeoGeometryUtils::TGeoShapeToTGeoTessellated(TGeoShape const* shape)
 {
   auto& buf = shape->GetBuffer3D(TBuffer3D::kRawSizes | TBuffer3D::kRaw | TBuffer3D::kCore, false);

Detectors/Base/test/testO2Tessellated.cxx

@@ -0,0 +1,180 @@
+// Copyright 2019-2026 CERN and copyright holders of ALICE O2.
+// See https://alice-o2.web.cern.ch/copyright for details of the copyright holders.
+// All rights not expressly granted are reserved.
+//
+// This software is distributed under the terms of the GNU General Public
+// License v3 (GPL Version 3), copied verbatim in the file "COPYING".
+//
+// In applying this license CERN does not waive the privileges and immunities
+// granted to it by virtue of its status as an Intergovernmental Organization
+// or submit itself to any jurisdiction.
+/// \author Sandro Wenzel <sandro.wenzel@cern.ch>
+/// \since 2026-09
+
+#define BOOST_TEST_MODULE Test O2Tessellated class
+#define BOOST_TEST_MAIN
+#define BOOST_TEST_DYN_LINK
+#include <boost/test/unit_test.hpp>
+
+#include "DetectorsBase/O2Tessellated.h"
+
+#include "TGeoShape.h"
+
+#include <cmath>
+#include <limits>
+#include <vector>
+
+namespace
+{
+using o2::base::O2Tessellated;
+using Vertex_t = O2Tessellated::Vertex_t;
+
+/// A small deterministic generator, so a failing ray is reproducible from its seed alone.
+class Rng
+{
+ public:
+  explicit Rng(unsigned long long seed) : mState(seed) {}
+  double uniform(double low, double high)
+  {
+    mState = mState * 6364136223846793005ULL + 1442695040888963407ULL;
+    const double unit = static_cast<double>((mState >> 11) & ((1ULL << 53) - 1)) / static_cast<double>(1ULL << 53);
+    return low + unit * (high - low);
+  }
+
+ private:
+  unsigned long long mState;
+};
+
+/// Add the twelve outward-wound triangles of an axis-aligned box.
+void addBox(O2Tessellated& shape, double cx, double cy, double cz, double hx, double hy, double hz)
+{
+  const double x0 = cx - hx, x1 = cx + hx;
+  const double y0 = cy - hy, y1 = cy + hy;
+  const double z0 = cz - hz, z1 = cz + hz;
+  const Vertex_t corner[8] = {{x0, y0, z0}, {x1, y0, z0}, {x1, y1, z0}, {x0, y1, z0}, {x0, y0, z1}, {x1, y0, z1}, {x1, y1, z1}, {x0, y1, z1}};
+  // each quad is wound counter-clockwise seen from outside, so the facet normal points outward
+  const int quad[6][4] = {{0, 3, 2, 1}, {4, 5, 6, 7}, {0, 1, 5, 4}, {2, 3, 7, 6}, {1, 2, 6, 5}, {0, 4, 7, 3}};
+  for (const auto& face : quad) {
+    shape.AddFacet(corner[face[0]], corner[face[1]], corner[face[2]]);
+    shape.AddFacet(corner[face[0]], corner[face[2]], corner[face[3]]);
+  }
+}
+
+/// The Moeller-Trumbore distance used by O2Tessellated's leaf test, repeated here as the oracle.
+double rayTriangleReference(const double* origin, const double* dir, const Vertex_t& v0, const Vertex_t& v1,
+                            const Vertex_t& v2)
+{
+  constexpr double EPS = 1.e-8;
+  const double infinity = std::numeric_limits<double>::infinity();
+  const double e1[3] = {v1[0] - v0[0], v1[1] - v0[1], v1[2] - v0[2]};
+  const double e2[3] = {v2[0] - v0[0], v2[1] - v0[1], v2[2] - v0[2]};
+  const double p[3] = {dir[1] * e2[2] - dir[2] * e2[1], dir[2] * e2[0] - dir[0] * e2[2],
+                       dir[0] * e2[1] - dir[1] * e2[0]};
+  const double det = e1[0] * p[0] + e1[1] * p[1] + e1[2] * p[2];
+  if (std::abs(det) <= EPS) {
+    return infinity;
+  }
+  const double tvec[3] = {origin[0] - v0[0], origin[1] - v0[1], origin[2] - v0[2]};
+  const double invDet = 1.0 / det;
+  const double u = (tvec[0] * p[0] + tvec[1] * p[1] + tvec[2] * p[2]) * invDet;
+  if (u < 0.0 || u > 1.0) {
+    return infinity;
+  }
+  const double q[3] = {tvec[1] * e1[2] - tvec[2] * e1[1], tvec[2] * e1[0] - tvec[0] * e1[2],
+                       tvec[0] * e1[1] - tvec[1] * e1[0]};
+  const double v = (dir[0] * q[0] + dir[1] * q[1] + dir[2] * q[2]) * invDet;
+  if (v < 0.0 || u + v > 1.0) {
+    return infinity;
+  }
+  const double t = e2[0] * q[0] + e2[1] * q[1] + e2[2] * q[2];
+  return (t * invDet > 0.) ? t * invDet : infinity;
+}
+
+/// The unpruned answer: the nearest facet over every facet of the mesh, entering or exiting.
+double bruteForce(const O2Tessellated& shape, const double* origin, const double* dir, bool entering)
+{
+  double best = TGeoShape::Big();
+  for (int facet = 0; facet < shape.GetNfacets(); ++facet) {
+    const auto& description = shape.GetFacet(facet);
+    const Vertex_t& v0 = shape.GetVertex(description[0]);
+    const Vertex_t& v1 = shape.GetVertex(description[1]);
+    const Vertex_t& v2 = shape.GetVertex(description[2]);
+    const double e1[3] = {v1[0] - v0[0], v1[1] - v0[1], v1[2] - v0[2]};
+    const double e2[3] = {v2[0] - v0[0], v2[1] - v0[1], v2[2] - v0[2]};
+    const double normal[3] = {e1[1] * e2[2] - e1[2] * e2[1], e1[2] * e2[0] - e1[0] * e2[2],
+                              e1[0] * e2[1] - e1[1] * e2[0]};
+    const double along = normal[0] * dir[0] + normal[1] * dir[1] + normal[2] * dir[2];
+    // the same facing filter the shape applies: entering facets face the ray, exiting ones face away
+    if (entering ? (along > 0.) : (along <= 0.)) {
+      continue;
+    }
+    best = std::min(best, rayTriangleReference(origin, dir, v0, v1, v2));
+  }
+  return best;
+}
+
+/// Eight boxes in a row, so every axial ray meets sixteen facets and the BVH has many leaves.
+void buildRow(O2Tessellated& shape)
+{
+  for (int index = 0; index < 8; ++index) {
+    addBox(shape, -21. + 6. * index, 0., 0., 2., 3., 4.);
+  }
+  shape.CloseShape(true, false, false);
+}
+} // namespace
+
+BOOST_AUTO_TEST_CASE(PrunedRayQueriesEqualTheBruteForceMinimum)
+{
+  O2Tessellated shape("row");
+  buildRow(shape);
+  BOOST_CHECK_EQUAL(shape.GetNfacets(), 96);
+
+  Rng rng(20260912);
+  int outsideHits = 0;
+  int insideHits = 0;
+  for (int trial = 0; trial < 4000; ++trial) {
+    // origins inside the row and well outside it, so both directions are exercised
+    const double origin[3] = {rng.uniform(-40., 40.), rng.uniform(-12., 12.), rng.uniform(-12., 12.)};
+    double dir[3] = {rng.uniform(-1., 1.), rng.uniform(-1., 1.), rng.uniform(-1., 1.)};
+    const double norm = std::sqrt(dir[0] * dir[0] + dir[1] * dir[1] + dir[2] * dir[2]);
+    if (norm < 1.e-6) {
+      continue;
+    }
+    for (int index = 0; index < 3; ++index) {
+      dir[index] /= norm;
+    }
+
+    const double outside = shape.DistFromOutside(origin, dir, 1, TGeoShape::Big(), nullptr);
+    const double inside = shape.DistFromInside(origin, dir, 1, TGeoShape::Big(), nullptr);
+    const double outsideReference = bruteForce(shape, origin, dir, true);
+    const double insideReference = bruteForce(shape, origin, dir, false);
+
+    BOOST_CHECK_EQUAL(outside, outsideReference);
+    BOOST_CHECK_EQUAL(inside, insideReference);
+    outsideHits += outsideReference < TGeoShape::Big() ? 1 : 0;
+    insideHits += insideReference < TGeoShape::Big() ? 1 : 0;
+  }
+  // the case is only meaningful if the rays really hit the mesh; this sampling gives about 450
+  // entering and 3000 exiting hits
+  BOOST_CHECK_GT(outsideHits, 200);
+  BOOST_CHECK_GT(insideHits, 200);
+}
+
+BOOST_AUTO_TEST_CASE(APrunedRayFindsTheNearestOfManyFacetsAlongIt)
+{
+  O2Tessellated shape("row");
+  buildRow(shape);
+
+  // straight down the row: eight boxes, so sixteen entering and sixteen exiting facets are in line
+  const double origin[3] = {-40., 0., 0.};
+  const double dir[3] = {1., 0., 0.};
+  BOOST_CHECK_EQUAL(shape.DistFromOutside(origin, dir, 1, TGeoShape::Big(), nullptr), 17.);
+  BOOST_CHECK_EQUAL(shape.DistFromOutside(origin, dir, 1, TGeoShape::Big(), nullptr),
+                    bruteForce(shape, origin, dir, true));
+
+  // from inside the first box, the exit is its own far face and not a later box's
+  const double inner[3] = {-21., 0., 0.};
+  BOOST_CHECK_EQUAL(shape.DistFromInside(inner, dir, 1, TGeoShape::Big(), nullptr), 2.);
+  BOOST_CHECK_EQUAL(shape.DistFromInside(inner, dir, 1, TGeoShape::Big(), nullptr),
+                    bruteForce(shape, inner, dir, false));
+}