Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

feat: Add closest point and distance computation to SurfaceBounds #3990

Open
wants to merge 44 commits into
base: main
Choose a base branch
from
Open
Changes from 1 commit
Commits
Show all changes
44 commits
Select commit Hold shift + click to select a range
87106e9
tmp
andiwand Dec 12, 2024
4d5d10d
tmp
andiwand Dec 13, 2024
9bcc7f6
tmp
andiwand Dec 17, 2024
e6e0cc2
draft
andiwand Dec 17, 2024
bdf87b6
clean annulus
andiwand Dec 18, 2024
f856077
tmp
andiwand Dec 18, 2024
7f29993
tmp
andiwand Dec 18, 2024
a595467
tmp
andiwand Dec 18, 2024
793c950
includes
andiwand Dec 18, 2024
0c965e7
tmp
andiwand Dec 18, 2024
7a1c005
tmp
andiwand Dec 18, 2024
79ada20
make compile
andiwand Dec 18, 2024
17db365
tmp
andiwand Dec 18, 2024
289109a
jacobians
andiwand Dec 18, 2024
a034c42
Merge branch 'main' of github.com:acts-project/acts into tmp-distance…
andiwand Dec 18, 2024
edb6882
metric
andiwand Dec 19, 2024
a2bf8be
remove some unintended changes
andiwand Dec 19, 2024
43c7d23
default inside impl; fix compilation
andiwand Dec 19, 2024
26e68b3
Merge branch 'main' of github.com:acts-project/acts into tmp-distance…
andiwand Dec 19, 2024
7ea4115
deduplicate inside
andiwand Dec 19, 2024
890b472
remove bounds check helpers
andiwand Dec 19, 2024
c30f1cc
Merge branch 'main' of github.com:acts-project/acts into tmp-distance…
andiwand Dec 19, 2024
b9e2168
doc
andiwand Dec 19, 2024
cf9a545
rename tolerance mode
andiwand Dec 19, 2024
ef18a08
fix doc
andiwand Dec 19, 2024
c1de29d
rename tolerance mode
andiwand Dec 19, 2024
5032c8b
simplify distance
andiwand Dec 19, 2024
bf236b4
fix and docs
andiwand Dec 19, 2024
217b138
small fix
andiwand Dec 19, 2024
a0d665d
remove cartesian to bound jacobians
andiwand Dec 19, 2024
529a109
fix doc; tolerance less annulus inside check
andiwand Dec 19, 2024
97c568e
deal with different coordinate systems in annulus bounds
andiwand Dec 19, 2024
385c66d
try fix annulus closest point
andiwand Dec 19, 2024
cd13630
minor
andiwand Dec 19, 2024
c8be68f
minor
andiwand Dec 19, 2024
5727270
Merge branch 'main' of github.com:acts-project/acts into tmp-distance…
andiwand Dec 19, 2024
c6a4219
fixes and ref updates
andiwand Dec 19, 2024
c283860
Merge branch 'main' of github.com:acts-project/acts into tmp-distance…
andiwand Jan 21, 2025
d1c2343
adapt annulus bounds test; yet to check if correct
andiwand Jan 21, 2025
54bf281
Merge branch 'main' of github.com:acts-project/acts into tmp-distance…
andiwand Jan 21, 2025
9e80dbb
one more case
andiwand Jan 21, 2025
0e36076
after debugginRg
andiwand Jan 22, 2025
72f2fb5
revert
andiwand Jan 22, 2025
45d5f8d
update refs
andiwand Jan 22, 2025
File filter

Filter by extension

