|
| 1 | +// Copyright 2019-2026 CERN and copyright holders of ALICE O2. |
| 2 | +// See https://alice-o2.web.cern.ch/copyright for details of the copyright holders. |
| 3 | +// All rights not expressly granted are reserved. |
| 4 | +// |
| 5 | +// This software is distributed under the terms of the GNU General Public |
| 6 | +// License v3 (GPL Version 3), copied verbatim in the file "COPYING". |
| 7 | +// |
| 8 | +// In applying this license CERN does not waive the privileges and immunities |
| 9 | +// granted to it by virtue of its status as an Intergovernmental Organization |
| 10 | +// or submit itself to any jurisdiction. |
| 11 | +/// \author Sandro Wenzel <sandro.wenzel@cern.ch> |
| 12 | +/// \since 2026-09 |
| 13 | + |
| 14 | +#define BOOST_TEST_MODULE Test O2Tessellated class |
| 15 | +#define BOOST_TEST_MAIN |
| 16 | +#define BOOST_TEST_DYN_LINK |
| 17 | +#include <boost/test/unit_test.hpp> |
| 18 | + |
| 19 | +#include "DetectorsBase/O2Tessellated.h" |
| 20 | + |
| 21 | +#include "TGeoShape.h" |
| 22 | + |
| 23 | +#include <cmath> |
| 24 | +#include <limits> |
| 25 | +#include <vector> |
| 26 | + |
| 27 | +namespace |
| 28 | +{ |
| 29 | +using o2::base::O2Tessellated; |
| 30 | +using Vertex_t = O2Tessellated::Vertex_t; |
| 31 | + |
| 32 | +/// A small deterministic generator, so a failing ray is reproducible from its seed alone. |
| 33 | +class Rng |
| 34 | +{ |
| 35 | + public: |
| 36 | + explicit Rng(unsigned long long seed) : mState(seed) {} |
| 37 | + double uniform(double low, double high) |
| 38 | + { |
| 39 | + mState = mState * 6364136223846793005ULL + 1442695040888963407ULL; |
| 40 | + const double unit = static_cast<double>((mState >> 11) & ((1ULL << 53) - 1)) / static_cast<double>(1ULL << 53); |
| 41 | + return low + unit * (high - low); |
| 42 | + } |
| 43 | + |
| 44 | + private: |
| 45 | + unsigned long long mState; |
| 46 | +}; |
| 47 | + |
| 48 | +/// Add the twelve outward-wound triangles of an axis-aligned box. |
| 49 | +void addBox(O2Tessellated& shape, double cx, double cy, double cz, double hx, double hy, double hz) |
| 50 | +{ |
| 51 | + const double x0 = cx - hx, x1 = cx + hx; |
| 52 | + const double y0 = cy - hy, y1 = cy + hy; |
| 53 | + const double z0 = cz - hz, z1 = cz + hz; |
| 54 | + 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}}; |
| 55 | + // each quad is wound counter-clockwise seen from outside, so the facet normal points outward |
| 56 | + 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}}; |
| 57 | + for (const auto& face : quad) { |
| 58 | + shape.AddFacet(corner[face[0]], corner[face[1]], corner[face[2]]); |
| 59 | + shape.AddFacet(corner[face[0]], corner[face[2]], corner[face[3]]); |
| 60 | + } |
| 61 | +} |
| 62 | + |
| 63 | +/// The Moeller-Trumbore distance used by O2Tessellated's leaf test, repeated here as the oracle. |
| 64 | +double rayTriangleReference(const double* origin, const double* dir, const Vertex_t& v0, const Vertex_t& v1, |
| 65 | + const Vertex_t& v2) |
| 66 | +{ |
| 67 | + constexpr double EPS = 1.e-8; |
| 68 | + const double infinity = std::numeric_limits<double>::infinity(); |
| 69 | + const double e1[3] = {v1[0] - v0[0], v1[1] - v0[1], v1[2] - v0[2]}; |
| 70 | + const double e2[3] = {v2[0] - v0[0], v2[1] - v0[1], v2[2] - v0[2]}; |
| 71 | + const double p[3] = {dir[1] * e2[2] - dir[2] * e2[1], dir[2] * e2[0] - dir[0] * e2[2], |
| 72 | + dir[0] * e2[1] - dir[1] * e2[0]}; |
| 73 | + const double det = e1[0] * p[0] + e1[1] * p[1] + e1[2] * p[2]; |
| 74 | + if (std::abs(det) <= EPS) { |
| 75 | + return infinity; |
| 76 | + } |
| 77 | + const double tvec[3] = {origin[0] - v0[0], origin[1] - v0[1], origin[2] - v0[2]}; |
| 78 | + const double invDet = 1.0 / det; |
| 79 | + const double u = (tvec[0] * p[0] + tvec[1] * p[1] + tvec[2] * p[2]) * invDet; |
| 80 | + if (u < 0.0 || u > 1.0) { |
| 81 | + return infinity; |
| 82 | + } |
| 83 | + const double q[3] = {tvec[1] * e1[2] - tvec[2] * e1[1], tvec[2] * e1[0] - tvec[0] * e1[2], |
| 84 | + tvec[0] * e1[1] - tvec[1] * e1[0]}; |
| 85 | + const double v = (dir[0] * q[0] + dir[1] * q[1] + dir[2] * q[2]) * invDet; |
| 86 | + if (v < 0.0 || u + v > 1.0) { |
| 87 | + return infinity; |
| 88 | + } |
| 89 | + const double t = e2[0] * q[0] + e2[1] * q[1] + e2[2] * q[2]; |
| 90 | + return (t * invDet > 0.) ? t * invDet : infinity; |
| 91 | +} |
| 92 | + |
| 93 | +/// The unpruned answer: the nearest facet over every facet of the mesh, entering or exiting. |
| 94 | +double bruteForce(const O2Tessellated& shape, const double* origin, const double* dir, bool entering) |
| 95 | +{ |
| 96 | + double best = TGeoShape::Big(); |
| 97 | + for (int facet = 0; facet < shape.GetNfacets(); ++facet) { |
| 98 | + const auto& description = shape.GetFacet(facet); |
| 99 | + const Vertex_t& v0 = shape.GetVertex(description[0]); |
| 100 | + const Vertex_t& v1 = shape.GetVertex(description[1]); |
| 101 | + const Vertex_t& v2 = shape.GetVertex(description[2]); |
| 102 | + const double e1[3] = {v1[0] - v0[0], v1[1] - v0[1], v1[2] - v0[2]}; |
| 103 | + const double e2[3] = {v2[0] - v0[0], v2[1] - v0[1], v2[2] - v0[2]}; |
| 104 | + const double normal[3] = {e1[1] * e2[2] - e1[2] * e2[1], e1[2] * e2[0] - e1[0] * e2[2], |
| 105 | + e1[0] * e2[1] - e1[1] * e2[0]}; |
| 106 | + const double along = normal[0] * dir[0] + normal[1] * dir[1] + normal[2] * dir[2]; |
| 107 | + // the same facing filter the shape applies: entering facets face the ray, exiting ones face away |
| 108 | + if (entering ? (along > 0.) : (along <= 0.)) { |
| 109 | + continue; |
| 110 | + } |
| 111 | + best = std::min(best, rayTriangleReference(origin, dir, v0, v1, v2)); |
| 112 | + } |
| 113 | + return best; |
| 114 | +} |
| 115 | + |
| 116 | +/// Eight boxes in a row, so every axial ray meets sixteen facets and the BVH has many leaves. |
| 117 | +void buildRow(O2Tessellated& shape) |
| 118 | +{ |
| 119 | + for (int index = 0; index < 8; ++index) { |
| 120 | + addBox(shape, -21. + 6. * index, 0., 0., 2., 3., 4.); |
| 121 | + } |
| 122 | + shape.CloseShape(true, false, false); |
| 123 | +} |
| 124 | +} // namespace |
| 125 | + |
| 126 | +BOOST_AUTO_TEST_CASE(PrunedRayQueriesEqualTheBruteForceMinimum) |
| 127 | +{ |
| 128 | + O2Tessellated shape("row"); |
| 129 | + buildRow(shape); |
| 130 | + BOOST_CHECK_EQUAL(shape.GetNfacets(), 96); |
| 131 | + |
| 132 | + Rng rng(20260912); |
| 133 | + int outsideHits = 0; |
| 134 | + int insideHits = 0; |
| 135 | + for (int trial = 0; trial < 4000; ++trial) { |
| 136 | + // origins inside the row and well outside it, so both directions are exercised |
| 137 | + const double origin[3] = {rng.uniform(-40., 40.), rng.uniform(-12., 12.), rng.uniform(-12., 12.)}; |
| 138 | + double dir[3] = {rng.uniform(-1., 1.), rng.uniform(-1., 1.), rng.uniform(-1., 1.)}; |
| 139 | + const double norm = std::sqrt(dir[0] * dir[0] + dir[1] * dir[1] + dir[2] * dir[2]); |
| 140 | + if (norm < 1.e-6) { |
| 141 | + continue; |
| 142 | + } |
| 143 | + for (int index = 0; index < 3; ++index) { |
| 144 | + dir[index] /= norm; |
| 145 | + } |
| 146 | + |
| 147 | + const double outside = shape.DistFromOutside(origin, dir, 1, TGeoShape::Big(), nullptr); |
| 148 | + const double inside = shape.DistFromInside(origin, dir, 1, TGeoShape::Big(), nullptr); |
| 149 | + const double outsideReference = bruteForce(shape, origin, dir, true); |
| 150 | + const double insideReference = bruteForce(shape, origin, dir, false); |
| 151 | + |
| 152 | + BOOST_CHECK_EQUAL(outside, outsideReference); |
| 153 | + BOOST_CHECK_EQUAL(inside, insideReference); |
| 154 | + outsideHits += outsideReference < TGeoShape::Big() ? 1 : 0; |
| 155 | + insideHits += insideReference < TGeoShape::Big() ? 1 : 0; |
| 156 | + } |
| 157 | + // the case is only meaningful if the rays really hit the mesh; this sampling gives about 450 |
| 158 | + // entering and 3000 exiting hits |
| 159 | + BOOST_CHECK_GT(outsideHits, 200); |
| 160 | + BOOST_CHECK_GT(insideHits, 200); |
| 161 | +} |
| 162 | + |
| 163 | +BOOST_AUTO_TEST_CASE(APrunedRayFindsTheNearestOfManyFacetsAlongIt) |
| 164 | +{ |
| 165 | + O2Tessellated shape("row"); |
| 166 | + buildRow(shape); |
| 167 | + |
| 168 | + // straight down the row: eight boxes, so sixteen entering and sixteen exiting facets are in line |
| 169 | + const double origin[3] = {-40., 0., 0.}; |
| 170 | + const double dir[3] = {1., 0., 0.}; |
| 171 | + BOOST_CHECK_EQUAL(shape.DistFromOutside(origin, dir, 1, TGeoShape::Big(), nullptr), 17.); |
| 172 | + BOOST_CHECK_EQUAL(shape.DistFromOutside(origin, dir, 1, TGeoShape::Big(), nullptr), |
| 173 | + bruteForce(shape, origin, dir, true)); |
| 174 | + |
| 175 | + // from inside the first box, the exit is its own far face and not a later box's |
| 176 | + const double inner[3] = {-21., 0., 0.}; |
| 177 | + BOOST_CHECK_EQUAL(shape.DistFromInside(inner, dir, 1, TGeoShape::Big(), nullptr), 2.); |
| 178 | + BOOST_CHECK_EQUAL(shape.DistFromInside(inner, dir, 1, TGeoShape::Big(), nullptr), |
| 179 | + bruteForce(shape, inner, dir, false)); |
| 180 | +} |
0 commit comments