File indexing completed on 2025-07-11 07:53:38
0001
0002
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
0038 template <class T> inline constexpr T square(const T& x) { return x * x; }
0039 }
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
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
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);
0125 #else
0126 Acts::SeedFinderOrthogonal<eicrecon::SpacePoint> finder(
0127 m_seedFinderConfig);
0128 #endif
0129
0130 #if Acts_VERSION_MAJOR >= 37
0131
0132 Acts::SpacePointContainerConfig spConfig;
0133
0134
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
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
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
0216
0217 continue;
0218 }
0219
0220 auto slopeZ0 = lineFit(rzHitPositions);
0221 const auto xypos = findPCA(RX0Y0);
0222
0223
0224 int charge = determineCharge(xyHitPositions, xypos, RX0Y0);
0225
0226 float theta = atan(1. / std::get<0>(slopeZ0));
0227
0228 if (theta < 0) {
0229 theta += M_PI;
0230 }
0231 float eta = -log(tan(theta / 2.));
0232 float pt = R * m_cfg.bFieldInZ;
0233 float p = pt * cosh(eta);
0234 float qOverP = charge / p;
0235
0236
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
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);
0263 trackparam.setLoc({static_cast<float>(localpos(0)),
0264 static_cast<float>(localpos(1))});
0265 trackparam.setPhi(static_cast<float>(phi));
0266 trackparam.setTheta(theta);
0267 trackparam.setQOverP(qOverP);
0268 trackparam.setTime(10);
0269 edm4eic::Cov6f cov;
0270 cov(0, 0) = m_cfg.locaError / Acts::UnitConstants::mm;
0271 cov(1, 1) = m_cfg.locbError / Acts::UnitConstants::mm;
0272 cov(2, 2) = m_cfg.phiError / Acts::UnitConstants::rad;
0273 cov(3, 3) = m_cfg.thetaError / Acts::UnitConstants::rad;
0274 cov(4, 4) = m_cfg.qOverPError * Acts::UnitConstants::GeV;
0275 cov(5, 5) = m_cfg.timeError / Acts::UnitConstants::ns;
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
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
0323
0324
0325
0326
0327
0328
0329
0330
0331
0332
0333 std::tuple<float, float, float>
0334 eicrecon::TrackSeeding::circleFit(std::vector<std::pair<float, float>>& positions) {
0335
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
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;
0359 double Yi = y - meanY;
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
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
0389
0390
0391 static constexpr int iter_max = 99;
0392 double x = 0;
0393 double y = A0;
0394
0395
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
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
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;
0432 ysum = ysum + z;
0433 x2sum = x2sum + std::pow(r, 2);
0434 xysum = xysum + r * z;
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;
0440 const float b = (x2sum * ysum - xsum * xysum) / denominator;
0441 return std::make_tuple(a, b);
0442 }