Replace the TPC half-space cuts by bounded boxes

This replaces the fifteen TGeoHalfSpace cuts in the TPC support geometry by
equivalent bounded boxes and adds a unit test.

- TGeoHalfSpace is unbounded, so ROOT's GDML writer drops every volume that
  uses one and VGM cannot convert it to native Geant4 geometry.
- All fifteen occurrences in AliceO2 are TPC support structures, and every one
  of them is subtracted from a solid smaller than 10 cm in each direction.
- Each is now a cube of 100 cm half-size placed with one face on the cut plane,
  built by the new TGeoGeometryUtils::makeHalfSpaceBox.
- A subtraction takes its bounding box from the left-hand solid only, so the
  composite shapes are unchanged in extent.
- The cut term is written in parentheses. A trailing "shape:matrix" is unsafe:
  TGeoManager::Parse reads the last top-level ":" of an expression that already
  contains a top-level ")" as a transformation of the whole expression and then
  drops it, which would leave an unplaced cube swallowing the parent solid.
- The new test compares both forms point by point and along random rays, for
  the fifteen real planes, 200 arbitrary ones, and one compound expression of
  the shape that makes the parentheses necessary.

https://its.cern.ch/jira/browse/O2-4534

Co-Authored-By: Claude Opus 5 <noreply@anthropic.com>

Author: sawenzel

Detectors/Base/CMakeLists.txt

@@ -83,6 +83,13 @@ o2_target_root_dictionary(DetectorsBase
                                   include/DetectorsBase/O2Tessellated.h
                                   )
 
+o2_add_test(
+  HalfSpaceBox
+  SOURCES test/testHalfSpaceBox.cxx
+  COMPONENT_NAME DetectorsBase
+  PUBLIC_LINK_LIBRARIES O2::DetectorsBase
+  LABELS detectorsbase)
+
 if(BUILD_SIMULATION)
   if (NOT APPLE)
     o2_add_test(

Detectors/Base/include/DetectorsBase/TGeoGeometryUtils.h

@@ -30,6 +30,20 @@ class TGeoGeometryUtils
  public:
   ///< Transform any (primitive) TGeoShape to a tessellated representation
   static TGeoTessellated* TGeoShapeToTGeoTessellated(TGeoShape const*);
+
+  ///< Create a bounded stand-in for the half-space { x : (x - p) . n <= 0 }, which is what
+  ///< TGeoHalfSpace describes. Registers a cube of half-size `reach` under `name` and its
+  ///< placement under "<name>_tr". The stand-in agrees with the half-space everywhere within
+  ///< a distance `reach` of `p`, so `reach` must exceed the extent of the solid it is
+  ///< subtracted from. Unlike TGeoHalfSpace, the result can be exported to GDML and
+  ///< converted to native Geant4 geometry.
+  ///<
+  ///< Write the term in a composite expression **in parentheses**, as "-(<name>:<name>_tr)".
+  ///< A trailing "shape:matrix" is not safe: TGeoManager::Parse takes the last top-level ":"
+  ///< of an expression that already contains a top-level ")" to be a transformation of the
+  ///< whole expression, warns "no geometrical transformation allowed at this level" and then
+  ///< drops it - leaving an unplaced cube at the origin that swallows the parent solid.
+  static void makeHalfSpaceBox(const char* name, const double p[3], const double n[3], double reach);
 };
 
 } // namespace base

Detectors/Base/src/TGeoGeometryUtils.cxx

@@ -16,7 +16,11 @@
 #include <DetectorsBase/TGeoGeometryUtils.h>
 #include <TGeoShape.h>
 #include <TGeoTessellated.h>
+#include <TGeoBBox.h>
+#include <TGeoMatrix.h>
 #include <TBuffer3D.h>
+#include <TString.h>
+#include <cmath>
 #include <vector>
 
 namespace o2
@@ -140,5 +144,43 @@ TGeoTessellated* TGeoGeometryUtils::TGeoShapeToTGeoTessellated(TGeoShape const*
   return tes;
 }
 
