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/Plugins/EDM4hep/EDM4hepUtil.hpp"
0020 #include "Acts/Propagator/detail/CovarianceEngine.hpp"
0021 #include "Acts/Propagator/detail/JacobianEngine.hpp"
0022 #include "Acts/Surfaces/CurvilinearSurface.hpp"
0023 #include "Acts/Surfaces/PerigeeSurface.hpp"
0024 #include "Acts/Surfaces/PlaneSurface.hpp"
0025 #include "Acts/Surfaces/Surface.hpp"
0026 #include "Acts/Tests/CommonHelpers/FloatComparisons.hpp"
0027 #include "Acts/Utilities/Logger.hpp"
0028 #include "Acts/Utilities/TrackHelpers.hpp"
0029 #include "Acts/Utilities/Zip.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 BOOST_AUTO_TEST_SUITE(EDM4hepParameterConversion)
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[Acts::eBoundTheta]);
0052   double omega =
0053       par[Acts::eBoundQOverP] / std::sin(par[Acts::eBoundTheta]) * Bz;
0054
0055   auto J1 = EDM4hepUtil::detail::jacobianToEdm4hep(par[eBoundTheta],
0056                                                    par[eBoundQOverP], Bz);
0057
0058   BoundMatrix cov2 = J1 * cov * J1.transpose();
0059
0060   auto J2 = EDM4hepUtil::detail::jacobianFromEdm4hep(tanLambda, omega, Bz);
0061
0062   BoundMatrix cov3 = J2 * cov2 * J2.transpose();
0063
0064   CHECK_CLOSE_ABS(cov, cov3, 1e-9);
0065 }
0066
0067 BOOST_AUTO_TEST_CASE(ConvertTrackParametersToEdm4hepWithPerigee) {
0068   auto refSurface = Surface::makeShared<PerigeeSurface>(Vector3{50, 30, 20});
0069
0070   BoundVector par;
0071   par << 1_mm, 5_mm, 0, std::numbers::pi / 2., -1 / 1_GeV,
0072       5_ns;  // -> perpendicular to perigee and pointing right, should be PCA
0073
0074   BoundMatrix cov;
0075   cov.setIdentity();
0076   cov(5, 5) = 25_ns;
0077
0078   BoundTrackParameters boundPar{refSurface, par, cov,
0079                                 ParticleHypothesis::pion()};
0080
0081   double Bz = 2_T;
0082
0083   Acts::GeometryContext gctx;
0084
0085   EDM4hepUtil::detail::Parameters converted =
0086       EDM4hepUtil::detail::convertTrackParametersToEdm4hep(gctx, Bz, boundPar);
0087
0088   BOOST_CHECK(converted.covariance.has_value());
0089   BOOST_CHECK(converted.surface);
0090
0091   // input is already on perigee, should not be modified
0092   BOOST_CHECK_EQUAL(par.template head<2>(),
0093                     converted.values.template head<2>());
0094   BOOST_CHECK_EQUAL(
0095       (converted.covariance.value().template topLeftCorner<4, 4>()),
0096       ActsSquareMatrix<4>::Identity());
0097   BOOST_CHECK_GT(converted.covariance.value()(4, 4), 0);
0098   BOOST_CHECK_EQUAL(converted.covariance.value()(5, 5), 25_ns);
0099
0100   // convert back for roundtrip test
0101
0102   BoundTrackParameters roundtripPar =
0103       EDM4hepUtil::detail::convertTrackParametersFromEdm4hep(Bz, converted);
0104
0105   BOOST_CHECK(roundtripPar.parameters().isApprox(boundPar.parameters()));
0106   BOOST_CHECK(roundtripPar.covariance().value().isApprox(
0107       boundPar.covariance().value()));
0108 }
0109
0110 BOOST_AUTO_TEST_CASE(ConvertTrackParametersToEdm4hepWithOutPerigee) {
0111   auto planeSurface =
0112       CurvilinearSurface(Vector3{50, 30, 20}, Vector3{1, 1, 0.3}.normalized())
0113           .planeSurface();
0114
0115   BoundVector par;
0116   par << 1_mm, 5_mm, std::numbers::pi / 4., std::numbers::pi / 2. * 0.9,
0117       -1 / 1_GeV, 5_ns;
0118
0119   BoundMatrix cov;
0120   cov.setIdentity();
0121   cov(5, 5) = 25_ns;
0122
0123   BoundTrackParameters planePar{planeSurface, par, cov,
0124                                 ParticleHypothesis::pion()};
0125
0126   double Bz = 2_T;
0127
0128   Acts::GeometryContext gctx;
0129
0130   EDM4hepUtil::detail::Parameters converted =
0131       EDM4hepUtil::detail::convertTrackParametersToEdm4hep(gctx, Bz, planePar);
0132
0133   BOOST_CHECK(converted.covariance.has_value());
0134   BOOST_CHECK(converted.surface);
0135
0136   // input is not a perigee, so new params should be at 0, 0 on ad-hoc perigee
0137   BOOST_CHECK_EQUAL(converted.values.template head<2>(), (Vector2{0, 0}));
0138   CHECK_CLOSE_ABS(converted.values[2], par[2], 1e-6);
0139
0140   BOOST_CHECK_EQUAL(converted.covariance.value()(0, 0), 1);
0141
0142   BOOST_CHECK_LT(converted.covariance.value()(1, 1), 1.2);
0143   BOOST_CHECK_GT(converted.covariance.value()(1, 1), 1);
0144
0145   CHECK_CLOSE_ABS(converted.covariance.value()(2, 2), 1, 1e-6);
0146
0147   BOOST_CHECK_GT(converted.covariance.value()(3, 3), 1);
0148   BOOST_CHECK_LT(converted.covariance.value()(3, 3), 1.2);
0149
0150   BOOST_CHECK_GT(converted.covariance.value()(4, 4), 0);
0151   BOOST_CHECK_EQUAL(converted.covariance.value()(5, 5), 25_ns);
0152
0153   // convert back for roundtrip test
0154   BoundTrackParameters roundtripPar =
0155       EDM4hepUtil::detail::convertTrackParametersFromEdm4hep(Bz, converted);
0156
0157   BOOST_CHECK_NE(dynamic_cast<const Acts::PerigeeSurface*>(
0158                      &roundtripPar.referenceSurface()),
0159                  nullptr);
0160
0161   BOOST_CHECK((converted.covariance.value().topLeftCorner<3, 3>().isApprox(
0162       roundtripPar.covariance().value().topLeftCorner<3, 3>())));
0163   CHECK_CLOSE_ABS(roundtripPar.covariance().value()(3, 3), 1, 1e-6);
0164   CHECK_CLOSE_ABS(roundtripPar.covariance().value()(4, 4), 1, 1e-6);
0165   BOOST_CHECK_EQUAL(roundtripPar.covariance().value()(5, 5), 25_ns);
0166
0167   auto roundtripPlaneBoundParams =
0168       Acts::detail::boundToBoundConversion(gctx, roundtripPar, *planeSurface)
0169           .value();
0170
0171   BOOST_CHECK(roundtripPlaneBoundParams.parameters().isApprox(par));
0172
0173   CHECK_CLOSE_COVARIANCE(roundtripPlaneBoundParams.covariance().value(),
0174                          planePar.covariance().value(), 1e-3);
0175 }
0176
0177 BOOST_AUTO_TEST_CASE(ConvertTrackParametersToEdm4hepWithPerigeeNoCov) {
0178   auto refSurface = Surface::makeShared<PerigeeSurface>(Vector3{50, 30, 20});
0179
0180   BoundVector par;
0181   par << 1_mm, 5_mm, 0, std::numbers::pi / 2., -1 / 1_GeV,
0182       5_ns;  // -> perpendicular to perigee and pointing right, should be PCA
0183
0184   BoundTrackParameters boundPar{refSurface, par, std::nullopt,
0185                                 ParticleHypothesis::pion()};
0186
0187   double Bz = 2_T;
0188
0189   Acts::GeometryContext gctx;
0190
0191   EDM4hepUtil::detail::Parameters converted =
0192       EDM4hepUtil::detail::convertTrackParametersToEdm4hep(gctx, Bz, boundPar);
0193
0194   BOOST_CHECK(!converted.covariance.has_value());
0195   BOOST_CHECK(converted.surface);
0196
0197   // input is already on perigee, should not be modified
0198   BOOST_CHECK_EQUAL(par.template head<2>(),
0199                     converted.values.template head<2>());
0200
0201   // convert back for roundtrip test
0202
0203   BoundTrackParameters roundtripPar =
0204       EDM4hepUtil::detail::convertTrackParametersFromEdm4hep(Bz, converted);
0205
0206   BOOST_CHECK(roundtripPar.parameters().isApprox(boundPar.parameters()));
0207   BOOST_CHECK(!roundtripPar.covariance().has_value());
0208 }
0209
0210 BOOST_AUTO_TEST_CASE(ConvertTrackParametersToEdm4hepWithOutPerigeeNoCov) {
0211   auto refSurface =
0212       CurvilinearSurface(Vector3{50, 30, 20}, Vector3{1, 1, 0.3}.normalized())
0213           .planeSurface();
0214
0215   BoundVector par;
0216   par << 1_mm, 5_mm, std::numbers::pi / 4., std::numbers::pi / 2., -1 / 1_GeV,
0217       5_ns;
0218
0219   BoundTrackParameters boundPar{refSurface, par, std::nullopt,
0220                                 ParticleHypothesis::pion()};
0221
0222   double Bz = 2_T;
0223
0224   Acts::GeometryContext gctx;
0225
0226   EDM4hepUtil::detail::Parameters converted =
0227       EDM4hepUtil::detail::convertTrackParametersToEdm4hep(gctx, Bz, boundPar);
0228
0229   BOOST_CHECK(!converted.covariance.has_value());
0230   BOOST_CHECK(converted.surface);
0231
0232   // input is not a perigee, so new params should be at 0, 0 on ad-hoc perigee
0233   BOOST_CHECK_EQUAL(converted.values.template head<2>(), (Vector2{0, 0}));
0234   CHECK_CLOSE_ABS(converted.values[2], par[2], 1e-6);
0235
0236   // convert back for roundtrip test
0237   BoundTrackParameters roundtripPar =
0238       EDM4hepUtil::detail::convertTrackParametersFromEdm4hep(Bz, converted);
0239
0240   BOOST_CHECK_EQUAL(roundtripPar.parameters().template head<2>(),
0241                     (Vector2{0, 0}));
0242   BOOST_CHECK(roundtripPar.parameters().tail<4>().isApprox(par.tail<4>()));
0243   BOOST_CHECK(!roundtripPar.covariance().has_value());
0244 }
0245
0246 BOOST_AUTO_TEST_CASE(CovariancePacking) {
0247   BoundMatrix m;
0248   // clang-format off
0249   m << 1, 2, 3, 4, 5, 6,
0250        2, 2, 3, 4, 5, 6,
0251        3, 3, 3, 4, 5, 6,
0252        4, 4, 4, 4, 5, 6,
0253        5, 5, 5, 5, 5, 6,
0254        6, 6, 6, 6, 6, 6;
0255   // clang-format on
0256
0257   std::array<float, 21> values{};
0258   EDM4hepUtil::detail::packCovariance(m, values.data());
0259
0260   BoundMatrix m2;
0261   m2.setZero();
0262   EDM4hepUtil::detail::unpackCovariance(values.data(), m2);
0263
0264   CHECK_CLOSE_ABS(m, m2, 1e-9);
0265 }
0266
0267 BOOST_AUTO_TEST_CASE(RoundTripTests) {
0268   auto trackContainer = std::make_shared<Acts::VectorTrackContainer>();
0269   auto trackStateContainer = std::make_shared<Acts::VectorMultiTrajectory>();
0270   TrackContainer tracks(trackContainer, trackStateContainer);
0271
0272   using const_proxy_t = decltype(tracks)::ConstTrackProxy;
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         Acts::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(Acts::Surface::makeShared<PerigeeSurface>(pos));
0342     }
0343
0344     calculateTrackQuantities(track);
0345   }
0346
0347   edm4hep::TrackCollection edm4hepTracks;
0348
0349   Acts::GeometryContext gctx;
0350
0351   double Bz = 3_T;
0352
0353   auto logger = getDefaultLogger("EDM4hep", Logging::INFO);
0354
0355   for (const_proxy_t 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<Acts::VectorTrackContainer>(),
0374                             std::make_shared<Acts::VectorMultiTrajectory>());
0375
0376   for (const auto edm4hepTrack : edm4hepTracksConst) {
0377     EDM4hepUtil::readTrack(edm4hepTrack, readTracks.makeTrack(), Bz, *logger);
0378   }
0379
0380   BOOST_CHECK_EQUAL(tracks.size(), readTracks.size());
0381   std::size_t t = 0;
0382
0383   auto origTrackIt = tracks.begin();
0384   auto readTrackIt = readTracks.begin();
0385   while (origTrackIt != tracks.end() && readTrackIt != readTracks.end()) {
0386     BOOST_TEST_INFO_SCOPE("Track #" << t);
0387     auto orig = *origTrackIt;
0388     auto read = *readTrackIt;
0389
0390     CHECK_CLOSE_OR_SMALL(orig.parameters(), read.parameters(), 1e-5, 1e-8);
0391     CHECK_CLOSE_OR_SMALL(orig.covariance(), read.covariance(), 1e-5, 1e-8);
0392     BOOST_CHECK_EQUAL(orig.referenceSurface().center(gctx),
0393                       read.referenceSurface().center(gctx));
0394
0395     auto origTsIt = orig.trackStatesReversed().begin();
0396     auto readTsIt = read.trackStatesReversed().begin();
0397
0398     std::size_t tsi = 0;
0399     while (origTsIt != orig.trackStatesReversed().end() &&
0400            readTsIt != read.trackStatesReversed().end()) {
0401       BOOST_TEST_INFO_SCOPE("TS: #" << tsi);
0402       auto nextMeas = std::find_if(
0403           origTsIt, orig.trackStatesReversed().end(), [](const auto& ts) {
0404             return ts.typeFlags().test(TrackStateFlag::MeasurementFlag);
0405           });
0406       BOOST_CHECK(nextMeas != orig.trackStatesReversed().end());
0407       origTsIt = nextMeas;
0408       auto origTs = *origTsIt;
0409       auto readTs = *readTsIt;
0410
0411       BOOST_TEST_INFO_SCOPE(
0412           "orig parameters: " << origTs.parameters().transpose());
0413       BOOST_TEST_INFO_SCOPE(
0414           "read parameters: " << readTs.parameters().transpose());
0415       CHECK_CLOSE_OR_SMALL(origTs.smoothedCovariance(),
0416                            readTs.smoothedCovariance(), 1e-5, 1e-6);
0417       Vector3 newCenter = readTs.referenceSurface().center(
0418           gctx);  // new center is a perigee, but should be on the other
0419       // surface
0420       BOOST_CHECK(origTs.referenceSurface().isOnSurface(gctx, newCenter,
0421                                                         Vector3::Zero()));
0422
0423       // global hit positions should be the same
0424       Vector3 readGlobal = readTs.referenceSurface().localToGlobal(
0425           gctx, readTs.parameters().template head<2>(), Vector3::Zero());
0426       Vector3 origGlobal = origTs.referenceSurface().localToGlobal(
0427           gctx, origTs.parameters().template head<2>(), Vector3::Zero());
0428       CHECK_CLOSE_ABS(readGlobal, origGlobal, 1e-3);
0429       ++origTsIt;
0430       ++readTsIt;
0431       tsi++;
0432     }
0433     ++origTrackIt;
0434     ++readTrackIt;
0435
0436     t++;
0437   }
0438 }
0439
0440 BOOST_AUTO_TEST_SUITE_END()