Filter by extension

Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
Prev Previous commit
Next Next commit
tmp
andiwand committed Dec 18, 2024
commit 17db3653816b71179a6fae51e24e176d6c5e5e83
16 changes: 12 additions & 4 deletions Core/include/Acts/Surfaces/ConeBounds.hpp
Original file line number Diff line number Diff line change
@@ -74,13 +74,19 @@ class ConeBounds : public SurfaceBounds {
/// @param values The parameter array
ConeBounds(const std::array<double, eSize>& values) noexcept(false);

BoundsType type() const final;
BoundsType type() const final { return eCone; }

bool isCartesian() const final { return false; }
bool isCartesian() const final { return true; }

SquareMatrix2 boundToCartesianJacobian(const Vector2& lposition) const final;
SquareMatrix2 boundToCartesianJacobian(const Vector2& lposition) const final {
(void)lposition;
return SquareMatrix2::Identity();
}

SquareMatrix2 cartesianToBoundJacobian(const Vector2& lposition) const final;
SquareMatrix2 cartesianToBoundJacobian(const Vector2& lposition) const final {
(void)lposition;
return SquareMatrix2::Identity();
}

/// Return the bound values as dynamically sized vector
///
@@ -129,6 +135,8 @@ class ConeBounds : public SurfaceBounds {

/// Private helper function to shift a local 2D position
///
/// Shift r-phi coordinate to be centered around the average phi.
///
/// @param lposition The original local position
Vector2 shifted(const Vector2& lposition) const;
};
5 changes: 1 addition & 4 deletions Core/include/Acts/Surfaces/ConeSurface.hpp
Original file line number Diff line number Diff line change
@@ -11,7 +11,6 @@
#include "Acts/Definitions/Algebra.hpp"
#include "Acts/Definitions/Alignment.hpp"
#include "Acts/Definitions/Tolerance.hpp"
#include "Acts/Definitions/TrackParametrization.hpp"
#include "Acts/Geometry/GeometryContext.hpp"
#include "Acts/Geometry/Polyhedron.hpp"
#include "Acts/Surfaces/BoundaryTolerance.hpp"
@@ -23,8 +22,6 @@
#include "Acts/Utilities/Result.hpp"
#include "Acts/Utilities/detail/RealQuadraticEquation.hpp"

#include <cmath>
#include <cstddef>
#include <memory>
#include <numbers>
#include <string>
@@ -41,7 +38,7 @@ namespace Acts {
/// at the tip of the cone.
/// Propagations to a cone surface will be returned in
/// curvilinear coordinates.

///
class ConeSurface : public RegularSurface {
friend class Surface;

12 changes: 9 additions & 3 deletions Core/include/Acts/Surfaces/CylinderBounds.hpp
Original file line number Diff line number Diff line change
@@ -86,11 +86,17 @@ class CylinderBounds : public SurfaceBounds {

BoundsType type() const final { return eCylinder; }

bool isCartesian() const final { return false; }
bool isCartesian() const final { return true; }

SquareMatrix2 boundToCartesianJacobian(const Vector2& lposition) const final;
SquareMatrix2 boundToCartesianJacobian(const Vector2& lposition) const final {
(void)lposition;
return SquareMatrix2::Identity();
}

SquareMatrix2 cartesianToBoundJacobian(const Vector2& lposition) const final;
SquareMatrix2 cartesianToBoundJacobian(const Vector2& lposition) const final {
(void)lposition;
return SquareMatrix2::Identity();
}

/// Return the bound values as dynamically sized vector
///
3 changes: 0 additions & 3 deletions Core/include/Acts/Surfaces/CylinderSurface.hpp
Original file line number Diff line number Diff line change
@@ -11,7 +11,6 @@
#include "Acts/Definitions/Algebra.hpp"
#include "Acts/Definitions/Alignment.hpp"
#include "Acts/Definitions/Tolerance.hpp"
#include "Acts/Definitions/TrackParametrization.hpp"
#include "Acts/Geometry/GeometryContext.hpp"
#include "Acts/Geometry/Polyhedron.hpp"
#include "Acts/Surfaces/BoundaryTolerance.hpp"
@@ -24,8 +23,6 @@
#include "Acts/Utilities/Result.hpp"
#include "Acts/Utilities/detail/RealQuadraticEquation.hpp"

#include <cmath>
#include <cstddef>
#include <memory>
#include <numbers>
#include <string>
16 changes: 12 additions & 4 deletions Core/src/Surfaces/AnnulusBounds.cpp
Original file line number Diff line number Diff line change
@@ -160,14 +160,22 @@ std::vector<Vector2> AnnulusBounds::vertices(

SquareMatrix2 AnnulusBounds::boundToCartesianJacobian(
const Vector2& lposition) const {
(void)lposition;
return SquareMatrix2::Identity(); // TODO
SquareMatrix2 j;
j(0, 0) = std::cos(lposition[0]);
j(0, 1) = -lposition[1] * std::sin(lposition[0]);
j(1, 0) = std::sin(lposition[0]);
j(1, 1) = lposition[1] * std::cos(lposition[0]);
return j;
}

SquareMatrix2 AnnulusBounds::cartesianToBoundJacobian(
const Vector2& lposition) const {
(void)lposition;
return SquareMatrix2::Identity(); // TODO
SquareMatrix2 j;
j(0, 0) = std::cos(lposition[0]);
j(0, 1) = std::sin(lposition[0]);
j(1, 0) = -std::sin(lposition[0]) / lposition[1];
j(1, 1) = std::cos(lposition[0]) / lposition[1];
return j;
}

bool AnnulusBounds::inside(const Vector2& lposition) const {
17 changes: 0 additions & 17 deletions Core/src/Surfaces/ConeBounds.cpp
Original file line number Diff line number Diff line change
@@ -41,23 +41,6 @@ ConeBounds::ConeBounds(const std::array<double, eSize>& values) noexcept(false)
checkConsistency();
}

SurfaceBounds::BoundsType ConeBounds::type() const {
return SurfaceBounds::eCone;
}

SquareMatrix2 ConeBounds::boundToCartesianJacobian(
const Vector2& lposition) const {
(void)lposition;
return SquareMatrix2::Identity(); // TODO
}

SquareMatrix2 ConeBounds::cartesianToBoundJacobian(
const Vector2& lposition) const {
(void)lposition;
return SquareMatrix2::Identity(); // TODO
}

/// Shift r-phi coordinate to be centered around the average phi.
Vector2 ConeBounds::shifted(const Vector2& lposition) const {
using detail::radian_sym;

1 change: 0 additions & 1 deletion Core/src/Surfaces/ConeSurface.cpp
Original file line number Diff line number Diff line change
@@ -14,7 +14,6 @@
#include "Acts/Surfaces/detail/AlignmentHelper.hpp"
#include "Acts/Surfaces/detail/FacesHelper.hpp"
#include "Acts/Surfaces/detail/VerticesHelper.hpp"
#include "Acts/Utilities/Helpers.hpp"
#include "Acts/Utilities/Intersection.hpp"
#include "Acts/Utilities/ThrowAssert.hpp"
#include "Acts/Utilities/detail/RealQuadraticEquation.hpp"
26 changes: 2 additions & 24 deletions Core/src/Surfaces/CylinderBounds.cpp
Original file line number Diff line number Diff line change
@@ -9,7 +9,6 @@
#include "Acts/Surfaces/CylinderBounds.hpp"

#include "Acts/Definitions/Algebra.hpp"
#include "Acts/Definitions/TrackParametrization.hpp"
#include "Acts/Surfaces/BoundaryTolerance.hpp"
#include "Acts/Surfaces/detail/BoundaryCheckHelper.hpp"
#include "Acts/Surfaces/detail/VerticesHelper.hpp"
@@ -28,30 +27,9 @@ namespace Acts {
using VectorHelpers::perp;
using VectorHelpers::phi;

SquareMatrix2 CylinderBounds::boundToCartesianJacobian(
const Vector2& /*lposition*/) const {
SquareMatrix2 j;
j(0, eBoundLoc0) = get(eR);
j(0, eBoundLoc1) = 0;
j(1, eBoundLoc0) = 0;
j(1, eBoundLoc1) = 1;
return j;
}

SquareMatrix2 CylinderBounds::cartesianToBoundJacobian(
const Vector2& /*lposition*/) const {
SquareMatrix2 j;
j(0, eBoundLoc0) = 1 / get(eR);
j(0, eBoundLoc1) = 0;
j(1, eBoundLoc0) = 0;
j(1, eBoundLoc1) = 1;
return j;
}

Vector2 CylinderBounds::shifted(const Vector2& lposition) const {
return {
detail::radian_sym((lposition[eBoundLoc0] / get(eR)) - get(eAveragePhi)),
lposition[eBoundLoc1]};
return {detail::radian_sym((lposition[0] / get(eR)) - get(eAveragePhi)),
lposition[1]};
}

bool CylinderBounds::inside(const Vector2& lposition) const {
110 changes: 51 additions & 59 deletions Core/src/Surfaces/CylinderSurface.cpp
Original file line number Diff line number Diff line change
@@ -19,7 +19,6 @@
#include "Acts/Surfaces/detail/FacesHelper.hpp"
#include "Acts/Surfaces/detail/MergeHelper.hpp"
#include "Acts/Utilities/BinningType.hpp"
#include "Acts/Utilities/Helpers.hpp"
#include "Acts/Utilities/Intersection.hpp"
#include "Acts/Utilities/ThrowAssert.hpp"
#include "Acts/Utilities/detail/periodic.hpp"
@@ -34,46 +33,41 @@
#include <vector>

namespace Acts {
class DetectorElementBase;
} // namespace Acts

using Acts::VectorHelpers::perp;
using Acts::VectorHelpers::phi;
using VectorHelpers::perp;
using VectorHelpers::phi;

Acts::CylinderSurface::CylinderSurface(const CylinderSurface& other)
CylinderSurface::CylinderSurface(const CylinderSurface& other)
: GeometryObject(), RegularSurface(other), m_bounds(other.m_bounds) {}

Acts::CylinderSurface::CylinderSurface(const GeometryContext& gctx,
const CylinderSurface& other,
const Transform3& shift)
CylinderSurface::CylinderSurface(const GeometryContext& gctx,
const CylinderSurface& other,
const Transform3& shift)
: GeometryObject(),
RegularSurface(gctx, other, shift),
m_bounds(other.m_bounds) {}

Acts::CylinderSurface::CylinderSurface(const Transform3& transform,
double radius, double halfz,
double halfphi, double avphi,
double bevelMinZ, double bevelMaxZ)
CylinderSurface::CylinderSurface(const Transform3& transform, double radius,
double halfz, double halfphi, double avphi,
double bevelMinZ, double bevelMaxZ)
: RegularSurface(transform),
m_bounds(std::make_shared<const CylinderBounds>(
radius, halfz, halfphi, avphi, bevelMinZ, bevelMaxZ)) {}

Acts::CylinderSurface::CylinderSurface(
std::shared_ptr<const CylinderBounds> cbounds,
const DetectorElementBase& detelement)
CylinderSurface::CylinderSurface(std::shared_ptr<const CylinderBounds> cbounds,
const DetectorElementBase& detelement)
: RegularSurface(detelement), m_bounds(std::move(cbounds)) {
/// surfaces representing a detector element must have bounds
throw_assert(m_bounds, "CylinderBounds must not be nullptr");
}

Acts::CylinderSurface::CylinderSurface(
const Transform3& transform, std::shared_ptr<const CylinderBounds> cbounds)
CylinderSurface::CylinderSurface(const Transform3& transform,
std::shared_ptr<const CylinderBounds> cbounds)
: RegularSurface(transform), m_bounds(std::move(cbounds)) {
throw_assert(m_bounds, "CylinderBounds must not be nullptr");
}

Acts::CylinderSurface& Acts::CylinderSurface::operator=(
const CylinderSurface& other) {
CylinderSurface& CylinderSurface::operator=(const CylinderSurface& other) {
if (this != &other) {
Surface::operator=(other);
m_bounds = other.m_bounds;
@@ -82,11 +76,10 @@ Acts::CylinderSurface& Acts::CylinderSurface::operator=(
}

// return the binning position for ordering in the BinnedArray
Acts::Vector3 Acts::CylinderSurface::binningPosition(
const GeometryContext& gctx, BinningValue bValue) const {
Vector3 CylinderSurface::binningPosition(const GeometryContext& gctx,
BinningValue bValue) const {
// special binning type for R-type methods
if (bValue == Acts::BinningValue::binR ||
bValue == Acts::BinningValue::binRPhi) {
if (bValue == BinningValue::binR || bValue == BinningValue::binRPhi) {
double R = bounds().get(CylinderBounds::eR);
double phi = bounds().get(CylinderBounds::eAveragePhi);
return localToGlobal(gctx, Vector2{phi * R, 0}, Vector3{});
@@ -99,7 +92,7 @@ Acts::Vector3 Acts::CylinderSurface::binningPosition(
}

// return the measurement frame: it's the tangential plane
Acts::RotationMatrix3 Acts::CylinderSurface::referenceFrame(
RotationMatrix3 CylinderSurface::referenceFrame(
const GeometryContext& gctx, const Vector3& position,
const Vector3& /*direction*/) const {
RotationMatrix3 mFrame;
@@ -118,22 +111,22 @@ Acts::RotationMatrix3 Acts::CylinderSurface::referenceFrame(
return mFrame;
}

Acts::Surface::SurfaceType Acts::CylinderSurface::type() const {
Surface::SurfaceType CylinderSurface::type() const {
return Surface::Cylinder;
}

Acts::Vector3 Acts::CylinderSurface::localToGlobal(
const GeometryContext& gctx, const Vector2& lposition) const {
Vector3 CylinderSurface::localToGlobal(const GeometryContext& gctx,
const Vector2& lposition) const {
// create the position in the local 3d frame
double r = bounds().get(CylinderBounds::eR);
double phi = lposition[Acts::eBoundLoc0] / r;
Vector3 position(r * cos(phi), r * sin(phi), lposition[Acts::eBoundLoc1]);
double phi = lposition[0] / r;
Vector3 position(r * cos(phi), r * sin(phi), lposition[1]);
return transform(gctx) * position;
}

Acts::Result<Acts::Vector2> Acts::CylinderSurface::globalToLocal(
const GeometryContext& gctx, const Vector3& position,
double tolerance) const {
Result<Vector2> CylinderSurface::globalToLocal(const GeometryContext& gctx,
const Vector3& position,
double tolerance) const {
double inttol = tolerance;
if (tolerance == s_onSurfaceTolerance) {
// transform default value!
@@ -153,19 +146,19 @@ Acts::Result<Acts::Vector2> Acts::CylinderSurface::globalToLocal(
{bounds().get(CylinderBounds::eR) * phi(loc3Dframe), loc3Dframe.z()});
}

std::string Acts::CylinderSurface::name() const {
return "Acts::CylinderSurface";
std::string CylinderSurface::name() const {
return "CylinderSurface";
}

Acts::Vector3 Acts::CylinderSurface::normal(
const GeometryContext& gctx, const Acts::Vector2& lposition) const {
double phi = lposition[Acts::eBoundLoc0] / m_bounds->get(CylinderBounds::eR);
Vector3 CylinderSurface::normal(const GeometryContext& gctx,
const Vector2& lposition) const {
double phi = lposition[eBoundLoc0] / m_bounds->get(CylinderBounds::eR);
Vector3 localNormal(cos(phi), sin(phi), 0.);
return transform(gctx).linear() * localNormal;
}

Acts::Vector3 Acts::CylinderSurface::normal(
const GeometryContext& gctx, const Acts::Vector3& position) const {
Vector3 CylinderSurface::normal(const GeometryContext& gctx,
const Vector3& position) const {
const Transform3& sfTransform = transform(gctx);
// get it into the cylinder frame
Vector3 pos3D = sfTransform.inverse() * position;
@@ -175,19 +168,19 @@ Acts::Vector3 Acts::CylinderSurface::normal(
return sfTransform.linear() * pos3D.normalized();
}

double Acts::CylinderSurface::pathCorrection(
const GeometryContext& gctx, const Acts::Vector3& position,
const Acts::Vector3& direction) const {
double CylinderSurface::pathCorrection(const GeometryContext& gctx,
const Vector3& position,
const Vector3& direction) const {
Vector3 normalT = normal(gctx, position);
double cosAlpha = normalT.dot(direction);
return std::abs(1. / cosAlpha);
}

const Acts::CylinderBounds& Acts::CylinderSurface::bounds() const {
const CylinderBounds& CylinderSurface::bounds() const {
return (*m_bounds.get());
}

Acts::Polyhedron Acts::CylinderSurface::polyhedronRepresentation(
Polyhedron CylinderSurface::polyhedronRepresentation(
const GeometryContext& gctx, unsigned int quarterSegments) const {
auto ctrans = transform(gctx);

@@ -199,13 +192,12 @@ Acts::Polyhedron Acts::CylinderSurface::polyhedronRepresentation(
return Polyhedron(vertices, faces, triangularMesh, false);
}

Acts::Vector3 Acts::CylinderSurface::rotSymmetryAxis(
const GeometryContext& gctx) const {
Vector3 CylinderSurface::rotSymmetryAxis(const GeometryContext& gctx) const {
// fast access via transform matrix (and not rotation())
return transform(gctx).matrix().block<3, 1>(0, 2);
}

Acts::detail::RealQuadraticEquation Acts::CylinderSurface::intersectionSolver(
detail::RealQuadraticEquation CylinderSurface::intersectionSolver(
const Transform3& transform, const Vector3& position,
const Vector3& direction) const {
// Solve for radius R
@@ -227,7 +219,7 @@ Acts::detail::RealQuadraticEquation Acts::CylinderSurface::intersectionSolver(
return detail::RealQuadraticEquation(a, b, c);
}

Acts::SurfaceMultiIntersection Acts::CylinderSurface::intersect(
SurfaceMultiIntersection CylinderSurface::intersect(
const GeometryContext& gctx, const Vector3& position,
const Vector3& direction, const BoundaryTolerance& boundaryTolerance,
double tolerance) const {
@@ -294,7 +286,7 @@ Acts::SurfaceMultiIntersection Acts::CylinderSurface::intersect(
return {{second, first}, this};
}

Acts::AlignmentToPathMatrix Acts::CylinderSurface::alignmentToPathDerivative(
AlignmentToPathMatrix CylinderSurface::alignmentToPathDerivative(
const GeometryContext& gctx, const Vector3& position,
const Vector3& direction) const {
assert(isOnSurface(gctx, position, direction, BoundaryTolerance::Infinite()));
@@ -342,8 +334,7 @@ Acts::AlignmentToPathMatrix Acts::CylinderSurface::alignmentToPathDerivative(
return alignToPath;
}

Acts::ActsMatrix<2, 3>
Acts::CylinderSurface::localCartesianToBoundLocalDerivative(
ActsMatrix<2, 3> CylinderSurface::localCartesianToBoundLocalDerivative(
const GeometryContext& gctx, const Vector3& position) const {
using VectorHelpers::perp;
using VectorHelpers::phi;
@@ -363,11 +354,10 @@ Acts::CylinderSurface::localCartesianToBoundLocalDerivative(
return loc3DToLocBound;
}

std::pair<std::shared_ptr<Acts::CylinderSurface>, bool>
Acts::CylinderSurface::mergedWith(const CylinderSurface& other,
BinningValue direction, bool externalRotation,
const Logger& logger) const {
using namespace Acts::UnitLiterals;
std::pair<std::shared_ptr<CylinderSurface>, bool> CylinderSurface::mergedWith(
const CylinderSurface& other, BinningValue direction, bool externalRotation,
const Logger& logger) const {
using namespace UnitLiterals;

ACTS_VERBOSE("Merging cylinder surfaces in " << binningValueName(direction)
<< " direction");
@@ -458,7 +448,7 @@ Acts::CylinderSurface::mergedWith(const CylinderSurface& other,
double otherHlPhi = other.bounds().get(CylinderBounds::eHalfPhiSector);
double otherAvgPhi = other.bounds().get(CylinderBounds::eAveragePhi);

if (direction == Acts::BinningValue::binZ) {
if (direction == BinningValue::binZ) {
// z shift must match the bounds

if (std::abs(otherLocal.linear().col(eY)[eX]) >= tolerance &&
@@ -505,7 +495,7 @@ Acts::CylinderSurface::mergedWith(const CylinderSurface& other,
return {Surface::makeShared<CylinderSurface>(newTransform, newBounds),
zShift < 0};

} else if (direction == Acts::BinningValue::binRPhi) {
} else if (direction == BinningValue::binRPhi) {
// no z shift is allowed
if (std::abs(translation[2]) > tolerance) {
ACTS_ERROR(
@@ -568,3 +558,5 @@ Acts::CylinderSurface::mergedWith(const CylinderSurface& other,
binningValueName(direction));
}
}

} // namespace Acts