Back to home page

EIC code displayed by LXR
 
 

    

File indexing completed on 2026-10-11 08:14:32 

0001 // SPDX-License-Identifier: LGPL-3.0-or-later
0002 // Copyright (C) 2023 - 2026, Joe Osborn, Dmitry Romanov, Wouter Deconinck
0003 
0004 #include "TrackSeeding.h"
0005 
0006 #include <Acts/Definitions/Algebra.hpp>
0007 #include <Acts/Definitions/Units.hpp>
0008 #include <Acts/Seeding/SeedConfirmationRangeConfig.hpp>
0009 #include <Acts/Surfaces/PerigeeSurface.hpp>
0010 #include <Acts/Surfaces/Surface.hpp>
0011 #include <Acts/Utilities/Result.hpp>
0012 #include <edm4eic/Cov6f.h>
0013 #include <edm4eic/CovDiag3f.h>
0014 #include <edm4eic/EDM4eicVersion.h>
0015 #include <edm4hep/Vector2f.h>
0016 #include <edm4hep/Vector3f.h>
0017 #include <spdlog/common.h>
0018 #include <spdlog/logger.h>
0019 #include <spdlog/spdlog.h>
0020 #include <Eigen/Core>
0021 #include <Eigen/Geometry>
0022 #include <array>
0023 #include <cmath>
0024 #include <cstdint>
0025 #include <limits>
0026 #include <numbers>
0027 #include <span>
0028 #include <tuple>
0029 
0030 #include "extensions/spdlog/SpdlogToActs.h"
0031 
0032 // Acts version-specific includes
0033 #if TRACKSEEDING_HAS_SEEDING
0034 // Modern Seeding API includes
0035 #include <Acts/Definitions/Direction.hpp>
0036 #include <Acts/EventData/SeedContainer.hpp>
0037 #include <Acts/EventData/SeedProxy.hpp>
0038 #include <Acts/EventData/SpacePointColumns.hpp>
0039 #include <Acts/EventData/SpacePointContainer.hpp>
0040 #include <Acts/EventData/SpacePointProxy.hpp>
0041 #include <Acts/EventData/Types.hpp>
0042 #include <Acts/Geometry/Extent.hpp>
0043 #include <Acts/Seeding/BroadTripletSeedFilter.hpp>
0044 #include <Acts/Seeding/CylindricalSpacePointKDTree.hpp>
0045 #include <Acts/Seeding/DoubletSeedFinder.hpp>
0046 #include <Acts/Seeding/TripletSeedFinder.hpp>
0047 #include <Acts/Seeding/TripletSeeder.hpp>
0048 #include <Acts/Utilities/AxisDefinitions.hpp>
0049 #include <Acts/Utilities/Logger.hpp>
0050 
0051 namespace Acts {
0052 using SeedContainer2                             = Acts::SeedContainer;
0053 using SpacePointIndex2                           = Acts::SpacePointIndex;
0054 template <bool read_only> using SeedProxy2       = Acts::SeedProxy<read_only>;
0055 using SpacePointContainer2                       = Acts::SpacePointContainer;
0056 template <bool read_only> using SpacePointProxy2 = Acts::SpacePointProxy<read_only>;
0057 } // namespace Acts
0058 #endif
0059 
0060 #if TRACKSEEDING_HAS_SEEDING2
0061 // Modern Seeding2 API includes
0062 #include <Acts/Definitions/Direction.hpp>
0063 #include <Acts/EventData/SeedContainer2.hpp>
0064 #include <Acts/EventData/SeedProxy2.hpp>
0065 #include <Acts/EventData/SpacePointColumns.hpp>
0066 #include <Acts/EventData/SpacePointContainer2.hpp>
0067 #include <Acts/EventData/SpacePointProxy2.hpp>
0068 #include <Acts/EventData/Types.hpp>
0069 #include <Acts/Geometry/Extent.hpp>
0070 #include <Acts/Seeding2/BroadTripletSeedFilter.hpp>
0071 #include <Acts/Seeding2/CylindricalSpacePointKDTree.hpp>
0072 #include <Acts/Seeding2/DoubletSeedFinder.hpp>
0073 #include <Acts/Seeding2/TripletSeedFinder.hpp>
0074 #include <Acts/Seeding2/TripletSeeder.hpp>
0075 #include <Acts/Utilities/AxisDefinitions.hpp>
0076 #include <Acts/Utilities/Logger.hpp>
0077 #endif
0078 
0079 #if TRACKSEEDING_HAS_ORTHOGONAL
0080 // Orthogonal API includes
0081 #include <Acts/EventData/Seed.hpp>
0082 #include <Acts/EventData/SpacePointProxy.hpp>
0083 #include <Acts/Seeding/SeedFilter.hpp>
0084 #include <Acts/Seeding/SeedFilterConfig.hpp>
0085 #include <Acts/Seeding/SeedFinderConfig.hpp>
0086 #include <Acts/Seeding/SeedFinderOrthogonal.hpp>
0087 #include <Acts/Seeding/SeedFinderOrthogonalConfig.hpp>
0088 #include <Acts/Seeding/SeedFinderUtils.hpp>
0089 #endif
0090 
0091 namespace eicrecon {
0092 
0093 void TrackSeeding::init() {
0094   // Step 1: Resolve Auto to specific method based on Acts version
0095   m_resolvedMethod = m_cfg.seedingMethod;
0096   if (m_resolvedMethod == TrackSeedingConfig::SeedingMethod::Auto) {
0097 #if TRACKSEEDING_HAS_SEEDING2 || TRACKSEEDING_HAS_SEEDING
0098     // Prefer Seeding2 (modern API) when available
0099     m_resolvedMethod = TrackSeedingConfig::SeedingMethod::Seeding2;
0100 #elif TRACKSEEDING_HAS_ORTHOGONAL
0101     m_resolvedMethod = TrackSeedingConfig::SeedingMethod::Orthogonal;
0102 #else
0103 #error "No seeding method available - check Acts version compatibility"
0104 #endif
0105   }
0106 
0107   // Step 2: Validate method availability
0108 #if !(TRACKSEEDING_HAS_SEEDING2 || TRACKSEEDING_HAS_SEEDING)
0109   if (m_resolvedMethod == TrackSeedingConfig::SeedingMethod::Seeding2) {
0110     throw std::runtime_error("TrackSeeding: Seeding2 method not available in Acts " +
0111                              std::to_string(Acts_VERSION_MAJOR) + "." +
0112                              std::to_string(Acts_VERSION_MINOR) +
0113                              ". Use seedingMethod='auto' or 'orthogonal'.");
0114   }
0115 #endif
0116 #if !TRACKSEEDING_HAS_ORTHOGONAL
0117   if (m_resolvedMethod == TrackSeedingConfig::SeedingMethod::Orthogonal) {
0118     throw std::runtime_error("TrackSeeding: Orthogonal method not available in Acts " +
0119                              std::to_string(Acts_VERSION_MAJOR) + "." +
0120                              std::to_string(Acts_VERSION_MINOR) +
0121                              ". Use seedingMethod='auto' or 'seeding2'.");
0122   }
0123 #endif
0124 
0125   // Step 3: Initialize based on resolved method
0126 #if TRACKSEEDING_HAS_SEEDING2 || TRACKSEEDING_HAS_SEEDING
0127   if (m_resolvedMethod == TrackSeedingConfig::SeedingMethod::Seeding2) {
0128     // Initialize Seeding2
0129 #if TRACKSEEDING_HAS_SEEDING2 && TRACKSEEDING_HAS_ORTHOGONAL
0130     m_seedingData = SeedingData{};
0131     auto& data    = std::get<SeedingData>(m_seedingData);
0132 #else
0133     auto& data = m_seedingData;
0134 #endif
0135 
0136     // Create spdlog logger with matching log level for Acts
0137     auto spdlog_logger = spdlog::default_logger()->clone(std::string(this->name()));
0138     spdlog_logger->set_level(static_cast<spdlog::level::level_enum>(this->level()));
0139 
0140     data.actsLogger = eicrecon::getSpdlogLogger(std::string(this->name()), spdlog_logger);
0141 
0142     data.filterConfig.deltaInvHelixDiameter = m_cfg.deltaInvHelixDiameter;
0143     data.filterConfig.deltaRMin             = m_cfg.deltaRMin;
0144     data.filterConfig.compatSeedWeight      = m_cfg.compatSeedWeight;
0145     data.filterConfig.impactWeightFactor    = m_cfg.impactWeightFactor;
0146     data.filterConfig.zOriginWeightFactor   = m_cfg.zOriginWeightFactor;
0147     data.filterConfig.maxSeedsPerSpM        = m_cfg.maxSeedsPerSpM;
0148     data.filterConfig.compatSeedLimit       = m_cfg.compatSeedLimit;
0149     data.filterConfig.seedWeightIncrement   = m_cfg.seedWeightIncrement;
0150     data.filterConfig.seedConfirmation      = m_cfg.seedConfirmation;
0151 
0152     data.filterConfig.centralSeedConfirmationRange = Acts::SeedConfirmationRangeConfig{
0153         .zMinSeedConf            = m_cfg.zMinSeedConfCentral,
0154         .zMaxSeedConf            = m_cfg.zMaxSeedConfCentral,
0155         .rMaxSeedConf            = m_cfg.rMaxSeedConfCentral,
0156         .nTopForLargeR           = m_cfg.nTopForLargeRCentral,
0157         .nTopForSmallR           = m_cfg.nTopForSmallRCentral,
0158         .seedConfMinBottomRadius = m_cfg.seedConfMinBottomRadiusCentral,
0159         .seedConfMaxZOrigin      = m_cfg.seedConfMaxZOriginCentral,
0160         .minImpactSeedConf       = m_cfg.minImpactSeedConfCentral};
0161 
0162     data.filterConfig.forwardSeedConfirmationRange = Acts::SeedConfirmationRangeConfig{
0163         .zMinSeedConf            = m_cfg.zMinSeedConfForward,
0164         .zMaxSeedConf            = m_cfg.zMaxSeedConfForward,
0165         .rMaxSeedConf            = m_cfg.rMaxSeedConfForward,
0166         .nTopForLargeR           = m_cfg.nTopForLargeRForward,
0167         .nTopForSmallR           = m_cfg.nTopForSmallRForward,
0168         .seedConfMinBottomRadius = m_cfg.seedConfMinBottomRadiusForward,
0169         .seedConfMaxZOrigin      = m_cfg.seedConfMaxZOriginForward,
0170         .minImpactSeedConf       = m_cfg.minImpactSeedConfForward};
0171 
0172     data.seedFinder.emplace();
0173 
0174     // Configure bottom doublet finder
0175     Acts::DoubletSeedFinder::Config bottomDoubletFinderConfig;
0176     bottomDoubletFinderConfig.spacePointsSortedByRadius = false;
0177     bottomDoubletFinderConfig.candidateDirection        = Acts::Direction::Backward();
0178     bottomDoubletFinderConfig.deltaRMin                 = m_cfg.deltaRMinBottomSP;
0179     bottomDoubletFinderConfig.deltaRMax                 = m_cfg.deltaRMaxBottomSP;
0180     bottomDoubletFinderConfig.deltaZMin                 = m_cfg.deltaZMin;
0181     bottomDoubletFinderConfig.deltaZMax                 = m_cfg.deltaZMax;
0182     bottomDoubletFinderConfig.impactMax                 = m_cfg.impactMax;
0183     bottomDoubletFinderConfig.interactionPointCut       = m_cfg.interactionPointCut;
0184     bottomDoubletFinderConfig.collisionRegionMin        = m_cfg.collisionRegionMin;
0185     bottomDoubletFinderConfig.collisionRegionMax        = m_cfg.collisionRegionMax;
0186     bottomDoubletFinderConfig.cotThetaMax               = m_cfg.cotThetaMax;
0187     bottomDoubletFinderConfig.minPt                     = m_cfg.minPt;
0188     bottomDoubletFinderConfig.helixCutTolerance         = m_cfg.helixCutTolerance;
0189     data.bottomDoubletFinder                            = Acts::DoubletSeedFinder::create(
0190         Acts::DoubletSeedFinder::DerivedConfig(bottomDoubletFinderConfig, m_cfg.bFieldInZ));
0191 
0192     // Configure top doublet finder
0193     Acts::DoubletSeedFinder::Config topDoubletFinderConfig = bottomDoubletFinderConfig;
0194     topDoubletFinderConfig.candidateDirection              = Acts::Direction::Forward();
0195     topDoubletFinderConfig.deltaRMin                       = m_cfg.deltaRMinTopSP;
0196     topDoubletFinderConfig.deltaRMax                       = m_cfg.deltaRMaxTopSP;
0197     data.topDoubletFinder                                  = Acts::DoubletSeedFinder::create(
0198         Acts::DoubletSeedFinder::DerivedConfig(topDoubletFinderConfig, m_cfg.bFieldInZ));
0199 
0200     // Configure triplet finder
0201     Acts::TripletSeedFinder::Config tripletFinderConfig;
0202     tripletFinderConfig.useStripInfo      = false;
0203     tripletFinderConfig.sortedByCotTheta  = true;
0204     tripletFinderConfig.minPt             = m_cfg.minPt;
0205     tripletFinderConfig.sigmaScattering   = m_cfg.sigmaScattering;
0206     tripletFinderConfig.radLengthPerSeed  = m_cfg.radLengthPerSeed;
0207     tripletFinderConfig.impactMax         = m_cfg.impactMax;
0208     tripletFinderConfig.helixCutTolerance = m_cfg.helixCutTolerance;
0209     data.tripletFinder                    = Acts::TripletSeedFinder::create(
0210         Acts::TripletSeedFinder::DerivedConfig(tripletFinderConfig, m_cfg.bFieldInZ));
0211   }
0212 #endif
0213 
0214 #if TRACKSEEDING_HAS_ORTHOGONAL
0215   if (m_resolvedMethod == TrackSeedingConfig::SeedingMethod::Orthogonal) {
0216     // Initialize Orthogonal
0217 #if TRACKSEEDING_HAS_SEEDING2 && TRACKSEEDING_HAS_ORTHOGONAL
0218     m_seedingData = OrthogonalData{};
0219     auto& data    = std::get<OrthogonalData>(m_seedingData);
0220 #else
0221     auto& data = m_seedingData;
0222 #endif
0223 
0224     // Create spdlog logger with matching log level for Acts
0225     auto spdlog_logger = spdlog::default_logger()->clone(std::string(this->name()));
0226     spdlog_logger->set_level(static_cast<spdlog::level::level_enum>(this->level()));
0227 
0228     data.actsLogger = eicrecon::getSpdlogLogger(std::string(this->name()), spdlog_logger);
0229 
0230     data.seedFilterConfig.maxSeedsPerSpM        = m_cfg.maxSeedsPerSpM;
0231     data.seedFilterConfig.deltaRMin             = m_cfg.deltaRMin;
0232     data.seedFilterConfig.seedConfirmation      = m_cfg.seedConfirmation;
0233     data.seedFilterConfig.deltaInvHelixDiameter = m_cfg.deltaInvHelixDiameter;
0234     data.seedFilterConfig.impactWeightFactor    = m_cfg.impactWeightFactor;
0235     data.seedFilterConfig.zOriginWeightFactor   = m_cfg.zOriginWeightFactor;
0236     data.seedFilterConfig.compatSeedWeight      = m_cfg.compatSeedWeight;
0237     data.seedFilterConfig.compatSeedLimit       = m_cfg.compatSeedLimit;
0238     data.seedFilterConfig.seedWeightIncrement   = m_cfg.seedWeightIncrement;
0239 
0240     data.seedFilterConfig.centralSeedConfirmationRange = Acts::SeedConfirmationRangeConfig{
0241         .zMinSeedConf            = m_cfg.zMinSeedConfCentral,
0242         .zMaxSeedConf            = m_cfg.zMaxSeedConfCentral,
0243         .rMaxSeedConf            = m_cfg.rMaxSeedConfCentral,
0244         .nTopForLargeR           = m_cfg.nTopForLargeRCentral,
0245         .nTopForSmallR           = m_cfg.nTopForSmallRCentral,
0246         .seedConfMinBottomRadius = m_cfg.seedConfMinBottomRadiusCentral,
0247         .seedConfMaxZOrigin      = m_cfg.seedConfMaxZOriginCentral,
0248         .minImpactSeedConf       = m_cfg.minImpactSeedConfCentral};
0249 
0250     data.seedFilterConfig.forwardSeedConfirmationRange = Acts::SeedConfirmationRangeConfig{
0251         .zMinSeedConf            = m_cfg.zMinSeedConfForward,
0252         .zMaxSeedConf            = m_cfg.zMaxSeedConfForward,
0253         .rMaxSeedConf            = m_cfg.rMaxSeedConfForward,
0254         .nTopForLargeR           = m_cfg.nTopForLargeRForward,
0255         .nTopForSmallR           = m_cfg.nTopForSmallRForward,
0256         .seedConfMinBottomRadius = m_cfg.seedConfMinBottomRadiusForward,
0257         .seedConfMaxZOrigin      = m_cfg.seedConfMaxZOriginForward,
0258         .minImpactSeedConf       = m_cfg.minImpactSeedConfForward};
0259 
0260     data.seedFinderConfig.seedFilter =
0261         std::make_unique<Acts::SeedFilter<proxy_type>>(data.seedFilterConfig);
0262     data.seedFinderConfig.rMax               = m_cfg.rMax;
0263     data.seedFinderConfig.rMin               = m_cfg.rMin;
0264     data.seedFinderConfig.deltaRMinTopSP     = m_cfg.deltaRMinTopSP;
0265     data.seedFinderConfig.deltaRMaxTopSP     = m_cfg.deltaRMaxTopSP;
0266     data.seedFinderConfig.deltaRMinBottomSP  = m_cfg.deltaRMinBottomSP;
0267     data.seedFinderConfig.deltaRMaxBottomSP  = m_cfg.deltaRMaxBottomSP;
0268     data.seedFinderConfig.collisionRegionMin = m_cfg.collisionRegionMin;
0269     data.seedFinderConfig.collisionRegionMax = m_cfg.collisionRegionMax;
0270     data.seedFinderConfig.zMin               = m_cfg.zMin;
0271     data.seedFinderConfig.zMax               = m_cfg.zMax;
0272     data.seedFinderConfig.maxSeedsPerSpM     = m_cfg.maxSeedsPerSpM;
0273     data.seedFinderConfig.cotThetaMax        = m_cfg.cotThetaMax;
0274     data.seedFinderConfig.sigmaScattering    = m_cfg.sigmaScattering;
0275     data.seedFinderConfig.radLengthPerSeed   = m_cfg.radLengthPerSeed;
0276     data.seedFinderConfig.minPt              = m_cfg.minPt;
0277     data.seedFinderConfig.impactMax          = m_cfg.impactMax;
0278     data.seedFinderConfig.rMinMiddle         = m_cfg.rMinMiddle;
0279     data.seedFinderConfig.rMaxMiddle         = m_cfg.rMaxMiddle;
0280     data.seedFinderConfig.deltaPhiMax        = m_cfg.deltaPhiMax;
0281 
0282     data.seedFinderOptions.beamPos   = Acts::Vector2(m_cfg.beamPosX, m_cfg.beamPosY);
0283     data.seedFinderOptions.bFieldInZ = m_cfg.bFieldInZ;
0284 
0285     data.seedFinderConfig = data.seedFinderConfig.calculateDerivedQuantities();
0286     data.seedFinderOptions =
0287         data.seedFinderOptions.calculateDerivedQuantities(data.seedFinderConfig);
0288   }
0289 #endif
0290 }
0291 
0292 void TrackSeeding::process(const Input& input, const Output& output) const {
0293   const auto [trk_hits]        = input;
0294   auto [trk_seeds, trk_params] = output;
0295 
0296 #if TRACKSEEDING_HAS_SEEDING2 || TRACKSEEDING_HAS_SEEDING
0297   if (m_resolvedMethod == TrackSeedingConfig::SeedingMethod::Seeding2) {
0298     // Get Seeding2 data from variant or direct member
0299 #if TRACKSEEDING_HAS_SEEDING2 && TRACKSEEDING_HAS_ORTHOGONAL
0300     const auto& data = std::get<SeedingData>(m_seedingData);
0301 #else
0302     const auto& data = m_seedingData;
0303 #endif
0304 
0305     // ========== Seeding2 Processing ==========
0306     if (trk_hits->empty()) {
0307       return;
0308     }
0309     // Build SpacePointContainer2 from tracker hits
0310     Acts::SpacePointContainer2 spacePoints(
0311         Acts::SpacePointColumns::PackedXY | Acts::SpacePointColumns::PackedZR |
0312         Acts::SpacePointColumns::Phi | Acts::SpacePointColumns::VarianceZ |
0313         Acts::SpacePointColumns::VarianceR |
0314 #if TRACKSEEDING_HAS_SEEDING
0315         Acts::SpacePointColumns::CopiedFromIndex
0316 #else
0317         Acts::SpacePointColumns::CopyFromIndex
0318 #endif
0319     );
0320     spacePoints.reserve(trk_hits->size());
0321 
0322     Acts::Experimental::CylindricalSpacePointKDTreeBuilder kdTreeBuilder;
0323     kdTreeBuilder.reserve(trk_hits->size());
0324 
0325     Acts::Extent rRangeSPExtent;
0326 
0327     for (std::uint32_t i = 0; i < trk_hits->size(); ++i) {
0328       const auto& hit = (*trk_hits)[i];
0329       const float hx  = hit.getPosition()[0];
0330       const float hy  = hit.getPosition()[1];
0331       const float hz  = hit.getPosition()[2];
0332       const float hr  = std::hypot(hx, hy);
0333 
0334       if (hr < m_cfg.rMin || hr > m_cfg.rMax || hz < m_cfg.zMin || hz > m_cfg.zMax) {
0335         continue;
0336       }
0337 
0338       const float varR =
0339           (hx * hx * hit.getPositionError().xx + hy * hy * hit.getPositionError().yy) /
0340           (hx * hx + hy * hy + std::numeric_limits<float>::epsilon());
0341       const float varZ = hit.getPositionError().zz;
0342 
0343       Acts::SpacePointIndex2 spIdx = spacePoints.size();
0344       auto sp                      = spacePoints.createSpacePoint();
0345       sp.xy()                      = {hx, hy};
0346       sp.zr()                      = {hz, hr};
0347       sp.phi()                     = std::atan2(hy, hx);
0348       sp.varianceZ()               = varZ;
0349       sp.varianceR()               = varR;
0350 #if TRACKSEEDING_HAS_SEEDING
0351       sp.copiedFromIndex() = i;
0352 #else
0353       sp.copyFromIndex() = i;
0354 #endif
0355 
0356       kdTreeBuilder.insert(spIdx, sp.phi(), hr, hz);
0357       rRangeSPExtent.extend({hx, hy, hz});
0358     }
0359     if (spacePoints.empty()) {
0360       return;
0361     }
0362 
0363     // Build KD-tree
0364     Acts::Experimental::CylindricalSpacePointKDTree kdTree = kdTreeBuilder.build();
0365 
0366     // Configure KD-tree options for bottom and top candidate searches
0367     Acts::Experimental::CylindricalSpacePointKDTree::Options bottomOptions;
0368     bottomOptions.rMax   = m_cfg.rMax;
0369     bottomOptions.zMin   = m_cfg.zMin;
0370     bottomOptions.zMax   = m_cfg.zMax;
0371     bottomOptions.phiMin = m_cfg.phiMin;
0372     bottomOptions.phiMax = m_cfg.phiMax;
0373     // Use bottom-specific deltaR parameters for bottom region
0374     bottomOptions.deltaRMin          = m_cfg.deltaRMinBottomSP;
0375     bottomOptions.deltaRMax          = m_cfg.deltaRMaxBottomSP;
0376     bottomOptions.collisionRegionMin = m_cfg.collisionRegionMin;
0377     bottomOptions.collisionRegionMax = m_cfg.collisionRegionMax;
0378     bottomOptions.cotThetaMax        = m_cfg.cotThetaMax;
0379     bottomOptions.deltaPhiMax        = m_cfg.deltaPhiMax;
0380 
0381     Acts::Experimental::CylindricalSpacePointKDTree::Options topOptions = bottomOptions;
0382     // Use top-specific deltaR parameters for top region
0383     topOptions.deltaRMin = m_cfg.deltaRMinTopSP;
0384     topOptions.deltaRMax = m_cfg.deltaRMaxTopSP;
0385 
0386     // Configure doublet finders
0387     // (constructed once in init() and stored in data)
0388 
0389     const bool useVariableMiddleSPRange = m_cfg.useVariableMiddleSPRange;
0390     float rMiddleSPMin                  = 0.F;
0391     float rMiddleSPMax                  = 0.F;
0392     if (useVariableMiddleSPRange) {
0393       rMiddleSPMin = std::floor(rRangeSPExtent.min(Acts::AxisDirection::AxisR) / 2) * 2 +
0394                      m_cfg.deltaRMiddleMinSPRange;
0395       rMiddleSPMax = std::floor(rRangeSPExtent.max(Acts::AxisDirection::AxisR) / 2) * 2 -
0396                      m_cfg.deltaRMiddleMaxSPRange;
0397     }
0398 
0399     // Setup seed filter
0400     Acts::BroadTripletSeedFilter::State filterState;
0401     Acts::BroadTripletSeedFilter::Cache filterCache;
0402     Acts::BroadTripletSeedFilter seedFilter(data.filterConfig, filterState, filterCache,
0403                                             *data.actsLogger);
0404 
0405     // Thread-local storage for candidates and cache
0406     thread_local Acts::TripletSeeder::Cache seederCache;
0407     thread_local Acts::Experimental::CylindricalSpacePointKDTree::Candidates candidates;
0408 
0409     // Output seed container
0410     Acts::SeedContainer2 actsSeeds;
0411     actsSeeds.assignSpacePointContainer(spacePoints);
0412 
0413     // Run the seeding algorithm by iterating over all points in the tree
0414     for (const auto& middle : kdTree) {
0415       const auto spM = spacePoints.at(middle.second).asConst();
0416 
0417       // Cut: Ensure middle space point lies within valid r-region
0418       const float rM = spM.zr()[1];
0419       if (useVariableMiddleSPRange) {
0420         if (rM < rMiddleSPMin || rM > rMiddleSPMax) {
0421           continue;
0422         }
0423       } else {
0424         if (rM > m_cfg.rMaxMiddle || rM < m_cfg.rMinMiddle) {
0425           continue;
0426         }
0427       }
0428 
0429       // Remove middle SPs outside phi and z region of interest
0430       const float zM = spM.zr()[0];
0431       if (zM < m_cfg.zOutermostLayers.first || zM > m_cfg.zOutermostLayers.second) {
0432         continue;
0433       }
0434       if (const float phiM = spM.phi(); phiM > m_cfg.phiMax || phiM < m_cfg.phiMin) {
0435         continue;
0436       }
0437 
0438       // Determine number of top seeds for confirmation
0439       std::size_t nTopSeedConf = 0;
0440       if (m_cfg.seedConfirmation) {
0441         // Choose central or forward region based on z of middle space point
0442         bool useCentralRegion =
0443             (zM <= m_cfg.zMaxSeedConfCentral && zM >= m_cfg.zMinSeedConfCentral);
0444         float rMaxSeedConf =
0445             useCentralRegion ? m_cfg.rMaxSeedConfCentral : m_cfg.rMaxSeedConfForward;
0446         std::size_t nTopForLargeR =
0447             useCentralRegion ? m_cfg.nTopForLargeRCentral : m_cfg.nTopForLargeRForward;
0448         std::size_t nTopForSmallR =
0449             useCentralRegion ? m_cfg.nTopForSmallRCentral : m_cfg.nTopForSmallRForward;
0450         nTopSeedConf = rM > rMaxSeedConf ? nTopForLargeR : nTopForSmallR;
0451       }
0452 
0453       // Find valid tuples of (bottom, middle, top) candidates
0454       candidates.clear();
0455       kdTree.validTuples(topOptions, bottomOptions, spM, nTopSeedConf, candidates);
0456 
0457       // Process bottom-low-high and top-low-high combinations
0458       Acts::SpacePointContainer2::ConstSubset bottomSps =
0459           spacePoints.subset(candidates.bottom_lh_v).asConst();
0460       Acts::SpacePointContainer2::ConstSubset topSps =
0461           spacePoints.subset(candidates.top_lh_v).asConst();
0462       data.seedFinder->createSeedsFromGroup(seederCache, *data.bottomDoubletFinder,
0463                                             *data.topDoubletFinder, *data.tripletFinder, seedFilter,
0464                                             spacePoints, bottomSps, spM, topSps, actsSeeds);
0465 
0466       bottomSps = spacePoints.subset(candidates.bottom_hl_v).asConst();
0467       topSps    = spacePoints.subset(candidates.top_hl_v).asConst();
0468       data.seedFinder->createSeedsFromGroup(seederCache, *data.bottomDoubletFinder,
0469                                             *data.topDoubletFinder, *data.tripletFinder, seedFilter,
0470                                             spacePoints, bottomSps, spM, topSps, actsSeeds);
0471     }
0472 
0473     debug("Created {} track seeds from {} space points", actsSeeds.size(), spacePoints.size());
0474 
0475     // Convert Acts seeds to EDM4eic format
0476     for (const auto& actsSeed : actsSeeds) {
0477       const auto& spIndices = actsSeed.spacePointIndices();
0478 
0479       // Get the three space points (bottom, middle, top)
0480       std::array<std::array<float, 3>, 3> positions;
0481       for (std::size_t k = 0; k < 3; ++k) {
0482 #if TRACKSEEDING_HAS_SEEDING
0483         const std::uint32_t hitIdx = spacePoints.at(spIndices[k]).copiedFromIndex();
0484 #else
0485         const std::uint32_t hitIdx = spacePoints.at(spIndices[k]).copyFromIndex();
0486 #endif
0487         const auto hit = (*trk_hits)[hitIdx];
0488         positions[k]   = {hit.getPosition()[0], hit.getPosition()[1], hit.getPosition()[2]};
0489       }
0490 
0491       // Estimate track parameters from seed
0492       auto trackParams =
0493           estimateTrackParamsFromSeed(positions, actsSeed.vertexZ(), m_cfg.beamPosX, m_cfg.beamPosY,
0494                                       m_cfg.bFieldInZ, m_geoSvc, m_cfg);
0495       if (!trackParams.has_value()) {
0496         debug("Failed to estimate track parameters from seed");
0497         continue;
0498       }
0499       trk_params->push_back(trackParams.value());
0500 
0501       // Add seed to collection
0502       auto trk_seed = trk_seeds->create();
0503       trk_seed.setPerigee({0.F, 0.F, 0.F});
0504 #if EDM4EIC_VERSION_MAJOR > 8 || (EDM4EIC_VERSION_MAJOR == 8 && EDM4EIC_VERSION_MINOR > 5)
0505       trk_seed.setQuality(actsSeed.quality());
0506 #endif
0507       trk_seed.setParams(trackParams.value());
0508       for (std::size_t k = 0; k < 3; ++k) {
0509 #if TRACKSEEDING_HAS_SEEDING
0510         const std::uint32_t hitIdx = spacePoints.at(spIndices[k]).copiedFromIndex();
0511 #else
0512         const std::uint32_t hitIdx = spacePoints.at(spIndices[k]).copyFromIndex();
0513 #endif
0514         trk_seed.addToHits((*trk_hits)[hitIdx]);
0515       }
0516     }
0517   }
0518 #endif
0519 
0520 #if TRACKSEEDING_HAS_ORTHOGONAL
0521   if (m_resolvedMethod == TrackSeedingConfig::SeedingMethod::Orthogonal) {
0522     // Get Orthogonal data from variant or direct member
0523 #if TRACKSEEDING_HAS_SEEDING2 && TRACKSEEDING_HAS_ORTHOGONAL
0524     const auto& data = std::get<OrthogonalData>(m_seedingData);
0525 #else
0526     const auto& data = m_seedingData;
0527 #endif
0528 
0529     // ========== Orthogonal Processing ==========
0530     std::vector<const eicrecon::SpacePoint*> spacePoints = getSpacePoints(*trk_hits);
0531 
0532     Acts::SeedFinderOrthogonal<proxy_type> finder(data.seedFinderConfig, data.actsLogger->clone());
0533 
0534     Acts::SpacePointContainerConfig spConfig;
0535     Acts::SpacePointContainerOptions spOptions;
0536     spOptions.beamPos = {m_cfg.beamPosX, m_cfg.beamPosY};
0537 
0538     SpacePointContainerType container(spacePoints);
0539     Acts::SpacePointContainer<decltype(container), Acts::detail::RefHolder> spContainer(
0540         spConfig, spOptions, container);
0541 
0542     std::vector<Acts::Seed<proxy_type>> seeds =
0543         finder.createSeeds(data.seedFinderOptions, spContainer);
0544 
0545     // need to convert here from seed of proxies to seed of sps
0546     for (const auto& seed : seeds) {
0547       const auto& sps = seed.sp();
0548 
0549       auto seedToAdd = Acts::Seed<eicrecon::SpacePoint>(*sps[0]->externalSpacePoint(),
0550                                                         *sps[1]->externalSpacePoint(),
0551                                                         *sps[2]->externalSpacePoint());
0552       seedToAdd.setVertexZ(seed.z());
0553       seedToAdd.setQuality(seed.seedQuality());
0554 
0555       // Estimate track parameters
0556       auto trackParams = estimateTrackParamsFromSeed(seedToAdd);
0557       if (!trackParams.has_value()) {
0558         debug("Failed to estimate track parameters from seed");
0559         continue;
0560       }
0561       trk_params->push_back(trackParams.value());
0562 
0563       // Add seed to collection
0564       auto trk_seed = trk_seeds->create();
0565       trk_seed.setPerigee({0.F, 0.F, 0.F});
0566 #if EDM4EIC_VERSION_MAJOR > 8 || (EDM4EIC_VERSION_MAJOR == 8 && EDM4EIC_VERSION_MINOR > 5)
0567       trk_seed.setQuality(seedToAdd.seedQuality());
0568 #endif
0569       trk_seed.setParams(trackParams.value());
0570       trk_seed.addToHits(*sps[0]->externalSpacePoint());
0571       trk_seed.addToHits(*sps[1]->externalSpacePoint());
0572       trk_seed.addToHits(*sps[2]->externalSpacePoint());
0573     }
0574 
0575     for (auto& sp : spacePoints) {
0576       delete sp;
0577     }
0578   }
0579 #endif
0580 }
0581 
0582 // Shared core physics calculation for track parameter estimation
0583 std::optional<edm4eic::MutableTrackParameters> TrackSeeding::computeTrackParametersFromFit(
0584     const std::vector<std::pair<float, float>>& xyPositions,
0585     const std::vector<std::pair<float, float>>& rzPositions, float vertexZ, float bFieldInZ,
0586     const std::shared_ptr<const ActsGeometryProvider>& geoSvc, const TrackSeedingConfig& cfg) {
0587   // Make mutable copies for fitting functions
0588   auto xyPosCopy = xyPositions;
0589   auto rzPosCopy = rzPositions;
0590 
0591   auto RX0Y0 = circleFit(xyPosCopy);
0592   float R    = std::get<0>(RX0Y0);
0593   float X0   = std::get<1>(RX0Y0);
0594   float Y0   = std::get<2>(RX0Y0);
0595   if (!(std::isfinite(R) && std::isfinite(std::abs(X0)) && std::isfinite(std::abs(Y0)))) {
0596     // avoid float overflow for hits on a line
0597     return {};
0598   }
0599   if (std::hypot(X0, Y0) < std::numeric_limits<decltype(std::hypot(X0, Y0))>::epsilon() ||
0600       !std::isfinite(std::hypot(X0, Y0))) {
0601     return {};
0602   }
0603 
0604   auto slopeZ0     = lineFit(rzPosCopy);
0605   const auto xypos = findPCA(RX0Y0);
0606 
0607   // Determine charge
0608   int charge = determineCharge(xyPosCopy, xypos, RX0Y0);
0609 
0610   float theta = std::atan(1.F / std::get<0>(slopeZ0));
0611   // normalize to 0<theta<pi
0612   if (theta < 0) {
0613     theta += std::numbers::pi_v<float>;
0614   }
0615   float eta    = -std::log(std::tan(theta / 2.F));
0616   float pt     = R * bFieldInZ; // pt[GeV] = R[mm] * B[GeV/mm]
0617   float p      = pt * std::cosh(eta);
0618   float qOverP = static_cast<float>(charge) / p;
0619 
0620   // Calculate phi at xypos
0621   auto xpos  = xypos.first;
0622   auto ypos  = xypos.second;
0623   auto vxpos = -1. * charge * (ypos - Y0);
0624   auto vypos = charge * (xpos - X0);
0625   auto phi   = std::atan2(vypos, vxpos);
0626 
0627   auto perigee = Acts::Surface::makeShared<Acts::PerigeeSurface>(Acts::Vector3(0, 0, 0));
0628   Acts::Vector3 global(xypos.first, xypos.second, vertexZ);
0629 
0630   // Compute local position at PCA
0631   Acts::Vector3 direction(std::sin(theta) * std::cos(phi), std::sin(theta) * std::sin(phi),
0632                           std::cos(theta));
0633   auto local = perigee->globalToLocal(geoSvc->getActsGeometryContext(), global, direction);
0634   if (!local.ok()) {
0635     return {};
0636   }
0637   Acts::Vector2 localpos = local.value();
0638 
0639   auto trackparam = edm4eic::MutableTrackParameters();
0640   trackparam.setType(-1); // type --> seed(-1)
0641   trackparam.setLoc(
0642       {static_cast<float>(localpos(0)), static_cast<float>(localpos(1))}); // 2d location on surface
0643   trackparam.setPhi(static_cast<float>(phi));                              // phi [rad]
0644   trackparam.setTheta(theta);                                              // theta [rad]
0645   trackparam.setQOverP(qOverP);                                            // Q/p [e/GeV]
0646   trackparam.setTime(10);                                                  // time in ns
0647   edm4eic::Cov6f cov;
0648   cov(0, 0) = std::pow(cfg.locaError / Acts::UnitConstants::mm, 2);   // loc0 variance in mm^2
0649   cov(1, 1) = std::pow(cfg.locbError / Acts::UnitConstants::mm, 2);   // loc1 variance in mm^2
0650   cov(2, 2) = std::pow(cfg.phiError / Acts::UnitConstants::rad, 2);   // phi variance in rad^2
0651   cov(3, 3) = std::pow(cfg.thetaError / Acts::UnitConstants::rad, 2); // theta variance in rad^2
0652   cov(4, 4) =
0653       std::pow(cfg.qOverPError * Acts::UnitConstants::GeV, 2);      // qOverP variance in (e/GeV)^2
0654   cov(5, 5) = std::pow(cfg.timeError / Acts::UnitConstants::ns, 2); // time variance in ns^2
0655   trackparam.setCovariance(cov);
0656   return trackparam;
0657 }
0658 
0659 #if TRACKSEEDING_HAS_SEEDING2 || TRACKSEEDING_HAS_SEEDING
0660 std::optional<edm4eic::MutableTrackParameters> TrackSeeding::estimateTrackParamsFromSeed(
0661     const std::array<std::array<float, 3>, 3>& spPositions, float vertexZ,
0662     float beamPosX [[maybe_unused]], float beamPosY [[maybe_unused]], float bFieldInZ,
0663     const std::shared_ptr<const ActsGeometryProvider>& geoSvc, const TrackSeedingConfig& cfg) {
0664   // Note: beamPosX/beamPosY not currently used in track parameter estimation
0665   // but passed through for potential future improvements
0666   std::vector<std::pair<float, float>> xyPositions;
0667   std::vector<std::pair<float, float>> rzPositions;
0668   xyPositions.reserve(3);
0669   rzPositions.reserve(3);
0670   for (const auto& pos : spPositions) {
0671     xyPositions.emplace_back(pos[0], pos[1]);
0672     rzPositions.emplace_back(std::hypot(pos[0], pos[1]), pos[2]);
0673   }
0674 
0675   return computeTrackParametersFromFit(xyPositions, rzPositions, vertexZ, bFieldInZ, geoSvc, cfg);
0676 }
0677 #endif
0678 
0679 #if TRACKSEEDING_HAS_ORTHOGONAL
0680 std::vector<const eicrecon::SpacePoint*>
0681 TrackSeeding::getSpacePoints(const edm4eic::TrackerHitCollection& trk_hits) {
0682   std::vector<const eicrecon::SpacePoint*> spacepoints;
0683 
0684   for (const auto hit : trk_hits) {
0685     const eicrecon::SpacePoint* sp = new SpacePoint(hit);
0686     spacepoints.push_back(sp);
0687   }
0688 
0689   return spacepoints;
0690 }
0691 
0692 std::optional<edm4eic::MutableTrackParameters>
0693 TrackSeeding::estimateTrackParamsFromSeed(const Acts::Seed<SpacePoint>& seed) const {
0694   std::vector<std::pair<float, float>> xyHitPositions;
0695   std::vector<std::pair<float, float>> rzHitPositions;
0696   for (const auto& spptr : seed.sp()) {
0697     xyHitPositions.emplace_back(spptr->x(), spptr->y());
0698     rzHitPositions.emplace_back(spptr->r(), spptr->z());
0699   }
0700 
0701   return computeTrackParametersFromFit(xyHitPositions, rzHitPositions, seed.z(), m_cfg.bFieldInZ,
0702                                        m_geoSvc, m_cfg);
0703 }
0704 #endif
0705 
0706 std::pair<float, float> TrackSeeding::findPCA(std::tuple<float, float, float>& circleParams) {
0707   const float R  = std::get<0>(circleParams);
0708   const float X0 = std::get<1>(circleParams);
0709   const float Y0 = std::get<2>(circleParams);
0710 
0711   // Calculate point on circle closest to origin
0712   const double R0 = std::hypot(X0, Y0);
0713 
0714   const double xmin = X0 * (1. - R / R0);
0715   const double ymin = Y0 * (1. - R / R0);
0716 
0717   return std::make_pair(xmin, ymin);
0718 }
0719 
0720 int TrackSeeding::determineCharge(std::vector<std::pair<float, float>>& positions,
0721                                   const std::pair<float, float>& PCA,
0722                                   std::tuple<float, float, float>& RX0Y0) {
0723   const auto& firstpos = positions.at(0);
0724   auto hit_x           = firstpos.first;
0725   auto hit_y           = firstpos.second;
0726 
0727   auto xpos = PCA.first;
0728   auto ypos = PCA.second;
0729 
0730   float X0 = std::get<1>(RX0Y0);
0731   float Y0 = std::get<2>(RX0Y0);
0732 
0733   Acts::Vector3 B_z(0, 0, 1);
0734   Acts::Vector3 radial(X0 - xpos, Y0 - ypos, 0);
0735   Acts::Vector3 hit(hit_x - xpos, hit_y - ypos, 0);
0736 
0737   auto cross = radial.cross(hit);
0738 
0739   float dot = cross.dot(B_z);
0740 
0741   return copysign(1., -dot);
0742 }
0743 
0744 /**
0745    * Circle fit to a given set of data points (in 2D)
0746    * This is an algebraic fit, due to Taubin, based on the journal article
0747    * G. Taubin, "Estimation Of Planar Curves, Surfaces And Nonplanar
0748    * Space Curves Defined By Implicit Equations, With
0749    * Applications To Edge And Range Image Segmentation",
0750    * IEEE Trans. PAMI, Vol. 13, pages 1115-1138, (1991)
0751    * It works well whether data points are sampled along an entire circle or along a small arc.
0752    * It still has a small bias and its statistical accuracy is slightly lower than that of the geometric fit (minimizing geometric distances),
0753    * It provides a very good initial guess for a subsequent geometric fit.
0754    * Nikolai Chernov  (September 2012)
0755    */
0756 std::tuple<float, float, float>
0757 TrackSeeding::circleFit(std::vector<std::pair<float, float>>& positions) {
0758   // Compute x- and y- sample means
0759   double meanX  = 0;
0760   double meanY  = 0;
0761   double weight = 0;
0762 
0763   for (const auto& [x, y] : positions) {
0764     meanX += x;
0765     meanY += y;
0766     ++weight;
0767   }
0768   meanX /= weight;
0769   meanY /= weight;
0770 
0771   //     computing moments
0772 
0773   double Mxx = 0;
0774   double Myy = 0;
0775   double Mxy = 0;
0776   double Mxz = 0;
0777   double Myz = 0;
0778   double Mzz = 0;
0779 
0780   for (auto& [x, y] : positions) {
0781     double Xi = x - meanX; //  centered x-coordinates
0782     double Yi = y - meanY; //  centered y-coordinates
0783     double Zi = std::pow(Xi, 2) + std::pow(Yi, 2);
0784 
0785     Mxy += Xi * Yi;
0786     Mxx += Xi * Xi;
0787     Myy += Yi * Yi;
0788     Mxz += Xi * Zi;
0789     Myz += Yi * Zi;
0790     Mzz += Zi * Zi;
0791   }
0792   Mxx /= weight;
0793   Myy /= weight;
0794   Mxy /= weight;
0795   Mxz /= weight;
0796   Myz /= weight;
0797   Mzz /= weight;
0798 
0799   //  computing coefficients of the characteristic polynomial
0800   const double Mz     = Mxx + Myy;
0801   const double Cov_xy = Mxx * Myy - Mxy * Mxy;
0802   const double Var_z  = Mzz - Mz * Mz;
0803   const double A3     = 4 * Mz;
0804   const double A2     = -3 * Mz * Mz - Mzz;
0805   const double A1     = Var_z * Mz + 4 * Cov_xy * Mz - Mxz * Mxz - Myz * Myz;
0806   const double A0  = Mxz * (Mxz * Myy - Myz * Mxy) + Myz * (Myz * Mxx - Mxz * Mxy) - Var_z * Cov_xy;
0807   const double A22 = A2 + A2;
0808   const double A33 = A3 + A3 + A3;
0809 
0810   //    finding the root of the characteristic polynomial
0811   //    using Newton's method starting at x=0
0812   //    (it is guaranteed to converge to the right root)
0813   static constexpr int iter_max = 99;
0814   double x                      = 0;
0815   double y                      = A0;
0816 
0817   // usually, 4-6 iterations are enough
0818   for (int iter = 0; iter < iter_max; ++iter) {
0819     const double Dy   = A1 + x * (A22 + A33 * x);
0820     const double xnew = x - y / Dy;
0821     if ((xnew == x) || (!std::isfinite(xnew))) {
0822       break;
0823     }
0824 
0825     const double ynew = A0 + xnew * (A1 + xnew * (A2 + xnew * A3));
0826     if (std::abs(ynew) >= std::abs(y)) {
0827       break;
0828     }
0829 
0830     x = xnew;
0831     y = ynew;
0832   }
0833 
0834   //  computing parameters of the fitting circle
0835   const double DET     = std::pow(x, 2) - x * Mz + Cov_xy;
0836   const double Xcenter = (Mxz * (Myy - x) - Myz * Mxy) / DET / 2;
0837   const double Ycenter = (Myz * (Mxx - x) - Mxz * Mxy) / DET / 2;
0838 
0839   float X0 = Xcenter + meanX;
0840   float Y0 = Ycenter + meanY;
0841   float R  = std::sqrt(std::pow(Xcenter, 2) + std::pow(Ycenter, 2) + Mz);
0842 
0843   //  assembling the output
0844   return std::make_tuple(R, X0, Y0);
0845 }
0846 
0847 std::tuple<float, float> TrackSeeding::lineFit(std::vector<std::pair<float, float>>& positions) {
0848   double xsum  = 0;
0849   double x2sum = 0;
0850   double ysum  = 0;
0851   double xysum = 0;
0852   for (const auto& [r, z] : positions) {
0853     xsum  = xsum + r;
0854     ysum  = ysum + z;
0855     x2sum = x2sum + std::pow(r, 2);
0856     xysum = xysum + static_cast<double>(r) * static_cast<double>(z);
0857   }
0858 
0859   const auto npts          = positions.size();
0860   const double denominator = (x2sum * npts - std::pow(xsum, 2));
0861   const float a            = (xysum * npts - xsum * ysum) / denominator;
0862   const float b            = (x2sum * ysum - xsum * xysum) / denominator;
0863   return std::make_tuple(a, b);
0864 }
0865 
0866 } // namespace eicrecon
[ Source navigation ] [ Diff markup ] [ Identifier search ] [ general search ]

This page was automatically generated by the 2.3.7 LXR engine.

The LXR team
Valid CSS 2.1! Valid HTML 4.01!