Skip to content

Commit

Permalink
perf: sqrt over hypot in Core
Browse files Browse the repository at this point in the history
  • Loading branch information
andiwand committed Oct 5, 2024
1 parent 6cea6b2 commit f7da5c8
Show file tree
Hide file tree
Showing 13 changed files with 49 additions and 40 deletions.
8 changes: 4 additions & 4 deletions Core/include/Acts/EventData/GenericFreeTrackParameters.hpp
Original file line number Diff line number Diff line change
Expand Up @@ -154,12 +154,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 =
std::pow(m_params[eFreeDir0], 2) + std::pow(m_params[eFreeDir1], 2);
// absolute magnitude is f by construction
Scalar magnitude = std::hypot(transverseMagnitude, m_params[eFreeDir2]);
Scalar magnitude2 = transverseMagnitude2 + std::pow(m_params[eFreeDir2], 2);
// 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(); }
Expand Down
28 changes: 17 additions & 11 deletions Core/include/Acts/Propagator/EigenStepperDenseExtension.hpp
Original file line number Diff line number Diff line change
Expand Up @@ -108,7 +108,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] = std::sqrt(1 + std::pow(mass * qop[0], 2));
kQoP[0] = Lambdappi[0];
} else {
// Update parameters and check for momentum condition
Expand All @@ -122,7 +122,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] = std::sqrt(1 + std::pow(mass * qopNew, 2));
kQoP[i] = Lambdappi[i];
}
return true;
Expand Down Expand Up @@ -167,14 +167,16 @@ struct EigenStepperDenseExtension {
}

// Add derivative dlambda/ds = Lambda''
state.stepping.derivative(7) = -std::hypot(mass, newMomentum) * g /
(newMomentum * newMomentum * newMomentum);
state.stepping.derivative(7) =
-std::sqrt(std::pow(mass, 2) + std::pow(newMomentum, 2)) * 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) =
std::sqrt(1 + std::pow(mass / newMomentum, 2));
// Update time
state.stepping.pars[eFreeTime] +=
(h / 6.) * (tKi[0] + 2. * (tKi[1] + tKi[2]) + tKi[3]);
Expand Down Expand Up @@ -332,25 +334,29 @@ 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 / std::sqrt(1 + std::pow(qop[0] * mass, 2));
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 / std::sqrt(1 + std::pow(qopNew * mass, 2));
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 / std::sqrt(1 + std::pow(qopNew * mass, 2));
qopNew = qop[0] + half_h * Lambdappi[2];
double dtp4dl = qopNew * mass * mass / std::hypot(1, qopNew * mass);
double dtp4dl =
qopNew * mass * mass / std::sqrt(1 + std::pow(qopNew * mass, 2));

