Back to home page

EIC code displayed by LXR
 
 

    

File indexing completed on 2026-09-20 08:27:28 

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/Utilities/MathHelpers.hpp>
0009 #if Acts_VERSION_MAJOR >= 46
0010 #include <Acts/EventData/BoundTrackParameters.hpp>
0011 #else
0012 #include <Acts/EventData/GenericBoundTrackParameters.hpp>
0013 #endif
0014 #include <Acts/EventData/MultiTrajectoryHelpers.hpp>
0015 #include <Acts/EventData/VectorMultiTrajectory.hpp>
0016 #include <Acts/Geometry/GeometryIdentifier.hpp>
0017 #include <Acts/Geometry/TrackingGeometry.hpp>
0018 #include <Acts/MagneticField/MagneticFieldProvider.hpp>
0019 #include <Acts/Material/MaterialInteraction.hpp>
0020 #include <Acts/Propagator/ActorList.hpp>
0021 #include <Acts/Propagator/EigenStepper.hpp>
0022 #include <Acts/Propagator/MaterialInteractor.hpp>
0023 #include <Acts/Propagator/Navigator.hpp>
0024 #include <Acts/Propagator/Propagator.hpp>
0025 #include <Acts/Propagator/PropagatorResult.hpp>
0026 #include <Acts/Surfaces/CylinderBounds.hpp>
0027 #include <Acts/Surfaces/CylinderSurface.hpp>
0028 #include <Acts/Surfaces/DiscSurface.hpp>
0029 #include <Acts/Surfaces/RadialBounds.hpp>
0030 #include <Acts/Utilities/Logger.hpp>
0031 #include <DD4hep/Handle.h>
0032 #include <Evaluator/DD4hepUnits.h>
0033 #include <edm4eic/Cov2f.h>
0034 #include <edm4eic/Cov3f.h>
0035 #include <edm4hep/Vector3f.h>
0036 #include <edm4hep/utils/vector_utils.h>
0037 #include <spdlog/common.h>
0038 #include <Eigen/Core>
0039 #include <Eigen/Geometry>
0040 #include <algorithm>
0041 #include <any>
0042 #include <cmath>
0043 #include <cstdint>
0044 #include <functional>
0045 #include <iterator>
0046 #include <map>
0047 #include <optional>
0048 #include <stdexcept>
0049 #include <string>
0050 #include <tuple>
0051 #include <typeinfo>
0052 #include <utility>
0053 #include <variant>
0054 
0055 #include "algorithms/interfaces/detail/multilambda.h"
0056 #include "algorithms/tracking/ActsGeometryProvider.h"
0057 #include "algorithms/tracking/TrackPropagation.h"
0058 #include "algorithms/tracking/TrackPropagationConfig.h"
0059 #include "extensions/spdlog/SpdlogToActs.h"
0060 
0061 namespace eicrecon {
0062 
0063 void TrackPropagation::init() {
0064   const auto* detector = m_detector;
0065 
0066   std::map<uint32_t, std::size_t> system_id_layers;
0067 
0068   multilambda _toDouble = {
0069       [](const std::string& v) { return dd4hep::_toDouble(v); },
0070       [](const double& v) { return v; },
0071   };
0072 
0073   auto _toActsSurface =
0074       [&_toDouble, &detector, &system_id_layers](
0075           const std::variant<CylinderSurfaceConfig, DiscSurfaceConfig> surface_variant)
0076       -> std::shared_ptr<Acts::Surface> {
0077     if (std::holds_alternative<CylinderSurfaceConfig>(surface_variant)) {
0078       CylinderSurfaceConfig surface = std::get<CylinderSurfaceConfig>(surface_variant);
0079       const double rmin =
0080           std::visit(_toDouble, surface.rmin) / dd4hep::mm * Acts::UnitConstants::mm;
0081       const double zmin =
0082           std::visit(_toDouble, surface.zmin) / dd4hep::mm * Acts::UnitConstants::mm;
0083       const double zmax =
0084           std::visit(_toDouble, surface.zmax) / dd4hep::mm * Acts::UnitConstants::mm;
0085       const uint32_t system_id = detector->constant<uint32_t>(surface.id);
0086       auto bounds              = std::make_shared<Acts::CylinderBounds>(rmin, (zmax - zmin) / 2);
0087       auto t                   = Acts::Translation3(Acts::Vector3(0, 0, (zmax + zmin) / 2));
0088       auto tf                  = Acts::Transform3(t);
0089       auto acts_surface        = Acts::Surface::makeShared<Acts::CylinderSurface>(tf, bounds);
0090       acts_surface->assignGeometryId(
0091           Acts::GeometryIdentifier().withExtra(system_id).withLayer(++system_id_layers[system_id]));
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       acts_surface->assignGeometryId(
0108           Acts::GeometryIdentifier().withExtra(system_id).withLayer(++system_id_layers[system_id]));
0109       return acts_surface;
0110     }
0111     throw std::domain_error("Unknown surface type");
0112   };
0113   m_target_surfaces.resize(m_cfg.target_surfaces.size());
0114   std::ranges::transform(m_cfg.target_surfaces, m_target_surfaces.begin(), _toActsSurface);
0115   m_filter_surfaces.resize(m_cfg.filter_surfaces.size());
0116   std::ranges::transform(m_cfg.filter_surfaces, m_filter_surfaces.begin(), _toActsSurface);
0117 
0118   trace("Initialized");
0119 }
0120 
0121 void TrackPropagation::propagateToSurfaceList(const Input& input, const Output& output) const {
0122   const auto [tracks, track_states, tracks_acts] = input;
0123   auto [track_segments]                          = output;
0124 
0125   // logging
0126   trace("Propagate tracks: --------------------");
0127   trace("number of tracks: {}", tracks->size());
0128 
0129   // Construct ConstTrackContainer from underlying containers
0130   auto trackStateContainer = std::make_shared<Acts::ConstVectorMultiTrajectory>(*track_states);
0131   auto trackContainer      = std::make_shared<Acts::ConstVectorTrackContainer>(*tracks_acts);
0132   ActsExamples::ConstTrackContainer constTracks(trackContainer, trackStateContainer);
0133 
0134   // loop over input tracks
0135   std::size_t i = 0;
0136   for (const auto& track : constTracks) {
0137 
0138     // check if this track can be propagated to any filter surface
0139     bool track_reaches_filter_surface{false};
0140     for (const auto& filter_surface : m_filter_surfaces) {
0141       auto point = propagate(edm4eic::Track{}, track, constTracks, filter_surface);
0142       if (point) {
0143         track_reaches_filter_surface = true;
0144         break;
0145       }
0146     }
0147     if (!track_reaches_filter_surface) {
0148       ++i;
0149       continue;
0150     }
0151 
0152     // start a mutable TrackSegment
0153     auto track_segment = track_segments->create();
0154 
0155     // corresponding track
0156     if (tracks->size() == constTracks.size()) {
0157       trace("track segment connected to track {}", i);
0158       track_segment.setTrack((*tracks)[i]);
0159       ++i;
0160     }
0161 
0162     // zero measurements of segment length
0163     decltype(edm4eic::TrackSegmentData::length) length            = 0;
0164     decltype(edm4eic::TrackSegmentData::lengthError) length_error = 0;
0165 
0166     // loop over projection-target surfaces
0167     for (const auto& target_surface : m_target_surfaces) {
0168 
0169       // project the track to this surface
0170       auto point = propagate(edm4eic::Track{}, track, constTracks, target_surface);
0171       if (!point) {
0172         trace("<> Failed to propagate track to this plane");
0173         continue;
0174       }
0175 
0176       // logging
0177       trace("<> track: x=( {:>10.2f} {:>10.2f} {:>10.2f} )", point->position.x, point->position.y,
0178             point->position.z);
0179       trace("               p=( {:>10.2f} {:>10.2f} {:>10.2f} )", point->momentum.x,
0180             point->momentum.y, point->momentum.z);
0181 
0182       // track point cut
0183       if (!m_cfg.track_point_cut(*point)) {
0184         trace("                 => REJECTED by trackPointCut");
0185         if (m_cfg.skip_track_on_track_point_cut_failure) {
0186           break;
0187         }
0188         continue;
0189       }
0190 
0191       // update the `TrackSegment` length
0192       // FIXME: `length` and `length_error` are currently not used by any callers, and may not be correctly calculated here
0193       if (track_segment.points_size() > 0) {
0194         auto pos0 = point->position;
0195         auto pos1 = std::prev(track_segment.points_end())->position;
0196         auto dist = edm4hep::utils::magnitude(pos0 - pos1);
0197         length += dist;
0198         trace("               dist to previous point: {}", dist);
0199       }
0200 
0201       // add the `TrackPoint` to the `TrackSegment`
0202       track_segment.addToPoints(*point);
0203 
0204     } // end `targetSurfaces` loop
0205 
0206     // set final length and length error
0207     track_segment.setLength(length);
0208     track_segment.setLengthError(length_error);
0209 
0210   } // end loop over input tracks
0211 }
0212 
0213 std::unique_ptr<edm4eic::TrackPoint>
0214 TrackPropagation::propagate(const edm4eic::Track& /* track */,
0215                             const ActsExamples::ConstTrackProxy& acts_track,
0216                             const ActsExamples::ConstTrackContainer& trackContainer,
0217                             const std::shared_ptr<const Acts::Surface>& targetSurf) const {
0218 
0219   auto tipIndex = acts_track.tipIndex();
0220 
0221   trace("  Propagating track with tip index {}", tipIndex);
0222 
0223   // Collect the trajectory summary info
0224   auto trajState =
0225       Acts::MultiTrajectoryHelpers::trajectoryState(trackContainer.trackStateContainer(), tipIndex);
0226   int m_nMeasurements = trajState.nMeasurements;
0227   int m_nStates       = trajState.nStates;
0228 
0229   trace("  Num measurement in trajectory: {}", m_nMeasurements);
0230   trace("  Num states in trajectory     : {}", m_nStates);
0231 
0232   // Get track state at last measurement surface
0233   // For last measurement surface, filtered and smoothed results are equivalent
0234   auto trackState        = trackContainer.trackStateContainer().getTrackState(tipIndex);
0235   auto initSurface       = trackState.referenceSurface().getSharedPtr();
0236   const auto& initParams = trackState.filtered();
0237   const auto& initCov    = trackState.filteredCovariance();
0238 
0239   Acts::BoundTrackParameters initBoundParams(initSurface, initParams, initCov,
0240                                              acts_track.particleHypothesis());
0241 
0242   // Get pathlength of last track state with respect to perigee surface
0243   const auto initPathLength = trackState.pathLength();
0244 
0245   trace("    TrackPropagation. Propagating to surface # {}", typeid(targetSurf->type()).name());
0246 
0247   std::shared_ptr<const Acts::TrackingGeometry> trackingGeometry   = m_geoSvc->trackingGeometry();
0248   std::shared_ptr<const Acts::MagneticFieldProvider> magneticField = m_geoSvc->getFieldProvider();
0249 
0250   // Convert algorithm log level to Acts log level
0251   const auto spdlog_level = static_cast<spdlog::level::level_enum>(this->level());
0252   const auto acts_level   = eicrecon::SpdlogToActsLevel(spdlog_level);
0253   ACTS_LOCAL_LOGGER(Acts::getDefaultLogger("PROP", acts_level));
0254 
0255   using Propagator        = Acts::Propagator<Acts::EigenStepper<>, Acts::Navigator>;
0256   using PropagatorOptions = Propagator::template Options<Acts::ActorList<Acts::MaterialInteractor>>;
0257   Propagator propagator(Acts::EigenStepper<>(magneticField),
0258                         Acts::Navigator({.trackingGeometry = m_geoSvc->trackingGeometry()},
0259                                         logger().cloneWithSuffix("Navigator")),
0260                         logger().cloneWithSuffix("Propagator"));
0261 
0262   // Get run-scoped contexts from service
0263   const auto& gctx = m_geoSvc->getActsGeometryContext();
0264   const auto& mctx = m_geoSvc->getActsMagneticFieldContext();
0265 
0266   PropagatorOptions propagationOptions(gctx, mctx);
0267 
0268   // Some target surfaces (e.g. DIRC) may be inside the last measurement surface (e.g. BIC),
0269   // so we use a straight line intersection from the last measurement surface to the target
0270   // surface to determine if we have to propagate backwards.
0271   auto initPosition  = initBoundParams.position(gctx);
0272   auto initDirection = initBoundParams.direction();
0273   auto intersections = targetSurf->intersect(gctx, initPosition, initDirection);
0274 
0275   // Determine closest forward intersection (positive pathlength from perigee)
0276   auto intersection = intersections.closestForward();
0277   auto difference   = intersection.position() - initPosition;
0278   auto dot          = difference.dot(initBoundParams.direction());
0279 
0280   // Propagate forwards by default
0281   propagationOptions.direction = Acts::Direction::Forward();
0282 
0283   // but invert if the position difference is opposite to direction
0284   if (intersection.isValid() && dot < 0) {
0285 
0286     // The extra fields of the surface geometry ID contain the DD4hep system
0287     auto initSurfaceExtra   = initSurface->geometryId().extra();
0288     auto targetSurfaceExtra = targetSurf->geometryId().extra();
0289     debug("    inverting direction for propagator from surface {} to {}", initSurfaceExtra,
0290           targetSurfaceExtra);
0291     auto p1 = initBoundParams.position(gctx);
0292     debug("      initial position {} {} {}", p1.x(), p1.y(), p1.z());
0293     auto p2 = intersection.position();
0294     debug("      straight line intersection at {} {} {}", p2.x(), p2.y(), p2.z());
0295 
0296     // Propagate backwards
0297     propagationOptions.direction = Acts::Direction::Backward();
0298   }
0299 
0300   auto result = propagator.propagate(initBoundParams, *targetSurf, propagationOptions);
0301 
0302   // check propagation result
0303   if (!result.ok()) {
0304     trace("    propagation failed (!result.ok())");
0305     return nullptr;
0306   }
0307   trace("    propagation result is OK");
0308 
0309   // Pulling results to convenient variables
0310   auto trackStateParams  = *((*result).endParameters);
0311   const auto& parameter  = trackStateParams.parameters();
0312   const auto& covariance = *trackStateParams.covariance();
0313 
0314   // Path length
0315   const float pathLength      = initPathLength + (*result).pathLength;
0316   const float pathLengthError = 0;
0317   trace("    path len = {}", pathLength);
0318 
0319   // Position:
0320   auto projectionPos = trackStateParams.position(gctx);
0321   const decltype(edm4eic::TrackPoint::position) position{static_cast<float>(projectionPos(0)),
0322                                                          static_cast<float>(projectionPos(1)),
0323                                                          static_cast<float>(projectionPos(2))};
0324   const decltype(edm4eic::TrackPoint::positionError) positionError{0, 0, 0};
0325   trace("    pos x = {}", position.x);
0326   trace("    pos y = {}", position.y);
0327   trace("    pos z = {}", position.z);
0328 
0329   // Momentum
0330   const decltype(edm4eic::TrackPoint::momentum) momentum = edm4hep::utils::sphericalToVector(
0331       static_cast<float>(1.0 / std::abs(parameter[Acts::eBoundQOverP] * Acts::UnitConstants::GeV)),
0332       static_cast<float>(parameter[Acts::eBoundTheta] / Acts::UnitConstants::rad),
0333       static_cast<float>(parameter[Acts::eBoundPhi] / Acts::UnitConstants::rad));
0334   const decltype(edm4eic::TrackPoint::momentumError) momentumError{
0335       static_cast<float>(covariance(Acts::eBoundTheta, Acts::eBoundTheta) /
0336                          Acts::UnitConstants::rad / Acts::UnitConstants::rad),
0337       static_cast<float>(covariance(Acts::eBoundPhi, Acts::eBoundPhi) / Acts::UnitConstants::rad /
0338                          Acts::UnitConstants::rad),
0339       static_cast<float>(covariance(Acts::eBoundQOverP, Acts::eBoundQOverP) *
0340                          Acts::UnitConstants::GeV * Acts::UnitConstants::GeV),
0341       static_cast<float>(covariance(Acts::eBoundTheta, Acts::eBoundPhi) / Acts::UnitConstants::rad /
0342                          Acts::UnitConstants::rad),
0343       static_cast<float>(covariance(Acts::eBoundTheta, Acts::eBoundQOverP) /
0344                          Acts::UnitConstants::rad * Acts::UnitConstants::GeV),
0345       static_cast<float>(covariance(Acts::eBoundPhi, Acts::eBoundQOverP) /
0346                          Acts::UnitConstants::rad * Acts::UnitConstants::GeV)};
0347 
0348   // time
0349   const float time{static_cast<float>(parameter(Acts::eBoundTime) / Acts::UnitConstants::ns)};
0350   const float timeError{static_cast<float>(sqrt(covariance(Acts::eBoundTime, Acts::eBoundTime)) /
0351                                            Acts::UnitConstants::ns)};
0352 
0353   // Direction
0354   const float theta(parameter[Acts::eBoundTheta]);
0355   const float phi(parameter[Acts::eBoundPhi]);
0356   const decltype(edm4eic::TrackPoint::directionError) directionError{
0357       static_cast<float>(covariance(Acts::eBoundTheta, Acts::eBoundTheta) /
0358                          Acts::UnitConstants::rad / Acts::UnitConstants::rad),
0359       static_cast<float>(covariance(Acts::eBoundPhi, Acts::eBoundPhi) / Acts::UnitConstants::rad /
0360                          Acts::UnitConstants::rad),
0361       static_cast<float>(covariance(Acts::eBoundTheta, Acts::eBoundPhi) / Acts::UnitConstants::rad /
0362                          Acts::UnitConstants::rad)};
0363 
0364   // >oO debug print
0365   trace("    loc 0   = {:.4f} mm", parameter[Acts::eBoundLoc0] / Acts::UnitConstants::mm);
0366   trace("    loc 1   = {:.4f} mm", parameter[Acts::eBoundLoc1] / Acts::UnitConstants::mm);
0367   trace("    phi     = {:.4f} rad", parameter[Acts::eBoundPhi] / Acts::UnitConstants::rad);
0368   trace("    theta   = {:.4f} rad", parameter[Acts::eBoundTheta] / Acts::UnitConstants::rad);
0369   trace("    q/p     = {:.4f} / GeV", parameter[Acts::eBoundQOverP] * Acts::UnitConstants::GeV);
0370   trace("    p       = {:.4f} GeV", 1.0 / parameter[Acts::eBoundQOverP] / Acts::UnitConstants::GeV);
0371   trace("    err phi = {:.4f} rad", sqrt(covariance(Acts::eBoundPhi, Acts::eBoundPhi) /
0372                                          Acts::UnitConstants::rad / Acts::UnitConstants::rad));
0373   trace("    err th  = {:.4f} rad", sqrt(covariance(Acts::eBoundTheta, Acts::eBoundTheta) /
0374                                          Acts::UnitConstants::rad / Acts::UnitConstants::rad));
0375   trace("    err q/p = {:.4f} / GeV", sqrt(covariance(Acts::eBoundQOverP, Acts::eBoundQOverP) *
0376                                            Acts::UnitConstants::GeV * Acts::UnitConstants::GeV));
0377   trace("    chi2    = {:.4f}", trajState.chi2Sum);
0378   trace("    loc err = {:.4f} mm^2",
0379         static_cast<float>(covariance(Acts::eBoundLoc0, Acts::eBoundLoc0) /
0380                            Acts::UnitConstants::mm / Acts::UnitConstants::mm));
0381   trace("    loc err = {:.4f} mm^2",
0382         static_cast<float>(covariance(Acts::eBoundLoc1, Acts::eBoundLoc1) /
0383                            Acts::UnitConstants::mm / Acts::UnitConstants::mm));
0384   trace("    loc err = {:.4f} mm^2",
0385         static_cast<float>(covariance(Acts::eBoundLoc0, Acts::eBoundLoc1) /
0386                            Acts::UnitConstants::mm / Acts::UnitConstants::mm));
0387 
0388   uint64_t surface = targetSurf->geometryId().value();
0389   uint32_t system  = 0; // default value...will be set in TrackPropagation factory
0390 
0391   return std::make_unique<edm4eic::TrackPoint>(
0392       edm4eic::TrackPoint{.surface         = surface,
0393                           .system          = system,
0394                           .position        = position,
0395                           .positionError   = positionError,
0396                           .momentum        = momentum,
0397                           .momentumError   = momentumError,
0398                           .time            = time,
0399                           .timeError       = timeError,
0400                           .theta           = theta,
0401                           .phi             = phi,
0402                           .directionError  = directionError,
0403                           .pathlength      = pathLength,
0404                           .pathlengthError = pathLengthError});
0405 }
0406 
0407 } // 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!