File indexing completed on 2024-09-27 07:03:02
0001
0002
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
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
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
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);
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
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
0193
0194 continue;
0195 }
0196
0197 auto slopeZ0 = lineFit(rzHitPositions);
0198 const auto xypos = findPCA(RX0Y0);
0199
0200
0201
0202 int charge = determineCharge(xyHitPositions, xypos, RX0Y0);
0203
0204 float theta = atan(1./std::get<0>(slopeZ0));
0205
0206 if(theta < 0)
0207 { theta += M_PI; }
0208 float eta = -log(tan(theta/2.));
0209 float pt = R * m_cfg.bFieldInZ;
0210 float p = pt * cosh(eta);
0211 float qOverP = charge / p;
0212
0213
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
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);
0243 trackparam.setLoc({static_cast<float>(localpos(0)), static_cast<float>(localpos(1))});
0244 trackparam.setPhi(static_cast<float>(phi));
0245 trackparam.setTheta(theta);
0246 trackparam.setQOverP(qOverP);
0247 trackparam.setTime(10);
0248 edm4eic::Cov6f cov;
0249 cov(0,0) = m_cfg.locaError / Acts::UnitConstants::mm;
0250 cov(1,1) = m_cfg.locbError / Acts::UnitConstants::mm;
0251 cov(2,2) = m_cfg.phiError / Acts::UnitConstants::rad;
0252 cov(3,3) = m_cfg.thetaError / Acts::UnitConstants::rad;
0253 cov(4,4) = m_cfg.qOverPError * Acts::UnitConstants::GeV;
0254 cov(5,5) = m_cfg.timeError / Acts::UnitConstants::ns;
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
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
0302
0303
0304
0305
0306
0307
0308
0309
0310
0311
0312 std::tuple<float,float,float> eicrecon::TrackSeeding::circleFit(std::vector<std::pair<float,float>>& positions) const
0313 {
0314
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
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;
0340 double Yi = y - meanY;
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
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
0370
0371
0372 static constexpr int iter_max = 99;
0373 double x = 0;
0374 double y = A0;
0375
0376
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
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
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;
0411 ysum=ysum+z;
0412 x2sum=x2sum+std::pow(r,2);
0413 xysum=xysum+r*z;
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;
0419 const float b= (x2sum*ysum-xsum*xysum)/denominator;
0420 return std::make_tuple( a, b );
0421 }