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 // readTrackState is the inverse of writeTrackState for a single state, with
0445 // the field sampled at the same place on both sides.
0446 BOOST_AUTO_TEST_CASE(TrackStateRoundTrip) {
0447   auto gctx = GeometryContext::dangerouslyDefaultConstruct();
0448   MagneticFieldContext mctx{};
0449   ZGradientField field{2_T};
0450
0451   const Vector3 refPoint{1_mm, -2_mm, 300_mm};
0452   BoundVector par;
0453   par << 50_mm, 300_mm, 0.7, 1.1, 1. / 2_GeV, 0;
0454   BoundMatrix cov = BoundMatrix::Identity() * 1e-6;
0455   BoundTrackParameters params{Surface::makeShared<PerigeeSurface>(refPoint),
0456                               par, cov, ParticleHypothesis::pion()};
0457
0458   edm4hep::TrackState state = EDM4hepUtil::writeTrackState(
0459       gctx, mctx, edm4hep::TrackState::AtIP, params, field);
0460
0461   auto fieldCache = field.makeCache(mctx);
0462   BoundTrackParameters read =
0463       EDM4hepUtil::readTrackState(gctx, state, field, fieldCache);
0464
0465   CHECK_CLOSE_OR_SMALL(read.referenceSurface().center(gctx), refPoint, 1e-6,
0466                        1e-6);
0467   CHECK_CLOSE_OR_SMALL(read.parameters(), par, 1e-5, 1e-8);
0468   BOOST_REQUIRE(read.covariance().has_value());
0469   CHECK_CLOSE_OR_SMALL(read.covariance().value(), cov, 1e-5, 1e-8);
0470
0471   // The overload creating its own field cache gives the same result
0472   BoundTrackParameters readNoCache =
0473       EDM4hepUtil::readTrackState(gctx, mctx, state, field);
0474   BOOST_CHECK_EQUAL(readNoCache.parameters(), read.parameters());
0475   BOOST_CHECK_EQUAL(readNoCache.covariance().value(),
0476                     read.covariance().value());
0477 }
0478
0479 BOOST_AUTO_TEST_CASE(RoundTripTests) {
0480   auto trackContainer = std::make_shared<VectorTrackContainer>();
0481   auto trackStateContainer = std::make_shared<VectorMultiTrajectory>();
0482   TrackContainer tracks(trackContainer, trackStateContainer);
0483
0484   std::mt19937 rng{42};
0485   std::normal_distribution<double> gauss(0., 1.);
0486   std::uniform_real_distribution<double> f(-1, 1);
0487   std::uniform_real_distribution<double> r(0, 1);
0488   std::uniform_int_distribution<std::uint32_t> nTracks(2, 20);
0489   std::uniform_int_distribution<std::uint32_t> nTs(1, 20);
0490   std::uniform_real_distribution<double> phiDist(-std::numbers::pi,
0491                                                  std::numbers::pi);
0492   std::uniform_real_distribution<double> etaDist(-4, 4);
0493   std::uniform_real_distribution<double> ptDist(1_MeV, 10_GeV);
0494   std::uniform_real_distribution<double> qDist(0., 1.);
0495
0496   auto genParams = [&]() -> std::pair<BoundVector, BoundMatrix> {
0497     double d0 = 20_um * gauss(rng);
0498     double z0 = 20_mm * gauss(rng);
0499     double phi = phiDist(rng);
0500     double eta = etaDist(rng);
0501     double theta = 2 * std::atan(exp(-eta));
0502     double pt = ptDist(rng);
0503     double p = pt / std::sin(theta);
0504     double charge = qDist(rng) > 0.5 ? 1. : -1.;
0505     double qop = charge / p;
0506     double t = 5_ns * gauss(rng);
0507
0508     BoundVector par;
0509     par << d0, z0, phi, theta, qop, t;
0510     BoundMatrix cov;
0511     cov = BoundMatrix::Identity();
0512     cov.diagonal() << 20_um * 20_um, 20_mm * 20_mm, 0.1, 0.1, 1_GeV, 25_ns;
0513     return {par, cov};
0514   };
0515
0516   std::uint32_t numT = nTracks(rng);
0517   for (std::uint32_t t = 0; t < numT; t++) {
0518     auto track = tracks.makeTrack();
0519     {
0520       auto [par, cov] = genParams();
0521       track.parameters() = par;
0522       track.covariance() = cov;
0523     }
0524     track.setReferenceSurface(
0525         Surface::makeShared<PerigeeSurface>(Vector3{0, 0, 0}));
0526
0527     std::uint32_t numTs = nTs(rng);
0528     for (std::uint32_t i = 0; i < numTs; i++) {
0529       auto ts = track.appendTrackState(TrackStatePropMask::Smoothed);
0530       double crit = r(rng);
0531       if (crit < 0.1) {
0532         ts.typeFlags().setIsHole();
0533       } else if (crit < 0.2) {
0534         ts.typeFlags().setIsOutlier();
0535       } else if (crit < 0.3) {
0536         ts.typeFlags().setIsSharedHit();
0537       } else if (crit < 0.4) {
0538         ts.typeFlags().setIsMaterial();
0539       } else {
0540         ts.typeFlags().setIsMeasurement();
0541       }
0542
0543       auto [par, cov] = genParams();
0544       ts.smoothed() = par;
0545       ts.smoothedCovariance() = cov;
0546       Vector3 pos;
0547       pos << 1000 * f(rng), 1000 * f(rng), 3000 * f(rng);
0548       ts.setReferenceSurface(Surface::makeShared<PerigeeSurface>(pos));
0549     }
0550
0551     calculateTrackQuantities(track);
0552   }
0553
0554   edm4hep::TrackCollection edm4hepTracks;
0555
0556   auto gctx = GeometryContext::dangerouslyDefaultConstruct();
0557
0558   double Bz = 3_T;
0559
0560   auto logger = getDefaultLogger("EDM4hep", Logging::INFO);
0561
0562   for (const auto& track : tracks) {
0563     auto to = edm4hepTracks.create();
0564     EDM4hepUtil::writeTrack(gctx, track, to, Bz, {}, *logger);
0565   }
0566
0567   BOOST_CHECK_EQUAL(edm4hepTracks.size(), tracks.size());
0568
0569   auto tIt = tracks.begin();
0570   for (auto edm4hepTrack : edm4hepTracks) {
0571     auto track = *tIt;
0572     BOOST_CHECK_EQUAL(track.nMeasurements(),
0573                       edm4hepTrack.trackStates_size() - 1);
0574
0575     ++tIt;
0576   }
0577
0578   const edm4hep::TrackCollection& edm4hepTracksConst = edm4hepTracks;
0579
0580   TrackContainer readTracks(std::make_shared<VectorTrackContainer>(),
0581                             std::make_shared<VectorMultiTrajectory>());
0582
0583   for (const auto edm4hepTrack : edm4hepTracksConst) {
0584     auto track = readTracks.makeTrack();
0585     EDM4hepUtil::readTrack(gctx, edm4hepTrack, track, Bz, *logger);
0586   }
0587
0588   BOOST_CHECK_EQUAL(tracks.size(), readTracks.size());
0589   std::size_t t = 0;
0590
0591   auto origTrackIt = tracks.begin();
0592   auto readTrackIt = readTracks.begin();
0593   while (origTrackIt != tracks.end() && readTrackIt != readTracks.end()) {
0594     BOOST_TEST_INFO_SCOPE("Track #" << t);
0595     auto orig = *origTrackIt;
0596     auto read = *readTrackIt;
0597
0598     CHECK_CLOSE_OR_SMALL(orig.parameters(), read.parameters(), 1e-5, 1e-8);
0599     CHECK_CLOSE_OR_SMALL(orig.covariance(), read.covariance(), 1e-5, 1e-8);
0600     BOOST_CHECK_EQUAL(orig.referenceSurface().center(gctx),
0601                       read.referenceSurface().center(gctx));
0602
0603     // Track summary quantities have to survive the round trip as well
0604     CHECK_CLOSE_OR_SMALL(orig.chi2(), read.chi2(), 1e-5, 1e-8);
0605     BOOST_CHECK_EQUAL(orig.nDoF(), read.nDoF());
0606     BOOST_CHECK_EQUAL(orig.nHoles(), read.nHoles());
0607
0608     auto origTsIt = orig.trackStatesReversed().begin();
0609     auto readTsIt = read.trackStatesReversed().begin();
0610
0611     std::size_t tsi = 0;
0612     while (origTsIt != orig.trackStatesReversed().end() &&
0613            readTsIt != read.trackStatesReversed().end()) {
0614       BOOST_TEST_INFO_SCOPE("TS: #" << tsi);
0615       auto nextMeas = std::find_if(
0616           origTsIt, orig.trackStatesReversed().end(),
0617           [](const auto& ts) { return ts.typeFlags().isMeasurement(); });
0618       BOOST_CHECK(nextMeas != orig.trackStatesReversed().end());
0619       origTsIt = nextMeas;
0620       auto origTs = *origTsIt;
0621       auto readTs = *readTsIt;
0622
0623       BOOST_TEST_INFO_SCOPE(
0624           "orig parameters: " << origTs.parameters().transpose());
0625       BOOST_TEST_INFO_SCOPE(
0626           "read parameters: " << readTs.parameters().transpose());
0627       CHECK_CLOSE_OR_SMALL(origTs.smoothedCovariance(),
0628                            readTs.smoothedCovariance(), 1e-5, 1e-6);
0629       Vector3 newCenter = readTs.referenceSurface().center(
0630           gctx);  // new center is a perigee, but should be on the other
0631       // surface
0632       BOOST_CHECK(origTs.referenceSurface().isOnSurface(gctx, newCenter,
0633                                                         Vector3::Zero()));
0634
0635       // global hit positions should be the same
0636       Vector3 readGlobal = readTs.referenceSurface().localToGlobal(
0637           gctx, readTs.parameters().template head<2>(), Vector3::Zero());
0638       Vector3 origGlobal = origTs.referenceSurface().localToGlobal(
0639           gctx, origTs.parameters().template head<2>(), Vector3::Zero());
0640       CHECK_CLOSE_ABS(readGlobal, origGlobal, 1e-3);
0641       ++origTsIt;
0642       ++readTsIt;
0643       tsi++;
0644     }
0645     ++origTrackIt;
0646     ++readTrackIt;
0647
0648     t++;
0649   }
0650 }
0651
0652 BOOST_AUTO_TEST_SUITE_END()
0653
0654 }  // namespace ActsTests