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 <array>
0030 #include <numbers>
0031 #include <random>
0032
0033 #include <edm4hep/Constants.h>
0034 #include <edm4hep/TrackCollection.h>
0035 #include <edm4hep/TrackState.h>
0036 #include <edm4hep/TrackerHit3DCollection.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   auto gctx = GeometryContext::dangerouslyDefaultConstruct();
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       SquareMatrix<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   auto gctx = GeometryContext::dangerouslyDefaultConstruct();
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   // The time variance grows slightly: the path length from the plane surface
0155   // to the ad-hoc perigee varies with the other (uncertain) parameters, and
0156   // dt/ds couples that variation into the time.
0157   BOOST_CHECK_GT(converted.covariance.value()(5, 5), 25_ns);
0158   CHECK_CLOSE_REL(converted.covariance.value()(5, 5), 25_ns, 1e-4);
0159
0160   // convert back for roundtrip test
0161   BoundTrackParameters roundtripPar =
0162       EDM4hepUtil::detail::convertTrackParametersFromEdm4hep(Bz, converted);
0163
0164   BOOST_CHECK_NE(
0165       dynamic_cast<const PerigeeSurface*>(&roundtripPar.referenceSurface()),
0166       nullptr);
0167
0168   BOOST_CHECK((converted.covariance.value().topLeftCorner<3, 3>().isApprox(
0169       roundtripPar.covariance().value().topLeftCorner<3, 3>())));
0170   CHECK_CLOSE_ABS(roundtripPar.covariance().value()(3, 3), 1, 1e-6);
0171   CHECK_CLOSE_ABS(roundtripPar.covariance().value()(4, 4), 1, 1e-6);
0172   BOOST_CHECK_GT(roundtripPar.covariance().value()(5, 5), 25_ns);
0173   CHECK_CLOSE_REL(roundtripPar.covariance().value()(5, 5), 25_ns, 1e-4);
0174
0175   auto roundtripPlaneBoundParams =
0176       detail::boundToBoundConversion(gctx, roundtripPar, *planeSurface).value();
0177
0178   BOOST_CHECK(roundtripPlaneBoundParams.parameters().isApprox(par));
0179
0180   CHECK_CLOSE_COVARIANCE(roundtripPlaneBoundParams.covariance().value(),
0181                          planePar.covariance().value(), 1e-3);
0182 }
0183
0184 BOOST_AUTO_TEST_CASE(ConvertTrackParametersToEdm4hepWithPerigeeNoCov) {
0185   auto refSurface = Surface::makeShared<PerigeeSurface>(Vector3{50, 30, 20});
0186
0187   BoundVector par;
0188   par << 1_mm, 5_mm, 0, std::numbers::pi / 2., -1 / 1_GeV,
0189       5_ns;  // -> perpendicular to perigee and pointing right, should be PCA
0190
0191   BoundTrackParameters boundPar{refSurface, par, std::nullopt,
0192                                 ParticleHypothesis::pion()};
0193
0194   double Bz = 2_T;
0195
0196   auto gctx = GeometryContext::dangerouslyDefaultConstruct();
0197
0198   EDM4hepUtil::detail::Parameters converted =
0199       EDM4hepUtil::detail::convertTrackParametersToEdm4hep(gctx, Bz, boundPar);
0200
0201   BOOST_CHECK(!converted.covariance.has_value());
0202   BOOST_CHECK(converted.surface);
0203
0204   // input is already on perigee, should not be modified
0205   BOOST_CHECK_EQUAL(par.template head<2>(),
0206                     converted.values.template head<2>());
0207
0208   // convert back for roundtrip test
0209
0210   BoundTrackParameters roundtripPar =
0211       EDM4hepUtil::detail::convertTrackParametersFromEdm4hep(Bz, converted);
0212
0213   BOOST_CHECK(roundtripPar.parameters().isApprox(boundPar.parameters()));
0214   BOOST_CHECK(!roundtripPar.covariance().has_value());
0215 }
0216
0217 BOOST_AUTO_TEST_CASE(ConvertTrackParametersToEdm4hepWithOutPerigeeNoCov) {
0218   std::shared_ptr<PlaneSurface> refSurface =
0219       CurvilinearSurface(Vector3{50, 30, 20}, Vector3{1, 1, 0.3}.normalized())
0220           .planeSurface();
0221
0222   BoundVector par;
0223   par << 1_mm, 5_mm, std::numbers::pi / 4., std::numbers::pi / 2., -1 / 1_GeV,
0224       5_ns;
0225
0226   BoundTrackParameters boundPar{refSurface, par, std::nullopt,
0227                                 ParticleHypothesis::pion()};
0228
0229   double Bz = 2_T;
0230
0231   auto gctx = GeometryContext::dangerouslyDefaultConstruct();
0232
0233   EDM4hepUtil::detail::Parameters converted =
0234       EDM4hepUtil::detail::convertTrackParametersToEdm4hep(gctx, Bz, boundPar);
0235
0236   BOOST_CHECK(!converted.covariance.has_value());
0237   BOOST_CHECK(converted.surface);
0238
0239   // input is not a perigee, so new params should be at 0, 0 on ad-hoc perigee
0240   BOOST_CHECK_EQUAL(converted.values.template head<2>(), (Vector2{0, 0}));
0241   CHECK_CLOSE_ABS(converted.values[2], par[2], 1e-6);
0242
0243   // convert back for roundtrip test
0244   BoundTrackParameters roundtripPar =
0245       EDM4hepUtil::detail::convertTrackParametersFromEdm4hep(Bz, converted);
0246
0247   BOOST_CHECK_EQUAL(roundtripPar.parameters().template head<2>(),
0248                     (Vector2{0, 0}));
0249   BOOST_CHECK(roundtripPar.parameters().tail<4>().isApprox(par.tail<4>()));
0250   BOOST_CHECK(!roundtripPar.covariance().has_value());
0251 }
0252
0253 BOOST_AUTO_TEST_CASE(CovariancePacking) {
0254   BoundMatrix m;
0255   // clang-format off
0256   m << 1, 2, 3, 4, 5, 6,
0257        2, 2, 3, 4, 5, 6,
0258        3, 3, 3, 4, 5, 6,
0259        4, 4, 4, 4, 5, 6,
0260        5, 5, 5, 5, 5, 6,
0261        6, 6, 6, 6, 6, 6;
0262   // clang-format on
0263
0264   edm4hep::CovMatrix6f values{};
0265   EDM4hepUtil::detail::packCovariance(m, values);
0266
0267   BoundMatrix m2;
0268   m2.setZero();
0269   EDM4hepUtil::detail::unpackCovariance(values, m2);
0270
0271   CHECK_CLOSE_ABS(m, m2, 1e-9);
0272 }
0273
0274 // The packed covariance has to follow the EDM4hep parameter order
0275 // (edm4hep::TrackParams: d0, phi, omega, z0, tanLambda, time), which is a
0276 // different permutation from the order used internally by
0277 // detail::Parameters::values (d0, z0, phi, tanLambda, omega, time). A
0278 // pack/unpack round trip cannot catch a mismatch here, since it would be
0279 // symmetric, so check the absolute slots against EDM4hep's own accessor.
0280 BOOST_AUTO_TEST_CASE(CovariancePackingConvention) {
0281   // values index -> edm4hep parameter it represents
0282   constexpr std::array<edm4hep::TrackParams, 6> params{
0283       edm4hep::TrackParams::d0,         // values[0]
0284       edm4hep::TrackParams::z0,         // values[1]
0285       edm4hep::TrackParams::phi,        // values[2]
0286       edm4hep::TrackParams::tanLambda,  // values[3]
0287       edm4hep::TrackParams::omega,      // values[4]
0288       edm4hep::TrackParams::time,       // values[5]
0289   };
0290
0291   // Distinct, symmetric entries so every slot is identifiable
0292   SquareMatrix<6> m;
0293   for (std::size_t i = 0; i < 6; ++i) {
0294     for (std::size_t j = 0; j < 6; ++j) {
0295       m(i, j) = 10. * std::min(i, j) + std::max(i, j);
0296     }
0297   }
0298
0299   edm4hep::TrackState trackState{};
0300   EDM4hepUtil::detail::packCovariance(m, trackState.covMatrix);
0301
0302   for (std::size_t i = 0; i < 6; ++i) {
0303     for (std::size_t j = 0; j < 6; ++j) {
0304       BOOST_TEST_INFO_SCOPE("values index (" << i << ", " << j << ")");
0305       BOOST_CHECK_EQUAL(trackState.getCovMatrix(params.at(i), params.at(j)),
0306                         static_cast<float>(m(i, j)));
0307     }
0308   }
0309 }
0310
0311 // The hit lookup has to be invoked in the order the track states are emitted
0312 // (inside-out), not in the outside-in order they are iterated internally: a
0313 // stateful callback would otherwise associate its hits with the wrong states.
0314 BOOST_AUTO_TEST_CASE(HitLookupInvocationOrder) {
0315   TrackContainer tracks(std::make_shared<VectorTrackContainer>(),
0316                         std::make_shared<VectorMultiTrajectory>());
0317   auto track = tracks.makeTrack();
0318
0319   BoundVector par;
0320   par << 1_mm, 5_mm, 0.1, std::numbers::pi / 2. * 0.9, -1 / 1_GeV, 5_ns;
0321   track.parameters() = par;
0322   track.covariance() = BoundMatrix::Identity();
0323   track.setReferenceSurface(
0324       Surface::makeShared<PerigeeSurface>(Vector3{0, 0, 0}));
0325
0326   // Appended innermost first, so the internal reverse iteration visits them
0327   // outermost first.
0328   const std::array<double, 3> zs{100, 200, 300};
0329   for (double z : zs) {
0330     auto ts = track.appendTrackState(TrackStatePropMask::Smoothed);
0331     ts.typeFlags().setIsMeasurement();
0332     ts.smoothed() = par;
0333     ts.smoothedCovariance() = BoundMatrix::Identity();
0334     ts.setReferenceSurface(
0335         Surface::makeShared<PerigeeSurface>(Vector3{0, 0, z}));
0336   }
0337
0338   auto gctx = GeometryContext::dangerouslyDefaultConstruct();
0339   MagneticFieldContext mctx{};
0340   ConstantBField field{Vector3{0, 0, 2_T}};
0341
0342   edm4hep::TrackerHit3DCollection hitCollection;
0343   std::vector<edm4hep::TrackerHit3D> hitByZ;
0344   for (double z : zs) {
0345     auto hit = hitCollection.create();
0346     hit.setPosition({0, 0, z});
0347     hitByZ.push_back(hit);
0348   }
0349
0350   std::vector<double> seen;
0351   EDM4hepUtil::TrackerHitLookup lookup =
0352       [&](const AnyConstTrackStateProxy& state)
0353       -> std::optional<edm4hep::TrackerHit> {
0354     double z = state.referenceSurface().center(gctx).z();
0355     seen.push_back(z);
0356     for (std::size_t i = 0; i < zs.size(); ++i) {
0357       if (std::abs(z - zs[i]) < 1e-6) {
0358         return hitByZ[i];
0359       }
0360     }
0361     return std::nullopt;
0362   };
0363
0364   edm4hep::TrackCollection outTracks;
0365   auto to = outTracks.create();
0366   EDM4hepUtil::writeTrack(gctx, mctx, track, to, field, lookup);
0367
0368   // Invoked once per measurement state, inside-out
0369   BOOST_REQUIRE_EQUAL(seen.size(), zs.size());
0370   for (std::size_t i = 0; i < zs.size(); ++i) {
0371     BOOST_TEST_INFO_SCOPE("lookup call #" << i);
0372     CHECK_CLOSE_ABS(seen[i], zs[i], 1e-6);
0373   }
0374
0375   // ... and the hits are attached in that same order
0376   BOOST_REQUIRE_EQUAL(to.trackerHits_size(), zs.size());
0377   for (std::size_t i = 0; i < zs.size(); ++i) {
0378     BOOST_TEST_INFO_SCOPE("tracker hit #" << i);
0379     CHECK_CLOSE_ABS(to.getTrackerHits(i).getPosition().z, zs[i], 1e-6);
0380   }
0381 }
0382
0383 namespace {
0384 /// Bz grows with |z|, so the field at a track's point of closest approach
0385 /// differs from the field at its perigee reference point.
0386 class ZGradientField : public MagneticFieldProvider {
0387  public:
0388   explicit ZGradientField(double b0) : m_b0(b0) {}
0389
0390   Cache makeCache(const MagneticFieldContext& /*mctx*/) const override {
0391     return Cache(std::in_place_type<int>, 0);
0392   }
0393
0394   Result<Vector3> getField(const Vector3& position,
0395                            Cache& /*cache*/) const override {
0396     return Result<Vector3>::success(
0397         Vector3{0, 0, m_b0 * (1 + 0.001 * std::abs(position.z()))});
0398   }
0399
0400  private:
0401   double m_b0;
0402 };
0403 }  // namespace
0404
0405 // The field is sampled at the point of closest approach, which for parameters
0406 // already on a perigee is d0/z0 away from the reference point. Both sides have
0407 // to sample at the same place, or q/p comes back biased in a non-uniform field.
0408 BOOST_AUTO_TEST_CASE(NonUniformFieldRoundTrip) {
0409   auto gctx = GeometryContext::dangerouslyDefaultConstruct();
0410   MagneticFieldContext mctx{};
0411   ZGradientField field{2_T};
0412
0413   BoundVector par;
0414   par << 50_mm, 300_mm, 0.7, 1.1, 1. / 2_GeV, 0;
0415   BoundMatrix cov = BoundMatrix::Identity() * 1e-6;
0416
0417   TrackContainer tracks(std::make_shared<VectorTrackContainer>(),
0418                         std::make_shared<VectorMultiTrajectory>());
0419   auto track = tracks.makeTrack();
0420   track.parameters() = par;
0421   track.covariance() = cov;
0422   track.setReferenceSurface(
0423       Surface::makeShared<PerigeeSurface>(Vector3{0, 0, 0}));
0424
0425   auto ts = track.appendTrackState(TrackStatePropMask::Smoothed);
0426   ts.typeFlags().setIsMeasurement();
0427   ts.smoothed() = par;
0428   ts.smoothedCovariance() = cov;
0429   ts.setReferenceSurface(
0430       Surface::makeShared<PerigeeSurface>(Vector3{0, 0, 300_mm}));
0431
0432   edm4hep::TrackCollection edm4hepTracks;
0433   auto to = edm4hepTracks.create();
0434   EDM4hepUtil::writeTrack(gctx, mctx, track, to, field);
0435
0436   TrackContainer readTracks(std::make_shared<VectorTrackContainer>(),
0437                             std::make_shared<VectorMultiTrajectory>());
0438   auto read = readTracks.makeTrack();
0439   EDM4hepUtil::readTrack(gctx, mctx, to, read, field);
0440
0441   CHECK_CLOSE_OR_SMALL(track.parameters(), read.parameters(), 1e-5, 1e-8);
0442 }
0443
0444 BOOST_AUTO_TEST_CASE(RoundTripTests) {
0445   auto trackContainer = std::make_shared<VectorTrackContainer>();
0446   auto trackStateContainer = std::make_shared<VectorMultiTrajectory>();
0447   TrackContainer tracks(trackContainer, trackStateContainer);
0448
0449   std::mt19937 rng{42};
0450   std::normal_distribution<double> gauss(0., 1.);
0451   std::uniform_real_distribution<double> f(-1, 1);
0452   std::uniform_real_distribution<double> r(0, 1);
0453   std::uniform_int_distribution<std::uint32_t> nTracks(2, 20);
0454   std::uniform_int_distribution<std::uint32_t> nTs(1, 20);
0455   std::uniform_real_distribution<double> phiDist(-std::numbers::pi,
0456                                                  std::numbers::pi);
0457   std::uniform_real_distribution<double> etaDist(-4, 4);
0458   std::uniform_real_distribution<double> ptDist(1_MeV, 10_GeV);
0459   std::uniform_real_distribution<double> qDist(0., 1.);
0460
0461   auto genParams = [&]() -> std::pair<BoundVector, BoundMatrix> {
0462     double d0 = 20_um * gauss(rng);
0463     double z0 = 20_mm * gauss(rng);
0464     double phi = phiDist(rng);
0465     double eta = etaDist(rng);
0466     double theta = 2 * std::atan(exp(-eta));
0467     double pt = ptDist(rng);
0468     double p = pt / std::sin(theta);
0469     double charge = qDist(rng) > 0.5 ? 1. : -1.;
0470     double qop = charge / p;
0471     double t = 5_ns * gauss(rng);
0472
0473     BoundVector par;
0474     par << d0, z0, phi, theta, qop, t;
0475     BoundMatrix cov;
0476     cov = BoundMatrix::Identity();
0477     cov.diagonal() << 20_um * 20_um, 20_mm * 20_mm, 0.1, 0.1, 1_GeV, 25_ns;
0478     return {par, cov};
0479   };
0480
0481   std::uint32_t numT = nTracks(rng);
0482   for (std::uint32_t t = 0; t < numT; t++) {
0483     auto track = tracks.makeTrack();
0484     {
0485       auto [par, cov] = genParams();
0486       track.parameters() = par;
0487       track.covariance() = cov;
0488     }
0489     track.setReferenceSurface(
0490         Surface::makeShared<PerigeeSurface>(Vector3{0, 0, 0}));
0491
0492     std::uint32_t numTs = nTs(rng);
0493     for (std::uint32_t i = 0; i < numTs; i++) {
0494       auto ts = track.appendTrackState(TrackStatePropMask::Smoothed);
0495       double crit = r(rng);
0496       if (crit < 0.1) {
0497         ts.typeFlags().setIsHole();
0498       } else if (crit < 0.2) {
0499         ts.typeFlags().setIsOutlier();
0500       } else if (crit < 0.3) {
0501         ts.typeFlags().setIsSharedHit();
0502       } else if (crit < 0.4) {
0503         ts.typeFlags().setIsMaterial();
0504       } else {
0505         ts.typeFlags().setIsMeasurement();
0506       }
0507
0508       auto [par, cov] = genParams();
0509       ts.smoothed() = par;
0510       ts.smoothedCovariance() = cov;
0511       Vector3 pos;
0512       pos << 1000 * f(rng), 1000 * f(rng), 3000 * f(rng);
0513       ts.setReferenceSurface(Surface::makeShared<PerigeeSurface>(pos));
0514     }
0515
0516     calculateTrackQuantities(track);
0517   }
0518
0519   edm4hep::TrackCollection edm4hepTracks;
0520
0521   auto gctx = GeometryContext::dangerouslyDefaultConstruct();
0522
0523   double Bz = 3_T;
0524
0525   auto logger = getDefaultLogger("EDM4hep", Logging::INFO);
0526
0527   for (const auto& track : tracks) {
0528     auto to = edm4hepTracks.create();
0529     EDM4hepUtil::writeTrack(gctx, track, to, Bz, {}, *logger);
0530   }
0531
0532   BOOST_CHECK_EQUAL(edm4hepTracks.size(), tracks.size());
0533
0534   auto tIt = tracks.begin();
0535   for (auto edm4hepTrack : edm4hepTracks) {
0536     auto track = *tIt;
0537     BOOST_CHECK_EQUAL(track.nMeasurements(),
0538                       edm4hepTrack.trackStates_size() - 1);
0539
0540     ++tIt;
0541   }
0542
0543   const edm4hep::TrackCollection& edm4hepTracksConst = edm4hepTracks;
0544
0545   TrackContainer readTracks(std::make_shared<VectorTrackContainer>(),
0546                             std::make_shared<VectorMultiTrajectory>());
0547
0548   for (const auto edm4hepTrack : edm4hepTracksConst) {
0549     auto track = readTracks.makeTrack();
0550     EDM4hepUtil::readTrack(gctx, edm4hepTrack, track, Bz, *logger);
0551   }
0552
0553   BOOST_CHECK_EQUAL(tracks.size(), readTracks.size());
0554   std::size_t t = 0;
0555
0556   auto origTrackIt = tracks.begin();
0557   auto readTrackIt = readTracks.begin();
0558   while (origTrackIt != tracks.end() && readTrackIt != readTracks.end()) {
0559     BOOST_TEST_INFO_SCOPE("Track #" << t);
0560     auto orig = *origTrackIt;
0561     auto read = *readTrackIt;
0562
0563     CHECK_CLOSE_OR_SMALL(orig.parameters(), read.parameters(), 1e-5, 1e-8);
0564     CHECK_CLOSE_OR_SMALL(orig.covariance(), read.covariance(), 1e-5, 1e-8);
0565     BOOST_CHECK_EQUAL(orig.referenceSurface().center(gctx),
0566                       read.referenceSurface().center(gctx));
0567
0568     // Track summary quantities have to survive the round trip as well
0569     CHECK_CLOSE_OR_SMALL(orig.chi2(), read.chi2(), 1e-5, 1e-8);
0570     BOOST_CHECK_EQUAL(orig.nDoF(), read.nDoF());
0571     BOOST_CHECK_EQUAL(orig.nHoles(), read.nHoles());
0572
0573     auto origTsIt = orig.trackStatesReversed().begin();
0574     auto readTsIt = read.trackStatesReversed().begin();
0575
0576     std::size_t tsi = 0;
0577     while (origTsIt != orig.trackStatesReversed().end() &&
0578            readTsIt != read.trackStatesReversed().end()) {
0579       BOOST_TEST_INFO_SCOPE("TS: #" << tsi);
0580       auto nextMeas = std::find_if(
0581           origTsIt, orig.trackStatesReversed().end(),
0582           [](const auto& ts) { return ts.typeFlags().isMeasurement(); });
0583       BOOST_CHECK(nextMeas != orig.trackStatesReversed().end());
0584       origTsIt = nextMeas;
0585       auto origTs = *origTsIt;
0586       auto readTs = *readTsIt;
0587
0588       BOOST_TEST_INFO_SCOPE(
0589           "orig parameters: " << origTs.parameters().transpose());
0590       BOOST_TEST_INFO_SCOPE(
0591           "read parameters: " << readTs.parameters().transpose());
0592       CHECK_CLOSE_OR_SMALL(origTs.smoothedCovariance(),
0593                            readTs.smoothedCovariance(), 1e-5, 1e-6);
0594       Vector3 newCenter = readTs.referenceSurface().center(
0595           gctx);  // new center is a perigee, but should be on the other
0596       // surface
0597       BOOST_CHECK(origTs.referenceSurface().isOnSurface(gctx, newCenter,
0598                                                         Vector3::Zero()));
0599
0600       // global hit positions should be the same
0601       Vector3 readGlobal = readTs.referenceSurface().localToGlobal(
0602           gctx, readTs.parameters().template head<2>(), Vector3::Zero());
0603       Vector3 origGlobal = origTs.referenceSurface().localToGlobal(
0604           gctx, origTs.parameters().template head<2>(), Vector3::Zero());
0605       CHECK_CLOSE_ABS(readGlobal, origGlobal, 1e-3);
0606       ++origTsIt;
0607       ++readTsIt;
0608       tsi++;
0609     }
0610     ++origTrackIt;
0611     ++readTrackIt;
0612
0613     t++;
0614   }
0615 }
0616
0617 BOOST_AUTO_TEST_SUITE_END()
0618
0619 }  // namespace ActsTests