+
+///< Bounded stand-in for a TGeoHalfSpace
+void TGeoGeometryUtils::makeHalfSpaceBox(const char* name, const double p[3], const double n[3], double reach)
+{
+  // TGeoHalfSpace contains the points x with (p - x) . n >= 0, and normalizes n itself.
+  double nn[3] = {n[0], n[1], n[2]};
+  const double norm = std::sqrt(nn[0] * nn[0] + nn[1] * nn[1] + nn[2] * nn[2]);
+  for (auto& c : nn) {
+    c /= norm;
+  }
+
+  // an orthonormal triad (u, v, nn); the seed is chosen to stay away from nn
+  double a[3] = {1., 0., 0.};
+  if (std::abs(nn[0]) > 0.9) {
+    a[0] = 0.;
+    a[1] = 1.;
+  }
+  double u[3] = {a[1] * nn[2] - a[2] * nn[1], a[2] * nn[0] - a[0] * nn[2], a[0] * nn[1] - a[1] * nn[0]};
+  const double unorm = std::sqrt(u[0] * u[0] + u[1] * u[1] + u[2] * u[2]);
+  for (auto& c : u) {
+    c /= unorm;
+  }
+  const double v[3] = {nn[1] * u[2] - nn[2] * u[1], nn[2] * u[0] - nn[0] * u[2], nn[0] * u[1] - nn[1] * u[0]};
+
+  // rotation taking the local z axis onto nn (TGeoRotation stores the matrix row-wise,
+  // so the images of the local axes are its columns)
+  const double m[9] = {u[0], v[0], nn[0], u[1], v[1], nn[1], u[2], v[2], nn[2]};
+  auto* rot = new TGeoRotation(TString::Format("%s_rot", name));
+  rot->SetMatrix(m);
+
+  // centre the cube one half-size behind the plane, so its +z face lies on the plane
+  auto* tr = new TGeoCombiTrans(TString::Format("%s_tr", name), p[0] - reach * nn[0], p[1] - reach * nn[1],
+                                p[2] - reach * nn[2], rot);
+  tr->RegisterYourself();
+
+  new TGeoBBox(name, reach, reach, reach);
+}
+
 } // namespace base
 } // namespace o2

Detectors/Base/test/testHalfSpaceBox.cxx

