algorithms/tracking/TrackSeeding.cc

0001 // SPDX-License-Identifier: LGPL-3.0-or-later
0002 // Copyright (C) 2023  - 2025 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/EventData/Seed.hpp>
0009 #include <Acts/EventData/SpacePointProxy.hpp>
0010 #include <Acts/Seeding/SeedConfirmationRangeConfig.hpp>
0011 #include <Acts/Seeding/SeedFilter.hpp>
0012 #include <Acts/Seeding/SeedFilterConfig.hpp>
0013 #include <Acts/Seeding/SeedFinderConfig.hpp>
0014 #include <Acts/Seeding/SeedFinderOrthogonal.hpp>
0015 #include <Acts/Seeding/SeedFinderOrthogonalConfig.hpp>
0016 #include <Acts/Seeding/SeedFinderUtils.hpp>
0017 #include <Acts/Surfaces/PerigeeSurface.hpp>
0018 #include <Acts/Surfaces/Surface.hpp>
0019 #include <Acts/Utilities/KDTree.hpp> // IWYU pragma: keep FIXME KDTree missing in SeedFinderOrthogonal.hpp until Acts v23.0.0
0020 #include <Acts/Utilities/Result.hpp>
0021 #include <edm4eic/Cov6f.h>
0022 #include <edm4eic/EDM4eicVersion.h>
0023 #include <edm4hep/Vector2f.h>
0024 #include <edm4hep/Vector3f.h>
0025 #include <Eigen/Core>
0026 #include <Eigen/Geometry>
0027 #include <array>
0028 #include <cmath>
0029 #include <limits>
0030 #include <tuple>
0031
0032 namespace eicrecon {
0033
0034 void TrackSeeding::init() {
0035
0036   // Filter parameters
0037   m_seedFilterConfig.maxSeedsPerSpM        = m_cfg.maxSeedsPerSpM_filter;
0038   m_seedFilterConfig.deltaRMin             = m_cfg.deltaRMin;
0039   m_seedFilterConfig.seedConfirmation      = m_cfg.seedConfirmation;
0040   m_seedFilterConfig.deltaInvHelixDiameter = m_cfg.deltaInvHelixDiameter;
0041   m_seedFilterConfig.impactWeightFactor    = m_cfg.impactWeightFactor;
0042   m_seedFilterConfig.zOriginWeightFactor   = m_cfg.zOriginWeightFactor;
0043   m_seedFilterConfig.compatSeedWeight      = m_cfg.compatSeedWeight;
0044   m_seedFilterConfig.compatSeedLimit       = m_cfg.compatSeedLimit;
0045   m_seedFilterConfig.seedWeightIncrement   = m_cfg.seedWeightIncrement;
0046
0047   m_seedFilterConfig.centralSeedConfirmationRange = Acts::SeedConfirmationRangeConfig{
0048       .zMinSeedConf            = m_cfg.zMinSeedConfCentral,
0049       .zMaxSeedConf            = m_cfg.zMaxSeedConfCentral,
0050       .rMaxSeedConf            = m_cfg.rMaxSeedConfCentral,
0051       .nTopForLargeR           = m_cfg.nTopForLargeRCentral,
0052       .nTopForSmallR           = m_cfg.nTopForSmallRCentral,
0053       .seedConfMinBottomRadius = m_cfg.seedConfMinBottomRadiusCentral,
0054       .seedConfMaxZOrigin      = m_cfg.seedConfMaxZOriginCentral,
0055       .minImpactSeedConf       = m_cfg.minImpactSeedConfCentral};
0056
0057   m_seedFilterConfig.forwardSeedConfirmationRange = Acts::SeedConfirmationRangeConfig{
0058       .zMinSeedConf            = m_cfg.zMinSeedConfForward,
0059       .zMaxSeedConf            = m_cfg.zMaxSeedConfForward,
0060       .rMaxSeedConf            = m_cfg.rMaxSeedConfForward,
0061       .nTopForLargeR           = m_cfg.nTopForLargeRForward,
0062       .nTopForSmallR           = m_cfg.nTopForSmallRForward,
0063       .seedConfMinBottomRadius = m_cfg.seedConfMinBottomRadiusForward,
0064       .seedConfMaxZOrigin      = m_cfg.seedConfMaxZOriginForward,
0065       .minImpactSeedConf       = m_cfg.minImpactSeedConfForward};
0066
0067 #if Acts_VERSION_MAJOR < 42
0068   m_seedFilterConfig = m_seedFilterConfig.toInternalUnits();
0069 #endif
0070
0071   // Finder parameters
0072   m_seedFinderConfig.seedFilter =
0073       std::make_unique<Acts::SeedFilter<proxy_type>>(m_seedFilterConfig);
0074   m_seedFinderConfig.rMax               = m_cfg.rMax;
0075   m_seedFinderConfig.rMin               = m_cfg.rMin;
0076   m_seedFinderConfig.deltaRMinTopSP     = m_cfg.deltaRMinTopSP;
0077   m_seedFinderConfig.deltaRMaxTopSP     = m_cfg.deltaRMaxTopSP;
0078   m_seedFinderConfig.deltaRMinBottomSP  = m_cfg.deltaRMinBottomSP;
0079   m_seedFinderConfig.deltaRMaxBottomSP  = m_cfg.deltaRMaxBottomSP;
0080   m_seedFinderConfig.collisionRegionMin = m_cfg.collisionRegionMin;
0081   m_seedFinderConfig.collisionRegionMax = m_cfg.collisionRegionMax;
0082   m_seedFinderConfig.zMin               = m_cfg.zMin;
0083   m_seedFinderConfig.zMax               = m_cfg.zMax;
0084   m_seedFinderConfig.maxSeedsPerSpM     = m_cfg.maxSeedsPerSpM;
0085   m_seedFinderConfig.cotThetaMax        = m_cfg.cotThetaMax;
0086   m_seedFinderConfig.sigmaScattering    = m_cfg.sigmaScattering;
0087   m_seedFinderConfig.radLengthPerSeed   = m_cfg.radLengthPerSeed;
0088   m_seedFinderConfig.minPt              = m_cfg.minPt;
0089   m_seedFinderConfig.impactMax          = m_cfg.impactMax;
0090   m_seedFinderConfig.rMinMiddle         = m_cfg.rMinMiddle;
0091   m_seedFinderConfig.rMaxMiddle         = m_cfg.rMaxMiddle;
0092   m_seedFinderConfig.deltaPhiMax        = m_cfg.deltaPhiMax;
0093
0094   m_seedFinderOptions.beamPos   = Acts::Vector2(m_cfg.beamPosX, m_cfg.beamPosY);
0095   m_seedFinderOptions.bFieldInZ = m_cfg.bFieldInZ;
0096
0097   m_seedFinderConfig = m_seedFinderConfig
0098 #if Acts_VERSION_MAJOR < 42
0099                            .toInternalUnits()
0100 #endif
0101                            .calculateDerivedQuantities();
0102   m_seedFinderOptions = m_seedFinderOptions
0103 #if Acts_VERSION_MAJOR < 42
0104                             .toInternalUnits()
0105 #endif
0106                             .calculateDerivedQuantities(m_seedFinderConfig);
0107 }
0108
0109 void TrackSeeding::process(const Input& input, const Output& output) const {
0110
0111   const auto [trk_hits]        = input;
0112   auto [trk_seeds, trk_params] = output;
0113
0114   std::vector<const eicrecon::SpacePoint*> spacePoints = getSpacePoints(*trk_hits);
0115
0116   Acts::SeedFinderOrthogonal<proxy_type> finder(m_seedFinderConfig); // FIXME move into class scope
0117
0118   // Config
0119   Acts::SpacePointContainerConfig spConfig;
0120
0121   // Options
0122   Acts::SpacePointContainerOptions spOptions;
0123   spOptions.beamPos = {0., 0.};
0124
0125   SpacePointContainerType container(spacePoints);
0126   Acts::SpacePointContainer<decltype(container), Acts::detail::RefHolder> spContainer(
0127       spConfig, spOptions, container);
0128
0129   std::vector<Acts::Seed<proxy_type>> seeds = finder.createSeeds(m_seedFinderOptions, spContainer);
0130
0131   // need to convert here from seed of proxies to seed of sps
0132   for (const auto& seed : seeds) {
0133     const auto& sps = seed.sp();
0134
0135     auto seedToAdd = Acts::Seed<eicrecon::SpacePoint>(*sps[0]->externalSpacePoint(),
0136                                                       *sps[1]->externalSpacePoint(),
0137                                                       *sps[2]->externalSpacePoint());
0138     seedToAdd.setVertexZ(seed.z());
0139     seedToAdd.setQuality(seed.seedQuality());
0140
0141     // Estimate track parameters
0142     auto trackParams = estimateTrackParamsFromSeed(seedToAdd);
0143     if (!trackParams.has_value()) {
0144       debug("Failed to estimate track parameters from seed");
0145       continue;
0146     }
0147     trk_params->push_back(trackParams.value());
0148
0149     // Add seed to collection
0150     auto trk_seed = trk_seeds->create();
0151     trk_seed.setPerigee({0.f, 0.f, 0.f});
0152 #if EDM4EIC_VERSION_MAJOR > 8 || (EDM4EIC_VERSION_MAJOR == 8 && EDM4EIC_VERSION_MINOR > 5)
0153     trk_seed.setQuality(seedToAdd.seedQuality());
0154 #endif
0155     trk_seed.setParams(trackParams.value());
0156     trk_seed.addToHits(*sps[0]->externalSpacePoint());
0157     trk_seed.addToHits(*sps[1]->externalSpacePoint());
0158     trk_seed.addToHits(*sps[2]->externalSpacePoint());
0159   }
0160
0161   for (auto& sp : spacePoints) {
0162     delete sp;
0163   }
0164 }
0165
0166 std::vector<const eicrecon::SpacePoint*>
0167 TrackSeeding::getSpacePoints(const edm4eic::TrackerHitCollection& trk_hits) {
0168   std::vector<const eicrecon::SpacePoint*> spacepoints;
0169
0170   for (const auto hit : trk_hits) {
0171     const eicrecon::SpacePoint* sp = new SpacePoint(hit);
0172     spacepoints.push_back(sp);
0173   }
0174
0175   return spacepoints;
0176 }
0177
0178 std::optional<edm4eic::MutableTrackParameters>
0179 TrackSeeding::estimateTrackParamsFromSeed(const Acts::Seed<SpacePoint>& seed) const {
0180   std::vector<std::pair<float, float>> xyHitPositions;
0181   std::vector<std::pair<float, float>> rzHitPositions;
0182   for (const auto& spptr : seed.sp()) {
0183     xyHitPositions.emplace_back(spptr->x(), spptr->y());
0184     rzHitPositions.emplace_back(spptr->r(), spptr->z());
0185   }
0186
0187   auto RX0Y0 = circleFit(xyHitPositions);
0188   float R    = std::get<0>(RX0Y0);
0189   float X0   = std::get<1>(RX0Y0);
0190   float Y0   = std::get<2>(RX0Y0);
0191   if (!(std::isfinite(R) && std::isfinite(std::abs(X0)) && std::isfinite(std::abs(Y0)))) {
0192     // avoid float overflow for hits on a line
0193     return {};
0194   }
0195   if (std::hypot(X0, Y0) < std::numeric_limits<decltype(std::hypot(X0, Y0))>::epsilon() ||
0196       !std::isfinite(std::hypot(X0, Y0))) {
0197     //Avoid center of circle at origin, where there is no point-of-closest approach
0198     //Also, avoid float overflow on circle center
0199     return {};
0200   }
0201
0202   auto slopeZ0     = lineFit(rzHitPositions);
0203   const auto xypos = findPCA(RX0Y0);
0204
0205   //Determine charge
0206   int charge = determineCharge(xyHitPositions, xypos, RX0Y0);
0207
0208   float theta = atan(1. / std::get<0>(slopeZ0));
0209   // normalize to 0<theta<pi
0210   if (theta < 0) {
0211     theta += M_PI;
0212   }
0213   float eta    = -log(tan(theta / 2.));
0214   float pt     = R * m_cfg.bFieldInZ; // pt[GeV] = R[mm] * B[GeV/mm]
0215   float p      = pt * cosh(eta);
0216   float qOverP = charge / p;
0217
0218   //Calculate phi at xypos
0219   auto xpos = xypos.first;
0220   auto ypos = xypos.second;
0221
0222   auto vxpos = -1. * charge * (ypos - Y0);
0223   auto vypos = charge * (xpos - X0);
0224
0225   auto phi = atan2(vypos, vxpos);
0226
0227   const float z0 = seed.z();
0228   auto perigee   = Acts::Surface::makeShared<Acts::PerigeeSurface>(Acts::Vector3(0, 0, 0));
0229   Acts::Vector3 global(xypos.first, xypos.second, z0);
0230
0231   //Compute local position at PCA
0232   Acts::Vector2 localpos;
0233   Acts::Vector3 direction(sin(theta) * cos(phi), sin(theta) * sin(phi), cos(theta));
0234
0235   auto local = perigee->globalToLocal(m_geoSvc->getActsGeometryContext(), global, direction);
0236
0237   if (!local.ok()) {
0238     return {};
0239   }
0240
0241   localpos = local.value();
0242
0243   auto trackparam = edm4eic::MutableTrackParameters();
0244   trackparam.setType(-1); // type --> seed(-1)
0245   trackparam.setLoc(
0246       {static_cast<float>(localpos(0)), static_cast<float>(localpos(1))}); // 2d location on surface
0247   trackparam.setPhi(static_cast<float>(phi));                              // phi [rad]
0248   trackparam.setTheta(theta);                                              //theta [rad]
0249   trackparam.setQOverP(qOverP);                                            // Q/p [e/GeV]
0250   trackparam.setTime(10);                                                  // time in ns
0251   edm4eic::Cov6f cov;
0252   cov(0, 0) = m_cfg.locaError / Acts::UnitConstants::mm;    // loc0
0253   cov(1, 1) = m_cfg.locbError / Acts::UnitConstants::mm;    // loc1
0254   cov(2, 2) = m_cfg.phiError / Acts::UnitConstants::rad;    // phi
0255   cov(3, 3) = m_cfg.thetaError / Acts::UnitConstants::rad;  // theta
0256   cov(4, 4) = m_cfg.qOverPError * Acts::UnitConstants::GeV; // qOverP
0257   cov(5, 5) = m_cfg.timeError / Acts::UnitConstants::ns;    // time
0258   trackparam.setCovariance(cov);
0259
0260   return trackparam;
0261 }
0262
0263 std::pair<float, float> TrackSeeding::findPCA(std::tuple<float, float, float>& circleParams) {
0264   const float R  = std::get<0>(circleParams);
0265   const float X0 = std::get<1>(circleParams);
0266   const float Y0 = std::get<2>(circleParams);
0267
0268   const double R0 = std::hypot(X0, Y0);
0269
0270   //Calculate point on circle closest to origin
0271   const double xmin = X0 * (1. - R / R0);
0272   const double ymin = Y0 * (1. - R / R0);
0273
0274   return std::make_pair(xmin, ymin);
0275 }
0276
0277 int TrackSeeding::determineCharge(std::vector<std::pair<float, float>>& positions,
0278                                   const std::pair<float, float>& PCA,
0279                                   std::tuple<float, float, float>& RX0Y0) {
0280
0281   const auto& firstpos = positions.at(0);
0282   auto hit_x           = firstpos.first;
0283   auto hit_y           = firstpos.second;
0284
0285   auto xpos = PCA.first;
0286   auto ypos = PCA.second;
0287
0288   float X0 = std::get<1>(RX0Y0);
0289   float Y0 = std::get<2>(RX0Y0);
0290
0291   Acts::Vector3 B_z(0, 0, 1);
0292   Acts::Vector3 radial(X0 - xpos, Y0 - ypos, 0);
0293   Acts::Vector3 hit(hit_x - xpos, hit_y - ypos, 0);
0294
0295   auto cross = radial.cross(hit);
0296
0297   float dot = cross.dot(B_z);
0298
0299   return copysign(1., -dot);
0300 }
0301
0302 /**
0303    * Circle fit to a given set of data points (in 2D)
0304    * This is an algebraic fit, due to Taubin, based on the journal article
0305    * G. Taubin, "Estimation Of Planar Curves, Surfaces And Nonplanar
0306    * Space Curves Defined By Implicit Equations, With
0307    * Applications To Edge And Range Image Segmentation",
0308    * IEEE Trans. PAMI, Vol. 13, pages 1115-1138, (1991)
0309    * It works well whether data points are sampled along an entire circle or along a small arc.
0310    * It still has a small bias and its statistical accuracy is slightly lower than that of the geometric fit (minimizing geometric distances),
0311    * It provides a very good initial guess for a subsequent geometric fit.
0312    * Nikolai Chernov  (September 2012)
0313    */
0314 std::tuple<float, float, float>
0315 TrackSeeding::circleFit(std::vector<std::pair<float, float>>& positions) {
0316   // Compute x- and y- sample means
0317   double meanX  = 0;
0318   double meanY  = 0;
0319   double weight = 0;
0320
0321   for (const auto& [x, y] : positions) {
0322     meanX += x;
0323     meanY += y;
0324     ++weight;
0325   }
0326   meanX /= weight;
0327   meanY /= weight;
0328
0329   //     computing moments
0330
0331   double Mxx = 0;
0332   double Myy = 0;
0333   double Mxy = 0;
0334   double Mxz = 0;
0335   double Myz = 0;
0336   double Mzz = 0;
0337
0338   for (auto& [x, y] : positions) {
0339     double Xi = x - meanX; //  centered x-coordinates
0340     double Yi = y - meanY; //  centered y-coordinates
0341     double Zi = std::pow(Xi, 2) + std::pow(Yi, 2);
0342
0343     Mxy += Xi * Yi;
0344     Mxx += Xi * Xi;
0345     Myy += Yi * Yi;
0346     Mxz += Xi * Zi;
0347     Myz += Yi * Zi;
0348     Mzz += Zi * Zi;
0349   }
0350   Mxx /= weight;
0351   Myy /= weight;
0352   Mxy /= weight;
0353   Mxz /= weight;
0354   Myz /= weight;
0355   Mzz /= weight;
0356
0357   //  computing coefficients of the characteristic polynomial
0358
0359   const double Mz     = Mxx + Myy;
0360   const double Cov_xy = Mxx * Myy - Mxy * Mxy;
0361   const double Var_z  = Mzz - Mz * Mz;
0362   const double A3     = 4 * Mz;
0363   const double A2     = -3 * Mz * Mz - Mzz;
0364   const double A1     = Var_z * Mz + 4 * Cov_xy * Mz - Mxz * Mxz - Myz * Myz;
0365   const double A0  = Mxz * (Mxz * Myy - Myz * Mxy) + Myz * (Myz * Mxx - Mxz * Mxy) - Var_z * Cov_xy;
0366   const double A22 = A2 + A2;
0367   const double A33 = A3 + A3 + A3;
0368
0369   //    finding the root of the characteristic polynomial
0370   //    using Newton's method starting at x=0
0371   //    (it is guaranteed to converge to the right root)
0372   static constexpr int iter_max = 99;
0373   double x                      = 0;
0374   double y                      = A0;
0375
0376   // usually, 4-6 iterations are enough
0377   for (int iter = 0; iter < iter_max; ++iter) {
0378     const double Dy   = A1 + x * (A22 + A33 * x);
0379     const double xnew = x - y / Dy;
0380     if ((xnew == x) || (!std::isfinite(xnew))) {
0381       break;
0382     }
0383
0384     const double ynew = A0 + xnew * (A1 + xnew * (A2 + xnew * A3));
0385     if (std::abs(ynew) >= std::abs(y)) {
0386       break;
0387     }
0388
0389     x = xnew;
0390     y = ynew;
0391   }
0392
0393   //  computing parameters of the fitting circle
0394   const double DET     = std::pow(x, 2) - x * Mz + Cov_xy;
0395   const double Xcenter = (Mxz * (Myy - x) - Myz * Mxy) / DET / 2;
0396   const double Ycenter = (Myz * (Mxx - x) - Mxz * Mxy) / DET / 2;
0397
0398   //  assembling the output
0399   float X0 = Xcenter + meanX;
0400   float Y0 = Ycenter + meanY;
0401   float R  = std::sqrt(std::pow(Xcenter, 2) + std::pow(Ycenter, 2) + Mz);
0402   return std::make_tuple(R, X0, Y0);
0403 }
0404
0405 std::tuple<float, float> TrackSeeding::lineFit(std::vector<std::pair<float, float>>& positions) {
0406   double xsum  = 0;
0407   double x2sum = 0;
0408   double ysum  = 0;
0409   double xysum = 0;
0410   for (const auto& [r, z] : positions) {
0411     xsum  = xsum + r;               //calculate sigma(xi)
0412     ysum  = ysum + z;               //calculate sigma(yi)
0413     x2sum = x2sum + std::pow(r, 2); //calculate sigma(x^2i)
0414     xysum = xysum + r * z;          //calculate sigma(xi*yi)
0415   }
0416
0417   const auto npts          = positions.size();
0418   const double denominator = (x2sum * npts - std::pow(xsum, 2));
0419   const float a            = (xysum * npts - xsum * ysum) / denominator;  //calculate slope
0420   const float b            = (x2sum * ysum - xsum * xysum) / denominator; //calculate intercept
0421   return std::make_tuple(a, b);
0422 }
0423
0424 } // namespace eicrecon