12 #include "GaudiKernel/IInterface.h"
29 #include "Acts/Definitions/Units.hpp"
30 #include "Acts/EventData/TrackParameters.hpp"
31 #include "Acts/EventData/VectorTrackContainer.hpp"
32 #include "Acts/EventData/TransformationHelpers.hpp"
33 #include "Acts/Geometry/TrackingGeometry.hpp"
34 #include "Acts/Propagator/detail/JacobianEngine.hpp"
35 #include "Acts/Surfaces/PerigeeSurface.hpp"
36 #include "Acts/Surfaces/Surface.hpp"
38 #include "Acts/EventData/TrackStatePropMask.hpp"
39 #include "Acts/EventData/SourceLink.hpp"
52 #pragma GCC diagnostic push
53 #pragma GCC diagnostic ignored "-Wunused-function"
54 static void ActsMeasurementCheck(
const Acts::GeometryContext &gctx,
56 const Acts::Surface &surface,
57 const Acts::BoundVector &loc);
58 #pragma GCC diagnostic pop
60 static void ActsTrackParameterCheck(
61 const Acts::BoundTrackParameters &actsParameter,
62 const Acts::GeometryContext &gctx,
const Acts::BoundSquareMatrix &covpc,
63 const Acts::BoundVector &targetPars,
const Acts::BoundSquareMatrix &targetCov,
69 const std::string &
type,
const std::string &
name,
const IInterface *
parent)
74 if (m_extractMuonSurfaces){
77 if (!m_trackingGeometryTool.empty()) {
78 ATH_CHECK(m_trackingGeometryTool.retrieve());
79 m_trackingGeometry = m_trackingGeometryTool->trackingGeometry();
81 m_trackingGeometry->visitSurfaces([&](
const Acts::Surface *surface) {
86 surface->associatedDetectorElement());
92 actsElement->upstreamDetectorElement()) !=
nullptr);
97 m_actsSurfaceMap.insert({actsElement->identify(), surface});
100 <<
" has two ACTS surfaces: "
101 <<
it->second->geometryId() <<
" and "
102 << surface->geometryId());
107 if (m_extractMuonSurfaces){
110 unsigned int mapSize = m_actsSurfaceMap.size();
113 for (
auto [
id, surf] : reSurfaces){
114 const Acts::Surface*
tmp = surf.get();
115 m_actsSurfaceMap.insert(std::pair<Identifier, const Acts::Surface*>(
id,
tmp));
118 ATH_MSG_VERBOSE(
"After adding muon surfaces, the map has grown from "<<mapSize<<
" to "<<m_actsSurfaceMap.size());
120 return StatusCode::SUCCESS;
124 const Acts::Surface &actsSurface)
const {
127 actsSurface.associatedDetectorElement());
131 throw std::domain_error(
"No ATLAS surface corresponding to to the Acts one");
138 auto it = m_actsSurfaceMap.find(atlasID);
139 if (
it != m_actsSurfaceMap.end()) {
142 ATH_MSG_ERROR(
"No Acts surface corresponding to this ATLAS surface:");
144 throw std::domain_error(
"No Acts surface corresponding to the ATLAS one");
148 [[maybe_unused]]
const Acts::GeometryContext& gctx,
151 return Acts::SourceLink(&measurement);
155 std::vector<Acts::SourceLink>
159 auto outliers =
track.outliersOnTrack();
161 std::vector<Acts::SourceLink> sourceLinks;
162 sourceLinks.reserve(
hits->size() + outliers->size());
165 sourceLinks.push_back(trkMeasurementToSourceLink(gctx, **
it));
167 for (
auto it = outliers->begin();
it != outliers->end(); ++
it) {
168 sourceLinks.push_back(trkMeasurementToSourceLink(gctx, **
it));
174 const Acts::BoundTrackParameters
178 using namespace Acts::UnitLiterals;
179 std::shared_ptr<const Acts::Surface> actsSurface;
189 ATH_MSG_ERROR(
"Could not find ACTS detector surface for this TrackParameter:");
197 "trkTrackParametersToActsParameters:: No associated surface found (owner: "<<atlasParameter.
associatedSurface().
owner()<<
198 "). Creating a free surface. Trk parameters:");
203 actsSurface = Acts::Surface::makeShared<const Acts::PlaneSurface>(
207 actsSurface = Acts::Surface::makeShared<const Acts::PerigeeSurface>(
212 ATH_MSG_WARNING(
"No surface type found for this Trk::Surface. Creating a perigee surface.");
213 actsSurface = Acts::Surface::makeShared<const Acts::PerigeeSurface>(
219 auto atlasParam = atlasParameter.parameters();
220 if (actsSurface->bounds().type() ==
221 Acts::SurfaceBounds::BoundsType::eAnnulus) {
224 auto position = atlasParameter.
position();
226 actsSurface->globalToLocal(gctx, position, atlasParameter.
momentum());
230 atlasParameter.
charge() / (atlasParameter.
momentum().mag() * 1_MeV),
234 "Unable to convert annulus surface - globalToLocal failed");
239 atlasParameter.
charge() / (atlasParameter.
momentum().mag() * 1_MeV), 0.;
242 Acts::BoundSquareMatrix
cov = Acts::BoundSquareMatrix::Identity();
243 if (atlasParameter.covariance()) {
244 cov.topLeftCorner(5, 5) = *atlasParameter.covariance();
249 for (
int i = 0;
i <
cov.rows();
i++) {
252 for (
int i = 0;
i <
cov.cols();
i++) {
259 Acts::PdgParticle absPdg = Acts::makeAbsolutePdgParticle(
260 static_cast<Acts::PdgParticle
>(
261 m_pdgToParticleHypothesis.convert(hypothesis, atlasParameter.
charge())));
263 absPdg,
mass, Acts::AnyCharge{std::abs(
static_cast<float>(atlasParameter.
charge()))}};
265 return Acts::BoundTrackParameters(actsSurface,
params,
266 cov, actsHypothesis);
269 std::unique_ptr<Trk::TrackParameters>
271 const Acts::BoundTrackParameters &actsParameter,
272 const Acts::GeometryContext &gctx)
const {
274 using namespace Acts::UnitLiterals;
276 if (actsParameter.covariance()) {
278 newcov(actsParameter.covariance()->topLeftCorner<5, 5>());
280 for (
int i = 0;
i < newcov.rows();
i++) {
281 newcov(
i, 4) = newcov(
i, 4) * 1_MeV;
283 for (
int i = 0;
i < newcov.cols();
i++) {
284 newcov(4,
i) = newcov(4,
i) * 1_MeV;
286 cov = std::optional<AmgSymMatrix(5)>(newcov);
289 const Acts::Surface &actsSurface = actsParameter.referenceSurface();
291 switch (actsSurface.type()) {
296 actsSurfaceToTrkSurface(actsSurface));
297 return std::make_unique<Trk::AtaCone>(
298 actsParameter.get<Acts::eBoundLoc0>(),
299 actsParameter.get<Acts::eBoundLoc1>(),
300 actsParameter.get<Acts::eBoundPhi>(),
301 actsParameter.get<Acts::eBoundTheta>(),
302 actsParameter.get<Acts::eBoundQOverP>() * 1_MeV, coneSurface,
cov);
309 actsSurfaceToTrkSurface(actsSurface));
310 return std::make_unique<Trk::AtaCylinder>(
311 actsParameter.get<Acts::eBoundLoc0>(),
312 actsParameter.get<Acts::eBoundLoc1>(),
313 actsParameter.get<Acts::eBoundPhi>(),
314 actsParameter.get<Acts::eBoundTheta>(),
315 actsParameter.get<Acts::eBoundQOverP>() * 1_MeV, cylSurface,
cov);
320 &actsSurfaceToTrkSurface(actsSurface));
321 discSurface !=
nullptr) {
322 return std::make_unique<Trk::AtaDisc>(
323 actsParameter.get<Acts::eBoundLoc0>(),
324 actsParameter.get<Acts::eBoundLoc1>(),
325 actsParameter.get<Acts::eBoundPhi>(),
326 actsParameter.get<Acts::eBoundTheta>(),
327 actsParameter.get<Acts::eBoundQOverP>() * 1_MeV, *discSurface,
cov);
328 }
else if (
const auto *planeSurface =
330 &actsSurfaceToTrkSurface(actsSurface));
331 planeSurface !=
nullptr) {
333 auto helperSurface = Acts::Surface::makeShared<Acts::PlaneSurface>(
334 planeSurface->transform());
336 auto covpc = actsParameter.covariance().value();
338 Acts::FreeVector freePars =
339 Acts::transformBoundToFreeParameters(
340 actsSurface, gctx, actsParameter.parameters());
342 Acts::BoundVector targetPars =
343 Acts::transformFreeToBoundParameters(freePars,
344 *helperSurface, gctx)
347 Acts::FreeMatrix freeTransportJacobian = Acts::FreeMatrix::Identity();
350 freeToPathDerivatives.head<3>() = freePars.segment<3>(Acts::eFreeDir0);
352 auto boundToFreeJacobian = actsSurface.boundToFreeJacobian(
353 gctx, freePars.segment<3>(Acts::eFreePos0),
354 freePars.segment<3>(Acts::eFreeDir0));
356 Acts::BoundMatrix boundToBoundJac = Acts::detail::boundToBoundTransportJacobian(
357 gctx, freePars, boundToFreeJacobian, freeTransportJacobian,
358 freeToPathDerivatives, *helperSurface);
360 Acts::BoundMatrix targetCov =
361 boundToBoundJac * covpc * boundToBoundJac.transpose();
363 auto pars = std::make_unique<Trk::AtaPlane>(
364 targetPars[Acts::eBoundLoc0], targetPars[Acts::eBoundLoc1],
365 targetPars[Acts::eBoundPhi], targetPars[Acts::eBoundTheta],
366 targetPars[Acts::eBoundQOverP] * 1_MeV, *planeSurface,
367 targetCov.topLeftCorner<5, 5>());
369 if (m_visualDebugOutput) {
370 ActsTrackParameterCheck(actsParameter, gctx, covpc, targetPars,
371 targetCov, planeSurface);
377 throw std::runtime_error{
378 "Acts::DiscSurface is not associated with ATLAS disc or plane "
383 case Acts::Surface::SurfaceType::Perigee: {
385 return std::make_unique<Trk::Perigee>(
386 actsParameter.get<Acts::eBoundLoc0>(),
387 actsParameter.get<Acts::eBoundLoc1>(),
388 actsParameter.get<Acts::eBoundPhi>(),
389 actsParameter.get<Acts::eBoundTheta>(),
390 actsParameter.get<Acts::eBoundQOverP>() * 1_MeV, perSurface,
cov);
397 actsSurfaceToTrkSurface(actsSurface));
398 return std::make_unique<Trk::AtaPlane>(
399 actsParameter.get<Acts::eBoundLoc0>(),
400 actsParameter.get<Acts::eBoundLoc1>(),
401 actsParameter.get<Acts::eBoundPhi>(),
402 actsParameter.get<Acts::eBoundTheta>(),
403 actsParameter.get<Acts::eBoundQOverP>() * 1_MeV, plaSurface,
cov);
409 actsSurfaceToTrkSurface(actsSurface));
410 return std::make_unique<Trk::AtaStraightLine>(
411 actsParameter.get<Acts::eBoundLoc0>(),
412 actsParameter.get<Acts::eBoundLoc1>(),
413 actsParameter.get<Acts::eBoundPhi>(),
414 actsParameter.get<Acts::eBoundTheta>(),
415 actsParameter.get<Acts::eBoundQOverP>() * 1_MeV, lineSurface,
cov);
419 return std::make_unique<Trk::CurvilinearParameters>(
420 actsParameter.position(gctx), actsParameter.get<Acts::eBoundPhi>(),
421 actsParameter.get<Acts::eBoundTheta>(),
422 actsParameter.get<Acts::eBoundQOverP>() * 1_MeV,
cov);
425 case Acts::Surface::SurfaceType::Other: {
430 throw std::domain_error(
"Surface type not found");
436 const Acts::GeometryContext & gctx)
const {
438 << trackColl.size() <<
" tracks.");
439 unsigned int trkCount = 0;
440 std::vector<Identifier> failedIds;
441 for (
auto trk : trackColl) {
444 trk->trackStateOnSurfaces();
450 auto actsTrack = tc.getTrack(tc.addTrack());
451 auto& trackStateContainer = tc.trackStateContainer();
454 <<
" track states on surfaces.");
456 actsTrack.chi2() = trk->fitQuality()->chiSquared();
457 actsTrack.nDoF() = trk->fitQuality()->numberDoF();
461 bool first_tsos =
true;
462 int measurementsCount = 0;
463 for (
auto tsos : *trackStates) {
467 if (tsos->measurementOnTrack()) {
468 mask |= Acts::TrackStatePropMask::Calibrated;
470 if (tsos->trackParameters()) {
471 mask |= Acts::TrackStatePropMask::Smoothed;
475 auto index = Acts::MultiTrajectoryTraits::kInvalid;
477 index = actsTrack.tipIndex();
479 auto actsTSOS = trackStateContainer.getTrackState(
480 trackStateContainer.addTrackState(
mask,
index));
482 <<
" TSOS index within trajectory: "
483 << actsTSOS.index());
484 actsTrack.tipIndex() = actsTSOS.index();
486 if (tsos->trackParameters()) {
492 trkTrackParametersToActsParameters(*(tsos->trackParameters()), gctx);
495 if (!actsTrackParameterPositionCheck(
parameters, *(tsos->trackParameters()), gctx)){
496 failedIds.push_back(tsos->trackParameters()->associatedSurface().associatedDetectorElementIdentifier());
502 actsTrack.parameters() =
parameters.parameters();
503 actsTrack.covariance() = *
parameters.covariance();
504 actsTrack.setReferenceSurface(
505 parameters.referenceSurface().getSharedPtr());
508 actsTSOS.setReferenceSurface(
parameters.referenceSurface().getSharedPtr());
510 actsTSOS.smoothed() =
parameters.parameters();
511 actsTSOS.smoothedCovariance() = *
parameters.covariance();
514 if (!(actsTSOS.hasSmoothed() && actsTSOS.hasReferenceSurface())) {
516 << actsTSOS.hasSmoothed()
517 <<
"] or reference surface ["
518 << actsTSOS.hasReferenceSurface() <<
"].");
520 ATH_MSG_VERBOSE(
"TrackState has smoothed state and reference surface.");
524 ATH_MSG_ERROR(
"Unable to convert TrackParameter with exception ["<<
e.what()<<
"]. Will be missing from ACTS track."
525 <<(*tsos->trackParameters()));
528 if (tsos->measurementOnTrack()) {
529 auto &measurement = *(tsos->measurementOnTrack());
537 int dim = measurement.localParameters().dimension();
538 actsTSOS.allocateCalibrated(
dim);
540 actsTSOS.calibrated<1>() = measurement.localParameters();
541 actsTSOS.calibratedCovariance<1>() = measurement.localCovariance();
542 }
else if (
dim == 2) {
543 actsTSOS.calibrated<2>() = measurement.localParameters();
544 actsTSOS.calibratedCovariance<2>() = measurement.localCovariance();
546 throw std::domain_error(
"Cannot handle measurement dim>2");
548 actsTSOS.setUncalibratedSourceLink(Acts::SourceLink(
static_cast<ActsTrk::ATLASSourceLink>(tsos->measurementOnTrack())));
552 actsTrack.nMeasurements() = measurementsCount;
554 <<
" track states on surfaces.");
556 ATH_MSG_VERBOSE(
"Finished converting " << trackColl.size() <<
" tracks.");
558 if (!failedIds.empty()){
559 ATH_MSG_WARNING(
"Failed to convert "<<failedIds.size()<<
" track parameters.");
560 for (
auto id : failedIds){
561 ATH_MSG_WARNING(
"-> Failed for Identifier "<<m_idHelperSvc->toString(
id));
564 ATH_MSG_VERBOSE(
"ACTS Track container has " << tc.size() <<
" tracks.");
570 const Acts::GeometryContext &gctx)
const {
578 if (std::fabs(actsPos.x() - trkparameters.
position().x()) >
580 std::fabs(actsPos.y() - trkparameters.
position().y()) >
582 std::fabs(actsPos.z() - trkparameters.
position().z()) >
585 << actsPos <<
" vs Trk \n"
598 #pragma GCC diagnostic push
599 #pragma GCC diagnostic ignored "-Wunused-function"
600 static void ActsTrk::ActsMeasurementCheck(
602 const Acts::Surface &surface,
const Acts::BoundVector &loc) {
609 if (bounds ==
nullptr) {
610 throw std::runtime_error{
"Annulus but not XY"};
619 if (
auto res = surface.globalToLocal(gctx,
global, Acts::Vector3{});
623 throw std::runtime_error{
"Global position not on target surface"};
629 Acts::Surface::makeShared<Acts::PlaneSurface>(surf.
transform());
630 Acts::BoundVector locxypar;
631 locxypar.head<2>() = locxy;
633 locxypar[3] = M_PI_2;
636 Acts::FreeVector globalxypar = Acts::transformBoundToFreeParameters(
637 *planeSurface, gctx, locxypar);
638 auto boundToFree = planeSurface->boundToFreeJacobian(
639 gctx, globalxypar.segment<3>(Acts::eFreePos0),
640 globalxypar.segment<3>(Acts::eFreeDir0));
641 Acts::ActsSquareMatrix<2> xyToXyzJac = boundToFree.topLeftCorner<2, 2>();
643 Acts::BoundVector locpcpar;
644 locpcpar.head<2>() = locpc;
646 locpcpar[3] = M_PI_2;
649 Acts::FreeVector globalpcpar = Acts::transformBoundToFreeParameters(
650 surface, gctx, locpcpar);
652 boundToFree = surface.boundToFreeJacobian(
653 gctx, globalpcpar.segment<3>(Acts::eFreePos0),
654 globalpcpar.segment<3>(Acts::eFreeDir0));
655 Acts::ActsSquareMatrix<2> pcToXyzJac = boundToFree.topLeftCorner<2, 2>();
656 Acts::ActsSquareMatrix<2> xyzToPcJac = pcToXyzJac.inverse();
659 Acts::ActsMatrix<2, 2> covpc = covxy;
660 covpc = xyToXyzJac * covpc * xyToXyzJac.transpose();
661 covpc = xyzToPcJac * covpc * xyzToPcJac.transpose();
663 std::mt19937
gen{42 + surface.geometryId().value()};
664 std::normal_distribution<double> normal{0, 1};
665 std::uniform_real_distribution<double>
uniform{-1, 1};
667 Acts::ActsMatrix<2, 2> lltxy = covxy.llt().matrixL();
668 Acts::ActsMatrix<2, 2> lltpc = covpc.llt().matrixL();
670 for (
size_t i = 0;
i < 1
e4;
i++) {
671 std::cout <<
"ANNULUS COV: ";
672 std::cout << surface.geometryId();
680 std::cout <<
"," << xy.x() <<
"," << xy.y();
681 std::cout <<
"," <<
xyz.x() <<
"," <<
xyz.y() <<
"," <<
xyz.z();
691 std::cout <<
"," << rt.x() <<
"," << rt.y();
692 std::cout <<
"," <<
xyz.x() <<
"," <<
xyz.y() <<
"," <<
xyz.z();
695 std::cout << std::endl;
698 #pragma GCC diagnostic pop
700 void ActsTrk::ActsTrackParameterCheck(
701 const Acts::BoundTrackParameters &actsParameter,
702 const Acts::GeometryContext &gctx,
const Acts::BoundSquareMatrix &covpc,
703 const Acts::BoundVector &targetPars,
const Acts::BoundSquareMatrix &targetCov,
706 std::cout <<
"ANNULUS PAR COV: ";
707 std::cout << actsParameter.referenceSurface().geometryId();
708 for (
unsigned int i = 0;
i < 5;
i++) {
709 for (
unsigned int j = 0; j < 5; j++) {
710 std::cout <<
"," << covpc(
i, j);
713 for (
unsigned int i = 0;
i < 5;
i++) {
714 for (
unsigned int j = 0; j < 5; j++) {
715 std::cout <<
"," << targetCov(
i, j);
718 std::cout << std::endl;
720 std::mt19937
gen{4242 +
721 actsParameter.referenceSurface().geometryId().value()};
722 std::normal_distribution<double> normal{0, 1};
724 Acts::ActsMatrix<2, 2> lltxy =
725 targetCov.topLeftCorner<2, 2>().llt().matrixL();
726 Acts::ActsMatrix<2, 2> lltpc = covpc.topLeftCorner<2, 2>().llt().matrixL();
728 for (
size_t i = 0;
i < 1
e4;
i++) {
729 std::cout <<
"ANNULUS PAR: ";
730 std::cout << actsParameter.referenceSurface().geometryId();
732 Acts::ActsVector<2> rnd;
733 rnd << normal(
gen), normal(
gen);
737 Acts::ActsVector<2> xy =
738 lltxy.topLeftCorner<2, 2>() * rnd + targetPars.head<2>();
742 for (
unsigned int i = 0;
i < 2;
i++) {
743 std::cout <<
"," << xy[
i];
745 std::cout <<
"," <<
xyz.x() <<
"," <<
xyz.y() <<
"," <<
xyz.z();
749 Acts::ActsVector<2> rt = lltpc.topLeftCorner<2, 2>() * rnd +
750 actsParameter.parameters().head<2>();
752 gctx, Acts::Vector2{rt.head<2>()}, Acts::Vector3{});
754 for (
unsigned int i = 0;
i < 2;
i++) {
755 std::cout <<
"," << rt[
i];
757 std::cout <<
"," <<
xyz.x() <<
"," <<
xyz.y() <<
"," <<
xyz.z();
760 std::cout << std::endl;