refactor: Change transformCoordinates to avoid recalculation of variables in Seedfinder (acts-project#1831)

This PR tries to optimise SeedFinder by removing unnecessary calculations

Author: LuisFelipeCoelho

Core/include/Acts/Seeding/SeedFilterConfig.hpp

@@ -56,7 +56,7 @@ struct SeedFilterConfig {
 
   // maximum number of lower quality seeds in seed confirmation
   std::size_t maxSeedsPerSpMConf = std::numeric_limits<std::size_t>::max();
-  // maximum number of quality seeds for each middle-bottom SP-duplet in seed
+  // maximum number of quality seeds for each middle-bottom SP-dublet in seed
   // confirmation if the limit is reached we check if there is a lower quality
   // seed to be replaced
   std::size_t maxQualitySeedsPerSpMConf =

Core/include/Acts/Seeding/SeedFinder.hpp

@@ -128,16 +128,36 @@ class SeedFinder {
       const sp_range_t& topSPs) const;
 
  private:
+  /// Iterates over dublets and tests the compatibility between them by applying
+  /// a series of cuts that can be tested with only two SPs
+  /// @param spacePointData object contaning the spacepoint data
+  /// @param options frequently changing configuration (like beam position)
+  /// @param grid spacepoint grid
+  /// @param otherSPsNeighbours inner or outer space points to be used in the dublet
+  /// @param mediumSP space point candidate to be used as middle SP in a seed
+  /// @param linCircleVec vector contining inner or outer SP parameters after reference frame transformation to the u-v space
+  /// @param outVec Output object containing top or bottom SPs that are compatible with a certain middle SPs
+  /// @param deltaRMinSP minimum allowed r-distance between dublet components
+  /// @param deltaRMaxSP maximum allowed r-distance between dublet components
+  /// @param isBottom wheter otherSPs contains outer or inner SPs
   template <typename out_range_t>
   void getCompatibleDoublets(
+      Acts::SpacePointData& spacePointData,
       const Acts::SeedFinderOptions& options,
       const Acts::SpacePointGrid<external_spacepoint_t>& grid,
-      boost::container::small_vector<Acts::Neighbour<external_spacepoint_t>, 9>&
-          otherSPs,
+      boost::container::small_vector<Neighbour<external_spacepoint_t>, 9>&
+          otherSPsNeighbours,
       const InternalSpacePoint<external_spacepoint_t>& mediumSP,
-      out_range_t& outVec, const float& deltaRMinSP, const float& deltaRMaxSP,
-      bool isBottom) const;
+      std::vector<LinCircle>& linCircleVec, out_range_t& outVec,
+      const float& deltaRMinSP, const float& deltaRMaxSP, bool isBottom) const;
 
+  /// Iterates over the seed candidates tests the compatibility between three
+  /// SPs and calls for the seed confirmation
+  /// @param spacePointData object contaning the spacepoint data
+  /// @param SpM space point candidate to be used as middle SP in a seed
+  /// @param options frequently changing configuration (like beam position)
+  /// @param seedFilterState State object that holds memory used in SeedFilter
+  /// @param state State object that holds memory used
   void filterCandidates(Acts::SpacePointData& spacePointData,
                         const InternalSpacePoint<external_spacepoint_t>& SpM,
                         const Acts::SeedFinderOptions& options,

Core/include/Acts/Seeding/SeedFinder.ipp

@@ -120,7 +120,9 @@ void SeedFinder<external_spacepoint_t, platform_t>::createSeedsForGroup(
       }
     }
 
-    getCompatibleDoublets(options, grid, state.topNeighbours, *spM.get(),
+    // Iterate over middle-top dublets
+    getCompatibleDoublets(state.spacePointData, options, grid,
+                          state.topNeighbours, *spM.get(), state.linCircleTop,
                           state.compatTopSP, m_config.deltaRMinTopSP,
                           m_config.deltaRMaxTopSP, false);
 
@@ -151,9 +153,11 @@ void SeedFinder<external_spacepoint_t, platform_t>::createSeedsForGroup(
       }
     }
 
-    getCompatibleDoublets(options, grid, state.bottomNeighbours, *spM.get(),
-                          state.compatBottomSP, m_config.deltaRMinBottomSP,
-                          m_config.deltaRMaxBottomSP, true);
+    // Iterate over middle-bottom dublets
+    getCompatibleDoublets(
+        state.spacePointData, options, grid, state.bottomNeighbours, *spM.get(),
+        state.linCircleBottom, state.compatBottomSP, m_config.deltaRMinBottomSP,
+        m_config.deltaRMaxBottomSP, true);
 
     // no bottom SP found -> try next spM
     if (state.compatBottomSP.empty()) {
@@ -175,23 +179,31 @@ template <typename external_spacepoint_t, typename platform_t>
 template <typename out_range_t>
 inline void
 SeedFinder<external_spacepoint_t, platform_t>::getCompatibleDoublets(
+    Acts::SpacePointData& spacePointData,
     const Acts::SeedFinderOptions& options,
     const Acts::SpacePointGrid<external_spacepoint_t>& grid,
     boost::container::small_vector<Neighbour<external_spacepoint_t>, 9>&
         otherSPsNeighbours,
     const InternalSpacePoint<external_spacepoint_t>& mediumSP,
-    out_range_t& outVec, const float& deltaRMinSP, const float& deltaRMaxSP,
-    bool isBottom) const {
+    std::vector<LinCircle>& linCircleVec, out_range_t& outVec,
+    const float& deltaRMinSP, const float& deltaRMaxSP, bool isBottom) const {
   const int sign = isBottom ? -1 : 1;
 
   outVec.clear();
+  linCircleVec.clear();
 
   const float& rM = mediumSP.radius();
   const float& xM = mediumSP.x();
   const float& yM = mediumSP.y();
   const float& zM = mediumSP.z();
-  const float ratio_xM_rM = xM / rM;
-  const float ratio_yM_rM = yM / rM;
+  const float& varianceRM = mediumSP.varianceR();
+  const float& varianceZM = mediumSP.varianceZ();
+  const float cosPhiM = xM / rM;
+  const float sinPhiM = yM / rM;
+  float vIPAbs = 0;
+  if (m_config.interactionPointCut) {
+    vIPAbs = m_config.impactMax / (rM * rM);
+  }
 
   for (auto& otherSPCol : otherSPsNeighbours) {
     const auto& otherSPs = grid.at(otherSPCol.index);
@@ -209,6 +221,7 @@ SeedFinder<external_spacepoint_t, platform_t>::getCompatibleDoublets(
       const float rO = otherSP->radius();
       float deltaR = sign * (rO - rM);
 
+      // if r-distance is too small, try next SP in bin
       if (deltaR < deltaRMinSP) {
         if (isBottom) {
           break;
@@ -251,34 +264,69 @@ SeedFinder<external_spacepoint_t, platform_t>::getCompatibleDoublets(
         continue;
       }
 
+      const float deltaX = otherSP->x() - xM;
+      const float deltaY = otherSP->y() - yM;
+
+      // calculate projection fraction of spM->sp vector pointing in same
+      // direction as
+      // vector origin->spM (x) and projection fraction of spM->sp vector
+      // pointing orthogonal to origin->spM (y)
+      const float xNewFrame = deltaX * cosPhiM + deltaY * sinPhiM;
+      const float yNewFrame = deltaY * cosPhiM - deltaX * sinPhiM;
+
+      const float deltaR2 = (deltaX * deltaX + deltaY * deltaY);
+      const float iDeltaR2 = 1. / deltaR2;
+
+      // conformal transformation u=x/(x²+y²) v=y/(x²+y²) transform the
+      // circle into straight lines in the u/v plane the line equation can
+      // be described in terms of aCoef and bCoef, where v = aCoef * u +
+      // bCoef
+      const float uT = xNewFrame * iDeltaR2;
+      const float vT = yNewFrame * iDeltaR2;
+
+      // continue if interactionPointCut is disabled
       if (not m_config.interactionPointCut) {
+        const float iDeltaR = std::sqrt(iDeltaR2);
+        cotTheta = deltaZ * iDeltaR;
+
+        // error term for sp-pair without correlation of middle space point
+        const float Er =
+            ((varianceZM + otherSP->varianceZ()) +
+             (cotTheta * cotTheta) * (varianceRM + otherSP->varianceR())) *
+            iDeltaR2;
+
+        // fill output vectors
+        linCircleVec.push_back(fillLineCircle(
+            {deltaZ * iDeltaR, iDeltaR, Er, uT, vT, xNewFrame, yNewFrame}));
+        spacePointData.setDeltaR(otherSP->index(),
+                                 std::sqrt(deltaR2 + (deltaZ * deltaZ)));
         outVec.push_back(otherSP.get());
         continue;
       }
 
-      const float xVal =
-          (otherSP->x() - xM) * ratio_xM_rM + (otherSP->y() - yM) * ratio_yM_rM;
-      const float yVal =
-          (otherSP->y() - yM) * ratio_xM_rM - (otherSP->x() - xM) * ratio_yM_rM;
-
-      if (std::abs(rM * yVal) <= sign * m_config.impactMax * xVal) {
+      if (std::abs(rM * yNewFrame) <= sign * m_config.impactMax * xNewFrame) {
+        const float iDeltaR = std::sqrt(iDeltaR2);
+        cotTheta = deltaZ * iDeltaR;
+
+        // error term for sp-pair without correlation of middle space point
+        const float Er =
+            ((varianceZM + otherSP->varianceZ()) +
+             (cotTheta * cotTheta) * (varianceRM + otherSP->varianceR())) *
+            iDeltaR2;
+
+        // fill output vectors
+        linCircleVec.push_back(fillLineCircle(
+            {deltaZ * iDeltaR, iDeltaR, Er, uT, vT, xNewFrame, yNewFrame}));
+        spacePointData.setDeltaR(otherSP->index(),
+                                 std::sqrt(deltaR2 + (deltaZ * deltaZ)));
         outVec.push_back(otherSP.get());
         continue;
       }
 
-      // conformal transformation u=x/(x²+y²) v=y/(x²+y²) transform the
-      // circle into straight lines in the u/v plane the line equation can
-      // be described in terms of aCoef and bCoef, where v = aCoef * u +
-      // bCoef
-      const float uT = xVal / (xVal * xVal + yVal * yVal);
-      const float vT = yVal / (xVal * xVal + yVal * yVal);
       // in the rotated frame the interaction point is positioned at x = -rM
       // and y ~= impactParam
       const float uIP = -1. / rM;
-      float vIP = m_config.impactMax / (rM * rM);
-      if (sign * yVal > 0.) {
-        vIP = -vIP;
-      }
+      const float vIP = (sign * yNewFrame > 0.) ? -vIPAbs : vIPAbs;
       // we can obtain aCoef as the slope dv/du of the linear function,
       // estimated using du and dv between the two SP bCoef is obtained by
       // inserting aCoef into the linear equation
@@ -290,13 +338,27 @@ SeedFinder<external_spacepoint_t, platform_t>::getCompatibleDoublets(
       if ((bCoef * bCoef) * options.minHelixDiameter2 > (1 + aCoef * aCoef)) {
         continue;
       }
+
+      const float iDeltaR = std::sqrt(iDeltaR2);
+      cotTheta = deltaZ * iDeltaR;
+
+      // error term for sp-pair without correlation of middle space point
+      const float Er =
+          ((varianceZM + otherSP->varianceZ()) +
+           (cotTheta * cotTheta) * (varianceRM + otherSP->varianceR())) *
+          iDeltaR2;
+
+      // fill output vectors
+      linCircleVec.push_back(fillLineCircle(
+          {deltaZ * iDeltaR, iDeltaR, Er, uT, vT, xNewFrame, yNewFrame}));
+      spacePointData.setDeltaR(otherSP->index(),
+                               std::sqrt(deltaR2 + (deltaZ * deltaZ)));
       outVec.push_back(otherSP.get());
     }
   }
 }
 
 template <typename external_spacepoint_t, typename platform_t>
-
 inline void SeedFinder<external_spacepoint_t, platform_t>::filterCandidates(
     Acts::SpacePointData& spacePointData,
     const InternalSpacePoint<external_spacepoint_t>& spM,
@@ -306,16 +368,8 @@ inline void SeedFinder<external_spacepoint_t, platform_t>::filterCandidates(
   float varianceRM = spM.varianceR();
   float varianceZM = spM.varianceZ();
 
-  state.linCircleBottom.clear();
-  state.linCircleTop.clear();
-
   std::size_t numTopSP = state.compatTopSP.size();
 
-  transformCoordinates(state.spacePointData, state.compatBottomSP, spM, true,
-                       state.linCircleBottom);
-  transformCoordinates(state.spacePointData, state.compatTopSP, spM, false,
-                       state.linCircleTop);
-
   // sort: make index vector
   std::vector<std::size_t> sorted_bottoms(state.linCircleBottom.size());
   for (std::size_t i(0); i < sorted_bottoms.size(); ++i) {
@@ -356,7 +410,6 @@ inline void SeedFinder<external_spacepoint_t, platform_t>::filterCandidates(
     }
 
     auto lb = state.linCircleBottom[b];
-    seedFilterState.zOrigin = lb.Zo;
     float cotThetaB = lb.cotTheta;
     float Vb = lb.V;
     float Ub = lb.U;
@@ -515,6 +568,7 @@ inline void SeedFinder<external_spacepoint_t, platform_t>::filterCandidates(
 
       float deltaCotTheta = cotThetaB - cotThetaT;
       float deltaCotTheta2 = deltaCotTheta * deltaCotTheta;
+
       // Apply a cut on the compatibility between the r-z slope of the two
       // seed segments. This is done by comparing the squared difference
       // between slopes, and comparing to the squared uncertainty in this
@@ -571,7 +625,6 @@ inline void SeedFinder<external_spacepoint_t, platform_t>::filterCandidates(
 
       // sqrt(S2)/B = 2 * helixradius
       // calculated radius must not be smaller than minimum radius
-
       if (S2 < B2 * options.minHelixDiameter2) {
         continue;
       }
@@ -629,6 +682,8 @@ inline void SeedFinder<external_spacepoint_t, platform_t>::filterCandidates(
       continue;
     }
 
+    seedFilterState.zOrigin = spM.z() - rM * lb.cotTheta;
+
     m_config.seedFilter->filterSeeds_2SpFixed(
         state.spacePointData, *state.compatBottomSP[b], spM, state.topSpVec,
         state.curvatures, state.impactParameters, seedFilterState,

Core/include/Acts/Seeding/SeedFinderOrthogonal.ipp

@@ -169,7 +169,7 @@ bool SeedFinderOrthogonal<external_spacepoint_t>::validTuple(
   }
 
   /*
-   * Cut: Ensure that the origin of the duplet (the intersection of the line
+   * Cut: Ensure that the origin of the dublet (the intersection of the line
    * between them with the z axis) lies within the collision region.
    */
   float zOrigin = zL - rL * cotTheta;
@@ -335,7 +335,6 @@ void SeedFinderOrthogonal<external_spacepoint_t>::filterCandidates(
     }
 
     auto lb = linCircleBottom[b];
-    seedFilterState.zOrigin = lb.Zo;
     float cotThetaB = lb.cotTheta;
     float Vb = lb.V;
     float Ub = lb.U;
@@ -476,6 +475,8 @@ void SeedFinderOrthogonal<external_spacepoint_t>::filterCandidates(
       continue;
     }
 
+    seedFilterState.zOrigin = middle.z() - rM * lb.cotTheta;
+
     m_config.seedFilter->filterSeeds_2SpFixed(
         spacePointData, *bottom[b], middle, top_valid, curvatures,
         impactParameters, seedFilterState, candidates_collector);
@@ -498,7 +499,7 @@ void SeedFinderOrthogonal<external_spacepoint_t>::processFromMiddleSP(
    * bottom_lh_v denotes the candidates bottom seed points, assuming that the
    * track has monotonically _increasing_ z position. bottom_hl_v denotes the
    * candidate bottom points assuming that the track has monotonically
-   * _decreaing_ z position. top_lh_v are the candidate top points for an
+   * _decreasing_ z position. top_lh_v are the candidate top points for an
    * increasing z track, and top_hl_v are the candidate top points for a
    * decreasing z track.
    */

Core/include/Acts/Seeding/SeedFinderUtils.hpp

@@ -16,16 +16,13 @@
 namespace Acts {
 /// @brief A partial description of a circle in u-v space.
 struct LinCircle {
-  float Zo;
   float cotTheta;
   float iDeltaR;
   float Er;
   float U;
   float V;
   float x;
   float y;
-  float z;
-  float r;
 };
 
 /// @brief Transform two spacepoints to a u-v space circle.
@@ -57,7 +54,6 @@ LinCircle transformCoordinates(const external_spacepoint_t& sp,
 /// @param[in] spM The middle spacepoint.
 /// @param[in] bottom Should be true if vec are bottom spacepoints.
 /// @param[out] linCircleVec The output vector to write to.
-/// @returns Vector of sorted indexes for the vectors (vec and linCircleVec)
 template <typename external_spacepoint_t>
 void transformCoordinates(
     Acts::SpacePointData& spacePointData,
@@ -72,6 +68,12 @@ void transformCoordinates(Acts::SpacePointData& spacePointData,
                           std::vector<LinCircle>& linCircleVec,
                           callable_t&& extractFunction);
 
+/// @brief Fills LineCircle object for a SP dublet in u-v frame.
+///
+/// @param[in] lineCircleVariables Vector contaning LineCircle variables
+inline LinCircle fillLineCircle(
+    const std::array<float, 7>& lineCircleVariables);
+
 /// @brief Check the compatibility of spacepoint coordinates in xyz assuming the Bottom-Middle direction with the strip meassument details
 ///
 /// @tparam external_spacepoint_t The external spacepoint type.