traccc/seeding/doublet_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/seeding_config.hpp"
0016 #include "traccc/seeding/detail/spacepoint_type.hpp"
0017
0018 namespace traccc {
0019
0020 // helper functions used for both cpu and gpu
0021 struct doublet_finding_helper {
0022   /// Check if two spacepoints form doublets
0023   ///
0024   /// @param sp1 is middle spacepoint
0025   /// @param sp2 is bottom or top spacepoint
0026   /// @param config is configuration parameter
0027   /// @tparam otherSpType is whether it is for middle-bottom or middle-top
0028   /// doublet
0029   ///
0030   /// @return boolean value for compatibility
0031   ///
0032   template <details::spacepoint_type otherSpType, typename T1, typename T2>
0033   static inline TRACCC_HOST_DEVICE bool isCompatible(
0034       const edm::spacepoint<T1>& sp1, const edm::spacepoint<T2>& sp2,
0035       const seedfinder_config& config);
0036
0037   /// Do the conformal transformation on doublet's coordinate
0038   ///
0039   /// @param sp1 is middle spacepoint
0040   /// @param sp2 is bottom or top spacepoint
0041   /// @tparam otherSpType is whether it is for middle-bottom or middle-top
0042   /// doublet
0043   ///
0044   /// @return lin_circle which contains the transformed coordinate information
0045   ///
0046   template <details::spacepoint_type otherSpType, typename T1, typename T2>
0047   static inline TRACCC_HOST_DEVICE lin_circle transform_coordinates(
0048       const edm::spacepoint<T1>& sp1, const edm::spacepoint<T2>& sp2);
0049 };
0050
0051 template <details::spacepoint_type otherSpType, typename T1, typename T2>
0052 bool TRACCC_HOST_DEVICE doublet_finding_helper::isCompatible(
0053     const edm::spacepoint<T1>& sp1, const edm::spacepoint<T2>& sp2,
0054     const seedfinder_config& config) {
0055   static_assert(otherSpType == details::spacepoint_type::bottom ||
0056                 otherSpType == details::spacepoint_type::top);
0057
0058   scalar deltaR, cotTheta, zOrigin;
0059   if constexpr (otherSpType == details::spacepoint_type::bottom) {
0060     // check if R distance is too small, because bins are not R-sorted
0061     deltaR = sp1.radius() - sp2.radius();
0062     // actually cotTheta * deltaR to avoid division by 0 statements
0063     cotTheta = sp1.z() - sp2.z();
0064     // actually zOrigin * deltaR to avoid division by 0 statements
0065     zOrigin = sp1.z() * deltaR - sp1.radius() * cotTheta;
0066   } else {
0067     // check if R distance is too small, because bins are not R-sorted
0068     deltaR = sp2.radius() - sp1.radius();
0069     // actually cotTheta * deltaR to avoid division by 0 statements
0070     cotTheta = (sp2.z() - sp1.z());
0071     // actually zOrigin * deltaR to avoid division by 0 statements
0072     zOrigin = sp1.z() * deltaR - sp1.radius() * cotTheta;
0073   }
0074
0075   if ((deltaR >= config.deltaRMax) || (deltaR <= config.deltaRMin) ||
0076       (math::fabs(cotTheta) >= config.cotThetaMax * deltaR) ||
0077       (zOrigin <= config.collisionRegionMin * deltaR) ||
0078       (zOrigin >= config.collisionRegionMax * deltaR) ||
0079       math::fabs(cotTheta) >= config.deltaZMax) {
0080     return false;
0081   }
0082
0083   /*
0084    * The following cut is capable of discriminating some doublets on the
0085    * basis that it is impossible to find a third spacepoint for the doublet
0086    * that will keep the resulting triplet inside the helix radius bound.
0087    * This explanation is enhanced with Geogebra commands with can be entered
0088    * into the application directly to provide a visual "proof" of why this
0089    * cut works.
0090    *
0091    * We will start by creating two spacepoints at arbitrary locations (they
0092    * can be moved as desired):
0093    *
0094    * ```
0095    * A = (2, 4)
0096    * B = (3, 12)
0097    * ```
0098    *
0099    * We will also define a radius $R$:
0100    *
0101    * ```
0102    * R = 10
0103    * ```
0104    *
0105    * Next, we consider the fact that two points and a radius define exactly
0106    * two circles through those points and with that radius. This makes sense
0107    * because three points (six degrees of freedom) precisely define a single
0108    * circle, and two points and a radius (five DoFs) define two circles. We
0109    * will construct those circles now, with the radius being the minimum
0110    * helix radius from the configuration.
0111    *
0112    * To find the midpoints of these circles, we will first find the
0113    * perpendicular bisector of points $A$ and $B$:
0114    *
0115    * ```
0116    * M = 0.5 * (A + B)
0117    * ```
0118    */
0119   scalar midX = 0.5f * (sp1.x() + sp2.x());
0120   scalar midY = 0.5f * (sp1.y() + sp2.y());
0121
0122   /*
0123    * Then we will compute the slope of the perpendicular bisector:
0124    *
0125    * ```
0126    * s = (y(B) - y(A)) / (x(B) - x(A))
0127    * ```
0128    */
0129   scalar slope = (sp2.y() - sp1.y()) / (sp2.x() - sp1.x());
0130
0131   /*
0132    * Next, we can simply place circle midpoints on our perpendicular
0133    * bisector, but we cannot simply use the radius $R$ as the distance from
0134    * the midpoint of $A$ and $B$, we have account for the length of the
0135    * sagitta:
0136    *
0137    * ```
0138    * dX = x(B) - x(A)
0139    * dY = y(B) - y(A)
0140    * dXY2 = dX * dX + dY * dY
0141    * q = sqrt(R * R - dXY2 / 4)
0142    * ```
0143    */
0144   scalar deltaX = sp2.x() - sp1.x();
0145   scalar deltaY = sp2.y() - sp1.y();
0146   scalar deltaXY2 = deltaX * deltaX + deltaY * deltaY;
0147   scalar sagittaLength = math::sqrt(
0148       config.minHelixRadius * config.minHelixRadius - deltaXY2 / 4.f);
0149
0150   /*
0151    * We then compute the central angle between the points $A$, $B$, and the
0152    * midpoints of the circles we want to construct. Naively, this can be
0153    * done using the trigonomic functions:
0154    *
0155    * ```
0156    * theta = atan(1 / slope)
0157    * ```
0158    *
0159    * After which we can compute the delta between the midpoint of $A$ and
0160    * $B$ and the midpoints of our circles:
0161    *
0162    * ```
0163    * mdX = (R - q) * cos(theta)
0164    * mdY = (R - q) * sin(theta)
0165    * ```
0166    *
0167    * However, some trigonomy allows us to reduce this:
0168    *
0169    * ```
0170    * denom = sqrt((s * s + 1) / (s * s))
0171    * cosTheta = 1 / denom
0172    * sinTheta = -1 / (s * denom)
0173    * mdX = q * cosTheta
0174    * mdY = q * sinTheta
0175    * ```
0176    */
0177   scalar denom = math::sqrt((slope * slope + 1) / (slope * slope));
0178   scalar cosCentralAngle = 1.f / denom;
0179   scalar sinCentralAngle = -1.f / (slope * denom);
0180
0181   scalar mpDeltaX = sagittaLength * cosCentralAngle;
0182   scalar mpDeltaY = sagittaLength * sinCentralAngle;
0183
0184   /*
0185    * We now easily find the midpoints of the required circles:
0186    *
0187    * ```
0188    * V = (x(M) + mdX, y(M) + mdY)
0189    * W = (x(M) - mdX, y(M) - mdY)
0190    * ```
0191    */
0192   scalar mp1X = midX + mpDeltaX;
0193   scalar mp2X = midX - mpDeltaX;
0194   scalar mp1Y = midY + mpDeltaY;
0195   scalar mp2Y = midY - mpDeltaY;
0196
0197   /*
0198    * Finally, we compute the radii of the circles. The crucial intuition
0199    * here is that if either of the newly constructed circles _completely_
0200    * contain a circle of radius $R$ around the origin, then we can never
0201    * find a third spacepoint to complete the triplet. Thus, we compute the
0202    * radii. We leave them squared in the C++ code to avoid an unnecessary
0203    * square root.
0204    */
0205   scalar mp1R2 = mp1X * mp1X + mp1Y * mp1Y;
0206   scalar mp2R2 = mp2X * mp2X + mp2Y * mp2Y;
0207
0208   if (math::min(mp1R2, mp2R2) <= ((config.minHelixRadius - config.impactMax) *
0209                                   (config.minHelixRadius - config.impactMax))) {
0210     return false;
0211   }
0212
0213   return true;
0214 }
0215
0216 template <details::spacepoint_type otherSpType, typename T1, typename T2>
0217 lin_circle TRACCC_HOST_DEVICE doublet_finding_helper::transform_coordinates(
0218     const edm::spacepoint<T1>& sp1, const edm::spacepoint<T2>& sp2) {
0219   static_assert(otherSpType == details::spacepoint_type::bottom ||
0220                 otherSpType == details::spacepoint_type::top);
0221
0222   const scalar& xM = sp1.x();
0223   const scalar& yM = sp1.y();
0224   const scalar& zM = sp1.z();
0225   const scalar& rM = sp1.radius();
0226   const scalar& varianceZM = sp1.z_variance();
0227   const scalar& varianceRM = sp1.radius_variance();
0228   scalar cosPhiM = xM / rM;
0229   scalar sinPhiM = yM / rM;
0230
0231   scalar deltaX = sp2.x() - xM;
0232   scalar deltaY = sp2.y() - yM;
0233   scalar deltaZ = sp2.z() - zM;
0234   // calculate projection fraction of spM->sp vector pointing in same
0235   // direction as
0236   // vector origin->spM (x) and projection fraction of spM->sp vector pointing
0237   // orthogonal to origin->spM (y)
0238   scalar x = deltaX * cosPhiM + deltaY * sinPhiM;
0239   scalar y = deltaY * cosPhiM - deltaX * sinPhiM;
0240   // 1/(length of M -> SP)
0241   scalar iDeltaR2 =
0242       static_cast<scalar>(1.) / (deltaX * deltaX + deltaY * deltaY);
0243   scalar iDeltaR = std::sqrt(iDeltaR2);
0244   // cot_theta = (deltaZ/deltaR)
0245   scalar cot_theta = deltaZ * iDeltaR;
0246   if constexpr (otherSpType == details::spacepoint_type::bottom) {
0247     cot_theta = -cot_theta;
0248   }
0249   // VERY frequent (SP^3) access
0250   lin_circle l;
0251   l.m_cotTheta = cot_theta;
0252   // location on z-axis of this SP-duplet
0253   l.m_Zo = zM - rM * cot_theta;
0254   l.m_iDeltaR = iDeltaR;
0255   // transformation of circle equation (x,y) into linear equation (u,v)
0256   // x^2 + y^2 - 2x_0*x - 2y_0*y = 0
0257   // is transformed into
0258   // 1 - 2x_0*u - 2y_0*v = 0
0259   // using the following m_U and m_V
0260   // (u = A + B*v); A and B are created later on
0261   l.m_U = x * iDeltaR2;
0262   l.m_V = y * iDeltaR2;
0263   // error term for sp-pair without correlation of middle space point
0264   l.m_Er = ((varianceZM + sp2.z_variance()) +
0265             (cot_theta * cot_theta) * (varianceRM + sp2.radius_variance())) *
0266            iDeltaR2;
0267
0268   return l;
0269 }
0270
0271 }  // namespace traccc