File indexing completed on 2025-08-28 08:12:01
0001
0002
0003
0004
0005
0006
0007
0008
0009 #include "Acts/Seeding2/TripletSeedFinder.hpp"
0010
0011 #include "Acts/EventData/SpacePointContainer2.hpp"
0012 #include "Acts/Utilities/MathHelpers.hpp"
0013 #include "Acts/Utilities/Zip.hpp"
0014
0015 #include <ranges>
0016
0017 #include <Eigen/Dense>
0018 #include <boost/mp11.hpp>
0019 #include <boost/mp11/algorithm.hpp>
0020
0021 namespace Acts::Experimental {
0022
0023 namespace {
0024
0025
0026
0027 bool stripCoordinateCheck(float tolerance, const ConstSpacePointProxy2& sp,
0028 const std::array<float, 3>& spacePointPosition,
0029 std::array<float, 3>& outputCoordinates) {
0030 const std::array<float, 3>& topStripVector = sp.topStripVector();
0031 const std::array<float, 3>& bottomStripVector = sp.bottomStripVector();
0032 const std::array<float, 3>& stripCenterDistance = sp.stripCenterDistance();
0033
0034
0035 const std::array<float, 3> d1 = {
0036 topStripVector[1] * spacePointPosition[2] -
0037 topStripVector[2] * spacePointPosition[1],
0038 topStripVector[2] * spacePointPosition[0] -
0039 topStripVector[0] * spacePointPosition[2],
0040 topStripVector[0] * spacePointPosition[1] -
0041 topStripVector[1] * spacePointPosition[0]};
0042
0043
0044 const float bd1 = bottomStripVector[0] * d1[0] +
0045 bottomStripVector[1] * d1[1] + bottomStripVector[2] * d1[2];
0046
0047
0048
0049 const float s1 = stripCenterDistance[0] * d1[0] +
0050 stripCenterDistance[1] * d1[1] +
0051 stripCenterDistance[2] * d1[2];
0052 if (std::abs(s1) > std::abs(bd1) * tolerance) {
0053 return false;
0054 }
0055
0056
0057 const std::array<float, 3> d0 = {
0058 bottomStripVector[1] * spacePointPosition[2] -
0059 bottomStripVector[2] * spacePointPosition[1],
0060 bottomStripVector[2] * spacePointPosition[0] -
0061 bottomStripVector[0] * spacePointPosition[2],
0062 bottomStripVector[0] * spacePointPosition[1] -
0063 bottomStripVector[1] * spacePointPosition[0]};
0064
0065
0066
0067 float s0 = stripCenterDistance[0] * d0[0] + stripCenterDistance[1] * d0[1] +
0068 stripCenterDistance[2] * d0[2];
0069 if (std::abs(s0) > std::abs(bd1) * tolerance) {
0070 return false;
0071 }
0072
0073
0074
0075
0076 const std::array<float, 3>& topStripCenter = sp.topStripCenter();
0077
0078
0079
0080 s0 = s0 / bd1;
0081 outputCoordinates[0] = topStripCenter[0] + topStripVector[0] * s0;
0082 outputCoordinates[1] = topStripCenter[1] + topStripVector[1] * s0;
0083 outputCoordinates[2] = topStripCenter[2] + topStripVector[2] * s0;
0084 return true;
0085 }
0086
0087 template <bool useStripInfo, bool sortedByCotTheta>
0088 class Impl final : public TripletSeedFinder {
0089 public:
0090 explicit Impl(const DerivedConfig& config) : m_cfg(config) {}
0091
0092 const DerivedConfig& config() const override { return m_cfg; }
0093
0094 template <typename TopDoublets>
0095 void createPixelTripletTopCandidates(
0096 const ConstSpacePointProxy2& spM,
0097 const DoubletsForMiddleSp::Proxy& bottomDoublet, TopDoublets& topDoublets,
0098 TripletTopCandidates& tripletTopCandidates) const {
0099 const float rM = spM.zr()[1];
0100 const float varianceZM = spM.varianceZ();
0101 const float varianceRM = spM.varianceR();
0102
0103
0104 tripletTopCandidates.reserve(tripletTopCandidates.size() +
0105 topDoublets.size());
0106
0107 const float cotThetaB = bottomDoublet.cotTheta();
0108 const float erB = bottomDoublet.er();
0109 const float iDeltaRB = bottomDoublet.iDeltaR();
0110 const float Ub = bottomDoublet.u();
0111 const float Vb = bottomDoublet.v();
0112
0113
0114 const float iSinTheta2 = 1 + cotThetaB * cotThetaB;
0115 const float sigmaSquaredPtDependent = iSinTheta2 * m_cfg.sigmapT2perRadius;
0116
0117
0118
0119
0120
0121
0122
0123
0124
0125 const float scatteringInRegion2 = m_cfg.multipleScattering2 * iSinTheta2;
0126
0127 std::size_t topDoubletOffset = 0;
0128 for (auto [topDoublet, topDoubletIndex] :
0129 zip(topDoublets, std::ranges::iota_view<std::size_t, std::size_t>(
0130 0, topDoublets.size()))) {
0131 const SpacePointIndex2 spT = topDoublet.spacePointIndex();
0132 const float cotThetaT = topDoublet.cotTheta();
0133
0134
0135 const float cotThetaAvg2 = cotThetaB * cotThetaT;
0136
0137
0138
0139 const float error2 = topDoublet.er() + erB +
0140 2 * (cotThetaAvg2 * varianceRM + varianceZM) *
0141 iDeltaRB * topDoublet.iDeltaR();
0142
0143 const float deltaCotTheta = cotThetaB - cotThetaT;
0144 const float deltaCotTheta2 = deltaCotTheta * deltaCotTheta;
0145
0146
0147
0148
0149
0150
0151
0152
0153
0154
0155
0156 if (deltaCotTheta2 > error2 + scatteringInRegion2) {
0157 if constexpr (sortedByCotTheta) {
0158
0159
0160
0161 if (cotThetaB < cotThetaT) {
0162 break;
0163 }
0164 topDoubletOffset = topDoubletIndex + 1;
0165 }
0166 continue;
0167 }
0168
0169 const float dU = topDoublet.u() - Ub;
0170
0171 if (dU == 0) {
0172 continue;
0173 }
0174
0175
0176 const float A = (topDoublet.v() - Vb) / dU;
0177 const float S2 = 1 + A * A;
0178 const float B = Vb - A * Ub;
0179 const float B2 = B * B;
0180
0181
0182
0183 if (S2 < B2 * m_cfg.minHelixDiameter2) {
0184 continue;
0185 }
0186
0187
0188
0189
0190 const float iHelixDiameter2 = B2 / S2;
0191
0192
0193 float p2scatterSigma = iHelixDiameter2 * sigmaSquaredPtDependent;
0194 if (!std::isinf(m_cfg.maxPtScattering)) {
0195
0196
0197
0198
0199 if (B2 == 0 || m_cfg.pTPerHelixRadius * std::sqrt(S2 / B2) >
0200 2 * m_cfg.maxPtScattering) {
0201 const float pTscatterSigma =
0202 (m_cfg.highland / m_cfg.maxPtScattering) * m_cfg.sigmaScattering;
0203 p2scatterSigma = pTscatterSigma * pTscatterSigma * iSinTheta2;
0204 }
0205 }
0206
0207
0208 if (deltaCotTheta2 > error2 + p2scatterSigma) {
0209 if constexpr (sortedByCotTheta) {
0210 if (cotThetaB < cotThetaT) {
0211 break;
0212 }
0213 topDoubletOffset = topDoubletIndex;
0214 }
0215 continue;
0216 }
0217
0218
0219
0220
0221 const float im = std::abs((A - B * rM) * rM);
0222 if (im > m_cfg.impactMax) {
0223 continue;
0224 }
0225
0226
0227
0228 tripletTopCandidates.emplace_back(spT, B / std::sqrt(S2), im);
0229 }
0230
0231 if constexpr (sortedByCotTheta) {
0232
0233 topDoublets = topDoublets.subrange(topDoubletOffset);
0234 }
0235 }
0236
0237 template <typename TopDoublets>
0238 void createStripTripletTopCandidates(
0239 const SpacePointContainer2& spacePoints, const ConstSpacePointProxy2& spM,
0240 const DoubletsForMiddleSp::Proxy& bottomDoublet,
0241 const TopDoublets& topDoublets,
0242 TripletTopCandidates& tripletTopCandidates) const {
0243 const float rM = spM.zr()[1];
0244 const float cosPhiM = spM.xy()[0] / rM;
0245 const float sinPhiM = spM.xy()[1] / rM;
0246 const float varianceZM = spM.varianceZ();
0247 const float varianceRM = spM.varianceR();
0248
0249
0250 tripletTopCandidates.reserve(tripletTopCandidates.size() +
0251 topDoublets.size());
0252
0253 float cotThetaB = bottomDoublet.cotTheta();
0254 const float erB = bottomDoublet.er();
0255 const float iDeltaRB = bottomDoublet.iDeltaR();
0256 const float Vb = bottomDoublet.v();
0257 const float Ub = bottomDoublet.u();
0258
0259
0260 const float iSinTheta2 = 1 + cotThetaB * cotThetaB;
0261 const float sigmaSquaredPtDependent = iSinTheta2 * m_cfg.sigmapT2perRadius;
0262
0263
0264
0265
0266
0267
0268
0269
0270
0271 const float scatteringInRegion2 = m_cfg.multipleScattering2 * iSinTheta2;
0272
0273 const float sinTheta = 1 / std::sqrt(iSinTheta2);
0274 const float cosTheta = cotThetaB * sinTheta;
0275
0276
0277
0278 const std::array<float, 2> rotationTermsUVtoXY = {cosPhiM * sinTheta,
0279 sinPhiM * sinTheta};
0280
0281 for (auto topDoublet : topDoublets) {
0282
0283 float dU = topDoublet.u() - Ub;
0284 if (dU == 0) {
0285 continue;
0286 }
0287
0288
0289 const float A0 = (topDoublet.v() - Vb) / dU;
0290
0291 const float zPositionMiddle = cosTheta * std::sqrt(1 + A0 * A0);
0292
0293
0294
0295 const std::array<float, 3> positionMiddle = {
0296 rotationTermsUVtoXY[0] - rotationTermsUVtoXY[1] * A0,
0297 rotationTermsUVtoXY[0] * A0 + rotationTermsUVtoXY[1],
0298 zPositionMiddle};
0299
0300 std::array<float, 3> rMTransf{};
0301 if (!stripCoordinateCheck(m_cfg.toleranceParam, spM, positionMiddle,
0302 rMTransf)) {
0303 continue;
0304 }
0305
0306
0307 const float B0 = 2 * (Vb - A0 * Ub);
0308 const float Cb = 1 - B0 * bottomDoublet.y();
0309 const float Sb = A0 + B0 * bottomDoublet.x();
0310 const std::array<float, 3> positionBottom = {
0311 rotationTermsUVtoXY[0] * Cb - rotationTermsUVtoXY[1] * Sb,
0312 rotationTermsUVtoXY[0] * Sb + rotationTermsUVtoXY[1] * Cb,
0313 zPositionMiddle};
0314
0315 const ConstSpacePointProxy2 spB =
0316 spacePoints[bottomDoublet.spacePointIndex()];
0317 std::array<float, 3> rBTransf{};
0318 if (!stripCoordinateCheck(m_cfg.toleranceParam, spB, positionBottom,
0319 rBTransf)) {
0320 continue;
0321 }
0322
0323
0324 const float Ct = 1 - B0 * topDoublet.y();
0325 const float St = A0 + B0 * topDoublet.x();
0326 const std::array<float, 3> positionTop = {
0327 rotationTermsUVtoXY[0] * Ct - rotationTermsUVtoXY[1] * St,
0328 rotationTermsUVtoXY[0] * St + rotationTermsUVtoXY[1] * Ct,
0329 zPositionMiddle};
0330
0331 const ConstSpacePointProxy2 spT =
0332 spacePoints[topDoublet.spacePointIndex()];
0333 std::array<float, 3> rTTransf{};
0334 if (!stripCoordinateCheck(m_cfg.toleranceParam, spT, positionTop,
0335 rTTransf)) {
0336 continue;
0337 }
0338
0339
0340 const float xB = rBTransf[0] - rMTransf[0];
0341 const float yB = rBTransf[1] - rMTransf[1];
0342 const float zB = rBTransf[2] - rMTransf[2];
0343 const float xT = rTTransf[0] - rMTransf[0];
0344 const float yT = rTTransf[1] - rMTransf[1];
0345 const float zT = rTTransf[2] - rMTransf[2];
0346
0347 const float iDeltaRB2 = 1 / (xB * xB + yB * yB);
0348 const float iDeltaRT2 = 1 / (xT * xT + yT * yT);
0349
0350 cotThetaB = -zB * std::sqrt(iDeltaRB2);
0351 const float cotThetaT = zT * std::sqrt(iDeltaRT2);
0352
0353
0354 const float averageCotTheta = 0.5f * (cotThetaB + cotThetaT);
0355 const float cotThetaAvg2 = averageCotTheta * averageCotTheta;
0356
0357
0358
0359 const float error2 = topDoublet.er() + erB +
0360 2 * (cotThetaAvg2 * varianceRM + varianceZM) *
0361 iDeltaRB * topDoublet.iDeltaR();
0362
0363 const float deltaCotTheta = cotThetaB - cotThetaT;
0364 const float deltaCotTheta2 = deltaCotTheta * deltaCotTheta;
0365
0366
0367
0368
0369
0370
0371
0372
0373
0374
0375
0376 if (deltaCotTheta2 > error2 + scatteringInRegion2) {
0377
0378 continue;
0379 }
0380
0381 const float rMxy =
0382 std::sqrt(rMTransf[0] * rMTransf[0] + rMTransf[1] * rMTransf[1]);
0383 const float irMxy = 1 / rMxy;
0384 const float Ax = rMTransf[0] * irMxy;
0385 const float Ay = rMTransf[1] * irMxy;
0386
0387 const float ub = (xB * Ax + yB * Ay) * iDeltaRB2;
0388 const float vb = (yB * Ax - xB * Ay) * iDeltaRB2;
0389 const float ut = (xT * Ax + yT * Ay) * iDeltaRT2;
0390 const float vt = (yT * Ax - xT * Ay) * iDeltaRT2;
0391
0392 dU = ut - ub;
0393
0394 if (dU == 0) {
0395 continue;
0396 }
0397 const float A = (vt - vb) / dU;
0398 const float S2 = 1 + A * A;
0399 const float B = vb - A * ub;
0400 const float B2 = B * B;
0401
0402
0403
0404 if (S2 < B2 * m_cfg.minHelixDiameter2) {
0405 continue;
0406 }
0407
0408
0409
0410
0411 const float iHelixDiameter2 = B2 / S2;
0412
0413
0414 float p2scatterSigma = iHelixDiameter2 * sigmaSquaredPtDependent;
0415 if (!std::isinf(m_cfg.maxPtScattering)) {
0416
0417
0418
0419
0420 if (B2 == 0 || m_cfg.pTPerHelixRadius * std::sqrt(S2 / B2) >
0421 2 * m_cfg.maxPtScattering) {
0422 float pTscatterSigma =
0423 (m_cfg.highland / m_cfg.maxPtScattering) * m_cfg.sigmaScattering;
0424 p2scatterSigma = pTscatterSigma * pTscatterSigma * iSinTheta2;
0425 }
0426 }
0427
0428
0429 if (deltaCotTheta2 > error2 + p2scatterSigma) {
0430 continue;
0431 }
0432
0433
0434
0435
0436 const float im = std::abs((A - B * rMxy) * rMxy);
0437 if (im > m_cfg.impactMax) {
0438 continue;
0439 }
0440
0441
0442
0443 tripletTopCandidates.emplace_back(spT.index(), B / std::sqrt(S2), im);
0444 }
0445 }
0446
0447 void createTripletTopCandidates(
0448 const SpacePointContainer2& spacePoints, const ConstSpacePointProxy2& spM,
0449 const DoubletsForMiddleSp::Proxy& bottomDoublet,
0450 DoubletsForMiddleSp::Range& topDoublets,
0451 TripletTopCandidates& tripletTopCandidates) const override {
0452 if constexpr (useStripInfo) {
0453 createStripTripletTopCandidates(spacePoints, spM, bottomDoublet,
0454 topDoublets, tripletTopCandidates);
0455 } else {
0456 createPixelTripletTopCandidates(spM, bottomDoublet, topDoublets,
0457 tripletTopCandidates);
0458 }
0459 }
0460
0461 void createTripletTopCandidates(
0462 const SpacePointContainer2& spacePoints, const ConstSpacePointProxy2& spM,
0463 const DoubletsForMiddleSp::Proxy& bottomDoublet,
0464 DoubletsForMiddleSp::Subset& topDoublets,
0465 TripletTopCandidates& tripletTopCandidates) const override {
0466 if constexpr (useStripInfo) {
0467 createStripTripletTopCandidates(spacePoints, spM, bottomDoublet,
0468 topDoublets, tripletTopCandidates);
0469 } else {
0470 createPixelTripletTopCandidates(spM, bottomDoublet, topDoublets,
0471 tripletTopCandidates);
0472 }
0473 }
0474
0475 void createTripletTopCandidates(
0476 const SpacePointContainer2& spacePoints, const ConstSpacePointProxy2& spM,
0477 const DoubletsForMiddleSp::Proxy& bottomDoublet,
0478 DoubletsForMiddleSp::Subset2& topDoublets,
0479 TripletTopCandidates& tripletTopCandidates) const override {
0480 if constexpr (useStripInfo) {
0481 createStripTripletTopCandidates(spacePoints, spM, bottomDoublet,
0482 topDoublets, tripletTopCandidates);
0483 } else {
0484 createPixelTripletTopCandidates(spM, bottomDoublet, topDoublets,
0485 tripletTopCandidates);
0486 }
0487 }
0488
0489 private:
0490 DerivedConfig m_cfg;
0491 };
0492
0493 }
0494
0495 TripletSeedFinder::DerivedConfig::DerivedConfig(const Config& config,
0496 float bFieldInZ_)
0497 : Config(config), bFieldInZ(bFieldInZ_) {
0498 using namespace Acts::UnitLiterals;
0499
0500
0501 {
0502 const double xOverX0 = radLengthPerSeed;
0503 const double q2OverBeta2 = 1;
0504
0505 const double t = std::sqrt(xOverX0 * q2OverBeta2);
0506
0507
0508 highland =
0509 static_cast<float>(13.6_MeV * t * (1.0 + 0.038 * 2 * std::log(t)));
0510 }
0511
0512 const float maxScatteringAngle = highland / minPt;
0513 const float maxScatteringAngle2 = maxScatteringAngle * maxScatteringAngle;
0514
0515
0516 pTPerHelixRadius = bFieldInZ;
0517 minHelixDiameter2 = square(minPt * 2 / pTPerHelixRadius) * helixCutTolerance;
0518 const float pT2perRadius = square(highland / pTPerHelixRadius);
0519 sigmapT2perRadius = pT2perRadius * square(2 * sigmaScattering);
0520 multipleScattering2 = maxScatteringAngle2 * square(sigmaScattering);
0521 }
0522
0523 std::unique_ptr<TripletSeedFinder> TripletSeedFinder::create(
0524 const DerivedConfig& config) {
0525 using BooleanOptions =
0526 boost::mp11::mp_list<std::bool_constant<false>, std::bool_constant<true>>;
0527
0528 using UseStripInfoOptions = BooleanOptions;
0529 using SortedByCotThetaOptions = BooleanOptions;
0530
0531 using TripletOptions =
0532 boost::mp11::mp_product<boost::mp11::mp_list, UseStripInfoOptions,
0533 SortedByCotThetaOptions>;
0534
0535 std::unique_ptr<TripletSeedFinder> result;
0536 boost::mp11::mp_for_each<TripletOptions>([&](auto option) {
0537 using OptionType = decltype(option);
0538
0539 using UseStripInfo = boost::mp11::mp_at_c<OptionType, 0>;
0540 using SortedByCotTheta = boost::mp11::mp_at_c<OptionType, 1>;
0541
0542 if (config.useStripInfo != UseStripInfo::value ||
0543 config.sortedByCotTheta != SortedByCotTheta::value) {
0544 return;
0545 }
0546
0547
0548 if (result != nullptr) {
0549 throw std::runtime_error(
0550 "TripletSeedFinder: Multiple implementations found for one "
0551 "configuration");
0552 }
0553
0554
0555 result =
0556 std::make_unique<Impl<UseStripInfo::value, SortedByCotTheta::value>>(
0557 config);
0558 });
0559 if (result == nullptr) {
0560 throw std::runtime_error(
0561 "TripletSeedFinder: No implementation found for the given "
0562 "configuration");
0563 }
0564 return result;
0565 }
0566
0567 }