@@ -0,0 +1,210 @@
+// Copyright 2019-2020 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.
+
+/// \file testHalfSpaceBox.cxx
+/// \author Sandro Wenzel (CERN)
+/// \brief Checks that TGeoGeometryUtils::makeHalfSpaceBox reproduces TGeoHalfSpace
+
+#define BOOST_TEST_MODULE Test HalfSpaceBox
+#define BOOST_TEST_MAIN
+#define BOOST_TEST_DYN_LINK
+#include <boost/test/unit_test.hpp>
+
+#include "DetectorsBase/TGeoGeometryUtils.h"
+#include "TGeoManager.h"
+#include "TGeoBBox.h"
+#include "TGeoTube.h"
+#include "TGeoMatrix.h"
+#include "TGeoHalfSpace.h"
+#include "TGeoCompositeShape.h"
+#include "TMath.h"
+#include "TRandom3.h"
+#include "TString.h"
+#include <cmath>
+#include <vector>
+
+namespace
+{
+struct Plane {
+  const char* label;
+  double p[3];
+  double n[3];
+};
+
+// The fifteen half-space cuts of the TPC support structures (Detectors/TPC/simulation/src/Detector.cxx).
+// Largest solid any of them is subtracted from is 1.65 x 1.85 x 8.9 cm, hence the 10 cm parent below.
+std::vector<Plane> tpcPlanes()
+{
+  const double slope = TMath::Tan(22. * TMath::DegToRad());
+  const double intp = 1.245;
+  const double b = slope * slope + 1.;
+  const double p1[3] = {intp * slope / b, -intp / b, 0.};
+  const double p2[3] = {-intp * slope / b, -intp / b, 0.};
+  return {
+    {"sp1", {p1[0], p1[1], 0.}, {-p1[0], -p1[1], 0.}},
+    {"sp2", {p2[0], p2[1], 0.}, {-p2[0], -p2[1], 0.}},
+    {"cutil1", {0., 0.105, 0.}, {0., 1., 0.}},
+    {"cutomh1", {0., -1.05, -3.4}, {0., -TMath::Tan(30. * TMath::DegToRad()), 1.}},
+    {"cutomh2", {0., -1.05, 3.4}, {0., -TMath::Tan(30. * TMath::DegToRad()), -1.}},
+    {"cutomh3", {-1.65, 0., -0.9}, {TMath::Tan(75. * TMath::DegToRad()), 0., 1.}},
+    {"cutomh4", {-1.65, 0., 0.9}, {TMath::Tan(75. * TMath::DegToRad()), 0., -1.}},
+    {"cutomh5", {1.65, -1.05, 0.}, {-1., -TMath::Tan(20. * TMath::DegToRad()), 0.}},
+    {"cutohs1", {0., -0.186, 0.}, {0., -1., 0.}},
+    {"cutmmh1", {-1.65, 0., -0.9}, {8., 0., 8. * TMath::Tan(13. * TMath::DegToRad())}},
+    {"cutmmh2", {-1.65, 0., 0.9}, {8., 0., -8. * TMath::Tan(13. * TMath::DegToRad())}},
+    {"cutmmh3", {0., 1.85, -2.8}, {0., -6.1, 6.1 * TMath::Tan(20. * TMath::DegToRad())}},
+    {"cutmmh4", {0., 1.85, 2.8}, {0., -6.1, -6.1 * TMath::Tan(20. * TMath::DegToRad())}},
+    {"cutmmh5", {0.75, 0., -8.9}, {2.4 * TMath::Tan(30. * TMath::DegToRad()), 0., 2.4}},
+    {"cutmmh6", {0.75, 0., 8.9}, {2.4 * TMath::Tan(30. * TMath::DegToRad()), 0., -2.4}}};
+}
+
+// Compares "parent - halfspace" against "parent - box:box_tr" on random points and random rays.
+// Points closer than kSurfaceBand to the plane are skipped: on the surface itself the two
+// implementations may legitimately round to different sides.
+void compare(const TString& tag, const double p[3], const double n[3], double parentHalfSize, double reach,
+             TRandom3& rnd, int nPoints, int nRays, double& maxDistDiff)
+{
+  const double nl = std::sqrt(n[0] * n[0] + n[1] * n[1] + n[2] * n[2]);
+  BOOST_REQUIRE(nl > 1e-6);
+  constexpr double kSurfaceBand = 1e-9;
+
+  new TGeoBBox(TString::Format("parent_%s", tag.Data()).Data(), parentHalfSize, parentHalfSize, parentHalfSize);
+  new TGeoHalfSpace(TString::Format("hs_%s", tag.Data()).Data(), const_cast<double*>(p), const_cast<double*>(n));
+  o2::base::TGeoGeometryUtils::makeHalfSpaceBox(TString::Format("bx_%s", tag.Data()).Data(), p, n, reach);
+
+  auto* ref = new TGeoCompositeShape(TString::Format("ref_%s", tag.Data()),
+                                     TString::Format("parent_%s-hs_%s", tag.Data(), tag.Data()));
+  auto* box = new TGeoCompositeShape(TString::Format("new_%s", tag.Data()),
+                                     TString::Format("parent_%s-(bx_%s:bx_%s_tr)", tag.Data(), tag.Data(), tag.Data()));
+
+  const double range = 1.2 * parentHalfSize;
+  for (int k = 0; k < nPoints; ++k) {
+    double x[3];
+    for (int i = 0; i < 3; ++i) {
+      x[i] = rnd.Uniform(-range, range);
+    }
+    const double d = ((x[0] - p[0]) * n[0] + (x[1] - p[1]) * n[1] + (x[2] - p[2]) * n[2]) / nl;
+    if (std::abs(d) < kSurfaceBand) {
+      continue;
+    }
+    if (ref->Contains(x) != box->Contains(x)) {
+      BOOST_REQUIRE_MESSAGE(false, "containment differs for " << tag.Data() << " at (" << x[0] << "," << x[1] << ","
+                                                             << x[2] << "), distance to plane " << d);
+    }
+  }
+
+  for (int k = 0; k < nRays; ++k) {
+    double x[3], dir[3];
+    for (int i = 0; i < 3; ++i) {
+      x[i] = rnd.Uniform(-3. * range, 3. * range);
+      dir[i] = rnd.Uniform(-1., 1.);
+    }
+    const double dn = std::sqrt(dir[0] * dir[0] + dir[1] * dir[1] + dir[2] * dir[2]);
+    if (dn < 1e-6) {
+      continue;
+    }
+    for (int i = 0; i < 3; ++i) {
+      dir[i] /= dn;
+    }
+    const bool inside = ref->Contains(x);
+    if (inside != box->Contains(x)) {
+      continue; // a point sitting on the surface; covered by the containment loop above
+    }
+    const double d1 = inside ? ref->DistFromInside(x, dir, 3) : ref->DistFromOutside(x, dir, 3);
+    const double d2 = inside ? box->DistFromInside(x, dir, 3) : box->DistFromOutside(x, dir, 3);
+    if (d1 > 1e15 && d2 > 1e15) {
+      continue; // both miss
+    }
+    maxDistDiff = std::max(maxDistDiff, std::abs(d1 - d2));
+  }
+}
+} // namespace
+
+BOOST_AUTO_TEST_CASE(HalfSpaceBox_reproduces_TGeoHalfSpace)
+{
+  auto* geom = new TGeoManager("halfspacetest", "half-space replacement test");
+  TRandom3 rnd(20240101);
+  double maxDistDiff = 0.;
+
+  // the real TPC cuts
+  for (const auto& pl : tpcPlanes()) {
+    compare(pl.label, pl.p, pl.n, 10., 100., rnd, 200000, 20000, maxDistDiff);
+  }
+
+  // and a spread of arbitrary planes, to pin the rotation for normals in every octant
+  for (int i = 0; i < 200; ++i) {
+    double p[3], n[3];
+    for (int k = 0; k < 3; ++k) {
+      p[k] = rnd.Uniform(-5., 5.);
+      n[k] = rnd.Uniform(-1., 1.);
+    }
+    if (std::sqrt(n[0] * n[0] + n[1] * n[1] + n[2] * n[2]) < 1e-3) {
+      continue;
+    }
+    compare(TString::Format("rnd%d", i), p, n, 10., 100., rnd, 20000, 2000, maxDistDiff);
+  }
+
+  // the two shapes are not bit-identical, but they must agree to double round-off
+  BOOST_CHECK_SMALL(maxDistDiff, 1e-9);
+  BOOST_TEST_MESSAGE("maximum ray-distance difference: " << maxDistDiff);
+  delete geom;
+}
+
+// The composite expressions of the TPC support structures are not all of the simple
+// "parent - cut" shape: tpcihs6 subtracts a union and two placed tubes first. That shape is
+// what makes a *trailing* "cut:matrix" term unsafe to write unparenthesised, so keep a case
+// with the same structure.
+BOOST_AUTO_TEST_CASE(HalfSpaceBox_in_a_compound_expression)
+{
+  auto* geom = new TGeoManager("halfspacetest2", "half-space replacement, compound expression");
+  const double shift[3] = {0., -0.175, 0.};
+  const double p[3] = {0., 0.105, 0.};
+  const double n[3] = {0., 1., 0.};
+
+  new TGeoBBox("tpcihs1", 4.7, 0.66, 2.35);
+  new TGeoBBox("tpcihs2", 4.7, 0.485, 1.0, const_cast<double*>(shift));
+  new TGeoBBox("tpcihs3", 1.5, 0.485, 2.35, const_cast<double*>(shift));
+  new TGeoTube("tpcihs4", 0.0, 2.38, 0.1);
+  auto* trans2 = new TGeoTranslation("trans2", 0.0, 2.84, 2.25);
+  trans2->RegisterYourself();
+  auto* trans3 = new TGeoTranslation("trans3", 0.0, 2.84, -2.25);
+  trans3->RegisterYourself();
+  new TGeoHalfSpace("cutil1", const_cast<double*>(p), const_cast<double*>(n));
+  o2::base::TGeoGeometryUtils::makeHalfSpaceBox("bcutil1", p, n, 100.);
+
+  auto* ref = new TGeoCompositeShape(
+    "ref_tpcihs6", "tpcihs1-(tpcihs2+tpcihs3)-(tpcihs4:trans2)-(tpcihs4:trans3)-cutil1");
+  auto* box = new TGeoCompositeShape(
+    "new_tpcihs6", "tpcihs1-(tpcihs2+tpcihs3)-(tpcihs4:trans2)-(tpcihs4:trans3)-(bcutil1:bcutil1_tr)");
+
+  TRandom3 rnd(20240102);
+  long inRef = 0, inBox = 0;
+  for (int k = 0; k < 2000000; ++k) {
+    double x[3];
+    for (int i = 0; i < 3; ++i) {
+      x[i] = rnd.Uniform(-6., 6.);
+    }
+    if (std::abs(x[1] - p[1]) < 1e-9) {
+      continue;
+    }
+    const bool a = ref->Contains(x);
+    const bool b = box->Contains(x);
+    inRef += a;
+    inBox += b;
+    if (a != b) {
+      BOOST_REQUIRE_MESSAGE(false, "containment differs at (" << x[0] << "," << x[1] << "," << x[2] << ")");
+    }
+  }
+  // guards against both shapes being empty, which would make the comparison vacuous
+  BOOST_CHECK_GT(inRef, 0);
+  BOOST_CHECK_EQUAL(inRef, inBox);
+  delete geom;
+}