Back to home page

EIC code displayed by LXR

 
 

    


File indexing completed on 2025-07-11 07:53:38

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