Plugins/EDM4hep/ConvertTrackEDM4hepTest.cpp

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