Plugins/EDM4hep/ConvertTrackEDM4hepTest.cpp

0001 // This file is part of the Acts project.
0002 //
0003 // Copyright (C) 2023 CERN for the benefit of the Acts project
0004 //
0005 // This Source Code Form is subject to the terms of the Mozilla Public
0006 // License, v. 2.0. If a copy of the MPL was not distributed with this
0007 // file, You can obtain one at http://mozilla.org/MPL/2.0/.
0008
0009 #include <boost/test/data/test_case.hpp>
0010 #include <boost/test/tools/old/interface.hpp>
0011 #include <boost/test/unit_test.hpp>
0012 #include <boost/test/unit_test_suite.hpp>
0013
0014 #include "Acts/Definitions/TrackParametrization.hpp"
0015 #include "Acts/EventData/Charge.hpp"
0016 #include "Acts/EventData/GenericBoundTrackParameters.hpp"
0017 #include "Acts/EventData/MultiTrajectory.hpp"
0018 #include "Acts/EventData/TrackHelpers.hpp"
0019 #include "Acts/EventData/TrackStatePropMask.hpp"
0020 #include "Acts/EventData/VectorMultiTrajectory.hpp"
0021 #include "Acts/EventData/VectorTrackContainer.hpp"
0022 #include "Acts/Geometry/GeometryContext.hpp"
0023 #include "Acts/Plugins/EDM4hep/EDM4hepUtil.hpp"
0024 #include "Acts/Propagator/detail/CovarianceEngine.hpp"
0025 #include "Acts/Propagator/detail/JacobianEngine.hpp"
0026 #include "Acts/Surfaces/PerigeeSurface.hpp"
0027 #include "Acts/Surfaces/PlaneSurface.hpp"
0028 #include "Acts/Surfaces/Surface.hpp"
0029 #include "Acts/Tests/CommonHelpers/FloatComparisons.hpp"
0030 #include "Acts/Utilities/Logger.hpp"
0031 #include "Acts/Utilities/Zip.hpp"
0032
0033 #include <algorithm>
0034 #include <iostream>
0035 #include <random>
0036
0037 #include <edm4hep/TrackCollection.h>
0038
0039 using namespace Acts;
0040 using namespace Acts::UnitLiterals;
0041 BOOST_AUTO_TEST_SUITE(EDM4hepParameterConversion)
0042
0043 BOOST_AUTO_TEST_CASE(JacobianRoundtrip) {
0044   BoundVector par;
0045   par << 1_mm, 5_mm, 0.1, M_PI_2 * 0.9, -1 / 1_GeV, 5_ns;
0046
0047   BoundMatrix cov;
0048   cov.setIdentity();
0049
0050   double Bz = 2_T;
0051
0052   double tanLambda = std::tan(M_PI_2 - par[Acts::eBoundTheta]);
0053   double omega =
0054       par[Acts::eBoundQOverP] / std::sin(par[Acts::eBoundTheta]) * Bz;
0055
0056   auto J1 = EDM4hepUtil::detail::jacobianToEdm4hep(par[eBoundTheta],
0057                                                    par[eBoundQOverP], Bz);
0058
0059   BoundMatrix cov2 = J1 * cov * J1.transpose();
0060
0061   auto J2 = EDM4hepUtil::detail::jacobianFromEdm4hep(tanLambda, omega, Bz);
0062
0063   BoundMatrix cov3 = J2 * cov2 * J2.transpose();
0064
0065   CHECK_CLOSE_ABS(cov, cov3, 1e-9);
0066 }
0067
0068 BOOST_AUTO_TEST_CASE(ConvertTrackParametersToEdm4hepWithPerigee) {
0069   auto refSurface = Surface::makeShared<PerigeeSurface>(Vector3{50, 30, 20});
0070
0071   BoundVector par;
0072   par << 1_mm, 5_mm, 0, M_PI_2, -1 / 1_GeV,
0073       5_ns;  // -> perpendicular to perigee and pointing right, should be PCA
0074
0075   BoundMatrix cov;
0076   cov.setIdentity();
0077   cov(5, 5) = 25_ns;
0078
0079   BoundTrackParameters boundPar{refSurface, par, cov,
0080                                 ParticleHypothesis::pion()};
0081
0082   double Bz = 2_T;
0083
0084   Acts::GeometryContext gctx;
0085
0086   EDM4hepUtil::detail::Parameters converted =
0087       EDM4hepUtil::detail::convertTrackParametersToEdm4hep(gctx, Bz, boundPar);
0088
0089   BOOST_CHECK(converted.covariance.has_value());
0090   BOOST_CHECK(converted.surface);
0091
0092   // input is already on perigee, should not be modified
0093   BOOST_CHECK_EQUAL(par.template head<2>(),
0094                     converted.values.template head<2>());
0095   BOOST_CHECK_EQUAL(
0096       (converted.covariance.value().template topLeftCorner<4, 4>()),
0097       ActsSquareMatrix<4>::Identity());
0098   BOOST_CHECK_GT(converted.covariance.value()(4, 4), 0);
0099   BOOST_CHECK_EQUAL(converted.covariance.value()(5, 5), 25_ns);
0100
0101   // convert back for roundtrip test
0102
0103   BoundTrackParameters roundtripPar =
0104       EDM4hepUtil::detail::convertTrackParametersFromEdm4hep(Bz, converted);
0105
0106   BOOST_CHECK(roundtripPar.parameters().isApprox(boundPar.parameters()));
0107   BOOST_CHECK(roundtripPar.covariance().value().isApprox(
0108       boundPar.covariance().value()));
0109 }
0110
0111 BOOST_AUTO_TEST_CASE(ConvertTrackParametersToEdm4hepWithOutPerigee) {
0112   auto planeSurface = Surface::makeShared<PlaneSurface>(
0113       Vector3{50, 30, 20}, Vector3{1, 1, 0.3}.normalized());
0114
0115   BoundVector par;
0116   par << 1_mm, 5_mm, M_PI / 4., M_PI_2 * 0.9, -1 / 1_GeV, 5_ns;
0117
0118   BoundMatrix cov;
0119   cov.setIdentity();
0120   cov(5, 5) = 25_ns;
0121
0122   BoundTrackParameters planePar{planeSurface, par, cov,
0123                                 ParticleHypothesis::pion()};
0124
0125   double Bz = 2_T;
0126
0127   Acts::GeometryContext gctx;
0128
0129   EDM4hepUtil::detail::Parameters converted =
0130       EDM4hepUtil::detail::convertTrackParametersToEdm4hep(gctx, Bz, planePar);
0131
0132   BOOST_CHECK(converted.covariance.has_value());
0133   BOOST_CHECK(converted.surface);
0134
0135   // input is not a perigee, so new params should be at 0, 0 on ad-hoc perigee
0136   BOOST_CHECK_EQUAL(converted.values.template head<2>(), (Vector2{0, 0}));
0137   CHECK_CLOSE_ABS(converted.values[2], par[2], 1e-6);
0138
0139   BOOST_CHECK_EQUAL(converted.covariance.value()(0, 0), 1);
0140
0141   BOOST_CHECK_LT(converted.covariance.value()(1, 1), 1.2);
0142   BOOST_CHECK_GT(converted.covariance.value()(1, 1), 1);
0143
0144   CHECK_CLOSE_ABS(converted.covariance.value()(2, 2), 1, 1e-6);
0145
0146   BOOST_CHECK_GT(converted.covariance.value()(3, 3), 1);
0147   BOOST_CHECK_LT(converted.covariance.value()(3, 3), 1.2);
0148
0149   BOOST_CHECK_GT(converted.covariance.value()(4, 4), 0);
0150   BOOST_CHECK_EQUAL(converted.covariance.value()(5, 5), 25_ns);
0151
0152   // convert back for roundtrip test
0153   BoundTrackParameters roundtripPar =
0154       EDM4hepUtil::detail::convertTrackParametersFromEdm4hep(Bz, converted);
0155
0156   BOOST_CHECK_NE(dynamic_cast<const Acts::PerigeeSurface*>(
0157                      &roundtripPar.referenceSurface()),
0158                  nullptr);
0159
0160   BOOST_CHECK((converted.covariance.value().topLeftCorner<3, 3>().isApprox(
0161       roundtripPar.covariance().value().topLeftCorner<3, 3>())));
0162   CHECK_CLOSE_ABS(roundtripPar.covariance().value()(3, 3), 1, 1e-6);
0163   CHECK_CLOSE_ABS(roundtripPar.covariance().value()(4, 4), 1, 1e-6);
0164   BOOST_CHECK_EQUAL(roundtripPar.covariance().value()(5, 5), 25_ns);
0165
0166   auto roundtripPlaneBoundParams =
0167       Acts::detail::boundToBoundConversion(gctx, roundtripPar, *planeSurface)
0168           .value();
0169
0170   BOOST_CHECK(roundtripPlaneBoundParams.parameters().isApprox(par));
0171
0172   CHECK_CLOSE_COVARIANCE(roundtripPlaneBoundParams.covariance().value(),
0173                          planePar.covariance().value(), 1e-3);
0174 }
0175
0176 BOOST_AUTO_TEST_CASE(ConvertTrackParametersToEdm4hepWithPerigeeNoCov) {
0177   auto refSurface = Surface::makeShared<PerigeeSurface>(Vector3{50, 30, 20});
0178
0179   BoundVector par;
0180   par << 1_mm, 5_mm, 0, M_PI_2, -1 / 1_GeV,
0181       5_ns;  // -> perpendicular to perigee and pointing right, should be PCA
0182
0183   BoundTrackParameters boundPar{refSurface, par, std::nullopt,
0184                                 ParticleHypothesis::pion()};
0185
0186   double Bz = 2_T;
0187
0188   Acts::GeometryContext gctx;
0189
0190   EDM4hepUtil::detail::Parameters converted =
0191       EDM4hepUtil::detail::convertTrackParametersToEdm4hep(gctx, Bz, boundPar);
0192
0193   BOOST_CHECK(!converted.covariance.has_value());
0194   BOOST_CHECK(converted.surface);
0195
0196   // input is already on perigee, should not be modified
0197   BOOST_CHECK_EQUAL(par.template head<2>(),
0198                     converted.values.template head<2>());
0199
0200   // convert back for roundtrip test
0201
0202   BoundTrackParameters roundtripPar =
0203       EDM4hepUtil::detail::convertTrackParametersFromEdm4hep(Bz, converted);
0204
0205   BOOST_CHECK(roundtripPar.parameters().isApprox(boundPar.parameters()));
0206   BOOST_CHECK(!roundtripPar.covariance().has_value());
0207 }
0208
0209 BOOST_AUTO_TEST_CASE(ConvertTrackParametersToEdm4hepWithOutPerigeeNoCov) {
0210   auto refSurface = Surface::makeShared<PlaneSurface>(
0211       Vector3{50, 30, 20}, Vector3{1, 1, 0.3}.normalized());
0212
0213   BoundVector par;
0214   par << 1_mm, 5_mm, M_PI / 4., M_PI_2, -1 / 1_GeV, 5_ns;
0215
0216   BoundTrackParameters boundPar{refSurface, par, std::nullopt,
0217                                 ParticleHypothesis::pion()};
0218
0219   double Bz = 2_T;
0220
0221   Acts::GeometryContext gctx;
0222
0223   EDM4hepUtil::detail::Parameters converted =
0224       EDM4hepUtil::detail::convertTrackParametersToEdm4hep(gctx, Bz, boundPar);
0225
0226   BOOST_CHECK(!converted.covariance.has_value());
0227   BOOST_CHECK(converted.surface);
0228
0229   // input is not a perigee, so new params should be at 0, 0 on ad-hoc perigee
0230   BOOST_CHECK_EQUAL(converted.values.template head<2>(), (Vector2{0, 0}));
0231   CHECK_CLOSE_ABS(converted.values[2], par[2], 1e-6);
0232
0233   // convert back for roundtrip test
0234   BoundTrackParameters roundtripPar =
0235       EDM4hepUtil::detail::convertTrackParametersFromEdm4hep(Bz, converted);
0236
0237   BOOST_CHECK_EQUAL(roundtripPar.parameters().template head<2>(),
0238                     (Vector2{0, 0}));
0239   BOOST_CHECK(roundtripPar.parameters().tail<4>().isApprox(par.tail<4>()));
0240   BOOST_CHECK(!roundtripPar.covariance().has_value());
0241 }
0242
0243 BOOST_AUTO_TEST_CASE(CovariancePacking) {
0244   BoundMatrix m;
0245   // clang-format off
0246   m << 1, 2, 3, 4, 5, 6,
0247        2, 2, 3, 4, 5, 6,
0248        3, 3, 3, 4, 5, 6,
0249        4, 4, 4, 4, 5, 6,
0250        5, 5, 5, 5, 5, 6,
0251        6, 6, 6, 6, 6, 6;
0252   // clang-format on
0253
0254   std::array<float, 21> values{};
0255   EDM4hepUtil::detail::packCovariance(m, values.data());
0256
0257   BoundMatrix m2;
0258   m2.setZero();
0259   EDM4hepUtil::detail::unpackCovariance(values.data(), m2);
0260
0261   CHECK_CLOSE_ABS(m, m2, 1e-9);
0262 }
0263
0264 BOOST_AUTO_TEST_CASE(RoundTripTests) {
0265   auto trackContainer = std::make_shared<Acts::VectorTrackContainer>();
0266   auto trackStateContainer = std::make_shared<Acts::VectorMultiTrajectory>();
0267   TrackContainer tracks(trackContainer, trackStateContainer);
0268
0269   using mutable_proxy_t = decltype(tracks)::TrackProxy;
0270   using const_proxy_t = decltype(tracks)::ConstTrackProxy;
0271
0272   std::mt19937 rng{42};
0273   std::normal_distribution<double> gauss(0., 1.);
0274   std::uniform_real_distribution<double> f(-1, 1);
0275   std::uniform_real_distribution<double> r(0, 1);
0276   std::uniform_int_distribution<std::uint32_t> nTracks(2, 20);
0277   std::uniform_int_distribution<std::uint32_t> nTs(1, 20);
0278   std::uniform_real_distribution<double> phiDist(-M_PI, M_PI);
0279   std::uniform_real_distribution<double> etaDist(-4, 4);
0280   std::uniform_real_distribution<double> ptDist(1_MeV, 10_GeV);
0281   std::uniform_real_distribution<double> qDist(0., 1.);
0282
0283   auto genParams = [&]() -> std::pair<BoundVector, BoundMatrix> {
0284     double d0 = 20_um * gauss(rng);
0285     double z0 = 20_mm * gauss(rng);
0286     double phi = phiDist(rng);
0287     double eta = etaDist(rng);
0288     double theta = 2 * atan(exp(-eta));
0289     double pt = ptDist(rng);
0290     double p = pt / sin(theta);
0291     double charge = qDist(rng) > 0.5 ? 1. : -1.;
0292     double qop = charge / p;
0293     double t = 5_ns * gauss(rng);
0294
0295     BoundVector par;
0296     par << d0, z0, phi, theta, qop, t;
0297     BoundMatrix cov;
0298     cov = BoundMatrix::Identity();
0299     cov.diagonal() << 20_um * 20_um, 20_mm * 20_mm, 0.1, 0.1, 1_GeV, 25_ns;
0300     return {par, cov};
0301   };
0302
0303   std::uint32_t numT = nTracks(rng);
0304   for (std::uint32_t t = 0; t < numT; t++) {
0305     auto track = tracks.makeTrack();
0306     {
0307       auto [par, cov] = genParams();
0308       track.parameters() = par;
0309       track.covariance() = cov;
0310     }
0311     track.setReferenceSurface(
0312         Acts::Surface::makeShared<PerigeeSurface>(Vector3{0, 0, 0}));
0313
0314     std::uint32_t numTs = nTs(rng);
0315     for (std::uint32_t i = 0; i < numTs; i++) {
0316       auto ts = track.appendTrackState(TrackStatePropMask::Smoothed);
0317       double crit = r(rng);
0318       if (crit < 0.1) {
0319         ts.typeFlags().set(TrackStateFlag::HoleFlag);
0320         continue;
0321       } else if (crit < 0.2) {
0322         ts.typeFlags().set(TrackStateFlag::OutlierFlag);
0323         continue;
0324       } else if (crit < 0.3) {
0325         ts.typeFlags().set(TrackStateFlag::SharedHitFlag);
0326       } else if (crit < 0.4) {
0327         ts.typeFlags().set(TrackStateFlag::MaterialFlag);
0328         continue;
0329       }
0330
0331       ts.typeFlags().set(TrackStateFlag::MeasurementFlag);
0332
0333       auto [par, cov] = genParams();
0334       ts.smoothed() = par;
0335       ts.smoothedCovariance() = cov;
0336       Vector3 pos;
0337       pos << 1000 * f(rng), 1000 * f(rng), 3000 * f(rng);
0338       ts.setReferenceSurface(Acts::Surface::makeShared<PerigeeSurface>(pos));
0339     }
0340
0341     calculateTrackQuantities(track);
0342   }
0343
0344   edm4hep::TrackCollection edm4hepTracks;
0345
0346   Acts::GeometryContext gctx;
0347
0348   double Bz = 3_T;
0349
0350   auto logger = getDefaultLogger("EDM4hep", Logging::INFO);
0351
0352   for (const_proxy_t track : tracks) {
0353     auto to = edm4hepTracks.create();
0354     EDM4hepUtil::writeTrack(gctx, track, to, Bz, *logger);
0355   }
0356
0357   BOOST_CHECK_EQUAL(edm4hepTracks.size(), tracks.size());
0358
0359   auto tIt = tracks.begin();
0360   for (auto edm4hepTrack : edm4hepTracks) {
0361     auto track = *tIt;
0362     BOOST_CHECK_EQUAL(track.nMeasurements(),
0363                       edm4hepTrack.trackStates_size() - 1);
0364
0365     ++tIt;
0366   }
0367
0368   const edm4hep::TrackCollection& edm4hepTracksConst = edm4hepTracks;
0369
0370   TrackContainer readTracks(std::make_shared<Acts::VectorTrackContainer>(),
0371                             std::make_shared<Acts::VectorMultiTrajectory>());
0372
0373   for (const auto edm4hepTrack : edm4hepTracksConst) {
0374     EDM4hepUtil::readTrack(edm4hepTrack, readTracks.makeTrack(), Bz, *logger);
0375   }
0376
0377   BOOST_CHECK_EQUAL(tracks.size(), readTracks.size());
0378   std::size_t t = 0;
0379
0380   auto origTrackIt = tracks.begin();
0381   auto readTrackIt = readTracks.begin();
0382   while (origTrackIt != tracks.end() && readTrackIt != readTracks.end()) {
0383     BOOST_TEST_INFO_SCOPE("Track #" << t);
0384     auto orig = *origTrackIt;
0385     auto read = *readTrackIt;
0386
0387     CHECK_CLOSE_OR_SMALL(orig.parameters(), read.parameters(), 1e-5, 1e-8);
0388     CHECK_CLOSE_OR_SMALL(orig.covariance(), read.covariance(), 1e-5, 1e-8);
0389     BOOST_CHECK_EQUAL(orig.referenceSurface().center(gctx),
0390                       read.referenceSurface().center(gctx));
0391
0392     auto origTsIt = orig.trackStatesReversed().begin();
0393     auto readTsIt = read.trackStatesReversed().begin();
0394
0395     std::size_t tsi = 0;
0396     while (origTsIt != orig.trackStatesReversed().end() &&
0397            readTsIt != read.trackStatesReversed().end()) {
0398       BOOST_TEST_INFO_SCOPE("TS: #" << tsi);
0399       auto nextMeas = std::find_if(
0400           origTsIt, orig.trackStatesReversed().end(), [](const auto& ts) {
0401             return ts.typeFlags().test(TrackStateFlag::MeasurementFlag);
0402           });
0403       BOOST_CHECK(nextMeas != orig.trackStatesReversed().end());
0404       origTsIt = nextMeas;
0405       auto origTs = *origTsIt;
0406       auto readTs = *readTsIt;
0407
0408       BOOST_TEST_INFO_SCOPE(
0409           "orig parameters: " << origTs.parameters().transpose());
0410       BOOST_TEST_INFO_SCOPE(
0411           "read parameters: " << readTs.parameters().transpose());
0412       CHECK_CLOSE_OR_SMALL(origTs.smoothedCovariance(),
0413                            readTs.smoothedCovariance(), 1e-5, 1e-6);
0414       Vector3 newCenter = readTs.referenceSurface().center(
0415           gctx);  // new center is a perigee, but should be on the other
0416       // surface
0417       BOOST_CHECK(origTs.referenceSurface().isOnSurface(gctx, newCenter,
0418                                                         Vector3::Zero()));
0419
0420       // global hit positions should be the same
0421       Vector3 readGlobal = readTs.referenceSurface().localToGlobal(
0422           gctx, readTs.parameters().template head<2>(), Vector3::Zero());
0423       Vector3 origGlobal = origTs.referenceSurface().localToGlobal(
0424           gctx, origTs.parameters().template head<2>(), Vector3::Zero());
0425       CHECK_CLOSE_ABS(readGlobal, origGlobal, 1e-3);
0426       ++origTsIt;
0427       ++readTsIt;
0428       tsi++;
0429     }
0430     ++origTrackIt;
0431     ++readTrackIt;
0432
0433     t++;
0434   }
0435 }
0436
0437 BOOST_AUTO_TEST_SUITE_END()