14#define BOOST_TEST_MODULE Test O2BVHSurfaceSolid class
15#define BOOST_TEST_MAIN
16#define BOOST_TEST_DYN_LINK
17#include <boost/test/unit_test.hpp>
23#include "../src/BoundedSurface.h"
27#include "TGeoBoolNode.h"
28#include "TGeoCompositeShape.h"
30#include "TGeoManager.h"
31#include "TGeoMaterial.h"
32#include "TGeoMatrix.h"
33#include "TGeoMedium.h"
36#include "TGeoSphere.h"
39#include "TGeoVolume.h"
48#include <initializer_list>
57using Point2D = SurfaceSolid::Point2D;
58using Point3D = SurfaceSolid::Point3D;
61std::vector<Point2D> rectangleWire(
double extentU,
double extentV)
63 return {{0., 0.}, {extentU, 0.}, {extentU, extentV}, {0., extentV}};
66using BoundaryCurve = SurfaceSolid::PlanarBoundaryCurve;
69std::vector<BoundaryCurve> circleWire(
double radius,
bool clockwise =
false)
76std::vector<BoundaryCurve> paramRectWire(
double uMin,
double uMax,
double vMin,
double vMax)
78 return {BoundaryCurve::makeLine({uMin, vMin}, {uMax, vMin}), BoundaryCurve::makeLine({uMax, vMin}, {uMax, vMax}),
79 BoundaryCurve::makeLine({uMax, vMax}, {uMin, vMax}), BoundaryCurve::makeLine({uMin, vMax}, {uMin, vMin})};
83std::vector<surf::Curve2D> paramRectWireCurves(
double uMin,
double uMax,
double vMin,
double vMax)
85 return {surf::Curve2D::makeLine({uMin, vMin}, {uMax, vMin}), surf::Curve2D::makeLine({uMax, vMin}, {uMax, vMax}),
86 surf::Curve2D::makeLine({uMax, vMax}, {uMin, vMax}), surf::Curve2D::makeLine({uMin, vMax}, {uMin, vMin})};
92 double radius,
double holeRadius = 0.)
94 std::vector<std::vector<BoundaryCurve>> inners;
95 if (holeRadius > 0.) {
96 inners.push_back(circleWire(holeRadius,
true));
98 return solid.AddCurvedPlanarSurface(
center, axisU, axisV, circleWire(radius), inners);
111FaceFrame boxFaceFrame(
int faceIndex,
double halfX,
double halfY,
double halfZ)
115 return {{halfX, -halfY, -halfZ}, {0., 1., 0.}, {0., 0., 1.}, 2. * halfY, 2. * halfZ};
117 return {{-halfX, -halfY, -halfZ}, {0., 0., 1.}, {0., 1., 0.}, 2. * halfZ, 2. * halfY};
119 return {{-halfX, halfY, -halfZ}, {0., 0., 1.}, {1., 0., 0.}, 2. * halfZ, 2. * halfX};
121 return {{-halfX, -halfY, -halfZ}, {1., 0., 0.}, {0., 0., 1.}, 2. * halfX, 2. * halfZ};
123 return {{-halfX, -halfY, halfZ}, {1., 0., 0.}, {0., 1., 0.}, 2. * halfX, 2. * halfY};
125 return {{-halfX, -halfY, -halfZ}, {0., 1., 0.}, {1., 0., 0.}, 2. * halfY, 2. * halfX};
132bool addBoxFace(SurfaceSolid& solid,
int faceIndex,
double halfX,
double halfY,
double halfZ,
bool reversed =
false,
135 FaceFrame frame = boxFaceFrame(
faceIndex, halfX, halfY, halfZ);
137 std::swap(frame.axisU, frame.axisV);
138 std::swap(frame.extentU, frame.extentV);
140 for (
int dimension = 0; dimension < 3; ++dimension) {
141 frame.origin[dimension] +=
center[dimension];
143 return solid.AddPlanarSurface(frame.origin, frame.axisU, frame.axisV, rectangleWire(frame.extentU, frame.extentV));
146void addBoxSurfaces(SurfaceSolid& solid,
double halfX,
double halfY,
double halfZ,
150 BOOST_REQUIRE(addBoxFace(solid,
faceIndex, halfX, halfY, halfZ,
false,
center));
161void checkClose(
double value,
double reference,
double tolerance = 1.e-9)
166std::array<double, 3> unitDirection(
double x,
double y,
double z)
174void compareContainsGrid(
const SurfaceSolid& solid,
const TGeoShape&
reference,
double extent,
int samples)
176 for (
int stepX = 0; stepX <
samples; ++stepX) {
177 for (
int stepY = 0; stepY <
samples; ++stepY) {
178 for (
int stepZ = 0; stepZ <
samples; ++stepZ) {
179 const double point[3] = {-extent + 2. * extent * (stepX + 0.517) /
samples,
180 -extent + 2. * extent * (stepY + 0.263) /
samples,
181 -extent + 2. * extent * (stepZ + 0.741) /
samples};
182 BOOST_TEST_CONTEXT(
"point = (" << point[0] <<
", " << point[1] <<
", " << point[2] <<
")")
192void compareDistance(
const SurfaceSolid& solid,
const TGeoShape&
reference,
const std::array<double, 3>& point,
193 const std::array<double, 3>& direction,
double tolerance = 1.e-9)
195 BOOST_TEST_CONTEXT(
"point = (" << point[0] <<
", " << point[1] <<
", " << point[2] <<
") direction = ("
196 << direction[0] <<
", " << direction[1] <<
", " << direction[2] <<
")")
198 const bool inside =
reference.Contains(point.data());
201 checkClose(solid.DistFromInside(point.data(), direction.data(), 3),
202 reference.DistFromInside(point.data(), direction.data(), 3), tolerance);
204 checkClose(solid.DistFromOutside(point.data(), direction.data(), 3),
205 reference.DistFromOutside(point.data(), direction.data(), 3), tolerance);
217std::unique_ptr<SurfaceSolid> makeBoxSolid(
const char*
name,
double halfX,
double halfY,
double halfZ)
219 auto solid = std::make_unique<SurfaceSolid>(
name);
220 addBoxSurfaces(*solid, halfX, halfY, halfZ);
226std::unique_ptr<SurfaceSolid> makeTubeSolid(
const char*
name,
double innerRadius,
double outerRadius,
229 auto solid = std::make_unique<SurfaceSolid>(
name);
230 BOOST_REQUIRE(solid->AddCylindricalSurface({0., 0., 0.}, {0., 0., 1.}, {1., 0., 0.}, outerRadius, -halfHeight,
232 if (innerRadius > 0.) {
233 BOOST_REQUIRE(solid->AddCylindricalSurface({0., 0., 0.}, {0., 0., 1.}, {1., 0., 0.}, innerRadius, -halfHeight,
236 BOOST_REQUIRE(addDiskSurface(*solid, {0., 0., halfHeight}, {1., 0., 0.}, {0., 1., 0.}, outerRadius, innerRadius));
237 BOOST_REQUIRE(addDiskSurface(*solid, {0., 0., -halfHeight}, {1., 0., 0.}, {0., -1., 0.}, outerRadius, innerRadius));
242std::unique_ptr<SurfaceSolid> makeConeSolid(
const char*
name,
double radiusAtBottom,
double radiusAtTop,
245 auto solid = std::make_unique<SurfaceSolid>(
name);
246 BOOST_REQUIRE(solid->AddConicalSurface({0., 0., 0.}, {0., 0., 1.}, {1., 0., 0.}, radiusAtBottom, radiusAtTop,
247 -halfHeight, halfHeight));
248 BOOST_REQUIRE(addDiskSurface(*solid, {0., 0., halfHeight}, {1., 0., 0.}, {0., 1., 0.}, radiusAtTop));
249 BOOST_REQUIRE(addDiskSurface(*solid, {0., 0., -halfHeight}, {1., 0., 0.}, {0., -1., 0.}, radiusAtBottom));
254std::unique_ptr<SurfaceSolid> makeSphereSolid(
const char*
name,
double radius)
256 auto solid = std::make_unique<SurfaceSolid>(
name);
257 BOOST_REQUIRE(solid->AddSphericalSurface({0., 0., 0.}, {0., 0., 1.}, {1., 0., 0.}, radius));
262std::unique_ptr<SurfaceSolid> makeTorusSolid(
const char*
name,
double majorRadius,
double minorRadius)
264 auto solid = std::make_unique<SurfaceSolid>(
name);
265 BOOST_REQUIRE(solid->AddToroidalSurface({0., 0., 0.}, {0., 0., 1.}, {1., 0., 0.}, majorRadius, minorRadius));
272std::unique_ptr<SurfaceSolid> makeCapsuleSolid(
const char*
name,
double radius,
double halfHeight)
274 auto solid = std::make_unique<SurfaceSolid>(
name);
275 BOOST_REQUIRE(solid->AddCylindricalSurface({0., 0., 0.}, {0., 0., 1.}, {1., 0., 0.}, radius, -halfHeight,
277 BOOST_REQUIRE(solid->AddSphericalSurface({0., 0., halfHeight}, {0., 0., 1.}, {1., 0., 0.}, radius, 0.,
279 BOOST_REQUIRE(solid->AddSphericalSurface({0., 0., -halfHeight}, {0., 0., -1.}, {1., 0., 0.}, radius, 0.,
289const std::vector<std::array<double, 3>>& probeDirections()
291 static const std::vector<std::array<double, 3>> directions{
292 {1., 0., 0.}, {-1., 0., 0.}, {0., 1., 0.}, {0., -1., 0.}, {0., 0., 1.}, {0., 0., -1.}, unitDirection(1., 1., 0.), unitDirection(1., 0., 1.), unitDirection(0., 1., 1.), unitDirection(1., 1., 1.), unitDirection(-1., 1., -1.), unitDirection(0.37, -0.82, 0.44), unitDirection(-0.91, 0.13, 0.39), unitDirection(0.21, 0.55, -0.81)};
299std::vector<std::array<double, 3>> probeGrid(
double extent,
int samples)
301 std::vector<std::array<double, 3>> points;
303 for (
int stepX = 0; stepX <
samples; ++stepX) {
304 for (
int stepY = 0; stepY <
samples; ++stepY) {
305 for (
int stepZ = 0; stepZ <
samples; ++stepZ) {
306 points.push_back({-extent + 2. * extent * (stepX + 0.517) /
samples,
307 -extent + 2. * extent * (stepY + 0.263) /
samples,
308 -extent + 2. * extent * (stepZ + 0.741) /
samples});
323void checkDistanceAgainstLoop(
const SurfaceSolid& solid,
const std::array<double, 3>& point,
324 const std::array<double, 3>& direction,
double stepmax = TGeoShape::Big())
326 const double loopOutside = solid.DistFromOutside_Loop(point.data(), direction.data(), stepmax);
327 const double loopInside = solid.DistFromInside_Loop(point.data(), direction.data(), stepmax);
328 for (
const bool pruning : {
true,
false}) {
329 SurfaceSolid::SetRayTMaxPruning(pruning);
330 BOOST_TEST_CONTEXT(
"point = (" << point[0] <<
", " << point[1] <<
", " << point[2] <<
") direction = ("
331 << direction[0] <<
", " << direction[1] <<
", " << direction[2]
332 <<
") stepmax = " << stepmax <<
" pruning = " << pruning)
334 BOOST_CHECK_EQUAL(solid.DistFromOutside(point.data(), direction.data(), 3, stepmax), loopOutside);
335 BOOST_CHECK_EQUAL(solid.DistFromInside(point.data(), direction.data(), 3, stepmax), loopInside);
338 SurfaceSolid::SetRayTMaxPruning(
true);
342void sweepDistanceAgainstLoop(
const SurfaceSolid& solid,
double extent,
int samples)
344 for (
const auto& point : probeGrid(extent,
samples)) {
345 for (
const auto& direction : probeDirections()) {
346 checkDistanceAgainstLoop(solid, point, direction);
355void sweepDistanceAgainstReference(
const SurfaceSolid& solid,
const TGeoShape&
reference,
double extent,
int samples,
356 double tolerance = 1.e-9,
double skin = 1.e-6)
358 for (
const auto& point : probeGrid(extent,
samples)) {
359 const bool inside =
reference.Contains(point.data());
360 if (
reference.Safety(point.data(), inside) < skin) {
363 for (
const auto& direction : probeDirections()) {
364 BOOST_TEST_CONTEXT(
"point = (" << point[0] <<
", " << point[1] <<
", " << point[2] <<
") direction = ("
365 << direction[0] <<
", " << direction[1] <<
", " << direction[2] <<
")")
368 checkClose(solid.DistFromInside(point.data(), direction.data(), 3),
369 reference.DistFromInside(point.data(), direction.data(), 3), tolerance);
371 checkClose(solid.DistFromOutside(point.data(), direction.data(), 3),
372 reference.DistFromOutside(point.data(), direction.data(), 3), tolerance);
383 constexpr double halfX = 1.;
384 constexpr double halfY = 2.;
385 constexpr double halfZ = 3.;
387 SurfaceSolid solid(
"planarBox");
388 addBoxSurfaces(solid, halfX, halfY, halfZ);
396 int meshVertices = 0;
397 int meshSegments = 0;
398 int meshPolygons = 0;
399 solid.GetMeshNumbers(meshVertices, meshSegments, meshPolygons);
404 const std::array<std::array<double, 3>, 5> insidePoints{{{0., 0., 0.}, {0.9, 0., 0.}, {1., 0., 0.}, {1., 2., 3.}, {-1., -2., -3.}}};
405 for (
const auto& point : insidePoints) {
406 BOOST_TEST_CONTEXT(
"point = (" << point[0] <<
", " << point[1] <<
", " << point[2] <<
")")
413 const std::array<std::array<double, 3>, 4> outsidePoints{{{1.1, 0., 0.}, {0., 2.1, 0.}, {0., 0., -3.1}, {2., 3., 4.}}};
414 for (
const auto& point : outsidePoints) {
415 BOOST_TEST_CONTEXT(
"point = (" << point[0] <<
", " << point[1] <<
", " << point[2] <<
")")
422 const double fromLeft[3] = {-3., 0., 0.};
423 const double toRight[3] = {1., 0., 0.};
424 checkClose(solid.DistFromOutside(fromLeft, toRight, 3),
reference.DistFromOutside(fromLeft, toRight, 3));
426 const double fromFront[3] = {0., -5., 0.};
427 const double toBack[3] = {0., 1., 0.};
428 checkClose(solid.DistFromOutside(fromFront, toBack, 3),
reference.DistFromOutside(fromFront, toBack, 3));
430 const double fromCenter[3] = {0., 0., 0.};
431 const double alongX[3] = {1., 0., 0.};
432 const double alongZ[3] = {0., 0., -1.};
433 checkClose(solid.DistFromInside(fromCenter, alongX, 3),
reference.DistFromInside(fromCenter, alongX, 3));
434 checkClose(solid.DistFromInside(fromCenter, alongZ, 3),
reference.DistFromInside(fromCenter, alongZ, 3));
437 const double outsideSafetyPoint[3] = {2.5, 0., 0.};
438 checkClose(solid.Safety(fromCenter, kTRUE), halfX);
439 checkClose(solid.Safety(outsideSafetyPoint, kFALSE), 2.5 - halfX);
441 const double normalPoint[3] = {halfX, 0., 0.};
442 double normal[3] = {0., 0., 0.};
443 solid.ComputeNormal(normalPoint, alongX, normal);
444 checkClose(normal[0], 1.);
445 checkClose(normal[1], 0.);
446 checkClose(normal[2], 0.);
448 checkClose(solid.Capacity(), 8. * halfX * halfY * halfZ);
460 WireStatus reversedStatus = WireStatus::Valid;
461 auto reversedOuter = makeWire({{0., 0.}, {0., 1.}, {1., 1.}, {1., 0.}}, WireRole::Outer, reversedStatus);
462 BOOST_CHECK(reversedStatus == WireStatus::Reversed);
463 BOOST_CHECK_GT(reversedOuter.signedArea(), 0.);
466 WireStatus outerStatus = WireStatus::Valid;
467 auto outerWire = makeWire({{0., 0.}, {1., 0.}, {1., 1.}, {0., 1.}}, WireRole::Outer, outerStatus);
469 BOOST_CHECK_GT(outerWire.signedArea(), 0.);
472 WireStatus innerStatus = WireStatus::Valid;
473 auto innerWire = makeWire({{0., 0.}, {1., 0.}, {1., 1.}, {0., 1.}}, WireRole::Inner, innerStatus);
475 BOOST_CHECK_LT(innerWire.signedArea(), 0.);
479 WireStatus status = WireStatus::Valid;
483 BOOST_CHECK(!scratch.initialize({{0., 0.}, {1., 0.}, {2., 0.}}, WireRole::Outer, status));
487 BOOST_CHECK(!scratch.initialize({{0., 0.}, {1., 0.}, {0., 0.}, {1., 1.}}, WireRole::Outer, status));
488 BOOST_CHECK(status == WireStatus::DegenerateVertex);
491 const std::vector<surf::SurfaceEdge> openEdges{{{0., 0.}, {1., 0.}}, {{1., 0.}, {1., 1.}}, {{1., 1.}, {0.5, 0.5}}};
492 BOOST_CHECK(!scratch.initializeFromEdges(openEdges, WireRole::Outer, status));
496 const std::vector<surf::SurfaceEdge> closedEdges{
497 {{0., 0.}, {1., 0.}}, {{1., 0.}, {1., 1.}}, {{1., 1.}, {0., 1.}}, {{0., 1.}, {0., 0.}}};
498 BOOST_CHECK(scratch.initializeFromEdges(closedEdges, WireRole::Outer, status));
501 BOOST_CHECK(outerWire.classify({0.5, 0.5}) == surf::WireClassification::Inside);
502 BOOST_CHECK(outerWire.classify({1.5, 0.5}) == surf::WireClassification::Outside);
503 BOOST_CHECK(outerWire.classify({0.5, 0.}) == surf::WireClassification::Boundary);
513 : mVertices{firstVertex, secondVertex, thirdVertex}
515 mNormal =
normalized(
cross(secondVertex - firstVertex, thirdVertex - firstVertex));
520 for (
const auto&
vertex : mVertices) {
536 double bestDistanceSq = std::numeric_limits<double>::infinity();
537 for (
int vertexIndex = 0; vertexIndex < 3; ++vertexIndex) {
539 mVertices[(vertexIndex + 1) % 3]));
541 return bestDistanceSq;
558 void appendDisplayMesh(std::vector<Vec3>& vertices, std::vector<std::array<int, 3>>& triangles)
const override
560 const int firstVertexIndex =
static_cast<int>(vertices.size());
561 for (
const auto&
vertex : mVertices) {
562 vertices.push_back(
vertex);
564 triangles.push_back({firstVertexIndex, firstVertexIndex + 1, firstVertexIndex + 2});
569 for (
int vertexIndex = 0; vertexIndex < 3; ++vertexIndex) {
570 edges.emplace_back(mVertices[vertexIndex], mVertices[(vertexIndex + 1) % 3]);
575 std::array<Vec3, 3> mVertices;
582 auto dummy = std::make_unique<surf::DummyBoundedSurface>(
surf::Vec3{0., 0., 0.},
surf::Vec3{1., 0., 0.},
588 checkClose(
lower.xCoord, 0.);
589 checkClose(
upper.xCoord, 1.);
590 checkClose(
upper.yCoord, 1.);
592 const surf::Vec3 normal = dummy->normalAt({0., 0., 0.});
593 checkClose(std::abs(normal.
zCoord), 1.);
596 std::vector<surf::Vec3> vertices;
597 std::vector<std::array<int, 3>> triangles;
598 dummy->appendDisplayMesh(vertices, triangles);
603 std::vector<std::unique_ptr<surf::BoundedSurface>> surfaces;
604 surfaces.emplace_back(std::move(dummy));
612 constexpr double halfX = 1.;
613 constexpr double halfY = 1.5;
614 constexpr double halfZ = 2.;
617 SurfaceSolid missing(
"missingFaceBox");
619 BOOST_REQUIRE(addBoxFace(missing,
faceIndex, halfX, halfY, halfZ));
621 missing.CloseShape(
false);
625 SurfaceSolid reversed(
"reversedFaceBox");
626 BOOST_REQUIRE(addBoxFace(reversed, 0, halfX, halfY, halfZ,
true));
628 BOOST_REQUIRE(addBoxFace(reversed,
faceIndex, halfX, halfY, halfZ));
630 reversed.CloseShape(
false);
641 using Reliability = SurfaceSolid::NavigationReliability;
642 constexpr double halfX = 1.;
643 constexpr double halfY = 1.5;
644 constexpr double halfZ = 2.;
647 SurfaceSolid fresh(
"freshBox");
649 BOOST_REQUIRE(addBoxFace(fresh,
faceIndex, halfX, halfY, halfZ));
651 BOOST_CHECK(fresh.GetNavigationReliability() == Reliability::Undetermined);
654 fresh.CloseShape(
false);
655 BOOST_CHECK(fresh.GetNavigationReliability() == Reliability::Reliable);
662 SurfaceSolid missing(
"missingFaceBoxState");
664 BOOST_REQUIRE(addBoxFace(missing,
faceIndex, halfX, halfY, halfZ));
666 missing.CloseShape(
false);
667 BOOST_CHECK(missing.GetNavigationReliability() == Reliability::OpenSurfaceSet);
672 SurfaceSolid reversed(
"reversedFaceBoxState");
673 BOOST_REQUIRE(addBoxFace(reversed, 0, halfX, halfY, halfZ,
true));
675 BOOST_REQUIRE(addBoxFace(reversed,
faceIndex, halfX, halfY, halfZ));
677 reversed.CloseShape(
false);
678 BOOST_CHECK(reversed.GetNavigationReliability() == Reliability::ReversedFaces);
684 SurfaceSolid duplicated(
"duplicatedFaceBox");
685 for (
int pass = 0; pass < 2; ++pass) {
687 BOOST_REQUIRE(addBoxFace(duplicated,
faceIndex, halfX, halfY, halfZ));
690 duplicated.CloseShape(
false);
691 BOOST_CHECK(duplicated.GetNavigationReliability() == Reliability::NonManifold);
693 BOOST_CHECK(duplicated.GetNonManifoldEdgeCount() > 0);
695 BOOST_CHECK_EQUAL(std::string(SurfaceSolid::GetNavigationReliabilityName(Reliability::Reliable)),
"reliable");
696 BOOST_CHECK_EQUAL(std::string(SurfaceSolid::GetNavigationReliabilityName(Reliability::OpenSurfaceSet)),
698 BOOST_CHECK_EQUAL(std::string(SurfaceSolid::GetNavigationReliabilityName(Reliability::NonManifold)),
"non-manifold");
708 WireStatus status = WireStatus::Valid;
709 auto square = makeWire({{0., 0.}, {1., 0.}, {1., 1.}, {0., 1.}}, WireRole::Outer, status);
710 BOOST_REQUIRE(status == WireStatus::Valid);
713 BOOST_CHECK(square.classify({0.5, 0.5 * surf::kTolerance}) == WireClassification::Boundary);
714 BOOST_CHECK(square.classify({0.5, -0.5 * surf::kTolerance}) == WireClassification::Boundary);
716 BOOST_CHECK(square.classify({0.5, 1.e3 * surf::kTolerance}) == WireClassification::Inside);
717 BOOST_CHECK(square.classify({0.5, -1.e3 * surf::kTolerance}) == WireClassification::Outside);
721 std::string planeError;
722 BOOST_REQUIRE(plane.
initialize({0., 0., 0.}, {1., 0., 0.}, {0., 1., 0.},
723 {{0., 0.}, {1., 0.}, {1., 1.}, {0., 1.}}, {}, planeError));
726 std::vector<surf::RayHit> hits;
730 plane.appendIntersections(
origin, grazing, 0., 1.e30, hits);
735 plane.appendIntersections(
origin, steep, 0., 1.e30, hits);
736 BOOST_REQUIRE_EQUAL(hits.size(), 1u);
741 plane.appendIntersections(
origin, steep, 2., 1.e30, hits);
756 WireStatus status = WireStatus::Valid;
757 auto square = makeWire({{0., 0.}, {2., 0.}, {2., 3.}, {0., 3.}}, WireRole::Outer, status);
758 BOOST_REQUIRE(status == WireStatus::Valid);
759 checkClose(square.signedArea(), 6.);
764 checkClose(
lower.uCoord, 0.);
765 checkClose(
lower.vCoord, 0.);
766 checkClose(
upper.uCoord, 2.);
767 checkClose(
upper.vCoord, 3.);
770 const auto samples = square.sampledBoundary();
777 double parameter = -1.;
779 checkClose(projected.
uCoord, 1.);
780 checkClose(projected.
vCoord, 0.);
781 checkClose(parameter, 0.5);
783 bottom.closestPoint({-5., 1.}, parameter);
784 checkClose(parameter, 0.);
785 bottom.closestPoint({5., 1.}, parameter);
786 checkClose(parameter, 1.);
787 checkClose(std::sqrt(
bottom.distanceSq({1., 4.})), 4.);
790 WireStatus reversedStatus = WireStatus::Valid;
791 auto reversed = makeWire({{0., 0.}, {0., 3.}, {2., 3.}, {2., 0.}}, WireRole::Outer, reversedStatus);
792 BOOST_CHECK(reversedStatus == WireStatus::Reversed);
793 checkClose(reversed.signedArea(), 6.);
796 reversed.parametricBounds(reversedLower, reversedUpper);
797 checkClose(reversedUpper.uCoord, 2.);
798 checkClose(reversedUpper.vCoord, 3.);
802 const std::vector<surf::SurfaceEdge> openEdges{{{0., 0.}, {2., 0.}}, {{2., 0.}, {2., 3.}}, {{2., 3.}, {1., 1.}}};
807 BOOST_CHECK(square.classify({1., 0.}) == WireClassification::Boundary);
808 BOOST_CHECK(square.classify({1., 1.5}) == WireClassification::Inside);
809 BOOST_CHECK(square.classify({3., 1.5}) == WireClassification::Outside);
813 std::string faceError;
814 const std::vector<surf::Vec2> outer{{0., 0.}, {4., 0.}, {4., 4.}, {0., 4.}};
815 const std::vector<std::vector<surf::Vec2>> holes{{{1., 1.}, {3., 1.}, {3., 3.}, {1., 3.}}};
816 BOOST_REQUIRE(holedFace.
initialize({0., 0., 0.}, {1., 0., 0.}, {0., 1., 0.}, outer, holes, faceError));
818 bool boundary =
false;
826 checkClose(holedFace.
area(), 12.);
838 const Curve2D line = Curve2D::makeLine({0., 0.}, {4., 0.});
839 checkClose(line.startPoint().uCoord, 0.);
840 checkClose(line.endPoint().uCoord, 4.);
841 const surf::Vec2 lineTangent = line.tangentAt(0.5);
842 checkClose(lineTangent.
uCoord, 1.);
843 checkClose(lineTangent.
vCoord, 0.);
845 double lineParameter = -1.;
846 const surf::Vec2 lineProjection = line.closestPoint({1., 5.}, lineParameter);
847 checkClose(lineProjection.
uCoord, 1.);
848 checkClose(lineProjection.
vCoord, 0.);
849 checkClose(lineParameter, 0.25);
850 checkClose(std::sqrt(line.distanceSq({1., 5.})), 5.);
854 const Curve2D quarter = Curve2D::makeArc({0., 0.}, 2., 0.,
surf::kHalfPi);
855 checkClose(quarter.startPoint().uCoord, 2.);
856 checkClose(quarter.startPoint().vCoord, 0.);
857 checkClose(quarter.endPoint().uCoord, 0.);
858 checkClose(quarter.endPoint().vCoord, 2.);
860 const surf::Vec2 arcTangent = quarter.tangentAt(0.);
861 checkClose(arcTangent.
uCoord, 0.);
862 checkClose(arcTangent.
vCoord, 1.);
867 quarter.extendBounds(arcLower, arcUpper);
868 checkClose(arcLower.uCoord, 0.);
869 checkClose(arcLower.vCoord, 0.);
870 checkClose(arcUpper.uCoord, 2.);
871 checkClose(arcUpper.vCoord, 2.);
874 double arcParameter = -1.;
875 const surf::Vec2 arcProjection = quarter.closestPoint({5., 5.}, arcParameter);
876 checkClose(std::hypot(arcProjection.
uCoord, arcProjection.
vCoord), 2.);
877 checkClose(arcParameter, 0.5);
880 const Curve2D circle = Curve2D::makeCircle({1., -1.}, 3.);
883 circle.extendBounds(circleLower, circleUpper);
884 checkClose(circleLower.uCoord, -2.);
885 checkClose(circleUpper.uCoord, 4.);
886 checkClose(circleLower.vCoord, -4.);
887 checkClose(circleUpper.vCoord, 2.);
890 WireStatus status = WireStatus::Valid;
892 BOOST_REQUIRE(disk.initialize({Curve2D::makeCircle({0., 0.}, 2.)}, WireRole::Outer, status));
895 checkClose(disk.signedArea(),
surf::kPi * 4., 1.e-9);
896 BOOST_CHECK(disk.classify({0., 0.}) == WireClassification::Inside);
897 BOOST_CHECK(disk.classify({1.5, 0.}) == WireClassification::Inside);
898 BOOST_CHECK(disk.classify({3., 0.}) == WireClassification::Outside);
899 BOOST_CHECK(disk.classify({0., 3.}) == WireClassification::Outside);
900 BOOST_CHECK(disk.classify({2., 0.}) == WireClassification::Boundary);
901 BOOST_CHECK(disk.classify({0., -2.}) == WireClassification::Boundary);
904 WireStatus reversedStatus = WireStatus::Valid;
905 CurveWire reversedDisk;
906 BOOST_REQUIRE(reversedDisk.initialize({Curve2D::makeCircle({0., 0.}, 2.,
true)}, WireRole::Outer, reversedStatus));
907 BOOST_CHECK(reversedStatus == WireStatus::Reversed);
908 checkClose(reversedDisk.signedArea(),
surf::kPi * 4., 1.e-9);
915 BOOST_REQUIRE(outerRing.initialize({Curve2D::makeCircle({0., 0.}, 3.)}, WireRole::Outer, outerStatus));
916 BOOST_REQUIRE(innerRing.initialize({Curve2D::makeCircle({0., 0.}, 1.)}, WireRole::Inner, holeStatus));
918 BOOST_CHECK_LT(innerRing.signedArea(), 0.);
920 checkClose(outerRing.signedArea() + innerRing.signedArea(),
surf::kPi * (9. - 1.), 1.e-9);
924 BOOST_CHECK(outerRing.classify(materialPoint) == WireClassification::Inside);
925 BOOST_CHECK(innerRing.classify(materialPoint) == WireClassification::Outside);
928 BOOST_CHECK(outerRing.classify(holePoint) == WireClassification::Inside);
929 BOOST_CHECK(innerRing.classify(holePoint) == WireClassification::Inside);
935 const std::vector<Curve2D> halfDiskCurves{Curve2D::makeLine({-2., 0.}, {2., 0.}),
936 Curve2D::makeArc({0., 0.}, 2., 0.,
surf::kPi)};
937 BOOST_REQUIRE(halfDisk.initialize(halfDiskCurves, WireRole::Outer, halfStatus));
939 checkClose(halfDisk.signedArea(), 0.5 *
surf::kPi * 4., 1.e-9);
940 BOOST_CHECK(halfDisk.classify({0., 1.}) == WireClassification::Inside);
941 BOOST_CHECK(halfDisk.classify({0., -1.}) == WireClassification::Outside);
942 BOOST_CHECK(halfDisk.classify({0., 0.}) == WireClassification::Boundary);
947 const std::vector<Curve2D> openCurves{Curve2D::makeLine({0., 0.}, {2., 0.}),
948 Curve2D::makeLine({2., 0.}, {2., 2.})};
949 BOOST_CHECK(!openWire.initialize(openCurves, WireRole::Outer, openStatus));
963 const std::vector<Vec2> poles{{0., 0.}, {1., 2.}, {2., -1.}, {3., 1.}, {4., 0.}};
964 const std::vector<double> knots{0., 0., 0., 0., 0.5, 1., 1., 1., 1.};
965 const Curve2D spline = Curve2D::makeBSpline(3, poles, {}, knots);
967 checkClose(spline.startPoint().uCoord, 0.);
968 checkClose(spline.startPoint().vCoord, 0.);
969 checkClose(spline.endPoint().uCoord, 4.);
970 checkClose(spline.endPoint().vCoord, 0.);
972 Vec2
lower{1.e30, 1.e30};
973 Vec2
upper{-1.e30, -1.e30};
975 checkClose(
lower.uCoord, 0.);
976 checkClose(
upper.uCoord, 4.);
977 checkClose(
lower.vCoord, -1.);
978 checkClose(
upper.vCoord, 2.);
981 const std::vector<Vec2> circlePoles{{1., 0.}, {1., 1.}, {0., 1.}};
982 const std::vector<double> circleWeights{1., std::sqrt(0.5), 1.};
983 const std::vector<double> circleKnots{0., 0., 0., 1., 1., 1.};
984 const Curve2D quarter = Curve2D::makeBSpline(2, circlePoles, circleWeights, circleKnots);
987 const Vec2 point = quarter.pointAt(
static_cast<double>(
index) / 8);
988 checkClose(std::hypot(point.uCoord, point.vCoord), 1., 1.e-9);
992 const std::vector<Curve2D> loop{spline, Curve2D::makeLine({4., 0.}, {4., -3.}),
993 Curve2D::makeLine({4., -3.}, {0., -3.}), Curve2D::makeLine({0., -3.}, {0., 0.})};
994 WireStatus status = WireStatus::Valid;
996 BOOST_REQUIRE(wire.initialize(loop, WireRole::Outer, status));
997 double referenceArea = 0.;
998 const auto boundarySamples = wire.sampledBoundary();
999 for (
size_t k = 0; k + 1 < boundarySamples.size(); ++k) {
1000 referenceArea += 0.5 * (boundarySamples[k].uCoord * boundarySamples[k + 1].vCoord -
1001 boundarySamples[k + 1].uCoord * boundarySamples[k].vCoord);
1005 checkClose(wire.signedArea(), std::abs(referenceArea), 1.e-4);
1008 BOOST_CHECK(wire.classify({2., -1.5}) == WireClassification::Inside);
1009 BOOST_CHECK(wire.classify({2., -2.9}) == WireClassification::Inside);
1010 BOOST_CHECK(wire.classify({-1., -1.}) == WireClassification::Outside);
1011 BOOST_CHECK(wire.classify({2., 5.}) == WireClassification::Outside);
1012 BOOST_CHECK(wire.classify({0., 0.}) == WireClassification::Boundary);
1013 BOOST_CHECK(wire.classify({2., -3.}) == WireClassification::Boundary);
1017 const Curve2D arch = Curve2D::makeBSpline(2, {{0., 0.}, {1., 2.}, {2., 0.}}, {}, {0., 0., 0., 1., 1., 1.});
1018 checkClose(arch.pointAt(0.5).vCoord, 1.);
1019 WireStatus archStatus = WireStatus::Valid;
1020 CurveWire archRegion;
1021 BOOST_REQUIRE(archRegion.initialize({arch, Curve2D::makeLine({2., 0.}, {0., 0.})}, WireRole::Outer, archStatus));
1022 BOOST_CHECK(archRegion.classify({1., 0.5}) == WireClassification::Inside);
1023 BOOST_CHECK(archRegion.classify({1., 1.5}) == WireClassification::Outside);
1026 BOOST_CHECK(archRegion.classify({-1., 1.}) == WireClassification::Outside);
1027 BOOST_CHECK(archRegion.classify({3., 1.}) == WireClassification::Outside);
1030 Curve2D reversed = spline;
1031 reversed.reverseInPlace();
1032 checkClose(reversed.startPoint().uCoord, 4.);
1033 checkClose(reversed.endPoint().uCoord, 0.);
1034 checkClose(reversed.pointAt(0.25).uCoord, spline.pointAt(0.75).uCoord, 1.e-9);
1035 checkClose(reversed.pointAt(0.25).vCoord, spline.pointAt(0.75).vCoord, 1.e-9);
1052 const Curve2D sCurve =
1053 Curve2D::makeBSpline(3, {{0., 0.}, {0.5, 1.}, {1.5, -1.}, {2., 0.}}, {}, {0., 0., 0., 0., 1., 1., 1., 1.});
1054 BOOST_REQUIRE(sCurve.valid());
1055 checkClose(sCurve.pointAt(0.5).uCoord, 1., 1.e-12);
1056 checkClose(sCurve.pointAt(0.5).vCoord, 0., 1.e-12);
1057 double worstOffChord = 0.;
1058 for (
int step = 0; step <= 64; ++step) {
1059 worstOffChord = std::max(worstOffChord, std::abs(sCurve.pointAt(
static_cast<double>(step) / 64).vCoord));
1064 sCurve.bsplineSampleInto(
samples);
1066 double worstSampleError = 0.;
1067 for (
int step = 0; step <= 64; ++step) {
1068 const Vec2 onCurve = sCurve.pointAt(
static_cast<double>(step) / 64);
1069 double nearest = 1.e30;
1073 worstSampleError = std::max(worstSampleError, std::sqrt(nearest));
1082 std::vector<Vec2> straightPoles;
1083 std::vector<double> uniformKnots{0., 0., 0., 0.};
1084 constexpr int spanCount = 8;
1086 straightPoles.push_back({
static_cast<double>(
index), 0.});
1089 uniformKnots.push_back(
static_cast<double>(
index) / spanCount);
1091 uniformKnots.insert(uniformKnots.end(), {1., 1., 1., 1.});
1092 const Curve2D straight = Curve2D::makeBSpline(3, straightPoles, {}, uniformKnots);
1093 BOOST_REQUIRE(straight.valid());
1094 std::vector<Vec2> straightSamples;
1095 straight.bsplineSampleInto(straightSamples);
1096 BOOST_CHECK(
static_cast<int>(straightSamples.size()) >= spanCount + 1);
1106 BOOST_REQUIRE(annulus.
initialize({0., 0., 0.}, {1., 0., 0.}, {0., 1., 0.},
1107 {Curve2D::makeCircle({0., 0.}, 2.)},
1108 {{Curve2D::makeCircle({0., 0.}, 1.,
true)}},
error));
1114 {Curve2D::makeCircle({0., 0.}, 2.)}, {},
error));
1117 BOOST_CHECK(annulus.containsPointOnSurface({1.5, 0., 0.}));
1118 BOOST_CHECK(!annulus.containsPointOnSurface({0.5, 0., 0.}));
1119 BOOST_CHECK(!annulus.containsPointOnSurface({3., 0., 0.}));
1120 BOOST_CHECK(!annulus.containsPointOnSurface({1.5, 0., 0.5}));
1123 std::vector<surf::RayHit> hits;
1124 annulus.appendIntersections({1.5, 0., 1.}, {0., 0., -1.}, 0., 1.e30, hits);
1125 BOOST_REQUIRE_EQUAL(hits.size(), 1u);
1126 checkClose(hits.front().distance, 1.);
1127 checkClose(hits.front().normal.zCoord, 1.);
1129 annulus.appendIntersections({0.5, 0., 1.}, {0., 0., -1.}, 0., 1.e30, hits);
1134 checkClose(annulus.distanceSqToPatch({0., 0., 0.}), 1.);
1135 checkClose(annulus.distanceSqToPatch({1.5, 0., 2.}), 4.);
1136 checkClose(annulus.distanceSqToPatch({4., 0., 0.}), 4.);
1137 checkClose(annulus.distanceSqToPatch({0.5, 0., 1.}), 1.25);
1141 BOOST_REQUIRE(offsetDisk.
initialize({0., 0., 2.}, {1., 0., 0.}, {0., 1., 0.},
1142 {Curve2D::makeCircle({0., 0.}, 1.5)}, {},
error));
1143 checkClose(offsetDisk.capacityContribution(), 2. *
surf::kPi * 1.5 * 1.5 / 3., 1.e-9);
1152 BOOST_REQUIRE(cylinder.
initialize({0., 0., 0.}, {0., 0., 1.}, {1., 0., 0.}, 2., -3., 3., 0.,
surf::kTwoPi,
1156 std::vector<surf::RayHit> hits;
1158 BOOST_REQUIRE_EQUAL(hits.size(), 2u);
1159 const double chordHalf = std::sqrt(4. - 0.25);
1160 checkClose(hits[0].
distance, 5. - chordHalf);
1161 checkClose(hits[1].
distance, 5. + chordHalf);
1163 BOOST_CHECK_LT(hits[0].normal.xCoord, 0.);
1164 BOOST_CHECK_GT(hits[1].normal.xCoord, 0.);
1183 BOOST_REQUIRE(halfCylinder.
initialize({0., 0., 0.}, {0., 0., 1.}, {1., 0., 0.}, 2., -3., 3., 0.,
surf::kPi,
1189 BOOST_REQUIRE_EQUAL(hits.size(), 1u);
1190 checkClose(hits.front().distance, 7.);
1195 constexpr double radius = 2.;
1196 constexpr double halfHeight = 3.;
1198 SurfaceSolid solid(
"closedCylinder");
1199 BOOST_REQUIRE(solid.AddCylindricalSurface({0., 0., 0.}, {0., 0., 1.}, {1., 0., 0.}, radius, -halfHeight,
1202 BOOST_REQUIRE(addDiskSurface(solid, {0., 0., halfHeight}, {1., 0., 0.}, {0., 1., 0.}, radius));
1203 BOOST_REQUIRE(addDiskSurface(solid, {0., 0., -halfHeight}, {1., 0., 0.}, {0., -1., 0.}, radius));
1209 TGeoTube
reference(
"referenceTube", 0., radius, halfHeight);
1210 compareContainsGrid(solid,
reference, 4., 9);
1212 compareDistance(solid,
reference, {5., 0., 0.}, {-1., 0., 0.});
1213 compareDistance(solid,
reference, {0., 0., 5.}, {0., 0., -1.});
1214 compareDistance(solid,
reference, {5., 0.5, 1.}, {-1., 0., 0.});
1215 compareDistance(solid,
reference, {-4., -1., -2.}, unitDirection(1., 0.3, 0.5));
1216 compareDistance(solid,
reference, {0., 0., 0.}, {1., 0., 0.});
1217 compareDistance(solid,
reference, {0., 0., 0.}, {0., 0., 1.});
1218 compareDistance(solid,
reference, {0., 0., 0.}, unitDirection(1., 1., 1.));
1219 compareDistance(solid,
reference, {1., 0.5, -2.}, unitDirection(0.3, -0.4, 0.5));
1220 compareDistance(solid,
reference, {5., 2.5, 0.}, {-1., 0., 0.});
1224 const double center[3] = {0., 0., 0.};
1225 const double insidePoint[3] = {1., 0.5, 1.};
1226 const double radialOutside[3] = {4., 0., 0.};
1227 const double axialOutside[3] = {0., 0., 5.};
1228 const double cornerOutside[3] = {4., 0., 5.};
1229 checkClose(solid.Safety(
center, kTRUE), radius);
1230 checkClose(solid.Safety(insidePoint, kTRUE), radius - std::sqrt(1.25));
1231 checkClose(solid.Safety(radialOutside, kFALSE), 2.);
1232 checkClose(solid.Safety(axialOutside, kFALSE), 2.);
1233 checkClose(solid.Safety(cornerOutside, kFALSE), std::sqrt(8.));
1236 double normal[3] = {0., 0., 0.};
1237 const double sidePoint[3] = {radius, 0., 1.};
1238 const double alongX[3] = {1., 0., 0.};
1239 solid.ComputeNormal(sidePoint, alongX, normal);
1240 checkClose(normal[0], 1.);
1241 checkClose(normal[1], 0.);
1242 checkClose(normal[2], 0.);
1243 const double capPoint[3] = {0.5, 0.5, halfHeight};
1244 const double alongZ[3] = {0., 0., 1.};
1245 solid.ComputeNormal(capPoint, alongZ, normal);
1246 checkClose(normal[2], 1.);
1248 checkClose(solid.Capacity(),
reference.Capacity(), 1.e-9);
1250 int meshVertices = 0;
1251 int meshSegments = 0;
1252 int meshPolygons = 0;
1253 solid.GetMeshNumbers(meshVertices, meshSegments, meshPolygons);
1254 BOOST_CHECK_GT(meshVertices, 0);
1255 BOOST_CHECK_GT(meshPolygons, 0);
1260 constexpr double innerRadius = 1.;
1261 constexpr double outerRadius = 2.;
1262 constexpr double halfHeight = 3.;
1264 SurfaceSolid solid(
"hollowCylinder");
1265 BOOST_REQUIRE(solid.AddCylindricalSurface({0., 0., 0.}, {0., 0., 1.}, {1., 0., 0.}, outerRadius, -halfHeight,
1267 BOOST_REQUIRE(solid.AddCylindricalSurface({0., 0., 0.}, {0., 0., 1.}, {1., 0., 0.}, innerRadius, -halfHeight,
1269 BOOST_REQUIRE(addDiskSurface(solid, {0., 0., halfHeight}, {1., 0., 0.}, {0., 1., 0.}, outerRadius,
1271 BOOST_REQUIRE(addDiskSurface(solid, {0., 0., -halfHeight}, {1., 0., 0.}, {0., -1., 0.}, outerRadius,
1278 TGeoTube
reference(
"referenceHollowTube", innerRadius, outerRadius, halfHeight);
1279 compareContainsGrid(solid,
reference, 4., 9);
1282 compareDistance(solid,
reference, {0., 0., 0.}, {1., 0., 0.});
1283 compareDistance(solid,
reference, {0., 0., 2.}, unitDirection(0.4, 0.2, -1.));
1285 compareDistance(solid,
reference, {1.5, 0., 0.}, {1., 0., 0.});
1286 compareDistance(solid,
reference, {1.5, 0., 0.}, {-1., 0., 0.});
1287 compareDistance(solid,
reference, {-1.2, 0.8, 1.}, unitDirection(-0.2, 0.9, 0.4));
1288 compareDistance(solid,
reference, {5., 0., 0.}, {-1., 0., 0.});
1291 const double materialPoint[3] = {1.5, 0., 0.};
1292 checkClose(solid.Safety(materialPoint, kTRUE), 0.5);
1294 checkClose(solid.Capacity(),
reference.Capacity(), 1.e-9);
1299 constexpr double radius = 2.5;
1301 SurfaceSolid solid(
"fullSphere");
1302 BOOST_REQUIRE(solid.AddSphericalSurface({0., 0., 0.}, {0., 0., 1.}, {1., 0., 0.}, radius));
1309 TGeoSphere
reference(
"referenceSphere", 0., radius);
1310 compareContainsGrid(solid,
reference, 3.5, 9);
1312 compareDistance(solid,
reference, {5., 0., 0.}, {-1., 0., 0.});
1313 compareDistance(solid,
reference, {0., 0., 0.}, {1., 0., 0.});
1314 compareDistance(solid,
reference, {0., 0., 0.}, unitDirection(1., 1., 1.));
1315 compareDistance(solid,
reference, {1., 1., 1.}, unitDirection(-0.3, 0.5, 0.8));
1316 compareDistance(solid,
reference, {-4., 0.5, 0.5}, {1., 0., 0.});
1317 compareDistance(solid,
reference, {-4., 2.6, 0.}, {1., 0., 0.});
1320 const double insidePoint[3] = {1., 0., 0.};
1321 const double outsidePoint[3] = {4., 0., 0.};
1322 checkClose(solid.Safety(insidePoint, kTRUE), radius - 1.);
1323 checkClose(solid.Safety(outsidePoint, kFALSE), 4. - radius);
1325 double normal[3] = {0., 0., 0.};
1326 const double surfacePoint[3] = {radius, 0., 0.};
1327 const double alongX[3] = {1., 0., 0.};
1328 solid.ComputeNormal(surfacePoint, alongX, normal);
1329 checkClose(normal[0], 1.);
1331 checkClose(solid.Capacity(),
reference.Capacity(), 1.e-9);
1333 int meshVertices = 0;
1334 int meshSegments = 0;
1335 int meshPolygons = 0;
1336 solid.GetMeshNumbers(meshVertices, meshSegments, meshPolygons);
1337 BOOST_CHECK_GT(meshVertices, 0);
1338 BOOST_CHECK_GT(meshPolygons, 0);
1346 BOOST_REQUIRE(hemisphere.
initialize({0., 0., 0.}, {0., 0., 1.}, {1., 0., 0.}, 2., 0.,
surf::kHalfPi, 0.,
1353 std::vector<surf::RayHit> hits;
1355 BOOST_REQUIRE_EQUAL(hits.size(), 1u);
1356 checkClose(hits.front().distance, 3.);
1357 checkClose(hits.front().normal.zCoord, 1.);
1381 constexpr double halfHeight = 3.;
1382 constexpr double radiusAtBottom = 2.;
1383 constexpr double radiusAtTop = 1.;
1385 SurfaceSolid solid(
"truncatedCone");
1386 BOOST_REQUIRE(solid.AddConicalSurface({0., 0., 0.}, {0., 0., 1.}, {1., 0., 0.}, radiusAtBottom, radiusAtTop,
1387 -halfHeight, halfHeight));
1388 BOOST_REQUIRE(addDiskSurface(solid, {0., 0., halfHeight}, {1., 0., 0.}, {0., 1., 0.}, radiusAtTop));
1389 BOOST_REQUIRE(addDiskSurface(solid, {0., 0., -halfHeight}, {1., 0., 0.}, {0., -1., 0.}, radiusAtBottom));
1395 TGeoCone
reference(
"referenceCone", halfHeight, 0., radiusAtBottom, 0., radiusAtTop);
1396 compareContainsGrid(solid,
reference, 3.5, 9);
1398 compareDistance(solid,
reference, {5., 0., 0.}, {-1., 0., 0.});
1399 compareDistance(solid,
reference, {0., 0., 5.}, {0., 0., -1.});
1400 compareDistance(solid,
reference, {0., 0., 0.}, {1., 0., 0.});
1401 compareDistance(solid,
reference, {0., 0., 0.}, {0., 0., 1.});
1402 compareDistance(solid,
reference, {0., 0., 0.}, {0., 0., -1.});
1403 compareDistance(solid,
reference, {0.5, -0.3, 1.}, unitDirection(0.6, 0.4, 0.2));
1404 compareDistance(solid,
reference, {-4., 0.2, -2.}, unitDirection(1., 0.05, 0.3));
1408 const double center[3] = {0., 0., 0.};
1409 checkClose(solid.Safety(
center, kTRUE), 9. / std::sqrt(37.));
1412 double normal[3] = {0., 0., 0.};
1413 double referenceNormal[3] = {0., 0., 0.};
1414 const double sidePoint[3] = {1.5, 0., 0.};
1415 const double alongX[3] = {1., 0., 0.};
1416 solid.ComputeNormal(sidePoint, alongX, normal);
1417 reference.ComputeNormal(sidePoint, alongX, referenceNormal);
1418 checkClose(normal[0], referenceNormal[0]);
1419 checkClose(normal[1], referenceNormal[1]);
1420 checkClose(normal[2], referenceNormal[2]);
1422 checkClose(solid.Capacity(),
reference.Capacity(), 1.e-9);
1428 constexpr double halfHeight = 1.5;
1429 constexpr double baseRadius = 3.;
1431 SurfaceSolid solid(
"apexCone");
1432 BOOST_REQUIRE(solid.AddConicalSurface({0., 0., 0.}, {0., 0., 1.}, {1., 0., 0.}, baseRadius, 0., -halfHeight,
1434 BOOST_REQUIRE(addDiskSurface(solid, {0., 0., -halfHeight}, {1., 0., 0.}, {0., -1., 0.}, baseRadius));
1442 const auto analyticInside = [&](
double x,
double y,
double z) {
1443 return std::abs(
z) < halfHeight && std::hypot(
x,
y) < halfHeight -
z;
1445 const std::array<std::array<double, 3>, 7> probePoints{{{0., 0., 0.},
1452 for (
const auto& probe : probePoints) {
1453 BOOST_TEST_CONTEXT(
"point = (" << probe[0] <<
", " << probe[1] <<
", " << probe[2] <<
")")
1455 BOOST_CHECK_EQUAL(solid.Contains(probe.data()), analyticInside(probe[0], probe[1], probe[2]));
1460 const double insidePoint[3] = {0., 0., -1.};
1461 const double alongX[3] = {1., 0., 0.};
1462 checkClose(solid.DistFromInside(insidePoint, alongX, 3), 2.5);
1465 const double center[3] = {0., 0., 0.};
1466 checkClose(solid.Safety(
center, kTRUE), halfHeight / std::sqrt(2.), 1.e-9);
1469 checkClose(solid.Capacity(),
surf::kPi * baseRadius * baseRadius * 2. * halfHeight / 3., 1.e-9);
1475 constexpr double majorR = 3.;
1476 constexpr double minorR = 1.;
1479 BOOST_REQUIRE(torus.
initialize({0., 0., 0.}, {0., 0., 1.}, {1., 0., 0.}, majorR, minorR, 0.,
surf::kTwoPi, 0.,
1484 std::vector<surf::RayHit> hits;
1486 BOOST_REQUIRE_EQUAL(hits.size(), 4u);
1487 std::sort(hits.begin(), hits.end(),
1489 checkClose(hits[0].
distance, 6., 1.e-7);
1490 checkClose(hits[1].
distance, 8., 1.e-7);
1491 checkClose(hits[2].
distance, 12., 1.e-7);
1492 checkClose(hits[3].
distance, 14., 1.e-7);
1494 BOOST_CHECK_LT(hits[0].normal.xCoord, 0.);
1495 BOOST_CHECK_GT(hits[1].normal.xCoord, 0.);
1496 BOOST_CHECK_LT(hits[2].normal.xCoord, 0.);
1497 BOOST_CHECK_GT(hits[3].normal.xCoord, 0.);
1511 checkClose(outerNormal.
xCoord, 1.);
1513 checkClose(innerNormal.
xCoord, -1.);
1516 checkClose(topNormal.
zCoord, 1.);
1520 checkClose(torus.
distanceSqToPatch({0., 0., 0.}), (majorR - minorR) * (majorR - minorR));
1535 BOOST_REQUIRE(quarterTube.
initialize({0., 0., 0.}, {0., 0., 1.}, {1., 0., 0.}, majorR, minorR, 0.,
surf::kTwoPi, 0.,
1542 BOOST_REQUIRE_EQUAL(hits.size(), 1u);
1543 checkClose(hits.front().distance, 10. + minorR, 1.e-7);
1548 constexpr double majorR = 3.;
1549 constexpr double minorR = 1.;
1551 SurfaceSolid solid(
"fullTorus");
1552 BOOST_REQUIRE(solid.AddToroidalSurface({0., 0., 0.}, {0., 0., 1.}, {1., 0., 0.}, majorR, minorR));
1560 TGeoTorus
reference(
"referenceTorus", majorR, 0., minorR);
1561 compareContainsGrid(solid,
reference, 4.5, 11);
1564 const double outsidePoint[3] = {-10., 0., 0.};
1565 const double alongX[3] = {1., 0., 0.};
1566 checkClose(solid.DistFromOutside(outsidePoint, alongX, 3), 6., 1.e-7);
1567 const double materialPoint[3] = {majorR + minorR - 0.25, 0., 0.};
1569 checkClose(solid.DistFromInside(materialPoint, alongX, 3), 0.25, 1.e-7);
1572 compareDistance(solid,
reference, {-10., 0.3, 0.2}, {1., 0., 0.}, 1.e-6);
1573 compareDistance(solid,
reference, {0., 0., 5.}, {0., 0., -1.}, 1.e-6);
1576 checkClose(solid.Capacity(),
reference.Capacity(), 1.e-7);
1577 checkClose(solid.Capacity(), 2. *
surf::kPi *
surf::kPi * majorR * minorR * minorR, 1.e-9);
1579 int meshVertices = 0;
1580 int meshSegments = 0;
1581 int meshPolygons = 0;
1582 solid.GetMeshNumbers(meshVertices, meshSegments, meshPolygons);
1583 BOOST_CHECK_GT(meshVertices, 0);
1584 BOOST_CHECK_GT(meshPolygons, 0);
1593 constexpr double majorR = 4.;
1594 constexpr double minorR = 1.5;
1595 constexpr double tubeLow = 0.2;
1596 constexpr double tubeHigh = 2.0;
1600 BOOST_REQUIRE(scalarSection.
initialize({0.2, -0.1, 0.3}, {0., 0., 1.}, {1., 0., 0.}, majorR, minorR, 0.,
1604 const auto wire = paramRectWireCurves(0.,
surf::kHalfPi, tubeLow, tubeHigh);
1605 BOOST_REQUIRE(wireSection.
initialize({0.2, -0.1, 0.3}, {0., 0., 1.}, {1., 0., 0.}, majorR, minorR, 0.,
surf::kHalfPi,
1606 tubeLow, tubeHigh - tubeLow,
false, wire, {}, error));
1610 for (
double ring : {0.1, 0.7, 1.2, 1.7, 2.5}) {
1611 for (
double tube : {0.3, 0.8, 1.5, 1.9, 2.6}) {
1613 BOOST_TEST_CONTEXT(
"ring = " << ring <<
" tube = " << tube)
1629 constexpr double halfX = 1.;
1630 constexpr double halfY = 2.;
1631 constexpr double halfZ = 3.;
1633 SurfaceSolid solid(
"bvhBox");
1634 addBoxSurfaces(solid, halfX, halfY, halfZ);
1637 BOOST_REQUIRE(solid.HasBVH());
1643 BOOST_REQUIRE(solid.GetBVHRootBounds(
lower,
upper));
1644 const Point3D exactLower{-halfX, -halfY, -halfZ};
1645 const Point3D exactUpper{halfX, halfY, halfZ};
1647 for (
int dimension = 0; dimension < 3; ++dimension) {
1648 BOOST_TEST_CONTEXT(
"dimension = " << dimension)
1658 BOOST_CHECK_GE(solid.CountBVHRayCandidates({-2., 0., 0.}, {1., 0., 0.}), 2);
1664 constexpr double half = 1.;
1665 constexpr double centerX = 3.;
1666 SurfaceSolid twoBoxes(
"twoBoxes");
1667 addBoxSurfaces(twoBoxes,
half,
half,
half, {-centerX, 0., 0.});
1668 addBoxSurfaces(twoBoxes,
half,
half,
half, {centerX, 0., 0.});
1669 twoBoxes.CloseShape();
1670 BOOST_REQUIRE(twoBoxes.HasBVH());
1675 const auto analyticInside = [&](
const double* point) {
1676 return (std::abs(std::abs(point[0]) - centerX) <
half) && std::abs(point[1]) <
half && std::abs(point[2]) <
half;
1679 constexpr double extent = 5.;
1680 for (
int stepX = 0; stepX <
samples; ++stepX) {
1681 for (
int stepY = 0; stepY <
samples; ++stepY) {
1682 for (
int stepZ = 0; stepZ <
samples; ++stepZ) {
1683 const double point[3] = {-extent + 2. * extent * (stepX + 0.517) /
samples,
1684 -extent + 2. * extent * (stepY + 0.263) /
samples,
1685 -extent + 2. * extent * (stepZ + 0.741) /
samples};
1686 BOOST_TEST_CONTEXT(
"point = (" << point[0] <<
", " << point[1] <<
", " << point[2] <<
")")
1688 const bool bvhInside = twoBoxes.Contains(point);
1699 constexpr double halfX = 1.;
1700 constexpr double halfY = 2.;
1701 constexpr double halfZ = 3.;
1703 SurfaceSolid
box(
"boundaryBox");
1704 addBoxSurfaces(
box, halfX, halfY, halfZ);
1709 const std::array<std::array<double, 3>, 6> boundaryPoints{{
1713 {-halfX, 0., halfZ},
1714 {halfX, halfY, halfZ},
1715 {-halfX, -halfY, -halfZ}
1717 for (
const auto& point : boundaryPoints) {
1718 BOOST_TEST_CONTEXT(
"point = (" << point[0] <<
", " << point[1] <<
", " << point[2] <<
")")
1727 constexpr double radius = 1.;
1728 constexpr double halfHeight = 1.5;
1729 SurfaceSolid capsule(
"capsule");
1730 BOOST_REQUIRE(capsule.AddCylindricalSurface({0., 0., 0.}, {0., 0., 1.}, {1., 0., 0.}, radius, -halfHeight,
1732 BOOST_REQUIRE(capsule.AddSphericalSurface({0., 0., halfHeight}, {0., 0., 1.}, {1., 0., 0.}, radius, 0.,
1734 BOOST_REQUIRE(capsule.AddSphericalSurface({0., 0., -halfHeight}, {0., 0., -1.}, {1., 0., 0.}, radius, 0.,
1736 capsule.CloseShape();
1737 BOOST_REQUIRE(capsule.HasBVH());
1741 const auto capsuleInside = [&](
const double* point) {
1742 const double axialDistance = std::max(0., std::abs(point[2]) - halfHeight);
1743 return std::hypot(point[0], point[1], axialDistance) < radius;
1746 constexpr double extent = 3.;
1747 for (
int stepX = 0; stepX <
samples; ++stepX) {
1748 for (
int stepY = 0; stepY <
samples; ++stepY) {
1749 for (
int stepZ = 0; stepZ <
samples; ++stepZ) {
1750 const double point[3] = {-extent + 2. * extent * (stepX + 0.517) /
samples,
1751 -extent + 2. * extent * (stepY + 0.263) /
samples,
1752 -extent + 2. * extent * (stepZ + 0.741) /
samples};
1753 BOOST_TEST_CONTEXT(
"point = (" << point[0] <<
", " << point[1] <<
", " << point[2] <<
")")
1755 const bool bvhInside = capsule.Contains(point);
1764 const double onBarrel[3] = {radius, 0., 0.5};
1765 const double onCapApex[3] = {0., 0., halfHeight + radius};
1766 const double aboveApex[3] = {0., 0., halfHeight + radius + 0.01};
1767 const double onRim[3] = {radius, 0., halfHeight};
1774 checkClose(capsule.Capacity(),
1775 surf::kPi * radius * radius * 2. * halfHeight + 4. *
surf::kPi * radius * radius * radius / 3., 1.e-9);
1784 const std::array<std::pair<std::unique_ptr<SurfaceSolid>,
double>, 7> fixtures{{
1785 {makeBoxSolid(
"loopBox", 1., 2., 3.), 4.},
1786 {makeTubeSolid(
"loopTube", 0., 2., 3.), 4.},
1787 {makeTubeSolid(
"loopHollowTube", 1., 2., 3.), 4.},
1788 {makeConeSolid(
"loopCone", 2., 1., 3.), 4.},
1789 {makeSphereSolid(
"loopSphere", 2.5), 3.5},
1790 {makeTorusSolid(
"loopTorus", 3., 1.), 4.5},
1791 {makeCapsuleSolid(
"loopCapsule", 1., 1.5), 3.},
1794 for (
const auto& [solid, extent] : fixtures) {
1795 BOOST_TEST_CONTEXT(
"fixture = " << solid->GetName())
1797 BOOST_REQUIRE(solid->HasBVH());
1798 sweepDistanceAgainstLoop(*solid, extent, 5);
1810 const auto box = makeBoxSolid(
"sweepBox", 1., 2., 3.);
1811 TGeoBBox boxReference(
"sweepBoxReference", 1., 2., 3.);
1812 sweepDistanceAgainstReference(*
box, boxReference, 4.,
samples);
1814 const auto tube = makeTubeSolid(
"sweepTube", 0., 2., 3.);
1815 TGeoTube tubeReference(
"sweepTubeReference", 0., 2., 3.);
1816 sweepDistanceAgainstReference(*tube, tubeReference, 4.,
samples);
1818 const auto hollowTube = makeTubeSolid(
"sweepHollowTube", 1., 2., 3.);
1819 TGeoTube hollowTubeReference(
"sweepHollowTubeReference", 1., 2., 3.);
1820 sweepDistanceAgainstReference(*hollowTube, hollowTubeReference, 4.,
samples);
1822 const auto cone = makeConeSolid(
"sweepCone", 2., 1., 3.);
1823 TGeoCone coneReference(
"sweepConeReference", 3., 0., 2., 0., 1.);
1824 sweepDistanceAgainstReference(*cone, coneReference, 4.,
samples);
1826 const auto sphere = makeSphereSolid(
"sweepSphere", 2.5);
1827 TGeoSphere sphereReference(
"sweepSphereReference", 0., 2.5);
1828 sweepDistanceAgainstReference(*sphere, sphereReference, 3.5,
samples);
1831 const auto torus = makeTorusSolid(
"sweepTorus", 3., 1.);
1832 TGeoTorus torusReference(
"sweepTorusReference", 3., 0., 1.);
1833 sweepDistanceAgainstReference(*torus, torusReference, 4.5,
samples, 1.e-6);
1838 constexpr double halfX = 1.;
1839 constexpr double halfY = 2.;
1840 constexpr double halfZ = 3.;
1841 const auto box = makeBoxSolid(
"hardCaseBox", halfX, halfY, halfZ);
1848 const std::array<std::array<double, 3>, 4> throughFeature{{
1850 {-5., halfY, halfZ},
1854 const std::array<std::array<double, 3>, 4> throughFeatureDirection{{
1857 unitDirection(0., 1., 1.5),
1858 unitDirection(1., 2., 3.),
1861 checkDistanceAgainstLoop(*
box, throughFeature[
index], throughFeatureDirection[
index]);
1867 const std::array<std::array<double, 3>, 3> grazing{{
1872 const std::array<std::array<double, 3>, 3> grazingDirection{{
1875 unitDirection(1., 1., 0.),
1878 checkDistanceAgainstLoop(*
box, grazing[
index], grazingDirection[
index]);
1881 const auto tube = makeTubeSolid(
"hardCaseTube", 0., 2., 3.);
1882 checkDistanceAgainstLoop(*tube, {-5., 2., 0.}, {1., 0., 0.});
1883 checkDistanceAgainstLoop(*tube, {-5., 2. - 1.e-7, 0.}, {1., 0., 0.});
1884 checkDistanceAgainstLoop(*tube, {-5., 2. + 1.e-7, 0.}, {1., 0., 0.});
1889 const std::array<std::array<double, 3>, 4> onSurface{{
1893 {halfX, halfY, halfZ},
1895 for (
const auto& point : onSurface) {
1896 for (
const auto& direction : probeDirections()) {
1897 checkDistanceAgainstLoop(*
box, point, direction);
1902 const std::array<double, 3> justOutside{halfX + 1.e-6, 0., 0.};
1903 const std::array<double, 3> justInside{halfX - 1.e-6, 0., 0.};
1904 const std::array<double, 3> inward{-1., 0., 0.};
1905 checkClose(
box->DistFromOutside(justOutside.data(), inward.data(), 3), 1.e-6, 1.e-12);
1906 checkClose(
box->DistFromInside(justInside.data(), inward.data(), 3), 2. * halfX - 1.e-6, 1.e-12);
1907 checkClose(
box->DistFromOutside(justOutside.data(), inward.data(), 3),
1908 reference.DistFromOutside(justOutside.data(), inward.data(), 3), 1.e-12);
1911 const std::array<double, 3> farOutside{-5., 0., 0.};
1912 const std::array<double, 3> alongX{1., 0., 0.};
1913 const double entryDistance =
box->DistFromOutside(farOutside.data(), alongX.data(), 3);
1914 checkClose(entryDistance, 4.);
1916 BOOST_CHECK_EQUAL(
box->DistFromOutside(farOutside.data(), alongX.data(), 3, entryDistance * 0.5),
1919 BOOST_CHECK_EQUAL(
box->DistFromOutside(farOutside.data(), alongX.data(), 3, entryDistance - 1.e-3),
1922 checkClose(
box->DistFromOutside(farOutside.data(), alongX.data(), 3, entryDistance + 1.e-3), entryDistance);
1923 checkClose(
box->DistFromOutside(farOutside.data(), alongX.data(), 3, 100.), entryDistance);
1925 const std::array<double, 3>
center{0., 0., 0.};
1926 const double exitDistance =
box->DistFromInside(
center.data(), alongX.data(), 3);
1927 checkClose(exitDistance, halfX);
1929 checkClose(
box->DistFromInside(
center.data(), alongX.data(), 3, exitDistance * 2.), exitDistance);
1931 for (
const double stepmax : {entryDistance * 0.5, entryDistance - 1.e-9, entryDistance, entryDistance + 1.e-9,
1932 entryDistance * 2.}) {
1933 checkDistanceAgainstLoop(*
box, farOutside, alongX, stepmax);
1934 checkDistanceAgainstLoop(*
box,
center, alongX, stepmax);
1938 const std::array<double, 3> wayOff{-1000., 0., 0.};
1939 BOOST_CHECK_EQUAL(
box->DistFromOutside(wayOff.data(), alongX.data(), 3, 10.), TGeoShape::Big());
1940 checkClose(
box->DistFromOutside(wayOff.data(), alongX.data(), 3), 999.);
1948 constexpr int boxCount = 8;
1949 constexpr double half = 0.5;
1950 constexpr double spacing = 4.;
1952 SurfaceSolid
row(
"prunedRow");
1953 for (
int boxIndex = 0; boxIndex < boxCount; ++boxIndex) {
1957 BOOST_REQUIRE(
row.HasBVH());
1961 const std::array<double, 3> beforeRow{-5., 0., 0.};
1962 const std::array<double, 3> alongRow{1., 0., 0.};
1966 SurfaceSolid::ResetRayCandidateCounter();
1967 const double prunedDistance =
row.DistFromOutside(beforeRow.data(), alongRow.data(), 3);
1968 const long long prunedCandidates = SurfaceSolid::GetRayCandidateCount();
1970 SurfaceSolid::SetRayTMaxPruning(
false);
1971 SurfaceSolid::ResetRayCandidateCounter();
1972 const double unprunedDistance =
row.DistFromOutside(beforeRow.data(), alongRow.data(), 3);
1973 const long long unprunedCandidates = SurfaceSolid::GetRayCandidateCount();
1974 SurfaceSolid::SetRayTMaxPruning(
true);
1978 checkClose(prunedDistance, 5. -
half);
1980 BOOST_CHECK_GT(prunedCandidates, 0);
1981 BOOST_CHECK_LT(prunedCandidates, unprunedCandidates);
1985 sweepDistanceAgainstLoop(
row, 1.2 * boxCount * spacing / 2., 4);
1988 SurfaceSolid::ResetRayCandidateCounter();
1989 row.DistFromOutside_Loop(beforeRow.data(), alongRow.data());
2018 constexpr double radius = 2.;
2019 constexpr double rayOffsetY = 1.9;
2020 const double sphereHitX = -std::sqrt(radius * radius - rayOffsetY * rayOffsetY);
2021 const std::array<double, 3> rayOrigin{-10., rayOffsetY, 0.};
2022 const std::array<double, 3> alongX{1., 0., 0.};
2023 const double sphereDistance = sphereHitX - rayOrigin[0];
2028 for (
const double relativeOffset : {1.e-5, 3.e-5, 1.e-4, 3.e-4, 1.e-3, 3.e-3, 1.e-2}) {
2029 const double patchX = sphereHitX - relativeOffset * sphereDistance;
2030 BOOST_TEST_CONTEXT(
"relativeOffset = " << relativeOffset <<
" patchX = " << patchX)
2032 SurfaceSolid solid(
"nearTie");
2033 BOOST_REQUIRE(solid.AddSphericalSurface({0., 0., 0.}, {0., 0., 1.}, {1., 0., 0.}, radius));
2035 BOOST_REQUIRE(solid.AddPlanarSurface({patchX, rayOffsetY - 0.1, -0.1}, {0., 0., 1.}, {0., 1., 0.},
2036 rectangleWire(0.2, 0.2)));
2037 solid.CloseShape(
false);
2038 BOOST_REQUIRE(solid.HasBVH());
2041 const double expected = patchX - rayOrigin[0];
2042 for (
const bool pruning : {
true,
false}) {
2043 SurfaceSolid::SetRayTMaxPruning(pruning);
2044 BOOST_TEST_CONTEXT(
"pruning = " << pruning)
2046 const double distance = solid.DistFromOutside(rayOrigin.data(), alongX.data(), 3);
2051 SurfaceSolid::SetRayTMaxPruning(
true);
2061 constexpr double halfLen = 3.;
2062 constexpr double radius = 2.;
2063 constexpr double halfHeight = 4.;
2067 const std::vector<BoundaryCurve> stadiumWire{
2068 BoundaryCurve::makeLine({-halfLen, -radius}, {halfLen, -radius}),
2070 BoundaryCurve::makeLine({halfLen, radius}, {-halfLen, radius}),
2073 SurfaceSolid solid(
"stadiumPrism");
2075 BOOST_REQUIRE(solid.AddCurvedPlanarSurface({0., 0., halfHeight}, {1., 0., 0.}, {0., 1., 0.}, stadiumWire));
2076 BOOST_REQUIRE(solid.AddCurvedPlanarSurface({0., 0., -halfHeight}, {1., 0., 0.}, {0., -1., 0.}, stadiumWire));
2078 BOOST_REQUIRE(solid.AddPlanarSurface({-halfLen, radius, -halfHeight}, {0., 0., 1.}, {1., 0., 0.},
2079 rectangleWire(2. * halfHeight, 2. * halfLen)));
2080 BOOST_REQUIRE(solid.AddPlanarSurface({-halfLen, -radius, -halfHeight}, {1., 0., 0.}, {0., 0., 1.},
2081 rectangleWire(2. * halfLen, 2. * halfHeight)));
2083 BOOST_REQUIRE(solid.AddCylindricalSurface({halfLen, 0., 0.}, {0., 0., 1.}, {1., 0., 0.}, radius, -halfHeight,
2085 BOOST_REQUIRE(solid.AddCylindricalSurface({-halfLen, 0., 0.}, {0., 0., 1.}, {1., 0., 0.}, radius, -halfHeight,
2093 checkClose(solid.Capacity(), (4. * halfLen * radius +
surf::kPi * radius * radius) * 2. * halfHeight, 1.e-6);
2095 const auto stadiumInside = [&](
const double* point) {
2096 if (std::abs(point[2]) > halfHeight) {
2099 const double ax = std::abs(point[0]);
2100 const double dx = ax > halfLen ? ax - halfLen : 0.;
2101 return dx * dx + point[1] * point[1] <= radius * radius;
2105 const double extentX = 6., extentY = 3.5, extentZ = 5.;
2106 for (
int stepX = 0; stepX <
samples; ++stepX) {
2107 for (
int stepY = 0; stepY <
samples; ++stepY) {
2108 for (
int stepZ = 0; stepZ <
samples; ++stepZ) {
2109 const double point[3] = {-extentX + 2. * extentX * (stepX + 0.517) /
samples,
2110 -extentY + 2. * extentY * (stepY + 0.263) /
samples,
2111 -extentZ + 2. * extentZ * (stepZ + 0.741) /
samples};
2112 BOOST_TEST_CONTEXT(
"point = (" << point[0] <<
", " << point[1] <<
", " << point[2] <<
")")
2123 constexpr double radius = 2.;
2124 constexpr double halfHeight = 3.;
2129 SurfaceSolid solid(
"wireTrimmedCylinder");
2130 BOOST_REQUIRE(solid.AddCylindricalSurface({0., 0., 0.}, {0., 0., 1.}, {1., 0., 0.}, radius, -halfHeight,
2132 paramRectWire(0.,
surf::kTwoPi, -halfHeight, halfHeight)));
2133 BOOST_REQUIRE(addDiskSurface(solid, {0., 0., halfHeight}, {1., 0., 0.}, {0., 1., 0.}, radius));
2134 BOOST_REQUIRE(addDiskSurface(solid, {0., 0., -halfHeight}, {1., 0., 0.}, {0., -1., 0.}, radius));
2140 TGeoTube
reference(
"wireTrimTube", 0., radius, halfHeight);
2141 compareContainsGrid(solid,
reference, 4., 9);
2142 compareDistance(solid,
reference, {5., 0., 0.}, {-1., 0., 0.});
2143 compareDistance(solid,
reference, {0., 0., 5.}, {0., 0., -1.});
2144 compareDistance(solid,
reference, {-4., -1., -2.}, unitDirection(1., 0.3, 0.5));
2145 compareDistance(solid,
reference, {0., 0., 0.}, unitDirection(1., 1., 1.));
2146 compareDistance(solid,
reference, {5., 2.5, 0.}, {-1., 0., 0.});
2150 checkClose(solid.Capacity(),
reference.Capacity(), 1.e-6);
2152 double normal[3] = {0., 0., 0.};
2153 const double sidePoint[3] = {radius, 0., 1.};
2154 const double alongX[3] = {1., 0., 0.};
2155 solid.ComputeNormal(sidePoint, alongX, normal);
2156 checkClose(normal[0], 1.);
2157 checkClose(normal[1], 0.);
2158 checkClose(normal[2], 0.);
2163 constexpr double halfHeight = 3.;
2164 constexpr double radiusAtBottom = 2.;
2165 constexpr double radiusAtTop = 1.;
2167 SurfaceSolid solid(
"wireTrimmedCone");
2168 BOOST_REQUIRE(solid.AddConicalSurface({0., 0., 0.}, {0., 0., 1.}, {1., 0., 0.}, radiusAtBottom, radiusAtTop,
2170 paramRectWire(0.,
surf::kTwoPi, -halfHeight, halfHeight)));
2171 BOOST_REQUIRE(addDiskSurface(solid, {0., 0., halfHeight}, {1., 0., 0.}, {0., 1., 0.}, radiusAtTop));
2172 BOOST_REQUIRE(addDiskSurface(solid, {0., 0., -halfHeight}, {1., 0., 0.}, {0., -1., 0.}, radiusAtBottom));
2178 TGeoCone
reference(
"wireTrimCone", halfHeight, 0., radiusAtBottom, 0., radiusAtTop);
2179 compareContainsGrid(solid,
reference, 3.5, 9);
2180 compareDistance(solid,
reference, {5., 0., 0.}, {-1., 0., 0.});
2181 compareDistance(solid,
reference, {0., 0., 0.}, {1., 0., 0.});
2182 compareDistance(solid,
reference, {-4., 0.2, -2.}, unitDirection(1., 0.05, 0.3));
2186 checkClose(solid.Capacity(),
reference.Capacity(), 1.e-9);
2207 constexpr double kRelative = 1.e-12;
2211 BOOST_TEST_CONTEXT(what)
2214 const double exact = rectangle.capacityContribution();
2215 const double contour = wired.capacityContribution();
2216 BOOST_CHECK_GT(std::abs(exact), 1.e-6);
2217 checkClose(contour, exact, kRelative * std::abs(exact));
2224 const double phiLow = 0.3, phiHigh = 2.4, hLow = -0.8, hHigh = 1.9;
2225 BOOST_REQUIRE(rectangle.
initialize(centre, axis,
reference, 1.7, hLow, hHigh, phiLow, phiHigh - phiLow,
false,
2227 BOOST_REQUIRE(wired.
initialize(centre, axis,
reference, 1.7, hLow, hHigh, phiLow, phiHigh - phiLow,
false,
2228 paramRectWireCurves(phiLow, phiHigh, hLow, hHigh), {}, error));
2229 compare(
"cylinder", rectangle, wired);
2234 const double phiLow = -0.4, phiHigh = 1.9, hLow = 0.2, hHigh = 2.1;
2235 BOOST_REQUIRE(rectangle.
initialize(centre, axis,
reference, 1.1, 2.3, hLow, hHigh, phiLow, phiHigh - phiLow,
2237 BOOST_REQUIRE(wired.
initialize(centre, axis,
reference, 1.1, 2.3, hLow, hHigh, phiLow, phiHigh - phiLow,
false,
2238 paramRectWireCurves(phiLow, phiHigh, hLow, hHigh), {}, error));
2239 compare(
"cone", rectangle, wired);
2244 const double phiLow = 0.2, phiHigh = 2.7, thetaLow = 0.4, thetaHigh = 2.3;
2245 BOOST_REQUIRE(rectangle.
initialize(centre, axis,
reference, 2.2, thetaLow, thetaHigh, phiLow, phiHigh - phiLow,
2247 BOOST_REQUIRE(wired.
initialize(centre, axis,
reference, 2.2, thetaLow, thetaHigh, phiLow, phiHigh - phiLow,
false,
2248 paramRectWireCurves(phiLow, phiHigh, thetaLow, thetaHigh), {}, error));
2249 compare(
"sphere", rectangle, wired);
2254 const double ringLow = 0.1, ringHigh = 2.2, tubeLow = -0.3, tubeHigh = 1.8;
2255 BOOST_REQUIRE(rectangle.
initialize(centre, axis,
reference, 4., 1.4, ringLow, ringHigh - ringLow, tubeLow,
2256 tubeHigh - tubeLow,
false, error));
2257 BOOST_REQUIRE(wired.
initialize(centre, axis,
reference, 4., 1.4, ringLow, ringHigh - ringLow, tubeLow,
2258 tubeHigh - tubeLow,
false,
2259 paramRectWireCurves(ringLow, ringHigh, tubeLow, tubeHigh), {}, error));
2260 compare(
"torus", rectangle, wired);
2270 const double radius = 1.7;
2271 const std::vector<surf::Curve2D> trim{surf::Curve2D::makeCircle({1.0, 0.2}, 0.6)};
2272 BOOST_REQUIRE(disk.
initialize(centre, axis,
reference, radius, -2., 2., 0.,
surf::kTwoPi,
false, trim, {},
2278 std::vector<surf::CurveWire> innerWires;
2283 const auto gridRelativeError = [&](
int samples) {
2285 outerWire, innerWires,
2286 [&disk, radius](
double phi,
double height) {
2291 return std::abs(grid - contour) / std::abs(contour);
2293 const double at128 = gridRelativeError(128);
2294 const double at512 = gridRelativeError(512);
2295 const double at2048 = gridRelativeError(2048);
2299 BOOST_CHECK_LT(at512, 1.e-4);
2300 BOOST_CHECK_LT(at2048, 1.e-4);
2306 BOOST_CHECK_GT(at128, 1.e-3);
2307 BOOST_CHECK_GT(at2048, at512);
2317 const auto onCylinder = [](
double phi,
double height) {
2318 return Vec3{2. * std::cos(phi), 2. * std::sin(phi),
height};
2323 const std::vector<Curve2D> outer{Curve2D::makeLine({0., -3.}, {
surf::kTwoPi, -3.}),
2326 Curve2D::makeLine({0., 3.}, {0., -3.})};
2327 const std::vector<Curve2D> hole{Curve2D::makeLine({2.0, -1.}, {2.5, -1.}), Curve2D::makeLine({2.5, -1.}, {2.5, 1.}),
2328 Curve2D::makeLine({2.5, 1.}, {2.0, 1.}), Curve2D::makeLine({2.0, 1.}, {2.0, -1.})};
2329 BOOST_REQUIRE(windowed.
initialize({0., 0., 0.}, {0., 0., 1.}, {1., 0., 0.}, 2., -3., 3., 0.,
surf::kTwoPi,
false,
2330 outer, {hole}, error));
2337 std::vector<surf::RayHit> hits;
2338 windowed.
appendIntersections({0., 0., 0.}, {std::cos(2.25), std::sin(2.25), 0.}, 0., 1.e30, hits);
2341 windowed.
appendIntersections({0., 0., 0.}, {std::cos(0.5), std::sin(0.5), 0.}, 0., 1.e30, hits);
2342 BOOST_REQUIRE_EQUAL(hits.size(), 1u);
2343 checkClose(hits.front().distance, 2.);
2347 const std::vector<Curve2D> arcOuter{Curve2D::makeCircle({
surf::kPi, 0.}, 0.5)};
2348 BOOST_REQUIRE(arcTrim.
initialize({0., 0., 0.}, {0., 0., 1.}, {1., 0., 0.}, 2., -1., 1., 0.,
surf::kTwoPi,
false,
2349 arcOuter, {}, error));
2356 const auto onSphere = [](
double theta,
double phi) {
2357 return Vec3{2. * std::sin(theta) * std::cos(phi), 2. * std::sin(theta) * std::sin(phi), 2. * std::cos(theta)};
2360 BOOST_REQUIRE(sphereWire.
initialize({0., 0., 0.}, {0., 0., 1.}, {1., 0., 0.}, 2.,
surf::kHalfPi / 2.,
surf::kHalfPi,
2365 BOOST_REQUIRE(sphereScalar.
initialize({0., 0., 0.}, {0., 0., 1.}, {1., 0., 0.}, 2.,
surf::kHalfPi / 2.,
2375 const std::vector<Curve2D> wideOuter{Curve2D::makeLine({0., -1.}, {7., -1.}), Curve2D::makeLine({7., -1.}, {7., 1.}),
2376 Curve2D::makeLine({7., 1.}, {0., 1.}), Curve2D::makeLine({0., 1.}, {0., -1.})};
2377 BOOST_CHECK(!tooWide.
initialize({0., 0., 0.}, {0., 0., 1.}, {1., 0., 0.}, 2., -1., 1., 0.,
surf::kTwoPi,
false,
2378 wideOuter, {}, error));
2387template <
typename Po
intAt>
2388void checkMetricAgainstFiniteDifference(
const surf::BoundedSurface& surface,
const PointAt& pointAt,
2389 double uCoord,
double vCoord,
double tolerance = 1.e-6)
2391 const double step = 1.e-5;
2392 const surf::Vec3 dU = (pointAt(uCoord + step, vCoord) - pointAt(uCoord - step, vCoord)) * (0.5 / step);
2393 const surf::Vec3 dV = (pointAt(uCoord, vCoord + step) - pointAt(uCoord, vCoord - step)) * (0.5 / step);
2399 checkClose(gUU,
dot(dU, dU), tolerance);
2400 checkClose(gUV,
dot(dU, dV), tolerance);
2401 checkClose(gVV,
dot(dV, dV), tolerance);
2416 const Vec3 axisU{2., 0., 0.};
2417 const Vec3 axisV{1., 3., 0.};
2418 const std::vector<Vec2> unitSquare{{0., 0.}, {1., 0.}, {1., 1.}, {0., 1.}};
2420 BOOST_REQUIRE(plane.
initialize({0.5, -1., 2.}, axisU, axisV, unitSquare, {}, error));
2421 checkMetricAgainstFiniteDifference(plane, [&](
double u,
double v) {
return plane.
toGlobal({u,
v}); }, 0.3, 0.7);
2427 checkClose(gUU, 4.);
2428 checkClose(gUV, 2.);
2429 checkClose(gVV, 10.);
2437 BOOST_REQUIRE(curvedPlane.
initialize({0., 0., 0.}, {1., 0., 0.}, {0., 1., 0.},
2438 {surf::Curve2D::makeCircle({0., 0.}, 1.)}, {},
error));
2439 checkMetricAgainstFiniteDifference(
2440 curvedPlane, [&](
double u,
double v) {
return curvedPlane.toGlobal({u,
v}); }, 0.2, -0.4);
2444 for (
const double radius : {0.01, 100.}) {
2446 BOOST_REQUIRE(cylinder.
initialize({0., 0., 0.}, {0., 0., 1.}, {1., 0., 0.}, radius, -1., 1., 0.,
surf::kTwoPi,
2448 checkMetricAgainstFiniteDifference(cylinder, [&](
double u,
double v) {
return cylinder.
pointAt(u,
v); }, 1.1, 0.3, 1.e-4 * radius * radius);
2450 checkClose(std::sqrt(cylinder.
parametricLengthSqAt({1.1, 0.3}, {2.e-5, 0.})), 2.e-5 * radius, 1.e-12);
2455 BOOST_REQUIRE(sphere.
initialize({1., 2., 3.}, {0., 0., 1.}, {1., 0., 0.}, 2.5, 0.,
surf::kPi, 0.,
surf::kTwoPi,
2457 checkMetricAgainstFiniteDifference(sphere, [&](
double u,
double v) {
return sphere.
pointAt(
v, u); }, 0.9, 1.2);
2464 checkClose(gUU, 0.);
2465 checkClose(gVV, 2.5 * 2.5);
2468 checkClose(gUU, 0.);
2471 checkClose(gUU, 2.5 * 2.5);
2477 BOOST_REQUIRE(cone.
initialize({0., 0., 0.}, {0., 0., 1.}, {1., 0., 0.}, 0., 4., 0., 2., 0.,
surf::kTwoPi,
false,
2479 checkMetricAgainstFiniteDifference(cone, [&](
double u,
double v) {
return cone.
pointAt(u,
v); }, 2.0, 1.3);
2485 checkClose(gUU, 0.);
2486 checkClose(gVV, 1. + 2. * 2.);
2489 checkClose(gUU, 16.);
2495 BOOST_REQUIRE(torus.
initialize({0., 0., 0.}, {0., 0., 1.}, {1., 0., 0.}, 5., 1.5, 0.,
surf::kTwoPi, 0.,
2497 checkMetricAgainstFiniteDifference(torus, [&](
double u,
double v) {
return torus.
pointAt(u,
v); }, 0.7, 2.1);
2503 checkClose(gUU, 6.5 * 6.5);
2504 checkClose(gVV, 1.5 * 1.5);
2506 checkClose(gUU, 3.5 * 3.5);
2520 const auto trimWithPhiDrift = [](
double drift) {
2521 return std::vector<Curve2D>{Curve2D::makeLine({0.2, -1.}, {1.2, -1.}), Curve2D::makeLine({1.2, -1.}, {1.2, 1.}),
2522 Curve2D::makeLine({1.2, 1.}, {0.2 + drift, 1.}),
2523 Curve2D::makeLine({0.2 + drift, 1.}, {0.2 + drift, -1.})};
2525 const auto acceptsDrift = [&](
double radius,
double drift) {
2527 return cylinder.
initialize({0., 0., 0.}, {0., 0., 1.}, {1., 0., 0.}, radius, -1., 1., 0.,
surf::kTwoPi,
false,
2528 trimWithPhiDrift(drift), {}, error);
2557 const auto polygonAcceptsGap = [](
double gap) {
2558 const std::vector<SurfaceEdge>
edges{{{0., 0.}, {1., 0.}},
2559 {{1., 0.}, {1., 1.}},
2560 {{1., 1.}, {0., 1.}},
2561 {{0., 1.}, {gap, 0.}}};
2563 WireStatus status = WireStatus::Valid;
2566 const auto curveAcceptsGap = [](
double gap) {
2567 const std::vector<surf::Curve2D> curves{surf::Curve2D::makeLine({0., 0.}, {1., 0.}),
2568 surf::Curve2D::makeLine({1., 0.}, {1., 1.}),
2569 surf::Curve2D::makeLine({1., 1.}, {0., 1.}),
2570 surf::Curve2D::makeLine({0., 1.}, {gap, 0.})};
2572 WireStatus status = WireStatus::Valid;
2573 return wire.
initialize(curves, WireRole::Outer, status, {});
2578 for (
const double gap : {0., 1.e-8, 1.e-7}) {
2583 for (
const double gap : {1.e-5, 1.e-4}) {
2594void appendU32(std::vector<char>& bytes, uint32_t
value)
2596 const char*
raw =
reinterpret_cast<const char*
>(&
value);
2597 bytes.insert(bytes.end(),
raw,
raw +
sizeof(
value));
2600void appendDoubles(std::vector<char>& bytes, std::initializer_list<double>
values)
2603 const char*
raw =
reinterpret_cast<const char*
>(&
value);
2604 bytes.insert(bytes.end(),
raw,
raw +
sizeof(
value));
2611void appendSidecarHeader(std::vector<char>& bytes, uint32_t nSurfaces, uint32_t
version = 1,
2612 double modelTolerance = 0., uint32_t nModelEdges = 0)
2614 bytes.insert(bytes.end(), {
'O',
'2',
'S',
'S'});
2616 appendU32(bytes, nSurfaces);
2617 appendU32(bytes, 0);
2619 appendDoubles(bytes, {modelTolerance});
2622 appendU32(bytes, nModelEdges);
2627void appendPlaneRecord(std::vector<char>& bytes,
const FaceFrame& frame)
2629 appendU32(bytes, 1);
2630 appendU32(bytes, 0);
2631 appendU32(bytes, 9);
2632 appendDoubles(bytes, {frame.origin[0], frame.origin[1], frame.origin[2], frame.axisU[0], frame.axisU[1],
2633 frame.axisU[2], frame.axisV[0], frame.axisV[1], frame.axisV[2]});
2634 appendU32(bytes, 1);
2635 appendU32(bytes, 0);
2636 appendU32(bytes, 4);
2637 const double extentU = frame.extentU;
2638 const double extentV = frame.extentV;
2639 const std::array<std::array<double, 4>, 4>
edges{{{0., 0., extentU, 0.},
2640 {extentU, 0., extentU, extentV},
2641 {extentU, extentV, 0., extentV},
2642 {0., extentV, 0., 0.}}};
2643 for (
const auto& edge :
edges) {
2644 appendU32(bytes, 0);
2645 appendU32(bytes, 4);
2646 appendDoubles(bytes, {edge[0], edge[1], edge[2], edge[3]});
2652void appendDiskPlaneRecord(std::vector<char>& bytes,
const Point3D&
center,
const Point3D& axisU,
2653 const Point3D& axisV,
double radius,
double holeRadius = 0.)
2655 appendU32(bytes, 1);
2656 appendU32(bytes, 0);
2657 appendU32(bytes, 9);
2658 appendDoubles(bytes, {
center[0],
center[1],
center[2], axisU[0], axisU[1], axisU[2], axisV[0], axisV[1], axisV[2]});
2659 const uint32_t nWires = holeRadius > 0. ? 2u : 1u;
2660 appendU32(bytes, nWires);
2661 appendU32(bytes, 0);
2662 appendU32(bytes, 1);
2663 appendU32(bytes, 1);
2664 appendU32(bytes, 5);
2665 appendDoubles(bytes, {0., 0., radius, 0., 2. *
surf::kPi});
2666 if (holeRadius > 0.) {
2667 appendU32(bytes, 1);
2668 appendU32(bytes, 1);
2669 appendU32(bytes, 1);
2670 appendU32(bytes, 5);
2671 appendDoubles(bytes, {0., 0., holeRadius, 0., -2. *
surf::kPi});
2675std::filesystem::path writeSidecarFile(
const std::string&
name,
const std::vector<char>& bytes)
2677 const auto path = std::filesystem::temp_directory_path() /
name;
2678 std::ofstream out(
path, std::ios::binary);
2679 out.write(bytes.data(),
static_cast<std::streamsize
>(bytes.size()));
2680 BOOST_REQUIRE(out.good());
2688 constexpr double halfX = 1.;
2689 constexpr double halfY = 2.;
2690 constexpr double halfZ = 3.;
2692 std::vector<char> boxBytes;
2693 appendSidecarHeader(boxBytes, 6);
2695 appendPlaneRecord(boxBytes, boxFaceFrame(
faceIndex, halfX, halfY, halfZ));
2697 const auto boxPath = writeSidecarFile(
"o2_sidecar_roundtrip_box.bin", boxBytes);
2699 SurfaceSolid
box(
"sidecarBox");
2701 std::filesystem::remove(boxPath);
2707 TGeoBBox referenceBox(
"referenceBox", halfX, halfY, halfZ);
2708 compareContainsGrid(
box, referenceBox, 4., 7);
2709 compareDistance(
box, referenceBox, {5., 0.5, 0.5}, {-1., 0., 0.});
2710 compareDistance(
box, referenceBox, {0., 0., 0.}, unitDirection(1., 1., 1.));
2711 checkClose(
box.Capacity(), referenceBox.Capacity(), 1.e-9);
2714 constexpr double radius = 2.;
2715 constexpr double halfHeight = 3.;
2717 std::vector<char> tubeBytes;
2718 appendSidecarHeader(tubeBytes, 3);
2719 appendU32(tubeBytes, 2);
2720 appendU32(tubeBytes, 0);
2721 appendU32(tubeBytes, 14);
2722 appendDoubles(tubeBytes, {0., 0., 0., 0., 0., 1., 1., 0., 0., radius, -halfHeight, halfHeight, 0., 2. *
surf::kPi});
2723 appendU32(tubeBytes, 0);
2725 appendDiskPlaneRecord(tubeBytes, {0., 0., halfHeight}, {1., 0., 0.}, {0., 1., 0.}, radius);
2726 appendDiskPlaneRecord(tubeBytes, {0., 0., -halfHeight}, {1., 0., 0.}, {0., -1., 0.}, radius);
2727 const auto tubePath = writeSidecarFile(
"o2_sidecar_roundtrip_tube.bin", tubeBytes);
2729 SurfaceSolid tube(
"sidecarTube");
2731 std::filesystem::remove(tubePath);
2737 TGeoTube referenceTube(
"referenceTube", 0., radius, halfHeight);
2738 compareContainsGrid(tube, referenceTube, 4., 7);
2739 compareDistance(tube, referenceTube, {5., 0.5, 1.}, {-1., 0., 0.});
2740 compareDistance(tube, referenceTube, {0., 0., 0.}, unitDirection(1., 1., 1.));
2741 checkClose(tube.Capacity(), referenceTube.Capacity(), 1.e-9);
2744 const auto badPath = writeSidecarFile(
"o2_sidecar_bad_magic.bin", {
'X',
'X',
'X',
'X', 0, 0, 0, 0});
2745 SurfaceSolid bad(
"sidecarBad");
2747 std::filesystem::remove(badPath);
2752 std::vector<char> truncatedBytes;
2753 appendSidecarHeader(truncatedBytes, 1);
2754 const auto truncatedPath = writeSidecarFile(
"o2_sidecar_truncated.bin", truncatedBytes);
2755 SurfaceSolid truncated(
"sidecarTruncated");
2757 std::filesystem::remove(truncatedPath);
2765 constexpr double halfX = 1.;
2766 constexpr double halfY = 2.;
2767 constexpr double halfZ = 3.;
2768 const auto boxBytesWithHeader = [&](uint32_t
version,
double modelTolerance) {
2769 std::vector<char> bytes;
2770 appendSidecarHeader(bytes, 6,
version, modelTolerance);
2772 appendPlaneRecord(bytes, boxFaceFrame(
faceIndex, halfX, halfY, halfZ));
2776 const auto loadFrom = [](
const char*
name,
const std::vector<char>& bytes, SurfaceSolid& solid) {
2777 const auto path = writeSidecarFile(
name, bytes);
2779 std::filesystem::remove(
path);
2784 SurfaceSolid
v2(
"sidecarV2");
2785 BOOST_REQUIRE(loadFrom(
"o2_sidecar_v2.bin", boxBytesWithHeader(2, 3.5e-5),
v2));
2787 checkClose(
v2.GetModelTolerance(), 3.5e-5, 1.e-18);
2789 checkClose(
v2.GetModelTolerance(), 3.5e-5, 1.e-18);
2792 SurfaceSolid v2Silent(
"sidecarV2Silent");
2793 BOOST_REQUIRE(loadFrom(
"o2_sidecar_v2_silent.bin", boxBytesWithHeader(2, 0.), v2Silent));
2797 SurfaceSolid
v1(
"sidecarV1");
2798 BOOST_REQUIRE(loadFrom(
"o2_sidecar_v1.bin", boxBytesWithHeader(1, 0.),
v1));
2800 checkClose(
v1.GetModelTolerance(), 1.e-6, 1.e-18);
2803 SurfaceSolid handBuilt(
"handBuilt");
2805 handBuilt.SetModelTolerance(1.e-4);
2806 checkClose(handBuilt.GetModelTolerance(), 1.e-4, 1.e-18);
2807 handBuilt.SetModelTolerance(-1.);
2808 checkClose(handBuilt.GetModelTolerance(), 1.e-4, 1.e-18);
2814 SurfaceSolid v4(
"sidecarV4");
2815 BOOST_CHECK(!loadFrom(
"o2_sidecar_v4.bin", boxBytesWithHeader(4, 1.e-5), v4));
2819 std::vector<char> stump;
2820 stump.insert(stump.end(), {
'O',
'2',
'S',
'S'});
2821 appendU32(stump, 2);
2822 appendU32(stump, 6);
2823 appendU32(stump, 0);
2824 SurfaceSolid stumped(
"sidecarV2Stump");
2825 BOOST_CHECK(!loadFrom(
"o2_sidecar_v2_stump.bin", stump, stumped));
2833 constexpr double radius = 2.;
2834 constexpr double halfHeight = 3.;
2836 std::vector<char> bytes;
2837 appendSidecarHeader(bytes, 3);
2838 appendU32(bytes, 2);
2839 appendU32(bytes, 0);
2840 appendU32(bytes, 14);
2841 appendDoubles(bytes, {0., 0., 0., 0., 0., 1., 1., 0., 0., radius, -halfHeight, halfHeight, 0., 2. *
surf::kPi});
2842 appendU32(bytes, 1);
2843 appendU32(bytes, 0);
2844 appendU32(bytes, 4);
2845 const std::array<std::array<double, 4>, 4>
edges{{{0., -halfHeight, 2. *
surf::kPi, -halfHeight},
2847 {2. *
surf::kPi, halfHeight, 0., halfHeight},
2848 {0., halfHeight, 0., -halfHeight}}};
2849 for (
const auto& edge :
edges) {
2850 appendU32(bytes, 0);
2851 appendU32(bytes, 4);
2852 appendDoubles(bytes, {edge[0], edge[1], edge[2], edge[3]});
2854 appendDiskPlaneRecord(bytes, {0., 0., halfHeight}, {1., 0., 0.}, {0., 1., 0.}, radius);
2855 appendDiskPlaneRecord(bytes, {0., 0., -halfHeight}, {1., 0., 0.}, {0., -1., 0.}, radius);
2856 const auto path = writeSidecarFile(
"o2_sidecar_wiretrim_cylinder.bin", bytes);
2858 SurfaceSolid solid(
"sidecarWireCylinder");
2860 std::filesystem::remove(
path);
2866 TGeoTube
reference(
"wireTrimSidecarTube", 0., radius, halfHeight);
2867 compareContainsGrid(solid,
reference, 4., 7);
2868 compareDistance(solid,
reference, {5., 0.5, 1.}, {-1., 0., 0.});
2869 compareDistance(solid,
reference, {0., 0., 0.}, unitDirection(1., 1., 1.));
2870 checkClose(solid.Capacity(),
reference.Capacity(), 1.e-6);
2877 constexpr double majorR = 3.;
2878 constexpr double minorR = 1.;
2880 std::vector<char> bytes;
2881 appendSidecarHeader(bytes, 1);
2882 appendU32(bytes, 5);
2883 appendU32(bytes, 0);
2884 appendU32(bytes, 15);
2885 appendDoubles(bytes, {0., 0., 0., 0., 0., 1., 1., 0., 0., majorR, minorR, 0., 2. *
surf::kPi, 0., 2. *
surf::kPi});
2886 appendU32(bytes, 0);
2887 const auto path = writeSidecarFile(
"o2_sidecar_roundtrip_torus.bin", bytes);
2889 SurfaceSolid solid(
"sidecarTorus");
2891 std::filesystem::remove(
path);
2897 TGeoTorus
reference(
"sidecarTorusRef", majorR, 0., minorR);
2898 compareContainsGrid(solid,
reference, 4.5, 9);
2899 checkClose(solid.Capacity(),
reference.Capacity(), 1.e-7);
2911 constexpr double radius = 2.;
2912 constexpr double halfHeight = 3.;
2915 constexpr double joinGap = 5.e-6;
2917 const auto cylinderBytes = [&](uint32_t
version,
double modelTolerance) {
2918 std::vector<char> bytes;
2919 appendSidecarHeader(bytes, 3,
version, modelTolerance);
2920 appendU32(bytes, 2);
2921 appendU32(bytes, 0);
2922 appendU32(bytes, 14);
2923 appendDoubles(bytes, {0., 0., 0., 0., 0., 1., 1., 0., 0., radius, -halfHeight, halfHeight, 0., 2. *
surf::kPi});
2924 appendU32(bytes, 1);
2925 appendU32(bytes, 0);
2926 appendU32(bytes, 5);
2931 const std::array<std::array<double, 4>, 5>
edges{
2932 {{0., -halfHeight,
surf::kPi, -halfHeight},
2935 {2. *
surf::kPi, halfHeight, 0., halfHeight},
2936 {0., halfHeight, 0., -halfHeight}}};
2937 for (
const auto& edge :
edges) {
2938 appendU32(bytes, 0);
2939 appendU32(bytes, 4);
2940 appendDoubles(bytes, {edge[0], edge[1], edge[2], edge[3]});
2942 appendDiskPlaneRecord(bytes, {0., 0., halfHeight}, {1., 0., 0.}, {0., 1., 0.}, radius);
2943 appendDiskPlaneRecord(bytes, {0., 0., -halfHeight}, {1., 0., 0.}, {0., -1., 0.}, radius);
2946 const auto loadFrom = [](
const char*
name,
const std::vector<char>& bytes, SurfaceSolid& solid) {
2947 const auto path = writeSidecarFile(
name, bytes);
2949 std::filesystem::remove(
path);
2955 SurfaceSolid declared(
"sidecarJoinDeclared");
2956 BOOST_REQUIRE(loadFrom(
"o2_sidecar_join_declared.bin", cylinderBytes(2, 1.e-4), declared));
2958 declared.CloseShape();
2961 TGeoTube
reference(
"declaredToleranceTube", 0., radius, halfHeight);
2962 compareContainsGrid(declared,
reference, 4., 7);
2963 compareDistance(declared,
reference, {5., 0.5, 1.}, {-1., 0., 0.});
2966 SurfaceSolid silent(
"sidecarJoinSilent");
2967 BOOST_CHECK(!loadFrom(
"o2_sidecar_join_silent.bin", cylinderBytes(1, 0.), silent));
2971 SurfaceSolid tight(
"sidecarJoinTight");
2972 BOOST_CHECK(!loadFrom(
"o2_sidecar_join_tight.bin", cylinderBytes(2, 2.e-6), tight));
2981void appendBSplineEdgePlaneRecord(std::vector<char>& bytes,
const FaceFrame& frame)
2983 appendU32(bytes, 1);
2984 appendU32(bytes, 0);
2985 appendU32(bytes, 9);
2986 appendDoubles(bytes, {frame.origin[0], frame.origin[1], frame.origin[2], frame.axisU[0], frame.axisU[1],
2987 frame.axisU[2], frame.axisV[0], frame.axisV[1], frame.axisV[2]});
2988 appendU32(bytes, 1);
2989 appendU32(bytes, 0);
2990 appendU32(bytes, 4);
2991 const double extentU = frame.extentU;
2992 const double extentV = frame.extentV;
2994 appendU32(bytes, 2);
2995 appendU32(bytes, 22);
2996 appendDoubles(bytes, {3., 4.,
2997 0., 0., extentU / 3., 0., 2. * extentU / 3., 0., extentU, 0.,
2999 0., 0., 0., 0., 1., 1., 1., 1.});
3000 const std::array<std::array<double, 4>, 3> lines{
3001 {{extentU, 0., extentU, extentV}, {extentU, extentV, 0., extentV}, {0., extentV, 0., 0.}}};
3002 for (
const auto& edge : lines) {
3003 appendU32(bytes, 0);
3004 appendU32(bytes, 4);
3005 appendDoubles(bytes, {edge[0], edge[1], edge[2], edge[3]});
3017 const double w = std::sqrt(0.5);
3018 const std::vector<surf::Vec2> poles{{cu +
r, cv}, {cu +
r, cv +
r}, {cu, cv +
r}, {cu -
r, cv +
r}, {cu -
r, cv}, {cu -
r, cv -
r}, {cu, cv -
r}, {cu +
r, cv -
r}, {cu +
r, cv}};
3019 return surf::Curve2D::makeBSpline(2, poles, {1.,
w, 1.,
w, 1.,
w, 1.,
w, 1.},
3020 {0., 0., 0., 1., 1., 2., 2., 3., 3., 4., 4., 4.});
3023surf::Curve2D quarterCircleBSpline(
double cu,
double cv,
double r,
double a0)
3026 const double aMid = 0.5 * (a0 + a1);
3027 const std::vector<surf::Vec2> poles{{cu +
r * std::cos(a0), cv +
r * std::sin(a0)},
3028 {cu +
r * std::sqrt(2.) * std::cos(aMid), cv +
r * std::sqrt(2.) * std::sin(aMid)},
3029 {cu +
r * std::cos(a1), cv +
r * std::sin(a1)}};
3030 return surf::Curve2D::makeBSpline(2, poles, {1., std::sqrt(0.5), 1.}, {0., 0., 0., 1., 1., 1.});
3046 WireStatus status = WireStatus::Valid;
3047 BOOST_REQUIRE(exactWire.
initialize({surf::Curve2D::makeCircle({0., 0.}, 1.)}, WireRole::Outer, status));
3052 BOOST_REQUIRE(splineWire.
initialize({fullCircleBSpline(0., 0., 1.)}, WireRole::Outer, status));
3059 BOOST_CHECK(splineWire.
classify({1. - justInsideTheBand, 0.}) == WireClassification::Boundary);
3060 BOOST_CHECK(splineWire.
classify({1. + justInsideTheBand, 0.}) == WireClassification::Boundary);
3065 BOOST_CHECK(exactWire.classify({1. - justInsideTheBand, 0.}) == WireClassification::Inside);
3070 BOOST_REQUIRE(bigCylinder.
initialize({0., 0., 0.}, {0., 0., 1.}, {1., 0., 0.}, 100., -1., 1., 0.,
surf::kTwoPi,
3074 BOOST_REQUIRE(squareWire.
initialize({surf::Curve2D::makeLine({0., -1.}, {1., -1.}),
3075 surf::Curve2D::makeLine({1., -1.}, {1., 1.}),
3076 surf::Curve2D::makeLine({1., 1.}, {0., 1.}),
3077 surf::Curve2D::makeLine({0., 1.}, {0., -1.})},
3078 WireRole::Outer, status, bigMetric));
3080 BOOST_CHECK(squareWire.classify({0.5, 1. - 1.e-10}, bigMetric) == WireClassification::Inside);
3082 BOOST_CHECK(squareWire.classify({0.5, 1. - 1.e-12}, bigMetric) == WireClassification::Boundary);
3087 const double seamDrift = 4.e-7;
3092 second.
poles.front() = {second.
poles.front().uCoord + seamDrift, second.
poles.front().vCoord};
3094 BOOST_REQUIRE(driftedWire.
initialize({first, second, third, fourth}, WireRole::Outer, status));
3095 for (
const auto& curve : driftedWire.curves) {
3100 const Vec2 thisEnd = driftedWire.
curves[
index].loopEnd();
3101 const Vec2 nextStart = driftedWire.
curves[(
index + 1) % driftedWire.
curves.size()].loopStart();
3106 for (
const auto& curve : driftedWire.curves) {
3107 const auto& polyline = curve.bsplineSamples();
3108 BOOST_REQUIRE(polyline.size() >= 2);
3144 const double trimRadius = 0.5;
3146 BOOST_REQUIRE(splineTrimmed.
initialize({0., 0., 0.}, {0., 0., 1.}, {1., 0., 0.}, 2., -1., 1., 0.,
surf::kTwoPi,
3148 {quarterCircleBSpline(centrePhi, 0., trimRadius, 0.),
3149 quarterCircleBSpline(centrePhi, 0., trimRadius, surf::kHalfPi),
3150 quarterCircleBSpline(centrePhi, 0., trimRadius, surf::kPi),
3151 quarterCircleBSpline(centrePhi, 0., trimRadius, 3. * surf::kHalfPi)},
3156 BOOST_REQUIRE(arcTrimmed.
initialize({0., 0., 0.}, {0., 0., 1.}, {1., 0., 0.}, 2., -1., 1., 0.,
surf::kTwoPi,
false,
3157 {Curve2D::makeCircle({centrePhi, 0.}, trimRadius)}, {},
error));
3161 std::vector<surf::RayHit> hits;
3162 surface.
appendIntersections({0., 0., 0.}, {std::cos(phi), std::sin(phi), 0.}, 0., 1.e30, hits);
3169 const double justInside = 0.5 * band;
3173 const auto hits = hitAt(splineTrimmed, centrePhi);
3174 BOOST_REQUIRE_EQUAL(hits.size(), 1u);
3179 const auto hits = hitAt(splineTrimmed, centrePhi + trimRadius - justInside);
3180 BOOST_REQUIRE_EQUAL(hits.size(), 1u);
3186 const auto hits = hitAt(splineTrimmed, centrePhi + trimRadius + justInside);
3187 BOOST_REQUIRE_EQUAL(hits.size(), 1u);
3192 BOOST_CHECK(hitAt(splineTrimmed, centrePhi + trimRadius + 100. * band).
empty());
3196 const auto inside = hitAt(arcTrimmed, centrePhi + trimRadius - justInside);
3197 BOOST_REQUIRE_EQUAL(inside.size(), 1u);
3199 BOOST_CHECK(hitAt(arcTrimmed, centrePhi + trimRadius + justInside).
empty());
3204 BOOST_REQUIRE(plain.
initialize({0., 0., 0.}, {0., 0., 1.}, {1., 0., 0.}, 2., -1., 1., 0.,
surf::kTwoPi,
false,
3206 const auto hits = hitAt(plain, centrePhi);
3207 BOOST_REQUIRE_EQUAL(hits.size(), 1u);
3216 constexpr double halfX = 1.;
3217 constexpr double halfY = 2.;
3218 constexpr double halfZ = 3.;
3220 std::vector<char> bytes;
3221 appendSidecarHeader(bytes, 6);
3222 appendBSplineEdgePlaneRecord(bytes, boxFaceFrame(0, halfX, halfY, halfZ));
3224 appendPlaneRecord(bytes, boxFaceFrame(
faceIndex, halfX, halfY, halfZ));
3226 const auto path = writeSidecarFile(
"o2_sidecar_bspline_box.bin", bytes);
3228 SurfaceSolid
box(
"sidecarBSplineBox");
3230 std::filesystem::remove(
path);
3238 compareDistance(
box,
reference, {5., 0.5, 0.5}, {-1., 0., 0.});
3239 compareDistance(
box,
reference, {0., 0., 0.}, unitDirection(1., 1., 1.));
3240 checkClose(
box.Capacity(),
reference.Capacity(), 1.e-6);
3249 const auto onCylinder = [](
double phi,
double height) {
3250 return Vec3{2. * std::cos(phi), 2. * std::sin(phi),
height};
3257 const double trimRadius = 0.5;
3259 const std::vector<Curve2D> bsplineOuter{quarterCircleBSpline(centrePhi, 0., trimRadius, 0.),
3260 quarterCircleBSpline(centrePhi, 0., trimRadius,
surf::kHalfPi),
3261 quarterCircleBSpline(centrePhi, 0., trimRadius,
surf::kPi),
3262 quarterCircleBSpline(centrePhi, 0., trimRadius, 3. *
surf::kHalfPi)};
3263 BOOST_REQUIRE(bsplineDisk.
initialize({0., 0., 0.}, {0., 0., 1.}, {1., 0., 0.}, 2., -1., 1., 0.,
surf::kTwoPi,
false,
3264 bsplineOuter, {}, error));
3268 BOOST_REQUIRE(arcDisk.
initialize({0., 0., 0.}, {0., 0., 1.}, {1., 0., 0.}, 2., -1., 1., 0.,
surf::kTwoPi,
false,
3269 {Curve2D::makeCircle({centrePhi, 0.}, trimRadius)}, {},
error));
3275 for (
int phiStep = -12; phiStep <= 12; ++phiStep) {
3276 const double phi = centrePhi + 0.09 * phiStep;
3277 for (
int hStep = -12; hStep <= 12; ++hStep) {
3278 const double height = 0.09 * hStep;
3280 const double distToCentre = std::hypot(phi - centrePhi,
height);
3281 if (std::abs(distToCentre - trimRadius) < 5.e-3) {
3284 const Vec3 point = onCylinder(phi,
height);
3285 BOOST_CHECK_EQUAL(bsplineDisk.containsPointOnSurface(point), arcDisk.containsPointOnSurface(point));
3289 BOOST_CHECK_GT(compared, 100);
3293 std::vector<surf::RayHit> hits;
3294 bsplineDisk.appendIntersections({0., 0., 0.}, {std::cos(centrePhi), std::sin(centrePhi), 0.}, 0., 1.e30, hits);
3295 BOOST_REQUIRE_EQUAL(hits.size(), 1u);
3296 checkClose(hits.front().distance, 2.);
3298 bsplineDisk.appendIntersections({0., 0., 0.}, {std::cos(0.), std::sin(0.), 0.}, 0., 1.e30, hits);
3325 const auto onCylinder = [](
double phi,
double height) {
3326 return Vec3{2. * std::cos(phi), 2. * std::sin(phi),
height};
3330 const std::vector<Curve2D> outer{Curve2D::makeLine({0., -3.}, {
surf::kTwoPi, -3.}),
3333 Curve2D::makeLine({0., 3.}, {0., -3.})};
3337 const double trimRadius = 0.5;
3338 const std::vector<Curve2D> bsplineHole{quarterCircleBSpline(centrePhi, 0., trimRadius, 0.),
3339 quarterCircleBSpline(centrePhi, 0., trimRadius,
surf::kHalfPi),
3340 quarterCircleBSpline(centrePhi, 0., trimRadius,
surf::kPi),
3341 quarterCircleBSpline(centrePhi, 0., trimRadius, 3. *
surf::kHalfPi)};
3344 BOOST_REQUIRE(bsplineHoled.
initialize({0., 0., 0.}, {0., 0., 1.}, {1., 0., 0.}, 2., -3., 3., 0.,
surf::kTwoPi,
false,
3345 outer, {bsplineHole}, error));
3347 BOOST_REQUIRE(arcHoled.
initialize({0., 0., 0.}, {0., 0., 1.}, {1., 0., 0.}, 2., -3., 3., 0.,
surf::kTwoPi,
false,
3348 outer, {{Curve2D::makeCircle({centrePhi, 0.}, trimRadius)}},
error));
3351 BOOST_CHECK(!arcHoled.containsPointOnSurface(onCylinder(centrePhi, 0.)));
3352 BOOST_CHECK(!bsplineHoled.containsPointOnSurface(onCylinder(centrePhi, 0.)));
3354 BOOST_CHECK(bsplineHoled.containsPointOnSurface(onCylinder(0.5, 0.)));
3355 BOOST_CHECK(bsplineHoled.containsPointOnSurface(onCylinder(centrePhi, 2.5)));
3359 for (
int phiStep = -12; phiStep <= 12; ++phiStep) {
3360 const double phi = centrePhi + 0.09 * phiStep;
3361 for (
int hStep = -12; hStep <= 12; ++hStep) {
3362 const double height = 0.09 * hStep;
3363 if (std::abs(std::hypot(phi - centrePhi,
height) - trimRadius) < 5.e-3) {
3366 const Vec3 point = onCylinder(phi,
height);
3367 BOOST_CHECK_EQUAL(bsplineHoled.containsPointOnSurface(point), arcHoled.containsPointOnSurface(point));
3371 BOOST_CHECK_GT(compared, 100);
3374 std::vector<surf::RayHit> hits;
3375 bsplineHoled.appendIntersections({0., 0., 0.}, {std::cos(centrePhi), std::sin(centrePhi), 0.}, 0., 1.e30, hits);
3378 bsplineHoled.appendIntersections({0., 0., 0.}, {std::cos(0.5), std::sin(0.5), 0.}, 0., 1.e30, hits);
3385 BOOST_REQUIRE(singleEdgeHoled.
initialize({0., 0., 0.}, {0., 0., 1.}, {1., 0., 0.}, 2., -3., 3., 0.,
surf::kTwoPi,
3386 false, outer, {{fullCircleBSpline(centrePhi, 0., trimRadius)}},
error));
3387 BOOST_CHECK(!singleEdgeHoled.containsPointOnSurface(onCylinder(centrePhi, 0.)));
3388 BOOST_CHECK(singleEdgeHoled.containsPointOnSurface(onCylinder(0.5, 0.)));
3389 BOOST_CHECK(singleEdgeHoled.containsPointOnSurface(onCylinder(centrePhi, 2.5)));
3390 for (
int phiStep = -12; phiStep <= 12; ++phiStep) {
3391 const double phi = centrePhi + 0.09 * phiStep;
3392 for (
int hStep = -12; hStep <= 12; ++hStep) {
3393 const double height = 0.09 * hStep;
3394 if (std::abs(std::hypot(phi - centrePhi,
height) - trimRadius) < 5.e-3) {
3397 const Vec3 point = onCylinder(phi,
height);
3398 BOOST_CHECK_EQUAL(singleEdgeHoled.containsPointOnSurface(point), arcHoled.containsPointOnSurface(point));
3406BoundaryCurve quarterCircleBoundaryCurve(
double cu,
double cv,
double r,
double a0)
3409 const double aMid = 0.5 * (a0 + a1);
3410 const std::vector<Point2D> poles{{cu +
r * std::cos(a0), cv +
r * std::sin(a0)},
3411 {cu +
r * std::sqrt(2.) * std::cos(aMid), cv +
r * std::sqrt(2.) * std::sin(aMid)},
3412 {cu +
r * std::cos(a1), cv +
r * std::sin(a1)}};
3413 return BoundaryCurve::makeBSpline(2, poles, {1., std::sqrt(0.5), 1.}, {0., 0., 0., 1., 1., 1.});
3421void checkSolidsIdentical(
const SurfaceSolid& solid,
const SurfaceSolid&
other,
double extent,
int samples)
3429 static_cast<int>(solid.GetNavigationReliability()));
3437 for (
int dimension = 0; dimension < 3; ++dimension) {
3442 for (
const auto& point : probeGrid(extent,
samples)) {
3443 BOOST_TEST_CONTEXT(
"point = (" << point[0] <<
", " << point[1] <<
", " << point[2] <<
")")
3447 solid.Safety(point.data(), solid.Contains(point.data())));
3448 for (
const auto& direction : probeDirections()) {
3450 solid.DistFromOutside(point.data(), direction.data(), 3));
3452 solid.DistFromInside(point.data(), direction.data(), 3));
3459std::unique_ptr<SurfaceSolid> writeAndReadBack(
const SurfaceSolid& solid)
3461 const auto path = std::filesystem::temp_directory_path() /
3462 (std::string(
"o2_bvhsurfacesolid_persist_") + solid.GetName() +
".root");
3464 TFile
file(
path.string().c_str(),
"RECREATE");
3465 BOOST_REQUIRE(!
file.IsZombie());
3467 file.WriteObject(
const_cast<SurfaceSolid*
>(&solid),
"solid");
3469 std::unique_ptr<SurfaceSolid> restored;
3471 TFile
file(
path.string().c_str(),
"READ");
3472 BOOST_REQUIRE(!
file.IsZombie());
3473 restored.reset(
file.Get<SurfaceSolid>(
"solid"));
3475 std::filesystem::remove(
path);
3494 const auto box = makeBoxSolid(
"persistBox", 1., 2., 3.);
3495 const auto tube = makeTubeSolid(
"persistTube", 1., 2., 3.);
3496 const auto cone = makeConeSolid(
"persistCone", 2., 1., 3.);
3497 const auto sphere = makeSphereSolid(
"persistSphere", 2.);
3498 const auto torus = makeTorusSolid(
"persistTorus", 3., 1.);
3499 const auto capsule = makeCapsuleSolid(
"persistCapsule", 2., 3.);
3501 for (
const auto* solid : {
box.get(), tube.get(), cone.get(), sphere.get(), torus.get(), capsule.get()}) {
3502 BOOST_TEST_CONTEXT(
"solid = " << solid->GetName())
3504 const auto restored = writeAndReadBack(*solid);
3505 BOOST_REQUIRE(restored !=
nullptr);
3506 checkSolidsIdentical(*solid, *restored, 4.5, 5);
3512 SurfaceSolid trimmed(
"persistTrimmed");
3513 constexpr double radius = 2.;
3514 constexpr double halfHeight = 3.;
3515 const std::vector<BoundaryCurve> window{quarterCircleBoundaryCurve(
surf::kPi, 0., 0.5, 0.),
3519 BOOST_REQUIRE(trimmed.AddCylindricalSurface({0., 0., 0.}, {0., 0., 1.}, {1., 0., 0.}, radius, -halfHeight,
3521 BOOST_REQUIRE(addDiskSurface(trimmed, {0., 0., halfHeight}, {1., 0., 0.}, {0., 1., 0.}, radius));
3522 BOOST_REQUIRE(addDiskSurface(trimmed, {0., 0., -halfHeight}, {1., 0., 0.}, {0., -1., 0.}, radius));
3523 trimmed.CloseShape(
false);
3526 const auto restoredTrimmed = writeAndReadBack(trimmed);
3527 BOOST_REQUIRE(restoredTrimmed !=
nullptr);
3528 checkSolidsIdentical(trimmed, *restoredTrimmed, 4.5, 5);
3532 SurfaceSolid toleranced(
"persistTolerance");
3534 BOOST_REQUIRE(addBoxFace(toleranced,
faceIndex, 1., 2., 3.));
3536 toleranced.SetModelTolerance(7.25e-5);
3537 toleranced.CloseShape(
false);
3538 const auto restoredToleranced = writeAndReadBack(toleranced);
3539 BOOST_REQUIRE(restoredToleranced !=
nullptr);
3540 checkClose(restoredToleranced->GetModelTolerance(), 7.25e-5, 1.e-18);
3544 SurfaceSolid openBox(
"persistOpenBox");
3546 BOOST_REQUIRE(addBoxFace(openBox,
faceIndex, 1., 2., 3.));
3548 openBox.CloseShape(
false);
3549 BOOST_REQUIRE(!openBox.IsNavigable());
3551 static_cast<int>(SurfaceSolid::NavigationReliability::OpenSurfaceSet));
3553 const auto restoredOpenBox = writeAndReadBack(openBox);
3554 BOOST_REQUIRE(restoredOpenBox !=
nullptr);
3556 checkSolidsIdentical(openBox, *restoredOpenBox, 4.5, 5);
3565 SurfaceSolid
empty(
"emptySolid");
3567 static_cast<int>(SurfaceSolid::NavigationReliability::Undetermined));
3571 empty.CloseShape(
false);
3574 static_cast<int>(SurfaceSolid::NavigationReliability::Undetermined));
3584std::vector<std::array<double, 3>> spiralDirections(
int count)
3586 std::vector<std::array<double, 3>> directions;
3587 directions.reserve(
count);
3589 const double cosTheta = 1. - 2. * (
index + 0.5) /
count;
3590 const double sinTheta = std::sqrt(1. - cosTheta * cosTheta);
3591 const double phi = 2.399963229728653 *
index;
3592 directions.push_back({
sinTheta * std::cos(phi),
sinTheta * std::sin(phi), cosTheta});
3608 const auto box = makeBoxSolid(
"dirBox", 1., 2., 3.);
3609 const auto tube = makeTubeSolid(
"dirTube", 1., 2., 3.);
3610 const auto cone = makeConeSolid(
"dirCone", 2., 1., 3.);
3611 const auto sphere = makeSphereSolid(
"dirSphere", 2.);
3612 const auto torus = makeTorusSolid(
"dirTorus", 3., 1.);
3613 const auto capsule = makeCapsuleSolid(
"dirCapsule", 2., 3.);
3615 const auto directions = spiralDirections(13);
3616 for (
const auto* solid : {
box.get(), tube.get(), cone.get(), sphere.get(), torus.get(), capsule.get()}) {
3617 BOOST_TEST_CONTEXT(
"solid = " << solid->GetName())
3619 BOOST_REQUIRE(solid->IsNavigable());
3620 for (
const auto& point : probeGrid(4.5, 7)) {
3621 const bool reference = solid->Contains(point.data());
3622 for (
const auto& direction : directions) {
3623 BOOST_TEST_CONTEXT(
"point = (" << point[0] <<
", " << point[1] <<
", " << point[2] <<
") direction = ("
3624 << direction[0] <<
", " << direction[1] <<
", " << direction[2] <<
")")
3647 constexpr double halfX = 1.;
3648 constexpr double halfY = 2.;
3649 constexpr double halfZ = 3.;
3650 constexpr double gap = 0.05;
3652 SurfaceSolid gapped(
"gappedBox");
3654 BOOST_REQUIRE(addBoxFace(gapped,
faceIndex, halfX, halfY, halfZ));
3658 BOOST_REQUIRE(gapped.AddPlanarSurface({halfX, -halfY, -halfZ}, {0., 1., 0.}, {0., 0., 1.},
3659 rectangleWire(2. * halfY, halfZ - gap)));
3660 BOOST_REQUIRE(gapped.AddPlanarSurface({halfX, -halfY, gap}, {0., 1., 0.}, {0., 0., 1.},
3661 rectangleWire(2. * halfY, halfZ - gap)));
3662 gapped.CloseShape(
false);
3665 BOOST_REQUIRE(!gapped.IsNavigable());
3667 static_cast<int>(SurfaceSolid::NavigationReliability::OpenSurfaceSet));
3670 const std::array<double, 3> insidePoint{0., 0.3, 0.};
3671 const std::array<double, 3> throughGap{1., 0., 0.};
3672 BOOST_CHECK(!gapped.ContainsAlongDirection(insidePoint.data(), throughGap.data()));
3674 BOOST_CHECK(gapped.Contains_Loop(insidePoint.data()));
3678 const auto intact = makeBoxSolid(
"intactBox", halfX, halfY, halfZ);
3679 BOOST_CHECK(intact->Contains(insidePoint.data()));
3680 BOOST_CHECK(intact->ContainsAlongDirection(insidePoint.data(), throughGap.data()));
3684 for (
const auto& point : probeGrid(4.5, 7)) {
3685 BOOST_TEST_CONTEXT(
"point = (" << point[0] <<
", " << point[1] <<
", " << point[2] <<
")")
3687 BOOST_CHECK_EQUAL(gapped.Contains(point.data()), gapped.Contains_Loop(point.data()));
3700std::unique_ptr<SurfaceSolid> makeLPrismSolid(
const char*
name,
double height = 1.)
3703 const std::vector<Point2D> footprint{{0., 0.}, {3., 0.}, {3., 1.}, {1., 1.}, {1., 2.}, {0., 2.}};
3705 auto solid = std::make_unique<SurfaceSolid>(
name);
3708 std::vector<Point2D> bottomWire;
3709 bottomWire.reserve(footprint.size());
3710 for (
const auto&
vertex : footprint) {
3713 BOOST_REQUIRE(solid->AddPlanarSurface({0., 0., 0.}, {0., 1., 0.}, {1., 0., 0.}, bottomWire));
3714 BOOST_REQUIRE(solid->AddPlanarSurface({0., 0., height}, {1., 0., 0.}, {0., 1., 0.}, footprint));
3720 const auto&
end = footprint[(
index + 1) % footprint.size()];
3721 const double deltaU =
end[0] -
start[0];
3722 const double deltaV =
end[1] -
start[1];
3723 const double length = std::hypot(deltaU, deltaV);
3724 BOOST_REQUIRE(solid->AddPlanarSurface({start[0], start[1], 0.}, {deltaU / length, deltaV / length, 0.},
3727 solid->CloseShape();
3738 SurfaceSolid
box(
"preCloseBox");
3739 addBoxSurfaces(
box, 1., 2., 3.);
3740 BOOST_REQUIRE(!
box.HasBVH());
3742 for (
const auto& point : probeGrid(4.5, 5)) {
3743 BOOST_TEST_CONTEXT(
"point = (" << point[0] <<
", " << point[1] <<
", " << point[2] <<
")")
3745 const bool inside = std::abs(point[0]) < 1. && std::abs(point[1]) < 2. && std::abs(point[2]) < 3.;
3758 const auto box = makeBoxSolid(
"boundaryPolicyBox", 1., 2., 3.);
3759 const std::array<double, 3> onFace{1., 0.5, 0.5};
3760 const std::array<double, 3> outward{1., 0., 0.};
3761 const std::array<double, 3> inward{-1., 0., 0.};
3768 checkClose(
box->DistFromInside(onFace.data(), outward.data(), 3),
3769 reference.DistFromInside(onFace.data(), outward.data(), 3));
3770 checkClose(
box->DistFromOutside(onFace.data(), inward.data(), 3),
3771 reference.DistFromOutside(onFace.data(), inward.data(), 3));
3774 checkClose(
box->DistFromInside(onFace.data(), inward.data(), 3), 2.);
3777 for (
const auto& direction : {outward, inward}) {
3778 checkDistanceAgainstLoop(*
box, onFace, direction);
3794 const double invSqrt2 = 1. / std::sqrt(2.);
3797 const auto box = makeBoxSolid(
"grazeBox", 1., 2., 3.);
3798 const std::array<double, 3> grazeDirection{invSqrt2, -invSqrt2, 0.};
3799 const std::array<double, 3> grazeOrigin{1. - 5. * invSqrt2, 2. + 5. * invSqrt2, 0.};
3800 BOOST_REQUIRE(!
box->Contains(grazeOrigin.data()));
3803 const std::array<double, 3> pastTouch{1. + 1.e-3 * invSqrt2, 2. - 1.e-3 * invSqrt2, 0.};
3804 BOOST_REQUIRE(!
box->Contains(pastTouch.data()));
3806 BOOST_CHECK_EQUAL(
box->DistFromOutside(grazeOrigin.data(), grazeDirection.data(), 3), TGeoShape::Big());
3807 BOOST_CHECK_EQUAL(
box->DistFromOutside_Loop(grazeOrigin.data(), grazeDirection.data()), TGeoShape::Big());
3810 const auto prism = makeLPrismSolid(
"grazePrism");
3811 BOOST_REQUIRE(prism->IsNavigable());
3812 checkClose(prism->Capacity(), 4.);
3816 const std::array<double, 3> reflexDirection{invSqrt2, -invSqrt2, 0.};
3817 const std::array<double, 3> reflexOrigin{1. - 0.5 * invSqrt2, 1. + 0.5 * invSqrt2, 0.5};
3818 BOOST_REQUIRE(prism->Contains(reflexOrigin.data()));
3819 const std::array<double, 3> pastReflex{1. + 1.e-3 * invSqrt2, 1. - 1.e-3 * invSqrt2, 0.5};
3820 BOOST_REQUIRE(prism->Contains(pastReflex.data()));
3822 const double touchDistance = 0.5;
3823 const double exitDistance = 0.5 + std::sqrt(2.);
3824 const double reported = prism->DistFromInside(reflexOrigin.data(), reflexDirection.data(), 3);
3825 BOOST_CHECK_GT(reported, touchDistance + 1.e-6);
3826 checkClose(reported, exitDistance);
3827 BOOST_CHECK_EQUAL(prism->DistFromInside_Loop(reflexOrigin.data(), reflexDirection.data()), reported);
3833 const auto box = makeBoxSolid(
"describeBox", 1., 2., 3.);
3834 const SurfaceSolid::Point3D inside{0.2, 0.3, 0.4};
3835 const double invSqrt2 = 1. / std::sqrt(2.);
3837 const std::vector<std::pair<SurfaceSolid::Point3D, double>> cases{
3838 {{1., 0., 0.}, 0.8}, {{invSqrt2, -invSqrt2, 0.}, 0.8 * std::sqrt(2.)}, {{0., 0., -2.}, 3.4}};
3839 for (
const auto& [direction, exitDistance] : cases) {
3840 BOOST_TEST_CONTEXT(
"direction = (" << direction[0] <<
", " << direction[1] <<
", " << direction[2] <<
")")
3842 std::vector<SurfaceSolid::ContainsCrossing> bvhCrossings;
3843 std::vector<SurfaceSolid::ContainsCrossing> loopCrossings;
3844 box->DescribeContainsCrossings(inside, direction, bvhCrossings, loopCrossings);
3846 BOOST_REQUIRE_EQUAL(bvhCrossings.size(), 1u);
3847 BOOST_REQUIRE_EQUAL(loopCrossings.size(), 1u);
3849 checkClose(bvhCrossings[0].
distance, exitDistance);
3850 BOOST_CHECK_GT(bvhCrossings[0].normalAlignment, 0.);
3852 BOOST_CHECK(
box->ContainsAlongDirection(inside.data(), direction.data()));
3862 const auto prism = makeLPrismSolid(
"sweepPrism");
3863 BOOST_REQUIRE(prism->IsNavigable());
3865 sweepDistanceAgainstLoop(*prism, 3.5, 5);
3867 const auto directions = spiralDirections(13);
3868 for (
const auto& point : probeGrid(3.5, 7)) {
3869 const bool inside = prism->Contains(point.data());
3870 BOOST_TEST_CONTEXT(
"point = (" << point[0] <<
", " << point[1] <<
", " << point[2] <<
")")
3873 for (
const auto& direction : directions) {
3874 BOOST_CHECK_EQUAL(prism->ContainsAlongDirection(point.data(), direction.data()), inside);
3877 const bool expected = point[2] > 0. && point[2] < 1. &&
3878 ((point[0] > 0. && point[0] < 3. && point[1] > 0. && point[1] < 1.) ||
3879 (point[0] > 0. && point[0] < 1. && point[1] >= 1. && point[1] < 2.));
3899 const std::vector<Vec2> poles{{0., 0.}, {1., 2.}, {3., 2.}, {4., 0.}};
3900 const Curve2D clamped = Curve2D::makeBSpline(3, poles, {}, {0., 0., 0., 0., 1., 1., 1., 1.});
3901 const Curve2D uniform = Curve2D::makeBSpline(3, poles, {}, {0., 1., 2., 3., 4., 5., 6., 7.});
3903 BOOST_REQUIRE(clamped.valid());
3904 BOOST_REQUIRE(uniform.valid());
3907 for (
const auto* curve : {&clamped, &uniform}) {
3908 const Vec2
start = curve->startPoint();
3909 const Vec2
end = curve->endPoint();
3910 const Vec2 evaluatedStart = curve->pointAt(0.);
3911 const Vec2 evaluatedEnd = curve->pointAt(1.);
3912 checkClose(
start.uCoord, evaluatedStart.uCoord);
3913 checkClose(
start.vCoord, evaluatedStart.vCoord);
3914 checkClose(
end.uCoord, evaluatedEnd.uCoord);
3915 checkClose(
end.vCoord, evaluatedEnd.vCoord);
3920 checkClose(clamped.startPoint().uCoord, 0.);
3921 checkClose(clamped.endPoint().uCoord, 4.);
3922 BOOST_CHECK_GT(std::hypot(uniform.startPoint().uCoord - poles.front().uCoord,
3923 uniform.startPoint().vCoord - poles.front().vCoord),
3928 const Vec2 uniformStart = uniform.startPoint();
3929 const Vec2 uniformEnd = uniform.endPoint();
3932 BOOST_CHECK(wire.
initialize({uniform, Curve2D::makeLine(uniformEnd, uniformStart)}, surf::WireRole::Outer, status));
3950 const double sweep = 350. *
surf::kPi / 180.;
3951 const double overshoot = 0.4;
3952 const std::vector<Curve2D> outer{
3953 Curve2D::makeBSpline(2, {{0., -1.}, {-overshoot, 0.}, {0.5 * sweep, 1.}}, {}, {0., 0., 0., 1., 1., 1.}),
3954 Curve2D::makeBSpline(2, {{0.5 * sweep, 1.}, {sweep + overshoot, 0.}, {sweep, -1.}}, {}, {0., 0., 0., 1., 1., 1.}),
3955 Curve2D::makeLine({sweep, -1.}, {0., -1.})};
3958 Vec2 hullLower{1.e300, 1.e300};
3959 Vec2 hullUpper{-1.e300, -1.e300};
3962 BOOST_REQUIRE(hullWire.
initialize(outer, surf::WireRole::Outer, hullStatus));
3964 BOOST_REQUIRE_GT(hullUpper.uCoord - hullLower.uCoord,
surf::kTwoPi);
3967 Vec2 tightLower{1.e300, 1.e300};
3968 Vec2 tightUpper{-1.e300, -1.e300};
3970 BOOST_CHECK_LT(tightUpper.uCoord - tightLower.uCoord,
surf::kTwoPi);
3974 BOOST_CHECK_MESSAGE(surface.
initialize({0., 0., 0.}, {0., 0., 1.}, {1., 0., 0.}, 2., -2., 2., 0.,
surf::kTwoPi,
3975 false, outer, {}, error),
3976 "near-full-turn trim rejected: " << error);
3980 const std::vector<Curve2D> overWrapped{Curve2D::makeLine({0., -1.}, {
surf::kTwoPi + 0.5, -1.}),
3983 Curve2D::makeLine({0., 1.}, {0., -1.})};
3984 BOOST_CHECK(!tooWide.
initialize({0., 0., 0.}, {0., 0., 1.}, {1., 0., 0.}, 2., -2., 2., 0.,
surf::kTwoPi,
false,
3985 overWrapped, {}, error));
3996 SurfaceSolid solid(
"rejectingSolid");
3997 BOOST_REQUIRE(addBoxFace(solid, 0, 1., 2., 3.));
3998 BOOST_REQUIRE_EQUAL(solid.GetNsurfaces(), 1);
4002 BOOST_CHECK(!solid.AddPlanarSurface({0., 0., 0.}, {1., 0., 0.}, {1., 0., 0.}, rectangleWire(1., 1.)));
4003 BOOST_CHECK(!solid.AddPlanarSurface({0., 0., 0.}, {1., 0., 0.}, {0., 1., 0.}, {{0., 0.}, {1., 0.}}));
4004 BOOST_CHECK(!solid.AddCylindricalSurface({0., 0., 0.}, {0., 0., 1.}, {1., 0., 0.}, 0., -1., 1.));
4009 std::vector<char> bytes;
4010 appendSidecarHeader(bytes, 1);
4011 appendU32(bytes, 2);
4012 appendU32(bytes, 0);
4013 appendU32(bytes, 14);
4014 appendDoubles(bytes, {0., 0., 0., 0., 0., 1., 1., 0., 0., 0. , -1., 1., 0., 2. *
surf::kPi});
4015 appendU32(bytes, 0);
4016 const auto path = writeSidecarFile(
"o2_sidecar_rejected_face.bin", bytes);
4017 SurfaceSolid loaded(
"loadedRejecting");
4019 std::filesystem::remove(
path);
4034 constexpr double halfX = 1.;
4035 constexpr double halfY = 1.5;
4036 constexpr double halfZ = 2.;
4039 SurfaceSolid closedBox(
"rimClosedBox");
4040 addBoxSurfaces(closedBox, halfX, halfY, halfZ);
4041 closedBox.CloseShape(
false);
4042 BOOST_REQUIRE(closedBox.IsNavigable());
4046 BOOST_CHECK_SMALL(closedBox.GetMaxRimIsolation(), 1.e-12);
4047 BOOST_CHECK_SMALL(closedBox.GetUnmatchedRimLength(), 1.e-12);
4049 BOOST_CHECK_SMALL(closedBox.GetRimChordResolution(), 1.e-12);
4051 BOOST_CHECK_CLOSE(closedBox.GetTotalRimLength(), 16. * (halfX + halfY + halfZ), 1.e-9);
4056 constexpr double delta = 1.e-3;
4057 SurfaceSolid shiftedBox(
"rimShiftedBox");
4060 BOOST_REQUIRE(addBoxFace(shiftedBox,
faceIndex, halfX, halfY, halfZ,
false,
center));
4062 shiftedBox.CloseShape(
false);
4063 BOOST_CHECK_CLOSE(shiftedBox.GetMaxRimIsolation(), delta, 1.e-6);
4065 BOOST_CHECK(shiftedBox.GetMaxRimIsolation() > shiftedBox.GetRimChordResolution());
4075 constexpr double halfX = 1.;
4076 constexpr double halfY = 1.5;
4077 constexpr double halfZ = 2.;
4079 SurfaceSolid resampled(
"rimResampledBox");
4081 BOOST_REQUIRE(addBoxFace(resampled,
faceIndex, halfX, halfY, halfZ));
4084 const FaceFrame frame = boxFaceFrame(5, halfX, halfY, halfZ);
4085 const double extentU = frame.extentU;
4086 const double extentV = frame.extentV;
4087 BOOST_REQUIRE(resampled.AddPlanarSurface(frame.origin, frame.axisU, frame.axisV,
4089 {0.5 * extentU, 0.},
4091 {extentU, 0.5 * extentV},
4093 {0.5 * extentU, extentV},
4095 {0., 0.5 * extentV}}));
4096 resampled.CloseShape(
false);
4101 BOOST_CHECK(resampled.GetBoundaryEdgeCount() > 0);
4109 BOOST_CHECK_SMALL(resampled.GetMaxRimIsolation(), 1.e-12);
4110 BOOST_CHECK_SMALL(resampled.GetUnmatchedRimLength(), 1.e-12);
4118 constexpr double halfX = 1.;
4119 constexpr double halfY = 1.5;
4120 constexpr double halfZ = 2.;
4121 constexpr double delta = 1.e-3;
4123 SurfaceSolid shiftedBox(
"rimReportBox");
4126 BOOST_REQUIRE(addBoxFace(shiftedBox,
faceIndex, halfX, halfY, halfZ,
false,
center));
4128 shiftedBox.CloseShape(
false);
4130 const auto& rims = shiftedBox.GetRimReports();
4131 BOOST_REQUIRE_EQUAL(
static_cast<int>(rims.size()), shiftedBox.GetRimCount());
4133 int boundaryRims = 0;
4134 auto worst = SurfaceSolid::NavigationReliability::Reliable;
4135 double openLength = 0.;
4136 for (
const auto& rim : rims) {
4137 BOOST_CHECK(rim.surface >= 0 && rim.surface < shiftedBox.GetNsurfaces());
4141 openLength += rim.unmatchedLength;
4142 if (rim.state == SurfaceSolid::NavigationReliability::OpenSurfaceSet) {
4149 worst = std::max(worst, rim.state);
4152 BOOST_CHECK(worst == shiftedBox.GetNavigationReliability());
4153 BOOST_CHECK_CLOSE(openLength, shiftedBox.GetUnmatchedRimLength(), 1.e-9);
4156 const auto lifted = std::find_if(rims.begin(), rims.end(), [](
const auto& rim) { return rim.surface == 4; });
4157 BOOST_REQUIRE(lifted != rims.end());
4158 BOOST_CHECK(lifted->state == SurfaceSolid::NavigationReliability::OpenSurfaceSet);
4159 BOOST_CHECK_CLOSE(lifted->maxIsolation, delta, 1.e-6);
4160 BOOST_CHECK(lifted->maxIsolationFace >= 0 && lifted->maxIsolationFace != 4);
4162 BOOST_CHECK_CLOSE(lifted->maxIsolationPoint[2], halfZ + delta, 1.e-6);
4169 constexpr double radius = 1.;
4170 constexpr double halfHeight = 2.;
4171 SurfaceSolid openTube(
"rimOpenTube");
4173 openTube.AddCylindricalSurface({0., 0., 0.}, {0., 0., 1.}, {1., 0., 0.}, radius, -halfHeight, halfHeight));
4174 openTube.CloseShape(
false);
4184 BOOST_CHECK_CLOSE(openTube.GetUnmatchedRimLength(), openTube.GetTotalRimLength(), 1.e-9);
4185 BOOST_CHECK_CLOSE(openTube.GetTotalRimLength(), 2. *
surf::kTwoPi * radius, 1.);
4191 BOOST_CHECK_CLOSE(openTube.GetRimChordResolution(), exactSagitta, 1.);
4203 SurfaceSolid unstated(
"rimToleranceUnstated");
4204 addBoxSurfaces(unstated, 1., 1., 1.);
4205 unstated.CloseShape(
false);
4209 SurfaceSolid stated(
"rimToleranceStated");
4210 addBoxSurfaces(stated, 1., 1., 1.);
4211 stated.SetModelTolerance(2.5e-7);
4212 stated.CloseShape(
false);
4213 BOOST_CHECK_CLOSE(stated.GetRimMatchTolerance(), 2.5e-7, 1.e-9);
4228 BOOST_CHECK_CLOSE(circle.
uVariation(0., 1.), 8., 1.e-9);
4231 std::vector<surf::Vec2> poles{{0., 0.}, {1., 1.}, {2., 0.}, {1., -1.}, {0., 0.}};
4232 std::vector<double> knots{0., 0., 0., 0.25, 0.5, 0.75, 1., 1., 1.};
4233 const surf::Curve2D loop = surf::Curve2D::makeBSpline(2, poles, {}, knots);
4235 BOOST_CHECK_GT(loop.
uVariation(0., 1.), 0.5);
4252 const double radius = 1.2;
4253 const double height = 1.5;
4254 const double holeRadius = 0.3;
4256 const double centreV = 0.5 *
height;
4258 std::vector<o2::cad::O2BVHSurfaceSolid::PlanarBoundaryCurve> outer{
4266 const double corner = std::sqrt(2.) / 2.;
4267 const std::array<std::array<double, 2>, 9> unit{{{1., 0.}, {1., 1.}, {0., 1.}, {-1., 1.}, {-1., 0.}, {-1., -1.}, {0., -1.}, {1., -1.}, {1., 0.}}};
4268 std::vector<o2::cad::O2BVHSurfaceSolid::Point2D> poles;
4269 for (
const auto& pole : unit) {
4270 poles.push_back({centreU + holeRadius * pole[0], centreV + holeRadius * pole[1]});
4272 const std::vector<double>
weights{1., corner, 1., corner, 1., corner, 1., corner, 1.};
4273 const std::vector<double> knots{0., 0., 0., 0.25, 0.25, 0.5, 0.5, 0.75, 0.75, 1., 1., 1.};
4274 std::vector<std::vector<o2::cad::O2BVHSurfaceSolid::PlanarBoundaryCurve>> holes{
4278 BOOST_REQUIRE(withHole.
AddCylindricalSurface({0., 0., 0.}, {0., 0., 1.}, {1., 0., 0.}, radius, 0.,
height, 0.,
4280 std::vector<double> contributions;
4282 BOOST_REQUIRE_EQUAL(contributions.size(), 1u);
4287 BOOST_CHECK_CLOSE(contributions[0], radius * radius * chartArea / 3., 1.e-6);
4296 std::vector<surf::Vec2> poles{{0., 0.}, {1., 2.}, {2., 2.}, {3., 0.}};
4297 std::vector<double> knots{0., 1., 2., 3., 4., 5., 6., 7.};
4298 const surf::Curve2D unclamped = surf::Curve2D::makeBSpline(3, poles, {}, knots);
4302 BOOST_CHECK_GT(std::abs(
start.uCoord - poles.front().uCoord) + std::abs(
start.vCoord - poles.front().vCoord), 1.e-3);
4303 BOOST_CHECK_GT(std::abs(
end.uCoord - poles.back().uCoord) + std::abs(
end.vCoord - poles.back().vCoord), 1.e-3);
4306 BOOST_CHECK_SMALL(std::abs(unclamped.
pointAt(1.).
vCoord -
end.vCoord), 1.e-12);
4325std::vector<std::pair<std::vector<unsigned int>, std::vector<unsigned char>>>
4326 boxEdgeIdentity(
double halfX,
double halfY,
double halfZ)
4328 using Key = std::array<long long, 6>;
4329 auto quantize = [](
double value) {
return static_cast<long long>(std::llround(
value * 1.e9)); };
4330 std::map<Key, unsigned int> edgeIds;
4331 std::vector<std::pair<std::vector<unsigned int>, std::vector<unsigned char>>> perFace(6);
4334 const FaceFrame frame = boxFaceFrame(
faceIndex, halfX, halfY, halfZ);
4335 const auto corners = rectangleWire(frame.extentU, frame.extentV);
4337 const auto& startUV = corners[
segment];
4338 const auto& endUV = corners[(
segment + 1) % corners.size()];
4339 auto toGlobal = [&](
const Point2D& uv) {
4340 return Point3D{frame.origin[0] + frame.axisU[0] * uv[0] + frame.axisV[0] * uv[1],
4341 frame.origin[1] + frame.axisU[1] * uv[0] + frame.axisV[1] * uv[1],
4342 frame.origin[2] + frame.axisU[2] * uv[0] + frame.axisV[2] * uv[1]};
4347 quantize(
end[0]), quantize(
end[1]), quantize(
end[2])};
4348 const Key backwardKey{forwardKey[3], forwardKey[4], forwardKey[5],
4349 forwardKey[0], forwardKey[1], forwardKey[2]};
4351 const bool reversed = backwardKey < forwardKey;
4352 const Key canonical = reversed ? backwardKey : forwardKey;
4353 const auto inserted = edgeIds.emplace(canonical,
static_cast<unsigned int>(edgeIds.size()));
4354 perFace[
faceIndex].first.push_back(inserted.first->second);
4356 static_cast<unsigned char>(SurfaceSolid::kEdgeAnchored | (reversed ? SurfaceSolid::kEdgeReversed : 0u)));
4359 BOOST_REQUIRE_EQUAL(edgeIds.size(), 12u);
4364std::unique_ptr<SurfaceSolid> makeIdentifiedBox(
const char*
name,
double halfX,
double halfY,
double halfZ)
4366 auto solid = std::make_unique<SurfaceSolid>(
name);
4367 addBoxSurfaces(*solid, halfX, halfY, halfZ);
4368 const auto identity = boxEdgeIdentity(halfX, halfY, halfZ);
4381 const auto box = makeIdentifiedBox(
"identityBox", 1., 2., 3.);
4382 box->CloseShape(
false);
4398 BOOST_CHECK_SMALL(
box->GetMaxSharedEdgeDeviation(), 1.e-15);
4406 SurfaceSolid openBox(
"identityOpenBox");
4408 BOOST_REQUIRE(addBoxFace(openBox,
faceIndex, 1., 2., 3.));
4410 const auto identity = boxEdgeIdentity(1., 2., 3.);
4414 openBox.CloseShape(
false);
4423 static_cast<int>(SurfaceSolid::NavigationReliability::OpenSurfaceSet));
4431 const auto box = makeIdentifiedBox(
"identityReversed", 1., 2., 3.);
4433 auto identity = boxEdgeIdentity(1., 2., 3.);
4434 for (
auto& flag : identity[0].second) {
4435 flag =
static_cast<unsigned char>(flag ^ SurfaceSolid::kEdgeReversed);
4437 BOOST_REQUIRE(
box->SetSurfaceBoundaryEdges(0, identity[0].first, identity[0].second));
4438 box->CloseShape(
false);
4443 static_cast<int>(SurfaceSolid::NavigationReliability::ReversedFaces));
4447 const auto box = makeIdentifiedBox(
"identityNonManifold", 1., 2., 3.);
4448 BOOST_REQUIRE(addBoxFace(*
box, 0, 1., 2., 3.));
4449 const auto identity = boxEdgeIdentity(1., 2., 3.);
4450 BOOST_REQUIRE(
box->SetSurfaceBoundaryEdges(6, identity[0].first, identity[0].second));
4451 box->CloseShape(
false);
4455 static_cast<int>(SurfaceSolid::NavigationReliability::NonManifold));
4466 constexpr double delta = 3.e-4;
4467 SurfaceSolid shifted(
"identityShifted");
4471 BOOST_REQUIRE(addBoxFace(shifted,
faceIndex, 1., 2., 3.,
false, centre));
4473 const auto identity = boxEdgeIdentity(1., 2., 3.);
4477 shifted.CloseShape(
false);
4483 checkClose(shifted.GetMaxSharedEdgeDeviation(), delta, 1.e-12);
4494 SurfaceSolid flipped(
"identityFlippedWinding");
4495 const auto identity = boxEdgeIdentity(1., 2., 3.);
4496 std::vector<std::vector<unsigned int>>
ids(6);
4497 std::vector<std::vector<unsigned char>>
flags(6);
4499 FaceFrame frame = boxFaceFrame(
faceIndex, 1., 2., 3.);
4500 auto corners = rectangleWire(frame.extentU, frame.extentV);
4503 std::reverse(corners.begin(), corners.end());
4504 BOOST_REQUIRE(flipped.AddPlanarSurface(frame.origin, frame.axisU, frame.axisV, corners));
4505 const size_t n = corners.size();
4508 for (
size_t j = 0;
j <
n; ++
j) {
4515 flipped.CloseShape(
false);
4521 BOOST_CHECK_SMALL(flipped.GetMaxSharedEdgeDeviation(), 1.e-12);
4529 SurfaceSolid partial(
"identityPartial");
4530 addBoxSurfaces(partial, 1., 2., 3.);
4531 const auto identity = boxEdgeIdentity(1., 2., 3.);
4535 partial.CloseShape(
false);
4541 BOOST_CHECK(!partial.SetSurfaceBoundaryEdges(6, identity[0].first, identity[0].second));
4542 BOOST_CHECK(!partial.SetSurfaceBoundaryEdges(0, {1u, 2u}, {0}));
4549 const auto box = makeIdentifiedBox(
"identityPersist", 1., 2., 3.);
4550 box->CloseShape(
false);
4551 BOOST_REQUIRE(
box->HasEdgeIdentity());
4553 const auto restored = writeAndReadBack(*
box);
4554 BOOST_REQUIRE(restored !=
nullptr);
4559 checkClose(restored->GetMaxSharedEdgeDeviation(),
box->GetMaxSharedEdgeDeviation(), 1.e-18);
4560 checkSolidsIdentical(*
box, *restored, 4.5, 5);
4568 constexpr double halfX = 1.;
4569 constexpr double halfY = 2.;
4570 constexpr double halfZ = 3.;
4571 const auto identity = boxEdgeIdentity(halfX, halfY, halfZ);
4573 const auto boxBytes = [&](uint32_t
version,
bool writeIdentity) {
4574 std::vector<char> bytes;
4575 appendSidecarHeader(bytes, 6,
version, 1.e-7, writeIdentity ? 12u : 0u);
4577 appendPlaneRecord(bytes, boxFaceFrame(
faceIndex, halfX, halfY, halfZ));
4581 appendU32(bytes, writeIdentity ?
static_cast<uint32_t
>(
ids.size()) : 0u);
4582 if (writeIdentity) {
4583 for (
size_t e = 0; e <
ids.size(); ++e) {
4584 appendU32(bytes,
ids[e]);
4585 bytes.push_back(
static_cast<char>(
flags[e]));
4593 SurfaceSolid
v3(
"sidecarV3Identity");
4594 const auto v3Path = writeSidecarFile(
"o2_sidecar_v3_identity.bin", boxBytes(3,
true));
4596 std::filesystem::remove(v3Path);
4597 v3.CloseShape(
false);
4602 BOOST_CHECK_SMALL(
v3.GetMaxSharedEdgeDeviation(), 1.e-15);
4606 checkClose(
v3.Capacity(),
reference.Capacity(), 1.e-9);
4610 SurfaceSolid v3Silent(
"sidecarV3Silent");
4611 const auto silentPath = writeSidecarFile(
"o2_sidecar_v3_silent.bin", boxBytes(3,
false));
4613 std::filesystem::remove(silentPath);
4614 v3Silent.CloseShape(
false);
4619 SurfaceSolid
v2(
"sidecarV2StillLoads");
4620 const auto v2Path = writeSidecarFile(
"o2_sidecar_v2_still_loads.bin", boxBytes(2,
false));
4622 std::filesystem::remove(v2Path);
4623 v2.CloseShape(
false);
4626 checkClose(
v2.GetModelTolerance(), 1.e-7, 1.e-18);
4630 std::vector<char> mislabelled;
4631 appendSidecarHeader(mislabelled, 6, 2, 1.e-7);
4634 appendPlaneRecord(mislabelled, boxFaceFrame(
faceIndex, halfX, halfY, halfZ));
4636 SurfaceSolid mislabelledSolid(
"sidecarV3Mislabelled");
4637 const auto badPath = writeSidecarFile(
"o2_sidecar_v3_mislabelled.bin", mislabelled);
4639 std::filesystem::remove(badPath);
4650 const auto box = makeIdentifiedBox(
"identitySampling", 1., 2., 3.);
4651 box->CloseShape(
false);
4652 const double deviation =
box->GetMaxSharedEdgeDeviation();
4653 const bool navigable =
box->IsNavigable();
4659 SurfaceSolid resampled(
"identitySamplingResampled");
4660 const auto identity = boxEdgeIdentity(1., 2., 3.);
4662 const FaceFrame frame = boxFaceFrame(
faceIndex, 1., 2., 3.);
4663 auto corners = rectangleWire(frame.extentU, frame.extentV);
4664 std::vector<Point2D> dense;
4665 std::vector<unsigned int>
ids;
4666 std::vector<unsigned char>
flags;
4669 const auto&
b = corners[(
segment + 1) % corners.size()];
4671 dense.push_back({0.5 * (
a[0] +
b[0]), 0.5 * (
a[1] +
b[1])});
4679 BOOST_REQUIRE(resampled.AddPlanarSurface(frame.origin, frame.axisU, frame.axisV, dense));
4682 resampled.CloseShape(
false);
4693 checkClose(
box->GetMaxSharedEdgeDeviation(), deviation, 1.e-18);
4709std::array<double, 5> torusRayQuartic(
double majorRadius,
double minorRadius,
4712 const double dirDotDir = dir[0] * dir[0] + dir[1] * dir[1] + dir[2] * dir[2];
4713 const double originDotDir =
origin[0] * dir[0] +
origin[1] * dir[1] +
origin[2] * dir[2];
4714 const double originDotOrigin =
4716 const double constantK = majorRadius * majorRadius - minorRadius * minorRadius;
4717 const double transverseE = dir[0] * dir[0] + dir[1] * dir[1];
4718 const double transverseF =
origin[0] * dir[0] +
origin[1] * dir[1];
4720 const double fourRSquared = 4. * majorRadius * majorRadius;
4721 return {dirDotDir * dirDotDir,
4722 4. * dirDotDir * originDotDir,
4723 4. * originDotDir * originDotDir + 2. * dirDotDir * (originDotOrigin + constantK) -
4724 fourRSquared * transverseE,
4725 4. * originDotDir * (originDotOrigin + constantK) - 2. * fourRSquared * transverseF,
4726 (originDotOrigin + constantK) * (originDotOrigin + constantK) -
4727 fourRSquared * transverseG};
4732std::vector<double> torusRootsAtScale(
double scale)
4734 const Point3D origin{2.094269422822338 * scale, 3.292530879918199 * scale,
4735 1.9347519602583996 * scale};
4736 const Point3D dir{-0.7547297076674779, -0.03154700875395883, -0.655276929704412};
4737 const auto c = torusRayQuartic(2.5 * scale, 0.8 * scale,
origin, dir);
4739 return {roots.begin(), roots.end()};
4748 const auto reference = torusRootsAtScale(1.);
4749 BOOST_REQUIRE_EQUAL(
reference.size(), 2u);
4750 for (
const double scale : {0.5, 2., 10.}) {
4751 const auto roots = torusRootsAtScale(scale);
4752 BOOST_REQUIRE_EQUAL(roots.size(),
reference.size());
4753 for (
size_t i = 0;
i < roots.size(); ++
i) {
4754 checkClose(roots[
i],
reference[
i] * scale, 1.e-12);
4771 const auto reference = torusRootsAtScale(1.);
4772 BOOST_REQUIRE_EQUAL(
reference.size(), 2u);
4773 for (
const double factor : {0.15, 0.12, 0.05, 0.01}) {
4774 const auto roots = torusRootsAtScale(factor);
4775 BOOST_REQUIRE_EQUAL(roots.size(),
reference.size());
4776 for (
size_t i = 0;
i < roots.size(); ++
i) {
4784 const double scale = 0.12;
4785 const Point3D origin{2.094269422822338 * scale, 3.292530879918199 * scale,
4786 1.9347519602583996 * scale};
4787 const Point3D dir{-0.7547297076674779, -0.03154700875395883, -0.655276929704412};
4788 const auto c = torusRayQuartic(2.5 * scale, 0.8 * scale,
origin, dir);
4789 const auto evaluate = [&
c](
double t) {
4790 return (((
c[0] * t +
c[1]) * t +
c[2]) * t +
c[3]) * t +
c[4];
4793 const double t = root * scale;
4794 const double span = 0.02 * t;
4795 BOOST_CHECK_LT(evaluate(t - span) * evaluate(t + span), 0.);
4816 const auto dir = std::filesystem::temp_directory_path();
4820 double origin[3] = {1.0, 1.5, 2.0};
4826 auto* cube =
new TGeoBBox(
"cube", 2.0, 2.0, 2.0);
4827 auto* drill =
new TGeoTube(
"drill", 0.0, 0.8, 2.5);
4828 TGeoCompositeShape composite(
"shape",
new TGeoSubtraction(cube, drill,
nullptr,
nullptr));
4830 const std::vector<Point3D> probes{{0.5, 0.5, 0.5}, {1.0, 1.5, 2.0}, {3.0, 1.5, 2.0}, {0.0, 0.0, 0.0}, {1.9, 0.0, 0.0}, {0.0, 0.0, 1.9}, {-1.5, -1.5, 1.0}, {0.79, 0.0, 0.0}, {0.81, 0.0, 0.0}};
4831 const std::vector<Point3D> directions{{1., 0., 0.}, {0., 1., 0.}, {0., 0., 1.}, {-1., 0., 0.}, {0.6, 0.8, 0.}};
4833 for (
const TGeoShape* original : {
static_cast<const TGeoShape*
>(&
box),
4834 static_cast<const TGeoShape*
>(&composite)}) {
4835 const std::string
path = (dir / (std::string(
"o2_shape_sidecar_") + original->ClassName() +
".root")).string();
4837 BOOST_REQUIRE_MESSAGE(harness::saveShapeToRootFile(
path, *original, &error), error);
4839 std::unique_ptr<TGeoShape> loaded(harness::loadShapeFromRootFile(
path, &error));
4840 BOOST_REQUIRE_MESSAGE(loaded !=
nullptr, error);
4841 BOOST_CHECK_EQUAL(std::string(loaded->ClassName()), std::string(original->ClassName()));
4845 BOOST_CHECK_CLOSE(loaded->Capacity(), original->Capacity(), 5.0);
4848 for (
const auto& p : probes) {
4851 original->Safety(p.data(), original->Contains(p.data())));
4852 for (
const auto& d : directions) {
4854 original->DistFromOutside(p.data(), d.data(), 3));
4856 original->DistFromInside(p.data(), d.data(), 3));
4859 std::filesystem::remove(
path);
4866 const auto dir = std::filesystem::temp_directory_path();
4869 BOOST_CHECK(harness::loadShapeFromRootFile((dir /
"o2_shape_absent.root").
string(), &error) ==
nullptr);
4875 const std::string
path = (dir /
"o2_shape_not_a_shape.root").
string();
4877 TFile out(
path.c_str(),
"RECREATE");
4878 TNamed impostor(
"shape",
"not a shape");
4879 out.WriteTObject(&impostor,
"shape");
4885 std::filesystem::remove(
path);
4895 constexpr double kR = 1.5;
4896 constexpr double kZ = 2.0;
4897 constexpr double kError = 0.05;
4898 const TGeoTube truth(
"truth", 0., kR, kZ);
4899 const TGeoTube wrong(
"wrong", 0., kR + kError, kZ);
4905 const auto trueContains = [](
const Point3D& p) {
4906 return (std::hypot(p[0], p[1]) <= kR && std::fabs(p[2]) <= kZ) ? 1 : 0;
4908 const auto trueBoundaryDistance = [](
const Point3D& p) {
4909 const double r = std::hypot(p[0], p[1]);
4910 const double dr = kR -
r;
4911 const double dz = kZ - std::fabs(p[2]);
4912 if (dr > 0. && dz > 0.) {
4913 return std::min(dr, dz);
4915 return std::hypot(std::max(
r - kR, 0.), std::max(std::fabs(p[2]) - kZ, 0.));
4918 std::vector<Point3D> points;
4919 std::vector<int> containsState;
4920 std::vector<double> boundaryDistance;
4921 for (
int ix = -6; ix <= 6; ++ix) {
4922 for (
int iy = -6; iy <= 6; ++iy) {
4923 for (
int iz = -4; iz <= 4; ++iz) {
4924 const Point3D p{0.31 * ix, 0.29 * iy, 0.53 * iz};
4928 if (std::fabs(trueBoundaryDistance(p)) < 1.e-3) {
4931 points.push_back(p);
4932 containsState.push_back(trueContains(p));
4933 boundaryDistance.push_back(trueBoundaryDistance(p));
4937 BOOST_REQUIRE_GT(points.size(), 500u);
4939 harness::ValidationOptions opt;
4940 opt.meshBand = 1.e-6;
4941 opt.distanceTolerance = 1.e-9;
4943 auto containsTruth = harness::validateContainsAgainstOracle(&truth, points, containsState,
4944 boundaryDistance, opt);
4945 auto containsWrong = harness::validateContainsAgainstOracle(&wrong, points, containsState,
4946 boundaryDistance, opt);
4947 BOOST_CHECK_EQUAL(containsTruth.nMismatchUnexplained + containsTruth.nMismatchMissedSurface, 0u);
4948 BOOST_CHECK_GT(containsWrong.nMismatchUnexplained + containsWrong.nMismatchMissedSurface, 0u);
4950 auto safetyTruth = harness::validateSafetyAgainstOracle(&truth, points, boundaryDistance, opt);
4951 auto safetyWrong = harness::validateSafetyAgainstOracle(&wrong, points, boundaryDistance, opt);
4952 BOOST_CHECK_EQUAL(safetyTruth.nMismatchUnexplained + safetyTruth.nMismatchMissedSurface, 0u);
4953 BOOST_CHECK_GT(safetyWrong.nMismatchUnexplained + safetyWrong.nMismatchMissedSurface, 0u);
4958 std::vector<harness::Ray> rays;
4959 std::vector<double> rayDistance;
4960 std::vector<int> originState;
4961 for (
int i = 0;
i < 400; ++
i) {
4962 const double phi = 0.0173 *
i;
4963 const double z = -1.9 + 0.0095 *
i;
4964 const Point3D target{0.9 * kR * std::cos(2.1 * phi), 0.9 * kR * std::sin(2.1 * phi),
z};
4965 const Point3D origin{5.0 * std::cos(phi), 5.0 * std::sin(phi), 3.0 - 0.01 *
i};
4967 const double norm = std::sqrt(dir[0] * dir[0] + dir[1] * dir[1] + dir[2] * dir[2]);
4968 for (
auto& component : dir) {
4971 rays.push_back(harness::Ray{
origin, dir});
4972 rayDistance.push_back(truth.DistFromOutside(
origin.data(), dir.data(), 3));
4973 originState.push_back(trueContains(
origin));
4976 auto distTruth = harness::validateDistanceAgainstOracle(&truth, rays, rayDistance,
4977 false, opt, originState);
4978 auto distWrong = harness::validateDistanceAgainstOracle(&wrong, rays, rayDistance,
4979 false, opt, originState);
4980 BOOST_CHECK_EQUAL(distTruth.nMismatchUnexplained + distTruth.nMismatchMissedSurface, 0u);
4981 BOOST_CHECK_GT(distWrong.nMismatchUnexplained + distWrong.nMismatchMissedSurface, 0u);
5002TGeoCompositeShape* makePlacedTube(
const char*
name,
double rmin,
double rmax,
double dz,
5003 TGeoMatrix* matrixA, TGeoMatrix* matrixB)
5005 auto*
left =
new TGeoTube(Form(
"%s_l",
name), rmin, rmax, dz);
5006 auto*
right =
new TGeoTube(Form(
"%s_r",
name), rmin, rmax, dz);
5008 return new TGeoCompositeShape(
name,
node);
5017 constexpr double kRmin = 0.4;
5018 constexpr double kRmax = 1.0;
5019 constexpr double kDz = 5.0;
5020 const TGeoTube
reference(
"reference", kRmin, kRmax, kDz);
5022 auto* rotation =
new TGeoRotation(
"csgRot", 0., 0., 0.);
5023 rotation->RotateX(30.);
5024 rotation->RotateZ(17.);
5025 auto* matrixA =
new TGeoCombiTrans(0.3, 5.916, 2.0, rotation);
5026 auto* matrixB =
new TGeoCombiTrans(0.3, 5.916, 2.0, rotation);
5027 const TGeoCombiTrans placement(0.3, 5.916, 2.0, rotation);
5028 std::unique_ptr<TGeoCompositeShape> placed(
5029 makePlacedTube(
"csgPlaced", kRmin, kRmax, kDz, matrixA, matrixB));
5031 std::size_t probes = 0;
5032 std::size_t inside = 0;
5033 std::size_t outside = 0;
5034 for (
int ix = -8; ix <= 8; ++ix) {
5035 for (
int iy = -8; iy <= 8; ++iy) {
5036 for (
int iz = -8; iz <= 8; ++iz) {
5037 const Point3D master{0.3 + 0.37 * ix, 5.916 + 0.41 * iy, 2.0 + 0.43 * iz};
5039 placement.MasterToLocal(master.data(), local.data());
5042 const double r = std::hypot(local[0], local[1]);
5043 if (std::fabs(
r - kRmin) < 1.e-9 || std::fabs(
r - kRmax) < 1.e-9 ||
5044 std::fabs(std::fabs(local[2]) - kDz) < 1.e-9) {
5048 const bool wanted =
reference.Contains(local.data());
5049 BOOST_REQUIRE_EQUAL(placed->Contains(master.data()), wanted);
5050 BOOST_REQUIRE_CLOSE_FRACTION(placed->Safety(master.data(), wanted),
5051 reference.Safety(local.data(), wanted), 1.e-12);
5052 wanted ? ++inside : ++outside;
5055 Point3D{0.5773502691896258, 0.5773502691896258, 0.5773502691896258}}) {
5057 placement.MasterToLocalVect(dir.data(), localDir.data());
5059 BOOST_REQUIRE_CLOSE_FRACTION(placed->DistFromInside(master.data(), dir.data(), 3),
5060 reference.DistFromInside(local.data(), localDir.data(), 3),
5063 const double got = placed->DistFromOutside(master.data(), dir.data(), 3);
5064 const double want =
reference.DistFromOutside(local.data(), localDir.data(), 3);
5066 BOOST_REQUIRE_GT(got, 1.e20);
5068 BOOST_REQUIRE_CLOSE_FRACTION(got, want, 1.e-12);
5076 BOOST_CHECK_GT(probes, 2000u);
5077 BOOST_CHECK_GT(inside, 100u);
5078 BOOST_CHECK_GT(outside, 100u);
5084 const TGeoTube wrong(
"wrongReference", kRmin, kRmax + 0.05, kDz);
5085 std::size_t disagreements = 0;
5086 std::size_t shellProbes = 0;
5087 for (
int iphi = 0; iphi < 24; ++iphi) {
5088 const double phi = 2. * M_PI * iphi / 24.;
5089 for (
int iz = -3; iz <= 3; ++iz) {
5090 const double radius = kRmax + 0.025;
5091 const Point3D local{radius * std::cos(phi), radius * std::sin(phi), 1.3 * iz};
5093 placement.LocalToMaster(local.data(), master.data());
5095 disagreements += (placed->Contains(master.data()) != wrong.Contains(local.data())) ? 1 : 0;
5109 constexpr double kEyeRmin = 0.7;
5110 constexpr double kEyeRmax = 1.2;
5111 constexpr double kEyeDz = 0.75;
5112 constexpr double kRodR = 0.6;
5113 constexpr double kRodDz = 3.5;
5114 constexpr double kRodCentre = 3.5;
5116 auto* eyeRotation =
new TGeoRotation(
"csgEyeRot", 90., 90., 0.);
5117 auto* eyeMatrix =
new TGeoCombiTrans(0., 0., 0., eyeRotation);
5118 auto* rodMatrix =
new TGeoTranslation(0., 0., kRodCentre);
5119 auto* eye =
new TGeoTube(
"csgEye", kEyeRmin, kEyeRmax, kEyeDz);
5120 auto* rod =
new TGeoTube(
"csgRod", 0., kRodR, kRodDz);
5121 auto*
node =
new TGeoUnion(eye, rod, eyeMatrix, rodMatrix);
5122 std::unique_ptr<TGeoCompositeShape> ram(
new TGeoCompositeShape(
"csgRam",
node));
5124 const auto closedForm = [&](
const Point3D& p) {
5125 const double rEye = std::hypot(p[1], p[2]);
5126 const bool inEye = rEye >= kEyeRmin && rEye <= kEyeRmax && std::fabs(p[0]) <= kEyeDz;
5127 const double rRod = std::hypot(p[0], p[1]);
5128 const bool inRod = rRod <= kRodR && p[2] >= 0. && p[2] <= 2. * kRodDz;
5129 return inEye || inRod;
5131 const auto nearWall = [&](
const Point3D& p) {
5132 const double rEye = std::hypot(p[1], p[2]);
5133 const double rRod = std::hypot(p[0], p[1]);
5134 return std::fabs(rEye - kEyeRmin) < 1.e-9 || std::fabs(rEye - kEyeRmax) < 1.e-9 ||
5135 std::fabs(std::fabs(p[0]) - kEyeDz) < 1.e-9 || std::fabs(rRod - kRodR) < 1.e-9 ||
5136 std::fabs(p[2]) < 1.e-9 || std::fabs(p[2] - 2. * kRodDz) < 1.e-9;
5139 const std::filesystem::path
path =
5140 std::filesystem::temp_directory_path() /
"o2_csg_ram_shape.root";
5143 BOOST_REQUIRE_MESSAGE(harness::saveShapeToRootFile(
path.string(), *ram, &error), error);
5144 std::unique_ptr<TGeoShape> loaded(harness::loadShapeFromRootFile(
path.string(), &error));
5145 BOOST_REQUIRE_MESSAGE(loaded !=
nullptr, error);
5146 BOOST_CHECK_EQUAL(std::string(loaded->ClassName()), std::string(
"TGeoCompositeShape"));
5148 std::size_t inside = 0;
5149 std::size_t outside = 0;
5150 for (
int ix = -6; ix <= 6; ++ix) {
5151 for (
int iy = -6; iy <= 6; ++iy) {
5152 for (
int iz = -4; iz <= 20; ++iz) {
5153 const Point3D p{0.23 * ix, 0.27 * iy, 0.41 * iz};
5157 const bool wanted = closedForm(p);
5158 BOOST_REQUIRE_EQUAL(ram->Contains(p.data()), wanted);
5159 BOOST_REQUIRE_EQUAL(loaded->Contains(p.data()), wanted);
5160 wanted ? ++inside : ++outside;
5164 BOOST_CHECK_GT(inside, 50u);
5165 BOOST_CHECK_GT(outside, 500u);
5174 double minCapacity = ram->Capacity();
5175 double maxCapacity = minCapacity;
5176 for (
int i = 0;
i < 5; ++
i) {
5177 const double sampled = ram->Capacity();
5178 minCapacity = std::min(minCapacity, sampled);
5179 maxCapacity = std::max(maxCapacity, sampled);
5181 BOOST_CHECK_GT(minCapacity, 0.);
5182 const double spread = (maxCapacity - minCapacity) / (0.5 * (maxCapacity + minCapacity));
5183 BOOST_CHECK_GT(spread, 1.e-4);
5185 std::filesystem::remove(
path);
5220 BOOST_REQUIRE_EQUAL(boxCrossings.size(), 2u);
5221 BOOST_CHECK_SMALL(boxCrossings[0].t - 4., 1.e-12);
5222 BOOST_CHECK_SMALL(boxCrossings[1].t - 6., 1.e-12);
5226 TGeoTube tube(
"xrayTube", 0.5, 1.0, 2.0);
5228 BOOST_REQUIRE_EQUAL(tubeCrossings.size(), 4u);
5229 const double expected[4] = {4.0, 4.5, 5.5, 6.0};
5230 const int senses[4] = {+1, -1, +1, -1};
5231 for (
int i = 0;
i < 4; ++
i) {
5232 BOOST_CHECK_SMALL(tubeCrossings[
i].t -
expected[
i], 1.e-12);
5250 std::array<std::pair<std::unique_ptr<SurfaceSolid>,
double>, 7> fixtures{{
5251 {makeBoxSolid(
"xrayLoopBox", 1., 2., 3.), 4.},
5252 {makeTubeSolid(
"xrayLoopTube", 0., 2., 3.), 4.},
5253 {makeTubeSolid(
"xrayLoopHollowTube", 1., 2., 3.), 4.},
5254 {makeConeSolid(
"xrayLoopCone", 2., 1., 3.), 4.},
5255 {makeSphereSolid(
"xrayLoopSphere", 2.5), 3.5},
5256 {makeTorusSolid(
"xrayLoopTorus", 3., 1.), 4.5},
5257 {makeCapsuleSolid(
"xrayLoopCapsule", 1., 1.5), 3.},
5259 size_t comparedRays = 0;
5260 size_t comparedCrossings = 0;
5261 for (
const auto& [solid, extent] : fixtures) {
5262 BOOST_TEST_CONTEXT(
"fixture = " << solid->GetName())
5264 BOOST_REQUIRE(solid->HasBVH());
5265 const XRayPoint lo{-extent, -extent, -extent};
5266 const XRayPoint hi{extent, extent, extent};
5270 for (
const auto& ray : raster.
rays) {
5274 stepWithShapeApi(solid.get(), ray.origin, ray.dir, ray.tMax, cfg, bvhStats);
5277 ray.origin, ray.dir, ray.tMax, cfg, loopStats,
5278 [&solid](
const double* p) { return solid->Contains_Loop(p); },
5279 [&solid](
const double* p,
const double* d) { return solid->DistFromOutside_Loop(p, d); },
5280 [&solid](
const double* p,
const double* d) { return solid->DistFromInside_Loop(p, d); });
5281 BOOST_REQUIRE_EQUAL(viaBVH.size(), viaLoop.size());
5282 for (
size_t i = 0;
i < viaBVH.size(); ++
i) {
5288 comparedCrossings += viaBVH.size();
5292 BOOST_CHECK_GT(comparedRays, 2000u);
5293 BOOST_CHECK_GT(comparedCrossings, 2000u);
5301 const std::vector<Crossing> truth{{4.0, +1}, {4.5, -1}, {5.5, +1}, {6.0, -1}};
5302 const double tolerance = 1.e-6;
5312 auto perturbed = truth;
5313 perturbed[2].t += 1.e-3;
5315 compareLists(perturbed, truth, {}, {}, tolerance, displaced);
5320 BOOST_CHECK_SMALL(displaced.
worstDeltaT - 1.e-3, 1.e-12);
5322 auto dropped = truth;
5323 dropped.erase(dropped.begin() + 1);
5329 auto doubled = truth;
5330 doubled.insert(doubled.begin() + 1, {4.2, -1});
5332 compareLists(doubled, truth, {}, {}, tolerance, spurious);
5336 auto flipped = truth;
5337 flipped[1].kind = +1;
5339 compareLists(flipped, truth, {}, {}, tolerance, sense);
5344 auto nudged = truth;
5345 nudged[0].t += 1.e-9;
5374 auditCrossingList({{1.0, +1}, {2.0, -1}, {4.0, +1}, {6.0, -1}}, &
box,
origin, dir, 10., cfg,
5387 for (
const int n : {5, 16, 41}) {
5390 std::vector<double> byBeam(raster.
beams.size(), 0.);
5391 for (
const auto& ray : raster.
rays) {
5392 const double before = stats.insideLength;
5394 auditCrossingList(crossings,
nullptr, ray.origin, ray.dir, ray.tMax, cfg, stats);
5395 byBeam[ray.beam] += stats.insideLength - before;
5397 BOOST_CHECK_SMALL(
chordVolume(raster, byBeam) - 24., 1.e-9);
5404 TGeoSphere sphere(
"xrayVolSphere", 0., 1.);
5405 const double exact = 4. / 3. * 3.14159265358979323846;
5406 for (
const int n : {24, 96}) {
5409 for (
const auto& ray : raster.
rays) {
5410 const auto crossings =
stepWithShapeApi(&sphere, ray.origin, ray.dir, ray.tMax, cfg, stats);
5411 auditCrossingList(crossings,
nullptr, ray.origin, ray.dir, ray.tMax, cfg, stats);
5413 const double volume = stats.insideLength * raster.
cellArea[0];
5414 BOOST_CHECK_LT(std::fabs(volume - exact) / exact, 2.e-3);
5427 for (
const auto& ray : raster.
rays) {
5428 BOOST_REQUIRE(!
box.Contains(ray.origin.data()));
5430 const double end[3] = {ray.origin[0] + ray.tMax * ray.dir[0],
5431 ray.origin[1] + ray.tMax * ray.dir[1],
5432 ray.origin[2] + ray.tMax * ray.dir[2]};
5433 BOOST_REQUIRE(!
box.Contains(
end));
5440 for (
const auto& beams : {axes, fan}) {
5441 for (
const auto& beam : beams) {
5442 BOOST_CHECK_SMALL(
dot3(beam.dir, beam.dir) - 1., 1.e-12);
5443 BOOST_CHECK_SMALL(
dot3(beam.u, beam.v), 1.e-12);
5444 BOOST_CHECK_SMALL(
dot3(beam.u, beam.dir), 1.e-12);
5445 BOOST_CHECK_SMALL(
dot3(beam.v, beam.dir), 1.e-12);
5450 double worstFanAlignment = -1.;
5451 for (
size_t i = 0;
i < fan.size(); ++
i) {
5452 for (
size_t j =
i + 1;
j < fan.size(); ++
j) {
5453 worstFanAlignment = std::max(worstFanAlignment, std::fabs(
dot3(fan[
i].dir, fan[
j].dir)));
5456 BOOST_CHECK_LT(worstFanAlignment, 0.999);
5491double quarticBackwardError(
const std::array<double, 5>& coefficients,
double x)
5493 double value = 0., magnitude = 0., power = 1.;
5494 for (
int i = 0;
i < 5; ++
i) {
5495 const double term = coefficients[4 -
i] * power;
5497 magnitude += std::abs(term);
5500 return magnitude > 0. ? std::abs(
value) / magnitude : std::abs(
value);
5504std::array<double, 5> quarticFromRoots(
double r1,
double r2,
double r3,
double r4)
5506 return {1., -(r1 + r2 + r3 + r4),
5507 r1 * r2 + r1 * r3 + r1 * r4 + r2 * r3 + r2 * r4 + r3 * r4,
5508 -(r1 * r2 * r3 + r1 * r2 * r4 + r1 * r3 * r4 + r2 * r3 * r4), r1 * r2 * r3 * r4};
5511std::vector<double> sortedRoots(
const std::array<double, 5>&
c,
surf::QuarticBranch* branch =
nullptr)
5514 std::vector<double> roots(found.begin(), found.end());
5515 std::sort(roots.begin(), roots.end());
5520void checkRootsAreRoots(
const std::array<double, 5>&
c,
const std::vector<double>& roots)
5522 for (
const double root : roots) {
5523 BOOST_CHECK_LT(quarticBackwardError(
c, root), 1.e-12);
5539 const std::array<double, 5>
c{1.0, -1501.7280000044018, 845808.25396968238, -211752288.545858,
5540 19882619385.616932};
5545 const double firstRoot = 375.3392295779947145;
5546 const double secondRoot = 375.5247704240909448;
5553 const auto roots = sortedRoots(
c, &branch);
5554 BOOST_REQUIRE_EQUAL(roots.size(), 2u);
5555 checkClose(roots[0], firstRoot, 1.e-6);
5556 checkClose(roots[1], secondRoot, 1.e-6);
5557 checkRootsAreRoots(
c, roots);
5560 BOOST_CHECK(branch == surf::QuarticBranch::Biquadratic);
5572 for (
const double k : {1.e6, 1.e3, 1., 1.e-1, 1.e-2, 1.e-3, 1.e-4, 1.e-5, 1.e-6, 1.e-8}) {
5573 const auto c = quarticFromRoots(1. * k, 2. * k, 3. * k, 7. * k);
5575 const auto roots = sortedRoots(
c, &branch);
5576 BOOST_REQUIRE_EQUAL(roots.size(), 4u);
5577 const double expected[4] = {1. * k, 2. * k, 3. * k, 7. * k};
5578 for (
int i = 0;
i < 4; ++
i) {
5581 checkRootsAreRoots(
c, roots);
5584 BOOST_CHECK(branch == surf::QuarticBranch::Resolvent);
5591 for (
const double k : {1., 1.e-2, 1.e-4, 1.e-6}) {
5592 const auto c = quarticFromRoots(-2. * k, -1. * k, 1. * k, 2.1 * k);
5593 const auto roots = sortedRoots(
c);
5594 BOOST_REQUIRE_EQUAL(roots.size(), 4u);
5595 const double expected[4] = {-2. * k, -1. * k, 1. * k, 2.1 * k};
5596 for (
int i = 0;
i < 4; ++
i) {
5599 checkRootsAreRoots(
c, roots);
5611 const std::array<double, 5>
c{1., 0., -5., 0., 4.};
5613 const auto roots = sortedRoots(
c, &branch);
5614 BOOST_CHECK(branch == surf::QuarticBranch::Biquadratic);
5615 BOOST_REQUIRE_EQUAL(roots.size(), 4u);
5616 const double expected[4] = {-2., -1., 1., 2.};
5617 for (
int i = 0;
i < 4; ++
i) {
5631 for (
const double centre : {0., 1., 100.}) {
5632 for (
const double k : {1., 1.e-3, 1.e3}) {
5633 const auto c = quarticFromRoots((centre - 2.) * k, (centre - 1.) * k, (centre + 1.) * k,
5636 const auto roots = sortedRoots(
c, &branch);
5637 BOOST_CHECK(branch == surf::QuarticBranch::Biquadratic);
5638 BOOST_REQUIRE_EQUAL(roots.size(), 4u);
5639 const double expected[4] = {(centre - 2.) * k, (centre - 1.) * k, (centre + 1.) * k,
5641 for (
int i = 0;
i < 4; ++
i) {
5642 BOOST_CHECK_LT(std::abs(roots[
i] -
expected[
i]), 1.e-9 * std::max(1.e-30, std::abs(
expected[
i])));
5644 checkRootsAreRoots(
c, roots);
5659 BOOST_CHECK(branch == surf::QuarticBranch::NotAQuartic);
5661 for (
const double k : {1.e4, 1., 1.e-4, 1.e-8}) {
5663 const double k2 = k * k;
5664 const std::array<double, 5> biquadratic{1., 0., 5. * k2, 0., 4. * k2 * k2};
5667 BOOST_CHECK(branch == surf::QuarticBranch::Biquadratic);
5671 const std::array<double, 5> asymmetric{1., 2. * k, 6. * k2, 2. * k2 * k, 5. * k2 * k2};
5673 BOOST_CHECK(branch == surf::QuarticBranch::Resolvent);
5686 for (
const double k : {1., 1.e-2, 1.e-4, 1.e-6}) {
5687 for (
const double delta : {1.e-1, 1.e-3, 1.e-6, 1.e-9, 1.e-12, 0.}) {
5688 const auto c = quarticFromRoots(-2. * k, -1. * k, 1. * k, (2. + delta) * k);
5689 const auto roots = sortedRoots(
c);
5690 BOOST_REQUIRE_EQUAL(roots.size(), 4u);
5691 const double expected[4] = {-2. * k, -1. * k, 1. * k, (2. + delta) * k};
5692 for (
int i = 0;
i < 4; ++
i) {
5695 checkRootsAreRoots(
c, roots);
5717 const Point3D centre{0.3, -0.7, 1.1};
5718 const Point3D axis{0., 0., 1.};
5719 const Point3D refU{1., 0., 0.};
5720 constexpr double radius = 2.5;
5721 constexpr double phi = 0.9;
5724 const double cx = centre[0] + radius * std::cos(phi);
5725 const double cy = centre[1] + radius * std::sin(phi);
5726 const Double_t onCylinder[3] = {cx, cy, centre[2] + 0.4};
5727 const Double_t awayFromAxis[3] = {std::cos(phi), std::sin(phi), 0.};
5729 for (
const bool innerWall : {
false,
true}) {
5730 SurfaceSolid solid(innerWall ?
"streamK_innerCyl" :
"streamK_outerCyl");
5731 BOOST_REQUIRE(solid.AddCylindricalSurface(centre, axis, refU, radius, -1., 1., 0.,
surf::kTwoPi, innerWall));
5732 Double_t
n[3] = {0., 0., 0.};
5733 solid.ComputeNormal(onCylinder,
nullptr,
n);
5734 const double alignment =
n[0] * awayFromAxis[0] +
n[1] * awayFromAxis[1] +
n[2] * awayFromAxis[2];
5736 BOOST_CHECK_CLOSE(alignment, innerWall ? -1. : 1., 1.e-9);
5743 for (
const bool innerWall : {
false,
true}) {
5744 SurfaceSolid solid(innerWall ?
"streamK_innerCone" :
"streamK_outerCone");
5746 BOOST_REQUIRE(solid.AddConicalSurface(centre, axis, refU, 1., 3., 0., 2., 0.,
surf::kTwoPi, innerWall));
5747 const double h = 1.0;
5748 const double r = 2.0;
5749 const Double_t onCone[3] = {centre[0] +
r * std::cos(phi), centre[1] +
r * std::sin(phi), centre[2] +
h};
5750 Double_t
n[3] = {0., 0., 0.};
5751 solid.ComputeNormal(onCone,
nullptr,
n);
5755 const double radial =
n[0] * std::cos(phi) +
n[1] * std::sin(phi);
5756 BOOST_CHECK_GT(innerWall ? -radial : radial, 0.5);
5759 for (
const bool innerWall : {
false,
true}) {
5760 SurfaceSolid solid(innerWall ?
"streamK_innerSph" :
"streamK_outerSph");
5761 BOOST_REQUIRE(solid.AddSphericalSurface(centre, axis, refU, radius, 0.,
surf::kPi, 0.,
surf::kTwoPi, innerWall));
5762 const double theta = 1.1;
5763 const Double_t onSphere[3] = {centre[0] + radius * std::sin(theta) * std::cos(phi),
5764 centre[1] + radius * std::sin(theta) * std::sin(phi),
5765 centre[2] + radius * std::cos(theta)};
5766 const double outward[3] = {std::sin(theta) * std::cos(phi), std::sin(theta) * std::sin(phi), std::cos(theta)};
5767 Double_t
n[3] = {0., 0., 0.};
5768 solid.ComputeNormal(onSphere,
nullptr,
n);
5769 const double alignment =
n[0] * outward[0] +
n[1] * outward[1] +
n[2] * outward[2];
5770 BOOST_CHECK_CLOSE(alignment, innerWall ? -1. : 1., 1.e-9);
5796TGeoCombiTrans* makeStreamNPlacement()
5798 auto* rotation =
new TGeoRotation(
"streamNRot", 0., 0., 0.);
5799 rotation->RotateX(30.);
5800 rotation->RotateZ(17.);
5801 rotation->RotateY(-41.);
5802 return new TGeoCombiTrans(0.3, 5.916, 2.0, rotation);
5809 const auto dir = std::filesystem::temp_directory_path();
5810 const std::string
path = (dir /
"o2_shape_placed.root").
string();
5812 const TGeoTube tube(
"shape", 0.4, 1.0, 5.0);
5813 std::unique_ptr<TGeoCombiTrans> placement(makeStreamNPlacement());
5816 BOOST_REQUIRE_MESSAGE(harness::saveShapeToRootFile(
path, tube, placement.get(), &error), error);
5818 std::unique_ptr<TGeoShape> loaded(harness::loadShapeFromRootFile(
path, &error));
5819 BOOST_REQUIRE_MESSAGE(loaded !=
nullptr, error);
5822 std::unique_ptr<TGeoHMatrix> back(harness::loadShapePlacementFromRootFile(
path));
5823 BOOST_REQUIRE(back !=
nullptr);
5824 for (
int i = 0;
i < 9; ++
i) {
5827 for (
int i = 0;
i < 3; ++
i) {
5832 const Point3D master{0.9, 6.2, 3.1};
5835 back->MasterToLocal(master.data(), viaFile.data());
5836 placement->MasterToLocal(master.data(), viaOriginal.data());
5837 for (
int i = 0;
i < 3; ++
i) {
5840 std::filesystem::remove(
path);
5846 const auto dir = std::filesystem::temp_directory_path();
5847 const TGeoTube tube(
"shape", 0.4, 1.0, 5.0);
5852 const std::string legacy = (dir /
"o2_shape_legacy.root").
string();
5853 BOOST_REQUIRE_MESSAGE(harness::saveShapeToRootFile(legacy, tube, &error), error);
5854 BOOST_CHECK(harness::loadShapePlacementFromRootFile(legacy) ==
nullptr);
5858 const std::string identity = (dir /
"o2_shape_identity.root").
string();
5859 TGeoHMatrix unit(
"unit");
5860 BOOST_REQUIRE_MESSAGE(harness::saveShapeToRootFile(identity, tube, &unit, &error), error);
5861 BOOST_CHECK(harness::loadShapePlacementFromRootFile(identity) ==
nullptr);
5865 BOOST_CHECK(harness::loadShapePlacementFromRootFile((dir /
"o2_shape_nothing.root").
string()) ==
5868 std::filesystem::remove(legacy);
5869 std::filesystem::remove(identity);
5876 constexpr double kRmin = 0.4;
5877 constexpr double kRmax = 1.0;
5878 constexpr double kDz = 5.0;
5880 std::unique_ptr<TGeoCombiTrans> placement(makeStreamNPlacement());
5881 const TGeoTube placedPrimitive(
"streamNTube", kRmin, kRmax, kDz);
5883 std::unique_ptr<TGeoCompositeShape> composite(
5884 makePlacedTube(
"streamNComposite", kRmin, kRmax, kDz,
new TGeoCombiTrans(*placement),
5885 new TGeoCombiTrans(*placement)));
5887 std::size_t probes = 0;
5888 std::size_t inside = 0;
5889 std::size_t disagreements = 0;
5892 const TGeoTube wrong(
"streamNWrong", kRmin, kRmax * 1.05, kDz);
5893 std::size_t controlDisagreements = 0;
5895 for (
int ix = -8; ix <= 8; ++ix) {
5896 for (
int iy = -8; iy <= 8; ++iy) {
5897 for (
int iz = -8; iz <= 8; ++iz) {
5898 const Point3D master{0.3 + 0.37 * ix, 5.916 + 0.41 * iy, 2.0 + 0.43 * iz};
5900 placement->MasterToLocal(master.data(), local.data());
5901 const double r = std::hypot(local[0], local[1]);
5902 if (std::fabs(
r - kRmin) < 1.e-9 || std::fabs(
r - kRmax) < 1.e-9 ||
5903 std::fabs(
r - kRmax * 1.05) < 1.e-9 || std::fabs(std::fabs(local[2]) - kDz) < 1.e-9) {
5907 const bool wanted = composite->Contains(master.data());
5908 if (placedPrimitive.Contains(local.data()) != wanted) {
5911 if (wrong.Contains(local.data()) != wanted) {
5912 ++controlDisagreements;
5917 BOOST_REQUIRE_CLOSE_FRACTION(placedPrimitive.Safety(local.data(), wanted),
5918 composite->Safety(master.data(), wanted), 1.e-12);
5920 Point3D{0.5773502691896258, 0.5773502691896258,
5921 0.5773502691896258}}) {
5923 placement->MasterToLocalVect(dir.data(), localDir.data());
5925 BOOST_REQUIRE_CLOSE_FRACTION(placedPrimitive.DistFromInside(local.data(),
5926 localDir.data(), 3),
5927 composite->DistFromInside(master.data(), dir.data(), 3),
5930 const double got = placedPrimitive.DistFromOutside(local.data(), localDir.data(), 3);
5931 const double want = composite->DistFromOutside(master.data(), dir.data(), 3);
5933 BOOST_REQUIRE_GT(got, 1.e20);
5935 BOOST_REQUIRE_CLOSE_FRACTION(got, want, 1.e-12);
5943 BOOST_CHECK_GT(controlDisagreements, 0u);
5944 BOOST_CHECK_GT(inside, 100u);
5945 BOOST_CHECK_GT(probes, 3000u);
5951 constexpr double kRmin = 0.4;
5952 constexpr double kRmax = 1.0;
5953 constexpr double kDz = 5.0;
5954 const double analytic = TMath::Pi() * (kRmax * kRmax - kRmin * kRmin) * 2. * kDz;
5956 const TGeoTube tube(
"streamNCapTube", kRmin, kRmax, kDz);
5957 BOOST_CHECK_CLOSE_FRACTION(tube.Capacity(), analytic, 1.e-14);
5961 std::unique_ptr<TGeoCombiTrans> placement(makeStreamNPlacement());
5962 std::unique_ptr<TGeoCompositeShape> composite(
5963 makePlacedTube(
"streamNCapComposite", kRmin, kRmax, kDz,
new TGeoCombiTrans(*placement),
5964 new TGeoCombiTrans(*placement)));
5969 const double first = composite->Capacity();
5970 const double second = composite->Capacity();
5971 BOOST_CHECK_NE(
first, second);
5972 BOOST_CHECK_GT(std::fabs(
first - analytic) / analytic, 1.e-6);
5983 constexpr double kRmin = 0.4;
5984 constexpr double kRmax = 1.0;
5985 constexpr double kDz = 5.0;
5987 std::unique_ptr<TGeoCombiTrans> shapePlacementOwned(makeStreamNPlacement());
5988 const TGeoHMatrix shapePlacement(*shapePlacementOwned);
5989 auto* partRotation =
new TGeoRotation(
"streamNPartRot", 37., 24., 61.);
5990 const TGeoCombiTrans partPlacement(-2.0, 7.0, 1.5, partRotation);
5995 partPlacement.MasterToLocal(master.data(), partFrame.data());
5996 shapePlacement.MasterToLocal(partFrame.data(), shapeFrame.data());
5997 const double r = std::hypot(shapeFrame[0], shapeFrame[1]);
5998 onWall = std::fabs(
r - kRmin) < 1.e-9 || std::fabs(
r - kRmax) < 1.e-9 ||
5999 std::fabs(std::fabs(shapeFrame[2]) - kDz) < 1.e-9;
6000 return r >= kRmin &&
r <= kRmax && std::fabs(shapeFrame[2]) <= kDz;
6005 TGeoHMatrix correct(partPlacement);
6006 correct.Multiply(&shapePlacement);
6007 TGeoHMatrix reversed(shapePlacement);
6008 reversed.Multiply(&partPlacement);
6009 TGeoHMatrix transposedRotation(shapePlacement);
6012 const double*
r = shapePlacement.GetRotationMatrix();
6013 for (
int i = 0;
i < 3; ++
i) {
6014 for (
int j = 0;
j < 3; ++
j) {
6015 rt[3 *
i +
j] =
r[3 *
j +
i];
6018 transposedRotation.SetRotation(rt);
6019 transposedRotation.SetBit(TGeoMatrix::kGeoRotation);
6021 TGeoHMatrix withTransposed(partPlacement);
6022 withTransposed.Multiply(&transposedRotation);
6023 const TGeoHMatrix partOnly(partPlacement);
6025 const std::vector<std::pair<std::string, const TGeoHMatrix*>> candidates{
6026 {
"part*shape", &correct},
6027 {
"shape*part", &reversed},
6028 {
"part*shape^T", &withTransposed},
6029 {
"part only", &partOnly}};
6035 const double* centre = correct.GetTranslation();
6037 std::vector<size_t> disagreements(candidates.size(), 0);
6039 size_t insideProbes = 0;
6041 for (
size_t c = 0;
c < candidates.size(); ++
c) {
6044 auto* manager =
new TGeoManager((
"streamN_" +
std::to_string(
c)).c_str(),
"composition order");
6045 auto* material =
new TGeoMaterial(
"Vacuum", 0., 0., 0.);
6046 auto* medium =
new TGeoMedium(
"Vacuum", 1, material);
6047 auto* world =
new TGeoVolume(
"TOP",
new TGeoBBox(
"streamNWorld", 30., 30., 30.), medium);
6048 auto* part =
new TGeoVolume(
"PART",
new TGeoTube(
"streamNNodeTube", kRmin, kRmax, kDz), medium);
6049 world->AddNode(part, 1,
new TGeoHMatrix(*candidates[
c].second));
6050 manager->SetTopVolume(world);
6051 manager->CloseGeometry();
6053 size_t localProbes = 0;
6054 size_t localInside = 0;
6055 for (
int ix = -14; ix <= 14; ++ix) {
6056 for (
int iy = -14; iy <= 14; ++iy) {
6057 for (
int iz = -14; iz <= 14; ++iz) {
6058 const Point3D master{centre[0] + 0.45 * ix, centre[1] + 0.47 * iy,
6059 centre[2] + 0.43 * iz};
6060 bool onWall =
false;
6061 const bool wanted =
reference(master, onWall);
6069 TGeoNode*
node = manager->FindNode(master[0], master[1], master[2]);
6070 const bool got =
node !=
nullptr && std::string(
node->GetVolume()->GetName()) ==
"PART";
6071 if (got != wanted) {
6077 probes = localProbes;
6078 insideProbes = localInside;
6080 gGeoManager =
nullptr;
6084 BOOST_CHECK_GT(probes, 5000u);
6085 BOOST_CHECK_GT(insideProbes, 200u);
6087 BOOST_CHECK_GT(disagreements[1], 0u);
6088 BOOST_CHECK_GT(disagreements[2], 0u);
6089 BOOST_CHECK_GT(disagreements[3], 0u);
6109 TGeoBBox fast(
"repBenchFast", 1., 1., 1.);
6116 BOOST_CHECK_GT(
samples.insidePoints, 100);
6117 BOOST_CHECK_LT(
samples.insidePoints,
static_cast<long long>(
samples.points.size()) - 100);
6137 TGeoTube tube(
"repBenchTube", 0.3, 1., 2.);
6145 BOOST_REQUIRE_EQUAL(
a.points.size(),
b.points.size());
6146 for (
size_t i = 0;
i <
a.points.size(); ++
i) {
6163 constexpr size_t kBytes = 64u << 20;
6164 auto block = std::make_unique<char[]>(kBytes);
6165 for (
size_t i = 0;
i < kBytes;
i += 4096) {
6166 block[
i] =
static_cast<char>(
i);
6172 BOOST_CHECK_LT((
readMemory() - before).heapInUseBytes, 8LL << 20);
6181 auto* manager =
new TGeoManager(
"repBenchLadder",
"ladder");
6182 for (
const int k : {2, 4, 8, 16, 32}) {
6197 BOOST_CHECK(
dynamic_cast<TGeoCompositeShape*
>(single) ==
nullptr);
6200 gGeoManager =
nullptr;
6218struct SampleStream {
6219 explicit SampleStream(std::uint64_t seed) { lcg.state =
seed | 1u; }
6221 double uniform() {
return lcg.
next(); }
6222 double symmetric(
double extent) {
return (2. * uniform() - 1.) * extent; }
6232std::vector<std::array<double, 3>> nearestPatchSample(
const SurfaceSolid& solid,
double extent,
int count)
6234 SampleStream
stream(0x5EAFE7Full);
6235 std::vector<std::array<double, 3>> points;
6236 points.reserve(
static_cast<size_t>(
count));
6238 std::array<double, 3> point{
stream.symmetric(extent),
stream.symmetric(extent),
stream.symmetric(extent)};
6239 switch (
index % 5) {
6243 std::array<double, 3> normal{0., 0., 0.};
6244 const double safety = solid.Safety_Loop(point.data(), solid.Contains(point.data()));
6245 solid.ComputeNormal_Loop(point.data(),
nullptr, normal.data());
6246 const double sign = solid.Contains(point.data()) ? 1. : -1.;
6247 for (
int dimension = 0; dimension < 3; ++dimension) {
6248 point[dimension] += sign * safety * normal[dimension];
6253 std::array<double, 3> normal{0., 0., 0.};
6254 const double safety = solid.Safety_Loop(point.data(), solid.Contains(point.data()));
6255 solid.ComputeNormal_Loop(point.data(),
nullptr, normal.data());
6256 const double sign = solid.Contains(point.data()) ? 1. : -1.;
6257 const double offset = safety - sign * 1.e-9 * std::max(1., extent);
6258 for (
int dimension = 0; dimension < 3; ++dimension) {
6259 point[dimension] += sign *
offset * normal[dimension];
6264 for (
auto& coordinate : point) {
6269 for (
auto& coordinate : point) {
6274 points.push_back(point);
6277 points.push_back({0., 0., 0.});
6283int countNearestPatchDisagreements(
const SurfaceSolid& solid,
const std::array<double, 3>& point,
6284 double* worstSafetyGap =
nullptr)
6286 int disagreements = 0;
6287 for (
const bool inside : {
true,
false}) {
6288 const double accelerated = solid.Safety(point.data(), inside);
6289 const double reference = solid.Safety_Loop(point.data(), inside);
6293 if (worstSafetyGap !=
nullptr) {
6294 *worstSafetyGap = std::max(*worstSafetyGap, accelerated -
reference);
6299 const std::array<double, 3> direction = unitDirection(0.37, -0.82, 0.44);
6300 for (
const double* dir : {
static_cast<const double*
>(
nullptr), direction.data()}) {
6301 std::array<double, 3> accelerated{0., 0., 0.};
6302 std::array<double, 3>
reference{0., 0., 0.};
6303 solid.ComputeNormal(point.data(), dir, accelerated.data());
6304 solid.ComputeNormal_Loop(point.data(), dir,
reference.data());
6305 for (
int dimension = 0; dimension < 3; ++dimension) {
6306 if (accelerated[dimension] !=
reference[dimension]) {
6311 return disagreements;
6316std::unique_ptr<SurfaceSolid> makeManyPatchSolid(
const char*
name,
int ringCount)
6318 auto solid = std::make_unique<SurfaceSolid>(
name);
6319 for (
int ring = 0; ring < ringCount; ++ring) {
6321 addBoxSurfaces(*solid, 0.4, 0.4, 0.4, {3. * std::cos(
angle), 3. * std::sin(
angle), 0.});
6323 solid->CloseShape();
6329BoundaryCurve quarterCircleBoundary(
double cu,
double cv,
double r,
double a0)
6332 const double aMid = 0.5 * (a0 + a1);
6333 std::vector<Point2D> poles{{cu +
r * std::cos(a0), cv +
r * std::sin(a0)},
6334 {cu +
r * std::sqrt(2.) * std::cos(aMid), cv +
r * std::sqrt(2.) * std::sin(aMid)},
6335 {cu +
r * std::cos(a1), cv +
r * std::sin(a1)}};
6336 return BoundaryCurve::makeBSpline(2, std::move(poles), {1., std::sqrt(0.5), 1.}, {0., 0., 0., 1., 1., 1.});
6343std::unique_ptr<SurfaceSolid> makeWireTrimmedSolid(
const char*
name)
6345 auto solid = std::make_unique<SurfaceSolid>(
name);
6346 for (
int window = 0; window < 4; ++window) {
6348 BOOST_REQUIRE(solid->AddCylindricalSurface(
6349 {0., 0., 0.}, {0., 0., 1.}, {1., 0., 0.}, 2., -1., 1., 0.,
surf::kTwoPi,
false,
6350 {quarterCircleBoundary(centrePhi, 0., 0.5, 0.), quarterCircleBoundary(centrePhi, 0., 0.5, surf::kHalfPi),
6351 quarterCircleBoundary(centrePhi, 0., 0.5, surf::kPi),
6352 quarterCircleBoundary(centrePhi, 0., 0.5, 3. * surf::kHalfPi)}));
6354 BOOST_REQUIRE(addDiskSurface(*solid, {0., 0., 1.}, {1., 0., 0.}, {0., 1., 0.}, 2.));
6355 BOOST_REQUIRE(addDiskSurface(*solid, {0., 0., -1.}, {1., 0., 0.}, {0., -1., 0.}, 2.));
6356 solid->CloseShape();
6360struct NearestPatchFixture {
6361 std::unique_ptr<SurfaceSolid> solid;
6365std::vector<NearestPatchFixture> nearestPatchFixtures()
6367 std::vector<NearestPatchFixture> fixtures;
6368 fixtures.push_back({makeBoxSolid(
"safetyBox", 1., 2., 3.), 4.});
6369 fixtures.push_back({makeTubeSolid(
"safetyTube", 0., 2., 3.), 4.});
6370 fixtures.push_back({makeTubeSolid(
"safetyHollowTube", 1., 2., 3.), 4.});
6371 fixtures.push_back({makeConeSolid(
"safetyCone", 2., 1., 3.), 4.});
6372 fixtures.push_back({makeSphereSolid(
"safetySphere", 2.5), 3.5});
6373 fixtures.push_back({makeTorusSolid(
"safetyTorus", 3., 1.), 4.5});
6374 fixtures.push_back({makeCapsuleSolid(
"safetyCapsule", 1., 1.5), 3.});
6375 fixtures.push_back({makeManyPatchSolid(
"safetyRing", 12), 4.5});
6376 fixtures.push_back({makeWireTrimmedSolid(
"safetyWireTrim"), 3.});
6403 const std::vector<Vec2> rectangle{{0., 0.}, {2., 0.}, {2., 3.}, {0., 3.}};
6404 BOOST_REQUIRE(polygon.
initialize({0., 0., 0.}, {1., 0., 0.}, {0., 1., 0.}, rectangle, {}, error));
6406 BOOST_REQUIRE(disk.
initialize({0., 0., 0.5}, {1., 0., 0.}, {0., 1., 0.},
6407 {surf::Curve2D::makeCircle({0., 0.}, 1.5)}, {},
error));
6412 BOOST_REQUIRE(splineFace.
initialize({0.2, -0.3, 1.1}, {1., 0., 0.}, {0., 1., 0.},
6413 {quarterCircleBSpline(0., 0., 1.2, 0.),
6414 quarterCircleBSpline(0., 0., 1.2, surf::kHalfPi),
6415 quarterCircleBSpline(0., 0., 1.2, surf::kPi),
6416 quarterCircleBSpline(0., 0., 1.2, 3. * surf::kHalfPi)},
6419 BOOST_REQUIRE(cylinder.
initialize({0.1, -0.2, 0.}, {0., 0., 1.}, {1., 0., 0.}, 2., -1., 1., 0.3,
6422 BOOST_REQUIRE(cone.
initialize({0., 0., 0.}, {0., 1., 0.}, {1., 0., 0.}, 2., 0.5, -1., 1., 0., 1.1 *
surf::kPi,
6425 BOOST_REQUIRE(sphere.
initialize({0.3, 0.4, -0.5}, {0., 0., 1.}, {1., 0., 0.}, 1.7, 0.2, 2.4, 0.,
6428 BOOST_REQUIRE(torus.
initialize({0., 0., 0.}, {0., 0., 1.}, {1., 0., 0.}, 3., 0.8, 0., 1.4 *
surf::kPi, 0.,
6431 const std::vector<const surf::BoundedSurface*> surfaces{&polygon, &disk, &splineFace,
6432 &cylinder, &cone, &sphere,
6435 SampleStream
stream(0xB0B0Dull);
6437 for (
const auto* surface : surfaces) {
6438 Vec3
lower{TGeoShape::Big(), TGeoShape::Big(), TGeoShape::Big()};
6439 Vec3
upper{-TGeoShape::Big(), -TGeoShape::Big(), -TGeoShape::Big()};
6441 const double boxLower[3] = {
lower.xCoord,
lower.yCoord,
lower.zCoord};
6442 const double boxUpper[3] = {
upper.xCoord,
upper.yCoord,
upper.zCoord};
6444 for (
int sample = 0; sample < 4000; ++sample) {
6445 const double scale = (sample % 4 == 3) ? 1.e6 : ((sample % 4 == 2) ? 20. : 5.);
6446 const Vec3 point{
stream.symmetric(scale),
stream.symmetric(scale),
stream.symmetric(scale)};
6447 const double coordinates[3] = {point.xCoord, point.yCoord, point.zCoord};
6448 double boxDistanceSq = 0.;
6449 for (
int dimension = 0; dimension < 3; ++dimension) {
6450 if (coordinates[dimension] < boxLower[dimension]) {
6451 const double gap = boxLower[dimension] - coordinates[dimension];
6452 boxDistanceSq += gap * gap;
6453 }
else if (coordinates[dimension] > boxUpper[dimension]) {
6454 const double gap = coordinates[dimension] - boxUpper[dimension];
6455 boxDistanceSq += gap * gap;
6462 BOOST_REQUIRE_LE(boxDistanceSq * (1. - 1.e-12), patchDistanceSq);
6472 size_t comparedPoints = 0;
6473 double worstSafetyGap = -std::numeric_limits<double>::infinity();
6474 for (
const auto& fixture : nearestPatchFixtures()) {
6475 BOOST_TEST_CONTEXT(
"fixture = " << fixture.solid->GetName())
6477 BOOST_REQUIRE(fixture.solid->HasBVH());
6478 for (
const auto& point : nearestPatchSample(*fixture.solid, fixture.extent, 2000)) {
6479 BOOST_TEST_CONTEXT(
"point = (" << point[0] <<
", " << point[1] <<
", " << point[2] <<
")")
6481 BOOST_REQUIRE_EQUAL(countNearestPatchDisagreements(*fixture.solid, point, &worstSafetyGap), 0);
6497 SurfaceSolid open(
"safetyBeforeClose");
6498 addBoxSurfaces(open, 1., 2., 3.);
6499 BOOST_REQUIRE(!open.HasBVH());
6500 const std::array<double, 3> probe{0.3, -1.1, 2.2};
6501 BOOST_CHECK_EQUAL(open.Safety(probe.data(), kTRUE), open.Safety_Loop(probe.data(), kTRUE));
6502 std::array<double, 3> viaBVH{0., 0., 0.};
6503 std::array<double, 3> viaLoop{0., 0., 0.};
6504 open.ComputeNormal(probe.data(),
nullptr, viaBVH.data());
6505 open.ComputeNormal_Loop(probe.data(),
nullptr, viaLoop.data());
6510 SurfaceSolid
empty(
"safetyEmpty");
6514 empty.ComputeNormal(probe.data(),
nullptr, viaBVH.data());
6515 empty.ComputeNormal_Loop(probe.data(),
nullptr, viaLoop.data());
6530 const auto box = makeBoxSolid(
"tieBreakBox", 2., 2., 2.);
6531 const std::array<double, 3> centre{0., 0., 0.};
6532 std::array<double, 3> viaBVH{0., 0., 0.};
6533 std::array<double, 3> viaLoop{0., 0., 0.};
6534 box->ComputeNormal(centre.data(),
nullptr, viaBVH.data());
6535 box->ComputeNormal_Loop(centre.data(),
nullptr, viaLoop.data());
6542 BOOST_CHECK_EQUAL(std::abs(viaLoop[0]) + std::abs(viaLoop[1]) + std::abs(viaLoop[2]), 1.);
6545 const auto tube = makeTubeSolid(
"tieBreakTube", 1., 2., 1.);
6546 tube->ComputeNormal(centre.data(),
nullptr, viaBVH.data());
6547 tube->ComputeNormal_Loop(centre.data(),
nullptr, viaLoop.data());
6565 BOOST_REQUIRE(!SurfaceSolid::GetSafetyBoundUnsoundForTest());
6567 const auto fixtures = nearestPatchFixtures();
6568 size_t caughtOnFixtures = 0;
6569 size_t prunableFixtures = 0;
6570 size_t sabotagedDisagreements = 0;
6571 size_t safetyTooLarge = 0;
6572 size_t safetyTooSmall = 0;
6573 for (
const auto& fixture : fixtures) {
6579 const auto points = nearestPatchSample(*fixture.solid, fixture.extent, 400);
6580 size_t disagreementsHere = 0;
6581 SurfaceSolid::SetSafetyBoundUnsoundForTest(
true);
6582 for (
const auto& point : points) {
6583 const double sabotaged = fixture.solid->Safety(point.data(), kTRUE);
6584 SurfaceSolid::SetSafetyBoundUnsoundForTest(
false);
6585 const double reference = fixture.solid->Safety_Loop(point.data(), kTRUE);
6586 SurfaceSolid::SetSafetyBoundUnsoundForTest(
true);
6588 ++disagreementsHere;
6589 (sabotaged >
reference) ? ++safetyTooLarge : ++safetyTooSmall;
6591 disagreementsHere +=
static_cast<size_t>(countNearestPatchDisagreements(*fixture.solid, point));
6593 SurfaceSolid::SetSafetyBoundUnsoundForTest(
false);
6594 sabotagedDisagreements += disagreementsHere;
6595 if (disagreementsHere > 0) {
6597 BOOST_TEST_MESSAGE(
"sabotaged bound caught on " << fixture.solid->GetName() <<
": " << disagreementsHere
6598 <<
" disagreements over " << points.size() <<
" points");
6600 BOOST_TEST_MESSAGE(
"sabotaged bound NOT caught on " << fixture.solid->GetName() <<
" ("
6601 << fixture.solid->GetNsurfaces() <<
" patches)");
6608 BOOST_CHECK_GT(sabotagedDisagreements, 100u);
6610 BOOST_CHECK_GT(safetyTooLarge, 0u);
6615 for (
const auto& fixture : fixtures) {
6616 for (
const auto& point : nearestPatchSample(*fixture.solid, fixture.extent, 400)) {
6617 BOOST_REQUIRE_EQUAL(countNearestPatchDisagreements(*fixture.solid, point), 0);
6620 BOOST_CHECK(!SurfaceSolid::GetSafetyBoundUnsoundForTest());
6629 const auto ring = makeManyPatchSolid(
"candidateRing", 24);
6630 BOOST_REQUIRE_EQUAL(ring->GetNsurfaces(), 144);
6631 const auto points = nearestPatchSample(*ring, 5., 500);
6633 SurfaceSolid::ResetSafetyCandidateCounter();
6634 for (
const auto& point : points) {
6635 ring->Safety(point.data(), kTRUE);
6637 const long long acceleratedCandidates = SurfaceSolid::GetSafetyCandidateCount();
6638 const double perCall =
static_cast<double>(acceleratedCandidates) / points.size();
6640 BOOST_TEST_MESSAGE(
"Safety candidates per call: " << perCall <<
" of " << ring->GetNsurfaces() <<
" patches");
6641 BOOST_CHECK_GT(acceleratedCandidates, 0);
6642 BOOST_CHECK_LT(perCall, 0.4 * ring->GetNsurfaces());
6645 SurfaceSolid::ResetSafetyCandidateCounter();
6646 for (
const auto& point : points) {
6647 ring->Safety_Loop(point.data(), kTRUE);
6648 std::array<double, 3> normal{0., 0., 0.};
6649 ring->ComputeNormal_Loop(point.data(),
nullptr, normal.data());
6655 SurfaceSolid::ResetSafetyCandidateCounter();
6656 for (
const auto& point : points) {
6657 std::array<double, 3> normal{0., 0., 0.};
6658 ring->ComputeNormal(point.data(),
nullptr, normal.data());
6660 const double normalPerCall =
static_cast<double>(SurfaceSolid::GetSafetyCandidateCount()) / points.size();
6661 BOOST_TEST_MESSAGE(
"ComputeNormal candidates per call: " << normalPerCall);
6662 BOOST_CHECK_LT(normalPerCall, 0.4 * ring->GetNsurfaces());
6663 BOOST_CHECK_GE(normalPerCall, perCall);
6678using CoverBox = surf::BoundedSurface::CoverBox;
6682double coverBoxDistanceSq(
const CoverBox&
box,
const surf::Vec3& point)
6686 const double lower[3] = {
box.first.xCoord,
box.first.yCoord,
box.first.zCoord};
6687 const double upper[3] = {
box.second.xCoord,
box.second.yCoord,
box.second.zCoord};
6688 for (
int dimension = 0; dimension < 3; ++dimension) {
6689 const double gap = std::max({
lower[dimension] - coordinates[dimension],
6690 coordinates[dimension] -
upper[dimension], 0.});
6696double minCoverBoxDistanceSq(
const std::vector<CoverBox>& boxes,
const surf::Vec3& point)
6698 double best = std::numeric_limits<double>::infinity();
6699 for (
const auto&
box : boxes) {
6700 best = std::min(best, coverBoxDistanceSq(
box, point));
6705bool anyCoverBoxContains(
const std::vector<CoverBox>& boxes,
const surf::Vec3& point,
double slack)
6707 for (
const auto&
box : boxes) {
6708 if (point.
xCoord >=
box.first.xCoord - slack && point.
xCoord <=
box.second.xCoord + slack &&
6709 point.
yCoord >=
box.first.yCoord - slack && point.
yCoord <=
box.second.yCoord + slack &&
6710 point.
zCoord >=
box.first.zCoord - slack && point.
zCoord <=
box.second.zCoord + slack) {
6721void checkCoverBoxProperties(
const surf::BoundedSurface& surface,
const std::vector<surf::Vec3>& patchPoints,
6724 std::vector<CoverBox> boxes;
6726 BOOST_TEST_CONTEXT(
"surface = " <<
label)
6728 BOOST_REQUIRE(!boxes.empty());
6729 for (
const auto& point : patchPoints) {
6730 BOOST_TEST_CONTEXT(
"patch point = (" << point.
xCoord <<
", " << point.
yCoord <<
", " << point.
zCoord <<
")")
6732 BOOST_CHECK(anyCoverBoxContains(boxes, point, 1.e-9));
6735 SampleStream
stream(0xC0FEB0C5ull);
6738 constexpr double kProbeScales[3] = {1.5, 8., 300.};
6740 const double scale = kProbeScales[
index % 3];
6743 const double boxDistanceSq = minCoverBoxDistanceSq(boxes, point);
6744 BOOST_TEST_CONTEXT(
"point = (" << point.
xCoord <<
", " << point.
yCoord <<
", " << point.
zCoord <<
")")
6746 BOOST_CHECK_LE(boxDistanceSq, patchDistanceSq * (1. + 1.e-9) + 1.e-18);
6756 std::vector<CoverBox> boxes;
6758 BOOST_TEST_CONTEXT(
"surface = " <<
label)
6760 BOOST_CHECK_GT(boxes.size(), 1u);
6774 const Vec3 skewCenter{0.4, -0.2, 0.1};
6775 const Vec3 skewAxis{0.2, 0.3, 1.};
6776 const Vec3 referenceU{1., 0., 0.};
6780 BOOST_REQUIRE(cylinder.
initialize(skewCenter, skewAxis, referenceU, 1.7, -0.8, 1.2, 0.4, 1.9,
false, error));
6781 std::vector<Vec3> patchPoints;
6782 for (
int stepPhi = 0; stepPhi <= 12; ++stepPhi) {
6783 for (
int stepH = 0; stepH <= 4; ++stepH) {
6784 patchPoints.push_back(cylinder.
pointAt(0.4 + 1.9 * stepPhi / 12., -0.8 + 2. * stepH / 4.));
6787 checkCoverBoxProperties(cylinder, patchPoints,
"partial cylinder, skew axis");
6788 checkEmitsSeveralCoverBoxes(cylinder,
"partial cylinder, skew axis");
6792 BOOST_REQUIRE(cylinder.
initialize({0., 0., 0.}, {0., 0., 1.}, referenceU, 2., -1., 1., 0.,
surf::kTwoPi,
6794 std::vector<Vec3> patchPoints;
6795 for (
int stepPhi = 0; stepPhi <= 24; ++stepPhi) {
6796 patchPoints.push_back(cylinder.
pointAt(
surf::kTwoPi * stepPhi / 24., -1. + 2. * (stepPhi % 5) / 4.));
6798 checkCoverBoxProperties(cylinder, patchPoints,
"full cylinder");
6799 checkEmitsSeveralCoverBoxes(cylinder,
"full cylinder");
6803 BOOST_REQUIRE(cylinder.
initialize({0., 0., 0.}, {0., 0., 1.}, referenceU, 2., -1., 1., 0.,
surf::kTwoPi,
false,
6804 paramRectWireCurves(0.3, 2.1, -0.5, 0.7), {}, error));
6805 std::vector<Vec3> patchPoints;
6806 for (
int stepPhi = 0; stepPhi <= 10; ++stepPhi) {
6807 for (
int stepH = 0; stepH <= 4; ++stepH) {
6808 const double phi = 0.3 + 1.8 * stepPhi / 10.;
6809 const double height = -0.5 + 1.2 * stepH / 4.;
6815 BOOST_REQUIRE(!patchPoints.empty());
6816 checkCoverBoxProperties(cylinder, patchPoints,
"wire-trimmed cylinder");
6822 std::vector<Vec3> patchPoints;
6823 for (
int stepTheta = 0; stepTheta <= 8; ++stepTheta) {
6824 for (
int stepPhi = 0; stepPhi < 16; ++stepPhi) {
6828 checkCoverBoxProperties(sphere, patchPoints,
"full sphere, skew frame");
6829 checkEmitsSeveralCoverBoxes(sphere,
"full sphere, skew frame");
6835 BOOST_REQUIRE(cap.
initialize({0., 0., 0.}, {0., 0., 1.}, referenceU, 2., 0., 0.6, 0.2, 1.1,
false, error));
6836 std::vector<Vec3> patchPoints;
6837 for (
int stepTheta = 0; stepTheta <= 4; ++stepTheta) {
6838 for (
int stepPhi = 0; stepPhi <= 6; ++stepPhi) {
6839 patchPoints.push_back(cap.
pointAt(0.6 * stepTheta / 4., 0.2 + 1.1 * stepPhi / 6.));
6842 checkCoverBoxProperties(cap, patchPoints,
"spherical cap");
6846 BOOST_REQUIRE(cone.
initialize(skewCenter, skewAxis, referenceU, 2., 0.5, -0.9, 1.1, 0.7, 2.3,
false, error));
6847 std::vector<Vec3> patchPoints;
6848 for (
int stepPhi = 0; stepPhi <= 10; ++stepPhi) {
6849 for (
int stepH = 0; stepH <= 4; ++stepH) {
6850 patchPoints.push_back(cone.
pointAt(0.7 + 2.3 * stepPhi / 10., -0.9 + 2. * stepH / 4.));
6853 checkCoverBoxProperties(cone, patchPoints,
"partial cone, skew axis");
6854 checkEmitsSeveralCoverBoxes(cone,
"partial cone, skew axis");
6858 BOOST_REQUIRE(torus.
initialize(skewCenter, skewAxis, referenceU, 2.4, 0.7, 0.3, 2.1, -0.4, 1.7,
false, error));
6859 std::vector<Vec3> patchPoints;
6860 for (
int stepRing = 0; stepRing <= 10; ++stepRing) {
6861 for (
int stepTube = 0; stepTube <= 6; ++stepTube) {
6862 patchPoints.push_back(torus.
pointAt(0.3 + 2.1 * stepRing / 10., -0.4 + 1.7 * stepTube / 6.));
6865 checkCoverBoxProperties(torus, patchPoints,
"partial torus, skew axis");
6866 checkEmitsSeveralCoverBoxes(torus,
"partial torus, skew axis");
6870 BOOST_REQUIRE(torus.
initialize({0., 0., 0.}, {0., 0., 1.}, referenceU, 3., 1., 0.,
surf::kTwoPi, 0.,
6872 std::vector<Vec3> patchPoints;
6873 for (
int stepRing = 0; stepRing < 16; ++stepRing) {
6874 for (
int stepTube = 0; stepTube < 8; ++stepTube) {
6878 checkCoverBoxProperties(torus, patchPoints,
"full torus");
6882 const std::vector<surf::Vec2> rectangle{{0., 0.}, {2., 0.}, {2., 3.}, {0., 3.}};
6883 BOOST_REQUIRE(polygon.
initialize({0.2, -0.4, 0.5}, {1., 0.2, 0.}, {-0.1, 1., 0.3}, rectangle, {}, error));
6884 std::vector<Vec3> patchPoints;
6885 patchPoints.push_back(polygon.
toGlobal({0.01, 0.01}));
6886 patchPoints.push_back(polygon.
toGlobal({1.9, 2.9}));
6887 checkCoverBoxProperties(polygon, patchPoints,
"planar polygon");
6898 const auto sphere = makeSphereSolid(
"subBoxSphere", 2.5);
6899 BOOST_CHECK_EQUAL(sphere->CountBVHRayCandidates({2.2, 2.2, -5.}, {0., 0., 1.}), 0);
6900 BOOST_CHECK_GE(sphere->CountBVHRayCandidates({0., 0., -5.}, {0., 0., 1.}), 1);
6903 const auto tube = makeTubeSolid(
"subBoxTube", 0., 2., 1.);
6907 const auto torus = makeTorusSolid(
"subBoxTorus", 3., 1.);
6908 BOOST_CHECK_EQUAL(torus->CountBVHRayCandidates({3.4, 3.4, -5.}, {0., 0., 1.}), 0);
6909 BOOST_CHECK_GE(torus->CountBVHRayCandidates({3., 0., -5.}, {0., 0., 1.}), 1);
6913 SurfaceSolid quarter(
"subBoxQuarterCylinder");
6914 BOOST_REQUIRE(quarter.AddCylindricalSurface({0., 0., 0.}, {0., 0., 1.}, {1., 0., 0.}, 2., -1., 1.,
6916 quarter.CloseShape(
false);
6917 BOOST_REQUIRE(quarter.HasBVH());
6918 BOOST_CHECK_EQUAL(quarter.CountBVHRayCandidates({-1.9, -5., 0.}, {0., 1., 0.}), 0);
6919 BOOST_CHECK_GE(quarter.CountBVHRayCandidates({5., 0., 0.}, {-1., 0., 0.}), 1);
6930 std::unique_ptr<SurfaceSolid> solid;
6933 std::vector<Fixture> fixtures;
6934 fixtures.push_back({makeSphereSolid(
"subBoxSweepSphere", 2.5), 3.5});
6935 fixtures.push_back({makeTorusSolid(
"subBoxSweepTorus", 3., 1.), 4.5});
6936 fixtures.push_back({makeCapsuleSolid(
"subBoxSweepCapsule", 1., 1.5), 3.});
6937 fixtures.push_back({makeConeSolid(
"subBoxSweepCone", 2., 1., 3.), 4.});
6938 fixtures.push_back({makeWireTrimmedSolid(
"subBoxSweepWireTrim"), 3.});
6940 for (
const auto& fixture : fixtures) {
6941 BOOST_TEST_CONTEXT(
"fixture = " << fixture.solid->GetName())
6943 sweepDistanceAgainstLoop(*fixture.solid, fixture.extent, 4);
6944 for (
const auto& point : probeGrid(fixture.extent, 4)) {
6945 BOOST_TEST_CONTEXT(
"point = (" << point[0] <<
", " << point[1] <<
", " << point[2] <<
")")
6947 BOOST_CHECK_EQUAL(fixture.solid->Contains(point.data()), fixture.solid->Contains_Loop(point.data()));
6948 std::vector<SurfaceSolid::ContainsCrossing> bvhCrossings;
6949 std::vector<SurfaceSolid::ContainsCrossing> loopCrossings;
6950 fixture.solid->DescribeContainsCrossings({point[0], point[1], point[2]}, bvhCrossings, loopCrossings);
6951 BOOST_REQUIRE_EQUAL(bvhCrossings.size(), loopCrossings.size());
header::DataOrigin origin
o2::monitoring::tags::Key Key
std::unique_ptr< expressions::Node > node
o2::raw::RawFileWriter * raw
Validation and timing harness for TGeoShape navigation, typed on plain TGeoShape*.
std::vector< SidecarEdge > edges
Per-call cost, memory and the synthetic boolean ladder: the measuring parts of the representation com...
The X-ray transport benchmark's algorithms: stepping, auditing and comparing ordered crossing lists.
Class for time synchronization of RawReader instances.
bool AddCylindricalSurface(const Point3D ¢erPoint, const Point3D &axis, const Point3D &referenceAxisU, double radius, double heightMin, double heightMax, double phiStart=0., double phiSweep=6.283185307179586, bool innerWall=false)
Add a cylindrical wall of radius around axis over a height range and a phi sweep; innerWall points th...
void GetSurfaceCapacityContributions(std::vector< double > &contributions) const
Each face's divergence-theorem contribution to Capacity(), in record order.
Abstract analytic surface patch: one support surface plus its trim, with the kernels the navigation n...
double parametricLengthSqAt(const Vec2 &uv, const Vec2 &delta) const
The 3D length squared spanned by a parametric displacement delta starting at uv.
virtual void parametricMetric(const Vec2 &uv, double &gUU, double &gUV, double &gVV) const =0
The first fundamental form at uv, turning parametric displacements into 3D lengths; it varies over th...
virtual void appendCoverBoxes(std::vector< CoverBox > &boxes) const
virtual double distanceSqToPatch(const Vec3 &point) const =0
Squared distance from a 3D point to the trimmed patch (used for Safety).
A cone whose radius varies linearly with height, trimmed as the cylinder; one radius may be zero (an ...
Vec3 pointAt(double phi, double height) const
bool initialize(const Vec3 ¢erPoint, const Vec3 &axis, const Vec3 &referenceAxisU, double radiusAtMin, double radiusAtMax, double heightMin, double heightMax, double phiStart, double phiSweep, bool innerWall, std::string &errorMessage)
void parametricMetric(const Vec2 &uv, double &gUU, double &gUV, double &gVV) const override
(u, v) = (phi[rad], h[cm]): the azimuthal scale is the local radius, and a step in h spans sqrt(1 + s...
A plane trimmed by curved (line/arc/B-spline) loops in an orthonormal frame: exact caps,...
bool initialize(const Vec3 &surfaceOrigin, const Vec3 &surfaceAxisU, const Vec3 &surfaceAxisV, const std::vector< Curve2D > &outerCurves, const std::vector< std::vector< Curve2D > > &innerCurves, std::string &errorMessage, double joinTolerance=kWireJoinTolerance)
A cylinder of given radius around an axis, trimmed to a (phi, h) rectangle or by curve wires; innerWa...
Vec3 pointAt(double phi, double height) const
Vec3 normalAt(const Vec3 &point) const override
Outward-oriented normal at (or nearest to) the given point.
bool capacityIsExact() const override
Whether capacityContribution() is analytically exact for this surface.
void appendIntersections(const Vec3 &rayOrigin, const Vec3 &rayDirection, double minDistance, double maxDistance, std::vector< RayHit > &hits) const override
Append every hit of the ray with the trimmed patch in [minDistance, maxDistance], with the outward no...
bool pointInTrim(double phi, double height, bool *boundary=nullptr) const
True if the (phi, h) point lies in the trim wire (phi unwrapped into the wire window).
void conservativeBounds(Vec3 &lower, Vec3 &upper) const override
Accumulate a conservative axis-aligned bounding box of the trimmed patch.
bool initialize(const Vec3 ¢erPoint, const Vec3 &axis, const Vec3 &referenceAxisU, double radius, double heightMin, double heightMax, double phiStart, double phiSweep, bool innerWall, std::string &errorMessage)
bool containsPointOnSurface(const Vec3 &point) const override
True if the 3D point lies on the trimmed patch within tolerance.
double distanceSqToPatch(const Vec3 &point) const override
Distance to the patch: exact for the parametric rectangle, a lower bound for a wire trim.
double capacityContribution() const override
Divergence-theorem contribution over the (phi, h) rectangle; a wire trim uses the contour form,...
A trivial bounded surface, a single 3D triangle, to exercise the BoundedSurface interface.
bool capacityIsExact() const override
Whether capacityContribution() is analytically exact for this surface.
DummyBoundedSurface(const Vec3 &firstVertex, const Vec3 &secondVertex, const Vec3 &thirdVertex)
void conservativeBounds(Vec3 &lower, Vec3 &upper) const override
Accumulate a conservative axis-aligned bounding box of the trimmed patch.
void appendDirectedEdges(std::vector< std::pair< Vec3, Vec3 > > &edges) const override
Append the 3D directed boundary edges of the patch, for solid-closure validation.
Vec3 normalAt(const Vec3 &) const override
Outward-oriented normal at (or nearest to) the given point.
double capacityContribution() const override
Signed divergence-theorem contribution to the enclosed volume.
bool containsPointOnSurface(const Vec3 &) const override
True if the 3D point lies on the trimmed patch within tolerance.
void appendDisplayMesh(std::vector< Vec3 > &vertices, std::vector< std::array< int, 3 > > &triangles) const override
Append this patch's visualization triangulation (navigation must never depend on it).
double distanceSqToPatch(const Vec3 &point) const override
Squared distance from a 3D point to the trimmed patch (used for Safety).
void parametricMetric(const Vec2 &, double &gUU, double &gUV, double &gVV) const override
A triangle carries no parametric domain, so the form is the identity.
void appendIntersections(const Vec3 &, const Vec3 &, double, double, std::vector< RayHit > &) const override
Append every hit of the ray with the trimmed patch in [minDistance, maxDistance], with the outward no...
bool initialize(const Vec3 &surfaceOrigin, const Vec3 &surfaceAxisU, const Vec3 &surfaceAxisV, const std::vector< Vec2 > &outerWireVertices, const std::vector< std::vector< Vec2 > > &innerWireVertices, std::string &errorMessage)
Vec3 toGlobal(const Vec2 &point) const
void parametricMetric(const Vec2 &, double &gUU, double &gUV, double &gVV) const override
Constant over the plane, with a cross term: the frame axes need be neither unit-length nor orthogonal...
bool containsLocal(const Vec2 &point, bool *boundary=nullptr) const
A sphere of given radius trimmed to a (theta, phi) rectangle or by curve wires; innerWall points the ...
bool initialize(const Vec3 ¢er, const Vec3 &polarAxis, const Vec3 &referenceAxisU, double radius, double thetaMin, double thetaMax, double phiStart, double phiSweep, bool innerWall, std::string &errorMessage)
Vec3 pointAt(double theta, double phi) const
void appendIntersections(const Vec3 &rayOrigin, const Vec3 &rayDirection, double minDistance, double maxDistance, std::vector< RayHit > &hits) const override
Append every hit of the ray with the trimmed patch in [minDistance, maxDistance], with the outward no...
double distanceSqToPatch(const Vec3 &point) const override
Distance to the patch: exact inside the trim, else the full-sphere distance, a lower bound.
void parametricMetric(const Vec2 &uv, double &gUU, double &gUV, double &gVV) const override
(u, v) = (phi[rad], theta[rad]); gUU vanishes at either pole.
bool containsPointOnSurface(const Vec3 &point) const override
True if the 3D point lies on the trimmed patch within tolerance.
double capacityContribution() const override
Divergence-theorem contribution over the (theta, phi) rectangle; a wire trim uses the contour form in...
bool containsPointOnSurface(const Vec3 &point) const override
True if the 3D point lies on the trimmed patch within tolerance.
Vec3 pointAt(double phiRing, double phiTube) const
bool initialize(const Vec3 ¢erPoint, const Vec3 &axis, const Vec3 &referenceAxisU, double majorRadius, double minorRadius, double phiStart, double phiSweep, double tubeStart, double tubeSweep, bool innerWall, std::string &errorMessage)
void appendIntersections(const Vec3 &rayOrigin, const Vec3 &rayDirection, double minDistance, double maxDistance, std::vector< RayHit > &hits) const override
Append every hit of the ray with the trimmed patch in [minDistance, maxDistance], with the outward no...
bool capacityIsExact() const override
Whether capacityContribution() is analytically exact for this surface.
Vec3 normalAt(const Vec3 &point) const override
Outward-oriented normal at (or nearest to) the given point.
double distanceSqToPatch(const Vec3 &point) const override
Distance to the patch: exact for the full torus by the meridian distance, a lower bound for a trimmed...
double capacityContribution() const override
Divergence-theorem contribution over the (phiRing, phiTube) rectangle; a wire trim uses the contour f...
void parametricMetric(const Vec2 &uv, double &gUU, double &gUV, double &gVV) const override
(u, v) = (phiRing[rad], phiTube[rad]): the tube scale is r, the ring scale the distance from the axis...
GLenum GLsizeiptr const void GLsizei faceIndex
GLsizei const GLuint const GLfloat * weights
GLuint const GLchar * name
GLint GLsizei GLsizei height
GLboolean GLboolean GLboolean b
GLsizei GLsizei GLchar * source
GLsizei GLsizei GLfloat distance
GLsizei const GLfloat * value
GLenum GLsizei GLsizei GLint * values
GLuint GLsizei GLsizei * length
GLuint GLsizei const GLchar * label
GLint GLenum GLboolean normalized
GLsizei const GLchar *const * path
GLfloat GLfloat GLfloat GLfloat v3
GLboolean GLboolean GLboolean GLboolean a
GLubyte GLubyte GLubyte GLubyte w
GLfloat GLfloat GLfloat v2
GLsizei const GLint * box
GLdouble GLdouble GLdouble z
void report(gsl::span< o2::InteractionTimeRecord > irs, int threshold, bool verbose)
QuerySamples buildQuerySamples(const TGeoShape *reference, const std::string &referenceName, const Point3D &bboxMin, const Point3D &bboxMax, int nPoints, int nRays, uint64_t seed=20260802ULL, double inflate=0.12)
TGeoShape * buildBooleanLadder(int leaves, LadderShape shape, const std::string &tag)
TimingStat timeContainsPass(const TGeoShape *shape, const QuerySamples &s, int warmup, int passes)
MemorySnapshot readMemory()
TimingStat timeDistInPass(const TGeoShape *shape, const QuerySamples &s, int warmup, int passes)
TimingStat timeDistOutPass(const TGeoShape *shape, const QuerySamples &s, int warmup, int passes)
BooleanTreeStats booleanTreeStats(const TGeoShape *shape)
std::array< double, 3 > Point3D
double distanceSq(const Vec2 &firstPoint, const Vec2 &secondPoint)
QuarticBranch
Which of solveQuarticReal's branches produced its roots, for the tests.
constexpr double kBSplineFlatness
Chord flatness of the adaptive B-spline sampler, in the curve's parametric units; a B-spline trim is ...
WireClassification
Classification of a parametric point against a closed wire.
double pointSegmentDistanceSq(const Vec2 &point, const Vec2 &segmentStart, const Vec2 &segmentEnd)
double dot(const Vec3 &firstVector, const Vec3 &secondVector)
double integrateOverCurveTrim(const CurveWire &outerWire, const std::vector< CurveWire > &innerWires, const Integrand &integrand, int samplesPerAxis=128)
Midpoint-rule integral of integrand over the trimmed region; kept as the independent check of the con...
bool sameIntersection(double firstDistance, double secondDistance)
ParametricMetric parametricMetricOf(const Surface &surface)
bool buildCurveTrim(const std::vector< Curve2D > &outerTrim, const std::vector< std::vector< Curve2D > > &innerTrims, CurveWire &outerWire, std::vector< CurveWire > &innerWires, Vec2 &lower, Vec2 &upper, std::string &errorMessage, const ParametricMetric &metric={}, double joinTolerance=kWireJoinTolerance)
Build validated outer and inner trim wires and the outer loop's parametric bounds; rejects a trim wid...
constexpr double kBVHBoxTolerance
Widening of the BVH leaf boxes before the outward float rounding; it dominates every navigation lengt...
ClosureReport validateClosure(const std::vector< std::unique_ptr< BoundedSurface > > &surfaces, double modelTolerance=0.)
Validate closure and orientation of surfaces by half-edges, measure the rims, and count edge identiti...
QuarticRoots solveQuarticReal(double a4, double a3, double a2, double a1, double a0, QuarticBranch *takenBranch=nullptr)
Vec3 cross(const Vec3 &firstVector, const Vec3 &secondVector)
constexpr double kRimMatchTolerance
Rim-matching distance in cm when the model states no tolerance: the extractor precision,...
Vec3 normalized(const Vec3 &vector)
constexpr double kTolerance
generic length tolerance
constexpr int kArcSamples
Chords per full-circle arc for display and rims, shared by all surfaces so shared rims match; divisib...
constexpr double kIntersectionTolerance
clustering of near-equal intersections
std::vector< Crossing > stepWithShapeApi(const TGeoShape *shape, const Point3D &origin, const Point3D &dir, double tMax, const StepConfig &cfg, Robustness &stats)
Mode (a): the same loop driven by the ordinary TGeoShape virtuals.
void compareLists(const std::vector< Crossing > &candidate, const std::vector< Crossing > &reference, const Point3D &origin, const Point3D &dir, double tolerance, ListComparison &out)
double chordVolume(const Raster &raster, const std::vector< double > &insideLengthPerBeam)
std::vector< Crossing > stepCrossingsWithKernels(const Point3D &origin, const Point3D &dir, double tMax, const StepConfig &cfg, Robustness &stats, ContainsFn contains, DistOutFn distFromOutside, DistInFn distFromInside)
void auditCrossingList(const std::vector< Crossing > &crossings, const TGeoShape *shape, const Point3D &origin, const Point3D &dir, double tMax, const StepConfig &cfg, Robustness &stats)
double dot3(const Point3D &a, const Point3D &b)
std::vector< Beam > buildFanBeams(int count)
Raster buildRaster(const Point3D &bboxMin, const Point3D &bboxMax, int n, const std::vector< Beam > &beams, double transverseMargin)
std::vector< Beam > buildBeams(const std::string &axesSpec, double tiltDegrees)
bool LoadSurfaceSolid(const std::string &file, O2BVHSurfaceSolid &solid)
int32_t const char * file
int const float const TrackSeed< NLayers > & seed
const bool const bool const int FollowDirection BestTrial TrackITSInternal< NLayers > & best
ROOT::Math::PositionVector2D< ROOT::Math::Cartesian2D< T >, ROOT::Math::DefaultCoordinateSystemTag > Point2D
std::string to_string(gsl::span< T, Size > span)
static PlanarBoundaryCurve makeLine(const Point2D &start, const Point2D &end)
static PlanarBoundaryCurve makeBSpline(int splineDegree, std::vector< Point2D > splinePoles, std::vector< double > splineWeights, std::vector< double > splineKnots)
long long nodes
TGeoCompositeShape / TGeoBoolNode pairs.
Whether a set of bounded surfaces forms a closed, consistently oriented 2-manifold,...
One trimmed boundary curve in a surface's (u, v) domain: a line segment, a circular arc or a clamped ...
std::vector< Vec2 > poles
double uVariation(double from, double to) const
An upper bound on how far u travels along the curve between from and to.
Vec2 pointAt(double parameter) const
Point at curve parameter parameter in [0, 1] (0 at the start, 1 at the end).
One closed, oriented boundary loop of Curve2D segments: outer loops wind counter-clockwise,...
double representationTolerance() const
The widest gap between the loop's representation and its boundary, in parametric units; 0 for lines a...
void tightParametricBounds(Vec2 &lower, Vec2 &upper) const
Add the loop's extent measured on the curves to a parametric bounding box; use it to reject a wire as...
void parametricBounds(Vec2 &lower, Vec2 &upper) const
Add the loop's conservative extent, a B-spline's pole hull included, to a parametric bounding box.
std::vector< Curve2D > curves
WireClassification classify(const Vec2 &point, double lengthFloor) const
Classify a point against the loop with band floor lengthFloor: Boundary within the band,...
bool initialize(const std::vector< Curve2D > &inputCurves, WireRole wireRole, WireStatus &status, const ParametricMetric &metric={}, double joinTolerance=kWireJoinTolerance)
Build and validate the wire from an ordered closed list of curves, joining within joinTolerance throu...
One ray/surface intersection: the ray parameter and the outward normal; a quadric patch can give seve...
One straight line segment of a polygon wire, in a surface's parametric (u, v) domain.
One closed, oriented polygon loop in a surface's parametric domain: outer loops wind counter-clockwis...
bool initializeFromEdges(const std::vector< SurfaceEdge > &edges, WireRole wireRole, WireStatus &status, const ParametricMetric &metric={}, double joinTolerance=kWireJoinTolerance)
Build and validate the wire from an ordered edge list, joining within joinTolerance through metric,...
bool initialize(const std::vector< Vec2 > &inputVertices, WireRole wireRole, WireStatus &status, const ParametricMetric &metric={})
Build and validate the wire from an implicitly closed vertex ring; metric turns separations into 3D l...
A 2D point/vector in a surface's parametric (u, v) domain.
A 3D point/vector in the solid's local frame.
long long displaced
same position in both lists, more than tolerance apart
long long raysIdentical
the whole ordered list matched, position and sense
double worstDeltaT
max |dt| over positionally matched crossings, cm
long long missing
in the reference, absent from the candidate
long long extra
in the candidate, absent from the reference
std::vector< RayDef > rays
std::vector< Beam > beams
std::vector< double > cellArea
long long oddCrossingLists
long long parityMismatchIntervals
double insideLength
summed inside-segment length, cm (the chord integral)
BOOST_AUTO_TEST_CASE(PlanarBoxNavigationMatchesTGeoBBox)
o2::cad::harness::Point3D XRayPoint
void compare(std::string_view s1, std::string_view s2)
std::map< std::string, ID > expected
VectorOfTObjectPtrs other
BOOST_CHECK_EQUAL(triggersD.size(), triggers.size())