Skip to content

Commit

Permalink
fix: EDM4hep track covariance conversion (acts-project#2823)
Browse files Browse the repository at this point in the history 
Blocked by:
- acts-project#2781 

Co-authored-by: Andreas Stefl <487211+andiwand@users.noreply.github.com>
  • Loading branch information
@paulgessinger @andiwand @LaraCalic
2 people authored and LaraCalic committed Feb 10, 2024
1 parent c360a47 commit b52ed05
Show file tree
Hide file tree
Showing 2 changed files with 76 additions and 29 deletions.
32 changes: 15 additions & 17 deletions Plugins/EDM4hep/src/EDM4hepUtil.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -9,9 +9,12 @@
#include "Acts/Plugins/EDM4hep/EDM4hepUtil.hpp"

#include "Acts/Definitions/Algebra.hpp"
#include "Acts/Definitions/TrackParametrization.hpp"
#include "Acts/Definitions/Units.hpp"
#include "Acts/EventData/MultiTrajectory.hpp"
#include "Acts/EventData/MultiTrajectoryHelpers.hpp"
#include "Acts/Propagator/detail/CovarianceEngine.hpp"
#include "Acts/Propagator/detail/JacobianEngine.hpp"

#include "edm4hep/TrackState.h"

Expand All @@ -37,7 +40,7 @@ ActsSquareMatrix<6> jacobianToEdm4hep(double theta, double qOverP, double Bz) {
// =
//
// [ 2 ]
// [- cot (theta) - 1 0 ]
// [- csc (theta) 0 ]
// [ ]
// [-B*q/p*cos(theta) B ]
// [------------------ ----------]
Expand All @@ -46,10 +49,9 @@ ActsSquareMatrix<6> jacobianToEdm4hep(double theta, double qOverP, double Bz) {

ActsSquareMatrix<6> J;
J.setIdentity();
double cotTheta = std::tan(M_PI_2 + theta);
J(3, 3) = -cotTheta * cotTheta - 1; // d(tanLambda) / dTheta
J(4, 4) = Bz / std::sin(theta); // dOmega / d(qop)
double sinTheta = std::sin(theta);
J(3, 3) = -1.0 / (sinTheta * sinTheta);
J(4, 4) = Bz / sinTheta; // dOmega / d(qop)
J(4, 3) = -Bz * qOverP * std::cos(theta) /
(sinTheta * sinTheta); // dOmega / dTheta
return J;
Expand Down Expand Up @@ -116,14 +118,8 @@ Parameters convertTrackParametersToEdm4hep(const Acts::GeometryContext& gctx,

std::shared_ptr<const Acts::Surface> refSurface =
params.referenceSurface().getSharedPtr();

Acts::BoundVector targetPars = params.parameters();
std::optional<Acts::FreeVector> freePars;

auto makeFreePars = [&]() {
return Acts::detail::transformBoundToFreeParameters(
params.referenceSurface(), gctx, params.parameters());
};
std::optional<Acts::BoundSquareMatrix> targetCov = params.covariance();

// If the reference surface is a perigee surface, we use that. Otherwise
// we create a new perigee surface at the global position of the track
Expand All @@ -134,20 +130,22 @@ Parameters convertTrackParametersToEdm4hep(const Acts::GeometryContext& gctx,
// We need to convert to the target parameters
// Keep the free parameters around we might need them for the covariance
// conversion
freePars = makeFreePars();
targetPars = Acts::detail::transformFreeToBoundParameters(freePars.value(),
*refSurface, gctx)
.value();

auto perigeeParams = Acts::detail::boundToBoundConversion(
gctx, params, *refSurface, Vector3{0, 0, Bz})
.value();
targetPars = perigeeParams.parameters();
targetCov = perigeeParams.covariance();
}

Parameters result;
result.surface = refSurface;

// Only run covariance conversion if we have a covariance input
if (params.covariance()) {
if (targetCov) {
Acts::ActsSquareMatrix<6> J = jacobianToEdm4hep(
targetPars[eBoundTheta], targetPars[eBoundQOverP], Bz);
result.covariance = J * params.covariance().value() * J.transpose();
result.covariance = J * targetCov.value() * J.transpose();
}

result.values[0] = targetPars[Acts::eBoundLoc0];
Expand Down
73 changes: 61 additions & 12 deletions Tests/UnitTests/Plugins/EDM4hep/ConvertTrackEDM4hepTest.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -21,6 +21,8 @@
#include "Acts/EventData/VectorTrackContainer.hpp"
#include "Acts/Geometry/GeometryContext.hpp"
#include "Acts/Plugins/EDM4hep/EDM4hepUtil.hpp"
#include "Acts/Propagator/detail/CovarianceEngine.hpp"
#include "Acts/Propagator/detail/JacobianEngine.hpp"
#include "Acts/Surfaces/PerigeeSurface.hpp"
#include "Acts/Surfaces/PlaneSurface.hpp"
#include "Acts/Surfaces/Surface.hpp"
Expand All @@ -29,6 +31,7 @@
#include "Acts/Utilities/Zip.hpp"

#include <algorithm>
#include <iostream>
#include <random>

#include <edm4hep/TrackCollection.h>
Expand All @@ -37,6 +40,31 @@ using namespace Acts;
using namespace Acts::UnitLiterals;
BOOST_AUTO_TEST_SUITE(EDM4hepParameterConversion)

BOOST_AUTO_TEST_CASE(JacobianRoundtrip) {
BoundVector par;
par << 1_mm, 5_mm, 0.1, M_PI_2 * 0.9, -1 / 1_GeV, 5_ns;

BoundMatrix cov;
cov.setIdentity();

double Bz = 2_T;

double tanLambda = std::tan(M_PI_2 - par[Acts::eBoundTheta]);
double omega =
par[Acts::eBoundQOverP] / std::sin(par[Acts::eBoundTheta]) * Bz;

auto J1 = EDM4hepUtil::detail::jacobianToEdm4hep(par[eBoundTheta],
par[eBoundQOverP], Bz);

BoundMatrix cov2 = J1 * cov * J1.transpose();

auto J2 = EDM4hepUtil::detail::jacobianFromEdm4hep(tanLambda, omega, Bz);

BoundMatrix cov3 = J2 * cov2 * J2.transpose();

CHECK_CLOSE_ABS(cov, cov3, 1e-9);
}

BOOST_AUTO_TEST_CASE(ConvertTrackParametersToEdm4hepWithPerigee) {
auto refSurface = Surface::makeShared<PerigeeSurface>(Vector3{50, 30, 20});

Expand Down Expand Up @@ -81,25 +109,25 @@ BOOST_AUTO_TEST_CASE(ConvertTrackParametersToEdm4hepWithPerigee) {
}

BOOST_AUTO_TEST_CASE(ConvertTrackParametersToEdm4hepWithOutPerigee) {
auto refSurface = Surface::makeShared<PlaneSurface>(
auto planeSurface = Surface::makeShared<PlaneSurface>(
Vector3{50, 30, 20}, Vector3{1, 1, 0.3}.normalized());

BoundVector par;
par << 1_mm, 5_mm, M_PI / 4., M_PI_2, -1 / 1_GeV, 5_ns;
par << 1_mm, 5_mm, M_PI / 4., M_PI_2 * 0.9, -1 / 1_GeV, 5_ns;

BoundMatrix cov;
cov.setIdentity();
cov(5, 5) = 25_ns;

BoundTrackParameters boundPar{refSurface, par, cov,
BoundTrackParameters planePar{planeSurface, par, cov,
ParticleHypothesis::pion()};

double Bz = 2_T;

Acts::GeometryContext gctx;

EDM4hepUtil::detail::Parameters converted =
EDM4hepUtil::detail::convertTrackParametersToEdm4hep(gctx, Bz, boundPar);
EDM4hepUtil::detail::convertTrackParametersToEdm4hep(gctx, Bz, planePar);

BOOST_CHECK(converted.covariance.has_value());
BOOST_CHECK(converted.surface);
Expand All @@ -108,20 +136,41 @@ BOOST_AUTO_TEST_CASE(ConvertTrackParametersToEdm4hepWithOutPerigee) {
BOOST_CHECK_EQUAL(converted.values.template head<2>(), (Vector2{0, 0}));
CHECK_CLOSE_ABS(converted.values[2], par[2], 1e-6);

BOOST_CHECK((converted.covariance.value().template topLeftCorner<4, 4>())
.isApprox(ActsSquareMatrix<4>::Identity()));
BOOST_CHECK_EQUAL(converted.covariance.value()(0, 0), 1);

BOOST_CHECK_LT(converted.covariance.value()(1, 1), 1.2);
BOOST_CHECK_GT(converted.covariance.value()(1, 1), 1);

CHECK_CLOSE_ABS(converted.covariance.value()(2, 2), 1, 1e-6);

BOOST_CHECK_GT(converted.covariance.value()(3, 3), 1);
BOOST_CHECK_LT(converted.covariance.value()(3, 3), 1.2);

BOOST_CHECK_GT(converted.covariance.value()(4, 4), 0);
BOOST_CHECK_EQUAL(converted.covariance.value()(5, 5), 25_ns);

// convert back for roundtrip test
BoundTrackParameters roundtripPar =
EDM4hepUtil::detail::convertTrackParametersFromEdm4hep(Bz, converted);

BOOST_CHECK_EQUAL(roundtripPar.parameters().template head<2>(),
(Vector2{0, 0}));
BOOST_CHECK(roundtripPar.parameters().tail<4>().isApprox(par.tail<4>()));
BOOST_CHECK(roundtripPar.covariance().value().isApprox(
boundPar.covariance().value()));
BOOST_CHECK_NE(dynamic_cast<const Acts::PerigeeSurface*>(
&roundtripPar.referenceSurface()),
nullptr);

BOOST_CHECK((converted.covariance.value().topLeftCorner<3, 3>().isApprox(
roundtripPar.covariance().value().topLeftCorner<3, 3>())));
CHECK_CLOSE_ABS(roundtripPar.covariance().value()(3, 3), 1, 1e-6);
CHECK_CLOSE_ABS(roundtripPar.covariance().value()(4, 4), 1, 1e-6);
BOOST_CHECK_EQUAL(roundtripPar.covariance().value()(5, 5), 25_ns);

auto roundtripPlaneBoundParams =
Acts::detail::boundToBoundConversion(gctx, roundtripPar, *planeSurface)
.value();

BOOST_CHECK(roundtripPlaneBoundParams.parameters().isApprox(par));

CHECK_CLOSE_COVARIANCE(roundtripPlaneBoundParams.covariance().value(),
planePar.covariance().value(), 1e-3);
}

BOOST_AUTO_TEST_CASE(ConvertTrackParametersToEdm4hepWithPerigeeNoCov) {
Expand Down Expand Up @@ -365,7 +414,7 @@ BOOST_AUTO_TEST_CASE(RoundTripTests) {
readTs.smoothedCovariance(), 1e-5, 1e-6);
Vector3 newCenter = readTs.referenceSurface().center(
gctx); // new center is a perigee, but should be on the other
// surface
// surface
BOOST_CHECK(origTs.referenceSurface().isOnSurface(gctx, newCenter,
Vector3::Zero()));

Expand Down

0 comments on commit b52ed05

Please sign in to comment.