refactor!: compile-time DetectorMeasurementInfo for SSS

Core/include/Acts/Seeding/SeedFinder.hpp

@@ -34,6 +34,8 @@ namespace Acts {
 
 enum class SpacePointCandidateType : short { eBottom, eTop };
 
+enum class DetectorMeasurementInfo : short { eDefault, eDetailed };
+
 template <typename external_spacepoint_t, typename platform_t = void*>
 class SeedFinder {
   ///////////////////////////////////////////////////////////////////
@@ -164,6 +166,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

@@ -177,8 +177,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::eDetailed>(
+          state.spacePointData, *spM.get(), options, seedFilterState, state);
+    } else {
+      filterCandidates<Acts::DetectorMeasurementInfo::eDefault>(
+          state.spacePointData, *spM.get(), options, seedFilterState, state);
+    }
 
     m_config.seedFilter->filterSeeds_1SpFixed(
         state.spacePointData, state.candidates_collector,
@@ -431,6 +436,7 @@ 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,
@@ -455,7 +461,8 @@ inline void SeedFinder<external_spacepoint_t, platform_t>::filterCandidates(
     sorted_tops[i] = i;
   }
 
-  if (not m_config.useDetailedDoubleMeasurementInfo) {
+  if constexpr (detailedMeasurement ==
+                Acts::DetectorMeasurementInfo::eDefault) {
     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 <
@@ -517,7 +524,8 @@ 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) {
+    if constexpr (detailedMeasurement ==
+                  Acts::DetectorMeasurementInfo::eDetailed) {
       rotationTermsUVtoXY[0] = cosPhiM * sinTheta;
       rotationTermsUVtoXY[1] = sinPhiM * sinTheta;
     }
@@ -545,8 +553,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::eDetailed) {
         // protects against division by 0
         float dU = lt.U - Ub;
         if (dU == 0.) {
@@ -565,7 +581,6 @@ inline void SeedFinder<external_spacepoint_t, platform_t>::filterCandidates(
             rotationTermsUVtoXY[0] * A0 + rotationTermsUVtoXY[1],
             zPositionMiddle};
 
-        double rMTransf[3];
         if (!xyzCoordinateCheck(spacePointData, m_config, spM, positionMiddle,
                                 rMTransf)) {
           continue;
@@ -603,33 +618,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]);
-        double irMxy = 1 / rMxy;
-        float Ax = static_cast<float>(rMTransf[0] * irMxy);
-        float Ay = static_cast<float>(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;
       }
 
       // use geometric average
       float cotThetaAvg2 = cotThetaB * cotThetaT;
-      if (m_config.arithmeticAverageCotTheta) {
+      if constexpr (detailedMeasurement ==
+                    Acts::DetectorMeasurementInfo::eDetailed) {
         // use arithmetic average
         float averageCotTheta = 0.5 * (cotThetaB + cotThetaT);
         cotThetaAvg2 = averageCotTheta * averageCotTheta;
@@ -658,7 +664,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 (m_config.useDetailedDoubleMeasurementInfo) {
+        if constexpr (detailedMeasurement ==
+                      Acts::DetectorMeasurementInfo::eDetailed) {
           continue;
         }
         // break if cotTheta from bottom SP < cotTheta from top SP because
@@ -670,13 +677,27 @@ inline void SeedFinder<external_spacepoint_t, platform_t>::filterCandidates(
         continue;
       }
 
+      if constexpr (detailedMeasurement ==
+                    Acts::DetectorMeasurementInfo::eDetailed) {
+        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::eDetailed) {
         dU = ut - ub;
         // protects against division by 0
         if (dU == 0.) {
@@ -729,7 +750,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 (m_config.useDetailedDoubleMeasurementInfo) {
+        if constexpr (detailedMeasurement ==
+                      Acts::DetectorMeasurementInfo::eDetailed) {
           continue;
         }
         if (cotThetaB - cotThetaT < 0) {
@@ -741,9 +763,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::eDetailed) {
+        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

@@ -68,10 +68,6 @@ 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;
 

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

@@ -390,7 +390,6 @@ def itkSeedingAlgConfig(
         deltaRMaxBottomSP = 150 * u.mm
         deltaZMax = float("inf") * u.mm
         interactionPointCut = True
-        arithmeticAverageCotTheta = False
         impactMax = 2 * u.mm
         zBinsCustomLooping = [
             1,
@@ -461,7 +460,6 @@ def itkSeedingAlgConfig(
         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]
         zBinNeighborsTop = [
@@ -533,7 +531,6 @@ def itkSeedingAlgConfig(
         impactMax=impactMax,
         deltaPhiMax=deltaPhiMax,
         interactionPointCut=interactionPointCut,
-        arithmeticAverageCotTheta=arithmeticAverageCotTheta,
         deltaZMax=deltaZMax,
         maxPtScattering=maxPtScattering,
         zBinEdges=zBinEdges,

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

@@ -36,7 +36,6 @@
         "impactMax",
         "deltaPhiMax",
         "interactionPointCut",
-        "arithmeticAverageCotTheta",
         "deltaZMax",
         "maxPtScattering",
         "zBinEdges",
@@ -57,7 +56,7 @@
         "z",  # (min,max)
         "zOutermostLayers",  # (min,max)
     ],
-    defaults=[None] * 20 + [(None, None)] * 8,
+    defaults=[None] * 19 + [(None, None)] * 8,
 )
 SeedFinderOptionsArg = namedtuple(
     "SeedFinderOptions", ["beamPos", "bFieldInZ"], defaults=[(None, None), None]
@@ -202,8 +201,10 @@ def addSeeding(
         Defaults specified in Examples/Algorithms/TrackFinding/include/ActsExamples/TrackFinding/TrackParamsEstimationAlgorithm.hpp
     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, deltaZMax, maxPtScattering, zBinEdges, zBinsCustomLooping, rRangeMiddleSP, useVariableMiddleSPRange, binSizeR, seedConfirmation, centralSeedConfirmationRange, forwardSeedConfirmationRange, deltaR, deltaRBottomSP, deltaRTopSP, deltaRMiddleSPRange, collisionRegion, r, z, zOutermostLayers)
         Defaults specified in Examples/Algorithms/TrackFinding/include/ActsExamples/TrackFinding/TrackParamsEstimationAlgorithm.hpp
-    seedFinderConfigArg : SeedFinderConfigArg(maxSeedsPerSpM, cotThetaMax, sigmaScattering, radLengthPerSeed, minPt, impactMax, deltaPhiMax, interactionPointCut, arithmeticAverageCotTheta, deltaZMax, maxPtScattering, zBinEdges, zBinsCustomLooping, skipZMiddleBinSearch, 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, skipZMiddleBinSearch, 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)
@@ -541,7 +542,6 @@ def addStandardSeeding(
             minPt=seedFinderConfigArg.minPt,
             impactMax=seedFinderConfigArg.impactMax,
             interactionPointCut=seedFinderConfigArg.interactionPointCut,
-            arithmeticAverageCotTheta=seedFinderConfigArg.arithmeticAverageCotTheta,
             deltaZMax=seedFinderConfigArg.deltaZMax,
             maxPtScattering=seedFinderConfigArg.maxPtScattering,
             zBinEdges=seedFinderConfigArg.zBinEdges,

Examples/Python/src/TrackFinding.cpp

@@ -134,7 +134,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);