refactor!: compile-time DetectorMeasurementInfo for SSS

Core/include/Acts/Seeding/SeedFinder.hpp

@@ -34,6 +34,8 @@ namespace Acts {
 
 enum class SpacePointCandidateType : short { BOTTOM, TOP };
 
+enum class DetectorMeasurementInfo : short { DEFAULT, DETAILED };
+
 template <typename external_spacepoint_t, typename platform_t = void*>
 class SeedFinder {
   ///////////////////////////////////////////////////////////////////
@@ -159,6 +161,7 @@ class SeedFinder {
   /// @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
+  template <Acts::DetectorMeasurementInfo detailedMeasurement>
   void filterCandidates(Acts::SpacePointData& spacePointData,
                         const InternalSpacePoint<external_spacepoint_t>& SpM,
                         const Acts::SeedFinderOptions& options,

Core/include/Acts/Seeding/SeedFinder.ipp

@@ -172,8 +172,13 @@ void SeedFinder<external_spacepoint_t, platform_t>::createSeedsForGroup(
     }
 
     // filter candidates
-    filterCandidates(state.spacePointData, *spM.get(), options, seedFilterState,
-                     state);
+    if (m_config.useDetailedDoubleMeasurementInfo) {
+      filterCandidates<Acts::DetectorMeasurementInfo::DETAILED>(
+          state.spacePointData, *spM.get(), options, seedFilterState, state);
+    } else {
+      filterCandidates<Acts::DetectorMeasurementInfo::DEFAULT>(
+          state.spacePointData, *spM.get(), options, seedFilterState, state);
+    }
 
     m_config.seedFilter->filterSeeds_1SpFixed(
         state.spacePointData, state.candidates_collector,
@@ -367,12 +372,15 @@ SeedFinder<external_spacepoint_t, platform_t>::getCompatibleDoublets(
 }
 
 template <typename external_spacepoint_t, typename platform_t>
+template <Acts::DetectorMeasurementInfo detailedMeasurement>
 inline void SeedFinder<external_spacepoint_t, platform_t>::filterCandidates(
     Acts::SpacePointData& spacePointData,
     const InternalSpacePoint<external_spacepoint_t>& spM,
     const Acts::SeedFinderOptions& options, SeedFilterState& seedFilterState,
     SeedingState& state) const {
   float rM = spM.radius();
+  float cosPhiM = spM.x() / rM;
+  float sinPhiM = spM.y() / rM;
   float varianceRM = spM.varianceR();
   float varianceZM = spM.varianceZ();
 
@@ -389,17 +397,19 @@ inline void SeedFinder<external_spacepoint_t, platform_t>::filterCandidates(
     sorted_tops[i] = i;
   }
 
-  std::sort(sorted_bottoms.begin(), sorted_bottoms.end(),
-            [&state](const std::size_t& a, const std::size_t& b) -> bool {
-              return state.linCircleBottom[a].cotTheta <
-                     state.linCircleBottom[b].cotTheta;
-            });
-
-  std::sort(sorted_tops.begin(), sorted_tops.end(),
-            [&state](const std::size_t& a, const std::size_t& b) -> bool {
-              return state.linCircleTop[a].cotTheta <
-                     state.linCircleTop[b].cotTheta;
-            });
+  if constexpr (detailedMeasurement == Acts::DetectorMeasurementInfo::DEFAULT) {
+    std::sort(sorted_bottoms.begin(), sorted_bottoms.end(),
+              [&state](const std::size_t& a, const std::size_t& b) -> bool {
+                return state.linCircleBottom[a].cotTheta <
+                       state.linCircleBottom[b].cotTheta;
+              });
+
+    std::sort(sorted_tops.begin(), sorted_tops.end(),
+              [&state](const std::size_t& a, const std::size_t& b) -> bool {
+                return state.linCircleTop[a].cotTheta <
+                       state.linCircleTop[b].cotTheta;
+              });
+  }
 
   // Reserve enough space, in case current capacity is too little
   state.topSpVec.reserve(numTopSP);
@@ -449,9 +459,10 @@ inline void SeedFinder<external_spacepoint_t, platform_t>::filterCandidates(
     // coordinate transformation and checks for middle spacepoint
     // x and y terms for the rotation from UV to XY plane
     float rotationTermsUVtoXY[2] = {0, 0};
-    if (m_config.useDetailedDoubleMeasurementInfo) {
-      rotationTermsUVtoXY[0] = spM.x() * sinTheta / spM.radius();
-      rotationTermsUVtoXY[1] = spM.y() * sinTheta / spM.radius();
+    if constexpr (detailedMeasurement ==
+                  Acts::DetectorMeasurementInfo::DETAILED) {
+      rotationTermsUVtoXY[0] = cosPhiM * sinTheta;
+      rotationTermsUVtoXY[1] = sinPhiM * sinTheta;
     }
 
     // minimum number of compatible top SPs to trigger the filter for a certain
@@ -477,8 +488,16 @@ inline void SeedFinder<external_spacepoint_t, platform_t>::filterCandidates(
       float vb = 0.;
       float ut = 0.;
       float vt = 0.;
-
-      if (m_config.useDetailedDoubleMeasurementInfo) {
+      double rMTransf[3];
+      float xB = 0.;
+      float yB = 0.;
+      float xT = 0.;
+      float yT = 0.;
+      float iDeltaRB2 = 0.;
+      float iDeltaRT2 = 0.;
+
+      if constexpr (detailedMeasurement ==
+                    Acts::DetectorMeasurementInfo::DETAILED) {
         // protects against division by 0
         float dU = lt.U - Ub;
         if (dU == 0.) {
@@ -488,14 +507,15 @@ inline void SeedFinder<external_spacepoint_t, platform_t>::filterCandidates(
         // x_0 and y_0
         float A0 = (lt.V - Vb) / dU;
 
+        float zPositionMiddle = cosTheta * std::sqrt(1 + A0 * A0);
+
         // position of Middle SP converted from UV to XY assuming cotTheta
         // evaluated from the Bottom and Middle SPs double
         double positionMiddle[3] = {
             rotationTermsUVtoXY[0] - rotationTermsUVtoXY[1] * A0,
             rotationTermsUVtoXY[0] * A0 + rotationTermsUVtoXY[1],
-            cosTheta * std::sqrt(1 + A0 * A0)};
+            zPositionMiddle};
 
-        double rMTransf[3];
         if (!xyzCoordinateCheck(spacePointData, m_config, spM, positionMiddle,
                                 rMTransf)) {
           continue;
@@ -508,7 +528,7 @@ inline void SeedFinder<external_spacepoint_t, platform_t>::filterCandidates(
         double positionBottom[3] = {
             rotationTermsUVtoXY[0] * Cb - rotationTermsUVtoXY[1] * Sb,
             rotationTermsUVtoXY[0] * Sb + rotationTermsUVtoXY[1] * Cb,
-            cosTheta * std::sqrt(1 + A0 * A0)};
+            zPositionMiddle};
 
         auto spB = state.compatBottomSP[b];
         double rBTransf[3];
@@ -523,7 +543,7 @@ inline void SeedFinder<external_spacepoint_t, platform_t>::filterCandidates(
         double positionTop[3] = {
             rotationTermsUVtoXY[0] * Ct - rotationTermsUVtoXY[1] * St,
             rotationTermsUVtoXY[0] * St + rotationTermsUVtoXY[1] * Ct,
-            cosTheta * std::sqrt(1 + A0 * A0)};
+            zPositionMiddle};
 
         auto spT = state.compatTopSP[t];
         double rTTransf[3];
@@ -533,32 +553,24 @@ inline void SeedFinder<external_spacepoint_t, platform_t>::filterCandidates(
         }
 
         // bottom and top coordinates in the spM reference frame
-        float xB = rBTransf[0] - rMTransf[0];
-        float yB = rBTransf[1] - rMTransf[1];
+        xB = rBTransf[0] - rMTransf[0];
+        yB = rBTransf[1] - rMTransf[1];
         float zB = rBTransf[2] - rMTransf[2];
-        float xT = rTTransf[0] - rMTransf[0];
-        float yT = rTTransf[1] - rMTransf[1];
+        xT = rTTransf[0] - rMTransf[0];
+        yT = rTTransf[1] - rMTransf[1];
         float zT = rTTransf[2] - rMTransf[2];
 
-        float iDeltaRB2 = 1. / (xB * xB + yB * yB);
-        float iDeltaRT2 = 1. / (xT * xT + yT * yT);
+        iDeltaRB2 = 1. / (xB * xB + yB * yB);
+        iDeltaRT2 = 1. / (xT * xT + yT * yT);
 
         cotThetaB = -zB * std::sqrt(iDeltaRB2);
         cotThetaT = zT * std::sqrt(iDeltaRT2);
-
-        rMxy = std::sqrt(rMTransf[0] * rMTransf[0] + rMTransf[1] * rMTransf[1]);
-        float Ax = rMTransf[0] / rMxy;
-        float Ay = rMTransf[1] / rMxy;
-
-        ub = (xB * Ax + yB * Ay) * iDeltaRB2;
-        vb = (yB * Ax - xB * Ay) * iDeltaRB2;
-        ut = (xT * Ax + yT * Ay) * iDeltaRT2;
-        vt = (yT * Ax - xT * Ay) * iDeltaRT2;
       }
 
       // use geometric average
       float cotThetaAvg2 = cotThetaB * cotThetaT;
-      if (m_config.arithmeticAverageCotTheta) {
+      if constexpr (detailedMeasurement ==
+                    Acts::DetectorMeasurementInfo::DETAILED) {
         // use arithmetic average
         float averageCotTheta = 0.5 * (cotThetaB + cotThetaT);
         cotThetaAvg2 = averageCotTheta * averageCotTheta;
@@ -587,7 +599,8 @@ inline void SeedFinder<external_spacepoint_t, platform_t>::filterCandidates(
       // fair for scattering and measurement uncertainties)
       if (deltaCotTheta2 > (error2 + scatteringInRegion2)) {
         // skip top SPs based on cotTheta sorting when producing triplets
-        if (not m_config.skipPreviousTopSP) {
+        if constexpr (detailedMeasurement ==
+                      Acts::DetectorMeasurementInfo::DETAILED) {
           continue;
         }
         // break if cotTheta from bottom SP < cotTheta from top SP because
@@ -599,13 +612,27 @@ inline void SeedFinder<external_spacepoint_t, platform_t>::filterCandidates(
         continue;
       }
 
+      if constexpr (detailedMeasurement ==
+                    Acts::DetectorMeasurementInfo::DETAILED) {
+        rMxy = std::sqrt(rMTransf[0] * rMTransf[0] + rMTransf[1] * rMTransf[1]);
+        double irMxy = 1 / rMxy;
+        float Ax = rMTransf[0] * irMxy;
+        float Ay = rMTransf[1] * irMxy;
+
+        ub = (xB * Ax + yB * Ay) * iDeltaRB2;
+        vb = (yB * Ax - xB * Ay) * iDeltaRB2;
+        ut = (xT * Ax + yT * Ay) * iDeltaRT2;
+        vt = (yT * Ax - xT * Ay) * iDeltaRT2;
+      }
+
       float dU = 0;
       float A = 0;
       float S2 = 0;
       float B = 0;
       float B2 = 0;
 
-      if (m_config.useDetailedDoubleMeasurementInfo) {
+      if constexpr (detailedMeasurement ==
+                    Acts::DetectorMeasurementInfo::DETAILED) {
         dU = ut - ub;
         // protects against division by 0
         if (dU == 0.) {
@@ -656,7 +683,8 @@ inline void SeedFinder<external_spacepoint_t, platform_t>::filterCandidates(
 
       // if deltaTheta larger than allowed scattering for calculated pT, skip
       if (deltaCotTheta2 > (error2 + p2scatterSigma)) {
-        if (not m_config.skipPreviousTopSP) {
+        if constexpr (detailedMeasurement ==
+                      Acts::DetectorMeasurementInfo::DETAILED) {
           continue;
         }
         if (cotThetaB - cotThetaT < 0) {
@@ -668,9 +696,13 @@ inline void SeedFinder<external_spacepoint_t, platform_t>::filterCandidates(
       // A and B allow calculation of impact params in U/V plane with linear
       // function
       // (in contrast to having to solve a quadratic function in x/y plane)
-      float Im = m_config.useDetailedDoubleMeasurementInfo
-                     ? std::abs((A - B * rMxy) * rMxy)
-                     : std::abs((A - B * rM) * rM);
+      float Im = 0;
+      if constexpr (detailedMeasurement ==
+                    Acts::DetectorMeasurementInfo::DETAILED) {
+        Im = std::abs((A - B * rMxy) * rMxy);
+      } else {
+        Im = std::abs((A - B * rM) * rM);
+      }
 
       if (Im > m_config.impactMax) {
         continue;

Core/include/Acts/Seeding/SeedFinderConfig.hpp

@@ -64,16 +64,9 @@ struct SeedFinderConfig {
   // enable cut on the compatibility between interaction point and SPs
   bool interactionPointCut = false;
 
-  // use arithmetic average in the calculation of the squared error on the
-  // difference in tan(theta)
-  bool arithmeticAverageCotTheta = false;
-
   // non equidistant binning in z
   std::vector<float> zBinEdges;
 
-  // skip top SPs based on cotTheta sorting when producing triplets
-  bool skipPreviousTopSP = false;
-
   // seed confirmation
   bool seedConfirmation = false;
   // parameters for central seed confirmation

Examples/Python/python/acts/examples/itk.py

@@ -387,7 +387,6 @@ def itkSeedingAlgConfig(inputSpacePointsType: InputSpacePointsType):
         deltaRMaxBottomSP = 150 * u.mm
         deltaZMax = float("inf") * u.mm
         interactionPointCut = True
-        arithmeticAverageCotTheta = False
         impactMax = 2 * u.mm
         zBinsCustomLooping = [
             1,
@@ -402,7 +401,6 @@ def itkSeedingAlgConfig(inputSpacePointsType: InputSpacePointsType):
             5,
             7,
         ]  # enable custom z looping when searching for SPs, must contain numbers from 1 to the total number of bin in zBinEdges
-        skipPreviousTopSP = True
         zBinNeighborsTop = [
             [0, 0],
             [-1, 0],
@@ -451,10 +449,8 @@ def itkSeedingAlgConfig(inputSpacePointsType: InputSpacePointsType):
         deltaRMaxBottomSP = deltaRMaxTopSP
         deltaZMax = 900 * u.mm
         interactionPointCut = False
-        arithmeticAverageCotTheta = True
         impactMax = 20 * u.mm
         zBinsCustomLooping = [6, 7, 5, 8, 4, 9, 3, 10, 2, 11, 1]
-        skipPreviousTopSP = False
         zBinNeighborsTop = [
             [0, 0],
             [-1, 0],
@@ -503,11 +499,9 @@ def itkSeedingAlgConfig(inputSpacePointsType: InputSpacePointsType):
         impactMax=impactMax,
         deltaPhiMax=deltaPhiMax,
         interactionPointCut=interactionPointCut,
-        arithmeticAverageCotTheta=arithmeticAverageCotTheta,
         deltaZMax=deltaZMax,
         maxPtScattering=maxPtScattering,
         zBinEdges=zBinEdges,
-        skipPreviousTopSP=skipPreviousTopSP,
         zBinsCustomLooping=zBinsCustomLooping,
         rRangeMiddleSP=rRangeMiddleSP,
         useVariableMiddleSPRange=useVariableMiddleSPRange,

Examples/Python/python/acts/examples/reconstruction.py

@@ -36,11 +36,9 @@
         "impactMax",
         "deltaPhiMax",
         "interactionPointCut",
-        "arithmeticAverageCotTheta",
         "deltaZMax",
         "maxPtScattering",
         "zBinEdges",
-        "skipPreviousTopSP",
         "zBinsCustomLooping",
         "rRangeMiddleSP",
         "useVariableMiddleSPRange",
@@ -57,7 +55,7 @@
         "z",  # (min,max)
         "zOutermostLayers",  # (min,max)
     ],
-    defaults=[None] * 20 + [(None, None)] * 8,
+    defaults=[None] * 18 + [(None, None)] * 8,
 )
 SeedFinderOptionsArg = namedtuple(
     "SeedFinderOptions", ["beamPos", "bFieldInZ"], defaults=[(None, None), None]
@@ -212,7 +210,7 @@ def addSeeding(
     initialVarInflation : list
         List of 6 scale factors to inflate the initial covariance matrix
         Defaults (all 1) specified in Examples/Algorithms/TruthTracking/ActsExamples/TruthTracking/ParticleSmearing.hpp
-    seedFinderConfigArg : SeedFinderConfigArg(maxSeedsPerSpM, cotThetaMax, sigmaScattering, radLengthPerSeed, minPt, impactMax, deltaPhiMax, interactionPointCut, arithmeticAverageCotTheta, deltaZMax, maxPtScattering, zBinEdges, skipPreviousTopSP, zBinsCustomLooping, rRangeMiddleSP, useVariableMiddleSPRange, binSizeR, seedConfirmation, centralSeedConfirmationRange, forwardSeedConfirmationRange, deltaR, deltaRBottomSP, deltaRTopSP, deltaRMiddleSPRange, collisionRegion, r, z, zOutermostLayers)
+    seedFinderConfigArg : SeedFinderConfigArg(maxSeedsPerSpM, cotThetaMax, sigmaScattering, radLengthPerSeed, minPt, impactMax, deltaPhiMax, interactionPointCut, deltaZMax, maxPtScattering, zBinEdges, zBinsCustomLooping, rRangeMiddleSP, useVariableMiddleSPRange, binSizeR, seedConfirmation, centralSeedConfirmationRange, forwardSeedConfirmationRange, deltaR, deltaRBottomSP, deltaRTopSP, deltaRMiddleSPRange, collisionRegion, r, z, zOutermostLayers)
         SeedFinderConfig settings. deltaR, deltaRBottomSP, deltaRTopSP, deltaRMiddleSPRange, collisionRegion, r, z, zOutermostLayers are ranges specified as a tuple of (min,max). beamPos is specified as (x,y).
         Defaults specified in Core/include/Acts/Seeding/SeedFinderConfig.hpp
     seedFinderOptionsArg :  SeedFinderOptionsArg(bFieldInZ, beamPos)
@@ -546,11 +544,9 @@ def addStandardSeeding(
             minPt=seedFinderConfigArg.minPt,
             impactMax=seedFinderConfigArg.impactMax,
             interactionPointCut=seedFinderConfigArg.interactionPointCut,
-            arithmeticAverageCotTheta=seedFinderConfigArg.arithmeticAverageCotTheta,
             deltaZMax=seedFinderConfigArg.deltaZMax,
             maxPtScattering=seedFinderConfigArg.maxPtScattering,
             zBinEdges=seedFinderConfigArg.zBinEdges,
-            skipPreviousTopSP=seedFinderConfigArg.skipPreviousTopSP,
             zBinsCustomLooping=seedFinderConfigArg.zBinsCustomLooping,
             rRangeMiddleSP=seedFinderConfigArg.rRangeMiddleSP,
             useVariableMiddleSPRange=seedFinderConfigArg.useVariableMiddleSPRange,

Examples/Python/src/TrackFinding.cpp

@@ -98,7 +98,6 @@ void addTrackFinding(Context& ctx) {
     ACTS_PYTHON_MEMBER(impactMax);
     ACTS_PYTHON_MEMBER(deltaZMax);
     ACTS_PYTHON_MEMBER(zBinEdges);
-    ACTS_PYTHON_MEMBER(skipPreviousTopSP);
     ACTS_PYTHON_MEMBER(interactionPointCut);
     ACTS_PYTHON_MEMBER(zBinsCustomLooping);
     ACTS_PYTHON_MEMBER(useVariableMiddleSPRange);
@@ -111,7 +110,6 @@ void addTrackFinding(Context& ctx) {
     ACTS_PYTHON_MEMBER(seedConfirmation);
     ACTS_PYTHON_MEMBER(centralSeedConfirmationRange);
     ACTS_PYTHON_MEMBER(forwardSeedConfirmationRange);
-    ACTS_PYTHON_MEMBER(arithmeticAverageCotTheta);
     ACTS_PYTHON_MEMBER(useDetailedDoubleMeasurementInfo);
     ACTS_PYTHON_STRUCT_END();
     patchKwargsConstructor(c);
@@ -153,7 +151,6 @@ void addTrackFinding(Context& ctx) {
     ACTS_PYTHON_MEMBER(rMin);
     ACTS_PYTHON_MEMBER(radLengthPerSeed);
     ACTS_PYTHON_MEMBER(deltaZMax);
-    ACTS_PYTHON_MEMBER(skipPreviousTopSP);
     ACTS_PYTHON_MEMBER(interactionPointCut);
     ACTS_PYTHON_MEMBER(deltaPhiMax);
     ACTS_PYTHON_MEMBER(highland);