//~ D(3, 7) = h * mass * mass * qop[0] /
//~ std::hypot(1., mass * qop[0])
Expand All @@ -376,7 +382,7 @@ struct EigenStepperDenseExtension {
PdgParticle absPdg = particleHypothesis.absolutePdg();
float absQ = particleHypothesis.absoluteCharge();

energy[0] = std::hypot(initialMomentum, mass);
energy[0] = std::sqrt(std::pow(initialMomentum, 2) + std::pow(mass, 2));
// 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.
Expand Down Expand Up @@ -430,7 +436,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] = std::sqrt(std::pow(currentMomentum, 2) + std::pow(mass, 2));
dPds[i] = g * energy[i] / currentMomentum;
qop[i] = stepper.charge(state.stepping) / currentMomentum;
// Calculate term for later error propagation
Expand Down
9 changes: 5 additions & 4 deletions Core/include/Acts/Propagator/StraightLineStepper.hpp
Original file line number Diff line number Diff line change
Expand Up @@ -335,9 +335,10 @@ class StraightLineStepper {
prop_state.stepping.derivative.template head<3>() =
direction(prop_state.stepping);
// dt / ds
prop_state.stepping.derivative(eFreeTime) =
std::hypot(1., prop_state.stepping.particleHypothesis.mass() /
absoluteMomentum(prop_state.stepping));
prop_state.stepping.derivative(eFreeTime) = std::sqrt(
1. + std::pow(prop_state.stepping.particleHypothesis.mass() /
absoluteMomentum(prop_state.stepping),
2));
// d (dr/ds) / ds : == 0
prop_state.stepping.derivative.template segment<3>(4) =
Acts::Vector3::Zero().transpose();
Expand Down Expand Up @@ -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 = std::sqrt(1. + std::pow(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;
Expand Down
Original file line number Diff line number Diff line change
Expand Up @@ -146,7 +146,8 @@ 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 =
std::sqrt(std::pow(mass, 2) + std::pow(momentum, 2)) - 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
Expand Down
7 changes: 3 additions & 4 deletions Core/include/Acts/Seeding/EstimateTrackParamsFromSeed.hpp
Original file line number Diff line number Diff line change
Expand Up @@ -242,12 +242,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();

Expand Down Expand Up @@ -277,7 +276,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 / std::sqrt(1. + std::pow(invTanTheta, 2));

if (params.hasNaN()) {
ACTS_ERROR(
Expand Down
2 changes: 1 addition & 1 deletion Core/include/Acts/Utilities/VectorHelpers.hpp
Original file line number Diff line number Diff line change
Expand Up @@ -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
Expand Down
4 changes: 2 additions & 2 deletions Core/src/Material/Interactions.cpp
Original file line number Diff line number Diff line change
Expand Up @@ -395,7 +395,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 = std::sqrt(std::pow(m, 2) + std::pow(momentum, 2));

float dEdx = computeBremsstrahlungLossMean(m, energy);

Expand Down Expand Up @@ -424,7 +424,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 = std::sqrt(std::pow(m, 2) + std::pow(momentum, 2));

// compute derivative w/ respect to energy.
float derE = deriveBremsstrahlungLossMeanE(m);
Expand Down
2 changes: 1 addition & 1 deletion Core/src/Surfaces/DiscSurface.cpp
Original file line number Diff line number Diff line change
Expand Up @@ -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]));
}

Expand Down
9 changes: 5 additions & 4 deletions Core/src/Utilities/SpacePointUtility.cpp
Original file line number Diff line number Diff line change
Expand Up @@ -76,7 +76,7 @@ SpacePointUtility::globalCoords(
//
auto x = globalPos[ePos0];
auto y = globalPos[ePos1];
auto scale = 2 / std::hypot(x, y);
auto scale = 2 / std::sqrt(std::pow(x, 2) + std::pow(y, 2));
ActsMatrix<2, 3> jacXyzToRhoZ = ActsMatrix<2, 3>::Zero();
jacXyzToRhoZ(0, ePos0) = scale * x;
jacXyzToRhoZ(0, ePos1) = scale * y;
Expand Down Expand Up @@ -112,8 +112,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 = std::sqrt(std::pow(var1, 2) + std::pow(var2, 2));
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);
Expand All @@ -138,7 +139,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;
Expand Down
2 changes: 1 addition & 1 deletion Core/src/Vertexing/HelicalTrackLinearizer.cpp
Original file line number Diff line number Diff line change
Expand Up @@ -85,7 +85,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 / std::sqrt(std::pow(p, 2) + std::pow(m0, 2));
// Transverse speed (i.e., speed in the x-y plane)
ActsScalar betaT = beta * sinTheta;

Expand Down
4 changes: 2 additions & 2 deletions Core/src/Vertexing/ImpactPointEstimator.cpp
Original file line number Diff line number Diff line change
Expand Up @@ -207,7 +207,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 / std::sqrt(std::pow(p, 2) + std::pow(m0, 2));

pcaStraightTrack[3] = timeOnTrack + distanceToPca / beta;
}
Expand Down Expand Up @@ -282,7 +282,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 / std::sqrt(std::pow(p, 2) + std::pow(m0, 2));

pca[3] = tP - rho / (beta * sinTheta) * (phi - phiP);
}
Expand Down
Original file line number Diff line number Diff line change
Expand Up @@ -182,9 +182,10 @@ 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::hypot(genParticle.px(), genParticle.py(), genParticle.pz()) *
1_GeV);
particle.setAbsoluteMomentum(std::sqrt(std::pow(genParticle.px(), 2) +
std::pow(genParticle.py(), 2) +
std::pow(genParticle.pz(), 2)) *
1_GeV);

particles.push_back(std::move(particle));
}
Expand Down
4 changes: 2 additions & 2 deletions Examples/Algorithms/TrackFinding/src/HoughTransformSeeder.cpp
Original file line number Diff line number Diff line change
Expand Up @@ -478,7 +478,7 @@ void ActsExamples::HoughTransformSeeder::addSpacePoints(
ACTS_DEBUG("Inserting " << spContainer.size() << " space points from "
<< isp->key());
for (auto& sp : spContainer) {
double r = std::hypot(sp.x(), sp.y());
double r = std::sqrt(std::pow(sp.x(), 2) + std::pow(sp.y(), 2));
double z = sp.z();
float phi = std::atan2(sp.y(), sp.x());
ResultUnsigned hitlayer = m_cfg.layerIDFinder(r).value();
Expand Down Expand Up @@ -548,7 +548,7 @@ void ActsExamples::HoughTransformSeeder::addMeasurements(
Acts::Vector3 globalFakeMom(1, 1, 1);
Acts::Vector3 globalPos =
surface->localToGlobal(ctx.geoContext, localPos, globalFakeMom);
double r = std::hypot(globalPos[Acts::ePos0], globalPos[Acts::ePos1]);
double r = globalPos.head<2>().norm();
double phi = std::atan2(globalPos[Acts::ePos1], globalPos[Acts::ePos0]);
double z = globalPos[Acts::ePos2];
ResultUnsigned hitlayer = m_cfg.layerIDFinder(r);
Expand Down

0 comments on commit f7da5c8

Please sign in to comment.