add helper

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>
@@ -155,9 +156,9 @@ class GenericFreeTrackParameters {
     // w/ f,sin(theta) positive, the transverse magnitude is then
     //   sqrt(f^2*sin^2(theta)) = f*sin(theta)
     Scalar transverseMagnitude2 =
-        std::pow(m_params[eFreeDir0], 2) + std::pow(m_params[eFreeDir1], 2);
+        square(m_params[eFreeDir0]) + square(m_params[eFreeDir1]);
     // absolute magnitude is f by construction
-    Scalar magnitude2 = transverseMagnitude2 + std::pow(m_params[eFreeDir2], 2);
+    Scalar magnitude2 = transverseMagnitude2 + square(m_params[eFreeDir2]);
     // such that we can extract sin(theta) = f*sin(theta) / f
     return std::sqrt(transverseMagnitude2 / magnitude2) * absoluteMomentum();
   }

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::sqrt(1 + std::pow(mass * qop[0], 2));
+      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::sqrt(1 + std::pow(mass * qopNew, 2));
+      tKi[i] = hypot(1, mass * qopNew);
       kQoP[i] = Lambdappi[i];
     }
     return true;
@@ -167,16 +168,14 @@ struct EigenStepperDenseExtension {
     }
 
     // Add derivative dlambda/ds = Lambda''
-    state.stepping.derivative(7) =
-        -std::sqrt(std::pow(mass, 2) + std::pow(newMomentum, 2)) * g /
-        (newMomentum * newMomentum * newMomentum);
+    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::sqrt(1 + std::pow(mass / newMomentum, 2));
+    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]);
@@ -334,29 +333,25 @@ struct EigenStepperDenseExtension {
     //~ (3. * g + qop[0] * dgdqop(energy[0], .mass,
     //~ absPdg, meanEnergyLoss));
 
-    double dtp1dl =
-        qop[0] * mass * mass / std::sqrt(1 + std::pow(qop[0] * mass, 2));
+    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::sqrt(1 + std::pow(qopNew * mass, 2));
+    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::sqrt(1 + std::pow(qopNew * mass, 2));
+    double dtp3dl = qopNew * mass * mass / hypot(1, qopNew * mass);
     qopNew = qop[0] + half_h * Lambdappi[2];
-    double dtp4dl =
-        qopNew * mass * mass / std::sqrt(1 + std::pow(qopNew * mass, 2));
+    double dtp4dl = qopNew * mass * mass / hypot(1, qopNew * mass);
 
     //~ D(3, 7) = h * mass * mass * qop[0] /
     //~ std::hypot(1., mass * qop[0])
@@ -382,7 +377,7 @@ struct EigenStepperDenseExtension {
     PdgParticle absPdg = particleHypothesis.absolutePdg();
     float absQ = particleHypothesis.absoluteCharge();
 
-    energy[0] = std::sqrt(std::pow(initialMomentum, 2) + std::pow(mass, 2));
+    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.
@@ -436,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::sqrt(std::pow(currentMomentum, 2) + std::pow(mass, 2));
+    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>
@@ -335,10 +336,9 @@ class StraightLineStepper {
       prop_state.stepping.derivative.template head<3>() =
           direction(prop_state.stepping);
       // dt / ds
-      prop_state.stepping.derivative(eFreeTime) = std::sqrt(
-          1. + std::pow(prop_state.stepping.particleHypothesis.mass() /
-                            absoluteMomentum(prop_state.stepping),
-                        2));
+      prop_state.stepping.derivative(eFreeTime) =
+          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();
@@ -425,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::sqrt(1. + std::pow(m / p, 2));
+    const auto dtds = hypot(1., m / p, 2);
     // 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,8 +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::sqrt(std::pow(mass, 2) + std::pow(momentum, 2)) - 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>
@@ -276,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::sqrt(1. + std::pow(invTanTheta, 2));
+  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/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::sqrt(std::pow(m, 2) + std::pow(momentum, 2));
+  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::sqrt(std::pow(m, 2) + std::pow(momentum, 2));
+  const float energy = hypot(m, momentum);
 
   // compute derivative w/ respect to energy.
   float derE = deriveBremsstrahlungLossMeanE(m);

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::sqrt(std::pow(x, 2) + std::pow(y, 2));
+  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,7 +113,7 @@ Vector2 SpacePointUtility::calcRhoZVars(
   const auto var2 = paramCovAccessor(slinkBack).second(0, 0);
 
   // strip1 and strip2 are tilted at +/- theta/2
-  double sigma = std::sqrt(std::pow(var1, 2) + std::pow(var2, 2));
+  double sigma = hypot(var1, var2);
   double sigma_x = sigma / (2 * sin(theta * 0.5));
   double sigma_y = sigma / (2 * cos(theta * 0.5));
 

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::sqrt(std::pow(p, 2) + std::pow(m0, 2));
+  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::sqrt(std::pow(p, 2) + std::pow(m0, 2));
+      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::sqrt(std::pow(p, 2) + std::pow(m0, 2));
+    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"
@@ -182,10 +183,9 @@ Pythia8Generator::operator()(RandomEngine& rng) {
     particle.setPosition4(pos4);
     // normalization/ units are not import for the direction
     particle.setDirection(genParticle.px(), genParticle.py(), genParticle.pz());
-    particle.setAbsoluteMomentum(std::sqrt(std::pow(genParticle.px(), 2) +
-                                           std::pow(genParticle.py(), 2) +
-                                           std::pow(genParticle.pz(), 2)) *
-                                 1_GeV);
+    particle.setAbsoluteMomentum(
+        Acts::hypot(genParticle.px(), genParticle.py(), genParticle.pz()) *
+        1_GeV);
 
     particles.push_back(std::move(particle));
   }

Examples/Algorithms/TrackFinding/src/HoughTransformSeeder.cpp

@@ -15,6 +15,7 @@
 #include "Acts/Geometry/TrackingGeometry.hpp"
 #include "Acts/Surfaces/Surface.hpp"
 #include "Acts/Utilities/Enumerate.hpp"
+#include "Acts/Utilities/MathHelpers.hpp"
 #include "ActsExamples/EventData/GeometryContainers.hpp"
 #include "ActsExamples/EventData/Index.hpp"
 #include "ActsExamples/EventData/IndexSourceLink.hpp"
@@ -478,7 +479,7 @@ void ActsExamples::HoughTransformSeeder::addSpacePoints(
     ACTS_DEBUG("Inserting " << spContainer.size() << " space points from "
                             << isp->key());
     for (auto& sp : spContainer) {
-      double r = std::sqrt(std::pow(sp.x(), 2) + std::pow(sp.y(), 2));
+      double r = Acts::hypot(sp.x(), sp.y());
       double z = sp.z();
       float phi = std::atan2(sp.y(), sp.x());
       ResultUnsigned hitlayer = m_cfg.layerIDFinder(r).value();