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   std::shared_ptr<PlaneSurface> 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   std::shared_ptr<PlaneSurface> 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   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(-std::numbers::pi,
0279                                                  std::numbers::pi);
0280   std::uniform_real_distribution<double> etaDist(-4, 4);
0281   std::uniform_real_distribution<double> ptDist(1_MeV, 10_GeV);
0282   std::uniform_real_distribution<double> qDist(0., 1.);
0283
0284   auto genParams = [&]() -> std::pair<BoundVector, BoundMatrix> {
0285     double d0 = 20_um * gauss(rng);
0286     double z0 = 20_mm * gauss(rng);
0287     double phi = phiDist(rng);
0288     double eta = etaDist(rng);
0289     double theta = 2 * atan(exp(-eta));
0290     double pt = ptDist(rng);
0291     double p = pt / sin(theta);
0292     double charge = qDist(rng) > 0.5 ? 1. : -1.;
0293     double qop = charge / p;
0294     double t = 5_ns * gauss(rng);
0295
0296     BoundVector par;
0297     par << d0, z0, phi, theta, qop, t;
0298     BoundMatrix cov;
0299     cov = BoundMatrix::Identity();
0300     cov.diagonal() << 20_um * 20_um, 20_mm * 20_mm, 0.1, 0.1, 1_GeV, 25_ns;
0301     return {par, cov};
0302   };
0303
0304   std::uint32_t numT = nTracks(rng);
0305   for (std::uint32_t t = 0; t < numT; t++) {
0306     auto track = tracks.makeTrack();
0307     {
0308       auto [par, cov] = genParams();
0309       track.parameters() = par;
0310       track.covariance() = cov;
0311     }
0312     track.setReferenceSurface(
0313         Acts::Surface::makeShared<PerigeeSurface>(Vector3{0, 0, 0}));
0314
0315     std::uint32_t numTs = nTs(rng);
0316     for (std::uint32_t i = 0; i < numTs; i++) {
0317       auto ts = track.appendTrackState(TrackStatePropMask::Smoothed);
0318       double crit = r(rng);
0319       if (crit < 0.1) {
0320         ts.typeFlags().set(TrackStateFlag::HoleFlag);
0321         continue;
0322       } else if (crit < 0.2) {
0323         ts.typeFlags().set(TrackStateFlag::OutlierFlag);
0324         continue;
0325       } else if (crit < 0.3) {
0326         ts.typeFlags().set(TrackStateFlag::SharedHitFlag);
0327       } else if (crit < 0.4) {
0328         ts.typeFlags().set(TrackStateFlag::MaterialFlag);
0329         continue;
0330       }
0331
0332       ts.typeFlags().set(TrackStateFlag::MeasurementFlag);
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(Acts::Surface::makeShared<PerigeeSurface>(pos));
0340     }
0341
0342     calculateTrackQuantities(track);
0343   }
0344
0345   edm4hep::TrackCollection edm4hepTracks;
0346
0347   Acts::GeometryContext gctx;
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<Acts::VectorTrackContainer>(),
0372                             std::make_shared<Acts::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(), [](const auto& ts) {
0403             return ts.typeFlags().test(TrackStateFlag::MeasurementFlag);
0404           });
0405       BOOST_CHECK(nextMeas != orig.trackStatesReversed().end());
0406       origTsIt = nextMeas;
0407       auto origTs = *origTsIt;
0408       auto readTs = *readTsIt;
0409
0410       BOOST_TEST_INFO_SCOPE(
0411           "orig parameters: " << origTs.parameters().transpose());
0412       BOOST_TEST_INFO_SCOPE(
0413           "read parameters: " << readTs.parameters().transpose());
0414       CHECK_CLOSE_OR_SMALL(origTs.smoothedCovariance(),
0415                            readTs.smoothedCovariance(), 1e-5, 1e-6);
0416       Vector3 newCenter = readTs.referenceSurface().center(
0417           gctx);  // new center is a perigee, but should be on the other
0418       // surface
0419       BOOST_CHECK(origTs.referenceSurface().isOnSurface(gctx, newCenter,
0420                                                         Vector3::Zero()));
0421
0422       // global hit positions should be the same
0423       Vector3 readGlobal = readTs.referenceSurface().localToGlobal(
0424           gctx, readTs.parameters().template head<2>(), Vector3::Zero());
0425       Vector3 origGlobal = origTs.referenceSurface().localToGlobal(
0426           gctx, origTs.parameters().template head<2>(), Vector3::Zero());
0427       CHECK_CLOSE_ABS(readGlobal, origGlobal, 1e-3);
0428       ++origTsIt;
0429       ++readTsIt;
0430       tsi++;
0431     }
0432     ++origTrackIt;
0433     ++readTrackIt;
0434
0435     t++;
0436   }
0437 }
0438
0439 BOOST_AUTO_TEST_SUITE_END()