Back to home page

EIC code displayed by LXR

 
 

    


File indexing completed on 2025-08-28 08:12:01

0001 // This file is part of the ACTS project.
0002 //
0003 // Copyright (C) 2016 CERN for the benefit of the ACTS project
0004 //
0005 // This Source Code Form is subject to the terms of the Mozilla Public
0006 // License, v. 2.0. If a copy of the MPL was not distributed with this
0007 // file, You can obtain one at https://mozilla.org/MPL/2.0/.
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 /// Check the compatibility of strip space point coordinates in xyz assuming
0026 /// the Bottom-Middle direction with the strip measurement details
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   // cross product between top strip vector and spacePointPosition
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   // scalar product between bottom strip vector and d1
0044   const float bd1 = bottomStripVector[0] * d1[0] +
0045                     bottomStripVector[1] * d1[1] + bottomStripVector[2] * d1[2];
0046 
0047   // compatibility check using distance between strips to evaluate if
0048   // spacepointPosition is inside the bottom detector element
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   // cross product between bottom strip vector and spacePointPosition
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   // compatibility check using distance between strips to evaluate if
0066   // spacePointPosition is inside the top detector element
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   // if arrive here spacePointPosition is compatible with strip directions and
0074   // detector elements
0075 
0076   const std::array<float, 3>& topStripCenter = sp.topStripCenter();
0077 
0078   // spacePointPosition corrected with respect to the top strip position and
0079   // direction and the distance between the strips
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     // Reserve enough space, in case current capacity is too little
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     // 1+(cot^2(theta)) = 1/sin^2(theta)
0114     const float iSinTheta2 = 1 + cotThetaB * cotThetaB;
0115     const float sigmaSquaredPtDependent = iSinTheta2 * m_cfg.sigmapT2perRadius;
0116     // calculate max scattering for min momentum at the seed's theta angle
0117     // scaling scatteringAngle^2 by sin^2(theta) to convert pT^2 to p^2
0118     // accurate would be taking 1/atan(thetaBottom)-1/atan(thetaTop) <
0119     // scattering
0120     // but to avoid trig functions we approximate cot by scaling by
0121     // 1/sin^4(theta)
0122     // resolving with pT to p scaling --> only divide by sin^2(theta)
0123     // max approximation error for allowed scattering angles of 0.04 rad at
0124     // eta=infinity: ~8.5%
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       // use geometric average
0135       const float cotThetaAvg2 = cotThetaB * cotThetaT;
0136 
0137       // add errors of spB-spM and spM-spT pairs and add the correlation term
0138       // for errors on spM
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       // Apply a cut on the compatibility between the r-z slope of the two
0147       // seed segments. This is done by comparing the squared difference
0148       // between slopes, and comparing to the squared uncertainty in this
0149       // difference - we keep a seed if the difference is compatible within
0150       // the assumed uncertainties. The uncertainties get contribution from
0151       // the  space-point-related squared error (error2) and a scattering term
0152       // calculated assuming the minimum pt we expect to reconstruct
0153       // (scatteringInRegion2). This assumes gaussian error propagation which
0154       // allows just adding the two errors if they are uncorrelated (which is
0155       // fair for scattering and measurement uncertainties)
0156       if (deltaCotTheta2 > error2 + scatteringInRegion2) {
0157         if constexpr (sortedByCotTheta) {
0158           // skip top SPs based on cotTheta sorting when producing triplets
0159           // break if cotTheta from bottom SP < cotTheta from top SP because
0160           // the SP are sorted by cotTheta
0161           if (cotThetaB < cotThetaT) {
0162             break;
0163           }
0164           topDoubletOffset = topDoubletIndex + 1;
0165         }
0166         continue;
0167       }
0168 
0169       const float dU = topDoublet.u() - Ub;
0170       // protects against division by 0
0171       if (dU == 0) {
0172         continue;
0173       }
0174       // A and B are evaluated as a function of the circumference parameters
0175       // x_0 and y_0
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       // sqrt(S2)/B = 2 * helixradius
0182       // calculated radius must not be smaller than minimum radius
0183       if (S2 < B2 * m_cfg.minHelixDiameter2) {
0184         continue;
0185       }
0186 
0187       // refinement of the cut on the compatibility between the r-z slope of
0188       // the two seed segments using a scattering term scaled by the actual
0189       // measured pT (p2scatterSigma)
0190       const float iHelixDiameter2 = B2 / S2;
0191       // convert p(T) to p scaling by sin^2(theta) AND scale by 1/sin^4(theta)
0192       // from rad to deltaCotTheta
0193       float p2scatterSigma = iHelixDiameter2 * sigmaSquaredPtDependent;
0194       if (!std::isinf(m_cfg.maxPtScattering)) {
0195         // if pT > maxPtScattering, calculate allowed scattering angle using
0196         // maxPtScattering instead of pt.
0197         // To avoid 0-divison the pT check is skipped in case of B2==0, and
0198         // p2scatterSigma is calculated directly from maxPtScattering
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       // if deltaTheta larger than allowed scattering for calculated pT, skip
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       // A and B allow calculation of impact params in U/V plane with linear
0219       // function
0220       // (in contrast to having to solve a quadratic function in x/y plane)
0221       const float im = std::abs((A - B * rM) * rM);
0222       if (im > m_cfg.impactMax) {
0223         continue;
0224       }
0225 
0226       // inverse diameter is signed depending on if the curvature is
0227       // positive/negative in phi
0228       tripletTopCandidates.emplace_back(spT, B / std::sqrt(S2), im);
0229     }  // loop on tops
0230 
0231     if constexpr (sortedByCotTheta) {
0232       // remove the top doublets that were skipped due to cotTheta sorting
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     // Reserve enough space, in case current capacity is too little
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     // 1+(cot^2(theta)) = 1/sin^2(theta)
0260     const float iSinTheta2 = 1 + cotThetaB * cotThetaB;
0261     const float sigmaSquaredPtDependent = iSinTheta2 * m_cfg.sigmapT2perRadius;
0262     // calculate max scattering for min momentum at the seed's theta angle
0263     // scaling scatteringAngle^2 by sin^2(theta) to convert pT^2 to p^2
0264     // accurate would be taking 1/atan(thetaBottom)-1/atan(thetaTop) <
0265     // scattering
0266     // but to avoid trig functions we approximate cot by scaling by
0267     // 1/sin^4(theta)
0268     // resolving with pT to p scaling --> only divide by sin^2(theta)
0269     // max approximation error for allowed scattering angles of 0.04 rad at
0270     // eta=infinity: ~8.5%
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     // coordinate transformation and checks for middle spacepoint
0277     // x and y terms for the rotation from UV to XY plane
0278     const std::array<float, 2> rotationTermsUVtoXY = {cosPhiM * sinTheta,
0279                                                       sinPhiM * sinTheta};
0280 
0281     for (auto topDoublet : topDoublets) {
0282       // protects against division by 0
0283       float dU = topDoublet.u() - Ub;
0284       if (dU == 0) {
0285         continue;
0286       }
0287       // A and B are evaluated as a function of the circumference parameters
0288       // x_0 and y_0
0289       const float A0 = (topDoublet.v() - Vb) / dU;
0290 
0291       const float zPositionMiddle = cosTheta * std::sqrt(1 + A0 * A0);
0292 
0293       // position of Middle SP converted from UV to XY assuming cotTheta
0294       // evaluated from the Bottom and Middle SPs double
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       // coordinate transformation and checks for bottom spacepoint
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       // coordinate transformation and checks for top spacepoint
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       // bottom and top coordinates in the spM reference frame
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       // use arithmetic average
0354       const float averageCotTheta = 0.5f * (cotThetaB + cotThetaT);
0355       const float cotThetaAvg2 = averageCotTheta * averageCotTheta;
0356 
0357       // add errors of spB-spM and spM-spT pairs and add the correlation term
0358       // for errors on spM
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       // Apply a cut on the compatibility between the r-z slope of the two
0367       // seed segments. This is done by comparing the squared difference
0368       // between slopes, and comparing to the squared uncertainty in this
0369       // difference - we keep a seed if the difference is compatible within
0370       // the assumed uncertainties. The uncertainties get contribution from
0371       // the  space-point-related squared error (error2) and a scattering term
0372       // calculated assuming the minimum pt we expect to reconstruct
0373       // (scatteringInRegion2). This assumes gaussian error propagation which
0374       // allows just adding the two errors if they are uncorrelated (which is
0375       // fair for scattering and measurement uncertainties)
0376       if (deltaCotTheta2 > error2 + scatteringInRegion2) {
0377         // skip top SPs based on cotTheta sorting when producing triplets
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       // protects against division by 0
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       // sqrt(S2)/B = 2 * helixradius
0403       // calculated radius must not be smaller than minimum radius
0404       if (S2 < B2 * m_cfg.minHelixDiameter2) {
0405         continue;
0406       }
0407 
0408       // refinement of the cut on the compatibility between the r-z slope of
0409       // the two seed segments using a scattering term scaled by the actual
0410       // measured pT (p2scatterSigma)
0411       const float iHelixDiameter2 = B2 / S2;
0412       // convert p(T) to p scaling by sin^2(theta) AND scale by 1/sin^4(theta)
0413       // from rad to deltaCotTheta
0414       float p2scatterSigma = iHelixDiameter2 * sigmaSquaredPtDependent;
0415       if (!std::isinf(m_cfg.maxPtScattering)) {
0416         // if pT > maxPtScattering, calculate allowed scattering angle using
0417         // maxPtScattering instead of pt.
0418         // To avoid 0-divison the pT check is skipped in case of B2==0, and
0419         // p2scatterSigma is calculated directly from maxPtScattering
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       // if deltaTheta larger than allowed scattering for calculated pT, skip
0429       if (deltaCotTheta2 > error2 + p2scatterSigma) {
0430         continue;
0431       }
0432 
0433       // A and B allow calculation of impact params in U/V plane with linear
0434       // function
0435       // (in contrast to having to solve a quadratic function in x/y plane)
0436       const float im = std::abs((A - B * rMxy) * rMxy);
0437       if (im > m_cfg.impactMax) {
0438         continue;
0439       }
0440 
0441       // inverse diameter is signed depending on if the curvature is
0442       // positive/negative in phi
0443       tripletTopCandidates.emplace_back(spT.index(), B / std::sqrt(S2), im);
0444     }  // loop on tops
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 }  // namespace
0494 
0495 TripletSeedFinder::DerivedConfig::DerivedConfig(const Config& config,
0496                                                 float bFieldInZ_)
0497     : Config(config), bFieldInZ(bFieldInZ_) {
0498   using namespace Acts::UnitLiterals;
0499 
0500   // similar to `theta0Highland` in `Core/src/Material/Interactions.cpp`
0501   {
0502     const double xOverX0 = radLengthPerSeed;
0503     const double q2OverBeta2 = 1;  // q^2=1, beta^2~1
0504     // RPP2018 eq. 33.15 (treats beta and q² consistently)
0505     const double t = std::sqrt(xOverX0 * q2OverBeta2);
0506     // log((x/X0) * (q²/beta²)) = log((sqrt(x/X0) * (q/beta))²)
0507     //                          = 2 * log(sqrt(x/X0) * (q/beta))
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   // bFieldInZ is in (pT/radius) natively, no need for conversion
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;  // skip if the configuration does not match
0545     }
0546 
0547     // check if we already have an implementation for this configuration
0548     if (result != nullptr) {
0549       throw std::runtime_error(
0550           "TripletSeedFinder: Multiple implementations found for one "
0551           "configuration");
0552     }
0553 
0554     // create the implementation for the given configuration
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 }  // namespace Acts::Experimental