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