Merge branch 'main' into addSimplePolygonBrepConv

.pre-commit-config.yaml

@@ -86,3 +86,10 @@ repos:
       language: system
       entry: docs/parse_cmake_options.py --write docs/getting_started.md
       files: ^CMakeLists.txt$
+
+  - repo: local
+    hooks:
+    - id: leftover_conflict_markers
+      name: Leftover conflict markers
+      language: system
+      entry: git diff --staged --check

Core/include/Acts/EventData/GenericFreeTrackParameters.hpp

@@ -13,6 +13,7 @@
 #include "Acts/Definitions/TrackParametrization.hpp"
 #include "Acts/EventData/TrackParametersConcept.hpp"
 #include "Acts/EventData/detail/PrintParameters.hpp"
+#include "Acts/Utilities/MathHelpers.hpp"
 #include "Acts/Utilities/UnitVectors.hpp"
 
 #include <cassert>
@@ -154,12 +155,12 @@ class GenericFreeTrackParameters {
     //   [f*sin(theta)*cos(phi), f*sin(theta)*sin(phi), f*cos(theta)]
     // w/ f,sin(theta) positive, the transverse magnitude is then
     //   sqrt(f^2*sin^2(theta)) = f*sin(theta)
-    Scalar transverseMagnitude =
-        std::hypot(m_params[eFreeDir0], m_params[eFreeDir1]);
+    Scalar transverseMagnitude2 =
+        square(m_params[eFreeDir0]) + square(m_params[eFreeDir1]);
     // absolute magnitude is f by construction
-    Scalar magnitude = std::hypot(transverseMagnitude, m_params[eFreeDir2]);
+    Scalar magnitude2 = transverseMagnitude2 + square(m_params[eFreeDir2]);
     // such that we can extract sin(theta) = f*sin(theta) / f
-    return (transverseMagnitude / magnitude) * absoluteMomentum();
+    return std::sqrt(transverseMagnitude2 / magnitude2) * absoluteMomentum();
   }
   /// Momentum three-vector.
   Vector3 momentum() const { return absoluteMomentum() * direction(); }

Core/include/Acts/Propagator/EigenStepperDenseExtension.hpp

@@ -17,6 +17,7 @@
 #include "Acts/Material/MaterialSlab.hpp"
 #include "Acts/Propagator/EigenStepperDefaultExtension.hpp"
 #include "Acts/Propagator/Propagator.hpp"
