diff --git a/CMakeLists.txt b/CMakeLists.txt index 45da9f68b..531ea5daf 100644 --- a/CMakeLists.txt +++ b/CMakeLists.txt @@ -71,6 +71,31 @@ find_package( atomic regex) +# Find Eigen3 (required by ACTS) +find_package(Eigen3 CONFIG REQUIRED) + +# GNUInstallDirs must be included before ACTS is configured as a sub-project; +# ACTS only calls include(GNUInstallDirs) when it is the top-level project, so +# CMAKE_INSTALL_INCLUDEDIR would otherwise be empty and its install() calls fail. +include(GNUInstallDirs) + +# Find or fetch ACTS - must be at top level so both Ecal and Tracking can use it +find_package(Acts 47.0.0 QUIET) +if (NOT Acts_FOUND) + message(STATUS "Did not find Acts >= 47.0.0, downloading and compiling v47.0.0") + include(FetchContent) + FetchContent_Declare( + Acts + URL https://github.com/acts-project/acts/archive/refs/tags/v47.0.0.tar.gz + URL_HASH MD5=833c02edee49827f25be74c80d8ee654 + DOWNLOAD_EXTRACT_TIMESTAMP TRUE + ) + FetchContent_MakeAvailable(Acts) + FetchContent_GetProperties(Acts) + include_directories(SYSTEM "${Acts_SOURCE_DIR}/Core/include") +else() + message(STATUS "Found Acts ${Acts_VERSION}") +endif() # FunctionalCoreTest.cxx doesnt comply with clang-tidy but that's ok set_source_files_properties(Framework/test/FunctionalCoreTest.cxx PROPERTIES COMPILE_OPTIONS "-Wno-null-dereference") diff --git a/Ecal/CMakeLists.txt b/Ecal/CMakeLists.txt index 59c57dfe4..aa0923be5 100644 --- a/Ecal/CMakeLists.txt +++ b/Ecal/CMakeLists.txt @@ -17,15 +17,13 @@ setup_library( find_package(ONNXRuntime 1.2.0) find_package(Eigen3 CONFIG REQUIRED) - -# Find ACTS for track finding in ECAL -find_package(Acts REQUIRED COMPONENTS Core) +# Acts is set up by the top-level CMakeLists.txt setup_library(module Ecal dependencies ROOT::Physics Framework::Framework Recon::Event Recon::Recon Tools::Tools DetDescr::DetDescr ONNXRuntime::Interface Ecal::Event Tracking::Event Tracking::Tracking SimCore::Event - Eigen3::Eigen ActsCore + Eigen3::Eigen Acts::Core ) setup_test(dependencies Ecal::Ecal) diff --git a/Ecal/include/Ecal/EcalTrackFinderProcessor.h b/Ecal/include/Ecal/EcalTrackFinderProcessor.h index 6318499a9..9f7ced6f8 100644 --- a/Ecal/include/Ecal/EcalTrackFinderProcessor.h +++ b/Ecal/include/Ecal/EcalTrackFinderProcessor.h @@ -35,6 +35,8 @@ #include "Acts/Surfaces/RectangleBounds.hpp" #include "Acts/Surfaces/Surface.hpp" #include "Acts/TrackFinding/CombinatorialKalmanFilter.hpp" +#include "Acts/TrackFinding/MeasurementSelector.hpp" +#include "Acts/TrackFinding/TrackStateCreator.hpp" #include "Acts/TrackFitting/GainMatrixUpdater.hpp" // C++ diff --git a/Ecal/src/Ecal/EcalTrackFinderProcessor.cxx b/Ecal/src/Ecal/EcalTrackFinderProcessor.cxx index 55c154699..6ad81a925 100644 --- a/Ecal/src/Ecal/EcalTrackFinderProcessor.cxx +++ b/Ecal/src/Ecal/EcalTrackFinderProcessor.cxx @@ -14,8 +14,8 @@ // ACTS #include "Acts/Definitions/Units.hpp" +#include "Acts/EventData/BoundTrackParameters.hpp" #include "Acts/EventData/MultiTrajectory.hpp" -#include "Acts/EventData/TrackParameters.hpp" #include "Acts/EventData/TransformationHelpers.hpp" #include "Acts/Geometry/CuboidVolumeBuilder.hpp" #include "Acts/Geometry/GeometryContext.hpp" @@ -23,8 +23,7 @@ #include "Acts/Geometry/TrackingGeometryBuilder.hpp" #include "Acts/Geometry/TrackingVolume.hpp" #include "Acts/MagneticField/MagneticFieldContext.hpp" -#include "Acts/Propagator/AbortList.hpp" -#include "Acts/Propagator/ActionList.hpp" +#include "Acts/Propagator/ActorList.hpp" #include "Acts/Propagator/MaterialInteractor.hpp" #include "Acts/Propagator/StandardAborters.hpp" #include "Acts/Propagator/detail/SteppingLogger.hpp" @@ -121,7 +120,9 @@ void EcalTrackFinderProcessor::onNewRun(const ldmx::RunHeader&) { Acts::Vector3 back_acts = tracking::sim::utils::ldmx2Acts(Acts::Vector3(0.0, 0.0, ecal_back_z)); - // Create layer configurations - one per ECAL layer + // Create layer configurations - one per ECAL layer. + // IMPORTANT: layer_configs must be in ascending x-order so LayerArrayCreator + // gets a monotone sequence for BinningType::arbitrary along AxisX. std::vector layer_configs; double clearance = 1.0; // mm envelope around each layer surface @@ -143,7 +144,8 @@ void EcalTrackFinderProcessor::onNewRun(const ldmx::RunHeader&) { ecal_vol_cfg.name = "EcalVolume"; ecal_vol_cfg.layerCfg = layer_configs; ecal_vol_cfg.volumeMaterial = - std::make_shared(Acts::Material()); + std::make_shared( + Acts::Material::Vacuum()); // Build the tracking geometry Acts::CuboidVolumeBuilder cvb; @@ -173,6 +175,11 @@ void EcalTrackFinderProcessor::onNewRun(const ldmx::RunHeader&) { // Only care about sensitive surfaces (those with a sensitive ID) if (surface->geometryId().sensitive() == 0) return; + // CuboidVolumeBuilder creates new PlaneSurface objects inside the geometry + // rather than using our originals. Mark them sensitive here so the CKF + // actor doesn't skip them as passive surfaces. + const_cast(surface)->assignIsSensitive(true); + // Match to ECAL layer by z position (LDMX frame) Acts::Vector3 center_ldmx = tracking::sim::utils::acts2Ldmx(surface->center(gctx)); @@ -197,16 +204,19 @@ void EcalTrackFinderProcessor::onNewRun(const ldmx::RunHeader&) { // Setup stepper with zero B-field const auto stepper = Acts::EigenStepper<>{zero_b_field}; + auto acts_logging_level = + debug_ ? Acts::Logging::VERBOSE : Acts::Logging::FATAL; + // Setup navigator with tracking geometry Acts::Navigator::Config nav_cfg{tracking_geometry_}; nav_cfg.resolveSensitive = true; nav_cfg.resolvePassive = false; nav_cfg.resolveMaterial = false; - const Acts::Navigator navigator(nav_cfg); + const Acts::Navigator navigator( + nav_cfg, + Acts::getDefaultLogger("ECAL_NAV", acts_logging_level)); // Create propagator - auto acts_logging_level = - debug_ ? Acts::Logging::VERBOSE : Acts::Logging::WARNING; propagator_ = std::make_unique( stepper, navigator, Acts::getDefaultLogger("ECAL_PROP", acts_logging_level)); @@ -239,10 +249,14 @@ void EcalTrackFinderProcessor::createEcalSurfaces() { auto surface = Acts::Surface::makeShared(transform, bounds); + // Mark as sensitive so the CKF actor creates track states here. + // PlaneSurfaces default to isSensitive()=false; without this flag the CKF + // treats them as passive material surfaces and creates no track states. + surface->assignIsSensitive(true); + // Assign a geometry ID (use layer as volume, 0 as layer in ACTS sense) - Acts::GeometryIdentifier geo_id; - geo_id.setVolume(layer); - geo_id.setLayer(0); + Acts::GeometryIdentifier geo_id = + Acts::GeometryIdentifier().withVolume(layer).withLayer(0); surface->assignGeometryId(geo_id); layer_surfaces_[layer] = surface; @@ -395,17 +409,20 @@ std::vector EcalTrackFinderProcessor::findSeeds( return seeds; } - // Create seed track at reference surface (ECAL front) - // Intersect line with reference surface + // Create seed track at the FIRST ECAL layer surface (layer 0). + // Using a surface that IS in the tracking geometry (has associatedLayer set) + // ensures the ACTS Navigator can use the fast initialization path and find + // all sensitive surfaces during CKF propagation. auto& gctx = getCondition( tracking::geo::GeometryContext::NAME) .get(); - // For a plane surface, normal is the third column of rotation (z-direction in - // local frame) + auto& seed_surface = layer_surfaces_.begin()->second; + + // For a plane surface, normal is the third column of rotation Acts::Vector3 ref_normal = - reference_surface_->transform(gctx).rotation().col(2); - Acts::Vector3 ref_center = reference_surface_->center(gctx); + seed_surface->localToGlobalTransform(gctx).rotation().col(2); + Acts::Vector3 ref_center = seed_surface->center(gctx); double t = (ref_center - position).dot(ref_normal) / direction.dot(ref_normal); @@ -416,14 +433,14 @@ std::vector EcalTrackFinderProcessor::findSeeds( Acts::Vector3 seed_mom = p_estimate * direction; // Charge (assume positive) - Acts::ActsScalar q = Acts::UnitConstants::e; + double q = Acts::UnitConstants::e; - // Convert to bound parameters at reference surface + // Convert to bound parameters at the first layer surface Acts::FreeVector seed_free = tracking::sim::utils::toFreeParameters(seed_pos, seed_mom, q); auto bound_params_result = Acts::transformFreeToBoundParameters( - seed_free, *reference_surface_, gctx); + seed_free, *seed_surface, gctx); if (!bound_params_result.ok()) { ldmx_log(warn) << "Failed to create bound parameters for seed"; @@ -441,12 +458,12 @@ std::vector EcalTrackFinderProcessor::findSeeds( stddev[Acts::eBoundQOverP] = 0.5 / p_estimate; // 50% uncertainty stddev[Acts::eBoundTime] = 10.0 * Acts::UnitConstants::ns; - Acts::BoundSquareMatrix bound_cov = stddev.cwiseProduct(stddev).asDiagonal(); + Acts::BoundMatrix bound_cov = stddev.cwiseProduct(stddev).asDiagonal(); // Create ldmx::Track seed ldmx::Track seed; - // Convert reference surface position to LDMX frame + // Convert layer 0 surface position to LDMX frame Acts::Vector3 ref_ldmx = tracking::sim::utils::acts2Ldmx(ref_center); seed.setPerigeeLocation(ref_ldmx[0], ref_ldmx[1], ref_ldmx[2]); @@ -486,7 +503,7 @@ EcalTrackFinderProcessor::makeGeoIdSourceLinkMap( // matches our source link map keys. auto it = layer_geo_ids_.find(layer); if (it == layer_geo_ids_.end()) { - ldmx_log(debug) << "No builder geometry ID for layer " << layer; + ldmx_log(warn) << "No builder geometry ID for layer " << layer; continue; } @@ -531,7 +548,9 @@ void EcalTrackFinderProcessor::produce(framework::Event& event) { // Create source link map auto geo_id_sl_map = makeGeoIdSourceLinkMap(measurements); - ldmx_log(debug) << "Source link map: " << geo_id_sl_map.size() << " entries"; + ldmx_log(info) << "Source link map: " << geo_id_sl_map.size() + << " entries from " << measurements.size() + << " measurements"; // Find seed tracks auto seed_tracks = findSeeds(measurements); @@ -568,18 +587,7 @@ void EcalTrackFinderProcessor::produce(framework::Event& event) { tracking::sim::LdmxMeasurementCalibrator calibrator{measurements}; - Acts::CombinatorialKalmanFilterExtensions ckf_extensions; - ckf_extensions.calibrator - .connect<&tracking::sim::LdmxMeasurementCalibrator::calibrate< - Acts::VectorMultiTrajectory>>(&calibrator); - ckf_extensions.updater.connect< - &Acts::GainMatrixUpdater::operator()>( - &kf_updater); - ckf_extensions.measurementSelector - .connect<&Acts::MeasurementSelector::select>( - &meas_sel); - - // Setup source link accessor + // Setup source link accessor iterator type and lambda struct SourceLinkAccIt { using BaseIt = decltype(geo_id_sl_map.begin()); BaseIt it_; @@ -613,11 +621,25 @@ void EcalTrackFinderProcessor::produce(framework::Event& event) { return {SourceLinkAccIt{begin}, SourceLinkAccIt{end}}; }; - Acts::SourceLinkAccessorDelegate - source_link_accessor_delegate; - source_link_accessor_delegate + // v46: calibrator and measurementSelector moved to TrackStateCreator + Acts::TrackStateCreator track_state_creator; + track_state_creator.sourceLinkAccessor .connect<&decltype(source_link_accessor)::operator(), decltype(source_link_accessor)>(&source_link_accessor); + track_state_creator.calibrator + .connect<&tracking::sim::LdmxMeasurementCalibrator::calibrate< + Acts::VectorMultiTrajectory>>(&calibrator); + track_state_creator.measurementSelector + .connect<&Acts::MeasurementSelector::select>( + &meas_sel); + + Acts::CombinatorialKalmanFilterExtensions ckf_extensions; + ckf_extensions.updater.connect< + &Acts::GainMatrixUpdater::operator()>( + &kf_updater); + ckf_extensions.createTrackStates.connect<&Acts::TrackStateCreator< + SourceLinkAccIt, TrackContainer>::createTrackStates>( + &track_state_creator); // Create track container Acts::VectorTrackContainer vtc; @@ -628,10 +650,10 @@ void EcalTrackFinderProcessor::produce(framework::Event& event) { for (size_t seed_idx = 0; seed_idx < seed_tracks.size(); ++seed_idx) { const auto& seed = seed_tracks[seed_idx]; - // Convert seed to BoundTrackParameters - // The seed was created with bound parameters on the reference PlaneSurface, - // so we must start the CKF from that same surface (NOT a PerigeeSurface, - // which interprets loc0/loc1 as d0/z0 instead of plane-local coordinates). + // Convert seed to BoundTrackParameters. + // Start from layer_surfaces_[0] which is part of the tracking geometry and + // has associatedLayer() set — this lets the Navigator use the fast + // initialization path and correctly traverse all 32 ECAL layers. Acts::BoundVector param_vec; param_vec << seed.getD0(), seed.getZ0(), seed.getPhi(), seed.getTheta(), seed.getQoP(), seed.getT(); @@ -640,18 +662,17 @@ void EcalTrackFinderProcessor::produce(framework::Event& event) { << " loc1=" << param_vec[1] << " phi=" << param_vec[2] << " theta=" << param_vec[3] << " qop=" << param_vec[4]; - Acts::BoundSquareMatrix cov_mat = + Acts::BoundMatrix cov_mat = tracking::sim::utils::unpackCov(seed.getPerigeeCov()); - auto part_hypo{Acts::SinglyChargedParticleHypothesis::electron()}; - Acts::BoundTrackParameters start_params(reference_surface_, param_vec, + auto part_hypo{Acts::ParticleHypothesis::electron()}; + auto& layer0_surface = layer_surfaces_.begin()->second; + Acts::BoundTrackParameters start_params(layer0_surface, param_vec, cov_mat, part_hypo); // Setup CKF options - const Acts::CombinatorialKalmanFilterOptions - ckf_options(gctx, mctx, cctx, source_link_accessor_delegate, - ckf_extensions, propagator_options); + const Acts::CombinatorialKalmanFilterOptions ckf_options( + gctx, mctx, cctx, ckf_extensions, propagator_options); // Run CKF auto results = ckf_->findTracks(start_params, ckf_options, tc); @@ -663,11 +684,28 @@ void EcalTrackFinderProcessor::produce(framework::Event& event) { } auto& tracks_from_seed = results.value(); - ldmx_log(debug) << "CKF returned " << tracks_from_seed.size() - << " tracks from seed " << seed_idx; + ldmx_log(info) << "CKF returned " << tracks_from_seed.size() + << " tracks from seed " << seed_idx; for (auto& track : tracks_from_seed) { + // Count track state types before smoothing + int n_meas = 0, n_holes = 0, n_outliers = 0, n_total = 0; + for (const auto& ts : track.trackStatesReversed()) { + ++n_total; + if (ts.typeFlags().isMeasurement()) ++n_meas; + if (ts.typeFlags().isHole()) ++n_holes; + if (ts.typeFlags().isOutlier()) ++n_outliers; + } + ldmx_log(info) << "Track states: total=" << n_total + << " meas=" << n_meas << " holes=" << n_holes + << " outliers=" << n_outliers; + // Smooth the track - Acts::smoothTrack(gctx, track); + auto smooth_result = Acts::smoothTrack(gctx, track); + if (!smooth_result.ok()) { + ldmx_log(warn) << "smoothTrack failed: " + << smooth_result.error().message(); + continue; + } // Create output track ldmx::Track trk; @@ -689,7 +727,7 @@ void EcalTrackFinderProcessor::produce(framework::Event& event) { } if (!found_smoothed) { - ldmx_log(warn) << "No smoothed track state found"; + ldmx_log(warn) << "No smoothed track state found after smoothing"; continue; } @@ -755,7 +793,7 @@ void EcalTrackFinderProcessor::produce(framework::Event& event) { trk.setPerigeeCov(cov_vec); Acts::Vector3 ref_loc_ldmx = - tracking::sim::utils::acts2Ldmx(reference_surface_->center(gctx)); + tracking::sim::utils::acts2Ldmx(layer_surfaces_.begin()->second->center(gctx)); trk.setPerigeeLocation(ref_loc_ldmx[0], ref_loc_ldmx[1], ref_loc_ldmx[2]); trk.setChi2(track.chi2()); @@ -766,11 +804,10 @@ void EcalTrackFinderProcessor::produce(framework::Event& event) { // Add measurement indices for (const auto ts : track.trackStatesReversed()) { - if (ts.typeFlags().test(Acts::TrackStateFlag::MeasurementFlag) && - ts.hasUncalibratedSourceLink()) { - const acts_examples::IndexSourceLink sl = - ts.getUncalibratedSourceLink() - .template get(); + if (ts.typeFlags().isMeasurement() && ts.hasUncalibratedSourceLink()) { + Acts::SourceLink usl = ts.getUncalibratedSourceLink(); + const acts_examples::IndexSourceLink& sl = + usl.get(); trk.addMeasurementIndex(sl.index()); } } @@ -832,11 +869,11 @@ void EcalTrackFinderProcessor::produce(framework::Event& event) { } else { // Fallback: sum on-track hit energies only for (const auto ts : track.trackStatesReversed()) { - if (ts.typeFlags().test(Acts::TrackStateFlag::MeasurementFlag) && + if (ts.typeFlags().isMeasurement() && ts.hasUncalibratedSourceLink()) { - const acts_examples::IndexSourceLink sl = - ts.getUncalibratedSourceLink() - .template get(); + Acts::SourceLink usl = ts.getUncalibratedSourceLink(); + const acts_examples::IndexSourceLink& sl = + usl.get(); track_energy += measurement_energies[sl.index()]; } } diff --git a/Tracking/CMakeLists.txt b/Tracking/CMakeLists.txt index 6a09400bd..3a6654d50 100644 --- a/Tracking/CMakeLists.txt +++ b/Tracking/CMakeLists.txt @@ -13,27 +13,7 @@ setup_library( ) find_package(Eigen3 CONFIG REQUIRED) -# if you want to try an alternate Acts, you can -# 1. comment out the stuff here and replace with add_subdirectory(acts) -# 2. git clone acts into this directory and choose your version -find_package(Acts QUIET) -if (NOT Acts_FOUND) - message(STATUS "Did not find Acts, downloading and compiling it here") - include(FetchContent) - FetchContent_Declare( - Acts - URL https://github.com/acts-project/acts/archive/refs/tags/v36.0.0.tar.gz - URL_HASH MD5=f543dd8ba030bea2e4ee2f1b07dbe7c0 - DOWNLOAD_EXTRACT_TIMESTAMP TRUE - ) - FetchContent_MakeAvailable(Acts) - FetchContent_GetProperties(Acts) - # Adding Acts as "SYSTEM" - # which will silence compiler warnings from these 3rd party softwares - include_directories(SYSTEM "${Acts_SOURCE_DIR}/Core/include") -else() - message(STATUS "Found Acts ${Acts_VERSION}") -endif() +# Acts is set up by the top-level CMakeLists.txt file(GLOB SRC_FILES CONFIGURE_DEPENDS ${PROJECT_SOURCE_DIR}/src/Tracking/Sim/[a-zA-z]*.cxx @@ -55,7 +35,7 @@ list(APPEND SRC_FILES setup_library(module Tracking dependencies Framework::Configure Framework::Framework - ActsCore + Acts::Core Geant4::Interface ROOT::Physics Tracking::Event diff --git a/Tracking/include/Tracking/Reco/ActsUtils.h b/Tracking/include/Tracking/Reco/ActsUtils.h index 1415448c9..99569a658 100644 --- a/Tracking/include/Tracking/Reco/ActsUtils.h +++ b/Tracking/include/Tracking/Reco/ActsUtils.h @@ -4,7 +4,7 @@ #include "Acts/Definitions/Algebra.hpp" #include "Acts/Definitions/TrackParametrization.hpp" #include "Acts/Definitions/Units.hpp" -#include "Acts/EventData/TrackParameters.hpp" +#include "Acts/EventData/BoundTrackParameters.hpp" #include "Acts/Surfaces/PerigeeSurface.hpp" namespace tracking { diff --git a/Tracking/include/Tracking/Reco/CKFProcessor.h b/Tracking/include/Tracking/Reco/CKFProcessor.h index 79b7574f9..8f8bb324b 100644 --- a/Tracking/include/Tracking/Reco/CKFProcessor.h +++ b/Tracking/include/Tracking/Reco/CKFProcessor.h @@ -20,7 +20,7 @@ #include "Acts/Definitions/Common.hpp" #include "Acts/Definitions/TrackParametrization.hpp" #include "Acts/Definitions/Units.hpp" -#include "Acts/EventData/TrackParameters.hpp" +#include "Acts/EventData/BoundTrackParameters.hpp" #include "Acts/Utilities/Logger.hpp" // geometry @@ -35,13 +35,13 @@ // propagation testing #include "Acts/MagneticField/ConstantBField.hpp" -#include "Acts/Propagator/AbortList.hpp" -#include "Acts/Propagator/ActionList.hpp" -#include "Acts/Propagator/DenseEnvironmentExtension.hpp" +#include "Acts/Propagator/ActorList.hpp" +#include "Acts/Propagator/EigenStepperDenseExtension.hpp" #include "Acts/Propagator/MaterialInteractor.hpp" #include "Acts/Propagator/Navigator.hpp" #include "Acts/Propagator/Propagator.hpp" #include "Acts/Propagator/StandardAborters.hpp" +#include "Acts/Propagator/VoidNavigator.hpp" #include "Acts/Propagator/detail/SteppingLogger.hpp" #include "Acts/Surfaces/PerigeeSurface.hpp" #include "Acts/Utilities/Logger.hpp" @@ -55,7 +55,7 @@ #include "Acts/Geometry/GeometryIdentifier.hpp" #include "Acts/TrackFinding/CombinatorialKalmanFilter.hpp" #include "Acts/TrackFinding/MeasurementSelector.hpp" -#include "Acts/TrackFitting/GainMatrixSmoother.hpp" +#include "Acts/TrackFinding/TrackStateCreator.hpp" #include "Acts/TrackFitting/GainMatrixUpdater.hpp" #include "Acts/Utilities/CalibrationContext.hpp" @@ -78,13 +78,14 @@ #include "Tracking/Sim/BFieldXYZUtils.h" // mg Aug 2024 not sure if these are needed... using Updater = Acts::GainMatrixUpdater; -using Smoother = Acts::GainMatrixSmoother; using ActionList = - Acts::ActionList; -using AbortList = Acts::AbortList; + Acts::ActorList; using CkfPropagator = Acts::Propagator, Acts::Navigator>; +using ExtrapPropagator = + Acts::Propagator, Acts::VoidNavigator>; using TrackContainer = Acts::TrackContainer; @@ -222,8 +223,10 @@ class CKFProcessor final : public TrackingGeometryUser { const Acts::CombinatorialKalmanFilter> ckf_; - // Track Extrapolator Tool - std::shared_ptr> + // Track Extrapolator Tool (uses VoidNavigator to propagate freely to any + // surface) + std::unique_ptr propagator_extrap_; + std::shared_ptr> trk_extrap_; // Zero-B CKF as fallback @@ -231,7 +234,8 @@ class CKFProcessor final : public TrackingGeometryUser { std::unique_ptr< const Acts::CombinatorialKalmanFilter> ckf_zero_b_; - std::shared_ptr> + std::unique_ptr propagator_extrap_zero_b_; + std::shared_ptr> trk_extrap_zero_b_; // Const-B (1.5T) CKF as fallback for tagger @@ -239,7 +243,8 @@ class CKFProcessor final : public TrackingGeometryUser { std::unique_ptr< const Acts::CombinatorialKalmanFilter> ckf_const_b_; - std::shared_ptr> + std::unique_ptr propagator_extrap_const_b_; + std::shared_ptr> trk_extrap_const_b_; /// n seeds and n tracks diff --git a/Tracking/include/Tracking/Reco/GSFProcessor.h b/Tracking/include/Tracking/Reco/GSFProcessor.h index e63e7223d..cd3d03252 100644 --- a/Tracking/include/Tracking/Reco/GSFProcessor.h +++ b/Tracking/include/Tracking/Reco/GSFProcessor.h @@ -19,7 +19,7 @@ #include "Acts/Definitions/Common.hpp" #include "Acts/Definitions/TrackParametrization.hpp" #include "Acts/Definitions/Units.hpp" -#include "Acts/EventData/TrackParameters.hpp" +#include "Acts/EventData/BoundTrackParameters.hpp" #include "Acts/Utilities/Logger.hpp" // geometry @@ -34,13 +34,13 @@ // propagation testing #include "Acts/MagneticField/ConstantBField.hpp" -#include "Acts/Propagator/AbortList.hpp" -#include "Acts/Propagator/ActionList.hpp" -#include "Acts/Propagator/DenseEnvironmentExtension.hpp" +#include "Acts/Propagator/ActorList.hpp" +#include "Acts/Propagator/EigenStepperDenseExtension.hpp" #include "Acts/Propagator/MaterialInteractor.hpp" #include "Acts/Propagator/Navigator.hpp" #include "Acts/Propagator/Propagator.hpp" #include "Acts/Propagator/StandardAborters.hpp" +#include "Acts/Propagator/VoidNavigator.hpp" #include "Acts/Propagator/detail/SteppingLogger.hpp" #include "Acts/Surfaces/PerigeeSurface.hpp" #include "Acts/Utilities/Logger.hpp" @@ -55,7 +55,6 @@ #include "Acts/Geometry/GeometryIdentifier.hpp" #include "Acts/TrackFinding/CombinatorialKalmanFilter.hpp" #include "Acts/TrackFinding/MeasurementSelector.hpp" -#include "Acts/TrackFitting/GainMatrixSmoother.hpp" #include "Acts/TrackFitting/GainMatrixUpdater.hpp" #include "Acts/Utilities/CalibrationContext.hpp" @@ -80,8 +79,8 @@ #include "Tracking/Sim/BFieldXYZUtils.h" using ActionList = - Acts::ActionList; -using AbortList = Acts::AbortList; + Acts::ActorList; // using GsfPropagator = Acts::Propagator< // Acts::MultiEigenStepperLoop< @@ -94,7 +93,8 @@ using AbortList = Acts::AbortList; using MultiStepper = Acts::MultiEigenStepperLoop<>; using Propagator = Acts::Propagator, Acts::Navigator>; using GsfPropagator = Acts::Propagator; -using BetheHeitlerApprox = Acts::AtlasBetheHeitlerApprox<6, 5>; +using GsfExtrapPropagator = + Acts::Propagator, Acts::VoidNavigator>; namespace tracking { namespace reco { @@ -164,8 +164,13 @@ class GSFProcessor final : public TrackingGeometryUser { // Processing time counter // double processing_time_{0.}; - /// Time profiling data for performance analysis - std::map profiling_map_; + int nevents_{0}; + int n_input_tracks_{0}; + int n_gsf_failed_{0}; + int n_output_tracks_{0}; + int n_target_extrap_failed_{0}; + int n_ecal_extrap_failed_{0}; + double processing_time_{0.}; // refitting of tracks // bool kf_refit_{false}; @@ -218,8 +223,8 @@ class GSFProcessor final : public TrackingGeometryUser { std::string seed_coll_name_{"seedTracks"}; /// Gaussian Sum Fitter instance for track refitting - std::unique_ptr> + std::unique_ptr< + const Acts::GaussianSumFitter> gsf_; /// Collection name for input tracks to be refit @@ -273,13 +278,19 @@ class GSFProcessor final : public TrackingGeometryUser { /// Layer ID to ACTS Surface mapping for hit surface lookup std::unordered_map layer_surface_map_; - // Track Extrapolator Tool - std::shared_ptr> + // Track Extrapolator Tool (VoidNavigator to reach surfaces outside geometry) + std::unique_ptr propagator_extrap_; + std::shared_ptr> trk_extrap_; - /// Beam origin surface at z=-700 mm (tagger track initialization) + /// Beam origin surface at z=-700 mm (tagger post-fit extrapolation via + /// VoidNavigator) std::shared_ptr beam_origin_surface_; + /// Tagger GSF start surface at x≈-617mm in ACTS (1mm inside tagger volume + /// outer boundary ~-618mm, 1.5mm upstream of L1 sensors at x=-615.5mm) + std::shared_ptr tagger_start_surface_; + /// Target surface at z=0 mm (recoil track initialization, perigee output) std::shared_ptr target_surface_; diff --git a/Tracking/include/Tracking/Reco/GreedyAmbiguitySolver.h b/Tracking/include/Tracking/Reco/GreedyAmbiguitySolver.h index d9264bcb3..01e531231 100644 --- a/Tracking/include/Tracking/Reco/GreedyAmbiguitySolver.h +++ b/Tracking/include/Tracking/Reco/GreedyAmbiguitySolver.h @@ -86,7 +86,14 @@ class GreedyAmbiguitySolver final : public TrackingGeometryUser { */ void produce(framework::Event& event) override; + void onProcessEnd() override; + private: + int nevents_{0}; + int n_input_tracks_{0}; + int n_output_tracks_{0}; + double processing_time_{0.}; + /// Maximum amount of shared hits per track. std::uint32_t maximum_shared_hits_{1}; /// Maximum number of iterations diff --git a/Tracking/include/Tracking/Reco/SeedFinderProcessor.h b/Tracking/include/Tracking/Reco/SeedFinderProcessor.h index 27ece431c..0bf9e53d3 100644 --- a/Tracking/include/Tracking/Reco/SeedFinderProcessor.h +++ b/Tracking/include/Tracking/Reco/SeedFinderProcessor.h @@ -22,9 +22,6 @@ #include "Acts/Definitions/Algebra.hpp" #include "Acts/MagneticField/MagneticFieldContext.hpp" #include "Acts/Seeding/EstimateTrackParamsFromSeed.hpp" -#include "Acts/Seeding/Seed.hpp" -#include "Acts/Seeding/SeedFilter.hpp" -#include "Acts/Seeding/SpacePointGrid.hpp" #include "Acts/Utilities/CalibrationContext.hpp" #include "Acts/Utilities/Intersection.hpp" @@ -87,8 +84,8 @@ class SeedFinderProcessor : public TrackingGeometryUser { const Acts::Vector3& perigee_location, const ldmx::Measurements& pmeas_tgt); - void lineParabolaToHelix(const Acts::ActsVector<5> parameters, - Acts::ActsVector<5>& helix_parameters, + void lineParabolaToHelix(const Acts::Vector<5> parameters, + Acts::Vector<5>& helix_parameters, Acts::Vector3 ref); Acts::Vector3 b_field_; diff --git a/Tracking/include/Tracking/Reco/TrackExtrapolatorTool.h b/Tracking/include/Tracking/Reco/TrackExtrapolatorTool.h index bc38a6b13..3aaa05b20 100644 --- a/Tracking/include/Tracking/Reco/TrackExtrapolatorTool.h +++ b/Tracking/include/Tracking/Reco/TrackExtrapolatorTool.h @@ -5,21 +5,22 @@ #include #include "Acts/Definitions/TrackParametrization.hpp" +#include "Acts/EventData/ParticleHypothesis.hpp" #include "Acts/EventData/TrackContainer.hpp" #include "Acts/EventData/TrackProxy.hpp" #include "Acts/Geometry/GeometryContext.hpp" #include "Acts/MagneticField/MagneticFieldContext.hpp" -#include "Acts/Propagator/AbortList.hpp" -#include "Acts/Propagator/ActionList.hpp" +#include "Acts/Propagator/ActorList.hpp" #include "Acts/Propagator/MaterialInteractor.hpp" #include "Acts/Propagator/Propagator.hpp" #include "Acts/Propagator/detail/SteppingLogger.hpp" +#include "Acts/Utilities/TrackHelpers.hpp" #include "Tracking/Event/Track.h" #include "Tracking/Sim/TrackingUtils.h" using ActionList = - Acts::ActionList; -using AbortList = Acts::AbortList; + Acts::ActorList; namespace tracking { namespace reco { @@ -31,10 +32,7 @@ class TrackExtrapolatorTool { TrackExtrapolatorTool(propagator_t propagator, const Acts::GeometryContext& gctx, const Acts::MagneticFieldContext& mctx) - : propagator_(std::move(propagator)) { - gctx_ = gctx; - mctx_ = mctx; - } + : propagator_(std::move(propagator)), gctx_(gctx), mctx_(mctx) {} /** * Turn on/off internal debug flag @@ -53,8 +51,7 @@ class TrackExtrapolatorTool { @return optional with BoundTrackParameters */ - using PropagatorOptions = - typename propagator_t::template Options; + using PropagatorOptions = typename propagator_t::template Options; std::optional extrapolate( const Acts::BoundTrackParameters pars, @@ -66,9 +63,9 @@ class TrackExtrapolatorTool { if (max_step_size_ > 0) p_options.stepping.maxStepSize = max_step_size_; if (path_limit_ > 0) p_options.pathLimit = path_limit_; - p_options.direction = intersection.intersections()[0].pathLength() >= 0 - ? Acts::Direction::Forward - : Acts::Direction::Backward; + p_options.direction = intersection[0].pathLength() >= 0 + ? Acts::Direction::Forward() + : Acts::Direction::Backward(); auto result = propagator_.propagate(pars, *target_surface, p_options); @@ -112,88 +109,54 @@ class TrackExtrapolatorTool { << target_surface.get() << std::endl; } - // get first and last track state on surface - if (debug_) - std::cout - << "[TrackExtrapolatorTool] Getting outermost track state...\n"; - auto outermost = *(track.trackStatesReversed().begin()); - if (debug_) - std::cout - << "[TrackExtrapolatorTool] Getting innermost track state...\n"; - auto begin = track.trackStatesReversed().begin(); - std::advance(begin, track.nTrackStates() - 1); - auto innermost = *begin; - if (debug_) - std::cout << "[TrackExtrapolatorTool] Got innermost and outermost " - "track states\n"; - - // I'm checking which track state is closer to the origin of the target - // surface to decide from where to start the extrapolation to the surface. I - // use the coordinate along the beam axis. - if (debug_) - std::cout << "[TrackExtrapolatorTool] Calculating distances...\n"; - double first_dis = std::abs( - innermost.referenceSurface().transform(gctx_).translation()(0) - - target_surface->transform(gctx_).translation()(0)); - - double last_dis = std::abs( - outermost.referenceSurface().transform(gctx_).translation()(0) - - target_surface->transform(gctx_).translation()(0)); - if (debug_) - std::cout << "[TrackExtrapolatorTool] first_dis = " << first_dis - << ", last_dis = " << last_dis << std::endl; - - // This is the track state to use for the extrapolation + if (track.nTrackStates() == 0) { + return std::nullopt; + } - const auto& ts = first_dis < last_dis ? innermost : outermost; - if (debug_) std::cout << "[TrackExtrapolatorTool] Selected track state\n"; + // Use ACTS's built-in helper to find the measurement track state + // (first or last) that is closest to the target surface. This correctly + // handles holes and material-only states which lack filtered parameters. + auto state_result = Acts::findTrackStateForExtrapolation( + gctx_, track, *target_surface, + Acts::TrackExtrapolationStrategy::firstOrLast); - // Get the BoundTrackStateParameters + if (!state_result.ok()) { + return std::nullopt; + } - if (debug_) - std::cout << "[TrackExtrapolatorTool] Getting reference surface...\n"; + const auto& ts = state_result->first; const auto& surface = ts.referenceSurface(); - if (debug_) - std::cout << "[TrackExtrapolatorTool] Checking hasSmoothed...\n"; - bool has_smoothed = ts.hasSmoothed(); - if (debug_) - std::cout << "[TrackExtrapolatorTool] has_smoothed = " << has_smoothed - << std::endl; - // Use smoothed parameters if available, otherwise use filtered Acts::BoundVector params; Acts::BoundMatrix cov; - if (has_smoothed) { + if (ts.hasSmoothed()) { if (debug_) - std::cout << "[TrackExtrapolatorTool] Using smoothed parameters...\n"; + std::cout << "[TrackExtrapolatorTool] Using smoothed parameters\n"; params = ts.smoothed(); cov = ts.smoothedCovariance(); - } else { + } else if (ts.hasFiltered()) { if (debug_) - std::cout << "[TrackExtrapolatorTool] Using filtered parameters...\n"; + std::cout << "[TrackExtrapolatorTool] Using filtered parameters\n"; params = ts.filtered(); cov = ts.filteredCovariance(); + } else { + return std::nullopt; } - if (debug_) - std::cout << "[TrackExtrapolatorTool] Got all track state components\n"; if (debug_) { - std::cout << "Surface::" << surface.transform(gctx_).translation() + std::cout << "Surface::" + << surface.localToGlobalTransform(gctx_).translation() << std::endl; - std::cout << "HasSmoothed::" << has_smoothed << std::endl; + std::cout << "HasSmoothed::" << ts.hasSmoothed() << std::endl; std::cout << "Parameters::" << params.transpose() << std::endl; } - // mg Aug 2024 ... v36 takes the particle...assume electron - if (debug_) - std::cout - << "[TrackExtrapolatorTool] Creating BoundTrackParameters...\n"; - auto part_hypo{Acts::SinglyChargedParticleHypothesis::electron()}; + + auto part_hypo{Acts::ParticleHypothesis::electron()}; Acts::BoundTrackParameters sp(surface.getSharedPtr(), params, cov, part_hypo); if (debug_) - std::cout << "[TrackExtrapolatorTool] BoundTrackParameters created, " - "calling extrapolate(BTP)...\n"; + std::cout << "[TrackExtrapolatorTool] calling extrapolate(BTP)...\n"; auto result = extrapolate(sp, target_surface); if (debug_) std::cout << "[TrackExtrapolatorTool] extrapolate DONE\n"; return result; @@ -224,7 +187,7 @@ class TrackExtrapolatorTool { // Get the BoundTrackStateParameters // assume electron for now - auto part_hypo{Acts::SinglyChargedParticleHypothesis::electron()}; + auto part_hypo{Acts::ParticleHypothesis::electron()}; Acts::BoundTrackParameters state_parameters(surface.getSharedPtr(), smoothed, cov, part_hypo); @@ -276,7 +239,8 @@ class TrackExtrapolatorTool { if (opt_pars) { if (debug_) { - Acts::Vector3 surf_loc = target_surface->transform(gctx_).translation(); + Acts::Vector3 surf_loc = + target_surface->localToGlobalTransform(gctx_).translation(); std::cout << "[TrackExtrapolatorTool] Surface location: (" << surf_loc(0) << ", " << surf_loc(1) << ", " << surf_loc(2) << ")\n"; diff --git a/Tracking/include/Tracking/Reco/TruthSeedProcessor.h b/Tracking/include/Tracking/Reco/TruthSeedProcessor.h index 07b25dcc5..92b514651 100644 --- a/Tracking/include/Tracking/Reco/TruthSeedProcessor.h +++ b/Tracking/include/Tracking/Reco/TruthSeedProcessor.h @@ -19,7 +19,7 @@ #include "Acts/Definitions/Algebra.hpp" #include "Acts/Definitions/TrackParametrization.hpp" -#include "Acts/EventData/TrackParameters.hpp" +#include "Acts/EventData/BoundTrackParameters.hpp" #include "Acts/Propagator/Navigator.hpp" #include "Acts/Propagator/Propagator.hpp" #include "Acts/Surfaces/PerigeeSurface.hpp" @@ -189,9 +189,6 @@ class TruthSeedProcessor : public TrackingGeometryUser { const std::shared_ptr& origin_surface, const std::shared_ptr& target_surface); - /// The ACTS geometry context properly - Acts::GeometryContext gctx_; - /// pdg_ids of the particles we want to select for the seeds std::vector pdg_ids_{11}; diff --git a/Tracking/include/Tracking/Reco/VertexProcessor.h b/Tracking/include/Tracking/Reco/VertexProcessor.h index 217b11d1e..d59ce0171 100644 --- a/Tracking/include/Tracking/Reco/VertexProcessor.h +++ b/Tracking/include/Tracking/Reco/VertexProcessor.h @@ -39,7 +39,9 @@ #include "TLorentzVector.h" // Propagator with void navigator -using VoidPropagator = Acts::Propagator>; +#include "Acts/Propagator/VoidNavigator.hpp" +using VoidPropagator = + Acts::Propagator, Acts::VoidNavigator>; namespace tracking { namespace reco { @@ -77,7 +79,6 @@ class VertexProcessor : public framework::Producer { private: /// The contexts - TODO: they should move to some global location, I guess - Acts::GeometryContext gctx_; Acts::MagneticFieldContext bctx_; // Event counter diff --git a/Tracking/include/Tracking/Reco/Vertexer.h b/Tracking/include/Tracking/Reco/Vertexer.h index 0a058222e..8ca4eceeb 100644 --- a/Tracking/include/Tracking/Reco/Vertexer.h +++ b/Tracking/include/Tracking/Reco/Vertexer.h @@ -40,7 +40,9 @@ #include "Acts/Surfaces/PerigeeSurface.hpp" // Propagator with void navigator -using VoidPropagator = Acts::Propagator>; +#include "Acts/Propagator/VoidNavigator.hpp" +using VoidPropagator = + Acts::Propagator, Acts::VoidNavigator>; namespace tracking { namespace reco { @@ -62,7 +64,6 @@ class Vertexer : public framework::Producer { const std::vector& recoil_tracks); private: - Acts::GeometryContext gctx_; Acts::MagneticFieldContext bctx_; int nevents_{0}; diff --git a/Tracking/include/Tracking/Sim/BFieldXYZUtils.h b/Tracking/include/Tracking/Sim/BFieldXYZUtils.h index 7c526677a..f3d8eeda5 100644 --- a/Tracking/include/Tracking/Sim/BFieldXYZUtils.h +++ b/Tracking/include/Tracking/Sim/BFieldXYZUtils.h @@ -8,7 +8,7 @@ #include "Acts/MagneticField/BFieldMapUtils.hpp" #include "Acts/MagneticField/InterpolatedBFieldMap.hpp" #include "Acts/MagneticField/MagneticFieldContext.hpp" -#include "Acts/Utilities/AxisFwd.hpp" +#include "Acts/Utilities/AxisDefinitions.hpp" #include "Acts/Utilities/Grid.hpp" #include "Acts/Utilities/Interpolation.hpp" #include "Acts/Utilities/Result.hpp" @@ -211,8 +211,8 @@ inline InterpolatedMagneticField3 makeMagneticFieldMapXyzFromText( localToGlobalBin, GenericTransformPos transformPosition, GenericTransformBField transformMagneticField, - const std::string& fieldMapFile, Acts::ActsScalar lengthUnit, - Acts::ActsScalar BFieldUnit, bool firstOctant, bool rotateAxes) { + const std::string& fieldMapFile, double lengthUnit, double BFieldUnit, + bool firstOctant, bool rotateAxes) { /// [1] Read in field map file // Grid position points in x, y and z std::vector x_pos; diff --git a/Tracking/include/Tracking/Sim/GeometryContainers.h b/Tracking/include/Tracking/Sim/GeometryContainers.h index b013a542f..7feddc9a0 100644 --- a/Tracking/include/Tracking/Sim/GeometryContainers.h +++ b/Tracking/include/Tracking/Sim/GeometryContainers.h @@ -105,11 +105,11 @@ template inline acts_examples::Range::const_iterator> selectVolume(const GeometryIdMultiset& container, Acts::GeometryIdentifier::Value volume) { - auto cmp = Acts::GeometryIdentifier().setVolume(volume); + auto cmp = Acts::GeometryIdentifier().withVolume(volume); auto beg = std::lower_bound(container.begin(), container.end(), cmp, detail::CompareGeometryId{}); // WARNING overflows to volume==0 if the input volume is the last one - cmp = Acts::GeometryIdentifier().setVolume(volume + 1u); + cmp = Acts::GeometryIdentifier().withVolume(volume + 1u); // optimize search by using the lower bound as start point. also handles // volume overflows since the geo id would be located before the start of // the upper edge search window. @@ -131,11 +131,11 @@ inline acts_examples::Range::const_iterator> selectLayer(const GeometryIdMultiset& container, Acts::GeometryIdentifier::Value volume, Acts::GeometryIdentifier::Value layer) { - auto cmp = Acts::GeometryIdentifier().setVolume(volume).setLayer(layer); + auto cmp = Acts::GeometryIdentifier().withVolume(volume).withLayer(layer); auto beg = std::lower_bound(container.begin(), container.end(), cmp, detail::CompareGeometryId{}); // WARNING resets to layer==0 if the input layer is the last one - cmp = Acts::GeometryIdentifier().setVolume(volume).setLayer(layer + 1u); + cmp = Acts::GeometryIdentifier().withVolume(volume).withLayer(layer + 1u); // optimize search by using the lower bound as start point. also handles // volume overflows since the geo id would be located before the start of // the upper edge search window. @@ -166,10 +166,10 @@ inline auto selectModule(const GeometryIdMultiset& container, Acts::GeometryIdentifier::Value volume, Acts::GeometryIdentifier::Value layer, Acts::GeometryIdentifier::Value module_) { - return selectModule( - container, - Acts::GeometryIdentifier().setVolume(volume).setLayer(layer).setSensitive( - module_)); + return selectModule(container, Acts::GeometryIdentifier() + .withVolume(volume) + .withLayer(layer) + .withSensitive(module_)); } /// Select all elements for the lowest non-zero identifier component. diff --git a/Tracking/include/Tracking/Sim/LdmxSpacePoint.h b/Tracking/include/Tracking/Sim/LdmxSpacePoint.h index 2b0ec252b..d68c2b317 100644 --- a/Tracking/include/Tracking/Sim/LdmxSpacePoint.h +++ b/Tracking/include/Tracking/Sim/LdmxSpacePoint.h @@ -114,7 +114,7 @@ class LdmxSpacePoint { // Projection matrix from the full space to the (u,v) space. // This can be expanded to (u,v,t) space in the case time needs to be added. - Acts::ActsMatrix<2, 6> projector_; + Acts::Matrix<2, 6> projector_; private: void initialize() { diff --git a/Tracking/include/Tracking/Sim/MeasurementCalibrator.h b/Tracking/include/Tracking/Sim/MeasurementCalibrator.h index 1c090d2ed..75ad1bc58 100644 --- a/Tracking/include/Tracking/Sim/MeasurementCalibrator.h +++ b/Tracking/include/Tracking/Sim/MeasurementCalibrator.h @@ -4,6 +4,7 @@ #include #include "Acts/Definitions/Algebra.hpp" +#include "Acts/Definitions/TrackParametrization.hpp" #include "Acts/EventData/MultiTrajectory.hpp" #include "Acts/EventData/SourceLink.hpp" #include "Acts/EventData/VectorMultiTrajectory.hpp" @@ -84,13 +85,9 @@ class LdmxMeasurementCalibrator { // make tsCalCov 2x2 block the local_cov we just set ts_cal_cov.block(0, 0, 2, 2) = local_cov; - Acts::ActsMatrix<2, 6> projector; - projector.setZero(); - projector(0, 0) = 1.; - projector(1, 1) = 1.; - - trackState.setProjector(projector); - trackState.setUncalibratedSourceLink(genericSourceLink); + trackState.setProjectorSubspaceIndices( + std::array{Acts::eBoundLoc0, Acts::eBoundLoc1}); + trackState.setUncalibratedSourceLink(Acts::SourceLink{genericSourceLink}); } /// Find the measurement corresponding to the source link. @@ -124,12 +121,9 @@ class LdmxMeasurementCalibrator { ts_cal_cov.setZero(); ts_cal_cov(0, 0) = (meas.getLocalCovariance())[0]; - Acts::ActsMatrix<2, 6> projector; - projector.setZero(); - projector(0, 0) = 1.; - projector(1, 1) = 1.; - trackState.setProjector(projector.row(0)); - trackState.setUncalibratedSourceLink(genericSourceLink); + trackState.setProjectorSubspaceIndices( + std::array{Acts::eBoundLoc0}); + trackState.setUncalibratedSourceLink(Acts::SourceLink{genericSourceLink}); } // Function to test the measurement calibrator diff --git a/Tracking/include/Tracking/Sim/SeedToTrackParamMaker.h b/Tracking/include/Tracking/Sim/SeedToTrackParamMaker.h index 9103ca24d..8560f2c14 100644 --- a/Tracking/include/Tracking/Sim/SeedToTrackParamMaker.h +++ b/Tracking/include/Tracking/Sim/SeedToTrackParamMaker.h @@ -2,7 +2,6 @@ #include "Acts/Definitions/Algebra.hpp" #include "Acts/Definitions/TrackParametrization.hpp" -#include "Acts/Seeding/Seed.hpp" // #include "Acts/Utilities/VectorHelpers.hpp" #include @@ -44,10 +43,11 @@ class SeedToTrackParamMaker { /// This resembles the method used in ATLAS for the seed fitting /// L811 /// https://acode-browser.usatlas.bnl.gov/lxr/source/athena/InnerDetector/InDetRecTools/SiTrackMakerTool_xk/src/SiTrackMaker_xk.cxx - template - bool fitSeedAtlas(const Acts::Seed& seed, - std::array& data, const Acts::Transform3& Tp, - const double& bFieldZ); + // Acts::Seed removed in v47 — fitSeedAtlas(Seed) unused, commented out + // template + // bool fitSeedAtlas(const Acts::Seed& seed, + // std::array& data, const Acts::Transform3& Tp, + // const double& bFieldZ); template bool fitSeedAtlas(const std::vector& sp, @@ -56,9 +56,10 @@ class SeedToTrackParamMaker { /// This is a simple Line and Parabola fit (from HPS reconstruction by Robert /// Johnson) - template - bool fitSeedLinPar(const Acts::Seed& seed, - std::vector& data); + // Acts::Seed removed in v47 — fitSeedLinPar(Seed) unused, commented out + // template + // bool fitSeedLinPar(const Acts::Seed& seed, + // std::vector& data); /// Estimate the full track parameters from three space points /// @@ -91,9 +92,8 @@ class SeedToTrackParamMaker { template std::optional estimateTrackParamsFromSeed( const Acts::Transform3& Tp, spacepoint_iterator_t spBegin, - spacepoint_iterator_t spEnd, Acts::Vector3 bField, - Acts::ActsScalar bFieldMin, - Acts::ActsScalar mass = 139.57018 * Acts::UnitConstants::MeV) { + spacepoint_iterator_t spEnd, Acts::Vector3 bField, double bFieldMin, + double mass = 139.57018 * Acts::UnitConstants::MeV) { // Check the number of provided space points size_t num_sp = std::distance(spBegin, spEnd); if (num_sp != 3) { @@ -102,8 +102,8 @@ class SeedToTrackParamMaker { } // Convert bField to Tesla - Acts::ActsScalar b_field_in_tesla = bField.norm() / Acts::UnitConstants::T; - Acts::ActsScalar b_field_min_in_tesla = bFieldMin / Acts::UnitConstants::T; + double b_field_in_tesla = bField.norm() / Acts::UnitConstants::T; + double b_field_min_in_tesla = bFieldMin / Acts::UnitConstants::T; // Check if magnetic field is too small if (b_field_in_tesla < b_field_min_in_tesla) { // @todo shall we use straight-line estimation and use default q/pt in @@ -156,8 +156,7 @@ class SeedToTrackParamMaker { // Lambda to transform the coordinates to the (u, v) space auto uv_transform = [](const Acts::Vector3& local) -> Acts::Vector2 { Acts::Vector2 uv; - Acts::ActsScalar denominator = - local.x() * local.x() + local.y() * local.y(); + double denominator = local.x() * local.x() + local.y() * local.y(); uv.x() = local.x() / denominator; uv.y() = local.y() / denominator; return uv; @@ -168,15 +167,15 @@ class SeedToTrackParamMaker { // A,B are slope and intercept of the straight line in the u,v plane // connecting the three points - Acts::ActsScalar a = (uv2.y() - uv1.y()) / (uv2.x() - uv1.x()); - Acts::ActsScalar b = uv2.y() - a * uv2.x(); + double a = (uv2.y() - uv1.y()) / (uv2.x() - uv1.x()); + double b = uv2.y() - a * uv2.x(); // Curvature (with a sign) estimate - Acts::ActsScalar rho = -2.0 * b / std::hypot(1., a); + double rho = -2.0 * b / std::hypot(1., a); // The projection of the top space point on the transverse plane of the new // frame - Acts::ActsScalar rn = local2.x() * local2.x() + local2.y() * local2.y(); + double rn = local2.x() * local2.x() + local2.y() * local2.y(); // The (1/tanTheta) of momentum in the new frame, - Acts::ActsScalar inv_tan_theta = + double inv_tan_theta = local2.z() * std::sqrt(1. / rn) / (1. + rho * rho * rn); // The momentum direction in the new frame (the center of the circle has the // coordinate (-1.*A/(2*B), 1./(2*B))) @@ -200,23 +199,23 @@ class SeedToTrackParamMaker { // The estimated q/pt in [GeV/c]^-1 (note that the pt is the projection of // momentum on the transverse plane of the new frame) - Acts::ActsScalar q_over_pt = + double q_over_pt = rho * (Acts::UnitConstants::m) / (0.3 * b_field_in_tesla); // The estimated q/p in [GeV/c]^-1 params[Acts::eBoundQOverP] = q_over_pt / std::hypot(1., inv_tan_theta); // The estimated momentum, and its projection along the magnetic field // diretion - Acts::ActsScalar p_in_ge_v = std::abs(1.0 / params[Acts::eBoundQOverP]); - Acts::ActsScalar pz_in_ge_v = 1.0 / std::abs(q_over_pt) * inv_tan_theta; - Acts::ActsScalar mass_in_ge_v = mass / Acts::UnitConstants::GeV; + double p_in_ge_v = std::abs(1.0 / params[Acts::eBoundQOverP]); + double pz_in_ge_v = 1.0 / std::abs(q_over_pt) * inv_tan_theta; + double mass_in_ge_v = mass / Acts::UnitConstants::GeV; // The estimated velocity, and its projection along the magnetic field // diretion - Acts::ActsScalar v = p_in_ge_v / std::hypot(p_in_ge_v, mass_in_ge_v); - Acts::ActsScalar vz = pz_in_ge_v / std::hypot(p_in_ge_v, mass_in_ge_v); + double v = p_in_ge_v / std::hypot(p_in_ge_v, mass_in_ge_v); + double vz = pz_in_ge_v / std::hypot(p_in_ge_v, mass_in_ge_v); // The z_ coordinate of the bottom space point along the magnetic field // direction - Acts::ActsScalar pathz = sp_global_positions[0].dot(bField) / bField.norm(); + double pathz = sp_global_positions[0].dot(bField) / bField.norm(); // The estimated time (use path length along magnetic field only if it's not // zero) if (pathz != 0) { diff --git a/Tracking/include/Tracking/Sim/SeedToTrackParamMaker.ipp b/Tracking/include/Tracking/Sim/SeedToTrackParamMaker.ipp index b48bcb2f0..72fec4ae9 100644 --- a/Tracking/include/Tracking/Sim/SeedToTrackParamMaker.ipp +++ b/Tracking/include/Tracking/Sim/SeedToTrackParamMaker.ipp @@ -175,12 +175,13 @@ bool SeedToTrackParamMaker::karimakiFit( /// see /// https://acode-browser.usatlas.bnl.gov/lxr/source/athena/InnerDetector/InDetRecTools/SiTrackMakerTool_xk/src/SiTrackMaker_xk.cxx -template -bool SeedToTrackParamMaker::fitSeedAtlas( - const Acts::Seed& seed, std::array& data, - const Acts::Transform3& Tp, const double& bFieldZ) { - return FitSeedAtlas(seed.sp(), data, Tp, bFieldZ); -} +// Acts::Seed removed in v47 — fitSeedAtlas(Seed) unused, commented out +// template +// bool SeedToTrackParamMaker::fitSeedAtlas( +// const Acts::Seed& seed, std::array& data, +// const Acts::Transform3& Tp, const double& bFieldZ) { +// return FitSeedAtlas(seed.sp(), data, Tp, bFieldZ); +// } // double H = 0.0015; //kTesla @@ -272,11 +273,12 @@ bool SeedToTrackParamMaker::fitSeedAtlas( return true; } -template -bool SeedToTrackParamMaker::fitSeedLinPar( - const Acts::Seed& seed, std::vector& data) { - return true; -} +// Acts::Seed removed in v47 — fitSeedLinPar(Seed) unused, commented out +// template +// bool SeedToTrackParamMaker::fitSeedLinPar( +// const Acts::Seed& seed, std::vector& data) { +// return true; +// } } // namespace sim } // namespace tracking diff --git a/Tracking/include/Tracking/Sim/TrackingUtils.h b/Tracking/include/Tracking/Sim/TrackingUtils.h index 139d98ff0..c7819f5b3 100644 --- a/Tracking/include/Tracking/Sim/TrackingUtils.h +++ b/Tracking/include/Tracking/Sim/TrackingUtils.h @@ -34,7 +34,7 @@ #include "Acts/Definitions/PdgParticle.hpp" #include "Acts/Definitions/TrackParametrization.hpp" #include "Acts/Definitions/Units.hpp" -#include "Acts/EventData/TrackParameters.hpp" +#include "Acts/EventData/BoundTrackParameters.hpp" #include "Acts/Surfaces/PerigeeSurface.hpp" #include "Acts/Surfaces/PlaneSurface.hpp" #include "Acts/Surfaces/Surface.hpp" @@ -61,9 +61,9 @@ ldmx::LdmxSpacePoint* convertSimHitToLdmxSpacePoint( // BoundSymMatrix doesn't exist in v36 .. use BoundSquareMatrix // have to change this everywhere .. I think using BoundSysMatrix was defined // exactly the same as BoundSquareMatrix is now in ACTs -void flatCov(Acts::BoundSquareMatrix cov, std::vector& v_cov); +void flatCov(Acts::BoundMatrix cov, std::vector& v_cov); -Acts::BoundSquareMatrix unpackCov(const std::vector& v_cov); +Acts::BoundMatrix unpackCov(const std::vector& v_cov); // Rotate LDMX global -> ACTS frame: z_ldmx->x_acts, x_ldmx->y_acts, // y_ldmx->z_acts (0 0 1) * (x,y,z)_ldmx = x_acts (1 0 0) * (x,y,z)_ldmx = @@ -80,7 +80,7 @@ Acts::Vector3 acts2Ldmx(Acts::Vector3 acts_v); // Transform position, momentum and charge to free parameters Acts::FreeVector toFreeParameters(Acts::Vector3 pos_, Acts::Vector3 mom, - Acts::ActsScalar q); + double q); // Pack the acts track parameters into something that is serializable for the // event bus diff --git a/Tracking/include/Tracking/geo/DetectorElement.h b/Tracking/include/Tracking/geo/DetectorElement.h index 93a7ae566..bbd3eecb0 100644 --- a/Tracking/include/Tracking/geo/DetectorElement.h +++ b/Tracking/include/Tracking/geo/DetectorElement.h @@ -3,11 +3,11 @@ #include #include "Acts/Definitions/Algebra.hpp" -#include "Acts/Geometry/DetectorElementBase.hpp" #include "Acts/Geometry/GeometryContext.hpp" #include "Acts/Geometry/GeometryIdentifier.hpp" #include "Acts/Material/HomogeneousSurfaceMaterial.hpp" #include "Acts/Surfaces/Surface.hpp" +#include "Acts/Surfaces/SurfacePlacementBase.hpp" #include "Framework/Exception/Exception.h" #include "Tracking/geo/GeoUtils.h" @@ -16,7 +16,7 @@ namespace tracking::geo { -class DetectorElement : public Acts::DetectorElementBase { +class DetectorElement : public Acts::SurfacePlacementBase { public: // The detector element is initialized with the surface initial transformation // created from the TrackingGeometry constructor/parser @@ -54,16 +54,18 @@ class DetectorElement : public Acts::DetectorElementBase { // corrections Could be interesting to cache the transformations and re-update // all of them when IoV changes - const Acts::Transform3& transform( + const Acts::Transform3& localToGlobalTransform( const Acts::GeometryContext& gctx) const override; const Acts::Surface& surface() const override; Acts::Surface& surface() override; + bool isSensitive() const override { return true; } + // The thickness of the detector element is taken from the center of the // associated surface - double thickness() const override; + double thickness() const; Acts::GeometryIdentifier geometryId() const { if (!m_surface_) diff --git a/Tracking/include/Tracking/geo/TrackingGeometry.h b/Tracking/include/Tracking/geo/TrackingGeometry.h index 48e5ca78b..1beb3ced7 100644 --- a/Tracking/include/Tracking/geo/TrackingGeometry.h +++ b/Tracking/include/Tracking/geo/TrackingGeometry.h @@ -6,6 +6,7 @@ #include #include #include +#include #include "Acts/Definitions/Units.hpp" #include "Acts/Material/HomogeneousSurfaceMaterial.hpp" diff --git a/Tracking/python/full_tracking_sequence.py b/Tracking/python/full_tracking_sequence.py index 132fde89e..6d82d72c1 100644 --- a/Tracking/python/full_tracking_sequence.py +++ b/Tracking/python/full_tracking_sequence.py @@ -164,7 +164,7 @@ def tagged(name): depletion_voltage=70.0, noise_electrons=1000.0, threshold_electrons=3000.0, - out_raw_collection=tagged("TaggerSimHits"), + out_raw_collection=tagged("TaggerSimSiStripHits"), ) digi_recoil = tracking.DigitizationProcessor( @@ -178,7 +178,7 @@ def tagged(name): depletion_voltage=70.0, noise_electrons=1000.0, threshold_electrons=3000.0, - out_raw_collection=tagged("RecoilSimHits"), + out_raw_collection=tagged("RecoilSimSiStripHits"), ) fit_tagger = tracking.StripFitProcessor( @@ -238,6 +238,9 @@ def tagged(name): instance_name=tagged("SeedTagger"), input_hits_collection=tagger_meas_collection, out_seed_collection=tagged("TaggerRecoSeeds"), + # Perigee upstream of all tagger sensors (ACTS x from -12.5 to -615.5 mm). + # World boundary is at ACTS x = -650 mm. + perigee_location=[-617.0, 0.0, 0.0], pmin=0.03, pmax=63.0, d0min=-36.9, diff --git a/Tracking/src/Tracking/Reco/CKFProcessor.cxx b/Tracking/src/Tracking/Reco/CKFProcessor.cxx index 8c9c29499..14bde6761 100644 --- a/Tracking/src/Tracking/Reco/CKFProcessor.cxx +++ b/Tracking/src/Tracking/Reco/CKFProcessor.cxx @@ -6,6 +6,7 @@ #include "Tracking/Event/Track.h" #include "Tracking/Reco/TruthMatchingTool.h" #include "Tracking/Sim/GeometryContainers.h" +#include "Tracking/geo/DetectorElement.h" //--- C++ StdLib ---// #include //std::vector reverse @@ -91,8 +92,13 @@ void CKFProcessor::onNewRun(const ldmx::RunHeader& rh) { // Setup the finder / fitters ckf_ = std::make_unique>( *propagator_, Acts::getDefaultLogger("CKF", acts_logging_level)); + // Extrapolation uses VoidNavigator so it can reach surfaces outside the + // tracking geometry (e.g. ECAL scoring plane) without being stopped at + // volume boundaries. + propagator_extrap_ = std::make_unique( + Acts::EigenStepper<>{map}, Acts::VoidNavigator{}); trk_extrap_ = std::make_shared>( - *propagator_, geometryContext(), magneticFieldContext()); + *propagator_extrap_, geometryContext(), magneticFieldContext()); // Setup zero-B CKF as fallback Acts::ConstantBField zero_b_field(Acts::Vector3(0., 0., 0.)); @@ -103,9 +109,14 @@ void CKFProcessor::onNewRun(const ldmx::RunHeader& rh) { ckf_zero_b_ = std::make_unique>( *propagator_zero_b_, Acts::getDefaultLogger("CKF_ZERO_B", acts_logging_level)); + propagator_extrap_zero_b_ = std::make_unique( + Acts::EigenStepper<>{ + std::make_shared(zero_b_field)}, + Acts::VoidNavigator{}); trk_extrap_zero_b_ = std::make_shared>( - *propagator_zero_b_, geometryContext(), magneticFieldContext()); + *propagator_extrap_zero_b_, geometryContext(), + magneticFieldContext()); // Setup const-B (1.5T) CKF as fallback for tagger propagator_const_b_ = @@ -113,9 +124,12 @@ void CKFProcessor::onNewRun(const ldmx::RunHeader& rh) { ckf_const_b_ = std::make_unique>( *propagator_const_b_, Acts::getDefaultLogger("CKF_CONST_B", acts_logging_level)); + propagator_extrap_const_b_ = std::make_unique( + Acts::EigenStepper<>{const_b_field}, Acts::VoidNavigator{}); trk_extrap_const_b_ = std::make_shared>( - *propagator_const_b_, geometryContext(), magneticFieldContext()); + *propagator_extrap_const_b_, geometryContext(), + magneticFieldContext()); } // end of CKFProcessor::onNewRun() void CKFProcessor::produce(framework::Event& event) { @@ -136,7 +150,7 @@ void CKFProcessor::produce(framework::Event& event) { // Move this at the start of the producer Acts::PropagatorOptions + Acts::NavigatorPlainOptions, ActionList> propagator_options(geometryContext(), magneticFieldContext()); propagator_options.pathLimit = std::numeric_limits::max(); @@ -146,14 +160,14 @@ void CKFProcessor::produce(framework::Event& event) { // Switch the material interaction on/off & eventually into logging mode auto& m_interactor = - propagator_options.actionList.get(); + propagator_options.actorList.get(); m_interactor.multipleScattering = true; m_interactor.energyLoss = true; m_interactor.recordInteractions = false; // The logger can be switched to sterile, e.g. for timing logging auto& s_logger = - propagator_options.actionList.get(); + propagator_options.actorList.get(); s_logger.sterile = true; // Set a maximum step size propagator_options.stepping.maxStepSize = @@ -230,7 +244,7 @@ void CKFProcessor::produce(framework::Event& event) { param_vec << seed.getD0(), seed.getZ0(), seed.getPhi(), seed.getTheta(), seed.getQoP(), seed.getT(); - Acts::BoundSquareMatrix cov_mat = + Acts::BoundMatrix cov_mat = tracking::sim::utils::unpackCov(seed.getPerigeeCov()); ldmx_log(debug) << " For seed index_ = " << seed_track_index @@ -246,7 +260,7 @@ void CKFProcessor::produce(framework::Event& event) { << cov_mat(1, 1) << ", " << cov_mat(2, 2) << ")"; // need to set particle hypothesis...set to electron for now... - auto part_hypo{Acts::SinglyChargedParticleHypothesis::electron()}; + auto part_hypo{Acts::ParticleHypothesis::electron()}; start_parameters.push_back(Acts::BoundTrackParameters( perigee_surface, param_vec, cov_mat, part_hypo)); @@ -274,29 +288,7 @@ void CKFProcessor::produce(framework::Event& event) { tracking::sim::LdmxMeasurementCalibrator calibrator{measurements}; - Acts::CombinatorialKalmanFilterExtensions ckf_extensions; - - if (use1_dmeasurements_) { - ckf_extensions.calibrator - .connect<&tracking::sim::LdmxMeasurementCalibrator::calibrate1d< - Acts::VectorMultiTrajectory>>(&calibrator); - } else { - ckf_extensions.calibrator - .connect<&tracking::sim::LdmxMeasurementCalibrator::calibrate< - Acts::VectorMultiTrajectory>>(&calibrator); - } - - ckf_extensions.updater.connect< - &Acts::GainMatrixUpdater::operator()>( - &kf_updater); - - ckf_extensions.measurementSelector - .connect<&Acts::MeasurementSelector::select>( - &meas_sel); - - ldmx_log(debug) << "SourceLinkAccessor..."; - - // Create source link accessor and connect delegate + // Create source link accessor iterator type and lambda struct SourceLinkAccIt { using BaseIt = decltype(geo_id_sl_map.begin()); BaseIt it_; @@ -306,7 +298,6 @@ void CKFProcessor::produce(framework::Event& event) { using difference_type = typename BaseIt::difference_type; using iterator_category = typename BaseIt::iterator_category; - // using value_type = typename BaseIt::value_type::second_type; using value_type = Acts::SourceLink; using pointer = typename BaseIt::pointer; using reference = value_type&; @@ -322,9 +313,6 @@ void CKFProcessor::produce(framework::Event& event) { bool operator!=(const SourceLinkAccIt& other) const { return !(*this == other); } - // const value_type& operator*() const { return it->second; } - - // by value value_type operator*() const { return value_type{it_->second}; } }; @@ -334,11 +322,31 @@ void CKFProcessor::produce(framework::Event& event) { return {SourceLinkAccIt{begin}, SourceLinkAccIt{end}}; }; - Acts::SourceLinkAccessorDelegate - source_link_accessor_delegate; - source_link_accessor_delegate + // v46: calibrator and measurementSelector moved to TrackStateCreator + Acts::TrackStateCreator track_state_creator; + track_state_creator.sourceLinkAccessor .connect<&decltype(source_link_accessor)::operator(), decltype(source_link_accessor)>(&source_link_accessor); + if (use1_dmeasurements_) { + track_state_creator.calibrator + .connect<&tracking::sim::LdmxMeasurementCalibrator::calibrate1d< + Acts::VectorMultiTrajectory>>(&calibrator); + } else { + track_state_creator.calibrator + .connect<&tracking::sim::LdmxMeasurementCalibrator::calibrate< + Acts::VectorMultiTrajectory>>(&calibrator); + } + track_state_creator.measurementSelector + .connect<&Acts::MeasurementSelector::select>( + &meas_sel); + + Acts::CombinatorialKalmanFilterExtensions ckf_extensions; + ckf_extensions.updater.connect< + &Acts::GainMatrixUpdater::operator()>( + &kf_updater); + ckf_extensions.createTrackStates.connect<&Acts::TrackStateCreator< + SourceLinkAccIt, TrackContainer>::createTrackStates>( + &track_state_creator); ldmx_log(debug) << "Setting up surfaces..."; @@ -364,14 +372,12 @@ void CKFProcessor::produce(framework::Event& event) { ldmx_log(debug) << "---------------------------"; ldmx_log(debug) << "Candidate Track ID = " << track_id; // Define the CKF options here: - const Acts::CombinatorialKalmanFilterOptions - ckf_options(TrackingGeometryUser::geometryContext(), - TrackingGeometryUser::magneticFieldContext(), - TrackingGeometryUser::calibrationContext(), - source_link_accessor_delegate, ckf_extensions, - propagator_options, true /* multiple scattering */, - false /* energy loss */); + const Acts::CombinatorialKalmanFilterOptions ckf_options( + TrackingGeometryUser::geometryContext(), + TrackingGeometryUser::magneticFieldContext(), + TrackingGeometryUser::calibrationContext(), ckf_extensions, + static_cast(propagator_options), + true /* multiple scattering */, false /* energy loss */); ldmx_log(debug) << " Checking options: multiple scattering = " << ckf_options.multipleScattering @@ -433,7 +439,11 @@ void CKFProcessor::produce(framework::Event& event) { // For now it seems this loop is only looping on a single element for (auto& track : tracks_from_seed) { // do the track smoothing...this is not done in the CKF code anymore - Acts::smoothTrack(geometryContext(), track); // from TrackHelpers + auto smooth_result = Acts::smoothTrack(geometryContext(), track); + if (!smooth_result.ok()) { + ldmx_log(warn) << "smoothTrack failed: " + << smooth_result.error().message(); + } // Build the output Track ldmx::Track trk; @@ -505,7 +515,8 @@ void CKFProcessor::produce(framework::Event& event) { } // Perigee location: target surface origin rotated to LDMX frame Acts::Vector3 target_loc_ldmx = tracking::sim::utils::acts2Ldmx( - target_surface_->transform(geometryContext()).translation()); + target_surface_->localToGlobalTransform(geometryContext()) + .translation()); trk.setPerigeeLocation(target_loc_ldmx[0], target_loc_ldmx[1], target_loc_ldmx[2]); @@ -536,7 +547,7 @@ void CKFProcessor::produce(framework::Event& event) { ldmx_log(debug) << " Checking Track State index_ = " << trk_state_index << " at location " << ts.referenceSurface() - .transform(geometryContext()) + .localToGlobalTransform(geometryContext()) .translation() .transpose(); @@ -550,11 +561,10 @@ void CKFProcessor::produce(framework::Event& event) { // Check if the track state is a measurement auto type_flags = ts.typeFlags(); - if (type_flags.test(Acts::TrackStateFlag::MeasurementFlag) && - ts.hasUncalibratedSourceLink()) { - const acts_examples::IndexSourceLink sl = - ts.getUncalibratedSourceLink() - .template get(); + if (type_flags.isMeasurement() && ts.hasUncalibratedSourceLink()) { + Acts::SourceLink usl = ts.getUncalibratedSourceLink(); + const acts_examples::IndexSourceLink& sl = + usl.get(); ldmx::Measurement ldmx_meas = measurements.at(sl.index()); ldmx_log(debug) << " Adding measurement to ldmx::track with " @@ -595,7 +605,7 @@ void CKFProcessor::produce(framework::Event& event) { // Get the local frame (transform from global to local) auto local_frame_transform = - meas_surface.transform(geometryContext()); + meas_surface.localToGlobalTransform(geometryContext()); Acts::Vector3 local_momentum = local_frame_transform.rotation().transpose() * global_momentum; @@ -612,8 +622,10 @@ void CKFProcessor::produce(framework::Event& event) { // cos(angle) = 1 / sqrt(1 + tan(angle)^2) // path_length = thickness / cos(angle) float sensor_thickness = 0.0f; - if (const auto* det_el = meas_surface.associatedDetectorElement()) { - sensor_thickness = static_cast(det_el->thickness()); + if (const auto* placement = meas_surface.surfacePlacement()) { + sensor_thickness = static_cast( + static_cast(placement) + ->thickness()); } else { ldmx_log(warn) << "No detector element for measurement surface" << " — skipping dE/dx for this hit"; diff --git a/Tracking/src/Tracking/Reco/DigitizationProcessor.cxx b/Tracking/src/Tracking/Reco/DigitizationProcessor.cxx index a066d50cc..1daea5a45 100644 --- a/Tracking/src/Tracking/Reco/DigitizationProcessor.cxx +++ b/Tracking/src/Tracking/Reco/DigitizationProcessor.cxx @@ -67,7 +67,7 @@ void DigitizationProcessor::onProcessStart() { std::ofstream csv(dump_geo_csv_); csv << "layer_id,cx,cy,cz,Ux,Uy,Uz,Vx,Vy,Vz,Wx,Wy,Wz\n"; for (const auto& [layer_id, surface] : geometry().layer_surface_map_) { - const auto& xf = surface->transform(geometryContext()); + const auto& xf = surface->localToGlobalTransform(geometryContext()); const auto ctr = xf.translation(); // centre [mm in Acts units] const auto r = xf.rotation(); const auto u = r.col(0); @@ -162,7 +162,7 @@ void DigitizationProcessor::buildLorentzCache() { // Sensor W-normal = 3rd column of the rotation matrix const Acts::Vector3 w_hat = - surface->transform(geometryContext()).rotation().col(2); + surface->localToGlobalTransform(geometryContext()).rotation().col(2); const double bw = b_t.dot(w_hat); // [T] @@ -356,10 +356,11 @@ std::vector DigitizationProcessor::digitizeHits( auto hit_surface{geometry().getSurface(layer_id)}; if (!hit_surface) continue; - ldmx_log(trace) << "Local to global\n" - << hit_surface->transform(geometryContext()).rotation() - << "\n" - << hit_surface->transform(geometryContext()).translation(); + ldmx_log(trace) + << "Local to global\n" + << hit_surface->localToGlobalTransform(geometryContext()).rotation() + << "\n" + << hit_surface->localToGlobalTransform(geometryContext()).translation(); // ----------------------------------------------------------------------- // Project global hit position onto the surface (2D local coords) @@ -392,18 +393,20 @@ std::vector DigitizationProcessor::digitizeHits( // ----------------------------------------------------------------------- if (use_charge_digitization_) { // Read sensor thickness from the geometry. - const auto* det_el = hit_surface->associatedDetectorElement(); - if (!det_el) { + const auto* placement = hit_surface->surfacePlacement(); + if (!placement) { ldmx_log(warn) << "No detector element for layer_id=" << layer_id << " — skipping hit"; continue; } - const double thickness = det_el->thickness(); + const double thickness = + static_cast(placement) + ->thickness(); strip_digitizer_->setThickness(thickness); // Build the full 3D local position and direction for charge simulation. const Acts::Transform3 surf_transform = - hit_surface->transform(geometryContext()); + hit_surface->localToGlobalTransform(geometryContext()); // 3D local position: apply the inverse surface transform to the global // hit position so that we know the depth (W) coordinate. diff --git a/Tracking/src/Tracking/Reco/GSFProcessor.cxx b/Tracking/src/Tracking/Reco/GSFProcessor.cxx index 88c0eeed6..7b72c9ca9 100644 --- a/Tracking/src/Tracking/Reco/GSFProcessor.cxx +++ b/Tracking/src/Tracking/Reco/GSFProcessor.cxx @@ -1,6 +1,8 @@ #include "Tracking/Reco/GSFProcessor.h" #include +#include +#include #include "Acts/EventData/SourceLink.hpp" #include "Tracking/Event/Track.h" @@ -13,6 +15,9 @@ GSFProcessor::GSFProcessor(const std::string& name, framework::Process& process) void GSFProcessor::onNewRun(const ldmx::RunHeader& rh) { beam_origin_surface_ = tracking::sim::utils::unboundSurface(-700); + // 1mm inside tagger ACTS volume outer boundary (~-618mm), 1.5mm upstream of + // L1 at x=-615.5mm + tagger_start_surface_ = tracking::sim::utils::unboundSurface(-617.); target_surface_ = tracking::sim::utils::unboundSurface(0.); ecal_surface_ = tracking::sim::utils::unboundSurface(240.5); @@ -56,7 +61,8 @@ void GSFProcessor::onNewRun(const ldmx::RunHeader& rh) { GsfPropagator(std::move(multi_stepper), std::move(navigator), Acts::getDefaultLogger("GSF_PROP", acts_logging_level)); - BetheHeitlerApprox bethe_heitler = Acts::makeDefaultBetheHeitlerApprox(); + auto bethe_heitler = std::make_shared( + Acts::makeDefaultBetheHeitlerApprox()); gsf_ = std::make_unique>( std::move(gsf_propagator), std::move(bethe_heitler), @@ -67,8 +73,10 @@ void GSFProcessor::onNewRun(const ldmx::RunHeader& rh) { stepper, navigator, Acts::getDefaultLogger("GSF_EXTRAP", acts_logging_level)); + propagator_extrap_ = std::make_unique( + Acts::EigenStepper<>{map}, Acts::VoidNavigator{}); trk_extrap_ = std::make_shared>( - *propagator_, geometryContext(), magneticFieldContext()); + *propagator_extrap_, geometryContext(), magneticFieldContext()); } void GSFProcessor::configure(framework::config::Parameters& parameters) { @@ -106,6 +114,8 @@ void GSFProcessor::configure(framework::config::Parameters& parameters) { } // end of configure() void GSFProcessor::produce(framework::Event& event) { + auto t_start = std::chrono::high_resolution_clock::now(); + // General Setup auto tg{geometry()}; @@ -154,7 +164,7 @@ void GSFProcessor::produce(framework::Event& event) { // Move this at the start of the producer Acts::PropagatorOptions + Acts::NavigatorPlainOptions, ActionList> propagator_options(geometryContext(), magneticFieldContext()); propagator_options.pathLimit = std::numeric_limits::max(); @@ -165,14 +175,14 @@ void GSFProcessor::produce(framework::Event& event) { // Switch the material interaction on/off & eventually into logging mode auto& m_interactor = - propagator_options.actionList.get(); + propagator_options.actorList.get(); m_interactor.multipleScattering = true; m_interactor.energyLoss = true; m_interactor.recordInteractions = false; // The logger can be switched to sterile, e.g. for timing logging auto& s_logger = - propagator_options.actionList.get(); + propagator_options.actorList.get(); s_logger.sterile = true; // Set a maximum step size propagator_options.stepping.maxStepSize = @@ -187,7 +197,8 @@ void GSFProcessor::produce(framework::Event& event) { Acts::GsfOptions gsf_options{ geometryContext(), magneticFieldContext(), calibrationContext()}; gsf_options.extensions = gsf_extensions; - gsf_options.propagatorPlainOptions = propagator_options; + gsf_options.propagatorPlainOptions = + static_cast(propagator_options); gsf_options.maxComponents = max_components_; gsf_options.weightCutoff = weight_cutoff_; gsf_options.abortOnError = abort_on_error_; @@ -201,7 +212,11 @@ void GSFProcessor::produce(framework::Event& event) { Acts::TrackContainer tc{vtc, mtj}; // Loop on tracks + n_input_tracks_ += static_cast(tracks.size()); + unsigned int itrk = 0; + int n_gsf_ok_evt = 0; + int n_tgt_fail_evt = 0; ldmx_log(debug) << "Starting GSF processing of " << tracks.size() << " tracks"; @@ -249,19 +264,24 @@ void GSFProcessor::produce(framework::Event& event) { Acts::BoundTrackParameters trk_btp = tracking::sim::utils::boundTrackParameters(track, perigee); - // GSF starting parameters: for tagger, extrapolate back to beam origin; - // for recoil, use target parameters directly. Acts::BoundTrackParameters trk_btp_fit_start = trk_btp; + // For tagger: backward-extrapolate (via VoidNavigator) from the target + // perigee (x=0mm) to just inside the tagger outer boundary (x≈-650mm), then + // run the GSF forward (+x) through L1→L7. The CKF stores tagger track + // perigees at the target, so we must back-propagate before handing off to + // the GSF. if (tagger_tracking_) { - auto opt_beam_origin = - trk_extrap_->extrapolate(trk_btp, beam_origin_surface_); - if (!opt_beam_origin) { - ldmx_log(warn) << "Failed extrapolating to beam origin for GSF start. " - "Skipping.."; + auto opt_tagger_start = + trk_extrap_->extrapolate(trk_btp, tagger_start_surface_); + if (!opt_tagger_start) { + ldmx_log(debug) + << " Failed pre-fit extrapolation to tagger start surface (itrk=" + << itrk << ")"; + ++n_gsf_failed_; continue; } - trk_btp_fit_start = *opt_beam_origin; + trk_btp_fit_start = *opt_tagger_start; } ldmx_log(debug) << " Perigee surface (acts): (" << track.getPerigeeX() @@ -287,9 +307,10 @@ void GSFProcessor::produce(framework::Event& event) { ldmx_log(debug) << " About to run GSF fit with " << fit_track_source_links.size() << " source links"; - // Update GSF reference surface for this track + // GSF reference surface: for tagger use the start surface (x=-648mm, inside + // geometry), for recoil use the target (x=0mm). if (tagger_tracking_) { - gsf_ref_surface = beam_origin_surface_; + gsf_ref_surface = tagger_start_surface_; } else { gsf_ref_surface = target_surface_; } @@ -300,30 +321,43 @@ void GSFProcessor::produce(framework::Event& event) { trk_btp_fit_start, gsf_options, tc); if (!gsf_refit_result.ok()) { - ldmx_log(warn) << "GSF re-fit failed: " - << gsf_refit_result.error().message(); + ldmx_log(debug) << " GSF re-fit failed (itrk=" << itrk + << "): " << gsf_refit_result.error().message(); + if (n_gsf_failed_ < 5) + ldmx_log(info) << " [GSF dbg] fit failed (first few): " + << gsf_refit_result.error().message(); + ++n_gsf_failed_; continue; } - ldmx_log(debug) << " GSF fit succeeded, tc.size() = " << tc.size(); - - if (tc.size() < 1) continue; + ++n_gsf_ok_evt; + ldmx_log(debug) << " GSF fit succeeded (itrk=" << itrk + << "), tc.size()=" << tc.size(); - auto gsftrk = tc.getTrack(0); + auto gsftrk = gsf_refit_result.value(); // calculateTrackQuantities(gsftrk); const Acts::BoundVector& perigee_pars = gsftrk.parameters(); const Acts::BoundMatrix& trk_cov = gsftrk.covariance(); const Acts::Surface& perigee_surface = gsftrk.referenceSurface(); - ldmx_log(debug) - << " Reference Surface (acts-x, acts-y, acts-z) = (" - << perigee_surface.transform(geometryContext()).translation()(0) << ", " - << perigee_surface.transform(geometryContext()).translation()(1) << ", " - << perigee_surface.transform(geometryContext()).translation()(2) << ")"; - - ldmx_log(debug) << " Found track has " << gsftrk.nTrackStates() - << " track states"; + ldmx_log(debug) << " Reference Surface (acts-x, acts-y, acts-z) = (" + << perigee_surface.localToGlobalTransform(geometryContext()) + .translation()(0) + << ", " + << perigee_surface.localToGlobalTransform(geometryContext()) + .translation()(1) + << ", " + << perigee_surface.localToGlobalTransform(geometryContext()) + .translation()(2) + << ")"; + + ldmx_log(debug) << " nTrackStates=" << gsftrk.nTrackStates() + << " nMeasurements=" << gsftrk.nMeasurements() + << " chi2=" << gsftrk.chi2(); + if (gsftrk.nTrackStates() == 0) + ldmx_log(info) << " [GSF dbg] track has 0 states after fit (itrk=" + << itrk << ");"; ldmx_log(debug) << " Track parameters (d0, z0, phi, theta, q/p)= (" << perigee_pars[Acts::eBoundLoc0] << ", " @@ -337,6 +371,7 @@ void GSFProcessor::produce(framework::Event& event) { // Extrapolate GSF track to target surface to get perigee parameters auto opt_target = trk_extrap_->extrapolate(gsftrk, target_surface_); + ldmx_log(debug) << " Extrapolating to target (itrk=" << itrk << ")"; if (opt_target) { ldmx_log(debug) << " GSF target extrapolation succeeded"; auto ts_at_target = tracking::sim::utils::makeTrackState( @@ -351,7 +386,8 @@ void GSFProcessor::produce(framework::Event& event) { trk.setPerigeeCov(cov_vec); } Acts::Vector3 target_loc_ldmx = tracking::sim::utils::acts2Ldmx( - target_surface_->transform(geometryContext()).translation()); + target_surface_->localToGlobalTransform(geometryContext()) + .translation()); trk.setPerigeeLocation(target_loc_ldmx[0], target_loc_ldmx[1], target_loc_ldmx[2]); @@ -363,8 +399,10 @@ void GSFProcessor::produce(framework::Event& event) { << opt_target->parameters()[Acts::eBoundTheta] << ", " << opt_target->parameters()[Acts::eBoundQOverP] << ")"; } else { - ldmx_log(debug) << " GSF target extrapolation failed, using GSF fit " - "parameters at reference surface"; + ++n_target_extrap_failed_; + ++n_tgt_fail_evt; + ldmx_log(debug) << " GSF target extrapolation failed (itrk=" << itrk + << "), using GSF fit parameters at reference surface"; trk.setPerigeeParameters( tracking::sim::utils::convertActsToLdmxPars(perigee_pars)); std::vector v_trk_cov; @@ -385,6 +423,8 @@ void GSFProcessor::produce(framework::Event& event) { if (opt_ecal) trk.addTrackState(tracking::sim::utils::makeTrackState( geometryContext(), *opt_ecal, ldmx::AtECAL)); + else + ++n_ecal_extrap_failed_; } trk.setChi2(gsftrk.chi2()); @@ -406,11 +446,34 @@ void GSFProcessor::produce(framework::Event& event) { } // loop on tracks + ldmx_log(debug) << "[GSF evt " << nevents_ << "] in=" << tracks.size() + << " gsf_ok=" << n_gsf_ok_evt + << " tgt_fail=" << n_tgt_fail_evt + << " out=" << out_tracks.size(); + + n_output_tracks_ += static_cast(out_tracks.size()); event.add(out_trk_collection_, out_tracks); + + auto t_end = std::chrono::high_resolution_clock::now(); + processing_time_ += + std::chrono::duration(t_end - t_start).count(); + ++nevents_; } // end of produce() -void GSFProcessor::onProcessStart() {}; -void GSFProcessor::onProcessEnd() {}; +void GSFProcessor::onProcessStart() {} + +void GSFProcessor::onProcessEnd() { + ldmx_log(info) << "--------------------------------- "; + ldmx_log(info) << "GSF: " << n_output_tracks_ << " output tracks / " + << n_input_tracks_ << " input tracks"; + ldmx_log(info) << "AVG Time/Event: " << std::fixed << std::setprecision(1) + << processing_time_ / nevents_ << " ms"; + ldmx_log(info) << "GSF Fit Failures: " << n_gsf_failed_; + ldmx_log(info) << "Extrapolation Failures::"; + ldmx_log(info) << " Target: " << n_target_extrap_failed_ << " times"; + if (!tagger_tracking_) + ldmx_log(info) << " ECAL: " << n_ecal_extrap_failed_ << " times"; +} } // namespace reco } // namespace tracking diff --git a/Tracking/src/Tracking/Reco/GreedyAmbiguitySolver.cxx b/Tracking/src/Tracking/Reco/GreedyAmbiguitySolver.cxx index 25b4452c6..57412f860 100644 --- a/Tracking/src/Tracking/Reco/GreedyAmbiguitySolver.cxx +++ b/Tracking/src/Tracking/Reco/GreedyAmbiguitySolver.cxx @@ -1,9 +1,11 @@ #include "Tracking/Reco/GreedyAmbiguitySolver.h" #include +#include +#include #include "Acts/EventData/SourceLink.hpp" -#include "Acts/EventData/TrackHelpers.hpp" +#include "Acts/Utilities/TrackHelpers.hpp" namespace tracking { namespace reco { @@ -177,6 +179,8 @@ void GreedyAmbiguitySolver::configure( } void GreedyAmbiguitySolver::produce(framework::Event& event) { + auto t_start = std::chrono::high_resolution_clock::now(); + GreedyAmbiguitySolver::State state; std::vector out_tracks; @@ -184,6 +188,7 @@ void GreedyAmbiguitySolver::produce(framework::Event& event) { if (!event.exists(track_collection_, input_pass_name_)) { ldmx_log(debug) << "Track collection not found, exiting"; + ++nevents_; return; } const auto& tracks = @@ -191,11 +196,14 @@ void GreedyAmbiguitySolver::produce(framework::Event& event) { if (!event.exists(meas_collection_, input_pass_name_)) { ldmx_log(debug) << "Measurement collection not found, exiting"; + ++nevents_; return; } const auto& measurements = event.getCollection( meas_collection_, input_pass_name_); + n_input_tracks_ += static_cast(tracks.size()); + computeInitialState(tracks, measurements, state, tg, tracking::sim::utils::sourceLinkHash, tracking::sim::utils::sourceLinkEquality); @@ -209,16 +217,33 @@ void GreedyAmbiguitySolver::produce(framework::Event& event) { } } + n_output_tracks_ += static_cast(out_tracks.size()); + event.add(out_trk_collection_, out_tracks); - // for (auto iTrack : initial_state.selectedTracks) { - // std::cout << event.getEventNumber() << " " << iTrack << " " << - // initial_state.trackChi2[iTrack] << " " << - // initial_state.measurementsPerTrack[iTrack].size() << std::endl; - // } + auto t_end = std::chrono::high_resolution_clock::now(); + processing_time_ += + std::chrono::duration(t_end - t_start).count(); + ++nevents_; +} - ldmx_log(info) << " Resolved to " << state.selected_tracks_.size() - << " tracks from " << " " << tracks.size(); +void GreedyAmbiguitySolver::onProcessEnd() { + double avg_in = + nevents_ > 0 ? static_cast(n_input_tracks_) / nevents_ : 0.; + double avg_out = + nevents_ > 0 ? static_cast(n_output_tracks_) / nevents_ : 0.; + double retention = + n_input_tracks_ > 0 ? 100.0 * n_output_tracks_ / n_input_tracks_ : 0.; + ldmx_log(info) << "--------------------------------- "; + ldmx_log(info) << "GAS: " << n_output_tracks_ << " output tracks / " + << n_input_tracks_ << " input tracks"; + ldmx_log(info) << "AVG Time/Event: " << std::fixed << std::setprecision(1) + << processing_time_ / nevents_ << " ms"; + ldmx_log(info) << "AVG tracks in/event: " << std::fixed + << std::setprecision(1) << avg_in; + ldmx_log(info) << "AVG tracks out/event: " << std::fixed + << std::setprecision(1) << avg_out << " (" << std::fixed + << std::setprecision(1) << retention << "% retained)"; } } // namespace reco diff --git a/Tracking/src/Tracking/Reco/SeedFinderProcessor.cxx b/Tracking/src/Tracking/Reco/SeedFinderProcessor.cxx index 47b572128..b9452b6bb 100644 --- a/Tracking/src/Tracking/Reco/SeedFinderProcessor.cxx +++ b/Tracking/src/Tracking/Reco/SeedFinderProcessor.cxx @@ -125,8 +125,7 @@ void SeedFinderProcessor::produce(framework::Event& event) { const auto& perigee_cov = tagtrk.getPerigeeCov(); if (!perigee_cov.empty()) { - Acts::BoundSquareMatrix cov = - tracking::sim::utils::unpackCov(perigee_cov); + Acts::BoundMatrix cov = tracking::sim::utils::unpackCov(perigee_cov); double locu = tagtrk.getD0(); double locv = tagtrk.getZ0(); double covuu = @@ -236,8 +235,8 @@ ldmx::Track SeedFinderProcessor::seedTracker( // In this way it's easier to incorporate the tagger track extrapolation to // the fit - Acts::ActsMatrix<5, 5> a = Acts::ActsMatrix<5, 5>::Zero(); - Acts::ActsVector<5> y = Acts::ActsVector<5>::Zero(); + Acts::Matrix<5, 5> a = Acts::Matrix<5, 5>::Zero(); + Acts::Vector<5> y = Acts::Vector<5>::Zero(); for (auto meas : vmeas) { double xmeas = meas.getGlobalPosition()[0] - xOrigin; @@ -246,8 +245,10 @@ ldmx::Track SeedFinderProcessor::seedTracker( const Acts::Surface* hit_surface = geometry().getSurface(meas.getLayerID()); // Get the global to local transformation - auto rot = hit_surface->transform(geometryContext()).rotation(); - auto tr = hit_surface->transform(geometryContext()).translation(); + auto rot = + hit_surface->localToGlobalTransform(geometryContext()).rotation(); + auto tr = + hit_surface->localToGlobalTransform(geometryContext()).translation(); auto rotl2g = rot.transpose(); @@ -258,7 +259,7 @@ ldmx::Track SeedFinderProcessor::seedTracker( yhit_.push_back(meas.getGlobalPosition()[1]); zhit_.push_back(meas.getGlobalPosition()[2]); - Acts::ActsMatrix<2, 5> a_i; + Acts::Matrix<2, 5> a_i; a_i(0, 0) = rotl2g(0, 1); a_i(0, 1) = rotl2g(0, 1) * xmeas; @@ -279,7 +280,7 @@ ldmx::Track SeedFinderProcessor::seedTracker( loc(0) = meas.getLocalPosition()[0]; loc(1) = 0.; // weight matrix - Acts::ActsMatrix<2, 2> w_i = Acts::ActsMatrix<2, 2>::Zero(); + Acts::Matrix<2, 2> w_i = Acts::Matrix<2, 2>::Zero(); w_i(0, 0) = 1. / (u_error_ * u_error_); w_i(1, 1) = 1. / (v_error_ * v_error_); @@ -287,11 +288,11 @@ ldmx::Track SeedFinderProcessor::seedTracker( Acts::Vector2 yprime_i = loc + offset - xoffset; y += (a_i.transpose()) * w_i * yprime_i; - Acts::ActsMatrix<2, 5> wa_i = (w_i * a_i); + Acts::Matrix<2, 5> wa_i = (w_i * a_i); a += a_i.transpose() * wa_i; } - Acts::ActsVector<5> b; + Acts::Vector<5> b; b = a.inverse() * y; b0_.push_back(b(0)); @@ -300,17 +301,22 @@ ldmx::Track SeedFinderProcessor::seedTracker( b3_.push_back(b(3)); b4_.push_back(b(4)); - // Acts::ActsVector<5> hlx = Acts::ActsVector<5>::Zero(); - Acts::ActsVector<3> ref{0., 0., 0.}; + // Acts::Vector<5> hlx = Acts::Vector<5>::Zero(); + Acts::Vector<3> ref{0., 0., 0.}; + // relative_perigee_x is the perigee position in the fit frame (fit-x = ACTS x + // - xOrigin). It is used only for evaluating the fitted curve (y, z, slopes). + // The PerigeeSurface and seed_pos must use the absolute ACTS x coordinate, + // which is perigee_location(0) directly. double relative_perigee_x = perigee_location(0) - xOrigin; std::shared_ptr seed_perigee = Acts::Surface::makeShared(Acts::Vector3( - relative_perigee_x, perigee_location(1), perigee_location(2))); + perigee_location(0), perigee_location(1), perigee_location(2))); - // in mm - Acts::Vector3 seed_pos{relative_perigee_x, + // in mm — x is absolute ACTS x; y and z evaluated at fit-x = + // relative_perigee_x + Acts::Vector3 seed_pos{perigee_location(0), b(0) + b(1) * relative_perigee_x + b(2) * relative_perigee_x * relative_perigee_x, b(3) + b(4) * relative_perigee_x}; @@ -324,7 +330,7 @@ ldmx::Track SeedFinderProcessor::seedTracker( // Convert it to MeV since that's what TrackUtils assumes Acts::Vector3 seed_mom = p * dir / Acts::UnitConstants::MeV; - Acts::ActsScalar q = + double q = b(2) < 0 ? -1 * Acts::UnitConstants::e : +1 * Acts::UnitConstants::e; // Linear intersection with the perigee line. TODO:: Use propagator instead @@ -344,7 +350,7 @@ ldmx::Track SeedFinderProcessor::seedTracker( (*seed_perigee).intersect(geometryContext(), seed_pos, dir); Acts::FreeVector seed_free = tracking::sim::utils::toFreeParameters( - intersection.intersections()[0].position(), seed_mom, q); + intersection[0].position(), seed_mom, q); auto bound_params = Acts::transformFreeToBoundParameters( seed_free, *seed_perigee, geometryContext()) @@ -370,11 +376,13 @@ ldmx::Track SeedFinderProcessor::seedTracker( ldmx_log(debug) << "Making covariance matrix as diagonal matrix with inflated terms"; - Acts::BoundSquareMatrix bound_cov = stddev.cwiseProduct(stddev).asDiagonal(); + Acts::BoundMatrix bound_cov = stddev.cwiseProduct(stddev).asDiagonal(); ldmx_log(debug) << "...now putting together the seed track ..."; ldmx::Track trk = ldmx::Track(); + // Store the perigee surface position (absolute ACTS coordinates) converted to + // LDMX frame so CKFProcessor can reconstruct the same surface. Acts::Vector3 perigee_ldmx = tracking::sim::utils::acts2Ldmx(perigee_location); trk.setPerigeeLocation(perigee_ldmx(0), perigee_ldmx(1), perigee_ldmx(2)); @@ -393,7 +401,7 @@ ldmx::Track SeedFinderProcessor::seedTracker( ldmx_log(debug) << "...making the ParticleHypothesis ...assume electron for now"; - auto part_hypo{Acts::SinglyChargedParticleHypothesis::electron()}; + auto part_hypo{Acts::ParticleHypothesis::electron()}; ldmx_log(debug) << "Making BoundTrackParameters seedParameters"; Acts::BoundTrackParameters seed_parameters( diff --git a/Tracking/src/Tracking/Reco/TruthSeedProcessor.cxx b/Tracking/src/Tracking/Reco/TruthSeedProcessor.cxx index 833eb7053..2f8fa8343 100644 --- a/Tracking/src/Tracking/Reco/TruthSeedProcessor.cxx +++ b/Tracking/src/Tracking/Reco/TruthSeedProcessor.cxx @@ -9,7 +9,6 @@ TruthSeedProcessor::TruthSeedProcessor(const std::string& name, : TrackingGeometryUser(name, process) {} void TruthSeedProcessor::onNewRun(const ldmx::RunHeader& rh) { - gctx_ = Acts::GeometryContext(); normal_ = std::make_shared>(0., 1.); // Custom transformation of the interpolated bfield map @@ -186,12 +185,11 @@ void TruthSeedProcessor::createTruthTrack( // " < v_seed_cov; tracking::sim::utils::flatCov(bound_cov, v_seed_cov); seed.setPerigeeParameters(v_seed_params); @@ -458,8 +454,7 @@ ldmx::Track TruthSeedProcessor::seedFromTruth(const ldmx::Track& tt, stddev[Acts::eBoundTheta] = 5 * Acts::UnitConstants::degree; stddev[Acts::eBoundQOverP] = (1. / p) * (1. / p) * sigma_p; - Acts::BoundSquareMatrix bound_cov = - stddev.cwiseProduct(stddev).asDiagonal(); + Acts::BoundMatrix bound_cov = stddev.cwiseProduct(stddev).asDiagonal(); std::vector v_seed_cov; tracking::sim::utils::flatCov(bound_cov, v_seed_cov); seed.setPerigeeParameters(v_seed_params); @@ -795,13 +790,13 @@ void TruthSeedProcessor::produce(framework::Event& event) { double q_ecal = phit.getCharge() * Acts::UnitConstants::e; auto ecal_free = tracking::sim::utils::toFreeParameters(ep, em, q_ecal); auto ecal_bound = Acts::transformFreeToBoundParameters( - ecal_free, *ecal_surface, gctx_); + ecal_free, *ecal_surface, geometryContext()); if (ecal_bound.ok()) { - auto part{Acts::GenericParticleHypothesis(Acts::ParticleHypothesis( - Acts::PdgParticle(particle_hypothesis_)))}; - Acts::BoundTrackParameters ecal_pars( - ecal_surface, ecal_bound.value(), - Acts::BoundSquareMatrix::Identity(), part); + auto part{Acts::ParticleHypothesis( + Acts::PdgParticle(particle_hypothesis_))}; + Acts::BoundTrackParameters ecal_pars(ecal_surface, ecal_bound.value(), + Acts::BoundMatrix::Identity(), + part); truth_recoil_track.addTrackState(tracking::sim::utils::makeTrackState( geometryContext(), ecal_pars, ldmx::AtECAL)); } diff --git a/Tracking/src/Tracking/Reco/VertexProcessor.cxx b/Tracking/src/Tracking/Reco/VertexProcessor.cxx index f04df9830..a0cd7b3e6 100644 --- a/Tracking/src/Tracking/Reco/VertexProcessor.cxx +++ b/Tracking/src/Tracking/Reco/VertexProcessor.cxx @@ -14,7 +14,6 @@ VertexProcessor::VertexProcessor(const std::string& name, : framework::Producer(name, process) {} void VertexProcessor::onProcessStart() { - gctx_ = Acts::GeometryContext(); bctx_ = Acts::MagneticFieldContext(); h_m_ = new TH1F("m", "m", 100, 0., 1.); @@ -63,28 +62,9 @@ void VertexProcessor::produce(framework::Event& event) { // Set up propagator with void navigator propagator_ = std::make_shared(stepper); - // Track linearizer in the proximity of the vertex location - using Linearizer = Acts::HelicalTrackLinearizer; - Linearizer::Config linearizer_config; - linearizer_config.bField = sp_interpolated_b_field_; - linearizer_config.propagator = propagator_; - Linearizer linearizer(linearizer_config); - - // Set up Billoir Vertex Fitter - using VertexFitter = Acts::FullBilloirVertexFitter; - - VertexFitter::Config vertex_fitter_cfg; - - VertexFitter billoir_fitter(vertex_fitter_cfg); - - // VertexFitter::State state(sp_interpolated_bField_->makeCache(bctx_)); - - // Unconstrained fit - // See - // https://github.com/acts-project/acts/blob/main/Tests/UnitTests/Core/Vertexing/FullBilloirVertexFitterTests.cpp#L149 - // For constraint implementation - - Acts::VertexingOptions vf_options(gctx_, bctx_); + // Note: FullBilloirVertexFitter setup commented out — fit() is not called yet + // and v46 Config now requires extractParameters/trackLinearizer delegates. + // Acts::VertexingOptions vf_options(gctx_, bctx_); // Retrieve the track collection const auto& tracks = @@ -115,10 +95,10 @@ void VertexProcessor::produce(framework::Event& event) { tracks.at(i_track).getPhi(), tracks.at(i_track).getTheta(), tracks.at(i_track).getQoP(), tracks.at(i_track).getT(); - Acts::BoundSquareMatrix cov_mat = + Acts::BoundMatrix cov_mat = tracking::sim::utils::unpackCov(tracks.at(i_track).getPerigeeCov()); - auto part{Acts::GenericParticleHypothesis(Acts::ParticleHypothesis( - Acts::PdgParticle(tracks.at(i_track).getPdgID())))}; + auto part{Acts::ParticleHypothesis( + Acts::PdgParticle(tracks.at(i_track).getPdgID()))}; billoir_tracks.push_back(Acts::BoundTrackParameters( perigee_surface, param_vec, std::move(cov_mat), part)); } @@ -158,13 +138,12 @@ void VertexProcessor::produce(framework::Event& event) { seeds.at(i_seed).getPhi(), seeds.at(i_seed).getTheta(), seeds.at(i_seed).getQoP(), seeds.at(i_seed).getT(); - Acts::BoundSquareMatrix cov_mat = + Acts::BoundMatrix cov_mat = tracking::sim::utils::unpackCov(seeds.at(i_seed).getPerigeeCov()); int pion_pdg_id = 211; // pi+ if (seeds.at(i_seed).getCharge() < 0) pion_pdg_id = -211; // BoundTrackParameters needs the particle hypothesis - auto part{Acts::GenericParticleHypothesis( - Acts::ParticleHypothesis(Acts::PdgParticle(pion_pdg_id)))}; + auto part{Acts::ParticleHypothesis(Acts::PdgParticle(pion_pdg_id))}; auto bound_seed_params = Acts::BoundTrackParameters( perigee_surface, param_vec, std::move(cov_mat), part); diff --git a/Tracking/src/Tracking/Reco/Vertexer.cxx b/Tracking/src/Tracking/Reco/Vertexer.cxx index a9567d25f..3f5ec63b4 100644 --- a/Tracking/src/Tracking/Reco/Vertexer.cxx +++ b/Tracking/src/Tracking/Reco/Vertexer.cxx @@ -44,7 +44,6 @@ void Vertexer::onProcessStart() { h_tz0_vs_rz0_ = new TH2F("h_tz0_vs_rz0", "h_tz0_vs_rz0", 100, -40, 40, 100, -40, 40); - gctx_ = Acts::GeometryContext(); bctx_ = Acts::MagneticFieldContext(); /* @@ -90,34 +89,9 @@ void Vertexer::produce(framework::Event& event) { nevents_++; // auto start = std::chrono::high_resolution_clock::now(); - // Track linearizer in the proximity of the vertex location - using Linearizer = Acts::HelicalTrackLinearizer; - Linearizer::Config linearizer_config; - linearizer_config.bField = b_field_; - linearizer_config.propagator = propagator_; - Linearizer linearizer(linearizer_config); - - // Set up Billoir Vertex Fitter - using VertexFitter = Acts::FullBilloirVertexFitter; - - // Alternatively one can use - // using VertexFitter = - // Acts::FullBilloirVertexFitter; - - VertexFitter::Config vertex_fitter_cfg; - VertexFitter billoir_fitter(vertex_fitter_cfg); - // mg Aug 2024 .. State doesn't exist in v36 and isn't used here anyway - // VertexFitter::State state(sp_interpolated_bField_->makeCache(bctx_)); - - // Unconstrained fit - // See - // https://github.com/acts-project/acts/blob/main/Tests/UnitTests/Core/Vertexing/FullBilloirVertexFitterTests.cpp#L149 - // For constraint implementation - - // Acts::VertexingOptions vfOptions(gctx_, - // bctx_); - // mg Aug 2024 ... VertexingOptions template change in v36 - Acts::VertexingOptions vf_options(gctx_, bctx_); + // Note: FullBilloirVertexFitter setup commented out — fit() is not called + // and v46 Config now requires extractParameters/trackLinearizer delegates. + // Acts::VertexingOptions vf_options(gctx_, bctx_); // Retrive the two track collections diff --git a/Tracking/src/Tracking/Sim/PropagatorStepWriter.cxx b/Tracking/src/Tracking/Sim/PropagatorStepWriter.cxx index 858232f91..2a3e10ee9 100644 --- a/Tracking/src/Tracking/Sim/PropagatorStepWriter.cxx +++ b/Tracking/src/Tracking/Sim/PropagatorStepWriter.cxx @@ -6,7 +6,6 @@ #include #include -#include "Acts/Geometry/DetectorElementBase.hpp" #include "Framework/Exception/Exception.h" // mg ... I don't think these are used, and they are not defined in acts v36 // #include "Acts/Plugins/Identification/IdentifiedDetectorElement.hpp" @@ -166,9 +165,10 @@ bool PropagatorStepWriter::writeSteps( // double accuracy = // step.stepSize.value(Acts::ConstrainedStep::accuracy()); double accuracy = step.stepSize.accuracy(); - double actor = step.stepSize.value(Acts::ConstrainedStep::actor); - double aborter = step.stepSize.value(Acts::ConstrainedStep::aborter); - double user = step.stepSize.value(Acts::ConstrainedStep::user); + double actor = step.stepSize.value(Acts::ConstrainedStep::Type::Actor); + double aborter = + step.stepSize.value(Acts::ConstrainedStep::Type::Navigator); + double user = step.stepSize.value(Acts::ConstrainedStep::Type::User); double act2 = actor * actor; double acc2 = accuracy * accuracy; double abo2 = aborter * aborter; diff --git a/Tracking/src/Tracking/Sim/TrackingUtils.cxx b/Tracking/src/Tracking/Sim/TrackingUtils.cxx index 48f23870a..06e7e3548 100644 --- a/Tracking/src/Tracking/Sim/TrackingUtils.cxx +++ b/Tracking/src/Tracking/Sim/TrackingUtils.cxx @@ -76,15 +76,15 @@ ldmx::LdmxSpacePoint* convertSimHitToLdmxSpacePoint( sigma_v * sigma_v, hit.getID()); } -void flatCov(Acts::BoundSquareMatrix cov, std::vector& v_cov) { +void flatCov(Acts::BoundMatrix cov, std::vector& v_cov) { v_cov.clear(); v_cov.reserve(cov.rows() * (cov.rows() + 1) / 2); for (int i = 0; i < cov.rows(); i++) for (int j = i; j < cov.cols(); j++) v_cov.push_back(cov(i, j)); } -Acts::BoundSquareMatrix unpackCov(const std::vector& v_cov) { - Acts::BoundSquareMatrix cov; +Acts::BoundMatrix unpackCov(const std::vector& v_cov) { + Acts::BoundMatrix cov; int e{0}; for (int i = 0; i < cov.rows(); i++) for (int j = i; j < cov.cols(); j++) { @@ -121,9 +121,9 @@ Acts::Vector3 acts2Ldmx(Acts::Vector3 acts_v) { // Transform position, momentum and charge to free parameters Acts::FreeVector toFreeParameters(Acts::Vector3 pos_, Acts::Vector3 mom, - Acts::ActsScalar q) { + double q) { Acts::FreeVector free_params; - Acts::ActsScalar p = mom.norm() * Acts::UnitConstants::MeV; + double p = mom.norm() * Acts::UnitConstants::MeV; free_params[Acts::eFreePos0] = pos_(Acts::ePos0) * Acts::UnitConstants::mm; free_params[Acts::eFreePos1] = pos_(Acts::ePos1) * Acts::UnitConstants::mm; @@ -133,7 +133,7 @@ Acts::FreeVector toFreeParameters(Acts::Vector3 pos_, Acts::Vector3 mom, free_params[Acts::eFreeDir1] = mom(1) / mom.norm(); free_params[Acts::eFreeDir2] = mom(2) / mom.norm(); free_params[Acts::eFreeQOverP] = - (q != Acts::ActsScalar(0)) ? (q / p) : 0.; // 1. / p instead? + (q != double(0)) ? (q / p) : 0.; // 1. / p instead? return free_params; } @@ -156,8 +156,8 @@ Acts::BoundVector boundState(const ldmx::Track& trk) { Acts::BoundTrackParameters boundTrackParameters( const ldmx::Track& trk, std::shared_ptr perigee) { Acts::BoundVector param_vec = boundState(trk); - Acts::BoundSquareMatrix cov_mat = unpackCov(trk.getPerigeeCov()); - auto part_hypo{Acts::SinglyChargedParticleHypothesis::electron()}; + Acts::BoundMatrix cov_mat = unpackCov(trk.getPerigeeCov()); + auto part_hypo{Acts::ParticleHypothesis::electron()}; return Acts::BoundTrackParameters(perigee, param_vec, std::move(cov_mat), part_hypo); } @@ -247,7 +247,7 @@ ldmx::Track::TrackState makeTrackState( const Acts::BoundToFreeMatrix j_btf = bound_pars.referenceSurface().boundToFreeJacobian(gctx, acts_pos, acts_dir); - const Acts::FreeSquareMatrix free_cov = + const Acts::FreeMatrix free_cov = j_btf * bound_cov.value() * j_btf.transpose(); // Step 2: Drop time row/col (eFreeTime = 3) -> 7x7 diff --git a/Tracking/src/Tracking/dqm/StraightTracksDQM.cxx b/Tracking/src/Tracking/dqm/StraightTracksDQM.cxx index 1bd72b05f..70f8a3669 100644 --- a/Tracking/src/Tracking/dqm/StraightTracksDQM.cxx +++ b/Tracking/src/Tracking/dqm/StraightTracksDQM.cxx @@ -135,17 +135,16 @@ void StraightTracksDQM::trackMonitoringUnique( double track_state_loc1_ecal = track.getEcalLayer1Y(); int track_pdg_id = track.getPdgID(); - double sigma_phi = phiAngleError(track.getSlopeX(), track.getCov()); + const std::vector cov = track.getCov(); + double sigma_phi = phiAngleError(track.getSlopeX(), cov); double sigma_theta = - thetaAngleError(track.getSlopeX(), track.getSlopeY(), track.getCov()); - double sigma_loc0_target = std::sqrt(track.getCov()[4]); - double sigma_loc1_target = std::sqrt(track.getCov()[9]); + thetaAngleError(track.getSlopeX(), track.getSlopeY(), cov); + double sigma_loc0_target = std::sqrt(cov[4]); + double sigma_loc1_target = std::sqrt(cov[9]); double sigma_loc0_ecal = - locError(track.getCov().at(0), track.getCov().at(4), - track.getCov().at(1), track.getEcalLayer1Z()); + locError(cov.at(0), cov.at(4), cov.at(1), track.getEcalLayer1Z()); double sigma_loc1_ecal = - locError(track.getCov().at(7), track.getCov().at(9), - track.getCov().at(8), track.getEcalLayer1Z()); + locError(cov.at(7), cov.at(9), cov.at(8), track.getEcalLayer1Z()); histograms_.fill(title + "phi", trk_phi); histograms_.fill(title + "theta", trk_theta); diff --git a/Tracking/src/Tracking/dqm/TrackingRecoDQM.cxx b/Tracking/src/Tracking/dqm/TrackingRecoDQM.cxx index 3ec5b2928..59ba3fbac 100644 --- a/Tracking/src/Tracking/dqm/TrackingRecoDQM.cxx +++ b/Tracking/src/Tracking/dqm/TrackingRecoDQM.cxx @@ -382,7 +382,7 @@ void TrackingRecoDQM::trackMonitoring( } // Covariance matrix - Acts::BoundSquareMatrix cov = + Acts::BoundMatrix cov = tracking::sim::utils::unpackCov(track.getPerigeeCov()); double sigmad0 = sqrt( diff --git a/Tracking/src/Tracking/geo/DetectorElement.cxx b/Tracking/src/Tracking/geo/DetectorElement.cxx index 3d335c318..d4dd6e107 100644 --- a/Tracking/src/Tracking/geo/DetectorElement.cxx +++ b/Tracking/src/Tracking/geo/DetectorElement.cxx @@ -7,7 +7,7 @@ namespace tracking::geo { DetectorElement::~DetectorElement() {}; -const Acts::Transform3& DetectorElement::transform( +const Acts::Transform3& DetectorElement::localToGlobalTransform( const Acts::GeometryContext& gctx) const { if (!m_surface_) EXCEPTION_RAISE("BadGeometry", diff --git a/Tracking/src/Tracking/geo/GeometryContext.cxx b/Tracking/src/Tracking/geo/GeometryContext.cxx index 3d76f15d3..689a9a57d 100644 --- a/Tracking/src/Tracking/geo/GeometryContext.cxx +++ b/Tracking/src/Tracking/geo/GeometryContext.cxx @@ -8,9 +8,9 @@ namespace tracking::geo { const std::string GeometryContext::NAME = "TrackingGeometryContext"; -GeometryContext::GeometryContext() : framework::ConditionsObject(NAME) { - acts_gc_ = this; -} +GeometryContext::GeometryContext() + : framework::ConditionsObject(NAME), + acts_gc_(Acts::GeometryContext(this)) {} const Acts::GeometryContext& GeometryContext::get() const { return acts_gc_; } @@ -20,8 +20,7 @@ void GeometryContext::loadTransformations(const tgSurfMap& surf_map) { for (auto entry : surf_map) { alignment_map_[entry.first] = - static_cast( - (entry.second)->associatedDetectorElement()) + static_cast((entry.second)->surfacePlacement()) ->uncorrectedTransform(); } } diff --git a/Tracking/src/Tracking/geo/TrackersTrackingGeometry.cxx b/Tracking/src/Tracking/geo/TrackersTrackingGeometry.cxx index 40a934954..5e1916c78 100644 --- a/Tracking/src/Tracking/geo/TrackersTrackingGeometry.cxx +++ b/Tracking/src/Tracking/geo/TrackersTrackingGeometry.cxx @@ -20,6 +20,16 @@ TrackersTrackingGeometry::TrackersTrackingGeometry( Acts::CuboidVolumeBuilder::VolumeConfig recoil_volume_cfg = buildVolumeConfig( recoil_, recoil_layout_, tracker_y_length, tracker_z_length, "Recoil"); + // Extend the recoil volume upstream so the low-x (ACTS) edge is at -1mm, + // placing the target (x=0) clearly inside the volume for the CKF Navigator. + { + double downstream_x = + recoil_volume_cfg.position[0] + recoil_volume_cfg.length[0] / 2.0; + constexpr double low_x = -1.0; // mm + recoil_volume_cfg.length[0] = downstream_x - low_x; + recoil_volume_cfg.position[0] = (downstream_x + low_x) / 2.0; + } + std::vector vol_builder_configs{ tagger_volume_cfg, recoil_volume_cfg}; @@ -366,7 +376,7 @@ std::shared_ptr TrackersTrackingGeometry::getSurfacePtr( // After this call each surface will use the underlying detectorElement // transformation which will take care of effectively reading the gctx - surface->assignDetectorElement(std::move(*det_element)); + surface->assignSurfacePlacement(*det_element); det_elements_.push_back(det_element); return surface; @@ -408,7 +418,7 @@ TrackersTrackingGeometry::buildVolumeConfig( sub_det_volume_config.name = volumeName; // Vacuum material - Acts::Material subdet_mat = Acts::Material(); + Acts::Material subdet_mat = Acts::Material::Vacuum(); sub_det_volume_config.volumeMaterial = std::make_shared(subdet_mat); diff --git a/Tracking/src/Tracking/geo/TrackingGeometry.cxx b/Tracking/src/Tracking/geo/TrackingGeometry.cxx index e20c44543..42e3f566e 100644 --- a/Tracking/src/Tracking/geo/TrackingGeometry.cxx +++ b/Tracking/src/Tracking/geo/TrackingGeometry.cxx @@ -191,11 +191,11 @@ void TrackingGeometry::dumpGeometry(const std::string& outputDir, size_t output_precision = 6; Acts::ObjVisualization3D obj_vis(output_precision, output_scalor); - Acts::ViewConfig container_view = Acts::ViewConfig({220, 220, 220}); - Acts::ViewConfig volume_view = Acts::ViewConfig({220, 220, 0}); - Acts::ViewConfig sensitive_view = Acts::ViewConfig({0, 180, 240}); - Acts::ViewConfig passive_view = Acts::ViewConfig({240, 280, 0}); - Acts::ViewConfig grid_view = Acts::ViewConfig({220, 0, 0}); + Acts::ViewConfig container_view{.color = {220, 220, 220}}; + Acts::ViewConfig volume_view{.color = {220, 220, 0}}; + Acts::ViewConfig sensitive_view{.color = {0, 180, 240}}; + Acts::ViewConfig passive_view{.color = {240, 280, 0}}; + Acts::ViewConfig grid_view{.color = {220, 0, 0}}; Acts::GeometryView3D::drawTrackingVolume( obj_vis, *(t_geometry_->highestTrackingVolume()), gctx, container_view,