Analysis Software
Documentation for sPHENIX simulation software
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
ConvertTrackEDM4hepTest.cpp
Go to the documentation of this file. Or view the newest version in sPHENIX GitHub for file ConvertTrackEDM4hepTest.cpp
1 // This file is part of the Acts project.
2 //
3 // Copyright (C) 2023 CERN for the benefit of the Acts project
4 //
5 // This Source Code Form is subject to the terms of the Mozilla Public
6 // License, v. 2.0. If a copy of the MPL was not distributed with this
7 // file, You can obtain one at http://mozilla.org/MPL/2.0/.
8 
9 #include <boost/test/data/test_case.hpp>
10 #include <boost/test/tools/old/interface.hpp>
11 #include <boost/test/unit_test.hpp>
12 #include <boost/test/unit_test_suite.hpp>
13 
14 #include "Acts/Definitions/TrackParametrization.hpp"
15 #include "Acts/EventData/Charge.hpp"
16 #include "Acts/EventData/GenericBoundTrackParameters.hpp"
17 #include "Acts/EventData/MultiTrajectory.hpp"
18 #include "Acts/EventData/TrackHelpers.hpp"
19 #include "Acts/EventData/TrackStatePropMask.hpp"
20 #include "Acts/EventData/VectorMultiTrajectory.hpp"
21 #include "Acts/EventData/VectorTrackContainer.hpp"
22 #include "Acts/Geometry/GeometryContext.hpp"
23 #include "Acts/Plugins/EDM4hep/EDM4hepUtil.hpp"
24 #include "Acts/Surfaces/PerigeeSurface.hpp"
25 #include "Acts/Surfaces/PlaneSurface.hpp"
26 #include "Acts/Surfaces/Surface.hpp"
27 #include "Acts/Tests/CommonHelpers/FloatComparisons.hpp"
28 #include "Acts/Utilities/Logger.hpp"
29 #include "Acts/Utilities/Zip.hpp"
30 
31 #include <algorithm>
32 #include <random>
33 
34 #include <edm4hep/TrackCollection.h>
35 
36 using namespace Acts;
37 using namespace Acts::UnitLiterals;
38 BOOST_AUTO_TEST_SUITE(EDM4hepParameterConversion)
39 
40 BOOST_AUTO_TEST_CASE(ConvertTrackParametersToEdm4hepWithPerigee) {
41  auto refSurface = Surface::makeShared<PerigeeSurface>(Vector3{50, 30, 20});
42 
43  BoundVector par;
44  par << 1_mm, 5_mm, 0, M_PI_2, -1 / 1_GeV,
45  5_ns; // -> perpendicular to perigee and pointing right, should be PCA
46 
47  BoundMatrix cov;
48  cov.setIdentity();
49  cov(5, 5) = 25_ns;
50 
51  BoundTrackParameters boundPar{refSurface, par, cov,
52  ParticleHypothesis::pion()};
53 
54  double Bz = 2_T;
55 
56  Acts::GeometryContext gctx;
57 
58  EDM4hepUtil::detail::Parameters converted =
59  EDM4hepUtil::detail::convertTrackParametersToEdm4hep(gctx, Bz, boundPar);
60 
61  BOOST_CHECK(converted.covariance.has_value());
62  BOOST_CHECK(converted.surface);
63 
64  // input is already on perigee, should not be modified
65  BOOST_CHECK_EQUAL(par.template head<2>(),
66  converted.values.template head<2>());
67  BOOST_CHECK_EQUAL(
68  (converted.covariance.value().template topLeftCorner<4, 4>()),
69  ActsSquareMatrix<4>::Identity());
70  BOOST_CHECK(converted.covariance.value()(4, 4) > 0);
71  BOOST_CHECK_EQUAL(converted.covariance.value()(5, 5), 25_ns);
72 
73  // convert back for roundtrip test
74 
75  BoundTrackParameters roundtripPar =
76  EDM4hepUtil::detail::convertTrackParametersFromEdm4hep(Bz, converted);
77 
78  BOOST_CHECK(roundtripPar.parameters().isApprox(boundPar.parameters()));
79  BOOST_CHECK(roundtripPar.covariance().value().isApprox(
80  boundPar.covariance().value()));
81 }
82 
83 BOOST_AUTO_TEST_CASE(ConvertTrackParametersToEdm4hepWithOutPerigee) {
84  auto refSurface = Surface::makeShared<PlaneSurface>(
85  Vector3{50, 30, 20}, Vector3{1, 1, 0.3}.normalized());
86 
87  BoundVector par;
88  par << 1_mm, 5_mm, M_PI / 4., M_PI_2, -1 / 1_GeV, 5_ns;
89 
90  BoundMatrix cov;
91  cov.setIdentity();
92  cov(5, 5) = 25_ns;
93 
94  BoundTrackParameters boundPar{refSurface, par, cov,
95  ParticleHypothesis::pion()};
96 
97  double Bz = 2_T;
98 
99  Acts::GeometryContext gctx;
100 
101  EDM4hepUtil::detail::Parameters converted =
102  EDM4hepUtil::detail::convertTrackParametersToEdm4hep(gctx, Bz, boundPar);
103 
104  BOOST_CHECK(converted.covariance.has_value());
105  BOOST_CHECK(converted.surface);
106 
107  // input is not a perigee, so new params should be at 0, 0 on ad-hoc perigee
108  BOOST_CHECK_EQUAL(converted.values.template head<2>(), (Vector2{0, 0}));
109  CHECK_CLOSE_ABS(converted.values[2], par[2], 1e-6);
110 
111  BOOST_CHECK((converted.covariance.value().template topLeftCorner<4, 4>())
112  .isApprox(ActsSquareMatrix<4>::Identity()));
113  BOOST_CHECK(converted.covariance.value()(4, 4) > 0);
114  BOOST_CHECK_EQUAL(converted.covariance.value()(5, 5), 25_ns);
115 
116  // convert back for roundtrip test
117  BoundTrackParameters roundtripPar =
118  EDM4hepUtil::detail::convertTrackParametersFromEdm4hep(Bz, converted);
119 
120  BOOST_CHECK_EQUAL(roundtripPar.parameters().template head<2>(),
121  (Vector2{0, 0}));
122  BOOST_CHECK(roundtripPar.parameters().tail<4>().isApprox(par.tail<4>()));
123  BOOST_CHECK(roundtripPar.covariance().value().isApprox(
124  boundPar.covariance().value()));
125 }
126 
127 BOOST_AUTO_TEST_CASE(ConvertTrackParametersToEdm4hepWithPerigeeNoCov) {
128  auto refSurface = Surface::makeShared<PerigeeSurface>(Vector3{50, 30, 20});
129 
130  BoundVector par;
131  par << 1_mm, 5_mm, 0, M_PI_2, -1 / 1_GeV,
132  5_ns; // -> perpendicular to perigee and pointing right, should be PCA
133 
134  BoundTrackParameters boundPar{refSurface, par, std::nullopt,
135  ParticleHypothesis::pion()};
136 
137  double Bz = 2_T;
138 
139  Acts::GeometryContext gctx;
140 
141  EDM4hepUtil::detail::Parameters converted =
142  EDM4hepUtil::detail::convertTrackParametersToEdm4hep(gctx, Bz, boundPar);
143 
144  BOOST_CHECK(!converted.covariance.has_value());
145  BOOST_CHECK(converted.surface);
146 
147  // input is already on perigee, should not be modified
148  BOOST_CHECK_EQUAL(par.template head<2>(),
149  converted.values.template head<2>());
150 
151  // convert back for roundtrip test
152 
153  BoundTrackParameters roundtripPar =
154  EDM4hepUtil::detail::convertTrackParametersFromEdm4hep(Bz, converted);
155 
156  BOOST_CHECK(roundtripPar.parameters().isApprox(boundPar.parameters()));
157  BOOST_CHECK(!roundtripPar.covariance().has_value());
158 }
159 
160 BOOST_AUTO_TEST_CASE(ConvertTrackParametersToEdm4hepWithOutPerigeeNoCov) {
161  auto refSurface = Surface::makeShared<PlaneSurface>(
162  Vector3{50, 30, 20}, Vector3{1, 1, 0.3}.normalized());
163 
164  BoundVector par;
165  par << 1_mm, 5_mm, M_PI / 4., M_PI_2, -1 / 1_GeV, 5_ns;
166 
167  BoundTrackParameters boundPar{refSurface, par, std::nullopt,
168  ParticleHypothesis::pion()};
169 
170  double Bz = 2_T;
171 
172  Acts::GeometryContext gctx;
173 
174  EDM4hepUtil::detail::Parameters converted =
175  EDM4hepUtil::detail::convertTrackParametersToEdm4hep(gctx, Bz, boundPar);
176 
177  BOOST_CHECK(!converted.covariance.has_value());
178  BOOST_CHECK(converted.surface);
179 
180  // input is not a perigee, so new params should be at 0, 0 on ad-hoc perigee
181  BOOST_CHECK_EQUAL(converted.values.template head<2>(), (Vector2{0, 0}));
182  CHECK_CLOSE_ABS(converted.values[2], par[2], 1e-6);
183 
184  // convert back for roundtrip test
185  BoundTrackParameters roundtripPar =
186  EDM4hepUtil::detail::convertTrackParametersFromEdm4hep(Bz, converted);
187 
188  BOOST_CHECK_EQUAL(roundtripPar.parameters().template head<2>(),
189  (Vector2{0, 0}));
190  BOOST_CHECK(roundtripPar.parameters().tail<4>().isApprox(par.tail<4>()));
191  BOOST_CHECK(!roundtripPar.covariance().has_value());
192 }
193 
194 BOOST_AUTO_TEST_CASE(CovariancePacking) {
195  BoundMatrix m;
196  // clang-format off
197  m << 1, 2, 3, 4, 5, 6,
198  2, 2, 3, 4, 5, 6,
199  3, 3, 3, 4, 5, 6,
200  4, 4, 4, 4, 5, 6,
201  5, 5, 5, 5, 5, 6,
202  6, 6, 6, 6, 6, 6;
203  // clang-format on
204 
205  std::array<float, 21> values{};
206  EDM4hepUtil::detail::packCovariance(m, values.data());
207 
208  BoundMatrix m2;
209  m2.setZero();
210  EDM4hepUtil::detail::unpackCovariance(values.data(), m2);
211 
212  CHECK_CLOSE_ABS(m, m2, 1e-9);
213 }
214 
215 BOOST_AUTO_TEST_CASE(RoundTripTests) {
216  auto trackContainer = std::make_shared<Acts::VectorTrackContainer>();
217  auto trackStateContainer = std::make_shared<Acts::VectorMultiTrajectory>();
218  TrackContainer tracks(trackContainer, trackStateContainer);
219 
220  using mutable_proxy_t = decltype(tracks)::TrackProxy;
221  using const_proxy_t = decltype(tracks)::ConstTrackProxy;
222 
223  std::mt19937 rng{42};
224  std::normal_distribution<double> gauss(0., 1.);
225  std::uniform_real_distribution<double> f(-1, 1);
226  std::uniform_real_distribution<double> r(0, 1);
227  std::uniform_int_distribution<unsigned int> nTracks(2, 20);
228  std::uniform_int_distribution<unsigned int> nTs(1, 20);
229  std::uniform_real_distribution<double> phiDist(-M_PI, M_PI);
230  std::uniform_real_distribution<double> etaDist(-4, 4);
231  std::uniform_real_distribution<double> ptDist(1_MeV, 10_GeV);
232  std::uniform_real_distribution<double> qDist(0., 1.);
233 
234  auto genParams = [&]() -> std::pair<BoundVector, BoundMatrix> {
235  double d0 = 20_um * gauss(rng);
236  double z0 = 20_mm * gauss(rng);
237  double phi = phiDist(rng);
238  double eta = etaDist(rng);
239  double theta = 2 * atan(exp(-eta));
240  double pt = ptDist(rng);
241  double p = pt / sin(theta);
242  double charge = qDist(rng) > 0.5 ? 1. : -1.;
243  double qop = charge / p;
244  double t = 5_ns * gauss(rng);
245 
246  BoundVector par;
247  par << d0, z0, phi, theta, qop, t;
248  BoundMatrix cov;
249  cov = BoundMatrix::Identity();
250  cov.diagonal() << 20_um * 20_um, 20_mm * 20_mm, 0.1, 0.1, 1_GeV, 25_ns;
251  return {par, cov};
252  };
253 
254  size_t numT = nTracks(rng);
255  for (size_t t = 0; t < numT; t++) {
256  auto track = tracks.getTrack(tracks.addTrack());
257  {
258  auto [par, cov] = genParams();
259  track.parameters() = par;
260  track.covariance() = cov;
261  }
262  track.setReferenceSurface(
263  Acts::Surface::makeShared<PerigeeSurface>(Vector3{0, 0, 0}));
264 
265  size_t numTs = nTs(rng);
266  for (size_t i = 0; i < numTs; i++) {
267  auto ts = track.appendTrackState(TrackStatePropMask::Smoothed);
268  double crit = r(rng);
269  if (crit < 0.1) {
270  ts.typeFlags().set(TrackStateFlag::HoleFlag);
271  continue;
272  } else if (crit < 0.2) {
273  ts.typeFlags().set(TrackStateFlag::OutlierFlag);
274  continue;
275  } else if (crit < 0.3) {
276  ts.typeFlags().set(TrackStateFlag::SharedHitFlag);
277  } else if (crit < 0.4) {
278  ts.typeFlags().set(TrackStateFlag::MaterialFlag);
279  continue;
280  }
281 
282  ts.typeFlags().set(TrackStateFlag::MeasurementFlag);
283 
284  auto [par, cov] = genParams();
285  ts.smoothed() = par;
286  ts.smoothedCovariance() = cov;
287  Vector3 pos;
288  pos << 1000 * f(rng), 1000 * f(rng), 3000 * f(rng);
289  ts.setReferenceSurface(Acts::Surface::makeShared<PerigeeSurface>(pos));
290  }
291 
292  calculateTrackQuantities(track);
293  }
294 
295  edm4hep::TrackCollection edm4hepTracks;
296 
297  Acts::GeometryContext gctx;
298 
299  double Bz = 3_T;
300 
301  auto logger = getDefaultLogger("EDM4hep", Logging::INFO);
302 
303  for (const_proxy_t track : tracks) {
304  auto to = edm4hepTracks.create();
305  EDM4hepUtil::writeTrack(gctx, track, to, Bz, *logger);
306  }
307 
308  BOOST_CHECK_EQUAL(edm4hepTracks.size(), tracks.size());
309 
310  auto tIt = tracks.begin();
311  for (auto edm4hepTrack : edm4hepTracks) {
312  auto track = *tIt;
313  BOOST_CHECK_EQUAL(track.nMeasurements(),
314  edm4hepTrack.trackStates_size() - 1);
315 
316  ++tIt;
317  }
318 
319  const edm4hep::TrackCollection& edm4hepTracksConst = edm4hepTracks;
320 
321  TrackContainer readTracks(std::make_shared<Acts::VectorTrackContainer>(),
322  std::make_shared<Acts::VectorMultiTrajectory>());
323 
324  for (const auto edm4hepTrack : edm4hepTracksConst) {
325  EDM4hepUtil::readTrack(
326  edm4hepTrack, readTracks.getTrack(readTracks.addTrack()), Bz, *logger);
327  }
328 
329  BOOST_CHECK_EQUAL(tracks.size(), readTracks.size());
330  size_t t = 0;
331 
332  auto origTrackIt = tracks.begin();
333  auto readTrackIt = readTracks.begin();
334  while (origTrackIt != tracks.end() && readTrackIt != readTracks.end()) {
335  BOOST_TEST_INFO_SCOPE("Track #" << t);
336  auto orig = *origTrackIt;
337  auto read = *readTrackIt;
338 
339  CHECK_CLOSE_OR_SMALL(orig.parameters(), read.parameters(), 1e-5, 1e-8);
340  CHECK_CLOSE_OR_SMALL(orig.covariance(), read.covariance(), 1e-5, 1e-8);
341  BOOST_CHECK_EQUAL(orig.referenceSurface().center(gctx),
342  read.referenceSurface().center(gctx));
343 
344  auto origTsIt = orig.trackStatesReversed().begin();
345  auto readTsIt = read.trackStatesReversed().begin();
346 
347  size_t tsi = 0;
348  while (origTsIt != orig.trackStatesReversed().end() &&
349  readTsIt != read.trackStatesReversed().end()) {
350  BOOST_TEST_INFO_SCOPE("TS: #" << tsi);
351  auto nextMeas = std::find_if(
352  origTsIt, orig.trackStatesReversed().end(), [](const auto& ts) {
353  return ts.typeFlags().test(TrackStateFlag::MeasurementFlag);
354  });
355  BOOST_CHECK(nextMeas != orig.trackStatesReversed().end());
356  origTsIt = nextMeas;
357  auto origTs = *origTsIt;
358  auto readTs = *readTsIt;
359 
360  BOOST_TEST_INFO_SCOPE(
361  "orig parameters: " << origTs.parameters().transpose());
362  BOOST_TEST_INFO_SCOPE(
363  "read parameters: " << readTs.parameters().transpose());
364  CHECK_CLOSE_OR_SMALL(origTs.smoothedCovariance(),
365  readTs.smoothedCovariance(), 1e-5, 1e-6);
366  Vector3 newCenter = readTs.referenceSurface().center(
367  gctx); // new center is a perigee, but should be on the other
368  // surface
369  BOOST_CHECK(origTs.referenceSurface().isOnSurface(gctx, newCenter,
370  Vector3::Zero()));
371 
372  // global hit positions should be the same
373  Vector3 readGlobal = readTs.referenceSurface().localToGlobal(
374  gctx, readTs.parameters().template head<2>(), Vector3::Zero());
375  Vector3 origGlobal = origTs.referenceSurface().localToGlobal(
376  gctx, origTs.parameters().template head<2>(), Vector3::Zero());
377  CHECK_CLOSE_ABS(readGlobal, origGlobal, 1e-3);
378  ++origTsIt;
379  ++readTsIt;
380  tsi++;
381  }
382  ++origTrackIt;
383  ++readTrackIt;
384 
385  t++;
386  }
387 }
388 
389 BOOST_AUTO_TEST_SUITE_END()