+#include "Acts/Utilities/MathHelpers.hpp"
 
 namespace Acts {
 
@@ -108,7 +109,7 @@ struct EigenStepperDenseExtension {
       // Evaluate k for the time propagation
       Lambdappi[0] = -qop[0] * qop[0] * qop[0] * g * energy[0] / (q * q);
       //~ tKi[0] = std::hypot(1, mass / initialMomentum);
-      tKi[0] = std::hypot(1, mass * qop[0]);
+      tKi[0] = fastHypot(1, mass * qop[0]);
       kQoP[0] = Lambdappi[0];
     } else {
       // Update parameters and check for momentum condition
@@ -122,7 +123,7 @@ struct EigenStepperDenseExtension {
       // Evaluate k_i for the time propagation
       auto qopNew = qop[0] + h * Lambdappi[i - 1];
       Lambdappi[i] = -qopNew * qopNew * qopNew * g * energy[i] / (q * q);
-      tKi[i] = std::hypot(1, mass * qopNew);
+      tKi[i] = fastHypot(1, mass * qopNew);
       kQoP[i] = Lambdappi[i];
     }
     return true;
@@ -167,14 +168,14 @@ struct EigenStepperDenseExtension {
     }
 
     // Add derivative dlambda/ds = Lambda''
-    state.stepping.derivative(7) = -std::hypot(mass, newMomentum) * g /
+    state.stepping.derivative(7) = -fastHypot(mass, newMomentum) * g /
                                    (newMomentum * newMomentum * newMomentum);
 
     // Update momentum
     state.stepping.pars[eFreeQOverP] =
         stepper.charge(state.stepping) / newMomentum;
     // Add derivative dt/ds = 1/(beta * c) = sqrt(m^2 * p^{-2} + c^{-2})
-    state.stepping.derivative(3) = std::hypot(1, mass / newMomentum);
+    state.stepping.derivative(3) = fastHypot(1, mass / newMomentum);
     // Update time
     state.stepping.pars[eFreeTime] +=
         (h / 6.) * (tKi[0] + 2. * (tKi[1] + tKi[2]) + tKi[3]);
@@ -332,25 +333,25 @@ struct EigenStepperDenseExtension {
     //~ (3. * g + qop[0] * dgdqop(energy[0], .mass,
     //~ absPdg, meanEnergyLoss));
 
-    double dtp1dl = qop[0] * mass * mass / std::hypot(1, qop[0] * mass);
+    double dtp1dl = qop[0] * mass * mass / fastHypot(1, qop[0] * mass);
     double qopNew = qop[0] + half_h * Lambdappi[0];
 
     //~ double dtpp2dl = -mass * mass * qopNew *
     //~ qopNew *
     //~ (3. * g * (1. + half_h * jdL[0]) +
     //~ qopNew * dgdqop(energy[1], mass, absPdgCode, meanEnergyLoss));
 
-    double dtp2dl = qopNew * mass * mass / std::hypot(1, qopNew * mass);
+    double dtp2dl = qopNew * mass * mass / fastHypot(1, qopNew * mass);
     qopNew = qop[0] + half_h * Lambdappi[1];
 
     //~ double dtpp3dl = -mass * mass * qopNew *
     //~ qopNew *
     //~ (3. * g * (1. + half_h * jdL[1]) +
     //~ qopNew * dgdqop(energy[2], mass, absPdg, meanEnergyLoss));
 
-    double dtp3dl = qopNew * mass * mass / std::hypot(1, qopNew * mass);
+    double dtp3dl = qopNew * mass * mass / fastHypot(1, qopNew * mass);
     qopNew = qop[0] + half_h * Lambdappi[2];
-    double dtp4dl = qopNew * mass * mass / std::hypot(1, qopNew * mass);
+    double dtp4dl = qopNew * mass * mass / fastHypot(1, qopNew * mass);
 
     //~ D(3, 7) = h * mass * mass * qop[0] /
     //~ std::hypot(1., mass * qop[0])
@@ -376,7 +377,7 @@ struct EigenStepperDenseExtension {
     PdgParticle absPdg = particleHypothesis.absolutePdg();
     float absQ = particleHypothesis.absoluteCharge();
 
-    energy[0] = std::hypot(initialMomentum, mass);
+    energy[0] = fastHypot(initialMomentum, mass);
     // use unit length as thickness to compute the energy loss per unit length
     MaterialSlab slab(material, 1);
     // Use the same energy loss throughout the step.
@@ -430,7 +431,7 @@ struct EigenStepperDenseExtension {
                         const stepper_t& stepper, const int i) {
     // Update parameters related to a changed momentum
     currentMomentum = initialMomentum + h * dPds[i - 1];
-    energy[i] = std::hypot(currentMomentum, mass);
+    energy[i] = fastHypot(currentMomentum, mass);
     dPds[i] = g * energy[i] / currentMomentum;
     qop[i] = stepper.charge(state.stepping) / currentMomentum;
     // Calculate term for later error propagation

Core/include/Acts/Propagator/StraightLineStepper.hpp

@@ -28,6 +28,7 @@
 #include "Acts/Surfaces/Surface.hpp"
 #include "Acts/Utilities/Intersection.hpp"
 #include "Acts/Utilities/Logger.hpp"
+#include "Acts/Utilities/MathHelpers.hpp"
 #include "Acts/Utilities/Result.hpp"
 
 #include <algorithm>
@@ -336,8 +337,8 @@ class StraightLineStepper {
           direction(prop_state.stepping);
       // dt / ds
       prop_state.stepping.derivative(eFreeTime) =
-          std::hypot(1., prop_state.stepping.particleHypothesis.mass() /
-                             absoluteMomentum(prop_state.stepping));
+          fastHypot(1., prop_state.stepping.particleHypothesis.mass() /
+                            absoluteMomentum(prop_state.stepping));
       // d (dr/ds) / ds : == 0
       prop_state.stepping.derivative.template segment<3>(4) =
           Acts::Vector3::Zero().transpose();
@@ -424,7 +425,7 @@ class StraightLineStepper {
     const auto m = state.stepping.particleHypothesis.mass();
     const auto p = absoluteMomentum(state.stepping);
     // time propagates along distance as 1/b = sqrt(1 + m²/p²)
-    const auto dtds = std::hypot(1., m / p);
+    const auto dtds = fastHypot(1., m / p);
     // Update the track parameters according to the equations of motion
     Vector3 dir = direction(state.stepping);
     state.stepping.pars.template segment<3>(eFreePos0) += h * dir;

Core/include/Acts/Propagator/detail/PointwiseMaterialInteraction.hpp

@@ -17,6 +17,7 @@
 #include "Acts/Material/ISurfaceMaterial.hpp"
 #include "Acts/Material/MaterialSlab.hpp"
 #include "Acts/Surfaces/Surface.hpp"
+#include "Acts/Utilities/MathHelpers.hpp"
 
 #include <algorithm>
 #include <cmath>
@@ -146,7 +147,7 @@ struct PointwiseMaterialInteraction {
     const auto& particleHypothesis = stepper.particleHypothesis(state.stepping);
     // in forward(backward) propagation, energy decreases(increases) and
     // variances increase(decrease)
-    const auto nextE = std::hypot(mass, momentum) - Eloss * navDir;
+    const auto nextE = fastHypot(mass, momentum) - Eloss * navDir;
     // put particle at rest if energy loss is too large
     nextP = (mass < nextE) ? std::sqrt(nextE * nextE - mass * mass) : 0;
     // minimum momentum below which we will not push particles via material

Core/include/Acts/Seeding/EstimateTrackParamsFromSeed.hpp

@@ -13,6 +13,7 @@
 #include "Acts/Geometry/GeometryContext.hpp"
 #include "Acts/Surfaces/Surface.hpp"
 #include "Acts/Utilities/Logger.hpp"
+#include "Acts/Utilities/MathHelpers.hpp"
 
 #include <array>
 #include <cmath>
@@ -240,12 +241,11 @@ std::optional<BoundVector> estimateTrackParamsFromSeed(
   int sign = ia > 0 ? -1 : 1;
   const ActsScalar R = circleCenter.norm();
   ActsScalar invTanTheta =
-      local2.z() /
-      (2.f * R * std::asin(std::hypot(local2.x(), local2.y()) / (2.f * R)));
+      local2.z() / (2.f * R * std::asin(local2.head<2>().norm() / (2.f * R)));
   // The momentum direction in the new frame (the center of the circle has the
   // coordinate (-1.*A/(2*B), 1./(2*B)))
   ActsScalar A = -circleCenter(0) / circleCenter(1);
-  Vector3 transDirection(1., A, std::hypot(1, A) * invTanTheta);
+  Vector3 transDirection(1., A, fastHypot(1, A) * invTanTheta);
   // Transform it back to the original frame
   Vector3 direction = rotation * transDirection.normalized();
 
@@ -275,7 +275,7 @@ std::optional<BoundVector> estimateTrackParamsFromSeed(
   // momentum on the transverse plane of the new frame)
   ActsScalar qOverPt = sign * (UnitConstants::m) / (0.3 * bFieldInTesla * R);
   // The estimated q/p in [GeV/c]^-1
-  params[eBoundQOverP] = qOverPt / std::hypot(1., invTanTheta);
+  params[eBoundQOverP] = qOverPt / fastHypot(1., invTanTheta);
 
   if (params.hasNaN()) {
     ACTS_ERROR(

Core/include/Acts/Surfaces/detail/BoundaryCheckHelper.hpp

@@ -10,6 +10,8 @@
 
 #include "Acts/Surfaces/detail/VerticesHelper.hpp"
 
+#include <span>
+
 namespace Acts::detail {
 
 /// Check if a point is inside a box.
@@ -72,8 +74,7 @@ inline bool insideAlignedBox(const Vector2& lowerLeft,
 /// @param jacobianOpt The Jacobian to transform the distance to Cartesian
 ///
 /// @return True if the point is inside the polygon.
-template <typename Vector2Container>
-inline bool insidePolygon(const Vector2Container& vertices,
+inline bool insidePolygon(std::span<const Vector2> vertices,
                           const BoundaryTolerance& tolerance,
                           const Vector2& point,
                           const std::optional<SquareMatrix2>& jacobianOpt) {

Core/include/Acts/Surfaces/detail/VerticesHelper.hpp

@@ -13,6 +13,7 @@
 
 #include <algorithm>
 #include <cmath>
+#include <span>
 #include <utility>
 #include <vector>
 
@@ -174,9 +175,8 @@ bool onHyperPlane(const std::vector<Vector3>& vertices,
                   ActsScalar tolerance = s_onSurfaceTolerance);
 
 /// Calculate the closest point on the polygon.
-template <typename Vector2Container>
 inline Vector2 computeClosestPointOnPolygon(const Vector2& point,
-                                            const Vector2Container& vertices,
+                                            std::span<const Vector2> vertices,
                                             const SquareMatrix2& metric) {
   auto squaredNorm = [&](const Vector2& x) {
     return (x.transpose() * metric * x).value();

Core/include/Acts/Utilities/MathHelpers.hpp

@@ -0,0 +1,30 @@
+// This file is part of the ACTS project.
+//
+// Copyright (C) 2016 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 https://mozilla.org/MPL/2.0/.
+
+#pragma once
+
+#include <cmath>
+
+namespace Acts {
+
+template <typename T>
+constexpr auto square(T x) {
+  return x * x;
+}
+
+template <typename... T>
+constexpr auto hypotSquare(T... args) {
+  return (square(args) + ...);
+}
+
+template <typename... T>
+constexpr auto fastHypot(T... args) {
+  return std::sqrt(hypotSquare(args...));
+}
+
+}  // namespace Acts

Core/include/Acts/Utilities/VectorHelpers.hpp

@@ -73,7 +73,7 @@ double perp(const Eigen::MatrixBase<Derived>& v) noexcept {
     assert(v.rows() >= 2 &&
            "Perp function not valid for vectors not at least 2D");
   }
-  return std::hypot(v[0], v[1]);
+  return v.template head<2>().norm();
 }
 
 /// Calculate the theta angle (longitudinal w.r.t. z axis) of a vector

Core/src/Material/Interactions.cpp

@@ -10,6 +10,7 @@
 
 #include "Acts/Definitions/PdgParticle.hpp"
 #include "Acts/Material/Material.hpp"
+#include "Acts/Utilities/MathHelpers.hpp"
 
 #include <cassert>
 #include <cmath>
@@ -51,7 +52,7 @@ struct RelativisticQuantities {
     betaGamma = pOverM;
     assert((betaGamma >= 0) && "Negative betaGamma");
     // gamma = sqrt(m² + p²)/m = sqrt(1 + (p/m)²)
-    gamma = std::sqrt(1.0f + pOverM * pOverM);
+    gamma = Acts::fastHypot(1.0f, pOverM);
   }
 };
 
@@ -395,7 +396,7 @@ float Acts::computeEnergyLossRadiative(const MaterialSlab& slab,
   // particle momentum and energy
   // do not need to care about the sign since it is only used squared
   const float momentum = absQ / qOverP;
-  const float energy = std::hypot(m, momentum);
+  const float energy = fastHypot(m, momentum);
 
   float dEdx = computeBremsstrahlungLossMean(m, energy);
 
@@ -424,7 +425,7 @@ float Acts::deriveEnergyLossRadiativeQOverP(const MaterialSlab& slab,
   // particle momentum and energy
   // do not need to care about the sign since it is only used squared
   const float momentum = absQ / qOverP;
-  const float energy = std::hypot(m, momentum);
+  const float energy = fastHypot(m, momentum);
 
   // compute derivative w/ respect to energy.
   float derE = deriveBremsstrahlungLossMeanE(m);

Core/src/Surfaces/ConvexPolygonBounds.cpp

@@ -8,6 +8,7 @@
 
 #include "Acts/Surfaces/ConvexPolygonBounds.hpp"
 
+#include "Acts/Definitions/Algebra.hpp"
 #include "Acts/Surfaces/BoundaryTolerance.hpp"
 #include "Acts/Surfaces/detail/BoundaryCheckHelper.hpp"
 
@@ -51,8 +52,9 @@ Acts::ConvexPolygonBounds<Acts::PolygonDynamic>::type() const {
 bool Acts::ConvexPolygonBounds<Acts::PolygonDynamic>::inside(
     const Acts::Vector2& lposition,
     const Acts::BoundaryTolerance& boundaryTolerance) const {
-  return detail::insidePolygon(m_vertices, boundaryTolerance, lposition,
-                               std::nullopt);
+  return detail::insidePolygon(
+      std::span<const Vector2>(m_vertices.data(), m_vertices.size()),
+      boundaryTolerance, lposition, std::nullopt);
 }
 
 std::vector<Acts::Vector2> Acts::ConvexPolygonBounds<

Core/src/Surfaces/DiamondBounds.cpp

@@ -8,6 +8,7 @@
 
 #include "Acts/Surfaces/DiamondBounds.hpp"
 
+#include "Acts/Definitions/Algebra.hpp"
 #include "Acts/Surfaces/BoundaryTolerance.hpp"
 #include "Acts/Surfaces/detail/BoundaryCheckHelper.hpp"
 
@@ -22,7 +23,17 @@ Acts::SurfaceBounds::BoundsType Acts::DiamondBounds::type() const {
 bool Acts::DiamondBounds::inside(
     const Acts::Vector2& lposition,
     const Acts::BoundaryTolerance& boundaryTolerance) const {
-  return detail::insidePolygon(vertices(), boundaryTolerance, lposition,
+  // Vertices starting at lower left (min rel. phi)
+  // counter-clockwise
+  double x1 = get(DiamondBounds::eHalfLengthXnegY);
+  double y1 = get(DiamondBounds::eHalfLengthYneg);
+  double x2 = get(DiamondBounds::eHalfLengthXzeroY);
+  double y2 = 0.;
+  double x3 = get(DiamondBounds::eHalfLengthXposY);
+  double y3 = get(DiamondBounds::eHalfLengthYpos);
+  Vector2 vertices[] = {{-x1, -y1}, {x1, -y1}, {x2, y2},
+                        {x3, y3},   {-x3, y3}, {-x2, y2}};
+  return detail::insidePolygon(vertices, boundaryTolerance, lposition,
                                std::nullopt);
 }
 

Core/src/Surfaces/DiscSurface.cpp

@@ -202,7 +202,7 @@ Acts::Vector2 Acts::DiscSurface::localPolarToCartesian(
 
 Acts::Vector2 Acts::DiscSurface::localCartesianToPolar(
     const Vector2& lcart) const {
-  return Vector2(std::hypot(lcart[eBoundLoc0], lcart[eBoundLoc1]),
+  return Vector2(lcart.norm(),
                  std::atan2(lcart[eBoundLoc1], lcart[eBoundLoc0]));
 }
 

Core/src/Surfaces/TrapezoidBounds.cpp

@@ -89,7 +89,7 @@ bool Acts::TrapezoidBounds::inside(
 
   // at this stage, the point can only be in the triangles
   // run slow-ish polygon check
-  std::vector<Acts::Vector2> vertices = {
+  Vector2 vertices[] = {
       {-hlXnY, -hlY}, {hlXnY, -hlY}, {hlXpY, hlY}, {-hlXpY, hlY}};
   return detail::insidePolygon(vertices, boundaryTolerance, extPosition,
                                std::nullopt);