perf: std::sqrt over std::hypot in Core

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] = hypot(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] = hypot(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) = -hypot(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) = hypot(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 / hypot(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 / hypot(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 / hypot(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 / hypot(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] = hypot(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] = hypot(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));
+          hypot(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 = hypot(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 = hypot(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

@@ -14,6 +14,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>
@@ -242,12 +243,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, std::sqrt(1 + std::pow(A, 2)) * invTanTheta);
   // Transform it back to the original frame
   Vector3 direction = rotation * transDirection.normalized();
 
@@ -277,7 +277,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 / hypot(1., invTanTheta);
 
   if (params.hasNaN()) {
     ACTS_ERROR(

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>
+inline auto square(T x) {
+  return std::pow(x, 2);
+}
+
+template <typename... T>
+inline auto hypot2(T... args) {
+  return (square(args) + ...);
+}
+
+template <typename... T>
+inline auto hypot(T... args) {
+  return std::sqrt(hypot2(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 = hypot(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 = hypot(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 = hypot(m, momentum);
 
   // compute derivative w/ respect to energy.
   float derE = deriveBremsstrahlungLossMeanE(m);

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/Utilities/SpacePointUtility.cpp

@@ -15,6 +15,7 @@
 #include "Acts/SpacePointFormation/SpacePointBuilderOptions.hpp"
 #include "Acts/Surfaces/Surface.hpp"
 #include "Acts/Utilities/Helpers.hpp"
+#include "Acts/Utilities/MathHelpers.hpp"
 
 #include <algorithm>
 #include <cmath>
@@ -76,7 +77,7 @@ SpacePointUtility::globalCoords(
   //
   auto x = globalPos[ePos0];
   auto y = globalPos[ePos1];
-  auto scale = 2 / std::hypot(x, y);
+  auto scale = 2 / hypot(x, y);
   ActsMatrix<2, 3> jacXyzToRhoZ = ActsMatrix<2, 3>::Zero();
   jacXyzToRhoZ(0, ePos0) = scale * x;
   jacXyzToRhoZ(0, ePos1) = scale * y;
@@ -112,8 +113,9 @@ Vector2 SpacePointUtility::calcRhoZVars(
   const auto var2 = paramCovAccessor(slinkBack).second(0, 0);
 
   // strip1 and strip2 are tilted at +/- theta/2
-  double sigma_x = std::hypot(var1, var2) / (2 * sin(theta * 0.5));
-  double sigma_y = std::hypot(var1, var2) / (2 * cos(theta * 0.5));
+  double sigma = hypot(var1, var2);
+  double sigma_x = sigma / (2 * sin(theta * 0.5));
+  double sigma_y = sigma / (2 * cos(theta * 0.5));
 
   // projection to the surface with strip1.
   double sig_x1 = sigma_x * cos(0.5 * theta) + sigma_y * sin(0.5 * theta);
@@ -138,7 +140,7 @@ Vector2 SpacePointUtility::rhoZCovariance(const GeometryContext& gctx,
 
   auto x = globalPos[ePos0];
   auto y = globalPos[ePos1];
-  auto scale = 2 / std::hypot(x, y);
+  auto scale = 2 / globalPos.head<2>().norm();
   ActsMatrix<2, 3> jacXyzToRhoZ = ActsMatrix<2, 3>::Zero();
   jacXyzToRhoZ(0, ePos0) = scale * x;
   jacXyzToRhoZ(0, ePos1) = scale * y;

Core/src/Vertexing/HelicalTrackLinearizer.cpp

@@ -10,6 +10,7 @@
 
 #include "Acts/Propagator/PropagatorOptions.hpp"
 #include "Acts/Surfaces/PerigeeSurface.hpp"
+#include "Acts/Utilities/MathHelpers.hpp"
 #include "Acts/Vertexing/LinearizerTrackParameters.hpp"
 
 Acts::Result<Acts::LinearizedTrack>
@@ -85,7 +86,7 @@ Acts::HelicalTrackLinearizer::linearizeTrack(
   ActsScalar p = params.particleHypothesis().extractMomentum(qOvP);
 
   // Speed in units of c
-  ActsScalar beta = p / std::hypot(p, m0);
+  ActsScalar beta = p / hypot(p, m0);
   // Transverse speed (i.e., speed in the x-y plane)
   ActsScalar betaT = beta * sinTheta;
 

Core/src/Vertexing/ImpactPointEstimator.cpp

@@ -13,6 +13,7 @@
 #include "Acts/Propagator/PropagatorOptions.hpp"
 #include "Acts/Surfaces/PerigeeSurface.hpp"
 #include "Acts/Surfaces/PlaneSurface.hpp"
+#include "Acts/Utilities/MathHelpers.hpp"
 #include "Acts/Vertexing/VertexingError.hpp"
 
 namespace Acts {
@@ -207,7 +208,7 @@ Result<std::pair<Vector4, Vector3>> getDistanceAndMomentumImpl(
       ActsScalar p = trkParams.particleHypothesis().extractMomentum(qOvP);
 
       // Speed in units of c
-      ActsScalar beta = p / std::hypot(p, m0);
+      ActsScalar beta = p / hypot(p, m0);
 
       pcaStraightTrack[3] = timeOnTrack + distanceToPca / beta;
     }
@@ -282,7 +283,7 @@ Result<std::pair<Vector4, Vector3>> getDistanceAndMomentumImpl(
     ActsScalar p = trkParams.particleHypothesis().extractMomentum(qOvP);
 
     // Speed in units of c
-    ActsScalar beta = p / std::hypot(p, m0);
+    ActsScalar beta = p / hypot(p, m0);
 
     pca[3] = tP - rho / (beta * sinTheta) * (phi - phiP);
   }

Examples/Algorithms/GeneratorsPythia8/ActsExamples/Generators/Pythia8ProcessGenerator.cpp

@@ -8,6 +8,7 @@
 
 #include "ActsExamples/Generators/Pythia8ProcessGenerator.hpp"
 
+#include "Acts/Utilities/MathHelpers.hpp"
 #include "ActsExamples/EventData/SimVertex.hpp"
 #include "ActsFatras/EventData/Barcode.hpp"
 #include "ActsFatras/EventData/Particle.hpp"
@@ -183,7 +184,7 @@ Pythia8Generator::operator()(RandomEngine& rng) {
     // normalization/ units are not import for the direction
     particle.setDirection(genParticle.px(), genParticle.py(), genParticle.pz());
     particle.setAbsoluteMomentum(
-        std::hypot(genParticle.px(), genParticle.py(), genParticle.pz()) *
+        Acts::hypot(genParticle.px(), genParticle.py(), genParticle.pz()) *
         1_GeV);
 
     particles.push_back(std::move(particle));