Back to home page

EIC code displayed by LXR
 
 

    

Warning, file /EICrecon/src/algorithms/tracking/TrackPropagation.cc was not indexed or was modified since last indexation (in which case cross-reference links may be missing, inaccurate or erroneous).

0001 // SPDX-License-Identifier: LGPL-3.0-or-later
0002 // Copyright (C) 2022, 2023 Wenqing Fan, Barak Schmookler, Whitney Armstrong, Sylvester Joosten, Dmitry Romanov, Christopher Dilks, Wouter Deconinck
0003 
0004 #include <Acts/Definitions/Algebra.hpp>
0005 #include <Acts/Definitions/Direction.hpp>
0006 #include <Acts/Definitions/TrackParametrization.hpp>
0007 #include <Acts/Definitions/Units.hpp>
0008 #include <Acts/EventData/GenericBoundTrackParameters.hpp>
0009 #include <Acts/EventData/MultiTrajectoryHelpers.hpp>
0010 #include <Acts/EventData/ParticleHypothesis.hpp>
0011 #include <Acts/Geometry/GeometryIdentifier.hpp>
0012 #include <Acts/Geometry/TrackingGeometry.hpp>
0013 #include <Acts/MagneticField/MagneticFieldProvider.hpp>
0014 #include <Acts/Propagator/EigenStepper.hpp>
0015 #include <Acts/Propagator/Propagator.hpp>
0016 #include <Acts/Surfaces/CylinderBounds.hpp>
0017 #include <Acts/Surfaces/CylinderSurface.hpp>
0018 #include <Acts/Surfaces/DiscSurface.hpp>
0019 #include <Acts/Surfaces/RadialBounds.hpp>
0020 #include <Acts/Utilities/Logger.hpp>
0021 #include <ActsExamples/EventData/Trajectories.hpp>
0022 #include <DD4hep/Handle.h>
0023 #include <Eigen/Core>
0024 #include <Eigen/Geometry>
0025 #include <Evaluator/DD4hepUnits.h>
0026 #include <algorithm>
0027 #include <boost/container/vector.hpp>
0028 #include <cmath>
0029 #include <cstdint>
0030 #include <edm4hep/Vector3f.h>
0031 #include <edm4hep/utils/vector_utils.h>
0032 #include <fmt/core.h>
0033 #include <functional>
0034 #include <iterator>
0035 #include <map>
0036 #include <optional>
0037 #include <stdexcept>
0038 #include <string>
0039 #include <tuple>
0040 #include <typeinfo>
0041 #include <variant>
0042 
0043 #include "algorithms/tracking/ActsGeometryProvider.h"
0044 #include "algorithms/tracking/TrackPropagation.h"
0045 #include "algorithms/tracking/TrackPropagationConfig.h"
0046 #include "extensions/spdlog/SpdlogToActs.h"
0047 
0048 namespace eicrecon {
0049 
0050 template <typename... L> struct multilambda : L... {
0051   using L::operator()...;
0052   constexpr multilambda(L... lambda) : L(std::move(lambda))... {}
0053 };
0054 
0055 void TrackPropagation::init(const dd4hep::Detector* detector,
0056                             std::shared_ptr<const ActsGeometryProvider> geo_svc,
0057                             std::shared_ptr<spdlog::logger> logger) {
0058   m_geoSvc = geo_svc;
0059   m_log    = logger;
0060 
0061   std::map<uint32_t, std::size_t> system_id_layers;
0062 
0063   multilambda _toDouble = {
0064       [](const std::string& v) { return dd4hep::_toDouble(v); },
0065       [](const double& v) { return v; },
0066   };
0067 
0068   auto _toActsSurface =
0069       [&_toDouble, &detector, &system_id_layers](
0070           const std::variant<CylinderSurfaceConfig, DiscSurfaceConfig> surface_variant)
0071       -> std::shared_ptr<Acts::Surface> {
0072     if (std::holds_alternative<CylinderSurfaceConfig>(surface_variant)) {
0073       CylinderSurfaceConfig surface = std::get<CylinderSurfaceConfig>(surface_variant);
0074       const double rmin =
0075           std::visit(_toDouble, surface.rmin) / dd4hep::mm * Acts::UnitConstants::mm;
0076       const double zmin =
0077           std::visit(_toDouble, surface.zmin) / dd4hep::mm * Acts::UnitConstants::mm;
0078       const double zmax =
0079           std::visit(_toDouble, surface.zmax) / dd4hep::mm * Acts::UnitConstants::mm;
0080       const uint32_t system_id = detector->constant<uint32_t>(surface.id);
0081       auto bounds              = std::make_shared<Acts::CylinderBounds>(rmin, (zmax - zmin) / 2);
0082       auto t                   = Acts::Translation3(Acts::Vector3(0, 0, (zmax + zmin) / 2));
0083       auto tf                  = Acts::Transform3(t);
0084       auto acts_surface        = Acts::Surface::makeShared<Acts::CylinderSurface>(tf, bounds);
0085 #if Acts_VERSION_MAJOR >= 40
0086       acts_surface->assignGeometryId(
0087           Acts::GeometryIdentifier().withExtra(system_id).withLayer(++system_id_layers[system_id]));
0088 #else
0089       acts_surface->assignGeometryId(
0090           Acts::GeometryIdentifier().setExtra(system_id).setLayer(++system_id_layers[system_id]));
0091 #endif
0092       return acts_surface;
0093     }
0094     if (std::holds_alternative<DiscSurfaceConfig>(surface_variant)) {
0095       DiscSurfaceConfig surface = std::get<DiscSurfaceConfig>(surface_variant);
0096       const double zmin =
0097           std::visit(_toDouble, surface.zmin) / dd4hep::mm * Acts::UnitConstants::mm;
0098       const double rmin =
0099           std::visit(_toDouble, surface.rmin) / dd4hep::mm * Acts::UnitConstants::mm;
0100       const double rmax =
0101           std::visit(_toDouble, surface.rmax) / dd4hep::mm * Acts::UnitConstants::mm;
0102       const uint32_t system_id = detector->constant<uint32_t>(surface.id);
0103       auto bounds              = std::make_shared<Acts::RadialBounds>(rmin, rmax);
0104       auto t                   = Acts::Translation3(Acts::Vector3(0, 0, zmin));
0105       auto tf                  = Acts::Transform3(t);
0106       auto acts_surface        = Acts::Surface::makeShared<Acts::DiscSurface>(tf, bounds);
0107 #if Acts_VERSION_MAJOR >= 40
0108       acts_surface->assignGeometryId(
0109           Acts::GeometryIdentifier().withExtra(system_id).withLayer(++system_id_layers[system_id]));
0110 #else
0111       acts_surface->assignGeometryId(
0112           Acts::GeometryIdentifier().setExtra(system_id).setLayer(++system_id_layers[system_id]));
0113 #endif
0114       return acts_surface;
0115     }
0116     throw std::domain_error("Unknown surface type");
0117   };
0118   m_target_surfaces.resize(m_cfg.target_surfaces.size());
0119   std::transform(m_cfg.target_surfaces.cbegin(), m_cfg.target_surfaces.cend(),
0120                  m_target_surfaces.begin(), _toActsSurface);
0121   m_filter_surfaces.resize(m_cfg.filter_surfaces.size());
0122   std::transform(m_cfg.filter_surfaces.cbegin(), m_cfg.filter_surfaces.cend(),
0123                  m_filter_surfaces.begin(), _toActsSurface);
0124 
0125   m_log->trace("Initialized");
0126 }
0127 
0128 void TrackPropagation::propagateToSurfaceList(
0129     const std::tuple<const edm4eic::TrackCollection&,
0130                      const std::vector<const ActsExamples::Trajectories*>,
0131                      const std::vector<const ActsExamples::ConstTrackContainer*>>
0132         input,
0133     const std::tuple<edm4eic::TrackSegmentCollection*> output) const {
0134   const auto [tracks, acts_trajectories, acts_tracks] = input;
0135   auto [track_segments]                               = output;
0136 
0137   // logging
0138   m_log->trace("Propagate trajectories: --------------------");
0139   m_log->trace("number of tracks: {}", tracks.size());
0140   m_log->trace("number of acts_trajectories: {}", acts_trajectories.size());
0141   m_log->trace("number of acts_tracks: {}", acts_tracks.size());
0142 
0143   // loop over input trajectories
0144   for (std::size_t i = 0; const auto& traj : acts_trajectories) {
0145 
0146     // check if this trajectory can be propagated to any filter surface
0147     bool trajectory_reaches_filter_surface{false};
0148     for (const auto& filter_surface : m_filter_surfaces) {
0149       auto point = propagate(edm4eic::Track{}, traj, filter_surface);
0150       if (point) {
0151         trajectory_reaches_filter_surface = true;
0152         break;
0153       }
0154     }
0155     if (!trajectory_reaches_filter_surface) {
0156       ++i;
0157       continue;
0158     }
0159 
0160     // start a mutable TrackSegment
0161     auto track_segment = track_segments->create();
0162 
0163     // corresponding track
0164     if (tracks.size() == acts_trajectories.size()) {
0165       m_log->trace("track segment connected to track {}", i);
0166       track_segment.setTrack(tracks[i]);
0167       ++i;
0168     }
0169 
0170     // zero measurements of segment length
0171     decltype(edm4eic::TrackSegmentData::length) length            = 0;
0172     decltype(edm4eic::TrackSegmentData::lengthError) length_error = 0;
0173 
0174     // loop over projection-target surfaces
0175     for (const auto& target_surface : m_target_surfaces) {
0176 
0177       // project the trajectory `traj` to this surface
0178       auto point = propagate(edm4eic::Track{}, traj, target_surface);
0179       if (!point) {
0180         m_log->trace("<> Failed to propagate trajectory to this plane");
0181         continue;
0182       }
0183 
0184       // logging
0185       m_log->trace("<> trajectory: x=( {:>10.2f} {:>10.2f} {:>10.2f} )", point->position.x,
0186                    point->position.y, point->position.z);
0187       m_log->trace("               p=( {:>10.2f} {:>10.2f} {:>10.2f} )", point->momentum.x,
0188                    point->momentum.y, point->momentum.z);
0189 
0190       // track point cut
0191       if (!m_cfg.track_point_cut(*point)) {
0192         m_log->trace("                 => REJECTED by trackPointCut");
0193         if (m_cfg.skip_track_on_track_point_cut_failure) {
0194           break;
0195         }
0196         continue;
0197       }
0198 
0199       // update the `TrackSegment` length
0200       // FIXME: `length` and `length_error` are currently not used by any callers, and may not be correctly calculated here
0201       if (track_segment.points_size() > 0) {
0202         auto pos0 = point->position;
0203         auto pos1 = std::prev(track_segment.points_end())->position;
0204         auto dist = edm4hep::utils::magnitude(pos0 - pos1);
0205         length += dist;
0206         m_log->trace("               dist to previous point: {}", dist);
0207       }
0208 
0209       // add the `TrackPoint` to the `TrackSegment`
0210       track_segment.addToPoints(*point);
0211 
0212     } // end `targetSurfaces` loop
0213 
0214     // set final length and length error
0215     track_segment.setLength(length);
0216     track_segment.setLengthError(length_error);
0217 
0218   } // end loop over input trajectories
0219 }
0220 
0221 std::unique_ptr<edm4eic::TrackPoint>
0222 TrackPropagation::propagate(const edm4eic::Track& /* track */,
0223                             const ActsExamples::Trajectories* acts_trajectory,
0224                             const std::shared_ptr<const Acts::Surface>& targetSurf) const {
0225 
0226   // Get the entry index for the single trajectory
0227   // The trajectory entry indices and the multiTrajectory
0228   const auto& mj        = acts_trajectory->multiTrajectory();
0229   const auto& trackTips = acts_trajectory->tips();
0230 
0231   m_log->trace("  Number of elements in trackTips {}", trackTips.size());
0232 
0233   // Skip empty
0234   if (trackTips.empty()) {
0235     m_log->trace("  Empty multiTrajectory.");
0236     return nullptr;
0237   }
0238   const auto& trackTip = trackTips.front();
0239 
0240   // Collect the trajectory summary info
0241   auto trajState      = Acts::MultiTrajectoryHelpers::trajectoryState(mj, trackTip);
0242   int m_nMeasurements = trajState.nMeasurements;
0243   int m_nStates       = trajState.nStates;
0244 
0245   m_log->trace("  Num measurement in trajectory: {}", m_nMeasurements);
0246   m_log->trace("  Num states in trajectory     : {}", m_nStates);
0247 
0248   // Get track state at last measurement surface
0249   // For last measurement surface, filtered and smoothed results are equivalent
0250   auto trackState        = mj.getTrackState(trackTip);
0251   auto initSurface       = trackState.referenceSurface().getSharedPtr();
0252   const auto& initParams = trackState.filtered();
0253   const auto& initCov    = trackState.filteredCovariance();
0254 
0255   Acts::BoundTrackParameters initBoundParams(initSurface, initParams, initCov,
0256                                              Acts::ParticleHypothesis::pion());
0257 
0258   // Get pathlength of last track state with respect to perigee surface
0259   const auto initPathLength = trackState.pathLength();
0260 
0261   m_log->trace("    TrackPropagation. Propagating to surface # {}",
0262                typeid(targetSurf->type()).name());
0263 
0264   std::shared_ptr<const Acts::TrackingGeometry> trackingGeometry   = m_geoSvc->trackingGeometry();
0265   std::shared_ptr<const Acts::MagneticFieldProvider> magneticField = m_geoSvc->getFieldProvider();
0266   using Stepper                                                    = Acts::EigenStepper<>;
0267   using Propagator                                                 = Acts::Propagator<Stepper>;
0268   Stepper stepper(magneticField);
0269   Propagator propagator(stepper);
0270 
0271   ACTS_LOCAL_LOGGER(eicrecon::getSpdlogLogger("PROP", m_log));
0272 
0273 #if Acts_VERSION_MAJOR >= 36
0274   Propagator::template Options<> options(m_geoContext, m_fieldContext);
0275 #else
0276   Acts::PropagatorOptions<> options(m_geoContext, m_fieldContext);
0277 #endif
0278 
0279   auto result = propagator.propagate(initBoundParams, *targetSurf, options);
0280 
0281   // check propagation result
0282   if (!result.ok()) {
0283     m_log->trace("    propagation failed (!result.ok())");
0284     return nullptr;
0285   }
0286   m_log->trace("    propagation result is OK");
0287 
0288   // Pulling results to convenient variables
0289   auto trackStateParams  = *((*result).endParameters);
0290   const auto& parameter  = trackStateParams.parameters();
0291   const auto& covariance = *trackStateParams.covariance();
0292 
0293   // Path length
0294   const float pathLength      = initPathLength + (*result).pathLength;
0295   const float pathLengthError = 0;
0296   m_log->trace("    path len = {}", pathLength);
0297 
0298   // Position:
0299   auto projectionPos = trackStateParams.position(m_geoContext);
0300   const decltype(edm4eic::TrackPoint::position) position{static_cast<float>(projectionPos(0)),
0301                                                          static_cast<float>(projectionPos(1)),
0302                                                          static_cast<float>(projectionPos(2))};
0303   const decltype(edm4eic::TrackPoint::positionError) positionError{0, 0, 0};
0304   m_log->trace("    pos x = {}", position.x);
0305   m_log->trace("    pos y = {}", position.y);
0306   m_log->trace("    pos z = {}", position.z);
0307 
0308   // Momentum
0309   const decltype(edm4eic::TrackPoint::momentum) momentum = edm4hep::utils::sphericalToVector(
0310       static_cast<float>(1.0 / std::abs(parameter[Acts::eBoundQOverP])),
0311       static_cast<float>(parameter[Acts::eBoundTheta]),
0312       static_cast<float>(parameter[Acts::eBoundPhi]));
0313   const decltype(edm4eic::TrackPoint::momentumError) momentumError{
0314       static_cast<float>(covariance(Acts::eBoundTheta, Acts::eBoundTheta)),
0315       static_cast<float>(covariance(Acts::eBoundPhi, Acts::eBoundPhi)),
0316       static_cast<float>(covariance(Acts::eBoundQOverP, Acts::eBoundQOverP)),
0317       static_cast<float>(covariance(Acts::eBoundTheta, Acts::eBoundPhi)),
0318       static_cast<float>(covariance(Acts::eBoundTheta, Acts::eBoundQOverP)),
0319       static_cast<float>(covariance(Acts::eBoundPhi, Acts::eBoundQOverP))};
0320 
0321   // time
0322   const float time{static_cast<float>(parameter(Acts::eBoundTime))};
0323   const float timeError{static_cast<float>(sqrt(covariance(Acts::eBoundTime, Acts::eBoundTime)))};
0324 
0325   // Direction
0326   const float theta(parameter[Acts::eBoundTheta]);
0327   const float phi(parameter[Acts::eBoundPhi]);
0328   const decltype(edm4eic::TrackPoint::directionError) directionError{
0329       static_cast<float>(covariance(Acts::eBoundTheta, Acts::eBoundTheta)),
0330       static_cast<float>(covariance(Acts::eBoundPhi, Acts::eBoundPhi)),
0331       static_cast<float>(covariance(Acts::eBoundTheta, Acts::eBoundPhi))};
0332 
0333   // >oO debug print
0334   m_log->trace("    loc 0   = {:.4f}", parameter[Acts::eBoundLoc0]);
0335   m_log->trace("    loc 1   = {:.4f}", parameter[Acts::eBoundLoc1]);
0336   m_log->trace("    phi     = {:.4f}", parameter[Acts::eBoundPhi]);
0337   m_log->trace("    theta   = {:.4f}", parameter[Acts::eBoundTheta]);
0338   m_log->trace("    q/p     = {:.4f}", parameter[Acts::eBoundQOverP]);
0339   m_log->trace("    p       = {:.4f}", 1.0 / parameter[Acts::eBoundQOverP]);
0340   m_log->trace("    err phi = {:.4f}", sqrt(covariance(Acts::eBoundPhi, Acts::eBoundPhi)));
0341   m_log->trace("    err th  = {:.4f}", sqrt(covariance(Acts::eBoundTheta, Acts::eBoundTheta)));
0342   m_log->trace("    err q/p = {:.4f}", sqrt(covariance(Acts::eBoundQOverP, Acts::eBoundQOverP)));
0343   m_log->trace("    chi2    = {:.4f}", trajState.chi2Sum);
0344   m_log->trace("    loc err = {:.4f}",
0345                static_cast<float>(covariance(Acts::eBoundLoc0, Acts::eBoundLoc0)));
0346   m_log->trace("    loc err = {:.4f}",
0347                static_cast<float>(covariance(Acts::eBoundLoc1, Acts::eBoundLoc1)));
0348   m_log->trace("    loc err = {:.4f}",
0349                static_cast<float>(covariance(Acts::eBoundLoc0, Acts::eBoundLoc1)));
0350 
0351   uint64_t surface = targetSurf->geometryId().value();
0352   uint32_t system  = 0; // default value...will be set in TrackPropagation factory
0353 
0354   return std::make_unique<edm4eic::TrackPoint>(
0355       edm4eic::TrackPoint{surface, system, position, positionError, momentum, momentumError, time,
0356                           timeError, theta, phi, directionError, pathLength, pathLengthError});
0357 }
0358 
0359 } // namespace eicrecon
[ Source navigation ] [ Diff markup ] [ Identifier search ] [ general search ]

This page was automatically generated by the 2.3.7 LXR engine.

The LXR team
Valid CSS 2.1! Valid HTML 4.01!