algorithms/tracking/TrackSeeding.cc

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