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   // The time variance grows slightly: the path length from the plane surface
0151   // to the ad-hoc perigee varies with the other (uncertain) parameters, and
0152   // dt/ds couples that variation into the time.
0153   BOOST_CHECK_GT(converted.covariance.value()(5, 5), 25_ns);
0154   CHECK_CLOSE_REL(converted.covariance.value()(5, 5), 25_ns, 1e-4);
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_GT(roundtripPar.covariance().value()(5, 5), 25_ns);
0169   CHECK_CLOSE_REL(roundtripPar.covariance().value()(5, 5), 25_ns, 1e-4);
0170
0171   auto roundtripPlaneBoundParams =
0172       detail::boundToBoundConversion(gctx, roundtripPar, *planeSurface).value();
0173
0174   BOOST_CHECK(roundtripPlaneBoundParams.parameters().isApprox(par));
0175
0176   CHECK_CLOSE_COVARIANCE(roundtripPlaneBoundParams.covariance().value(),
0177                          planePar.covariance().value(), 1e-3);
0178 }
0179
0180 BOOST_AUTO_TEST_CASE(ConvertTrackParametersToEdm4hepWithPerigeeNoCov) {
0181   auto refSurface = Surface::makeShared<PerigeeSurface>(Vector3{50, 30, 20});
0182
0183   BoundVector par;
0184   par << 1_mm, 5_mm, 0, std::numbers::pi / 2., -1 / 1_GeV,
0185       5_ns;  // -> perpendicular to perigee and pointing right, should be PCA
0186
0187   BoundTrackParameters boundPar{refSurface, par, std::nullopt,
0188                                 ParticleHypothesis::pion()};
0189
0190   double Bz = 2_T;
0191
0192   auto gctx = GeometryContext::dangerouslyDefaultConstruct();
0193
0194   EDM4hepUtil::detail::Parameters converted =
0195       EDM4hepUtil::detail::convertTrackParametersToEdm4hep(gctx, Bz, boundPar);
0196
0197   BOOST_CHECK(!converted.covariance.has_value());
0198   BOOST_CHECK(converted.surface);
0199
0200   // input is already on perigee, should not be modified
0201   BOOST_CHECK_EQUAL(par.template head<2>(),
0202                     converted.values.template head<2>());
0203
0204   // convert back for roundtrip test
0205
0206   BoundTrackParameters roundtripPar =
0207       EDM4hepUtil::detail::convertTrackParametersFromEdm4hep(Bz, converted);
0208
0209   BOOST_CHECK(roundtripPar.parameters().isApprox(boundPar.parameters()));
0210   BOOST_CHECK(!roundtripPar.covariance().has_value());
0211 }
0212
0213 BOOST_AUTO_TEST_CASE(ConvertTrackParametersToEdm4hepWithOutPerigeeNoCov) {
0214   std::shared_ptr<PlaneSurface> refSurface =
0215       CurvilinearSurface(Vector3{50, 30, 20}, Vector3{1, 1, 0.3}.normalized())
0216           .planeSurface();
0217
0218   BoundVector par;
0219   par << 1_mm, 5_mm, std::numbers::pi / 4., std::numbers::pi / 2., -1 / 1_GeV,
0220       5_ns;
0221
0222   BoundTrackParameters boundPar{refSurface, par, std::nullopt,
0223                                 ParticleHypothesis::pion()};
0224
0225   double Bz = 2_T;
0226
0227   auto gctx = GeometryContext::dangerouslyDefaultConstruct();
0228
0229   EDM4hepUtil::detail::Parameters converted =
0230       EDM4hepUtil::detail::convertTrackParametersToEdm4hep(gctx, Bz, boundPar);
0231
0232   BOOST_CHECK(!converted.covariance.has_value());
0233   BOOST_CHECK(converted.surface);
0234
0235   // input is not a perigee, so new params should be at 0, 0 on ad-hoc perigee
0236   BOOST_CHECK_EQUAL(converted.values.template head<2>(), (Vector2{0, 0}));
0237   CHECK_CLOSE_ABS(converted.values[2], par[2], 1e-6);
0238
0239   // convert back for roundtrip test
0240   BoundTrackParameters roundtripPar =
0241       EDM4hepUtil::detail::convertTrackParametersFromEdm4hep(Bz, converted);
0242
0243   BOOST_CHECK_EQUAL(roundtripPar.parameters().template head<2>(),
0244                     (Vector2{0, 0}));
0245   BOOST_CHECK(roundtripPar.parameters().tail<4>().isApprox(par.tail<4>()));
0246   BOOST_CHECK(!roundtripPar.covariance().has_value());
0247 }
0248
0249 BOOST_AUTO_TEST_CASE(CovariancePacking) {
0250   BoundMatrix m;
0251   // clang-format off
0252   m << 1, 2, 3, 4, 5, 6,
0253        2, 2, 3, 4, 5, 6,
0254        3, 3, 3, 4, 5, 6,
0255        4, 4, 4, 4, 5, 6,
0256        5, 5, 5, 5, 5, 6,
0257        6, 6, 6, 6, 6, 6;
0258   // clang-format on
0259
0260   std::array<float, 21> values{};
0261   EDM4hepUtil::detail::packCovariance(m, values.data());
0262
0263   BoundMatrix m2;
0264   m2.setZero();
0265   EDM4hepUtil::detail::unpackCovariance(values.data(), m2);
0266
0267   CHECK_CLOSE_ABS(m, m2, 1e-9);
0268 }
0269
0270 BOOST_AUTO_TEST_CASE(RoundTripTests) {
0271   auto trackContainer = std::make_shared<VectorTrackContainer>();
0272   auto trackStateContainer = std::make_shared<VectorMultiTrajectory>();
0273   TrackContainer tracks(trackContainer, trackStateContainer);
0274
0275   std::mt19937 rng{42};
0276   std::normal_distribution<double> gauss(0., 1.);
0277   std::uniform_real_distribution<double> f(-1, 1);
0278   std::uniform_real_distribution<double> r(0, 1);
0279   std::uniform_int_distribution<std::uint32_t> nTracks(2, 20);
0280   std::uniform_int_distribution<std::uint32_t> nTs(1, 20);
0281   std::uniform_real_distribution<double> phiDist(-std::numbers::pi,
0282                                                  std::numbers::pi);
0283   std::uniform_real_distribution<double> etaDist(-4, 4);
0284   std::uniform_real_distribution<double> ptDist(1_MeV, 10_GeV);
0285   std::uniform_real_distribution<double> qDist(0., 1.);
0286
0287   auto genParams = [&]() -> std::pair<BoundVector, BoundMatrix> {
0288     double d0 = 20_um * gauss(rng);
0289     double z0 = 20_mm * gauss(rng);
0290     double phi = phiDist(rng);
0291     double eta = etaDist(rng);
0292     double theta = 2 * std::atan(exp(-eta));
0293     double pt = ptDist(rng);
0294     double p = pt / std::sin(theta);
0295     double charge = qDist(rng) > 0.5 ? 1. : -1.;
0296     double qop = charge / p;
0297     double t = 5_ns * gauss(rng);
0298
0299     BoundVector par;
0300     par << d0, z0, phi, theta, qop, t;
0301     BoundMatrix cov;
0302     cov = BoundMatrix::Identity();
0303     cov.diagonal() << 20_um * 20_um, 20_mm * 20_mm, 0.1, 0.1, 1_GeV, 25_ns;
0304     return {par, cov};
0305   };
0306
0307   std::uint32_t numT = nTracks(rng);
0308   for (std::uint32_t t = 0; t < numT; t++) {
0309     auto track = tracks.makeTrack();
0310     {
0311       auto [par, cov] = genParams();
0312       track.parameters() = par;
0313       track.covariance() = cov;
0314     }
0315     track.setReferenceSurface(
0316         Surface::makeShared<PerigeeSurface>(Vector3{0, 0, 0}));
0317
0318     std::uint32_t numTs = nTs(rng);
0319     for (std::uint32_t i = 0; i < numTs; i++) {
0320       auto ts = track.appendTrackState(TrackStatePropMask::Smoothed);
0321       double crit = r(rng);
0322       if (crit < 0.1) {
0323         ts.typeFlags().setIsHole();
0324       } else if (crit < 0.2) {
0325         ts.typeFlags().setIsOutlier();
0326       } else if (crit < 0.3) {
0327         ts.typeFlags().setIsSharedHit();
0328       } else if (crit < 0.4) {
0329         ts.typeFlags().setIsMaterial();
0330       } else {
0331         ts.typeFlags().setIsMeasurement();
0332       }
0333
0334       auto [par, cov] = genParams();
0335       ts.smoothed() = par;
0336       ts.smoothedCovariance() = cov;
0337       Vector3 pos;
0338       pos << 1000 * f(rng), 1000 * f(rng), 3000 * f(rng);
0339       ts.setReferenceSurface(Surface::makeShared<PerigeeSurface>(pos));
0340     }
0341
0342     calculateTrackQuantities(track);
0343   }
0344
0345   edm4hep::TrackCollection edm4hepTracks;
0346
0347   auto gctx = GeometryContext::dangerouslyDefaultConstruct();
0348
0349   double Bz = 3_T;
0350
0351   auto logger = getDefaultLogger("EDM4hep", Logging::INFO);
0352
0353   for (const auto& track : tracks) {
0354     auto to = edm4hepTracks.create();
0355     EDM4hepUtil::writeTrack(gctx, track, to, Bz, *logger);
0356   }
0357
0358   BOOST_CHECK_EQUAL(edm4hepTracks.size(), tracks.size());
0359
0360   auto tIt = tracks.begin();
0361   for (auto edm4hepTrack : edm4hepTracks) {
0362     auto track = *tIt;
0363     BOOST_CHECK_EQUAL(track.nMeasurements(),
0364                       edm4hepTrack.trackStates_size() - 1);
0365
0366     ++tIt;
0367   }
0368
0369   const edm4hep::TrackCollection& edm4hepTracksConst = edm4hepTracks;
0370
0371   TrackContainer readTracks(std::make_shared<VectorTrackContainer>(),
0372                             std::make_shared<VectorMultiTrajectory>());
0373
0374   for (const auto edm4hepTrack : edm4hepTracksConst) {
0375     auto track = readTracks.makeTrack();
0376     EDM4hepUtil::readTrack(edm4hepTrack, track, Bz, *logger);
0377   }
0378
0379   BOOST_CHECK_EQUAL(tracks.size(), readTracks.size());
0380   std::size_t t = 0;
0381
0382   auto origTrackIt = tracks.begin();
0383   auto readTrackIt = readTracks.begin();
0384   while (origTrackIt != tracks.end() && readTrackIt != readTracks.end()) {
0385     BOOST_TEST_INFO_SCOPE("Track #" << t);
0386     auto orig = *origTrackIt;
0387     auto read = *readTrackIt;
0388
0389     CHECK_CLOSE_OR_SMALL(orig.parameters(), read.parameters(), 1e-5, 1e-8);
0390     CHECK_CLOSE_OR_SMALL(orig.covariance(), read.covariance(), 1e-5, 1e-8);
0391     BOOST_CHECK_EQUAL(orig.referenceSurface().center(gctx),
0392                       read.referenceSurface().center(gctx));
0393
0394     auto origTsIt = orig.trackStatesReversed().begin();
0395     auto readTsIt = read.trackStatesReversed().begin();
0396
0397     std::size_t tsi = 0;
0398     while (origTsIt != orig.trackStatesReversed().end() &&
0399            readTsIt != read.trackStatesReversed().end()) {
0400       BOOST_TEST_INFO_SCOPE("TS: #" << tsi);
0401       auto nextMeas = std::find_if(
0402           origTsIt, orig.trackStatesReversed().end(),
0403           [](const auto& ts) { return ts.typeFlags().isMeasurement(); });
0404       BOOST_CHECK(nextMeas != orig.trackStatesReversed().end());
0405       origTsIt = nextMeas;
0406       auto origTs = *origTsIt;
0407       auto readTs = *readTsIt;
0408
0409       BOOST_TEST_INFO_SCOPE(
0410           "orig parameters: " << origTs.parameters().transpose());
0411       BOOST_TEST_INFO_SCOPE(
0412           "read parameters: " << readTs.parameters().transpose());
0413       CHECK_CLOSE_OR_SMALL(origTs.smoothedCovariance(),
0414                            readTs.smoothedCovariance(), 1e-5, 1e-6);
0415       Vector3 newCenter = readTs.referenceSurface().center(
0416           gctx);  // new center is a perigee, but should be on the other
0417       // surface
0418       BOOST_CHECK(origTs.referenceSurface().isOnSurface(gctx, newCenter,
0419                                                         Vector3::Zero()));
0420
0421       // global hit positions should be the same
0422       Vector3 readGlobal = readTs.referenceSurface().localToGlobal(
0423           gctx, readTs.parameters().template head<2>(), Vector3::Zero());
0424       Vector3 origGlobal = origTs.referenceSurface().localToGlobal(
0425           gctx, origTs.parameters().template head<2>(), Vector3::Zero());
0426       CHECK_CLOSE_ABS(readGlobal, origGlobal, 1e-3);
0427       ++origTsIt;
0428       ++readTsIt;
0429       tsi++;
0430     }
0431     ++origTrackIt;
0432     ++readTrackIt;
0433
0434     t++;
0435   }
0436 }
0437
0438 BOOST_AUTO_TEST_SUITE_END()
0439
0440 }  // namespace ActsTests