diff --git a/Core/include/Acts/Geometry/detail/Layer.ipp b/Core/include/Acts/Geometry/detail/Layer.ipp
index ed3d3ed1d3b35c5de8168048d88edceeddd6299d..ebbfeb6300a1e0544cf8dedbeaef5ad99d3dd4ad 100644
--- a/Core/include/Acts/Geometry/detail/Layer.ipp
+++ b/Core/include/Acts/Geometry/detail/Layer.ipp
@@ -246,16 +246,17 @@ const SurfaceIntersection Layer::surfaceOnApproach(
sIntersection.intersection.pathLength *=
std::copysign(1., options.navDir);
return sIntersection;
- } else if (sIntersection.alternatives.size() > 0.) {
+ } else if (sIntersection.alternative.status >=
+ Intersection::Status::reachable) {
// Test the alternative
- cLimit = sIntersection.alternatives[0].pathLength;
+ cLimit = sIntersection.alternative.pathLength;
withinLimit = (cLimit > oLimit and
cLimit * cLimit <= pLimit * pLimit + s_onSurfaceTolerance);
- if (sIntersection.alternatives[0] and withinLimit) {
+ if (sIntersection.alternative and withinLimit) {
// Set the right sign for the path length
- sIntersection.alternatives[0].pathLength *=
+ sIntersection.alternative.pathLength *=
std::copysign(1., options.navDir);
- return SurfaceIntersection(sIntersection.alternatives[0],
+ return SurfaceIntersection(sIntersection.alternative,
sIntersection.object);
}
}
diff --git a/Core/include/Acts/Geometry/detail/TrackingVolume.ipp b/Core/include/Acts/Geometry/detail/TrackingVolume.ipp
index a15dd93143347823208f44373938735c327fa917..16a5042eb13fe45869ecbbe7995ccbf49eb54526 100644
--- a/Core/include/Acts/Geometry/detail/TrackingVolume.ipp
+++ b/Core/include/Acts/Geometry/detail/TrackingVolume.ipp
@@ -110,15 +110,15 @@ std::vector<BoundaryIntersection> TrackingVolume::compatibleBoundaries(
sIntersection.object);
}
// Check the alternative
- if (!sIntersection.alternatives.empty()) {
+ if (sIntersection.alternative) {
// Test the alternative
- cLimit = sIntersection.alternatives[0].pathLength;
+ cLimit = sIntersection.alternative.pathLength;
withinLimit = (cLimit > oLimit and
cLimit * cLimit <= pLimit * pLimit + s_onSurfaceTolerance);
- if (sIntersection.alternatives[0] and withinLimit) {
- sIntersection.alternatives[0].pathLength *=
+ if (sIntersection.alternative and withinLimit) {
+ sIntersection.alternative.pathLength *=
std::copysign(1., options.navDir);
- return BoundaryIntersection(sIntersection.alternatives[0], bSurface,
+ return BoundaryIntersection(sIntersection.alternative, bSurface,
sIntersection.object);
}
}
diff --git a/Core/include/Acts/Propagator/DirectNavigator.hpp b/Core/include/Acts/Propagator/DirectNavigator.hpp
index 93ef7986bb58a96a6f8a242d5c3dcbbefe59eaec..a71a1a41002f89dab8ba2a40782d6904507d0425 100644
--- a/Core/include/Acts/Propagator/DirectNavigator.hpp
+++ b/Core/include/Acts/Propagator/DirectNavigator.hpp
@@ -172,9 +172,8 @@ class DirectNavigator {
// Intersect the next surface and go
double navStep = nextIntersection.intersection.pathLength;
double overstepLimit = stepper.overstepLimit(state.stepping);
- if (navStep < overstepLimit and !nextIntersection.alternatives.empty() and
- nextIntersection.alternatives[0]) {
- navStep = nextIntersection.alternatives[0].pathLength;
+ if (navStep < overstepLimit and nextIntersection.alternative) {
+ navStep = nextIntersection.alternative.pathLength;
}
// Set the step size - this has to be outsourced to the Stepper
diff --git a/Core/include/Acts/Propagator/detail/StandardAborters.hpp b/Core/include/Acts/Propagator/detail/StandardAborters.hpp
index efbf57a4caed4760fae385c0bcdc6ab5ebd9eac9..93ec9c9d85378b373abe8bc18ca441612d3ba460 100644
--- a/Core/include/Acts/Propagator/detail/StandardAborters.hpp
+++ b/Core/include/Acts/Propagator/detail/StandardAborters.hpp
@@ -187,11 +187,9 @@ struct SurfaceReached {
// Target is not reached, update the step size
const double overstepLimit = stepper.overstepLimit(state.stepping);
// Check the alternative solution
- if (distance < overstepLimit and
- sIntersection.alternatives.size() == 1 and
- sIntersection.alternatives[0]) {
+ if (distance < overstepLimit and sIntersection.alternative) {
// Update the distance to the alternative solution
- distance = sIntersection.alternatives[0].pathLength;
+ distance = sIntersection.alternative.pathLength;
}
state.stepping.stepSize.update(state.stepping.navDir * distance,
ConstrainedStep::aborter);
diff --git a/Core/include/Acts/Surfaces/BoundaryCheck.hpp b/Core/include/Acts/Surfaces/BoundaryCheck.hpp
index 3b53203bb85acbce678eb6b4f5525b4b73b78293..7ac3a9e7f59145f219bf2875ae39a37b9b07bfee 100644
--- a/Core/include/Acts/Surfaces/BoundaryCheck.hpp
+++ b/Core/include/Acts/Surfaces/BoundaryCheck.hpp
@@ -114,13 +114,19 @@ class BoundaryCheck {
double distance(const Vector2D& point, const Vector2D& lowerLeft,
const Vector2D& upperRight) const;
- private:
enum class Type {
eNone, ///< disable boundary check
eAbsolute, ///< absolute cut
eChi2 ///< chi2-based cut with full correlations
};
+ /// Broadcast the type
+ Type type() const;
+
+ // Broadcast the tolerance
+ const Vector2D& tolerance() const;
+
+ private:
/// Return a new BoundaryCheck with updated covariance.
/// @param jacobian Tranform Jacobian for the covariance
/// @warning This currently only transforms the covariance and does not work
@@ -171,6 +177,14 @@ class BoundaryCheck {
} // namespace Acts
+inline Acts::BoundaryCheck::Type Acts::BoundaryCheck::type() const {
+ return m_type;
+}
+
+inline const Acts::Vector2D& Acts::BoundaryCheck::tolerance() const {
+ return m_tolerance;
+}
+
inline Acts::BoundaryCheck::BoundaryCheck(bool check)
: m_weight(ActsSymMatrixD<2>::Identity()),
m_tolerance(0, 0),
diff --git a/Core/include/Acts/Surfaces/CylinderBounds.hpp b/Core/include/Acts/Surfaces/CylinderBounds.hpp
index 2c8570a1fda5caee3d165a913ddcd22c1f07b347..2f3e77d4169e448be560ee541b8337ddb584d6ab 100644
--- a/Core/include/Acts/Surfaces/CylinderBounds.hpp
+++ b/Core/include/Acts/Surfaces/CylinderBounds.hpp
@@ -120,6 +120,9 @@ class CylinderBounds : public SurfaceBounds {
/// This method returns the halflengthZ
double halflengthZ() const;
+ /// Returns true for full phi coverage
+ bool coversFullAzimuth() const;
+
private:
/// the bound radius, average, half phi and half Z
double m_radius, m_avgPhi, m_halfPhi, m_halfZ;
@@ -145,4 +148,9 @@ inline double CylinderBounds::halfPhiSector() const {
inline double CylinderBounds::halflengthZ() const {
return m_halfZ;
}
+
+inline bool CylinderBounds::coversFullAzimuth() const {
+ return (m_halfPhi == M_PI);
+}
+
} // namespace Acts
\ No newline at end of file
diff --git a/Core/include/Acts/Surfaces/CylinderSurface.hpp b/Core/include/Acts/Surfaces/CylinderSurface.hpp
index 36ededdc0ea31f6be79e059a3f5d46ba8efa0c68..3d9381a9341e967f5ddaec5ebca2b0d3b97e0136 100644
--- a/Core/include/Acts/Surfaces/CylinderSurface.hpp
+++ b/Core/include/Acts/Surfaces/CylinderSurface.hpp
@@ -281,7 +281,7 @@ class CylinderSurface : public Surface {
///
/// @return the quadratic equation
detail::RealQuadraticEquation intersectionSolver(
- const GeometryContext& gctx, const Vector3D& position,
+ const Transform3D& transform, const Vector3D& position,
const Vector3D& direction) const;
};
diff --git a/Core/include/Acts/Surfaces/DiscBounds.hpp b/Core/include/Acts/Surfaces/DiscBounds.hpp
index 496a2d6364c5fbc4e5ce2bcd87166c62f8e2dce9..2a0b0a1c654b005ee93d7b3da4f93a7471e6a5ba 100644
--- a/Core/include/Acts/Surfaces/DiscBounds.hpp
+++ b/Core/include/Acts/Surfaces/DiscBounds.hpp
@@ -20,6 +20,13 @@ namespace Acts {
/// common base class for all bounds that are in a r/phi frame
/// - simply introduced to avoid wrong bound assigments to surfaces
-class DiscBounds : public SurfaceBounds {};
+class DiscBounds : public SurfaceBounds {
+ public:
+ /// Returns true for full phi coverage
+ virtual bool coversFullAzimuth() const = 0;
+
+ /// Checks if it's inside the radius
+ virtual bool insideRadialBounds(double R, double tolerance = 0.) const = 0;
+};
} // namespace Acts
\ No newline at end of file
diff --git a/Core/include/Acts/Surfaces/DiscSurface.hpp b/Core/include/Acts/Surfaces/DiscSurface.hpp
index 7e34d6715939f1aae3cf720c64b1d4232d59e54f..8b6ff52144e32273ea4fc3e4b221c6f278bba205 100644
--- a/Core/include/Acts/Surfaces/DiscSurface.hpp
+++ b/Core/include/Acts/Surfaces/DiscSurface.hpp
@@ -20,6 +20,7 @@
#include "Acts/Surfaces/Surface.hpp"
#include "Acts/Surfaces/detail/PlanarHelper.hpp"
#include "Acts/Utilities/Definitions.hpp"
+#include "Acts/Utilities/Helpers.hpp"
namespace Acts {
diff --git a/Core/include/Acts/Surfaces/DiscTrapezoidalBounds.hpp b/Core/include/Acts/Surfaces/DiscTrapezoidalBounds.hpp
index 54e6244d040b79e642bb9b6c318abb8b46ddba98..ce2259138a90384f2a33806145d8d4936db13afd 100644
--- a/Core/include/Acts/Surfaces/DiscTrapezoidalBounds.hpp
+++ b/Core/include/Acts/Surfaces/DiscTrapezoidalBounds.hpp
@@ -108,6 +108,13 @@ class DiscTrapezoidalBounds : public DiscBounds {
/// This method returns the halflength in Y (this is Rmax -Rmin)
double halflengthY() const;
+ /// Returns true for full phi coverage - obviously false here
+ bool coversFullAzimuth() const final;
+
+ /// Checks if this is inside the radial coverage
+ /// given the a tolerance
+ bool insideRadialBounds(double R, double tolerance = 0.) const final;
+
private:
double m_rMin, m_rMax, m_minHalfX, m_maxHalfX, m_avgPhi;
double m_stereo; // TODO 2017-04-09 msmk: what is this good for?
@@ -167,4 +174,13 @@ inline double DiscTrapezoidalBounds::halflengthY() const {
return (hmax - hmin) / 2.0;
}
+inline bool DiscTrapezoidalBounds::coversFullAzimuth() const {
+ return false;
+}
+
+inline bool DiscTrapezoidalBounds::insideRadialBounds(double R,
+ double tolerance) const {
+ return (R + tolerance > m_rMin and R - tolerance < m_rMax);
+}
+
} // namespace Acts
\ No newline at end of file
diff --git a/Core/include/Acts/Surfaces/RadialBounds.hpp b/Core/include/Acts/Surfaces/RadialBounds.hpp
index 956a4eb82d7f02a9928b34220e9636553093333a..b458991b7d8915b97e2ecadb7345868618cfbf10 100644
--- a/Core/include/Acts/Surfaces/RadialBounds.hpp
+++ b/Core/include/Acts/Surfaces/RadialBounds.hpp
@@ -100,6 +100,13 @@ class RadialBounds : public DiscBounds {
/// Return method for the half phi span
double halfPhiSector() const;
+ /// Returns true for full phi coverage
+ bool coversFullAzimuth() const final;
+
+ /// Checks if this is inside the radial coverage
+ /// given the a tolerance
+ bool insideRadialBounds(double R, double tolerance = 0.) const final;
+
private:
double m_rMin, m_rMax, m_avgPhi, m_halfPhi;
@@ -126,4 +133,12 @@ inline double RadialBounds::halfPhiSector() const {
return m_halfPhi;
}
+inline bool RadialBounds::coversFullAzimuth() const {
+ return (m_halfPhi == M_PI);
+}
+
+inline bool RadialBounds::insideRadialBounds(double R, double tolerance) const {
+ return (R + tolerance > m_rMin and R - tolerance < m_rMax);
+}
+
} // namespace Acts
\ No newline at end of file
diff --git a/Core/include/Acts/Surfaces/detail/ConeSurface.ipp b/Core/include/Acts/Surfaces/detail/ConeSurface.ipp
index d6dd7122ee5c48f1d2a8223702f67f2f99d8bc69..0dc389f2a4522b41cfe89fe5bdb229928db8d2cf 100644
--- a/Core/include/Acts/Surfaces/detail/ConeSurface.ipp
+++ b/Core/include/Acts/Surfaces/detail/ConeSurface.ipp
@@ -108,12 +108,12 @@ inline SurfaceIntersection ConeSurface::intersect(
status2 == Intersection::Status::missed) {
// And add the alternative
if (qe.solutions > 1) {
- cIntersection.alternatives = {second};
+ cIntersection.alternative = second;
}
} else {
// And add the alternative
if (qe.solutions > 1) {
- cIntersection.alternatives = {first};
+ cIntersection.alternative = first;
}
cIntersection.intersection = second;
}
diff --git a/Core/include/Acts/Surfaces/detail/CylinderSurface.ipp b/Core/include/Acts/Surfaces/detail/CylinderSurface.ipp
index f5bf742e23b8b9597763c2a9168702e9fe898c85..c3f6a8a6bd343b95e64a1d5e1d29dce259addd02 100644
--- a/Core/include/Acts/Surfaces/detail/CylinderSurface.ipp
+++ b/Core/include/Acts/Surfaces/detail/CylinderSurface.ipp
@@ -17,13 +17,13 @@ inline const Vector3D CylinderSurface::rotSymmetryAxis(
}
inline detail::RealQuadraticEquation CylinderSurface::intersectionSolver(
- const GeometryContext& gctx, const Vector3D& position,
+ const Transform3D& transform, const Vector3D& position,
const Vector3D& direction) const {
// Solve for radius R
double R = bounds().r();
// Get the transformation matrtix
- const auto& tMatrix = transform(gctx).matrix();
+ const auto& tMatrix = transform.matrix();
Vector3D caxis = tMatrix.block<3, 1>(0, 2).transpose();
Vector3D ccenter = tMatrix.block<3, 1>(0, 3).transpose();
@@ -41,8 +41,9 @@ inline detail::RealQuadraticEquation CylinderSurface::intersectionSolver(
inline Intersection CylinderSurface::intersectionEstimate(
const GeometryContext& gctx, const Vector3D& position,
const Vector3D& direction, const BoundaryCheck& bcheck) const {
+ const auto& gctxTransform = transform(gctx);
// Solve the quadratic equation
- auto qe = intersectionSolver(gctx, position, direction);
+ auto qe = intersectionSolver(gctxTransform, position, direction);
// If no valid solution return a non-valid intersection
if (qe.solutions == 0) {
@@ -70,8 +71,10 @@ inline Intersection CylinderSurface::intersectionEstimate(
inline SurfaceIntersection CylinderSurface::intersect(
const GeometryContext& gctx, const Vector3D& position,
const Vector3D& direction, const BoundaryCheck& bcheck) const {
+ const auto& gctxTransform = transform(gctx);
+
// Solve the quadratic equation
- auto qe = intersectionSolver(gctx, position, direction);
+ auto qe = intersectionSolver(gctxTransform, position, direction);
// If no valid solution return a non-valid surfaceIntersection
if (qe.solutions == 0) {
@@ -84,23 +87,48 @@ inline SurfaceIntersection CylinderSurface::intersect(
qe.first * qe.first < s_onSurfaceTolerance * s_onSurfaceTolerance
? Intersection::Status::onSurface
: Intersection::Status::reachable;
- if (bcheck and not isOnSurface(gctx, solution1, direction, bcheck)) {
- status1 = Intersection::Status::missed;
- }
+ // Helper method for boundary check
+ auto boundaryCheck =
+ [&](const Vector3D& solution,
+ Intersection::Status status) -> Intersection::Status {
+ // No check to be done, return current status
+ if (!bcheck)
+ return status;
+ const auto& cBounds = bounds();
+ if (cBounds.coversFullAzimuth() and
+ bcheck.type() == BoundaryCheck::Type::eAbsolute) {
+ // Project out the current Z value via local z axis
+ // Built-in local to global for speed reasons
+ const auto& tMatrix = gctxTransform.matrix();
+ // Create the reference vector in local
+ const Vector3D vecLocal(solution - tMatrix.block<3, 1>(0, 3));
+ double cZ = vecLocal.dot(tMatrix.block<3, 1>(0, 2));
+ double tolerance = s_onSurfaceTolerance + bcheck.tolerance()[eLOC_Z];
+ double hZ = cBounds.halflengthZ() + tolerance;
+ return (cZ * cZ < hZ * hZ) ? status : Intersection::Status::missed;
+ }
+ return (isOnSurface(gctx, solution, direction, bcheck)
+ ? status
+ : Intersection::Status::missed);
+ };
+ // Check first solution for boundary compatiblity
+ status1 = boundaryCheck(solution1, status1);
+ // Set the intersection
+ Intersection first(solution1, qe.first, status1);
+ SurfaceIntersection cIntersection(first, this);
+ if (qe.solutions == 1) {
+ return cIntersection;
+ }
// Check the validity of the second solution
Vector3D solution2 = position + qe.second * direction;
Intersection::Status status2 =
qe.second * qe.second < s_onSurfaceTolerance * s_onSurfaceTolerance
? Intersection::Status::onSurface
: Intersection::Status::reachable;
- if (bcheck and not isOnSurface(gctx, solution2, direction, bcheck)) {
- status2 = Intersection::Status::missed;
- }
- // Set the intersection
- Intersection first(solution1, qe.first, status1);
+ // Check first solution for boundary compatiblity
+ status2 = boundaryCheck(solution2, status2);
Intersection second(solution2, qe.second, status2);
- SurfaceIntersection cIntersection(first, this);
// Check one if its valid or neither is valid
bool check1 = status1 != Intersection::Status::missed or
(status1 == Intersection::Status::missed and
@@ -109,14 +137,10 @@ inline SurfaceIntersection CylinderSurface::intersect(
if ((check1 and qe.first * qe.first < qe.second * qe.second) or
status2 == Intersection::Status::missed) {
// And add the alternative
- if (qe.solutions > 1) {
- cIntersection.alternatives = {second};
- }
+ cIntersection.alternative = second;
} else {
// And add the alternative
- if (qe.solutions > 1) {
- cIntersection.alternatives = {first};
- }
+ cIntersection.alternative = first;
cIntersection.intersection = second;
}
return cIntersection;
diff --git a/Core/include/Acts/Surfaces/detail/DiscSurface.ipp b/Core/include/Acts/Surfaces/detail/DiscSurface.ipp
index 7ac3ee399221ac9a2cd2ced487b843ab08a6c682..ed36924334331e085fa910a9defa87dcc148f508 100644
--- a/Core/include/Acts/Surfaces/detail/DiscSurface.ipp
+++ b/Core/include/Acts/Surfaces/detail/DiscSurface.ipp
@@ -110,15 +110,29 @@ inline const RotationMatrix3D DiscSurface::initJacobianToLocal(
inline Intersection DiscSurface::intersectionEstimate(
const GeometryContext& gctx, const Vector3D& position,
const Vector3D& direction, const BoundaryCheck& bcheck) const {
+ // Get the contextual transform
+ auto gctxTransform = transform(gctx);
// Use the intersection helper for planar surfaces
auto intersection =
- PlanarHelper::intersectionEstimate(transform(gctx), position, direction);
- // Evaluate (if necessary in terms of boundaries)
- // @todo: speed up isOnSurface - we know that it is on surface
- // all we need is to check if it's inside bounds
+ PlanarHelper::intersectionEstimate(gctxTransform, position, direction);
+ // Evaluate boundary check if requested (and reachable)
if (intersection.status != Intersection::Status::unreachable and bcheck and
- not isOnSurface(gctx, intersection.position, direction, bcheck)) {
- intersection.status = Intersection::Status::missed;
+ m_bounds != nullptr) {
+ // Built-in local to global for speed reasons
+ const auto& tMatrix = gctxTransform.matrix();
+ const Vector3D vecLocal(intersection.position - tMatrix.block<3, 1>(0, 3));
+ const Vector2D lcartesian =
+ tMatrix.block<3, 2>(0, 0).transpose() * vecLocal;
+ if (bcheck.type() == BoundaryCheck::Type::eAbsolute and
+ m_bounds->coversFullAzimuth()) {
+ double tolerance = s_onSurfaceTolerance + bcheck.tolerance()[eLOC_R];
+ if (not m_bounds->insideRadialBounds(VectorHelpers::perp(lcartesian),
+ tolerance)) {
+ intersection.status = Intersection::Status::missed;
+ }
+ } else if (not insideBounds(localCartesianToPolar(lcartesian), bcheck)) {
+ intersection.status = Intersection::Status::missed;
+ }
}
return intersection;
}
diff --git a/Core/include/Acts/Surfaces/detail/LineSurface.ipp b/Core/include/Acts/Surfaces/detail/LineSurface.ipp
index 036179c1903b2760d27f45425003a596f522db54..b120eb18cbc370bf656f240192b21a77eb3ff99a 100644
--- a/Core/include/Acts/Surfaces/detail/LineSurface.ipp
+++ b/Core/include/Acts/Surfaces/detail/LineSurface.ipp
@@ -114,14 +114,19 @@ inline Intersection LineSurface::intersectionEstimate(
status = (u * u < s_onSurfaceTolerance * s_onSurfaceTolerance)
? Intersection::Status::onSurface
: Intersection::Status::reachable;
- // evaluate validaty in terms of bounds
Vector3D result = (ma + u * ea);
- // it just needs to be a insideBounds() check
- // @todo there should be a faster check possible
- if (bcheck and not isOnSurface(gctx, result, direction, bcheck)) {
- status = Intersection::Status::missed;
+ // Evaluate the boundary check if requested
+ // m_bounds == nullptr prevents unecessary calulations for PerigeeSurface
+ if (bcheck and m_bounds) {
+ // At closest approach: check inside R or and inside Z
+ const Vector3D vecLocal(result - mb);
+ double cZ = vecLocal.dot(eb);
+ double hZ = m_bounds->halflengthZ() + s_onSurfaceTolerance;
+ if ((cZ * cZ > hZ * hZ) or ((vecLocal - cZ * eb).norm() >
+ m_bounds->r() + s_onSurfaceTolerance)) {
+ status = Intersection::Status::missed;
+ }
}
- // return the result with validity
return Intersection(result, u, status);
}
// return a false intersection
@@ -203,4 +208,4 @@ inline const BoundRowVector LineSurface::derivativeFactors(
// build the combined normal & longitudinal components
return (norm * (s_vec - pc * (long_mat * d_vec.asDiagonal() -
jacobian.block<3, BoundParsDim>(4, 0))));
-}
\ No newline at end of file
+}
diff --git a/Core/include/Acts/Surfaces/detail/PlaneSurface.ipp b/Core/include/Acts/Surfaces/detail/PlaneSurface.ipp
index 16cbcbc98c8c84ddeef3261c69a04b01766c40c4..f7565fec88ccd8dbf741ce6ba396b33bfd98d8a5 100644
--- a/Core/include/Acts/Surfaces/detail/PlaneSurface.ipp
+++ b/Core/include/Acts/Surfaces/detail/PlaneSurface.ipp
@@ -25,22 +25,28 @@ inline const Vector3D PlaneSurface::binningPosition(
inline double PlaneSurface::pathCorrection(const GeometryContext& gctx,
const Vector3D& position,
const Vector3D& direction) const {
- /// We can ignore the global position here
+ // We can ignore the global position here
return 1. / std::abs(Surface::normal(gctx, position).dot(direction));
}
inline Intersection PlaneSurface::intersectionEstimate(
const GeometryContext& gctx, const Vector3D& position,
const Vector3D& direction, const BoundaryCheck& bcheck) const {
+ // Get the contextual transform
+ const auto& gctxTransform = transform(gctx);
// Use the intersection helper for planar surfaces
auto intersection =
- PlanarHelper::intersectionEstimate(transform(gctx), position, direction);
- // Evaluate (if necessary in terms of boundaries)
- // @todo: speed up isOnSurface - we know that it is on surface
- // all we need is to check if it's inside bounds
- if (intersection.status != Intersection::Status::unreachable and bcheck and
- not isOnSurface(gctx, intersection.position, direction, bcheck)) {
- intersection.status = Intersection::Status::missed;
+ PlanarHelper::intersectionEstimate(gctxTransform, position, direction);
+ // Evaluate boundary check if requested (and reachable)
+ if (intersection.status != Intersection::Status::unreachable and bcheck) {
+ // Built-in local to global for speed reasons
+ const auto& tMatrix = gctxTransform.matrix();
+ // Create the reference vector in local
+ const Vector3D vecLocal(intersection.position - tMatrix.block<3, 1>(0, 3));
+ if (not insideBounds(tMatrix.block<3, 2>(0, 0).transpose() * vecLocal,
+ bcheck)) {
+ intersection.status = Intersection::Status::missed;
+ }
}
return intersection;
}
diff --git a/Core/include/Acts/Utilities/Intersection.hpp b/Core/include/Acts/Utilities/Intersection.hpp
index b8b42f2eece9dd22546fdc257d6f022aaff4fde9..8c300e8ea10731287145598e940295c6c9c703a8 100644
--- a/Core/include/Acts/Utilities/Intersection.hpp
+++ b/Core/include/Acts/Utilities/Intersection.hpp
@@ -92,7 +92,7 @@ class ObjectIntersection {
const representation_t* representation{nullptr};
/// The alternative intersections
- std::vector<Intersection> alternatives = {};
+ Intersection alternative{};
/// Default constructor
ObjectIntersection() = default;
@@ -123,7 +123,7 @@ class ObjectIntersection {
/// @brief Smaller operator for ordering & sorting
///
- /// This operator will ignore the alternatives, but simply
+ /// This operator will ignore the alternative, but simply
/// order the representing intersection
///
/// @param oi is the source intersection for comparison
@@ -134,7 +134,7 @@ class ObjectIntersection {
/// @brief Greater operator for ordering & sorting
///
- /// This operator will ignore the alternatives, but simply
+ /// This operator will ignore the alternative, but simply
/// order the representing intersection
///
/// @param oi is the source intersection for comparison
diff --git a/Tests/Core/Surfaces/CMakeLists.txt b/Tests/Core/Surfaces/CMakeLists.txt
index d81e0ea276f2a6fc62a8bdcdd4082b41f5f6def4..90dd6847b6f658447eb9d1db36628e1078c2d5e7 100644
--- a/Tests/Core/Surfaces/CMakeLists.txt
+++ b/Tests/Core/Surfaces/CMakeLists.txt
@@ -47,6 +47,10 @@ add_executable (SurfaceIntersectionTests SurfaceIntersectionTests.cpp)
target_link_libraries (SurfaceIntersectionTests PRIVATE ActsCore ActsTestsCommonHelpers ${Boost_UNIT_TEST_FRAMEWORK_LIBRARY})
add_test (NAME SurfaceIntersectionUnitTest COMMAND SurfaceIntersectionTests)
+add_executable (SurfaceIntersectionBenchmark SurfaceIntersectionBenchmark.cpp)
+target_link_libraries (SurfaceIntersectionBenchmark PRIVATE ActsCore ActsTestsCommonHelpers ${Boost_UNIT_TEST_FRAMEWORK_LIBRARY})
+add_test (NAME SurfaceIntersectionBenchmarkTest COMMAND SurfaceIntersectionBenchmark)
+
add_executable (InfiniteBoundsTests InfiniteBoundsTests.cpp)
target_link_libraries (InfiniteBoundsTests PRIVATE ActsCore ActsTestsCommonHelpers ${Boost_UNIT_TEST_FRAMEWORK_LIBRARY})
add_test (NAME InfiniteBoundsUnitTest COMMAND InfiniteBoundsTests)
diff --git a/Tests/Core/Surfaces/CylinderSurfaceTests.cpp b/Tests/Core/Surfaces/CylinderSurfaceTests.cpp
index 4748e3e1790ead75e537a6c10b5bd05c62463661..78c5a8b08a13c3f08bb2072981582d7370fd0166 100644
--- a/Tests/Core/Surfaces/CylinderSurfaceTests.cpp
+++ b/Tests/Core/Surfaces/CylinderSurfaceTests.cpp
@@ -191,11 +191,10 @@ BOOST_AUTO_TEST_CASE(CylinderSurfaceProperties) {
testContext, offSurface, direction, false);
BOOST_CHECK(bool(surfaceIntersect));
// there is a second solution & and it should be valid
- BOOST_CHECK(surfaceIntersect.alternatives.size() == 1);
- BOOST_CHECK(bool(*surfaceIntersect.alternatives.begin()));
+ BOOST_CHECK(surfaceIntersect.alternative);
// And it's path should be further away then the primary solution
double pn = surfaceIntersect.intersection.pathLength;
- double pa = surfaceIntersect.alternatives[0].pathLength;
+ double pa = surfaceIntersect.alternative.pathLength;
BOOST_CHECK(pn * pn < pa * pa);
//
/// Test pathCorrection
diff --git a/Tests/Core/Surfaces/SurfaceIntersectionBenchmark.cpp b/Tests/Core/Surfaces/SurfaceIntersectionBenchmark.cpp
new file mode 100644
index 0000000000000000000000000000000000000000..7b58fa37e45020b9f640cb005d5589fda4f074a5
--- /dev/null
+++ b/Tests/Core/Surfaces/SurfaceIntersectionBenchmark.cpp
@@ -0,0 +1,126 @@
+// This file is part of the Acts project.
+//
+// Copyright (C) 2019 CERN for the benefit of the Acts project
+//
+// This Source Code Form is subject to the terms of the Mozilla Public
+// License, v. 2.0. If a copy of the MPL was not distributed with this
+// file, You can obtain one at http://mozilla.org/MPL/2.0/.
+
+#define BOOST_TEST_MODULE Surface Intersection Benchmark
+
+#define BOOST_TEST_DYN_LINK
+#include <boost/test/unit_test.hpp>
+// leave blank line
+
+#include <boost/test/data/test_case.hpp>
+// leave blank line
+
+#include <boost/test/output_test_stream.hpp>
+// leave blank line
+
+#include <cmath>
+#include "Acts/Geometry/GeometryContext.hpp"
+#include "Acts/Surfaces/CylinderBounds.hpp"
+#include "Acts/Surfaces/CylinderSurface.hpp"
+#include "Acts/Surfaces/DiscSurface.hpp"
+#include "Acts/Surfaces/PlaneSurface.hpp"
+#include "Acts/Surfaces/RadialBounds.hpp"
+#include "Acts/Surfaces/RectangleBounds.hpp"
+#include "Acts/Surfaces/StrawSurface.hpp"
+#include "Acts/Utilities/Units.hpp"
+
+namespace bdata = boost::unit_test::data;
+namespace tt = boost::test_tools;
+using namespace Acts::UnitLiterals;
+
+namespace Acts {
+namespace Test {
+
+// Some randomness & number crunching
+unsigned int ntests = 100;
+unsigned int nrepts = 1000;
+const bool boundaryCheck = false;
+const bool testPlane = true;
+const bool testDisc = true;
+const bool testCylinder = true;
+const bool testStraw = true;
+
+// Create a test context
+GeometryContext tgContext = GeometryContext();
+
+// Create a test plane in 10 m distance
+// Some random transform
+Transform3D at = Transform3D::Identity() * Translation3D(0_m, 0_m, 10_m) *
+ AngleAxis3D(0.15, Vector3D(1.2, 1.2, 0.12).normalized());
+
+// Define the Plane surface
+auto rb = std::make_shared<RectangleBounds>(1_m, 1_m);
+auto aPlane = Surface::makeShared<PlaneSurface>(
+ std::make_shared<Transform3D>(at), std::move(rb));
+
+// Define the Disc surface
+auto db = std::make_shared<RadialBounds>(0.2_m, 1.2_m);
+auto aDisc = Surface::makeShared<DiscSurface>(std::make_shared<Transform3D>(at),
+ std::move(db));
+
+// Define a Cylinder surface
+auto cb = std::make_shared<CylinderBounds>(10_m, 100_m);
+auto aCylinder = Surface::makeShared<CylinderSurface>(
+ std::make_shared<Transform3D>(at), std::move(cb));
+
+// Define a Straw surface
+auto aStraw = Surface::makeShared<StrawSurface>(
+ std::make_shared<Transform3D>(at), 50_cm, 2_m);
+
+// The orgin of our attempts for plane, disc and cylinder
+Vector3D origin(0., 0., 0.);
+
+// The origin for straw/line attempts
+Vector3D originStraw(0.3_m, -0.2_m, 11_m);
+
+template <typename surface_t>
+void intersectionTest(const surface_t& surface, double phi, double theta) {
+ // Shoot at it
+ double cosPhi = std::cos(phi);
+ double sinPhi = std::sin(phi);
+ double cosTheta = std::cos(theta);
+ double sinTheta = std::sin(theta);
+
+ Vector3D direction(cosPhi * sinTheta, sinPhi * sinTheta, cosTheta);
+
+ for (unsigned int ir = 0; ir < nrepts; ++ir) {
+ auto intersect =
+ surface.intersect(tgContext, origin, direction, boundaryCheck);
+
+ (void)intersect;
+ }
+}
+
+BOOST_DATA_TEST_CASE(
+ benchmark_surface_intersections,
+ bdata::random(
+ (bdata::seed = 21,
+ bdata::distribution = std::uniform_real_distribution<>(-M_PI, M_PI))) ^
+ bdata::random((bdata::seed = 22,
+ bdata::distribution =
+ std::uniform_real_distribution<>(-0.3, 0.3))) ^
+ bdata::xrange(ntests),
+ phi, theta, index) {
+ (void)index;
+
+ if (testPlane) {
+ intersectionTest<PlaneSurface>(*aPlane, phi, theta);
+ }
+ if (testDisc) {
+ intersectionTest<DiscSurface>(*aDisc, phi, theta);
+ }
+ if (testCylinder) {
+ intersectionTest<CylinderSurface>(*aCylinder, phi, theta);
+ }
+ if (testStraw) {
+ intersectionTest<StrawSurface>(*aStraw, phi, theta + M_PI);
+ }
+}
+
+} // namespace Test
+} // namespace Acts
diff --git a/Tests/Core/Surfaces/SurfaceIntersectionTests.cpp b/Tests/Core/Surfaces/SurfaceIntersectionTests.cpp
index 33b4444a3475625b2d4fde7de46fa771a5a737b9..9c4afb65ce4822d49969065de208fbf5affe9089 100644
--- a/Tests/Core/Surfaces/SurfaceIntersectionTests.cpp
+++ b/Tests/Core/Surfaces/SurfaceIntersectionTests.cpp
@@ -79,12 +79,12 @@ BOOST_AUTO_TEST_CASE(CylinderIntersectionTests) {
BOOST_CHECK(aIntersection.intersection.status ==
Intersection::Status::onSurface);
// There MUST be a second solution
- BOOST_CHECK(aIntersection.alternatives.size() == 1);
+ BOOST_CHECK(aIntersection.alternative);
// The other intersection MUST be reachable
- BOOST_CHECK(aIntersection.alternatives[0].status ==
+ BOOST_CHECK(aIntersection.alternative.status ==
Intersection::Status::reachable);
// The other intersection is at 2 meter distance
- CHECK_CLOSE_ABS(aIntersection.alternatives[0].pathLength, -2_m,
+ CHECK_CLOSE_ABS(aIntersection.alternative.pathLength, -2_m,
s_onSurfaceTolerance);
// Intersect from the the center
@@ -97,12 +97,12 @@ BOOST_AUTO_TEST_CASE(CylinderIntersectionTests) {
BOOST_CHECK(cIntersection.intersection.status ==
Intersection::Status::reachable);
// There MUST be a second solution
- BOOST_CHECK(cIntersection.alternatives.size() == 1);
+ BOOST_CHECK(cIntersection.alternative);
// The other intersection MUST be reachable
- BOOST_CHECK(cIntersection.alternatives[0].status ==
+ BOOST_CHECK(cIntersection.alternative.status ==
Intersection::Status::reachable);
// There MUST be one forward one backwards solution
- BOOST_CHECK(cIntersection.alternatives[0].pathLength *
+ BOOST_CHECK(cIntersection.alternative.pathLength *
cIntersection.intersection.pathLength <
0);
@@ -116,12 +116,12 @@ BOOST_AUTO_TEST_CASE(CylinderIntersectionTests) {
BOOST_CHECK(oIntersection.intersection.status ==
Intersection::Status::reachable);
// There MUST be a second solution
- BOOST_CHECK(oIntersection.alternatives.size() == 1);
+ BOOST_CHECK(oIntersection.alternative);
// The other intersection MUST be reachable
- BOOST_CHECK(oIntersection.alternatives[0].status ==
+ BOOST_CHECK(oIntersection.alternative.status ==
Intersection::Status::reachable);
// There MUST be one forward one backwards solution
- BOOST_CHECK(oIntersection.alternatives[0].pathLength *
+ BOOST_CHECK(oIntersection.alternative.pathLength *
oIntersection.intersection.pathLength >
0);
@@ -144,27 +144,27 @@ BOOST_AUTO_TEST_CASE(CylinderIntersectionTests) {
BOOST_CHECK(eIntersection.intersection.status ==
Intersection::Status::reachable);
// There MUST be a second solution
- BOOST_CHECK(eIntersection.alternatives.size() == 1);
+ BOOST_CHECK(eIntersection.alternative);
// The other intersection MUST be reachable
- BOOST_CHECK(eIntersection.alternatives[0].status ==
+ BOOST_CHECK(eIntersection.alternative.status ==
Intersection::Status::reachable);
// And be the negative one
- BOOST_CHECK(eIntersection.alternatives[0].pathLength < 0.);
+ BOOST_CHECK(eIntersection.alternative.pathLength < 0.);
// Now re-do with boundary check
eIntersection = aCylinder->intersect(tgContext, atEdge, transTZ, true);
- // This should be the positive one
+ // This should be the negative one
BOOST_CHECK(eIntersection.intersection.pathLength < 0.);
// The status of this one should be reachable
BOOST_CHECK(eIntersection.intersection.status ==
Intersection::Status::reachable);
// There MUST be a second solution
- BOOST_CHECK(eIntersection.alternatives.size() == 1);
- // The other intersection MUST be reachable
- BOOST_CHECK(eIntersection.alternatives[0].status ==
+ BOOST_CHECK(!eIntersection.alternative);
+ // The other intersection MUST NOT be reachable
+ BOOST_CHECK(eIntersection.alternative.status ==
Intersection::Status::missed);
- // And be the negative one
- BOOST_CHECK(eIntersection.alternatives[0].pathLength > 0.);
+ // And be the positive one
+ BOOST_CHECK(eIntersection.alternative.pathLength > 0.);
};
// In a nominal world
@@ -205,12 +205,12 @@ BOOST_AUTO_TEST_CASE(ConeIntersectionTest) {
Intersection::Status::onSurface);
// There MUST be a second solution
- BOOST_CHECK(aIntersection.alternatives.size() == 1);
+ BOOST_CHECK(aIntersection.alternative);
// The other intersection MUST be reachable
- BOOST_CHECK(aIntersection.alternatives[0].status ==
+ BOOST_CHECK(aIntersection.alternative.status ==
Intersection::Status::reachable);
// The other intersection is at 2 meter distance
- CHECK_CLOSE_ABS(aIntersection.alternatives[0].pathLength, -4.,
+ CHECK_CLOSE_ABS(aIntersection.alternative.pathLength, -4.,
s_onSurfaceTolerance);
// Intersection from outside without chance of hitting the cylinder
@@ -256,6 +256,7 @@ BOOST_AUTO_TEST_CASE(PlanarIntersectionTest) {
// Intersect forward
auto fIntersection = aPlane->intersect(tgContext, before, direction, true);
+
// The intersection MUST be valid
BOOST_CHECK(fIntersection);
// The intersection MUST be reachable
@@ -264,7 +265,7 @@ BOOST_AUTO_TEST_CASE(PlanarIntersectionTest) {
// The path length MUST be positive
BOOST_CHECK(fIntersection.intersection.pathLength > 0.);
// The intersection MUST be unique
- BOOST_CHECK(fIntersection.alternatives.size() == 0);
+ BOOST_CHECK(!fIntersection.alternative);
// On surface intersection
auto oIntersection = aPlane->intersect(tgContext, onit, direction, true);
@@ -277,7 +278,7 @@ BOOST_AUTO_TEST_CASE(PlanarIntersectionTest) {
BOOST_CHECK(std::abs(oIntersection.intersection.pathLength) <
s_onSurfaceTolerance);
// The intersection MUST be unique
- BOOST_CHECK(oIntersection.alternatives.size() == 0);
+ BOOST_CHECK(!oIntersection.alternative);
// Intersect backwards
auto bIntersection = aPlane->intersect(tgContext, after, direction, true);
@@ -289,7 +290,7 @@ BOOST_AUTO_TEST_CASE(PlanarIntersectionTest) {
// The path length MUST be negative
BOOST_CHECK(bIntersection.intersection.pathLength < 0.);
// The intersection MUST be unique
- BOOST_CHECK(bIntersection.alternatives.size() == 0);
+ BOOST_CHECK(!bIntersection.alternative);
// An out of bounds attempt: missed
auto mIntersection = aPlane->intersect(tgContext, outside, direction, true);
@@ -301,7 +302,7 @@ BOOST_AUTO_TEST_CASE(PlanarIntersectionTest) {
// The path length MUST be negative
BOOST_CHECK(mIntersection.intersection.pathLength > 0.);
// The intersection MUST be unique
- BOOST_CHECK(mIntersection.alternatives.size() == 0);
+ BOOST_CHECK(!mIntersection.alternative);
// An invalid attempt
auto iIntersection = aPlane->intersect(tgContext, before, parallel, true);
@@ -311,7 +312,7 @@ BOOST_AUTO_TEST_CASE(PlanarIntersectionTest) {
BOOST_CHECK(iIntersection.intersection.status ==
Intersection::Status::unreachable);
// The intersection MUST be unique
- BOOST_CHECK(iIntersection.alternatives.size() == 0);
+ BOOST_CHECK(!iIntersection.alternative);
};
// In a nominal world
@@ -358,7 +359,7 @@ BOOST_AUTO_TEST_CASE(LineIntersectionTest) {
// The path length MUST be positive
BOOST_CHECK(fIntersection.intersection.pathLength > 0.);
// The intersection MUST be unique
- BOOST_CHECK(fIntersection.alternatives.size() == 0);
+ BOOST_CHECK(!fIntersection.alternative);
// On surface intersection - on the straw with random direction
auto oIntersection = aLine->intersect(tgContext, onit1, direction, true);
@@ -371,7 +372,7 @@ BOOST_AUTO_TEST_CASE(LineIntersectionTest) {
BOOST_CHECK(std::abs(oIntersection.intersection.pathLength) <
s_onSurfaceTolerance);
// The intersection MUST be unique
- BOOST_CHECK(oIntersection.alternatives.size() == 0);
+ BOOST_CHECK(!oIntersection.alternative);
// On surface intersecion - on the surface with normal vector
oIntersection = aLine->intersect(tgContext, onitP, normalP, true);
@@ -384,7 +385,7 @@ BOOST_AUTO_TEST_CASE(LineIntersectionTest) {
BOOST_CHECK(std::abs(oIntersection.intersection.pathLength) <
s_onSurfaceTolerance);
// The intersection MUST be unique
- BOOST_CHECK(oIntersection.alternatives.size() == 0);
+ BOOST_CHECK(!oIntersection.alternative);
// Intersect backwards
auto bIntersection = aLine->intersect(tgContext, after, direction, true);
@@ -396,7 +397,7 @@ BOOST_AUTO_TEST_CASE(LineIntersectionTest) {
// The path length MUST be negative
BOOST_CHECK(bIntersection.intersection.pathLength < 0.);
// The intersection MUST be unique
- BOOST_CHECK(bIntersection.alternatives.size() == 0);
+ BOOST_CHECK(!bIntersection.alternative);
// An out of bounds attempt: missed
auto mIntersection = aLine->intersect(tgContext, outside, direction, true);
@@ -408,7 +409,7 @@ BOOST_AUTO_TEST_CASE(LineIntersectionTest) {
// The path length MUST be negative
BOOST_CHECK(mIntersection.intersection.pathLength < 0.);
// The intersection MUST be unique
- BOOST_CHECK(mIntersection.alternatives.size() == 0);
+ BOOST_CHECK(!mIntersection.alternative);
// An invalid attempt
auto iIntersection = aLine->intersect(tgContext, before, parallel, true);
@@ -418,7 +419,7 @@ BOOST_AUTO_TEST_CASE(LineIntersectionTest) {
BOOST_CHECK(iIntersection.intersection.status ==
Intersection::Status::unreachable);
// The intersection MUST be unique
- BOOST_CHECK(iIntersection.alternatives.size() == 0);
+ BOOST_CHECK(!iIntersection.alternative);
};
// In a nominal world