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();
779 checkClose(projected.
uCoord, 1.);
780 checkClose(projected.
vCoord, 0.);
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));
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.};
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));
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(),
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,
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));
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);
2194 double boundarySag = 0.)
2196 std::vector<surf::Vec3> vertices;
2197 std::vector<std::array<int, 3>> triangles;
2202 double worstSag = 0.;
2203 for (
const auto& triangle : triangles) {
2205 const surf::Vec3& second = vertices[triangle[1]];
2206 const surf::Vec3& third = vertices[triangle[2]];
2210 (second + third) * 0.5, (third +
first) * 0.5}) {
2214 BOOST_CHECK_LE(worstSag, maxSag);
2222 constexpr double phiMin = 1.6;
2223 constexpr double phiMax = 4.683;
2227 BOOST_REQUIRE(halfCylinder.
initialize({0., 0., 0.}, {0., 0., 1.}, {1., 0., 0.}, 10., -20., 20., phiMin,
2228 phiMax - phiMin,
false, paramRectWireCurves(phiMin, phiMax, -20., 20.), {},
2230 checkDisplayMeshFollowsPatch(halfCylinder, 10., 26);
2234 BOOST_REQUIRE(quarterTorus.
initialize({0., 0., 0.}, {0., 0., 1.}, {1., 0., 0.}, 10., 0.3, 0.,
surf::kHalfPi, 0.,
2237 checkDisplayMeshFollowsPatch(quarterTorus, 10.3, 14, 0.3 * (1. - std::cos(
surf::kPi / 6.)));
2241 BOOST_REQUIRE(halfCone.
initialize({0., 0., 0.}, {0., 0., 1.}, {1., 0., 0.}, 8., 12., -10., 10., phiMin,
2242 phiMax - phiMin,
false, paramRectWireCurves(phiMin, phiMax, -10., 10.), {},
2244 checkDisplayMeshFollowsPatch(halfCone, 12., 26);
2256 const double determinant =
dot(edgeA, normalB);
2257 if (std::abs(determinant) < 1.e-12) {
2261 const double firstWeight =
dot(
offset, normalB) / determinant;
2262 const double secondWeight =
dot(direction,
cross(
offset, edgeA)) / determinant;
2264 return firstWeight >= 0. && secondWeight >= 0. && firstWeight + secondWeight <= 1. &&
distance > 0.;
2269void checkDisplayMeshLeavesHolesOpen(
const surf::BoundedSurface& patch,
int holeCount,
double expectedArea,
2272 std::vector<surf::Vec3> vertices;
2273 std::vector<std::array<int, 3>> triangles;
2277 int trianglesOverHole = 0;
2278 for (
const auto& triangle : triangles) {
2280 const surf::Vec3& second = vertices[triangle[1]];
2281 const surf::Vec3& third = vertices[triangle[2]];
2283 trianglesOverHole += rayHitsTriangle(rayOrigin, rayDirection,
first, second, third);
2285 if (expectedArea >= 0.) {
2286 checkClose(area, expectedArea);
2300 BOOST_REQUIRE(square.
initialize({0., 0., 0.}, {1., 0., 0.}, {0., 1., 0.}, {{-2., -2.}, {2., -2.}, {2., 2.}, {-2., 2.}},
2301 {{{-1., -1.}, {-1., 1.}, {1., 1.}, {1., -1.}}},
error));
2302 checkDisplayMeshLeavesHolesOpen(square, 1, 12., {0.013, 0.007, 1.}, down);
2306 BOOST_REQUIRE(annulus.initialize({0., 0., 0.}, {1., 0., 0.}, {0., 1., 0.}, {surf::Curve2D::makeCircle({0., 0.}, 2.)},
2307 {{surf::Curve2D::makeCircle({0., 0.}, 1.,
true)}},
error));
2309 checkDisplayMeshLeavesHolesOpen(annulus, 1, polygonFactor * (4. - 1.), {0.013, 0.007, 1.}, down);
2313 BOOST_REQUIRE(window.initialize({0., 0., 0.}, {0., 0., 1.}, {1., 0., 0.}, 10., -20., 20., 0.,
surf::kPi,
false,
2314 paramRectWireCurves(0.,
surf::kPi, -20., 20.),
2315 {paramRectWireCurves(1., 2., -5., 5.)},
error));
2316 checkDisplayMeshLeavesHolesOpen(window, 1, -1., {0., 0., 0.013}, {std::cos(1.5), std::sin(1.5), 0.});
2324 std::vector<surf::Vec3> vertices;
2325 std::vector<std::array<int, 3>> triangles;
2328 for (
const auto& triangle : triangles) {
2329 area += 0.5 *
norm(
cross(vertices[triangle[1]] - vertices[triangle[0]], vertices[triangle[2]] - vertices[triangle[0]]));
2339 const std::vector<surf::Vec2> outer{{-2., -2.}, {2., -2.}, {2., 2.}, {-2., 2.}};
2342 BOOST_REQUIRE(plain.
initialize({0., 0., 0.}, {1., 0., 0.}, {0., 1., 0.}, outer, {},
error));
2346 BOOST_REQUIRE(touching.
initialize({0., 0., 0.}, {1., 0., 0.}, {0., 1., 0.}, outer, {{{2., 2.}, {1., 0.}, {0., 1.}}},
error));
2347 checkClose(displayArea(touching), displayArea(plain));
2351 BOOST_REQUIRE(sharing.
initialize({0., 0., 0.}, {1., 0., 0.}, {0., 1., 0.}, outer,
2352 {{{-1., -1.}, {-1., 0.}, {0., 0.}, {0., -1.}}, {{0., 0.}, {0., 1.}, {1., 1.}, {1., 0.}}},
error));
2353 checkClose(displayArea(sharing), displayArea(plain));
2357 BOOST_REQUIRE(fullCylinder.
initialize({0., 0., 0.}, {0., 0., 1.}, {1., 0., 0.}, 10., -20., 20., 0.,
surf::kTwoPi,
false,
2360 BOOST_REQUIRE(seamWindow.
initialize({0., 0., 0.}, {0., 0., 1.}, {1., 0., 0.}, 10., -20., 20., 0.,
surf::kTwoPi,
false,
2362 {paramRectWireCurves(-0.5, 0.5, -5., 5.)},
error));
2363 checkClose(displayArea(seamWindow), displayArea(fullCylinder));
2384 constexpr double kRelative = 1.e-12;
2388 BOOST_TEST_CONTEXT(what)
2391 const double exact = rectangle.capacityContribution();
2392 const double contour = wired.capacityContribution();
2393 BOOST_CHECK_GT(std::abs(exact), 1.e-6);
2394 checkClose(contour, exact, kRelative * std::abs(exact));
2401 const double phiLow = 0.3, phiHigh = 2.4, hLow = -0.8, hHigh = 1.9;
2402 BOOST_REQUIRE(rectangle.
initialize(centre, axis,
reference, 1.7, hLow, hHigh, phiLow, phiHigh - phiLow,
false,
2404 BOOST_REQUIRE(wired.
initialize(centre, axis,
reference, 1.7, hLow, hHigh, phiLow, phiHigh - phiLow,
false,
2405 paramRectWireCurves(phiLow, phiHigh, hLow, hHigh), {},
error));
2406 compare(
"cylinder", rectangle, wired);
2411 const double phiLow = -0.4, phiHigh = 1.9, hLow = 0.2, hHigh = 2.1;
2412 BOOST_REQUIRE(rectangle.
initialize(centre, axis,
reference, 1.1, 2.3, hLow, hHigh, phiLow, phiHigh - phiLow,
2414 BOOST_REQUIRE(wired.
initialize(centre, axis,
reference, 1.1, 2.3, hLow, hHigh, phiLow, phiHigh - phiLow,
false,
2415 paramRectWireCurves(phiLow, phiHigh, hLow, hHigh), {},
error));
2416 compare(
"cone", rectangle, wired);
2421 const double phiLow = 0.2, phiHigh = 2.7, thetaLow = 0.4, thetaHigh = 2.3;
2422 BOOST_REQUIRE(rectangle.
initialize(centre, axis,
reference, 2.2, thetaLow, thetaHigh, phiLow, phiHigh - phiLow,
2424 BOOST_REQUIRE(wired.
initialize(centre, axis,
reference, 2.2, thetaLow, thetaHigh, phiLow, phiHigh - phiLow,
false,
2425 paramRectWireCurves(phiLow, phiHigh, thetaLow, thetaHigh), {},
error));
2426 compare(
"sphere", rectangle, wired);
2431 const double ringLow = 0.1, ringHigh = 2.2, tubeLow = -0.3, tubeHigh = 1.8;
2432 BOOST_REQUIRE(rectangle.
initialize(centre, axis,
reference, 4., 1.4, ringLow, ringHigh - ringLow, tubeLow,
2433 tubeHigh - tubeLow,
false,
error));
2434 BOOST_REQUIRE(wired.
initialize(centre, axis,
reference, 4., 1.4, ringLow, ringHigh - ringLow, tubeLow,
2435 tubeHigh - tubeLow,
false,
2436 paramRectWireCurves(ringLow, ringHigh, tubeLow, tubeHigh), {},
error));
2437 compare(
"torus", rectangle, wired);
2447 const double radius = 1.7;
2448 const std::vector<surf::Curve2D> trim{surf::Curve2D::makeCircle({1.0, 0.2}, 0.6)};
2449 BOOST_REQUIRE(disk.
initialize(centre, axis,
reference,
radius, -2., 2., 0.,
surf::kTwoPi,
false, trim, {},
2455 std::vector<surf::CurveWire> innerWires;
2460 const auto gridRelativeError = [&](
int samples) {
2462 outerWire, innerWires,
2468 return std::abs(grid - contour) / std::abs(contour);
2470 const double at128 = gridRelativeError(128);
2471 const double at512 = gridRelativeError(512);
2472 const double at2048 = gridRelativeError(2048);
2476 BOOST_CHECK_LT(at512, 1.e-4);
2477 BOOST_CHECK_LT(at2048, 1.e-4);
2483 BOOST_CHECK_GT(at128, 1.e-3);
2484 BOOST_CHECK_GT(at2048, at512);
2494 const auto onCylinder = [](
double phi,
double height) {
2495 return Vec3{2. * std::cos(phi), 2. * std::sin(phi),
height};
2500 const std::vector<Curve2D> outer{Curve2D::makeLine({0., -3.}, {
surf::kTwoPi, -3.}),
2503 Curve2D::makeLine({0., 3.}, {0., -3.})};
2504 const std::vector<Curve2D> hole{Curve2D::makeLine({2.0, -1.}, {2.5, -1.}), Curve2D::makeLine({2.5, -1.}, {2.5, 1.}),
2505 Curve2D::makeLine({2.5, 1.}, {2.0, 1.}), Curve2D::makeLine({2.0, 1.}, {2.0, -1.})};
2506 BOOST_REQUIRE(windowed.
initialize({0., 0., 0.}, {0., 0., 1.}, {1., 0., 0.}, 2., -3., 3., 0.,
surf::kTwoPi,
false,
2507 outer, {hole},
error));
2514 std::vector<surf::RayHit> hits;
2515 windowed.
appendIntersections({0., 0., 0.}, {std::cos(2.25), std::sin(2.25), 0.}, 0., 1.e30, hits);
2518 windowed.
appendIntersections({0., 0., 0.}, {std::cos(0.5), std::sin(0.5), 0.}, 0., 1.e30, hits);
2519 BOOST_REQUIRE_EQUAL(hits.size(), 1u);
2520 checkClose(hits.front().distance, 2.);
2524 const std::vector<Curve2D> arcOuter{Curve2D::makeCircle({
surf::kPi, 0.}, 0.5)};
2525 BOOST_REQUIRE(arcTrim.
initialize({0., 0., 0.}, {0., 0., 1.}, {1., 0., 0.}, 2., -1., 1., 0.,
surf::kTwoPi,
false,
2526 arcOuter, {},
error));
2533 const auto onSphere = [](
double theta,
double phi) {
2534 return Vec3{2. * std::sin(theta) * std::cos(phi), 2. * std::sin(theta) * std::sin(phi), 2. * std::cos(theta)};
2537 BOOST_REQUIRE(sphereWire.
initialize({0., 0., 0.}, {0., 0., 1.}, {1., 0., 0.}, 2.,
surf::kHalfPi / 2.,
surf::kHalfPi,
2542 BOOST_REQUIRE(sphereScalar.
initialize({0., 0., 0.}, {0., 0., 1.}, {1., 0., 0.}, 2.,
surf::kHalfPi / 2.,
2552 const std::vector<Curve2D> wideOuter{Curve2D::makeLine({0., -1.}, {7., -1.}), Curve2D::makeLine({7., -1.}, {7., 1.}),
2553 Curve2D::makeLine({7., 1.}, {0., 1.}), Curve2D::makeLine({0., 1.}, {0., -1.})};
2554 BOOST_CHECK(!tooWide.
initialize({0., 0., 0.}, {0., 0., 1.}, {1., 0., 0.}, 2., -1., 1., 0.,
surf::kTwoPi,
false,
2555 wideOuter, {},
error));
2564template <
typename Po
intAt>
2566 double uCoord,
double vCoord,
double tolerance = 1.e-6)
2568 const double step = 1.e-5;
2569 const surf::Vec3 dU = (pointAt(uCoord + step, vCoord) - pointAt(uCoord - step, vCoord)) * (0.5 / step);
2570 const surf::Vec3 dV = (pointAt(uCoord, vCoord + step) - pointAt(uCoord, vCoord - step)) * (0.5 / step);
2575 surface.parametricMetric({uCoord, vCoord}, gUU, gUV, gVV);
2576 checkClose(gUU,
dot(dU, dU), tolerance);
2577 checkClose(gUV,
dot(dU, dV), tolerance);
2578 checkClose(gVV,
dot(dV, dV), tolerance);
2593 const Vec3 axisU{2., 0., 0.};
2594 const Vec3 axisV{1., 3., 0.};
2595 const std::vector<Vec2> unitSquare{{0., 0.}, {1., 0.}, {1., 1.}, {0., 1.}};
2597 BOOST_REQUIRE(plane.
initialize({0.5, -1., 2.}, axisU, axisV, unitSquare, {},
error));
2598 checkMetricAgainstFiniteDifference(plane, [&](
double u,
double v) {
return plane.
toGlobal({u,
v}); }, 0.3, 0.7);
2604 checkClose(gUU, 4.);
2605 checkClose(gUV, 2.);
2606 checkClose(gVV, 10.);
2614 BOOST_REQUIRE(curvedPlane.
initialize({0., 0., 0.}, {1., 0., 0.}, {0., 1., 0.},
2615 {surf::Curve2D::makeCircle({0., 0.}, 1.)}, {},
error));
2616 checkMetricAgainstFiniteDifference(
2617 curvedPlane, [&](
double u,
double v) {
return curvedPlane.toGlobal({u,
v}); }, 0.2, -0.4);
2621 for (
const double radius : {0.01, 100.}) {
2623 BOOST_REQUIRE(cylinder.
initialize({0., 0., 0.}, {0., 0., 1.}, {1., 0., 0.},
radius, -1., 1., 0.,
surf::kTwoPi,
2625 checkMetricAgainstFiniteDifference(cylinder, [&](
double u,
double v) {
return cylinder.
pointAt(u,
v); }, 1.1, 0.3, 1.e-4 *
radius *
radius);
2632 BOOST_REQUIRE(sphere.
initialize({1., 2., 3.}, {0., 0., 1.}, {1., 0., 0.}, 2.5, 0.,
surf::kPi, 0.,
surf::kTwoPi,
2634 checkMetricAgainstFiniteDifference(sphere, [&](
double u,
double v) {
return sphere.
pointAt(
v, u); }, 0.9, 1.2);
2641 checkClose(gUU, 0.);
2642 checkClose(gVV, 2.5 * 2.5);
2645 checkClose(gUU, 0.);
2648 checkClose(gUU, 2.5 * 2.5);
2654 BOOST_REQUIRE(cone.
initialize({0., 0., 0.}, {0., 0., 1.}, {1., 0., 0.}, 0., 4., 0., 2., 0.,
surf::kTwoPi,
false,
2656 checkMetricAgainstFiniteDifference(cone, [&](
double u,
double v) {
return cone.
pointAt(u,
v); }, 2.0, 1.3);
2662 checkClose(gUU, 0.);
2663 checkClose(gVV, 1. + 2. * 2.);
2666 checkClose(gUU, 16.);
2672 BOOST_REQUIRE(torus.
initialize({0., 0., 0.}, {0., 0., 1.}, {1., 0., 0.}, 5., 1.5, 0.,
surf::kTwoPi, 0.,
2674 checkMetricAgainstFiniteDifference(torus, [&](
double u,
double v) {
return torus.
pointAt(u,
v); }, 0.7, 2.1);
2680 checkClose(gUU, 6.5 * 6.5);
2681 checkClose(gVV, 1.5 * 1.5);
2683 checkClose(gUU, 3.5 * 3.5);
2697 const auto trimWithPhiDrift = [](
double drift) {
2698 return std::vector<Curve2D>{Curve2D::makeLine({0.2, -1.}, {1.2, -1.}), Curve2D::makeLine({1.2, -1.}, {1.2, 1.}),
2699 Curve2D::makeLine({1.2, 1.}, {0.2 + drift, 1.}),
2700 Curve2D::makeLine({0.2 + drift, 1.}, {0.2 + drift, -1.})};
2702 const auto acceptsDrift = [&](
double radius,
double drift) {
2704 return cylinder.
initialize({0., 0., 0.}, {0., 0., 1.}, {1., 0., 0.},
radius, -1., 1., 0.,
surf::kTwoPi,
false,
2705 trimWithPhiDrift(drift), {},
error);
2734 const auto polygonAcceptsGap = [](
double gap) {
2735 const std::vector<SurfaceEdge>
edges{{{0., 0.}, {1., 0.}},
2736 {{1., 0.}, {1., 1.}},
2737 {{1., 1.}, {0., 1.}},
2738 {{0., 1.}, {gap, 0.}}};
2740 WireStatus status = WireStatus::Valid;
2743 const auto curveAcceptsGap = [](
double gap) {
2744 const std::vector<surf::Curve2D> curves{surf::Curve2D::makeLine({0., 0.}, {1., 0.}),
2745 surf::Curve2D::makeLine({1., 0.}, {1., 1.}),
2746 surf::Curve2D::makeLine({1., 1.}, {0., 1.}),
2747 surf::Curve2D::makeLine({0., 1.}, {gap, 0.})};
2749 WireStatus status = WireStatus::Valid;
2750 return wire.
initialize(curves, WireRole::Outer, status, {});
2755 for (
const double gap : {0., 1.e-8, 1.e-7}) {
2760 for (
const double gap : {1.e-5, 1.e-4}) {
2773 const char*
raw =
reinterpret_cast<const char*
>(&
value);
2777void appendDoubles(std::vector<char>&
bytes, std::initializer_list<double>
values)
2780 const char*
raw =
reinterpret_cast<const char*
>(&
value);
2789 double modelTolerance = 0.,
uint32_t nModelEdges = 0)
2791 bytes.insert(
bytes.end(), {
'O',
'2',
'S',
'S'});
2793 appendU32(
bytes, nSurfaces);
2794 appendU32(
bytes, 0);
2796 appendDoubles(
bytes, {modelTolerance});
2799 appendU32(
bytes, nModelEdges);
2804void appendPlaneRecord(std::vector<char>&
bytes,
const FaceFrame& frame)
2806 appendU32(
bytes, 1);
2807 appendU32(
bytes, 0);
2808 appendU32(
bytes, 9);
2809 appendDoubles(
bytes, {frame.origin[0], frame.origin[1], frame.origin[2], frame.axisU[0], frame.axisU[1],
2810 frame.axisU[2], frame.axisV[0], frame.axisV[1], frame.axisV[2]});
2811 appendU32(
bytes, 1);
2812 appendU32(
bytes, 0);
2813 appendU32(
bytes, 4);
2814 const double extentU = frame.extentU;
2815 const double extentV = frame.extentV;
2816 const std::array<std::array<double, 4>, 4>
edges{{{0., 0., extentU, 0.},
2817 {extentU, 0., extentU, extentV},
2818 {extentU, extentV, 0., extentV},
2819 {0., extentV, 0., 0.}}};
2820 for (
const auto& edge :
edges) {
2821 appendU32(
bytes, 0);
2822 appendU32(
bytes, 4);
2823 appendDoubles(
bytes, {edge[0], edge[1], edge[2], edge[3]});
2830 const Point3D& axisV,
double radius,
double holeRadius = 0.)
2832 appendU32(
bytes, 1);
2833 appendU32(
bytes, 0);
2834 appendU32(
bytes, 9);
2835 appendDoubles(
bytes, {
center[0],
center[1],
center[2], axisU[0], axisU[1], axisU[2], axisV[0], axisV[1], axisV[2]});
2836 const uint32_t nWires = holeRadius > 0. ? 2u : 1u;
2837 appendU32(
bytes, nWires);
2838 appendU32(
bytes, 0);
2839 appendU32(
bytes, 1);
2840 appendU32(
bytes, 1);
2841 appendU32(
bytes, 5);
2843 if (holeRadius > 0.) {
2844 appendU32(
bytes, 1);
2845 appendU32(
bytes, 1);
2846 appendU32(
bytes, 1);
2847 appendU32(
bytes, 5);
2852std::filesystem::path writeSidecarFile(
const std::string&
name,
const std::vector<char>&
bytes)
2854 const auto path = std::filesystem::temp_directory_path() /
name;
2855 std::ofstream out(
path, std::ios::binary);
2856 out.write(
bytes.data(),
static_cast<std::streamsize
>(
bytes.size()));
2857 BOOST_REQUIRE(out.good());
2865 constexpr double halfX = 1.;
2866 constexpr double halfY = 2.;
2867 constexpr double halfZ = 3.;
2869 std::vector<char> boxBytes;
2870 appendSidecarHeader(boxBytes, 6);
2872 appendPlaneRecord(boxBytes, boxFaceFrame(
faceIndex, halfX, halfY, halfZ));
2874 const auto boxPath = writeSidecarFile(
"o2_sidecar_roundtrip_box.bin", boxBytes);
2876 SurfaceSolid
box(
"sidecarBox");
2878 std::filesystem::remove(boxPath);
2884 TGeoBBox referenceBox(
"referenceBox", halfX, halfY, halfZ);
2885 compareContainsGrid(
box, referenceBox, 4., 7);
2886 compareDistance(
box, referenceBox, {5., 0.5, 0.5}, {-1., 0., 0.});
2887 compareDistance(
box, referenceBox, {0., 0., 0.}, unitDirection(1., 1., 1.));
2888 checkClose(
box.Capacity(), referenceBox.Capacity(), 1.e-9);
2891 constexpr double radius = 2.;
2892 constexpr double halfHeight = 3.;
2894 std::vector<char> tubeBytes;
2895 appendSidecarHeader(tubeBytes, 3);
2896 appendU32(tubeBytes, 2);
2897 appendU32(tubeBytes, 0);
2898 appendU32(tubeBytes, 14);
2899 appendDoubles(tubeBytes, {0., 0., 0., 0., 0., 1., 1., 0., 0.,
radius, -halfHeight, halfHeight, 0., 2. *
surf::kPi});
2900 appendU32(tubeBytes, 0);
2902 appendDiskPlaneRecord(tubeBytes, {0., 0., halfHeight}, {1., 0., 0.}, {0., 1., 0.},
radius);
2903 appendDiskPlaneRecord(tubeBytes, {0., 0., -halfHeight}, {1., 0., 0.}, {0., -1., 0.},
radius);
2904 const auto tubePath = writeSidecarFile(
"o2_sidecar_roundtrip_tube.bin", tubeBytes);
2906 SurfaceSolid tube(
"sidecarTube");
2908 std::filesystem::remove(tubePath);
2914 TGeoTube referenceTube(
"referenceTube", 0.,
radius, halfHeight);
2915 compareContainsGrid(tube, referenceTube, 4., 7);
2916 compareDistance(tube, referenceTube, {5., 0.5, 1.}, {-1., 0., 0.});
2917 compareDistance(tube, referenceTube, {0., 0., 0.}, unitDirection(1., 1., 1.));
2918 checkClose(tube.Capacity(), referenceTube.Capacity(), 1.e-9);
2921 const auto badPath = writeSidecarFile(
"o2_sidecar_bad_magic.bin", {
'X',
'X',
'X',
'X', 0, 0, 0, 0});
2922 SurfaceSolid bad(
"sidecarBad");
2924 std::filesystem::remove(badPath);
2929 std::vector<char> truncatedBytes;
2930 appendSidecarHeader(truncatedBytes, 1);
2931 const auto truncatedPath = writeSidecarFile(
"o2_sidecar_truncated.bin", truncatedBytes);
2932 SurfaceSolid truncated(
"sidecarTruncated");
2934 std::filesystem::remove(truncatedPath);
2942 constexpr double halfX = 1.;
2943 constexpr double halfY = 2.;
2944 constexpr double halfZ = 3.;
2945 const auto boxBytesWithHeader = [&](
uint32_t version,
double modelTolerance) {
2946 std::vector<char>
bytes;
2947 appendSidecarHeader(
bytes, 6,
version, modelTolerance);
2949 appendPlaneRecord(
bytes, boxFaceFrame(
faceIndex, halfX, halfY, halfZ));
2953 const auto loadFrom = [](
const char*
name,
const std::vector<char>&
bytes, SurfaceSolid& solid) {
2956 std::filesystem::remove(
path);
2961 SurfaceSolid
v2(
"sidecarV2");
2962 BOOST_REQUIRE(loadFrom(
"o2_sidecar_v2.bin", boxBytesWithHeader(2, 3.5e-5),
v2));
2964 checkClose(
v2.GetModelTolerance(), 3.5e-5, 1.e-18);
2966 checkClose(
v2.GetModelTolerance(), 3.5e-5, 1.e-18);
2969 SurfaceSolid v2Silent(
"sidecarV2Silent");
2970 BOOST_REQUIRE(loadFrom(
"o2_sidecar_v2_silent.bin", boxBytesWithHeader(2, 0.), v2Silent));
2974 SurfaceSolid
v1(
"sidecarV1");
2975 BOOST_REQUIRE(loadFrom(
"o2_sidecar_v1.bin", boxBytesWithHeader(1, 0.),
v1));
2977 checkClose(
v1.GetModelTolerance(), 1.e-6, 1.e-18);
2980 SurfaceSolid handBuilt(
"handBuilt");
2982 handBuilt.SetModelTolerance(1.e-4);
2983 checkClose(handBuilt.GetModelTolerance(), 1.e-4, 1.e-18);
2984 handBuilt.SetModelTolerance(-1.);
2985 checkClose(handBuilt.GetModelTolerance(), 1.e-4, 1.e-18);
2991 SurfaceSolid v4(
"sidecarV4");
2992 BOOST_CHECK(!loadFrom(
"o2_sidecar_v4.bin", boxBytesWithHeader(4, 1.e-5), v4));
2996 std::vector<char> stump;
2997 stump.insert(stump.end(), {
'O',
'2',
'S',
'S'});
2998 appendU32(stump, 2);
2999 appendU32(stump, 6);
3000 appendU32(stump, 0);
3001 SurfaceSolid stumped(
"sidecarV2Stump");
3002 BOOST_CHECK(!loadFrom(
"o2_sidecar_v2_stump.bin", stump, stumped));
3010 constexpr double radius = 2.;
3011 constexpr double halfHeight = 3.;
3013 std::vector<char>
bytes;
3014 appendSidecarHeader(
bytes, 3);
3015 appendU32(
bytes, 2);
3016 appendU32(
bytes, 0);
3017 appendU32(
bytes, 14);
3018 appendDoubles(
bytes, {0., 0., 0., 0., 0., 1., 1., 0., 0.,
radius, -halfHeight, halfHeight, 0., 2. *
surf::kPi});
3019 appendU32(
bytes, 1);
3020 appendU32(
bytes, 0);
3021 appendU32(
bytes, 4);
3022 const std::array<std::array<double, 4>, 4>
edges{{{0., -halfHeight, 2. *
surf::kPi, -halfHeight},
3024 {2. *
surf::kPi, halfHeight, 0., halfHeight},
3025 {0., halfHeight, 0., -halfHeight}}};
3026 for (
const auto& edge :
edges) {
3027 appendU32(
bytes, 0);
3028 appendU32(
bytes, 4);
3029 appendDoubles(
bytes, {edge[0], edge[1], edge[2], edge[3]});
3031 appendDiskPlaneRecord(
bytes, {0., 0., halfHeight}, {1., 0., 0.}, {0., 1., 0.},
radius);
3032 appendDiskPlaneRecord(
bytes, {0., 0., -halfHeight}, {1., 0., 0.}, {0., -1., 0.},
radius);
3033 const auto path = writeSidecarFile(
"o2_sidecar_wiretrim_cylinder.bin",
bytes);
3035 SurfaceSolid solid(
"sidecarWireCylinder");
3037 std::filesystem::remove(
path);
3044 compareContainsGrid(solid,
reference, 4., 7);
3045 compareDistance(solid,
reference, {5., 0.5, 1.}, {-1., 0., 0.});
3046 compareDistance(solid,
reference, {0., 0., 0.}, unitDirection(1., 1., 1.));
3047 checkClose(solid.Capacity(),
reference.Capacity(), 1.e-6);
3054 constexpr double majorR = 3.;
3055 constexpr double minorR = 1.;
3057 std::vector<char>
bytes;
3058 appendSidecarHeader(
bytes, 1);
3059 appendU32(
bytes, 5);
3060 appendU32(
bytes, 0);
3061 appendU32(
bytes, 15);
3062 appendDoubles(
bytes, {0., 0., 0., 0., 0., 1., 1., 0., 0., majorR, minorR, 0., 2. *
surf::kPi, 0., 2. *
surf::kPi});
3063 appendU32(
bytes, 0);
3064 const auto path = writeSidecarFile(
"o2_sidecar_roundtrip_torus.bin",
bytes);
3066 SurfaceSolid solid(
"sidecarTorus");
3068 std::filesystem::remove(
path);
3074 TGeoTorus
reference(
"sidecarTorusRef", majorR, 0., minorR);
3075 compareContainsGrid(solid,
reference, 4.5, 9);
3076 checkClose(solid.Capacity(),
reference.Capacity(), 1.e-7);
3088 constexpr double radius = 2.;
3089 constexpr double halfHeight = 3.;
3092 constexpr double joinGap = 5.e-6;
3094 const auto cylinderBytes = [&](
uint32_t version,
double modelTolerance) {
3095 std::vector<char>
bytes;
3096 appendSidecarHeader(
bytes, 3,
version, modelTolerance);
3097 appendU32(
bytes, 2);
3098 appendU32(
bytes, 0);
3099 appendU32(
bytes, 14);
3100 appendDoubles(
bytes, {0., 0., 0., 0., 0., 1., 1., 0., 0.,
radius, -halfHeight, halfHeight, 0., 2. *
surf::kPi});
3101 appendU32(
bytes, 1);
3102 appendU32(
bytes, 0);
3103 appendU32(
bytes, 5);
3108 const std::array<std::array<double, 4>, 5>
edges{
3109 {{0., -halfHeight,
surf::kPi, -halfHeight},
3112 {2. *
surf::kPi, halfHeight, 0., halfHeight},
3113 {0., halfHeight, 0., -halfHeight}}};
3114 for (
const auto& edge :
edges) {
3115 appendU32(
bytes, 0);
3116 appendU32(
bytes, 4);
3117 appendDoubles(
bytes, {edge[0], edge[1], edge[2], edge[3]});
3119 appendDiskPlaneRecord(
bytes, {0., 0., halfHeight}, {1., 0., 0.}, {0., 1., 0.},
radius);
3120 appendDiskPlaneRecord(
bytes, {0., 0., -halfHeight}, {1., 0., 0.}, {0., -1., 0.},
radius);
3123 const auto loadFrom = [](
const char*
name,
const std::vector<char>&
bytes, SurfaceSolid& solid) {
3126 std::filesystem::remove(
path);
3132 SurfaceSolid declared(
"sidecarJoinDeclared");
3133 BOOST_REQUIRE(loadFrom(
"o2_sidecar_join_declared.bin", cylinderBytes(2, 1.e-4), declared));
3135 declared.CloseShape();
3139 compareContainsGrid(declared,
reference, 4., 7);
3140 compareDistance(declared,
reference, {5., 0.5, 1.}, {-1., 0., 0.});
3143 SurfaceSolid silent(
"sidecarJoinSilent");
3144 BOOST_CHECK(!loadFrom(
"o2_sidecar_join_silent.bin", cylinderBytes(1, 0.), silent));
3148 SurfaceSolid tight(
"sidecarJoinTight");
3149 BOOST_CHECK(!loadFrom(
"o2_sidecar_join_tight.bin", cylinderBytes(2, 2.e-6), tight));
3158void appendBSplineEdgePlaneRecord(std::vector<char>&
bytes,
const FaceFrame& frame)
3160 appendU32(
bytes, 1);
3161 appendU32(
bytes, 0);
3162 appendU32(
bytes, 9);
3163 appendDoubles(
bytes, {frame.origin[0], frame.origin[1], frame.origin[2], frame.axisU[0], frame.axisU[1],
3164 frame.axisU[2], frame.axisV[0], frame.axisV[1], frame.axisV[2]});
3165 appendU32(
bytes, 1);
3166 appendU32(
bytes, 0);
3167 appendU32(
bytes, 4);
3168 const double extentU = frame.extentU;
3169 const double extentV = frame.extentV;
3171 appendU32(
bytes, 2);
3172 appendU32(
bytes, 22);
3173 appendDoubles(
bytes, {3., 4.,
3174 0., 0., extentU / 3., 0., 2. * extentU / 3., 0., extentU, 0.,
3176 0., 0., 0., 0., 1., 1., 1., 1.});
3177 const std::array<std::array<double, 4>, 3> lines{
3178 {{extentU, 0., extentU, extentV}, {extentU, extentV, 0., extentV}, {0., extentV, 0., 0.}}};
3179 for (
const auto& edge : lines) {
3180 appendU32(
bytes, 0);
3181 appendU32(
bytes, 4);
3182 appendDoubles(
bytes, {edge[0], edge[1], edge[2], edge[3]});
3194 const double w = std::sqrt(0.5);
3195 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}};
3196 return surf::Curve2D::makeBSpline(2, poles, {1.,
w, 1.,
w, 1.,
w, 1.,
w, 1.},
3197 {0., 0., 0., 1., 1., 2., 2., 3., 3., 4., 4., 4.});
3200surf::Curve2D quarterCircleBSpline(
double cu,
double cv,
double r,
double a0)
3203 const double aMid = 0.5 * (a0 + a1);
3204 const std::vector<surf::Vec2> poles{{cu +
r * std::cos(a0), cv +
r * std::sin(a0)},
3205 {cu +
r * std::sqrt(2.) * std::cos(aMid), cv +
r * std::sqrt(2.) * std::sin(aMid)},
3206 {cu +
r * std::cos(a1), cv +
r * std::sin(a1)}};
3207 return surf::Curve2D::makeBSpline(2, poles, {1., std::sqrt(0.5), 1.}, {0., 0., 0., 1., 1., 1.});
3223 WireStatus status = WireStatus::Valid;
3224 BOOST_REQUIRE(exactWire.
initialize({surf::Curve2D::makeCircle({0., 0.}, 1.)}, WireRole::Outer, status));
3229 BOOST_REQUIRE(splineWire.
initialize({fullCircleBSpline(0., 0., 1.)}, WireRole::Outer, status));
3236 BOOST_CHECK(splineWire.
classify({1. - justInsideTheBand, 0.}) == WireClassification::Boundary);
3237 BOOST_CHECK(splineWire.
classify({1. + justInsideTheBand, 0.}) == WireClassification::Boundary);
3242 BOOST_CHECK(exactWire.classify({1. - justInsideTheBand, 0.}) == WireClassification::Inside);
3247 BOOST_REQUIRE(bigCylinder.
initialize({0., 0., 0.}, {0., 0., 1.}, {1., 0., 0.}, 100., -1., 1., 0.,
surf::kTwoPi,
3251 BOOST_REQUIRE(squareWire.
initialize({surf::Curve2D::makeLine({0., -1.}, {1., -1.}),
3252 surf::Curve2D::makeLine({1., -1.}, {1., 1.}),
3253 surf::Curve2D::makeLine({1., 1.}, {0., 1.}),
3254 surf::Curve2D::makeLine({0., 1.}, {0., -1.})},
3255 WireRole::Outer, status, bigMetric));
3257 BOOST_CHECK(squareWire.classify({0.5, 1. - 1.e-10}, bigMetric) == WireClassification::Inside);
3259 BOOST_CHECK(squareWire.classify({0.5, 1. - 1.e-12}, bigMetric) == WireClassification::Boundary);
3264 const double seamDrift = 4.e-7;
3269 second.
poles.front() = {second.
poles.front().uCoord + seamDrift, second.
poles.front().vCoord};
3271 BOOST_REQUIRE(driftedWire.
initialize({first, second, third, fourth}, WireRole::Outer, status));
3272 for (
const auto& curve : driftedWire.curves) {
3277 const Vec2 thisEnd = driftedWire.
curves[
index].loopEnd();
3278 const Vec2 nextStart = driftedWire.
curves[(
index + 1) % driftedWire.
curves.size()].loopStart();
3283 for (
const auto& curve : driftedWire.curves) {
3284 const auto& polyline = curve.bsplineSamples();
3285 BOOST_REQUIRE(polyline.size() >= 2);
3321 const double trimRadius = 0.5;
3323 BOOST_REQUIRE(splineTrimmed.
initialize({0., 0., 0.}, {0., 0., 1.}, {1., 0., 0.}, 2., -1., 1., 0.,
surf::kTwoPi,
3325 {quarterCircleBSpline(centrePhi, 0., trimRadius, 0.),
3326 quarterCircleBSpline(centrePhi, 0., trimRadius, surf::kHalfPi),
3327 quarterCircleBSpline(centrePhi, 0., trimRadius, surf::kPi),
3328 quarterCircleBSpline(centrePhi, 0., trimRadius, 3. * surf::kHalfPi)},
3333 BOOST_REQUIRE(arcTrimmed.
initialize({0., 0., 0.}, {0., 0., 1.}, {1., 0., 0.}, 2., -1., 1., 0.,
surf::kTwoPi,
false,
3334 {Curve2D::makeCircle({centrePhi, 0.}, trimRadius)}, {},
error));
3338 std::vector<surf::RayHit> hits;
3339 surface.appendIntersections({0., 0., 0.}, {std::cos(phi), std::sin(phi), 0.}, 0., 1.e30, hits);
3346 const double justInside = 0.5 * band;
3350 const auto hits = hitAt(splineTrimmed, centrePhi);
3351 BOOST_REQUIRE_EQUAL(hits.size(), 1u);
3356 const auto hits = hitAt(splineTrimmed, centrePhi + trimRadius - justInside);
3357 BOOST_REQUIRE_EQUAL(hits.size(), 1u);
3363 const auto hits = hitAt(splineTrimmed, centrePhi + trimRadius + justInside);
3364 BOOST_REQUIRE_EQUAL(hits.size(), 1u);
3369 BOOST_CHECK(hitAt(splineTrimmed, centrePhi + trimRadius + 100. * band).
empty());
3373 const auto inside = hitAt(arcTrimmed, centrePhi + trimRadius - justInside);
3374 BOOST_REQUIRE_EQUAL(inside.size(), 1u);
3376 BOOST_CHECK(hitAt(arcTrimmed, centrePhi + trimRadius + justInside).
empty());
3381 BOOST_REQUIRE(plain.
initialize({0., 0., 0.}, {0., 0., 1.}, {1., 0., 0.}, 2., -1., 1., 0.,
surf::kTwoPi,
false,
3383 const auto hits = hitAt(plain, centrePhi);
3384 BOOST_REQUIRE_EQUAL(hits.size(), 1u);
3393 constexpr double halfX = 1.;
3394 constexpr double halfY = 2.;
3395 constexpr double halfZ = 3.;
3397 std::vector<char>
bytes;
3398 appendSidecarHeader(
bytes, 6);
3399 appendBSplineEdgePlaneRecord(
bytes, boxFaceFrame(0, halfX, halfY, halfZ));
3401 appendPlaneRecord(
bytes, boxFaceFrame(
faceIndex, halfX, halfY, halfZ));
3403 const auto path = writeSidecarFile(
"o2_sidecar_bspline_box.bin",
bytes);
3405 SurfaceSolid
box(
"sidecarBSplineBox");
3407 std::filesystem::remove(
path);
3415 compareDistance(
box,
reference, {5., 0.5, 0.5}, {-1., 0., 0.});
3416 compareDistance(
box,
reference, {0., 0., 0.}, unitDirection(1., 1., 1.));
3417 checkClose(
box.Capacity(),
reference.Capacity(), 1.e-6);
3426 const auto onCylinder = [](
double phi,
double height) {
3427 return Vec3{2. * std::cos(phi), 2. * std::sin(phi),
height};
3434 const double trimRadius = 0.5;
3436 const std::vector<Curve2D> bsplineOuter{quarterCircleBSpline(centrePhi, 0., trimRadius, 0.),
3437 quarterCircleBSpline(centrePhi, 0., trimRadius,
surf::kHalfPi),
3438 quarterCircleBSpline(centrePhi, 0., trimRadius,
surf::kPi),
3439 quarterCircleBSpline(centrePhi, 0., trimRadius, 3. *
surf::kHalfPi)};
3440 BOOST_REQUIRE(bsplineDisk.
initialize({0., 0., 0.}, {0., 0., 1.}, {1., 0., 0.}, 2., -1., 1., 0.,
surf::kTwoPi,
false,
3441 bsplineOuter, {},
error));
3445 BOOST_REQUIRE(arcDisk.
initialize({0., 0., 0.}, {0., 0., 1.}, {1., 0., 0.}, 2., -1., 1., 0.,
surf::kTwoPi,
false,
3446 {Curve2D::makeCircle({centrePhi, 0.}, trimRadius)}, {},
error));
3452 for (
int phiStep = -12; phiStep <= 12; ++phiStep) {
3453 const double phi = centrePhi + 0.09 * phiStep;
3454 for (
int hStep = -12; hStep <= 12; ++hStep) {
3455 const double height = 0.09 * hStep;
3457 const double distToCentre = std::hypot(phi - centrePhi,
height);
3458 if (std::abs(distToCentre - trimRadius) < 5.e-3) {
3461 const Vec3 point = onCylinder(phi,
height);
3462 BOOST_CHECK_EQUAL(bsplineDisk.containsPointOnSurface(point), arcDisk.containsPointOnSurface(point));
3466 BOOST_CHECK_GT(compared, 100);
3470 std::vector<surf::RayHit> hits;
3471 bsplineDisk.appendIntersections({0., 0., 0.}, {std::cos(centrePhi), std::sin(centrePhi), 0.}, 0., 1.e30, hits);
3472 BOOST_REQUIRE_EQUAL(hits.size(), 1u);
3473 checkClose(hits.front().distance, 2.);
3475 bsplineDisk.appendIntersections({0., 0., 0.}, {std::cos(0.), std::sin(0.), 0.}, 0., 1.e30, hits);
3502 const auto onCylinder = [](
double phi,
double height) {
3503 return Vec3{2. * std::cos(phi), 2. * std::sin(phi),
height};
3507 const std::vector<Curve2D> outer{Curve2D::makeLine({0., -3.}, {
surf::kTwoPi, -3.}),
3510 Curve2D::makeLine({0., 3.}, {0., -3.})};
3514 const double trimRadius = 0.5;
3515 const std::vector<Curve2D> bsplineHole{quarterCircleBSpline(centrePhi, 0., trimRadius, 0.),
3516 quarterCircleBSpline(centrePhi, 0., trimRadius,
surf::kHalfPi),
3517 quarterCircleBSpline(centrePhi, 0., trimRadius,
surf::kPi),
3518 quarterCircleBSpline(centrePhi, 0., trimRadius, 3. *
surf::kHalfPi)};
3521 BOOST_REQUIRE(bsplineHoled.
initialize({0., 0., 0.}, {0., 0., 1.}, {1., 0., 0.}, 2., -3., 3., 0.,
surf::kTwoPi,
false,
3522 outer, {bsplineHole},
error));
3524 BOOST_REQUIRE(arcHoled.
initialize({0., 0., 0.}, {0., 0., 1.}, {1., 0., 0.}, 2., -3., 3., 0.,
surf::kTwoPi,
false,
3525 outer, {{Curve2D::makeCircle({centrePhi, 0.}, trimRadius)}},
error));
3528 BOOST_CHECK(!arcHoled.containsPointOnSurface(onCylinder(centrePhi, 0.)));
3529 BOOST_CHECK(!bsplineHoled.containsPointOnSurface(onCylinder(centrePhi, 0.)));
3531 BOOST_CHECK(bsplineHoled.containsPointOnSurface(onCylinder(0.5, 0.)));
3532 BOOST_CHECK(bsplineHoled.containsPointOnSurface(onCylinder(centrePhi, 2.5)));
3536 for (
int phiStep = -12; phiStep <= 12; ++phiStep) {
3537 const double phi = centrePhi + 0.09 * phiStep;
3538 for (
int hStep = -12; hStep <= 12; ++hStep) {
3539 const double height = 0.09 * hStep;
3540 if (std::abs(std::hypot(phi - centrePhi,
height) - trimRadius) < 5.e-3) {
3543 const Vec3 point = onCylinder(phi,
height);
3544 BOOST_CHECK_EQUAL(bsplineHoled.containsPointOnSurface(point), arcHoled.containsPointOnSurface(point));
3548 BOOST_CHECK_GT(compared, 100);
3551 std::vector<surf::RayHit> hits;
3552 bsplineHoled.appendIntersections({0., 0., 0.}, {std::cos(centrePhi), std::sin(centrePhi), 0.}, 0., 1.e30, hits);
3555 bsplineHoled.appendIntersections({0., 0., 0.}, {std::cos(0.5), std::sin(0.5), 0.}, 0., 1.e30, hits);
3562 BOOST_REQUIRE(singleEdgeHoled.
initialize({0., 0., 0.}, {0., 0., 1.}, {1., 0., 0.}, 2., -3., 3., 0.,
surf::kTwoPi,
3563 false, outer, {{fullCircleBSpline(centrePhi, 0., trimRadius)}},
error));
3564 BOOST_CHECK(!singleEdgeHoled.containsPointOnSurface(onCylinder(centrePhi, 0.)));
3565 BOOST_CHECK(singleEdgeHoled.containsPointOnSurface(onCylinder(0.5, 0.)));
3566 BOOST_CHECK(singleEdgeHoled.containsPointOnSurface(onCylinder(centrePhi, 2.5)));
3567 for (
int phiStep = -12; phiStep <= 12; ++phiStep) {
3568 const double phi = centrePhi + 0.09 * phiStep;
3569 for (
int hStep = -12; hStep <= 12; ++hStep) {
3570 const double height = 0.09 * hStep;
3571 if (std::abs(std::hypot(phi - centrePhi,
height) - trimRadius) < 5.e-3) {
3574 const Vec3 point = onCylinder(phi,
height);
3575 BOOST_CHECK_EQUAL(singleEdgeHoled.containsPointOnSurface(point), arcHoled.containsPointOnSurface(point));
3583BoundaryCurve quarterCircleBoundaryCurve(
double cu,
double cv,
double r,
double a0)
3586 const double aMid = 0.5 * (a0 + a1);
3587 const std::vector<Point2D> poles{{cu +
r * std::cos(a0), cv +
r * std::sin(a0)},
3588 {cu +
r * std::sqrt(2.) * std::cos(aMid), cv +
r * std::sqrt(2.) * std::sin(aMid)},
3589 {cu +
r * std::cos(a1), cv +
r * std::sin(a1)}};
3590 return BoundaryCurve::makeBSpline(2, poles, {1., std::sqrt(0.5), 1.}, {0., 0., 0., 1., 1., 1.});
3598void checkSolidsIdentical(
const SurfaceSolid& solid,
const SurfaceSolid&
other,
double extent,
int samples)
3606 static_cast<int>(solid.GetNavigationReliability()));
3614 for (
int dimension = 0; dimension < 3; ++dimension) {
3619 for (
const auto& point : probeGrid(extent,
samples)) {
3620 BOOST_TEST_CONTEXT(
"point = (" << point[0] <<
", " << point[1] <<
", " << point[2] <<
")")
3624 solid.Safety(point.data(), solid.Contains(point.data())));
3625 for (
const auto& direction : probeDirections()) {
3627 solid.DistFromOutside(point.data(),
direction.data(), 3));
3629 solid.DistFromInside(point.data(),
direction.data(), 3));
3636std::unique_ptr<SurfaceSolid> writeAndReadBack(
const SurfaceSolid& solid)
3638 const auto path = std::filesystem::temp_directory_path() /
3639 (std::string(
"o2_bvhsurfacesolid_persist_") + solid.GetName() +
".root");
3641 TFile
file(
path.string().c_str(),
"RECREATE");
3642 BOOST_REQUIRE(!
file.IsZombie());
3644 file.WriteObject(
const_cast<SurfaceSolid*
>(&solid),
"solid");
3646 std::unique_ptr<SurfaceSolid> restored;
3648 TFile
file(
path.string().c_str(),
"READ");
3649 BOOST_REQUIRE(!
file.IsZombie());
3650 restored.reset(
file.Get<SurfaceSolid>(
"solid"));
3652 std::filesystem::remove(
path);
3671 const auto box = makeBoxSolid(
"persistBox", 1., 2., 3.);
3672 const auto tube = makeTubeSolid(
"persistTube", 1., 2., 3.);
3673 const auto cone = makeConeSolid(
"persistCone", 2., 1., 3.);
3674 const auto sphere = makeSphereSolid(
"persistSphere", 2.);
3675 const auto torus = makeTorusSolid(
"persistTorus", 3., 1.);
3676 const auto capsule = makeCapsuleSolid(
"persistCapsule", 2., 3.);
3678 for (
const auto* solid : {
box.get(), tube.get(), cone.get(), sphere.get(), torus.get(), capsule.get()}) {
3679 BOOST_TEST_CONTEXT(
"solid = " << solid->GetName())
3681 const auto restored = writeAndReadBack(*solid);
3682 BOOST_REQUIRE(restored !=
nullptr);
3683 checkSolidsIdentical(*solid, *restored, 4.5, 5);
3689 SurfaceSolid trimmed(
"persistTrimmed");
3690 constexpr double radius = 2.;
3691 constexpr double halfHeight = 3.;
3692 const std::vector<BoundaryCurve> window{quarterCircleBoundaryCurve(
surf::kPi, 0., 0.5, 0.),
3696 BOOST_REQUIRE(trimmed.AddCylindricalSurface({0., 0., 0.}, {0., 0., 1.}, {1., 0., 0.},
radius, -halfHeight,
3698 BOOST_REQUIRE(addDiskSurface(trimmed, {0., 0., halfHeight}, {1., 0., 0.}, {0., 1., 0.},
radius));
3699 BOOST_REQUIRE(addDiskSurface(trimmed, {0., 0., -halfHeight}, {1., 0., 0.}, {0., -1., 0.},
radius));
3700 trimmed.CloseShape(
false);
3703 const auto restoredTrimmed = writeAndReadBack(trimmed);
3704 BOOST_REQUIRE(restoredTrimmed !=
nullptr);
3705 checkSolidsIdentical(trimmed, *restoredTrimmed, 4.5, 5);
3709 SurfaceSolid toleranced(
"persistTolerance");
3711 BOOST_REQUIRE(addBoxFace(toleranced,
faceIndex, 1., 2., 3.));
3713 toleranced.SetModelTolerance(7.25e-5);
3714 toleranced.CloseShape(
false);
3715 const auto restoredToleranced = writeAndReadBack(toleranced);
3716 BOOST_REQUIRE(restoredToleranced !=
nullptr);
3717 checkClose(restoredToleranced->GetModelTolerance(), 7.25e-5, 1.e-18);
3721 SurfaceSolid openBox(
"persistOpenBox");
3723 BOOST_REQUIRE(addBoxFace(openBox,
faceIndex, 1., 2., 3.));
3725 openBox.CloseShape(
false);
3726 BOOST_REQUIRE(!openBox.IsNavigable());
3728 static_cast<int>(SurfaceSolid::NavigationReliability::OpenSurfaceSet));
3730 const auto restoredOpenBox = writeAndReadBack(openBox);
3731 BOOST_REQUIRE(restoredOpenBox !=
nullptr);
3733 checkSolidsIdentical(openBox, *restoredOpenBox, 4.5, 5);
3742 SurfaceSolid
empty(
"emptySolid");
3744 static_cast<int>(SurfaceSolid::NavigationReliability::Undetermined));
3748 empty.CloseShape(
false);
3751 static_cast<int>(SurfaceSolid::NavigationReliability::Undetermined));
3761std::vector<std::array<double, 3>> spiralDirections(
int count)
3763 std::vector<std::array<double, 3>> directions;
3764 directions.reserve(
count);
3766 const double cosTheta = 1. - 2. * (
index + 0.5) /
count;
3767 const double sinTheta = std::sqrt(1. - cosTheta * cosTheta);
3768 const double phi = 2.399963229728653 *
index;
3769 directions.push_back({
sinTheta * std::cos(phi),
sinTheta * std::sin(phi), cosTheta});
3785 const auto box = makeBoxSolid(
"dirBox", 1., 2., 3.);
3786 const auto tube = makeTubeSolid(
"dirTube", 1., 2., 3.);
3787 const auto cone = makeConeSolid(
"dirCone", 2., 1., 3.);
3788 const auto sphere = makeSphereSolid(
"dirSphere", 2.);
3789 const auto torus = makeTorusSolid(
"dirTorus", 3., 1.);
3790 const auto capsule = makeCapsuleSolid(
"dirCapsule", 2., 3.);
3792 const auto directions = spiralDirections(13);
3793 for (
const auto* solid : {
box.get(), tube.get(), cone.get(), sphere.get(), torus.get(), capsule.get()}) {
3794 BOOST_TEST_CONTEXT(
"solid = " << solid->GetName())
3796 BOOST_REQUIRE(solid->IsNavigable());
3797 for (
const auto& point : probeGrid(4.5, 7)) {
3798 const bool reference = solid->Contains(point.data());
3799 for (
const auto& direction : directions) {
3800 BOOST_TEST_CONTEXT(
"point = (" << point[0] <<
", " << point[1] <<
", " << point[2] <<
") direction = ("
3801 << direction[0] <<
", " << direction[1] <<
", " << direction[2] <<
")")
3824 constexpr double halfX = 1.;
3825 constexpr double halfY = 2.;
3826 constexpr double halfZ = 3.;
3827 constexpr double gap = 0.05;
3829 SurfaceSolid gapped(
"gappedBox");
3831 BOOST_REQUIRE(addBoxFace(gapped,
faceIndex, halfX, halfY, halfZ));
3835 BOOST_REQUIRE(gapped.AddPlanarSurface({halfX, -halfY, -halfZ}, {0., 1., 0.}, {0., 0., 1.},
3836 rectangleWire(2. * halfY, halfZ - gap)));
3837 BOOST_REQUIRE(gapped.AddPlanarSurface({halfX, -halfY, gap}, {0., 1., 0.}, {0., 0., 1.},
3838 rectangleWire(2. * halfY, halfZ - gap)));
3839 gapped.CloseShape(
false);
3842 BOOST_REQUIRE(!gapped.IsNavigable());
3844 static_cast<int>(SurfaceSolid::NavigationReliability::OpenSurfaceSet));
3847 const std::array<double, 3> insidePoint{0., 0.3, 0.};
3848 const std::array<double, 3> throughGap{1., 0., 0.};
3849 BOOST_CHECK(!gapped.ContainsAlongDirection(insidePoint.data(), throughGap.data()));
3851 BOOST_CHECK(gapped.Contains_Loop(insidePoint.data()));
3855 const auto intact = makeBoxSolid(
"intactBox", halfX, halfY, halfZ);
3856 BOOST_CHECK(intact->Contains(insidePoint.data()));
3857 BOOST_CHECK(intact->ContainsAlongDirection(insidePoint.data(), throughGap.data()));
3861 for (
const auto& point : probeGrid(4.5, 7)) {
3862 BOOST_TEST_CONTEXT(
"point = (" << point[0] <<
", " << point[1] <<
", " << point[2] <<
")")
3864 BOOST_CHECK_EQUAL(gapped.Contains(point.data()), gapped.Contains_Loop(point.data()));
3877std::unique_ptr<SurfaceSolid> makeLPrismSolid(
const char*
name,
double height = 1.)
3880 const std::vector<Point2D> footprint{{0., 0.}, {3., 0.}, {3., 1.}, {1., 1.}, {1., 2.}, {0., 2.}};
3882 auto solid = std::make_unique<SurfaceSolid>(
name);
3885 std::vector<Point2D> bottomWire;
3886 bottomWire.reserve(footprint.size());
3887 for (
const auto&
vertex : footprint) {
3890 BOOST_REQUIRE(solid->AddPlanarSurface({0., 0., 0.}, {0., 1., 0.}, {1., 0., 0.}, bottomWire));
3891 BOOST_REQUIRE(solid->AddPlanarSurface({0., 0., height}, {1., 0., 0.}, {0., 1., 0.}, footprint));
3897 const auto&
end = footprint[(
index + 1) % footprint.size()];
3898 const double deltaU =
end[0] -
start[0];
3899 const double deltaV =
end[1] -
start[1];
3900 const double length = std::hypot(deltaU, deltaV);
3901 BOOST_REQUIRE(solid->AddPlanarSurface({start[0], start[1], 0.}, {deltaU / length, deltaV / length, 0.},
3904 solid->CloseShape();
3915 SurfaceSolid
box(
"preCloseBox");
3916 addBoxSurfaces(
box, 1., 2., 3.);
3917 BOOST_REQUIRE(!
box.HasBVH());
3919 for (
const auto& point : probeGrid(4.5, 5)) {
3920 BOOST_TEST_CONTEXT(
"point = (" << point[0] <<
", " << point[1] <<
", " << point[2] <<
")")
3922 const bool inside = std::abs(point[0]) < 1. && std::abs(point[1]) < 2. && std::abs(point[2]) < 3.;
3935 const auto box = makeBoxSolid(
"boundaryPolicyBox", 1., 2., 3.);
3936 const std::array<double, 3> onFace{1., 0.5, 0.5};
3937 const std::array<double, 3> outward{1., 0., 0.};
3938 const std::array<double, 3> inward{-1., 0., 0.};
3945 checkClose(
box->DistFromInside(onFace.data(), outward.data(), 3),
3946 reference.DistFromInside(onFace.data(), outward.data(), 3));
3947 checkClose(
box->DistFromOutside(onFace.data(), inward.data(), 3),
3948 reference.DistFromOutside(onFace.data(), inward.data(), 3));
3951 checkClose(
box->DistFromInside(onFace.data(), inward.data(), 3), 2.);
3954 for (
const auto& direction : {outward, inward}) {
3955 checkDistanceAgainstLoop(*
box, onFace, direction);
3971 const double invSqrt2 = 1. / std::sqrt(2.);
3974 const auto box = makeBoxSolid(
"grazeBox", 1., 2., 3.);
3975 const std::array<double, 3> grazeDirection{invSqrt2, -invSqrt2, 0.};
3976 const std::array<double, 3> grazeOrigin{1. - 5. * invSqrt2, 2. + 5. * invSqrt2, 0.};
3977 BOOST_REQUIRE(!
box->Contains(grazeOrigin.data()));
3980 const std::array<double, 3> pastTouch{1. + 1.e-3 * invSqrt2, 2. - 1.e-3 * invSqrt2, 0.};
3981 BOOST_REQUIRE(!
box->Contains(pastTouch.data()));
3983 BOOST_CHECK_EQUAL(
box->DistFromOutside(grazeOrigin.data(), grazeDirection.data(), 3), TGeoShape::Big());
3984 BOOST_CHECK_EQUAL(
box->DistFromOutside_Loop(grazeOrigin.data(), grazeDirection.data()), TGeoShape::Big());
3987 const auto prism = makeLPrismSolid(
"grazePrism");
3988 BOOST_REQUIRE(prism->IsNavigable());
3989 checkClose(prism->Capacity(), 4.);
3993 const std::array<double, 3> reflexDirection{invSqrt2, -invSqrt2, 0.};
3994 const std::array<double, 3> reflexOrigin{1. - 0.5 * invSqrt2, 1. + 0.5 * invSqrt2, 0.5};
3995 BOOST_REQUIRE(prism->Contains(reflexOrigin.data()));
3996 const std::array<double, 3> pastReflex{1. + 1.e-3 * invSqrt2, 1. - 1.e-3 * invSqrt2, 0.5};
3997 BOOST_REQUIRE(prism->Contains(pastReflex.data()));
3999 const double touchDistance = 0.5;
4000 const double exitDistance = 0.5 + std::sqrt(2.);
4001 const double reported = prism->DistFromInside(reflexOrigin.data(), reflexDirection.data(), 3);
4002 BOOST_CHECK_GT(reported, touchDistance + 1.e-6);
4003 checkClose(reported, exitDistance);
4004 BOOST_CHECK_EQUAL(prism->DistFromInside_Loop(reflexOrigin.data(), reflexDirection.data()), reported);
4010 const auto box = makeBoxSolid(
"describeBox", 1., 2., 3.);
4011 const SurfaceSolid::Point3D inside{0.2, 0.3, 0.4};
4012 const double invSqrt2 = 1. / std::sqrt(2.);
4014 const std::vector<std::pair<SurfaceSolid::Point3D, double>> cases{
4015 {{1., 0., 0.}, 0.8}, {{invSqrt2, -invSqrt2, 0.}, 0.8 * std::sqrt(2.)}, {{0., 0., -2.}, 3.4}};
4016 for (
const auto& [direction, exitDistance] : cases) {
4017 BOOST_TEST_CONTEXT(
"direction = (" << direction[0] <<
", " << direction[1] <<
", " << direction[2] <<
")")
4019 std::vector<SurfaceSolid::ContainsCrossing> bvhCrossings;
4020 std::vector<SurfaceSolid::ContainsCrossing> loopCrossings;
4021 box->DescribeContainsCrossings(inside, direction, bvhCrossings, loopCrossings);
4023 BOOST_REQUIRE_EQUAL(bvhCrossings.size(), 1u);
4024 BOOST_REQUIRE_EQUAL(loopCrossings.size(), 1u);
4026 checkClose(bvhCrossings[0].
distance, exitDistance);
4027 BOOST_CHECK_GT(bvhCrossings[0].normalAlignment, 0.);
4029 BOOST_CHECK(
box->ContainsAlongDirection(inside.data(), direction.data()));
4039 const auto prism = makeLPrismSolid(
"sweepPrism");
4040 BOOST_REQUIRE(prism->IsNavigable());
4042 sweepDistanceAgainstLoop(*prism, 3.5, 5);
4044 const auto directions = spiralDirections(13);
4045 for (
const auto& point : probeGrid(3.5, 7)) {
4046 const bool inside = prism->Contains(point.data());
4047 BOOST_TEST_CONTEXT(
"point = (" << point[0] <<
", " << point[1] <<
", " << point[2] <<
")")
4050 for (
const auto& direction : directions) {
4051 BOOST_CHECK_EQUAL(prism->ContainsAlongDirection(point.data(), direction.data()), inside);
4054 const bool expected = point[2] > 0. && point[2] < 1. &&
4055 ((point[0] > 0. && point[0] < 3. && point[1] > 0. && point[1] < 1.) ||
4056 (point[0] > 0. && point[0] < 1. && point[1] >= 1. && point[1] < 2.));
4076 const std::vector<Vec2> poles{{0., 0.}, {1., 2.}, {3., 2.}, {4., 0.}};
4077 const Curve2D clamped = Curve2D::makeBSpline(3, poles, {}, {0., 0., 0., 0., 1., 1., 1., 1.});
4078 const Curve2D uniform = Curve2D::makeBSpline(3, poles, {}, {0., 1., 2., 3., 4., 5., 6., 7.});
4080 BOOST_REQUIRE(clamped.valid());
4081 BOOST_REQUIRE(uniform.valid());
4084 for (
const auto* curve : {&clamped, &uniform}) {
4085 const Vec2
start = curve->startPoint();
4086 const Vec2
end = curve->endPoint();
4087 const Vec2 evaluatedStart = curve->pointAt(0.);
4088 const Vec2 evaluatedEnd = curve->pointAt(1.);
4089 checkClose(
start.uCoord, evaluatedStart.uCoord);
4090 checkClose(
start.vCoord, evaluatedStart.vCoord);
4091 checkClose(
end.uCoord, evaluatedEnd.uCoord);
4092 checkClose(
end.vCoord, evaluatedEnd.vCoord);
4097 checkClose(clamped.startPoint().uCoord, 0.);
4098 checkClose(clamped.endPoint().uCoord, 4.);
4099 BOOST_CHECK_GT(std::hypot(uniform.startPoint().uCoord - poles.front().uCoord,
4100 uniform.startPoint().vCoord - poles.front().vCoord),
4105 const Vec2 uniformStart = uniform.startPoint();
4106 const Vec2 uniformEnd = uniform.endPoint();
4109 BOOST_CHECK(wire.
initialize({uniform, Curve2D::makeLine(uniformEnd, uniformStart)}, surf::WireRole::Outer, status));
4127 const double sweep = 350. *
surf::kPi / 180.;
4128 const double overshoot = 0.4;
4129 const std::vector<Curve2D> outer{
4130 Curve2D::makeBSpline(2, {{0., -1.}, {-overshoot, 0.}, {0.5 * sweep, 1.}}, {}, {0., 0., 0., 1., 1., 1.}),
4131 Curve2D::makeBSpline(2, {{0.5 * sweep, 1.}, {sweep + overshoot, 0.}, {sweep, -1.}}, {}, {0., 0., 0., 1., 1., 1.}),
4132 Curve2D::makeLine({sweep, -1.}, {0., -1.})};
4135 Vec2 hullLower{1.e300, 1.e300};
4136 Vec2 hullUpper{-1.e300, -1.e300};
4139 BOOST_REQUIRE(hullWire.
initialize(outer, surf::WireRole::Outer, hullStatus));
4141 BOOST_REQUIRE_GT(hullUpper.uCoord - hullLower.uCoord,
surf::kTwoPi);
4144 Vec2 tightLower{1.e300, 1.e300};
4145 Vec2 tightUpper{-1.e300, -1.e300};
4147 BOOST_CHECK_LT(tightUpper.uCoord - tightLower.uCoord,
surf::kTwoPi);
4151 BOOST_CHECK_MESSAGE(
surface.initialize({0., 0., 0.}, {0., 0., 1.}, {1., 0., 0.}, 2., -2., 2., 0.,
surf::kTwoPi,
4152 false, outer, {},
error),
4153 "near-full-turn trim rejected: " <<
error);
4157 const std::vector<Curve2D> overWrapped{Curve2D::makeLine({0., -1.}, {
surf::kTwoPi + 0.5, -1.}),
4160 Curve2D::makeLine({0., 1.}, {0., -1.})};
4161 BOOST_CHECK(!tooWide.
initialize({0., 0., 0.}, {0., 0., 1.}, {1., 0., 0.}, 2., -2., 2., 0.,
surf::kTwoPi,
false,
4162 overWrapped, {},
error));
4173 SurfaceSolid solid(
"rejectingSolid");
4174 BOOST_REQUIRE(addBoxFace(solid, 0, 1., 2., 3.));
4175 BOOST_REQUIRE_EQUAL(solid.GetNsurfaces(), 1);
4179 BOOST_CHECK(!solid.AddPlanarSurface({0., 0., 0.}, {1., 0., 0.}, {1., 0., 0.}, rectangleWire(1., 1.)));
4180 BOOST_CHECK(!solid.AddPlanarSurface({0., 0., 0.}, {1., 0., 0.}, {0., 1., 0.}, {{0., 0.}, {1., 0.}}));
4181 BOOST_CHECK(!solid.AddCylindricalSurface({0., 0., 0.}, {0., 0., 1.}, {1., 0., 0.}, 0., -1., 1.));
4186 std::vector<char>
bytes;
4187 appendSidecarHeader(
bytes, 1);
4188 appendU32(
bytes, 2);
4189 appendU32(
bytes, 0);
4190 appendU32(
bytes, 14);
4191 appendDoubles(
bytes, {0., 0., 0., 0., 0., 1., 1., 0., 0., 0. , -1., 1., 0., 2. *
surf::kPi});
4192 appendU32(
bytes, 0);
4193 const auto path = writeSidecarFile(
"o2_sidecar_rejected_face.bin",
bytes);
4194 SurfaceSolid loaded(
"loadedRejecting");
4196 std::filesystem::remove(
path);
4211 constexpr double halfX = 1.;
4212 constexpr double halfY = 1.5;
4213 constexpr double halfZ = 2.;
4216 SurfaceSolid closedBox(
"rimClosedBox");
4217 addBoxSurfaces(closedBox, halfX, halfY, halfZ);
4218 closedBox.CloseShape(
false);
4219 BOOST_REQUIRE(closedBox.IsNavigable());
4223 BOOST_CHECK_SMALL(closedBox.GetMaxRimIsolation(), 1.e-12);
4224 BOOST_CHECK_SMALL(closedBox.GetUnmatchedRimLength(), 1.e-12);
4226 BOOST_CHECK_SMALL(closedBox.GetRimChordResolution(), 1.e-12);
4228 BOOST_CHECK_CLOSE(closedBox.GetTotalRimLength(), 16. * (halfX + halfY + halfZ), 1.e-9);
4233 constexpr double delta = 1.e-3;
4234 SurfaceSolid shiftedBox(
"rimShiftedBox");
4237 BOOST_REQUIRE(addBoxFace(shiftedBox,
faceIndex, halfX, halfY, halfZ,
false,
center));
4239 shiftedBox.CloseShape(
false);
4240 BOOST_CHECK_CLOSE(shiftedBox.GetMaxRimIsolation(), delta, 1.e-6);
4242 BOOST_CHECK(shiftedBox.GetMaxRimIsolation() > shiftedBox.GetRimChordResolution());
4252 constexpr double halfX = 1.;
4253 constexpr double halfY = 1.5;
4254 constexpr double halfZ = 2.;
4256 SurfaceSolid resampled(
"rimResampledBox");
4258 BOOST_REQUIRE(addBoxFace(resampled,
faceIndex, halfX, halfY, halfZ));
4261 const FaceFrame frame = boxFaceFrame(5, halfX, halfY, halfZ);
4262 const double extentU = frame.extentU;
4263 const double extentV = frame.extentV;
4264 BOOST_REQUIRE(resampled.AddPlanarSurface(frame.origin, frame.axisU, frame.axisV,
4266 {0.5 * extentU, 0.},
4268 {extentU, 0.5 * extentV},
4270 {0.5 * extentU, extentV},
4272 {0., 0.5 * extentV}}));
4273 resampled.CloseShape(
false);
4278 BOOST_CHECK(resampled.GetBoundaryEdgeCount() > 0);
4286 BOOST_CHECK_SMALL(resampled.GetMaxRimIsolation(), 1.e-12);
4287 BOOST_CHECK_SMALL(resampled.GetUnmatchedRimLength(), 1.e-12);
4295 constexpr double halfX = 1.;
4296 constexpr double halfY = 1.5;
4297 constexpr double halfZ = 2.;
4298 constexpr double delta = 1.e-3;
4300 SurfaceSolid shiftedBox(
"rimReportBox");
4303 BOOST_REQUIRE(addBoxFace(shiftedBox,
faceIndex, halfX, halfY, halfZ,
false,
center));
4305 shiftedBox.CloseShape(
false);
4307 const auto& rims = shiftedBox.GetRimReports();
4308 BOOST_REQUIRE_EQUAL(
static_cast<int>(rims.size()), shiftedBox.GetRimCount());
4310 int boundaryRims = 0;
4311 auto worst = SurfaceSolid::NavigationReliability::Reliable;
4312 double openLength = 0.;
4313 for (
const auto& rim : rims) {
4314 BOOST_CHECK(rim.surface >= 0 && rim.surface < shiftedBox.GetNsurfaces());
4318 openLength += rim.unmatchedLength;
4319 if (rim.state == SurfaceSolid::NavigationReliability::OpenSurfaceSet) {
4326 worst = std::max(worst, rim.state);
4329 BOOST_CHECK(worst == shiftedBox.GetNavigationReliability());
4330 BOOST_CHECK_CLOSE(openLength, shiftedBox.GetUnmatchedRimLength(), 1.e-9);
4333 const auto lifted = std::find_if(rims.begin(), rims.end(), [](
const auto& rim) { return rim.surface == 4; });
4334 BOOST_REQUIRE(lifted != rims.end());
4335 BOOST_CHECK(lifted->state == SurfaceSolid::NavigationReliability::OpenSurfaceSet);
4336 BOOST_CHECK_CLOSE(lifted->maxIsolation, delta, 1.e-6);
4337 BOOST_CHECK(lifted->maxIsolationFace >= 0 && lifted->maxIsolationFace != 4);
4339 BOOST_CHECK_CLOSE(lifted->maxIsolationPoint[2], halfZ + delta, 1.e-6);
4346 constexpr double radius = 1.;
4347 constexpr double halfHeight = 2.;
4348 SurfaceSolid openTube(
"rimOpenTube");
4350 openTube.AddCylindricalSurface({0., 0., 0.}, {0., 0., 1.}, {1., 0., 0.},
radius, -halfHeight, halfHeight));
4351 openTube.CloseShape(
false);
4361 BOOST_CHECK_CLOSE(openTube.GetUnmatchedRimLength(), openTube.GetTotalRimLength(), 1.e-9);
4368 BOOST_CHECK_CLOSE(openTube.GetRimChordResolution(), exactSagitta, 1.);
4380 SurfaceSolid unstated(
"rimToleranceUnstated");
4381 addBoxSurfaces(unstated, 1., 1., 1.);
4382 unstated.CloseShape(
false);
4386 SurfaceSolid stated(
"rimToleranceStated");
4387 addBoxSurfaces(stated, 1., 1., 1.);
4388 stated.SetModelTolerance(2.5e-7);
4389 stated.CloseShape(
false);
4390 BOOST_CHECK_CLOSE(stated.GetRimMatchTolerance(), 2.5e-7, 1.e-9);
4405 BOOST_CHECK_CLOSE(circle.
uVariation(0., 1.), 8., 1.e-9);
4408 std::vector<surf::Vec2> poles{{0., 0.}, {1., 1.}, {2., 0.}, {1., -1.}, {0., 0.}};
4409 std::vector<double> knots{0., 0., 0., 0.25, 0.5, 0.75, 1., 1., 1.};
4410 const surf::Curve2D loop = surf::Curve2D::makeBSpline(2, poles, {}, knots);
4412 BOOST_CHECK_GT(loop.
uVariation(0., 1.), 0.5);
4429 const double radius = 1.2;
4430 const double height = 1.5;
4431 const double holeRadius = 0.3;
4433 const double centreV = 0.5 *
height;
4435 std::vector<o2::cad::O2BVHSurfaceSolid::PlanarBoundaryCurve> outer{
4443 const double corner = std::sqrt(2.) / 2.;
4444 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.}}};
4445 std::vector<o2::cad::O2BVHSurfaceSolid::Point2D> poles;
4446 for (
const auto& pole : unit) {
4447 poles.push_back({centreU + holeRadius * pole[0], centreV + holeRadius * pole[1]});
4449 const std::vector<double>
weights{1., corner, 1., corner, 1., corner, 1., corner, 1.};
4450 const std::vector<double> knots{0., 0., 0., 0.25, 0.25, 0.5, 0.5, 0.75, 0.75, 1., 1., 1.};
4451 std::vector<std::vector<o2::cad::O2BVHSurfaceSolid::PlanarBoundaryCurve>> holes{
4455 BOOST_REQUIRE(withHole.
AddCylindricalSurface({0., 0., 0.}, {0., 0., 1.}, {1., 0., 0.},
radius, 0.,
height, 0.,
4457 std::vector<double> contributions;
4459 BOOST_REQUIRE_EQUAL(contributions.size(), 1u);
4464 BOOST_CHECK_CLOSE(contributions[0],
radius *
radius * chartArea / 3., 1.e-6);
4473 std::vector<surf::Vec2> poles{{0., 0.}, {1., 2.}, {2., 2.}, {3., 0.}};
4474 std::vector<double> knots{0., 1., 2., 3., 4., 5., 6., 7.};
4475 const surf::Curve2D unclamped = surf::Curve2D::makeBSpline(3, poles, {}, knots);
4479 BOOST_CHECK_GT(std::abs(
start.uCoord - poles.front().uCoord) + std::abs(
start.vCoord - poles.front().vCoord), 1.e-3);
4480 BOOST_CHECK_GT(std::abs(
end.uCoord - poles.back().uCoord) + std::abs(
end.vCoord - poles.back().vCoord), 1.e-3);
4483 BOOST_CHECK_SMALL(std::abs(unclamped.
pointAt(1.).
vCoord -
end.vCoord), 1.e-12);
4502std::vector<std::pair<std::vector<unsigned int>, std::vector<unsigned char>>>
4503 boxEdgeIdentity(
double halfX,
double halfY,
double halfZ)
4505 using Key = std::array<long long, 6>;
4506 auto quantize = [](
double value) {
return static_cast<long long>(std::llround(
value * 1.e9)); };
4507 std::map<Key, unsigned int> edgeIds;
4508 std::vector<std::pair<std::vector<unsigned int>, std::vector<unsigned char>>> perFace(6);
4511 const FaceFrame frame = boxFaceFrame(
faceIndex, halfX, halfY, halfZ);
4512 const auto corners = rectangleWire(frame.extentU, frame.extentV);
4514 const auto& startUV = corners[
segment];
4515 const auto& endUV = corners[(
segment + 1) % corners.size()];
4516 auto toGlobal = [&](
const Point2D& uv) {
4517 return Point3D{frame.origin[0] + frame.axisU[0] * uv[0] + frame.axisV[0] * uv[1],
4518 frame.origin[1] + frame.axisU[1] * uv[0] + frame.axisV[1] * uv[1],
4519 frame.origin[2] + frame.axisU[2] * uv[0] + frame.axisV[2] * uv[1]};
4524 quantize(
end[0]), quantize(
end[1]), quantize(
end[2])};
4525 const Key backwardKey{forwardKey[3], forwardKey[4], forwardKey[5],
4526 forwardKey[0], forwardKey[1], forwardKey[2]};
4528 const bool reversed = backwardKey < forwardKey;
4529 const Key canonical = reversed ? backwardKey : forwardKey;
4530 const auto inserted = edgeIds.emplace(canonical,
static_cast<unsigned int>(edgeIds.size()));
4531 perFace[
faceIndex].first.push_back(inserted.first->second);
4533 static_cast<unsigned char>(SurfaceSolid::kEdgeAnchored | (reversed ? SurfaceSolid::kEdgeReversed : 0u)));
4536 BOOST_REQUIRE_EQUAL(edgeIds.size(), 12u);
4541std::unique_ptr<SurfaceSolid> makeIdentifiedBox(
const char*
name,
double halfX,
double halfY,
double halfZ)
4543 auto solid = std::make_unique<SurfaceSolid>(
name);
4544 addBoxSurfaces(*solid, halfX, halfY, halfZ);
4545 const auto identity = boxEdgeIdentity(halfX, halfY, halfZ);
4558 const auto box = makeIdentifiedBox(
"identityBox", 1., 2., 3.);
4559 box->CloseShape(
false);
4575 BOOST_CHECK_SMALL(
box->GetMaxSharedEdgeDeviation(), 1.e-15);
4583 SurfaceSolid openBox(
"identityOpenBox");
4585 BOOST_REQUIRE(addBoxFace(openBox,
faceIndex, 1., 2., 3.));
4587 const auto identity = boxEdgeIdentity(1., 2., 3.);
4591 openBox.CloseShape(
false);
4600 static_cast<int>(SurfaceSolid::NavigationReliability::OpenSurfaceSet));
4608 const auto box = makeIdentifiedBox(
"identityReversed", 1., 2., 3.);
4610 auto identity = boxEdgeIdentity(1., 2., 3.);
4611 for (
auto& flag : identity[0].second) {
4612 flag =
static_cast<unsigned char>(flag ^ SurfaceSolid::kEdgeReversed);
4614 BOOST_REQUIRE(
box->SetSurfaceBoundaryEdges(0, identity[0].first, identity[0].second));
4615 box->CloseShape(
false);
4620 static_cast<int>(SurfaceSolid::NavigationReliability::ReversedFaces));
4624 const auto box = makeIdentifiedBox(
"identityNonManifold", 1., 2., 3.);
4625 BOOST_REQUIRE(addBoxFace(*
box, 0, 1., 2., 3.));
4626 const auto identity = boxEdgeIdentity(1., 2., 3.);
4627 BOOST_REQUIRE(
box->SetSurfaceBoundaryEdges(6, identity[0].first, identity[0].second));
4628 box->CloseShape(
false);
4632 static_cast<int>(SurfaceSolid::NavigationReliability::NonManifold));
4643 constexpr double delta = 3.e-4;
4644 SurfaceSolid shifted(
"identityShifted");
4648 BOOST_REQUIRE(addBoxFace(shifted,
faceIndex, 1., 2., 3.,
false, centre));
4650 const auto identity = boxEdgeIdentity(1., 2., 3.);
4654 shifted.CloseShape(
false);
4660 checkClose(shifted.GetMaxSharedEdgeDeviation(), delta, 1.e-12);
4671 SurfaceSolid flipped(
"identityFlippedWinding");
4672 const auto identity = boxEdgeIdentity(1., 2., 3.);
4673 std::vector<std::vector<unsigned int>>
ids(6);
4674 std::vector<std::vector<unsigned char>>
flags(6);
4676 FaceFrame frame = boxFaceFrame(
faceIndex, 1., 2., 3.);
4677 auto corners = rectangleWire(frame.extentU, frame.extentV);
4680 std::reverse(corners.begin(), corners.end());
4681 BOOST_REQUIRE(flipped.AddPlanarSurface(frame.origin, frame.axisU, frame.axisV, corners));
4682 const size_t n = corners.size();
4685 for (
size_t j = 0;
j <
n; ++
j) {
4692 flipped.CloseShape(
false);
4698 BOOST_CHECK_SMALL(flipped.GetMaxSharedEdgeDeviation(), 1.e-12);
4706 SurfaceSolid partial(
"identityPartial");
4707 addBoxSurfaces(partial, 1., 2., 3.);
4708 const auto identity = boxEdgeIdentity(1., 2., 3.);
4712 partial.CloseShape(
false);
4718 BOOST_CHECK(!partial.SetSurfaceBoundaryEdges(6, identity[0].first, identity[0].second));
4719 BOOST_CHECK(!partial.SetSurfaceBoundaryEdges(0, {1u, 2u}, {0}));
4726 const auto box = makeIdentifiedBox(
"identityPersist", 1., 2., 3.);
4727 box->CloseShape(
false);
4728 BOOST_REQUIRE(
box->HasEdgeIdentity());
4730 const auto restored = writeAndReadBack(*
box);
4731 BOOST_REQUIRE(restored !=
nullptr);
4736 checkClose(restored->GetMaxSharedEdgeDeviation(),
box->GetMaxSharedEdgeDeviation(), 1.e-18);
4737 checkSolidsIdentical(*
box, *restored, 4.5, 5);
4745 constexpr double halfX = 1.;
4746 constexpr double halfY = 2.;
4747 constexpr double halfZ = 3.;
4748 const auto identity = boxEdgeIdentity(halfX, halfY, halfZ);
4751 std::vector<char>
bytes;
4752 appendSidecarHeader(
bytes, 6,
version, 1.e-7, writeIdentity ? 12u : 0u);
4754 appendPlaneRecord(
bytes, boxFaceFrame(
faceIndex, halfX, halfY, halfZ));
4759 if (writeIdentity) {
4760 for (
size_t e = 0;
e <
ids.size(); ++
e) {
4770 SurfaceSolid
v3(
"sidecarV3Identity");
4771 const auto v3Path = writeSidecarFile(
"o2_sidecar_v3_identity.bin", boxBytes(3,
true));
4773 std::filesystem::remove(v3Path);
4774 v3.CloseShape(
false);
4779 BOOST_CHECK_SMALL(
v3.GetMaxSharedEdgeDeviation(), 1.e-15);
4783 checkClose(
v3.Capacity(),
reference.Capacity(), 1.e-9);
4787 SurfaceSolid v3Silent(
"sidecarV3Silent");
4788 const auto silentPath = writeSidecarFile(
"o2_sidecar_v3_silent.bin", boxBytes(3,
false));
4790 std::filesystem::remove(silentPath);
4791 v3Silent.CloseShape(
false);
4796 SurfaceSolid
v2(
"sidecarV2StillLoads");
4797 const auto v2Path = writeSidecarFile(
"o2_sidecar_v2_still_loads.bin", boxBytes(2,
false));
4799 std::filesystem::remove(v2Path);
4800 v2.CloseShape(
false);
4803 checkClose(
v2.GetModelTolerance(), 1.e-7, 1.e-18);
4807 std::vector<char> mislabelled;
4808 appendSidecarHeader(mislabelled, 6, 2, 1.e-7);
4811 appendPlaneRecord(mislabelled, boxFaceFrame(
faceIndex, halfX, halfY, halfZ));
4813 SurfaceSolid mislabelledSolid(
"sidecarV3Mislabelled");
4814 const auto badPath = writeSidecarFile(
"o2_sidecar_v3_mislabelled.bin", mislabelled);
4816 std::filesystem::remove(badPath);
4827 const auto box = makeIdentifiedBox(
"identitySampling", 1., 2., 3.);
4828 box->CloseShape(
false);
4829 const double deviation =
box->GetMaxSharedEdgeDeviation();
4830 const bool navigable =
box->IsNavigable();
4836 SurfaceSolid resampled(
"identitySamplingResampled");
4837 const auto identity = boxEdgeIdentity(1., 2., 3.);
4839 const FaceFrame frame = boxFaceFrame(
faceIndex, 1., 2., 3.);
4840 auto corners = rectangleWire(frame.extentU, frame.extentV);
4841 std::vector<Point2D> dense;
4842 std::vector<unsigned int>
ids;
4843 std::vector<unsigned char>
flags;
4846 const auto&
b = corners[(
segment + 1) % corners.size()];
4848 dense.push_back({0.5 * (
a[0] +
b[0]), 0.5 * (
a[1] +
b[1])});
4856 BOOST_REQUIRE(resampled.AddPlanarSurface(frame.origin, frame.axisU, frame.axisV, dense));
4859 resampled.CloseShape(
false);
4870 checkClose(
box->GetMaxSharedEdgeDeviation(), deviation, 1.e-18);
4886std::array<double, 5> torusRayQuartic(
double majorRadius,
double minorRadius,
4889 const double dirDotDir = dir[0] * dir[0] + dir[1] * dir[1] + dir[2] * dir[2];
4890 const double originDotDir =
origin[0] * dir[0] +
origin[1] * dir[1] +
origin[2] * dir[2];
4891 const double originDotOrigin =
4893 const double constantK = majorRadius * majorRadius - minorRadius * minorRadius;
4894 const double transverseE = dir[0] * dir[0] + dir[1] * dir[1];
4895 const double transverseF =
origin[0] * dir[0] +
origin[1] * dir[1];
4897 const double fourRSquared = 4. * majorRadius * majorRadius;
4898 return {dirDotDir * dirDotDir,
4899 4. * dirDotDir * originDotDir,
4900 4. * originDotDir * originDotDir + 2. * dirDotDir * (originDotOrigin + constantK) -
4901 fourRSquared * transverseE,
4902 4. * originDotDir * (originDotOrigin + constantK) - 2. * fourRSquared * transverseF,
4903 (originDotOrigin + constantK) * (originDotOrigin + constantK) -
4904 fourRSquared * transverseG};
4909std::vector<double> torusRootsAtScale(
double scale)
4911 const Point3D origin{2.094269422822338 * scale, 3.292530879918199 * scale,
4912 1.9347519602583996 * scale};
4913 const Point3D dir{-0.7547297076674779, -0.03154700875395883, -0.655276929704412};
4914 const auto c = torusRayQuartic(2.5 * scale, 0.8 * scale,
origin, dir);
4916 return {roots.begin(), roots.end()};
4925 const auto reference = torusRootsAtScale(1.);
4926 BOOST_REQUIRE_EQUAL(
reference.size(), 2u);
4927 for (
const double scale : {0.5, 2., 10.}) {
4928 const auto roots = torusRootsAtScale(scale);
4929 BOOST_REQUIRE_EQUAL(roots.size(),
reference.size());
4930 for (
size_t i = 0;
i < roots.size(); ++
i) {
4931 checkClose(roots[
i],
reference[
i] * scale, 1.e-12);
4948 const auto reference = torusRootsAtScale(1.);
4949 BOOST_REQUIRE_EQUAL(
reference.size(), 2u);
4950 for (
const double factor : {0.15, 0.12, 0.05, 0.01}) {
4951 const auto roots = torusRootsAtScale(
factor);
4952 BOOST_REQUIRE_EQUAL(roots.size(),
reference.size());
4953 for (
size_t i = 0;
i < roots.size(); ++
i) {
4961 const double scale = 0.12;
4962 const Point3D origin{2.094269422822338 * scale, 3.292530879918199 * scale,
4963 1.9347519602583996 * scale};
4964 const Point3D dir{-0.7547297076674779, -0.03154700875395883, -0.655276929704412};
4965 const auto c = torusRayQuartic(2.5 * scale, 0.8 * scale,
origin, dir);
4966 const auto evaluate = [&
c](
double t) {
4967 return (((
c[0] *
t +
c[1]) *
t +
c[2]) *
t +
c[3]) *
t +
c[4];
4970 const double t = root * scale;
4971 const double span = 0.02 *
t;
4972 BOOST_CHECK_LT(evaluate(
t - span) * evaluate(
t + span), 0.);
4993 const auto dir = std::filesystem::temp_directory_path();
4997 double origin[3] = {1.0, 1.5, 2.0};
5003 auto* cube =
new TGeoBBox(
"cube", 2.0, 2.0, 2.0);
5004 auto* drill =
new TGeoTube(
"drill", 0.0, 0.8, 2.5);
5005 TGeoCompositeShape composite(
"shape",
new TGeoSubtraction(cube, drill,
nullptr,
nullptr));
5007 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}};
5008 const std::vector<Point3D> directions{{1., 0., 0.}, {0., 1., 0.}, {0., 0., 1.}, {-1., 0., 0.}, {0.6, 0.8, 0.}};
5010 for (
const TGeoShape* original : {
static_cast<const TGeoShape*
>(&
box),
5011 static_cast<const TGeoShape*
>(&composite)}) {
5012 const std::string
path = (dir / (std::string(
"o2_shape_sidecar_") + original->ClassName() +
".root")).string();
5014 BOOST_REQUIRE_MESSAGE(harness::saveShapeToRootFile(
path, *original, &
error),
error);
5016 std::unique_ptr<TGeoShape> loaded(harness::loadShapeFromRootFile(
path, &
error));
5017 BOOST_REQUIRE_MESSAGE(loaded !=
nullptr,
error);
5018 BOOST_CHECK_EQUAL(std::string(loaded->ClassName()), std::string(original->ClassName()));
5022 BOOST_CHECK_CLOSE(loaded->Capacity(), original->Capacity(), 5.0);
5025 for (
const auto& p : probes) {
5028 original->Safety(p.data(), original->Contains(p.data())));
5029 for (
const auto&
d : directions) {
5031 original->DistFromOutside(p.data(),
d.data(), 3));
5033 original->DistFromInside(p.data(),
d.data(), 3));
5036 std::filesystem::remove(
path);
5043 const auto dir = std::filesystem::temp_directory_path();
5046 BOOST_CHECK(harness::loadShapeFromRootFile((dir /
"o2_shape_absent.root").
string(), &
error) ==
nullptr);
5052 const std::string
path = (dir /
"o2_shape_not_a_shape.root").
string();
5054 TFile out(
path.c_str(),
"RECREATE");
5055 TNamed impostor(
"shape",
"not a shape");
5056 out.WriteTObject(&impostor,
"shape");
5062 std::filesystem::remove(
path);
5072 constexpr double kR = 1.5;
5073 constexpr double kZ = 2.0;
5074 constexpr double kError = 0.05;
5075 const TGeoTube truth(
"truth", 0., kR, kZ);
5076 const TGeoTube wrong(
"wrong", 0., kR + kError, kZ);
5082 const auto trueContains = [](
const Point3D& p) {
5083 return (std::hypot(p[0], p[1]) <= kR && std::fabs(p[2]) <= kZ) ? 1 : 0;
5085 const auto trueBoundaryDistance = [](
const Point3D& p) {
5086 const double r = std::hypot(p[0], p[1]);
5087 const double dr = kR -
r;
5088 const double dz = kZ - std::fabs(p[2]);
5089 if (dr > 0. && dz > 0.) {
5090 return std::min(dr, dz);
5092 return std::hypot(std::max(
r - kR, 0.), std::max(std::fabs(p[2]) - kZ, 0.));
5095 std::vector<Point3D> points;
5096 std::vector<int> containsState;
5097 std::vector<double> boundaryDistance;
5098 for (
int ix = -6; ix <= 6; ++ix) {
5099 for (
int iy = -6; iy <= 6; ++iy) {
5100 for (
int iz = -4; iz <= 4; ++iz) {
5101 const Point3D p{0.31 * ix, 0.29 * iy, 0.53 * iz};
5105 if (std::fabs(trueBoundaryDistance(p)) < 1.e-3) {
5108 points.push_back(p);
5109 containsState.push_back(trueContains(p));
5110 boundaryDistance.push_back(trueBoundaryDistance(p));
5114 BOOST_REQUIRE_GT(points.size(), 500u);
5116 harness::ValidationOptions opt;
5117 opt.meshBand = 1.e-6;
5118 opt.distanceTolerance = 1.e-9;
5120 auto containsTruth = harness::validateContainsAgainstOracle(&truth, points, containsState,
5121 boundaryDistance, opt);
5122 auto containsWrong = harness::validateContainsAgainstOracle(&wrong, points, containsState,
5123 boundaryDistance, opt);
5124 BOOST_CHECK_EQUAL(containsTruth.nMismatchUnexplained + containsTruth.nMismatchMissedSurface, 0u);
5125 BOOST_CHECK_GT(containsWrong.nMismatchUnexplained + containsWrong.nMismatchMissedSurface, 0u);
5127 auto safetyTruth = harness::validateSafetyAgainstOracle(&truth, points, boundaryDistance, opt);
5128 auto safetyWrong = harness::validateSafetyAgainstOracle(&wrong, points, boundaryDistance, opt);
5129 BOOST_CHECK_EQUAL(safetyTruth.nMismatchUnexplained + safetyTruth.nMismatchMissedSurface, 0u);
5130 BOOST_CHECK_GT(safetyWrong.nMismatchUnexplained + safetyWrong.nMismatchMissedSurface, 0u);
5135 std::vector<harness::Ray> rays;
5136 std::vector<double> rayDistance;
5137 std::vector<int> originState;
5138 for (
int i = 0;
i < 400; ++
i) {
5139 const double phi = 0.0173 *
i;
5140 const double z = -1.9 + 0.0095 *
i;
5141 const Point3D target{0.9 * kR * std::cos(2.1 * phi), 0.9 * kR * std::sin(2.1 * phi),
z};
5142 const Point3D origin{5.0 * std::cos(phi), 5.0 * std::sin(phi), 3.0 - 0.01 *
i};
5144 const double norm = std::sqrt(dir[0] * dir[0] + dir[1] * dir[1] + dir[2] * dir[2]);
5145 for (
auto& component : dir) {
5148 rays.push_back(harness::Ray{
origin, dir});
5149 rayDistance.push_back(truth.DistFromOutside(
origin.data(), dir.data(), 3));
5150 originState.push_back(trueContains(
origin));
5153 auto distTruth = harness::validateDistanceAgainstOracle(&truth, rays, rayDistance,
5154 false, opt, originState);
5155 auto distWrong = harness::validateDistanceAgainstOracle(&wrong, rays, rayDistance,
5156 false, opt, originState);
5157 BOOST_CHECK_EQUAL(distTruth.nMismatchUnexplained + distTruth.nMismatchMissedSurface, 0u);
5158 BOOST_CHECK_GT(distWrong.nMismatchUnexplained + distWrong.nMismatchMissedSurface, 0u);
5179TGeoCompositeShape* makePlacedTube(
const char*
name,
double rmin,
double rmax,
double dz,
5180 TGeoMatrix* matrixA, TGeoMatrix* matrixB)
5182 auto*
left =
new TGeoTube(Form(
"%s_l",
name), rmin, rmax, dz);
5183 auto*
right =
new TGeoTube(Form(
"%s_r",
name), rmin, rmax, dz);
5185 return new TGeoCompositeShape(
name,
node);
5194 constexpr double kRmin = 0.4;
5195 constexpr double kRmax = 1.0;
5196 constexpr double kDz = 5.0;
5197 const TGeoTube
reference(
"reference", kRmin, kRmax, kDz);
5199 auto*
rotation =
new TGeoRotation(
"csgRot", 0., 0., 0.);
5202 auto* matrixA =
new TGeoCombiTrans(0.3, 5.916, 2.0,
rotation);
5203 auto* matrixB =
new TGeoCombiTrans(0.3, 5.916, 2.0,
rotation);
5204 const TGeoCombiTrans placement(0.3, 5.916, 2.0,
rotation);
5205 std::unique_ptr<TGeoCompositeShape> placed(
5206 makePlacedTube(
"csgPlaced", kRmin, kRmax, kDz, matrixA, matrixB));
5208 std::size_t probes = 0;
5209 std::size_t inside = 0;
5210 std::size_t outside = 0;
5211 for (
int ix = -8; ix <= 8; ++ix) {
5212 for (
int iy = -8; iy <= 8; ++iy) {
5213 for (
int iz = -8; iz <= 8; ++iz) {
5214 const Point3D master{0.3 + 0.37 * ix, 5.916 + 0.41 * iy, 2.0 + 0.43 * iz};
5216 placement.MasterToLocal(master.data(), local.data());
5219 const double r = std::hypot(local[0], local[1]);
5220 if (std::fabs(
r - kRmin) < 1.e-9 || std::fabs(
r - kRmax) < 1.e-9 ||
5221 std::fabs(std::fabs(local[2]) - kDz) < 1.e-9) {
5225 const bool wanted =
reference.Contains(local.data());
5226 BOOST_REQUIRE_EQUAL(placed->Contains(master.data()), wanted);
5227 BOOST_REQUIRE_CLOSE_FRACTION(placed->Safety(master.data(), wanted),
5228 reference.Safety(local.data(), wanted), 1.e-12);
5229 wanted ? ++inside : ++outside;
5232 Point3D{0.5773502691896258, 0.5773502691896258, 0.5773502691896258}}) {
5234 placement.MasterToLocalVect(dir.data(), localDir.data());
5236 BOOST_REQUIRE_CLOSE_FRACTION(placed->DistFromInside(master.data(), dir.data(), 3),
5237 reference.DistFromInside(local.data(), localDir.data(), 3),
5240 const double got = placed->DistFromOutside(master.data(), dir.data(), 3);
5241 const double want =
reference.DistFromOutside(local.data(), localDir.data(), 3);
5243 BOOST_REQUIRE_GT(got, 1.e20);
5245 BOOST_REQUIRE_CLOSE_FRACTION(got, want, 1.e-12);
5253 BOOST_CHECK_GT(probes, 2000u);
5254 BOOST_CHECK_GT(inside, 100u);
5255 BOOST_CHECK_GT(outside, 100u);
5261 const TGeoTube wrong(
"wrongReference", kRmin, kRmax + 0.05, kDz);
5262 std::size_t disagreements = 0;
5263 std::size_t shellProbes = 0;
5264 for (
int iphi = 0; iphi < 24; ++iphi) {
5265 const double phi = 2. * M_PI * iphi / 24.;
5266 for (
int iz = -3; iz <= 3; ++iz) {
5267 const double radius = kRmax + 0.025;
5270 placement.LocalToMaster(local.data(), master.data());
5272 disagreements += (placed->Contains(master.data()) != wrong.Contains(local.data())) ? 1 : 0;
5286 constexpr double kEyeRmin = 0.7;
5287 constexpr double kEyeRmax = 1.2;
5288 constexpr double kEyeDz = 0.75;
5289 constexpr double kRodR = 0.6;
5290 constexpr double kRodDz = 3.5;
5291 constexpr double kRodCentre = 3.5;
5293 auto* eyeRotation =
new TGeoRotation(
"csgEyeRot", 90., 90., 0.);
5294 auto* eyeMatrix =
new TGeoCombiTrans(0., 0., 0., eyeRotation);
5295 auto* rodMatrix =
new TGeoTranslation(0., 0., kRodCentre);
5296 auto* eye =
new TGeoTube(
"csgEye", kEyeRmin, kEyeRmax, kEyeDz);
5297 auto* rod =
new TGeoTube(
"csgRod", 0., kRodR, kRodDz);
5298 auto*
node =
new TGeoUnion(eye, rod, eyeMatrix, rodMatrix);
5299 std::unique_ptr<TGeoCompositeShape> ram(
new TGeoCompositeShape(
"csgRam",
node));
5301 const auto closedForm = [&](
const Point3D& p) {
5302 const double rEye = std::hypot(p[1], p[2]);
5303 const bool inEye = rEye >= kEyeRmin && rEye <= kEyeRmax && std::fabs(p[0]) <= kEyeDz;
5304 const double rRod = std::hypot(p[0], p[1]);
5305 const bool inRod = rRod <= kRodR && p[2] >= 0. && p[2] <= 2. * kRodDz;
5306 return inEye || inRod;
5308 const auto nearWall = [&](
const Point3D& p) {
5309 const double rEye = std::hypot(p[1], p[2]);
5310 const double rRod = std::hypot(p[0], p[1]);
5311 return std::fabs(rEye - kEyeRmin) < 1.e-9 || std::fabs(rEye - kEyeRmax) < 1.e-9 ||
5312 std::fabs(std::fabs(p[0]) - kEyeDz) < 1.e-9 || std::fabs(rRod - kRodR) < 1.e-9 ||
5313 std::fabs(p[2]) < 1.e-9 || std::fabs(p[2] - 2. * kRodDz) < 1.e-9;
5316 const std::filesystem::path
path =
5317 std::filesystem::temp_directory_path() /
"o2_csg_ram_shape.root";
5320 BOOST_REQUIRE_MESSAGE(harness::saveShapeToRootFile(
path.string(), *ram, &
error),
error);
5321 std::unique_ptr<TGeoShape> loaded(harness::loadShapeFromRootFile(
path.string(), &
error));
5322 BOOST_REQUIRE_MESSAGE(loaded !=
nullptr,
error);
5323 BOOST_CHECK_EQUAL(std::string(loaded->ClassName()), std::string(
"TGeoCompositeShape"));
5325 std::size_t inside = 0;
5326 std::size_t outside = 0;
5327 for (
int ix = -6; ix <= 6; ++ix) {
5328 for (
int iy = -6; iy <= 6; ++iy) {
5329 for (
int iz = -4; iz <= 20; ++iz) {
5330 const Point3D p{0.23 * ix, 0.27 * iy, 0.41 * iz};
5334 const bool wanted = closedForm(p);
5335 BOOST_REQUIRE_EQUAL(ram->Contains(p.data()), wanted);
5336 BOOST_REQUIRE_EQUAL(loaded->Contains(p.data()), wanted);
5337 wanted ? ++inside : ++outside;
5341 BOOST_CHECK_GT(inside, 50u);
5342 BOOST_CHECK_GT(outside, 500u);
5351 double minCapacity = ram->Capacity();
5352 double maxCapacity = minCapacity;
5353 for (
int i = 0;
i < 5; ++
i) {
5354 const double sampled = ram->Capacity();
5355 minCapacity = std::min(minCapacity, sampled);
5356 maxCapacity = std::max(maxCapacity, sampled);
5358 BOOST_CHECK_GT(minCapacity, 0.);
5359 const double spread = (maxCapacity - minCapacity) / (0.5 * (maxCapacity + minCapacity));
5360 BOOST_CHECK_GT(spread, 1.e-4);
5362 std::filesystem::remove(
path);
5397 BOOST_REQUIRE_EQUAL(boxCrossings.size(), 2u);
5398 BOOST_CHECK_SMALL(boxCrossings[0].
t - 4., 1.e-12);
5399 BOOST_CHECK_SMALL(boxCrossings[1].
t - 6., 1.e-12);
5403 TGeoTube tube(
"xrayTube", 0.5, 1.0, 2.0);
5405 BOOST_REQUIRE_EQUAL(tubeCrossings.size(), 4u);
5406 const double expected[4] = {4.0, 4.5, 5.5, 6.0};
5407 const int senses[4] = {+1, -1, +1, -1};
5408 for (
int i = 0;
i < 4; ++
i) {
5409 BOOST_CHECK_SMALL(tubeCrossings[
i].
t -
expected[
i], 1.e-12);
5427 std::array<std::pair<std::unique_ptr<SurfaceSolid>,
double>, 7> fixtures{{
5428 {makeBoxSolid(
"xrayLoopBox", 1., 2., 3.), 4.},
5429 {makeTubeSolid(
"xrayLoopTube", 0., 2., 3.), 4.},
5430 {makeTubeSolid(
"xrayLoopHollowTube", 1., 2., 3.), 4.},
5431 {makeConeSolid(
"xrayLoopCone", 2., 1., 3.), 4.},
5432 {makeSphereSolid(
"xrayLoopSphere", 2.5), 3.5},
5433 {makeTorusSolid(
"xrayLoopTorus", 3., 1.), 4.5},
5434 {makeCapsuleSolid(
"xrayLoopCapsule", 1., 1.5), 3.},
5436 size_t comparedRays = 0;
5437 size_t comparedCrossings = 0;
5438 for (
const auto& [solid, extent] : fixtures) {
5439 BOOST_TEST_CONTEXT(
"fixture = " << solid->GetName())
5441 BOOST_REQUIRE(solid->HasBVH());
5442 const XRayPoint lo{-extent, -extent, -extent};
5443 const XRayPoint hi{extent, extent, extent};
5447 for (
const auto& ray : raster.
rays) {
5451 stepWithShapeApi(solid.get(), ray.origin, ray.dir, ray.tMax, cfg, bvhStats);
5454 ray.origin, ray.dir, ray.tMax, cfg, loopStats,
5455 [&solid](
const double* p) { return solid->Contains_Loop(p); },
5456 [&solid](
const double* p,
const double*
d) { return solid->DistFromOutside_Loop(p, d); },
5457 [&solid](
const double* p,
const double*
d) { return solid->DistFromInside_Loop(p, d); });
5458 BOOST_REQUIRE_EQUAL(viaBVH.size(), viaLoop.size());
5459 for (
size_t i = 0;
i < viaBVH.size(); ++
i) {
5465 comparedCrossings += viaBVH.size();
5469 BOOST_CHECK_GT(comparedRays, 2000u);
5470 BOOST_CHECK_GT(comparedCrossings, 2000u);
5478 const std::vector<Crossing> truth{{4.0, +1}, {4.5, -1}, {5.5, +1}, {6.0, -1}};
5479 const double tolerance = 1.e-6;
5489 auto perturbed = truth;
5490 perturbed[2].t += 1.e-3;
5492 compareLists(perturbed, truth, {}, {}, tolerance, displaced);
5497 BOOST_CHECK_SMALL(displaced.
worstDeltaT - 1.e-3, 1.e-12);
5499 auto dropped = truth;
5500 dropped.erase(dropped.begin() + 1);
5506 auto doubled = truth;
5507 doubled.insert(doubled.begin() + 1, {4.2, -1});
5509 compareLists(doubled, truth, {}, {}, tolerance, spurious);
5513 auto flipped = truth;
5514 flipped[1].kind = +1;
5516 compareLists(flipped, truth, {}, {}, tolerance, sense);
5521 auto nudged = truth;
5522 nudged[0].t += 1.e-9;
5551 auditCrossingList({{1.0, +1}, {2.0, -1}, {4.0, +1}, {6.0, -1}}, &
box,
origin, dir, 10., cfg,
5564 for (
const int n : {5, 16, 41}) {
5567 std::vector<double> byBeam(raster.
beams.size(), 0.);
5568 for (
const auto& ray : raster.
rays) {
5569 const double before = stats.insideLength;
5571 auditCrossingList(crossings,
nullptr, ray.origin, ray.dir, ray.tMax, cfg, stats);
5572 byBeam[ray.beam] += stats.insideLength - before;
5574 BOOST_CHECK_SMALL(
chordVolume(raster, byBeam) - 24., 1.e-9);
5581 TGeoSphere sphere(
"xrayVolSphere", 0., 1.);
5582 const double exact = 4. / 3. * 3.14159265358979323846;
5583 for (
const int n : {24, 96}) {
5586 for (
const auto& ray : raster.
rays) {
5587 const auto crossings =
stepWithShapeApi(&sphere, ray.origin, ray.dir, ray.tMax, cfg, stats);
5588 auditCrossingList(crossings,
nullptr, ray.origin, ray.dir, ray.tMax, cfg, stats);
5590 const double volume = stats.insideLength * raster.
cellArea[0];
5591 BOOST_CHECK_LT(std::fabs(volume - exact) / exact, 2.e-3);
5604 for (
const auto& ray : raster.
rays) {
5605 BOOST_REQUIRE(!
box.Contains(ray.origin.data()));
5607 const double end[3] = {ray.origin[0] + ray.tMax * ray.dir[0],
5608 ray.origin[1] + ray.tMax * ray.dir[1],
5609 ray.origin[2] + ray.tMax * ray.dir[2]};
5610 BOOST_REQUIRE(!
box.Contains(
end));
5617 for (
const auto& beams : {axes, fan}) {
5618 for (
const auto& beam : beams) {
5619 BOOST_CHECK_SMALL(
dot3(beam.dir, beam.dir) - 1., 1.e-12);
5620 BOOST_CHECK_SMALL(
dot3(beam.u, beam.v), 1.e-12);
5621 BOOST_CHECK_SMALL(
dot3(beam.u, beam.dir), 1.e-12);
5622 BOOST_CHECK_SMALL(
dot3(beam.v, beam.dir), 1.e-12);
5627 double worstFanAlignment = -1.;
5628 for (
size_t i = 0;
i < fan.size(); ++
i) {
5629 for (
size_t j =
i + 1;
j < fan.size(); ++
j) {
5630 worstFanAlignment = std::max(worstFanAlignment, std::fabs(
dot3(fan[
i].dir, fan[
j].dir)));
5633 BOOST_CHECK_LT(worstFanAlignment, 0.999);
5668double quarticBackwardError(
const std::array<double, 5>& coefficients,
double x)
5670 double value = 0., magnitude = 0., power = 1.;
5671 for (
int i = 0;
i < 5; ++
i) {
5672 const double term = coefficients[4 -
i] * power;
5674 magnitude += std::abs(term);
5677 return magnitude > 0. ? std::abs(
value) / magnitude : std::abs(
value);
5681std::array<double, 5> quarticFromRoots(
double r1,
double r2,
double r3,
double r4)
5683 return {1., -(r1 + r2 + r3 + r4),
5684 r1 * r2 + r1 * r3 + r1 * r4 + r2 * r3 + r2 * r4 + r3 * r4,
5685 -(r1 * r2 * r3 + r1 * r2 * r4 + r1 * r3 * r4 + r2 * r3 * r4), r1 * r2 * r3 * r4};
5688std::vector<double> sortedRoots(
const std::array<double, 5>&
c,
surf::QuarticBranch* branch =
nullptr)
5691 std::vector<double> roots(found.begin(), found.end());
5692 std::sort(roots.begin(), roots.end());
5697void checkRootsAreRoots(
const std::array<double, 5>&
c,
const std::vector<double>& roots)
5699 for (
const double root : roots) {
5700 BOOST_CHECK_LT(quarticBackwardError(
c, root), 1.e-12);
5716 const std::array<double, 5>
c{1.0, -1501.7280000044018, 845808.25396968238, -211752288.545858,
5717 19882619385.616932};
5722 const double firstRoot = 375.3392295779947145;
5723 const double secondRoot = 375.5247704240909448;
5730 const auto roots = sortedRoots(
c, &branch);
5731 BOOST_REQUIRE_EQUAL(roots.size(), 2u);
5732 checkClose(roots[0], firstRoot, 1.e-6);
5733 checkClose(roots[1], secondRoot, 1.e-6);
5734 checkRootsAreRoots(
c, roots);
5737 BOOST_CHECK(branch == surf::QuarticBranch::Biquadratic);
5749 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}) {
5750 const auto c = quarticFromRoots(1. * k, 2. * k, 3. * k, 7. * k);
5752 const auto roots = sortedRoots(
c, &branch);
5753 BOOST_REQUIRE_EQUAL(roots.size(), 4u);
5754 const double expected[4] = {1. * k, 2. * k, 3. * k, 7. * k};
5755 for (
int i = 0;
i < 4; ++
i) {
5758 checkRootsAreRoots(
c, roots);
5761 BOOST_CHECK(branch == surf::QuarticBranch::Resolvent);
5768 for (
const double k : {1., 1.e-2, 1.e-4, 1.e-6}) {
5769 const auto c = quarticFromRoots(-2. * k, -1. * k, 1. * k, 2.1 * k);
5770 const auto roots = sortedRoots(
c);
5771 BOOST_REQUIRE_EQUAL(roots.size(), 4u);
5772 const double expected[4] = {-2. * k, -1. * k, 1. * k, 2.1 * k};
5773 for (
int i = 0;
i < 4; ++
i) {
5776 checkRootsAreRoots(
c, roots);
5788 const std::array<double, 5>
c{1., 0., -5., 0., 4.};
5790 const auto roots = sortedRoots(
c, &branch);
5791 BOOST_CHECK(branch == surf::QuarticBranch::Biquadratic);
5792 BOOST_REQUIRE_EQUAL(roots.size(), 4u);
5793 const double expected[4] = {-2., -1., 1., 2.};
5794 for (
int i = 0;
i < 4; ++
i) {
5808 for (
const double centre : {0., 1., 100.}) {
5809 for (
const double k : {1., 1.e-3, 1.e3}) {
5810 const auto c = quarticFromRoots((centre - 2.) * k, (centre - 1.) * k, (centre + 1.) * k,
5813 const auto roots = sortedRoots(
c, &branch);
5814 BOOST_CHECK(branch == surf::QuarticBranch::Biquadratic);
5815 BOOST_REQUIRE_EQUAL(roots.size(), 4u);
5816 const double expected[4] = {(centre - 2.) * k, (centre - 1.) * k, (centre + 1.) * k,
5818 for (
int i = 0;
i < 4; ++
i) {
5819 BOOST_CHECK_LT(std::abs(roots[
i] -
expected[
i]), 1.e-9 * std::max(1.e-30, std::abs(
expected[
i])));
5821 checkRootsAreRoots(
c, roots);
5836 BOOST_CHECK(branch == surf::QuarticBranch::NotAQuartic);
5838 for (
const double k : {1.e4, 1., 1.e-4, 1.e-8}) {
5840 const double k2 = k * k;
5841 const std::array<double, 5> biquadratic{1., 0., 5. * k2, 0., 4. * k2 * k2};
5844 BOOST_CHECK(branch == surf::QuarticBranch::Biquadratic);
5848 const std::array<double, 5> asymmetric{1., 2. * k, 6. * k2, 2. * k2 * k, 5. * k2 * k2};
5850 BOOST_CHECK(branch == surf::QuarticBranch::Resolvent);
5863 for (
const double k : {1., 1.e-2, 1.e-4, 1.e-6}) {
5864 for (
const double delta : {1.e-1, 1.e-3, 1.e-6, 1.e-9, 1.e-12, 0.}) {
5865 const auto c = quarticFromRoots(-2. * k, -1. * k, 1. * k, (2. + delta) * k);
5866 const auto roots = sortedRoots(
c);
5867 BOOST_REQUIRE_EQUAL(roots.size(), 4u);
5868 const double expected[4] = {-2. * k, -1. * k, 1. * k, (2. + delta) * k};
5869 for (
int i = 0;
i < 4; ++
i) {
5872 checkRootsAreRoots(
c, roots);
5894 const Point3D centre{0.3, -0.7, 1.1};
5895 const Point3D axis{0., 0., 1.};
5896 const Point3D refU{1., 0., 0.};
5897 constexpr double radius = 2.5;
5898 constexpr double phi = 0.9;
5901 const double cx = centre[0] +
radius * std::cos(phi);
5902 const double cy = centre[1] +
radius * std::sin(phi);
5903 const Double_t onCylinder[3] = {cx, cy, centre[2] + 0.4};
5904 const Double_t awayFromAxis[3] = {std::cos(phi), std::sin(phi), 0.};
5906 for (
const bool innerWall : {
false,
true}) {
5907 SurfaceSolid solid(innerWall ?
"streamK_innerCyl" :
"streamK_outerCyl");
5908 BOOST_REQUIRE(solid.AddCylindricalSurface(centre, axis, refU,
radius, -1., 1., 0.,
surf::kTwoPi, innerWall));
5909 Double_t
n[3] = {0., 0., 0.};
5910 solid.ComputeNormal(onCylinder,
nullptr,
n);
5911 const double alignment =
n[0] * awayFromAxis[0] +
n[1] * awayFromAxis[1] +
n[2] * awayFromAxis[2];
5913 BOOST_CHECK_CLOSE(alignment, innerWall ? -1. : 1., 1.e-9);
5920 for (
const bool innerWall : {
false,
true}) {
5921 SurfaceSolid solid(innerWall ?
"streamK_innerCone" :
"streamK_outerCone");
5923 BOOST_REQUIRE(solid.AddConicalSurface(centre, axis, refU, 1., 3., 0., 2., 0.,
surf::kTwoPi, innerWall));
5924 const double h = 1.0;
5925 const double r = 2.0;
5926 const Double_t onCone[3] = {centre[0] +
r * std::cos(phi), centre[1] +
r * std::sin(phi), centre[2] +
h};
5927 Double_t
n[3] = {0., 0., 0.};
5928 solid.ComputeNormal(onCone,
nullptr,
n);
5932 const double radial =
n[0] * std::cos(phi) +
n[1] * std::sin(phi);
5933 BOOST_CHECK_GT(innerWall ? -radial : radial, 0.5);
5936 for (
const bool innerWall : {
false,
true}) {
5937 SurfaceSolid solid(innerWall ?
"streamK_innerSph" :
"streamK_outerSph");
5939 const double theta = 1.1;
5940 const Double_t onSphere[3] = {centre[0] +
radius * std::sin(theta) * std::cos(phi),
5941 centre[1] +
radius * std::sin(theta) * std::sin(phi),
5942 centre[2] +
radius * std::cos(theta)};
5943 const double outward[3] = {std::sin(theta) * std::cos(phi), std::sin(theta) * std::sin(phi), std::cos(theta)};
5944 Double_t
n[3] = {0., 0., 0.};
5945 solid.ComputeNormal(onSphere,
nullptr,
n);
5946 const double alignment =
n[0] * outward[0] +
n[1] * outward[1] +
n[2] * outward[2];
5947 BOOST_CHECK_CLOSE(alignment, innerWall ? -1. : 1., 1.e-9);
5973TGeoCombiTrans* makeStreamNPlacement()
5975 auto*
rotation =
new TGeoRotation(
"streamNRot", 0., 0., 0.);
5979 return new TGeoCombiTrans(0.3, 5.916, 2.0,
rotation);
5986 const auto dir = std::filesystem::temp_directory_path();
5987 const std::string
path = (dir /
"o2_shape_placed.root").
string();
5989 const TGeoTube tube(
"shape", 0.4, 1.0, 5.0);
5990 std::unique_ptr<TGeoCombiTrans> placement(makeStreamNPlacement());
5993 BOOST_REQUIRE_MESSAGE(harness::saveShapeToRootFile(
path, tube, placement.get(), &
error),
error);
5995 std::unique_ptr<TGeoShape> loaded(harness::loadShapeFromRootFile(
path, &
error));
5996 BOOST_REQUIRE_MESSAGE(loaded !=
nullptr,
error);
5999 std::unique_ptr<TGeoHMatrix> back(harness::loadShapePlacementFromRootFile(
path));
6000 BOOST_REQUIRE(back !=
nullptr);
6001 for (
int i = 0;
i < 9; ++
i) {
6004 for (
int i = 0;
i < 3; ++
i) {
6009 const Point3D master{0.9, 6.2, 3.1};
6012 back->MasterToLocal(master.data(), viaFile.data());
6013 placement->MasterToLocal(master.data(), viaOriginal.data());
6014 for (
int i = 0;
i < 3; ++
i) {
6017 std::filesystem::remove(
path);
6023 const auto dir = std::filesystem::temp_directory_path();
6024 const TGeoTube tube(
"shape", 0.4, 1.0, 5.0);
6029 const std::string legacy = (dir /
"o2_shape_legacy.root").
string();
6030 BOOST_REQUIRE_MESSAGE(harness::saveShapeToRootFile(legacy, tube, &
error),
error);
6031 BOOST_CHECK(harness::loadShapePlacementFromRootFile(legacy) ==
nullptr);
6035 const std::string identity = (dir /
"o2_shape_identity.root").
string();
6036 TGeoHMatrix unit(
"unit");
6037 BOOST_REQUIRE_MESSAGE(harness::saveShapeToRootFile(identity, tube, &unit, &
error),
error);
6038 BOOST_CHECK(harness::loadShapePlacementFromRootFile(identity) ==
nullptr);
6042 BOOST_CHECK(harness::loadShapePlacementFromRootFile((dir /
"o2_shape_nothing.root").
string()) ==
6045 std::filesystem::remove(legacy);
6046 std::filesystem::remove(identity);
6053 constexpr double kRmin = 0.4;
6054 constexpr double kRmax = 1.0;
6055 constexpr double kDz = 5.0;
6057 std::unique_ptr<TGeoCombiTrans> placement(makeStreamNPlacement());
6058 const TGeoTube placedPrimitive(
"streamNTube", kRmin, kRmax, kDz);
6060 std::unique_ptr<TGeoCompositeShape> composite(
6061 makePlacedTube(
"streamNComposite", kRmin, kRmax, kDz,
new TGeoCombiTrans(*placement),
6062 new TGeoCombiTrans(*placement)));
6064 std::size_t probes = 0;
6065 std::size_t inside = 0;
6066 std::size_t disagreements = 0;
6069 const TGeoTube wrong(
"streamNWrong", kRmin, kRmax * 1.05, kDz);
6070 std::size_t controlDisagreements = 0;
6072 for (
int ix = -8; ix <= 8; ++ix) {
6073 for (
int iy = -8; iy <= 8; ++iy) {
6074 for (
int iz = -8; iz <= 8; ++iz) {
6075 const Point3D master{0.3 + 0.37 * ix, 5.916 + 0.41 * iy, 2.0 + 0.43 * iz};
6077 placement->MasterToLocal(master.data(), local.data());
6078 const double r = std::hypot(local[0], local[1]);
6079 if (std::fabs(
r - kRmin) < 1.e-9 || std::fabs(
r - kRmax) < 1.e-9 ||
6080 std::fabs(
r - kRmax * 1.05) < 1.e-9 || std::fabs(std::fabs(local[2]) - kDz) < 1.e-9) {
6084 const bool wanted = composite->Contains(master.data());
6085 if (placedPrimitive.Contains(local.data()) != wanted) {
6088 if (wrong.Contains(local.data()) != wanted) {
6089 ++controlDisagreements;
6094 BOOST_REQUIRE_CLOSE_FRACTION(placedPrimitive.Safety(local.data(), wanted),
6095 composite->Safety(master.data(), wanted), 1.e-12);
6097 Point3D{0.5773502691896258, 0.5773502691896258,
6098 0.5773502691896258}}) {
6100 placement->MasterToLocalVect(dir.data(), localDir.data());
6102 BOOST_REQUIRE_CLOSE_FRACTION(placedPrimitive.DistFromInside(local.data(),
6103 localDir.data(), 3),
6104 composite->DistFromInside(master.data(), dir.data(), 3),
6107 const double got = placedPrimitive.DistFromOutside(local.data(), localDir.data(), 3);
6108 const double want = composite->DistFromOutside(master.data(), dir.data(), 3);
6110 BOOST_REQUIRE_GT(got, 1.e20);
6112 BOOST_REQUIRE_CLOSE_FRACTION(got, want, 1.e-12);
6120 BOOST_CHECK_GT(controlDisagreements, 0u);
6121 BOOST_CHECK_GT(inside, 100u);
6122 BOOST_CHECK_GT(probes, 3000u);
6128 constexpr double kRmin = 0.4;
6129 constexpr double kRmax = 1.0;
6130 constexpr double kDz = 5.0;
6131 const double analytic = TMath::Pi() * (kRmax * kRmax - kRmin * kRmin) * 2. * kDz;
6133 const TGeoTube tube(
"streamNCapTube", kRmin, kRmax, kDz);
6134 BOOST_CHECK_CLOSE_FRACTION(tube.Capacity(), analytic, 1.e-14);
6138 std::unique_ptr<TGeoCombiTrans> placement(makeStreamNPlacement());
6139 std::unique_ptr<TGeoCompositeShape> composite(
6140 makePlacedTube(
"streamNCapComposite", kRmin, kRmax, kDz,
new TGeoCombiTrans(*placement),
6141 new TGeoCombiTrans(*placement)));
6146 const double first = composite->Capacity();
6147 const double second = composite->Capacity();
6148 BOOST_CHECK_NE(
first, second);
6149 BOOST_CHECK_GT(std::fabs(
first - analytic) / analytic, 1.e-6);
6160 constexpr double kRmin = 0.4;
6161 constexpr double kRmax = 1.0;
6162 constexpr double kDz = 5.0;
6164 std::unique_ptr<TGeoCombiTrans> shapePlacementOwned(makeStreamNPlacement());
6165 const TGeoHMatrix shapePlacement(*shapePlacementOwned);
6166 auto* partRotation =
new TGeoRotation(
"streamNPartRot", 37., 24., 61.);
6167 const TGeoCombiTrans partPlacement(-2.0, 7.0, 1.5, partRotation);
6172 partPlacement.MasterToLocal(master.data(), partFrame.data());
6173 shapePlacement.MasterToLocal(partFrame.data(), shapeFrame.data());
6174 const double r = std::hypot(shapeFrame[0], shapeFrame[1]);
6175 onWall = std::fabs(
r - kRmin) < 1.e-9 || std::fabs(
r - kRmax) < 1.e-9 ||
6176 std::fabs(std::fabs(shapeFrame[2]) - kDz) < 1.e-9;
6177 return r >= kRmin &&
r <= kRmax && std::fabs(shapeFrame[2]) <= kDz;
6182 TGeoHMatrix correct(partPlacement);
6183 correct.Multiply(&shapePlacement);
6184 TGeoHMatrix reversed(shapePlacement);
6185 reversed.Multiply(&partPlacement);
6186 TGeoHMatrix transposedRotation(shapePlacement);
6189 const double*
r = shapePlacement.GetRotationMatrix();
6190 for (
int i = 0;
i < 3; ++
i) {
6191 for (
int j = 0;
j < 3; ++
j) {
6192 rt[3 *
i +
j] =
r[3 *
j +
i];
6195 transposedRotation.SetRotation(rt);
6196 transposedRotation.SetBit(TGeoMatrix::kGeoRotation);
6198 TGeoHMatrix withTransposed(partPlacement);
6199 withTransposed.Multiply(&transposedRotation);
6200 const TGeoHMatrix partOnly(partPlacement);
6202 const std::vector<std::pair<std::string, const TGeoHMatrix*>> candidates{
6203 {
"part*shape", &correct},
6204 {
"shape*part", &reversed},
6205 {
"part*shape^T", &withTransposed},
6206 {
"part only", &partOnly}};
6212 const double* centre = correct.GetTranslation();
6214 std::vector<size_t> disagreements(candidates.size(), 0);
6216 size_t insideProbes = 0;
6218 for (
size_t c = 0;
c < candidates.size(); ++
c) {
6221 auto* manager =
new TGeoManager((
"streamN_" +
std::to_string(
c)).c_str(),
"composition order");
6222 auto* material =
new TGeoMaterial(
"Vacuum", 0., 0., 0.);
6223 auto* medium =
new TGeoMedium(
"Vacuum", 1, material);
6224 auto* world =
new TGeoVolume(
"TOP",
new TGeoBBox(
"streamNWorld", 30., 30., 30.), medium);
6225 auto* part =
new TGeoVolume(
"PART",
new TGeoTube(
"streamNNodeTube", kRmin, kRmax, kDz), medium);
6226 world->AddNode(part, 1,
new TGeoHMatrix(*candidates[
c].second));
6227 manager->SetTopVolume(world);
6228 manager->CloseGeometry();
6230 size_t localProbes = 0;
6231 size_t localInside = 0;
6232 for (
int ix = -14; ix <= 14; ++ix) {
6233 for (
int iy = -14; iy <= 14; ++iy) {
6234 for (
int iz = -14; iz <= 14; ++iz) {
6235 const Point3D master{centre[0] + 0.45 * ix, centre[1] + 0.47 * iy,
6236 centre[2] + 0.43 * iz};
6237 bool onWall =
false;
6238 const bool wanted =
reference(master, onWall);
6246 TGeoNode*
node = manager->FindNode(master[0], master[1], master[2]);
6247 const bool got =
node !=
nullptr && std::string(
node->GetVolume()->GetName()) ==
"PART";
6248 if (got != wanted) {
6254 probes = localProbes;
6255 insideProbes = localInside;
6257 gGeoManager =
nullptr;
6261 BOOST_CHECK_GT(probes, 5000u);
6262 BOOST_CHECK_GT(insideProbes, 200u);
6264 BOOST_CHECK_GT(disagreements[1], 0u);
6265 BOOST_CHECK_GT(disagreements[2], 0u);
6266 BOOST_CHECK_GT(disagreements[3], 0u);
6286 TGeoBBox fast(
"repBenchFast", 1., 1., 1.);
6293 BOOST_CHECK_GT(
samples.insidePoints, 100);
6294 BOOST_CHECK_LT(
samples.insidePoints,
static_cast<long long>(
samples.points.size()) - 100);
6314 TGeoTube tube(
"repBenchTube", 0.3, 1., 2.);
6322 BOOST_REQUIRE_EQUAL(
a.points.size(),
b.points.size());
6323 for (
size_t i = 0;
i <
a.points.size(); ++
i) {
6340 constexpr size_t kBytes = 64u << 20;
6341 auto block = std::make_unique<char[]>(kBytes);
6342 for (
size_t i = 0;
i < kBytes;
i += 4096) {
6343 block[
i] =
static_cast<char>(
i);
6349 BOOST_CHECK_LT((
readMemory() - before).heapInUseBytes, 8LL << 20);
6358 auto* manager =
new TGeoManager(
"repBenchLadder",
"ladder");
6359 for (
const int k : {2, 4, 8, 16, 32}) {
6374 BOOST_CHECK(
dynamic_cast<TGeoCompositeShape*
>(single) ==
nullptr);
6377 gGeoManager =
nullptr;
6395struct SampleStream {
6396 explicit SampleStream(std::uint64_t seed) { lcg.state =
seed | 1u; }
6398 double uniform() {
return lcg.
next(); }
6399 double symmetric(
double extent) {
return (2. * uniform() - 1.) * extent; }
6409std::vector<std::array<double, 3>> nearestPatchSample(
const SurfaceSolid& solid,
double extent,
int count)
6411 SampleStream
stream(0x5EAFE7Full);
6412 std::vector<std::array<double, 3>> points;
6413 points.reserve(
static_cast<size_t>(
count));
6415 std::array<double, 3> point{
stream.symmetric(extent),
stream.symmetric(extent),
stream.symmetric(extent)};
6416 switch (
index % 5) {
6420 std::array<double, 3> normal{0., 0., 0.};
6421 const double safety = solid.Safety_Loop(point.data(), solid.Contains(point.data()));
6422 solid.ComputeNormal_Loop(point.data(),
nullptr, normal.data());
6423 const double sign = solid.Contains(point.data()) ? 1. : -1.;
6424 for (
int dimension = 0; dimension < 3; ++dimension) {
6425 point[dimension] += sign * safety * normal[dimension];
6430 std::array<double, 3> normal{0., 0., 0.};
6431 const double safety = solid.Safety_Loop(point.data(), solid.Contains(point.data()));
6432 solid.ComputeNormal_Loop(point.data(),
nullptr, normal.data());
6433 const double sign = solid.Contains(point.data()) ? 1. : -1.;
6434 const double offset = safety - sign * 1.e-9 * std::max(1., extent);
6435 for (
int dimension = 0; dimension < 3; ++dimension) {
6436 point[dimension] += sign *
offset * normal[dimension];
6451 points.push_back(point);
6454 points.push_back({0., 0., 0.});
6460int countNearestPatchDisagreements(
const SurfaceSolid& solid,
const std::array<double, 3>& point,
6461 double* worstSafetyGap =
nullptr)
6463 int disagreements = 0;
6464 for (
const bool inside : {
true,
false}) {
6465 const double accelerated = solid.Safety(point.data(), inside);
6466 const double reference = solid.Safety_Loop(point.data(), inside);
6470 if (worstSafetyGap !=
nullptr) {
6471 *worstSafetyGap = std::max(*worstSafetyGap, accelerated -
reference);
6476 const std::array<double, 3>
direction = unitDirection(0.37, -0.82, 0.44);
6477 for (
const double* dir : {
static_cast<const double*
>(
nullptr),
direction.data()}) {
6478 std::array<double, 3> accelerated{0., 0., 0.};
6479 std::array<double, 3>
reference{0., 0., 0.};
6480 solid.ComputeNormal(point.data(), dir, accelerated.data());
6481 solid.ComputeNormal_Loop(point.data(), dir,
reference.data());
6482 for (
int dimension = 0; dimension < 3; ++dimension) {
6483 if (accelerated[dimension] !=
reference[dimension]) {
6488 return disagreements;
6493std::unique_ptr<SurfaceSolid> makeManyPatchSolid(
const char*
name,
int ringCount)
6495 auto solid = std::make_unique<SurfaceSolid>(
name);
6496 for (
int ring = 0; ring < ringCount; ++ring) {
6498 addBoxSurfaces(*solid, 0.4, 0.4, 0.4, {3. * std::cos(
angle), 3. * std::sin(
angle), 0.});
6500 solid->CloseShape();
6506BoundaryCurve quarterCircleBoundary(
double cu,
double cv,
double r,
double a0)
6509 const double aMid = 0.5 * (a0 + a1);
6510 std::vector<Point2D> poles{{cu +
r * std::cos(a0), cv +
r * std::sin(a0)},
6511 {cu +
r * std::sqrt(2.) * std::cos(aMid), cv +
r * std::sqrt(2.) * std::sin(aMid)},
6512 {cu +
r * std::cos(a1), cv +
r * std::sin(a1)}};
6513 return BoundaryCurve::makeBSpline(2, std::move(poles), {1., std::sqrt(0.5), 1.}, {0., 0., 0., 1., 1., 1.});
6520std::unique_ptr<SurfaceSolid> makeWireTrimmedSolid(
const char*
name)
6522 auto solid = std::make_unique<SurfaceSolid>(
name);
6523 for (
int window = 0; window < 4; ++window) {
6525 BOOST_REQUIRE(solid->AddCylindricalSurface(
6526 {0., 0., 0.}, {0., 0., 1.}, {1., 0., 0.}, 2., -1., 1., 0.,
surf::kTwoPi,
false,
6527 {quarterCircleBoundary(centrePhi, 0., 0.5, 0.), quarterCircleBoundary(centrePhi, 0., 0.5, surf::kHalfPi),
6528 quarterCircleBoundary(centrePhi, 0., 0.5, surf::kPi),
6529 quarterCircleBoundary(centrePhi, 0., 0.5, 3. * surf::kHalfPi)}));
6531 BOOST_REQUIRE(addDiskSurface(*solid, {0., 0., 1.}, {1., 0., 0.}, {0., 1., 0.}, 2.));
6532 BOOST_REQUIRE(addDiskSurface(*solid, {0., 0., -1.}, {1., 0., 0.}, {0., -1., 0.}, 2.));
6533 solid->CloseShape();
6537struct NearestPatchFixture {
6538 std::unique_ptr<SurfaceSolid> solid;
6542std::vector<NearestPatchFixture> nearestPatchFixtures()
6544 std::vector<NearestPatchFixture> fixtures;
6545 fixtures.push_back({makeBoxSolid(
"safetyBox", 1., 2., 3.), 4.});
6546 fixtures.push_back({makeTubeSolid(
"safetyTube", 0., 2., 3.), 4.});
6547 fixtures.push_back({makeTubeSolid(
"safetyHollowTube", 1., 2., 3.), 4.});
6548 fixtures.push_back({makeConeSolid(
"safetyCone", 2., 1., 3.), 4.});
6549 fixtures.push_back({makeSphereSolid(
"safetySphere", 2.5), 3.5});
6550 fixtures.push_back({makeTorusSolid(
"safetyTorus", 3., 1.), 4.5});
6551 fixtures.push_back({makeCapsuleSolid(
"safetyCapsule", 1., 1.5), 3.});
6552 fixtures.push_back({makeManyPatchSolid(
"safetyRing", 12), 4.5});
6553 fixtures.push_back({makeWireTrimmedSolid(
"safetyWireTrim"), 3.});
6580 const std::vector<Vec2> rectangle{{0., 0.}, {2., 0.}, {2., 3.}, {0., 3.}};
6581 BOOST_REQUIRE(polygon.
initialize({0., 0., 0.}, {1., 0., 0.}, {0., 1., 0.}, rectangle, {},
error));
6583 BOOST_REQUIRE(disk.
initialize({0., 0., 0.5}, {1., 0., 0.}, {0., 1., 0.},
6584 {surf::Curve2D::makeCircle({0., 0.}, 1.5)}, {},
error));
6589 BOOST_REQUIRE(splineFace.
initialize({0.2, -0.3, 1.1}, {1., 0., 0.}, {0., 1., 0.},
6590 {quarterCircleBSpline(0., 0., 1.2, 0.),
6591 quarterCircleBSpline(0., 0., 1.2, surf::kHalfPi),
6592 quarterCircleBSpline(0., 0., 1.2, surf::kPi),
6593 quarterCircleBSpline(0., 0., 1.2, 3. * surf::kHalfPi)},
6596 BOOST_REQUIRE(cylinder.
initialize({0.1, -0.2, 0.}, {0., 0., 1.}, {1., 0., 0.}, 2., -1., 1., 0.3,
6599 BOOST_REQUIRE(cone.
initialize({0., 0., 0.}, {0., 1., 0.}, {1., 0., 0.}, 2., 0.5, -1., 1., 0., 1.1 *
surf::kPi,
6602 BOOST_REQUIRE(sphere.
initialize({0.3, 0.4, -0.5}, {0., 0., 1.}, {1., 0., 0.}, 1.7, 0.2, 2.4, 0.,
6605 BOOST_REQUIRE(torus.
initialize({0., 0., 0.}, {0., 0., 1.}, {1., 0., 0.}, 3., 0.8, 0., 1.4 *
surf::kPi, 0.,
6608 const std::vector<const surf::BoundedSurface*>
surfaces{&polygon, &disk, &splineFace,
6609 &cylinder, &cone, &sphere,
6612 SampleStream
stream(0xB0B0Dull);
6615 Vec3
lower{TGeoShape::Big(), TGeoShape::Big(), TGeoShape::Big()};
6616 Vec3
upper{-TGeoShape::Big(), -TGeoShape::Big(), -TGeoShape::Big()};
6618 const double boxLower[3] = {
lower.xCoord,
lower.yCoord,
lower.zCoord};
6619 const double boxUpper[3] = {
upper.xCoord,
upper.yCoord,
upper.zCoord};
6621 for (
int sample = 0; sample < 4000; ++sample) {
6622 const double scale = (sample % 4 == 3) ? 1.e6 : ((sample % 4 == 2) ? 20. : 5.);
6623 const Vec3 point{
stream.symmetric(scale),
stream.symmetric(scale),
stream.symmetric(scale)};
6624 const double coordinates[3] = {point.xCoord, point.yCoord, point.zCoord};
6625 double boxDistanceSq = 0.;
6626 for (
int dimension = 0; dimension < 3; ++dimension) {
6627 if (coordinates[dimension] < boxLower[dimension]) {
6628 const double gap = boxLower[dimension] - coordinates[dimension];
6629 boxDistanceSq += gap * gap;
6630 }
else if (coordinates[dimension] > boxUpper[dimension]) {
6631 const double gap = coordinates[dimension] - boxUpper[dimension];
6632 boxDistanceSq += gap * gap;
6635 const double patchDistanceSq =
surface->distanceSqToPatch(point);
6639 BOOST_REQUIRE_LE(boxDistanceSq * (1. - 1.e-12), patchDistanceSq);
6649 size_t comparedPoints = 0;
6650 double worstSafetyGap = -std::numeric_limits<double>::infinity();
6651 for (
const auto& fixture : nearestPatchFixtures()) {
6652 BOOST_TEST_CONTEXT(
"fixture = " << fixture.solid->GetName())
6654 BOOST_REQUIRE(fixture.solid->HasBVH());
6655 for (
const auto& point : nearestPatchSample(*fixture.solid, fixture.extent, 2000)) {
6656 BOOST_TEST_CONTEXT(
"point = (" << point[0] <<
", " << point[1] <<
", " << point[2] <<
")")
6658 BOOST_REQUIRE_EQUAL(countNearestPatchDisagreements(*fixture.solid, point, &worstSafetyGap), 0);
6674 SurfaceSolid open(
"safetyBeforeClose");
6675 addBoxSurfaces(open, 1., 2., 3.);
6676 BOOST_REQUIRE(!open.HasBVH());
6677 const std::array<double, 3> probe{0.3, -1.1, 2.2};
6678 BOOST_CHECK_EQUAL(open.Safety(probe.data(), kTRUE), open.Safety_Loop(probe.data(), kTRUE));
6679 std::array<double, 3> viaBVH{0., 0., 0.};
6680 std::array<double, 3> viaLoop{0., 0., 0.};
6681 open.ComputeNormal(probe.data(),
nullptr, viaBVH.data());
6682 open.ComputeNormal_Loop(probe.data(),
nullptr, viaLoop.data());
6687 SurfaceSolid
empty(
"safetyEmpty");
6691 empty.ComputeNormal(probe.data(),
nullptr, viaBVH.data());
6692 empty.ComputeNormal_Loop(probe.data(),
nullptr, viaLoop.data());
6707 const auto box = makeBoxSolid(
"tieBreakBox", 2., 2., 2.);
6708 const std::array<double, 3> centre{0., 0., 0.};
6709 std::array<double, 3> viaBVH{0., 0., 0.};
6710 std::array<double, 3> viaLoop{0., 0., 0.};
6711 box->ComputeNormal(centre.data(),
nullptr, viaBVH.data());
6712 box->ComputeNormal_Loop(centre.data(),
nullptr, viaLoop.data());
6719 BOOST_CHECK_EQUAL(std::abs(viaLoop[0]) + std::abs(viaLoop[1]) + std::abs(viaLoop[2]), 1.);
6722 const auto tube = makeTubeSolid(
"tieBreakTube", 1., 2., 1.);
6723 tube->ComputeNormal(centre.data(),
nullptr, viaBVH.data());
6724 tube->ComputeNormal_Loop(centre.data(),
nullptr, viaLoop.data());
6742 BOOST_REQUIRE(!SurfaceSolid::GetSafetyBoundUnsoundForTest());
6744 const auto fixtures = nearestPatchFixtures();
6745 size_t caughtOnFixtures = 0;
6746 size_t prunableFixtures = 0;
6747 size_t sabotagedDisagreements = 0;
6748 size_t safetyTooLarge = 0;
6749 size_t safetyTooSmall = 0;
6750 for (
const auto& fixture : fixtures) {
6756 const auto points = nearestPatchSample(*fixture.solid, fixture.extent, 400);
6757 size_t disagreementsHere = 0;
6758 SurfaceSolid::SetSafetyBoundUnsoundForTest(
true);
6759 for (
const auto& point : points) {
6760 const double sabotaged = fixture.solid->Safety(point.data(), kTRUE);
6761 SurfaceSolid::SetSafetyBoundUnsoundForTest(
false);
6762 const double reference = fixture.solid->Safety_Loop(point.data(), kTRUE);
6763 SurfaceSolid::SetSafetyBoundUnsoundForTest(
true);
6765 ++disagreementsHere;
6766 (sabotaged >
reference) ? ++safetyTooLarge : ++safetyTooSmall;
6768 disagreementsHere +=
static_cast<size_t>(countNearestPatchDisagreements(*fixture.solid, point));
6770 SurfaceSolid::SetSafetyBoundUnsoundForTest(
false);
6771 sabotagedDisagreements += disagreementsHere;
6772 if (disagreementsHere > 0) {
6774 BOOST_TEST_MESSAGE(
"sabotaged bound caught on " << fixture.solid->GetName() <<
": " << disagreementsHere
6775 <<
" disagreements over " << points.size() <<
" points");
6777 BOOST_TEST_MESSAGE(
"sabotaged bound NOT caught on " << fixture.solid->GetName() <<
" ("
6778 << fixture.solid->GetNsurfaces() <<
" patches)");
6785 BOOST_CHECK_GT(sabotagedDisagreements, 100u);
6787 BOOST_CHECK_GT(safetyTooLarge, 0u);
6792 for (
const auto& fixture : fixtures) {
6793 for (
const auto& point : nearestPatchSample(*fixture.solid, fixture.extent, 400)) {
6794 BOOST_REQUIRE_EQUAL(countNearestPatchDisagreements(*fixture.solid, point), 0);
6797 BOOST_CHECK(!SurfaceSolid::GetSafetyBoundUnsoundForTest());
6806 const auto ring = makeManyPatchSolid(
"candidateRing", 24);
6807 BOOST_REQUIRE_EQUAL(ring->GetNsurfaces(), 144);
6808 const auto points = nearestPatchSample(*ring, 5., 500);
6810 SurfaceSolid::ResetSafetyCandidateCounter();
6811 for (
const auto& point : points) {
6812 ring->Safety(point.data(), kTRUE);
6814 const long long acceleratedCandidates = SurfaceSolid::GetSafetyCandidateCount();
6815 const double perCall =
static_cast<double>(acceleratedCandidates) / points.size();
6817 BOOST_TEST_MESSAGE(
"Safety candidates per call: " << perCall <<
" of " << ring->GetNsurfaces() <<
" patches");
6818 BOOST_CHECK_GT(acceleratedCandidates, 0);
6819 BOOST_CHECK_LT(perCall, 0.4 * ring->GetNsurfaces());
6822 SurfaceSolid::ResetSafetyCandidateCounter();
6823 for (
const auto& point : points) {
6824 ring->Safety_Loop(point.data(), kTRUE);
6825 std::array<double, 3> normal{0., 0., 0.};
6826 ring->ComputeNormal_Loop(point.data(),
nullptr, normal.data());
6832 SurfaceSolid::ResetSafetyCandidateCounter();
6833 for (
const auto& point : points) {
6834 std::array<double, 3> normal{0., 0., 0.};
6835 ring->ComputeNormal(point.data(),
nullptr, normal.data());
6837 const double normalPerCall =
static_cast<double>(SurfaceSolid::GetSafetyCandidateCount()) / points.size();
6838 BOOST_TEST_MESSAGE(
"ComputeNormal candidates per call: " << normalPerCall);
6839 BOOST_CHECK_LT(normalPerCall, 0.4 * ring->GetNsurfaces());
6840 BOOST_CHECK_GE(normalPerCall, perCall);
6855using CoverBox = surf::BoundedSurface::CoverBox;
6859double coverBoxDistanceSq(
const CoverBox&
box,
const surf::Vec3& point)
6863 const double lower[3] = {
box.first.xCoord,
box.first.yCoord,
box.first.zCoord};
6864 const double upper[3] = {
box.second.xCoord,
box.second.yCoord,
box.second.zCoord};
6865 for (
int dimension = 0; dimension < 3; ++dimension) {
6866 const double gap = std::max({
lower[dimension] - coordinates[dimension],
6867 coordinates[dimension] -
upper[dimension], 0.});
6873double minCoverBoxDistanceSq(
const std::vector<CoverBox>& boxes,
const surf::Vec3& point)
6875 double best = std::numeric_limits<double>::infinity();
6876 for (
const auto&
box : boxes) {
6877 best = std::min(best, coverBoxDistanceSq(
box, point));
6882bool anyCoverBoxContains(
const std::vector<CoverBox>& boxes,
const surf::Vec3& point,
double slack)
6884 for (
const auto&
box : boxes) {
6885 if (point.
xCoord >=
box.first.xCoord - slack && point.
xCoord <=
box.second.xCoord + slack &&
6886 point.
yCoord >=
box.first.yCoord - slack && point.
yCoord <=
box.second.yCoord + slack &&
6887 point.
zCoord >=
box.first.zCoord - slack && point.
zCoord <=
box.second.zCoord + slack) {
6901 std::vector<CoverBox> boxes;
6902 surface.appendCoverBoxes(boxes);
6903 BOOST_TEST_CONTEXT(
"surface = " <<
label)
6905 BOOST_REQUIRE(!boxes.empty());
6906 for (
const auto& point : patchPoints) {
6907 BOOST_TEST_CONTEXT(
"patch point = (" << point.
xCoord <<
", " << point.
yCoord <<
", " << point.
zCoord <<
")")
6909 BOOST_CHECK(anyCoverBoxContains(boxes, point, 1.e-9));
6912 SampleStream
stream(0xC0FEB0C5ull);
6915 constexpr double kProbeScales[3] = {1.5, 8., 300.};
6917 const double scale = kProbeScales[
index % 3];
6919 const double patchDistanceSq =
surface.distanceSqToPatch(point);
6920 const double boxDistanceSq = minCoverBoxDistanceSq(boxes, point);
6921 BOOST_TEST_CONTEXT(
"point = (" << point.
xCoord <<
", " << point.
yCoord <<
", " << point.
zCoord <<
")")
6923 BOOST_CHECK_LE(boxDistanceSq, patchDistanceSq * (1. + 1.e-9) + 1.e-18);
6933 std::vector<CoverBox> boxes;
6934 surface.appendCoverBoxes(boxes);
6935 BOOST_TEST_CONTEXT(
"surface = " <<
label)
6937 BOOST_CHECK_GT(boxes.size(), 1u);
6951 const Vec3 skewCenter{0.4, -0.2, 0.1};
6952 const Vec3 skewAxis{0.2, 0.3, 1.};
6953 const Vec3 referenceU{1., 0., 0.};
6957 BOOST_REQUIRE(cylinder.
initialize(skewCenter, skewAxis, referenceU, 1.7, -0.8, 1.2, 0.4, 1.9,
false,
error));
6958 std::vector<Vec3> patchPoints;
6959 for (
int stepPhi = 0; stepPhi <= 12; ++stepPhi) {
6960 for (
int stepH = 0; stepH <= 4; ++stepH) {
6961 patchPoints.push_back(cylinder.
pointAt(0.4 + 1.9 * stepPhi / 12., -0.8 + 2. * stepH / 4.));
6964 checkCoverBoxProperties(cylinder, patchPoints,
"partial cylinder, skew axis");
6965 checkEmitsSeveralCoverBoxes(cylinder,
"partial cylinder, skew axis");
6969 BOOST_REQUIRE(cylinder.
initialize({0., 0., 0.}, {0., 0., 1.}, referenceU, 2., -1., 1., 0.,
surf::kTwoPi,
6971 std::vector<Vec3> patchPoints;
6972 for (
int stepPhi = 0; stepPhi <= 24; ++stepPhi) {
6973 patchPoints.push_back(cylinder.
pointAt(
surf::kTwoPi * stepPhi / 24., -1. + 2. * (stepPhi % 5) / 4.));
6975 checkCoverBoxProperties(cylinder, patchPoints,
"full cylinder");
6976 checkEmitsSeveralCoverBoxes(cylinder,
"full cylinder");
6980 BOOST_REQUIRE(cylinder.
initialize({0., 0., 0.}, {0., 0., 1.}, referenceU, 2., -1., 1., 0.,
surf::kTwoPi,
false,
6981 paramRectWireCurves(0.3, 2.1, -0.5, 0.7), {},
error));
6982 std::vector<Vec3> patchPoints;
6983 for (
int stepPhi = 0; stepPhi <= 10; ++stepPhi) {
6984 for (
int stepH = 0; stepH <= 4; ++stepH) {
6985 const double phi = 0.3 + 1.8 * stepPhi / 10.;
6986 const double height = -0.5 + 1.2 * stepH / 4.;
6992 BOOST_REQUIRE(!patchPoints.empty());
6993 checkCoverBoxProperties(cylinder, patchPoints,
"wire-trimmed cylinder");
6999 std::vector<Vec3> patchPoints;
7000 for (
int stepTheta = 0; stepTheta <= 8; ++stepTheta) {
7001 for (
int stepPhi = 0; stepPhi < 16; ++stepPhi) {
7005 checkCoverBoxProperties(sphere, patchPoints,
"full sphere, skew frame");
7006 checkEmitsSeveralCoverBoxes(sphere,
"full sphere, skew frame");
7012 BOOST_REQUIRE(cap.
initialize({0., 0., 0.}, {0., 0., 1.}, referenceU, 2., 0., 0.6, 0.2, 1.1,
false,
error));
7013 std::vector<Vec3> patchPoints;
7014 for (
int stepTheta = 0; stepTheta <= 4; ++stepTheta) {
7015 for (
int stepPhi = 0; stepPhi <= 6; ++stepPhi) {
7016 patchPoints.push_back(cap.
pointAt(0.6 * stepTheta / 4., 0.2 + 1.1 * stepPhi / 6.));
7019 checkCoverBoxProperties(cap, patchPoints,
"spherical cap");
7023 BOOST_REQUIRE(cone.
initialize(skewCenter, skewAxis, referenceU, 2., 0.5, -0.9, 1.1, 0.7, 2.3,
false,
error));
7024 std::vector<Vec3> patchPoints;
7025 for (
int stepPhi = 0; stepPhi <= 10; ++stepPhi) {
7026 for (
int stepH = 0; stepH <= 4; ++stepH) {
7027 patchPoints.push_back(cone.
pointAt(0.7 + 2.3 * stepPhi / 10., -0.9 + 2. * stepH / 4.));
7030 checkCoverBoxProperties(cone, patchPoints,
"partial cone, skew axis");
7031 checkEmitsSeveralCoverBoxes(cone,
"partial cone, skew axis");
7035 BOOST_REQUIRE(torus.
initialize(skewCenter, skewAxis, referenceU, 2.4, 0.7, 0.3, 2.1, -0.4, 1.7,
false,
error));
7036 std::vector<Vec3> patchPoints;
7037 for (
int stepRing = 0; stepRing <= 10; ++stepRing) {
7038 for (
int stepTube = 0; stepTube <= 6; ++stepTube) {
7039 patchPoints.push_back(torus.
pointAt(0.3 + 2.1 * stepRing / 10., -0.4 + 1.7 * stepTube / 6.));
7042 checkCoverBoxProperties(torus, patchPoints,
"partial torus, skew axis");
7043 checkEmitsSeveralCoverBoxes(torus,
"partial torus, skew axis");
7047 BOOST_REQUIRE(torus.
initialize({0., 0., 0.}, {0., 0., 1.}, referenceU, 3., 1., 0.,
surf::kTwoPi, 0.,
7049 std::vector<Vec3> patchPoints;
7050 for (
int stepRing = 0; stepRing < 16; ++stepRing) {
7051 for (
int stepTube = 0; stepTube < 8; ++stepTube) {
7055 checkCoverBoxProperties(torus, patchPoints,
"full torus");
7059 const std::vector<surf::Vec2> rectangle{{0., 0.}, {2., 0.}, {2., 3.}, {0., 3.}};
7060 BOOST_REQUIRE(polygon.
initialize({0.2, -0.4, 0.5}, {1., 0.2, 0.}, {-0.1, 1., 0.3}, rectangle, {},
error));
7061 std::vector<Vec3> patchPoints;
7062 patchPoints.push_back(polygon.
toGlobal({0.01, 0.01}));
7063 patchPoints.push_back(polygon.
toGlobal({1.9, 2.9}));
7064 checkCoverBoxProperties(polygon, patchPoints,
"planar polygon");
7075 const auto sphere = makeSphereSolid(
"subBoxSphere", 2.5);
7076 BOOST_CHECK_EQUAL(sphere->CountBVHRayCandidates({2.2, 2.2, -5.}, {0., 0., 1.}), 0);
7077 BOOST_CHECK_GE(sphere->CountBVHRayCandidates({0., 0., -5.}, {0., 0., 1.}), 1);
7080 const auto tube = makeTubeSolid(
"subBoxTube", 0., 2., 1.);
7084 const auto torus = makeTorusSolid(
"subBoxTorus", 3., 1.);
7085 BOOST_CHECK_EQUAL(torus->CountBVHRayCandidates({3.4, 3.4, -5.}, {0., 0., 1.}), 0);
7086 BOOST_CHECK_GE(torus->CountBVHRayCandidates({3., 0., -5.}, {0., 0., 1.}), 1);
7090 SurfaceSolid quarter(
"subBoxQuarterCylinder");
7091 BOOST_REQUIRE(quarter.AddCylindricalSurface({0., 0., 0.}, {0., 0., 1.}, {1., 0., 0.}, 2., -1., 1.,
7093 quarter.CloseShape(
false);
7094 BOOST_REQUIRE(quarter.HasBVH());
7095 BOOST_CHECK_EQUAL(quarter.CountBVHRayCandidates({-1.9, -5., 0.}, {0., 1., 0.}), 0);
7096 BOOST_CHECK_GE(quarter.CountBVHRayCandidates({5., 0., 0.}, {-1., 0., 0.}), 1);
7107 std::unique_ptr<SurfaceSolid> solid;
7110 std::vector<Fixture> fixtures;
7111 fixtures.push_back({makeSphereSolid(
"subBoxSweepSphere", 2.5), 3.5});
7112 fixtures.push_back({makeTorusSolid(
"subBoxSweepTorus", 3., 1.), 4.5});
7113 fixtures.push_back({makeCapsuleSolid(
"subBoxSweepCapsule", 1., 1.5), 3.});
7114 fixtures.push_back({makeConeSolid(
"subBoxSweepCone", 2., 1., 3.), 4.});
7115 fixtures.push_back({makeWireTrimmedSolid(
"subBoxSweepWireTrim"), 3.});
7117 for (
const auto& fixture : fixtures) {
7118 BOOST_TEST_CONTEXT(
"fixture = " << fixture.solid->GetName())
7120 sweepDistanceAgainstLoop(*fixture.solid, fixture.extent, 4);
7121 for (
const auto& point : probeGrid(fixture.extent, 4)) {
7122 BOOST_TEST_CONTEXT(
"point = (" << point[0] <<
", " << point[1] <<
", " << point[2] <<
")")
7124 BOOST_CHECK_EQUAL(fixture.solid->Contains(point.data()), fixture.solid->Contains_Loop(point.data()));
7125 std::vector<SurfaceSolid::ContainsCrossing> bvhCrossings;
7126 std::vector<SurfaceSolid::ContainsCrossing> loopCrossings;
7127 fixture.solid->DescribeContainsCrossings({point[0], point[1], point[2]}, bvhCrossings, loopCrossings);
7128 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
configStandalone Light OpenGL display display_light Apply smoothing to points Enable Z buffer Highlight links Highlight seeds Highlight sector tracks Highlight final tracks Show individual Show related Highlight grid Show only individual collision Color clusters belonging to track objects Set line Enable drawing TPC data Enable drawing TOF data Enable drawing magnetic field Set field line step count Invert colors configStandalone Heavy OpenGL display display_heavy Show tracks with ITS contribution Show tracks with TOF contribution Enable propagation of loopers Show clusters at nominal x of pad row for early transformed data Split CE tracks when they cross the central Mark clusters Mark adjacent clusters Hide rejected tracks Apply filter on tracks to be displayed Separate increase X h configStandalone Camera window renderer OpenGL display display_camera Orient the camera such that the y axis is always upwards Full Screen MultiSample Anti Downsampling Anti Enable Vertical Sync Target frame or if the backend cannot limit to the refresh Use OpenGL indirect draws to reduce number of draw calls Display h configStandalone Vulkan display display_vulkan Have uniform buffers in host accessible device memory h configStandalone Internal display display_internal Number of updates the Qt GUI may send ahead of the last one the display has before it waits for the display Opacity of the console background Number of console lines the history is shortened to this when it reaches twice the size Time in seconds an info message fades out after setInfoTime h OpenGL display display Font size Do not use Freetype for font Renderer for event Frontend to use for event ROOT macros used as track filter Console script file
configStandalone Light OpenGL display display_light Apply smoothing to points Enable Z buffer Highlight links Highlight seeds Highlight sector tracks Highlight final tracks Show individual Show related Highlight grid Show only individual collision Color clusters belonging to track objects Set line Enable drawing TPC data Enable drawing TOF data Enable drawing magnetic field Set field line step count Invert colors configStandalone Heavy OpenGL display display_heavy Show tracks with ITS contribution Show tracks with TOF contribution Enable propagation of loopers Show clusters at nominal x of pad row for early transformed data Split CE tracks when they cross the central Mark clusters Mark adjacent clusters Hide rejected tracks Apply filter on tracks to be displayed Separate increase X coordinate
Timeframe settings Average distance in cm of events merged into Randomize distance around average distance of merged events Simulate correct bunch interactions instead of placing only the average number of events A value[n] sets TFlen for[n] collisions in Number of bunches per trainsort Number of bunch trains Instantaneous interaction rate Do not simulate border Total number of collisions to be placed in the interior of all time Do not select random but walk over array of events in stride steps h Event generator settings h configStandalone Use CPU for GPU Force usage of the specified GPU device no CPU fallback Continue processing after an error First event to Number of events to Number of iterations to Number of iterations to Number of initial iterations excluded from const char e
Timeframe settings Average distance in cm of events merged into Randomize distance around average distance of merged events Simulate correct bunch interactions instead of placing only the average number of events A value[n] sets TFlen for[n] collisions in Number of bunches per trainsort Number of bunch trains Instantaneous interaction rate Do not simulate border Total number of collisions to be placed in the interior of all time Do not select random but walk over array of events in stride steps h Event generator settings h configStandalone Use CPU for GPU Force usage of the specified GPU device no CPU fallback Continue processing after an error First event to Number of events to Number of iterations to Number of iterations to Number of initial iterations excluded from const char Directory with events to Events directory is and not inside events d
configStandalone Reconstruction rec_tpc multiplier for the combined cluster track error during track following area in cm for the search of for z only used if decrease removal by factor
configStandalone Reconstruction rec_tpc multiplier for the combined cluster track error during track following area in cm for the search of for z only used if decrease removal by when above this local occupancy rowx Shift of TPC revert cluster flag edge bit distance to edge exceeds this parameter(fixed during CTF decoding)") AddOptionRTC(defaultZOffsetOverR
char error[CURL_ERROR_SIZE]
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...
Acts::RotationMatrix3 rotation
const Acts::Surface * surface
std::vector< std::shared_ptr< Surface > > surfaces
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 appendDisplayMesh(std::vector< Vec3 > &vertices, std::vector< std::array< int, 3 > > &triangles) const =0
Append this patch's visualization triangulation (navigation must never depend on it).
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).
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...
double normSq(const Vec3 &vector)
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
double norm(const Vec3 &vector)
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)
int const float const TrackSeed< NLayers > & seed
const bool const bool const int FollowDirection BestTrial TrackITSInternal< NLayers > & best
const float const float const float const float std::array< float, 3 > std::array< float, 3 > & direction
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())