traccc/seeding/triplet_finding_helper.hpp

0001 /** TRACCC library, part of the ACTS project (R&D line)
0002  *
0003  * (c) 2021-2025 CERN for the benefit of the ACTS project
0004  *
0005  * Mozilla Public License Version 2.0
0006  */
0007
0008 #pragma once
0009
0010 // Local include(s).
0011 #include "traccc/definitions/math.hpp"
0012 #include "traccc/edm/spacepoint_collection.hpp"
0013 #include "traccc/seeding/detail/doublet.hpp"
0014 #include "traccc/seeding/detail/lin_circle.hpp"
0015 #include "traccc/seeding/detail/triplet.hpp"
0016
0017 namespace traccc {
0018
0019 // helper function used for both cpu and gpu
0020 struct triplet_finding_helper {
0021   /// Check if two doublets with common middle spacepoint can form a triplet
0022   ///
0023   /// @param spM is middle spacepoint
0024   /// @param lb is transformed coordinate of middle-bottom doublet
0025   /// @param lt is transformed coordinate of middle-top doublet
0026   /// @param config is configuration parameter
0027   /// @param iSinTheta2 is the square of sin of pitch angle
0028   /// @param scatteringInRegion2 is the threshold for scattering angle for the
0029   /// lower pT cut
0030   /// @param curvature is curvature of triplet
0031   /// @param impact_parameter is impact parameter of triplet
0032   ///
0033   /// @return boolean value for compatibility
0034   template <typename T>
0035   static inline TRACCC_HOST_DEVICE bool isCompatible(
0036       const edm::spacepoint<T>& spM, const lin_circle& lb, const lin_circle& lt,
0037       const seedfinder_config& config, const scalar& iSinTheta2,
0038       const scalar& scatteringInRegion2, scalar& curvature,
0039       scalar& impact_parameter);
0040 };
0041
0042 template <typename T>
0043 bool TRACCC_HOST_DEVICE triplet_finding_helper::isCompatible(
0044     const edm::spacepoint<T>& spM, const lin_circle& lb, const lin_circle& lt,
0045     const seedfinder_config& config, const scalar& iSinTheta2,
0046     const scalar& scatteringInRegion2, scalar& curvature,
0047     scalar& impact_parameter) {
0048   // add errors of spB-spM and spM-spT pairs and add the correlation term
0049   // for errors on spM
0050   scalar error2 = lt.Er() + lb.Er() +
0051                   static_cast<scalar>(2.f) *
0052                       (lb.cotTheta() * lt.cotTheta() * spM.radius_variance() +
0053                        spM.z_variance()) *
0054                       lb.iDeltaR() * lt.iDeltaR();
0055
0056   scalar deltaCotTheta = lb.cotTheta() - lt.cotTheta();
0057   scalar deltaCotTheta2 = deltaCotTheta * deltaCotTheta;
0058   scalar error{0.f};
0059   scalar dCotThetaMinusError2{0.f};
0060
0061   // if the error is larger than the difference in theta, no need to
0062   // compare with scattering
0063   if (deltaCotTheta2 - error2 > 0) {
0064     deltaCotTheta = math::fabs(deltaCotTheta);
0065     // if deltaTheta larger than the scattering for the lower pT cut, skip
0066     error = std::sqrt(error2);
0067     dCotThetaMinusError2 = deltaCotTheta2 + error2 -
0068                            static_cast<scalar>(2.) * deltaCotTheta * error;
0069     // avoid taking root of scatteringInRegion
0070     // if left side of ">" is positive, both sides of unequality can be
0071     // squared
0072     // (scattering is always positive)
0073     if (dCotThetaMinusError2 > scatteringInRegion2) {
0074       return false;
0075     }
0076   }
0077
0078   // protects against division by 0
0079   scalar dU = lt.U() - lb.U();
0080   if (dU == static_cast<scalar>(0.)) {
0081     return false;
0082   }
0083
0084   // A and B are evaluated as a function of the circumference parameters
0085   // x_0 and y_0
0086   scalar A = (lt.V() - lb.V()) / dU;
0087   scalar S2 = static_cast<scalar>(1.) + A * A;
0088   scalar B = lb.V() - A * lb.U();
0089   scalar B2 = B * B;
0090   // sqrt(S2)/B = 2 * helixradius
0091   // calculated radius must not be smaller than minimum radius
0092   if (S2 < B2 * config.minHelixDiameter2) {
0093     return false;
0094   }
0095
0096   // 1/helixradius: (B/sqrt(S2))*2 (we leave everything squared)
0097   scalar iHelixDiameter2 = B2 / S2;
0098   // calculate scattering for p(T) calculated from seed curvature
0099   scalar pT2scatter =
0100       static_cast<scalar>(4.) * iHelixDiameter2 * config.pT2perRadius;
0101   // if pT > maxPtScattering, calculate allowed scattering angle using
0102   // maxPtScattering instead of pt.
0103   scalar pT =
0104       config.pTPerHelixRadius * std::sqrt(S2 / B2) / static_cast<scalar>(2.);
0105   if (pT > config.maxPtScattering) {
0106     scalar pTscatter = config.highland / config.maxPtScattering;
0107     pT2scatter = pTscatter * pTscatter;
0108   }
0109   // convert p(T) to p scaling by sin^2(theta) AND scale by 1/sin^4(theta)
0110   // from rad to deltaCotTheta
0111   scalar p2scatter = pT2scatter * iSinTheta2;
0112   // if deltaTheta larger than allowed scattering for calculated pT, skip
0113   if ((deltaCotTheta2 - error2 > static_cast<scalar>(0.)) &&
0114       (dCotThetaMinusError2 >
0115        p2scatter * config.sigmaScattering * config.sigmaScattering)) {
0116     return false;
0117   }
0118
0119   // calculate curvature
0120   curvature = B / std::sqrt(S2);
0121
0122   // A and B allow calculation of impact params in U/V plane with linear
0123   // function
0124   // (in contrast to having to solve a quadratic function in x/y plane)
0125   impact_parameter = math::fabs((A - B * spM.radius()) * spM.radius());
0126
0127   if (impact_parameter > config.impactMax) {
0128     return false;
0129   }
0130
0131   return true;
0132 }
0133
0134 }  // namespace traccc