feat!: Updates to GBTS for Athena Implementation (#3882)

Changes to GBTS seeding in Core and Examples needed to implement in Athena and make it work with any space point type.

Core/include/Acts/Seeding/GbtsDataStorage.hpp

@@ -29,20 +29,20 @@ struct GbtsSP {
   const space_point_t *SP;  // want inside to have pointer
   int gbtsID;
   int combined_ID;
-  GbtsSP(const space_point_t *sp, int id, int combined_id)
-      : SP(sp), gbtsID(id), combined_ID{combined_id} {
-    if (SP->sourceLinks().size() == 1) {  // pixels have 1 SL
-      m_isPixel = true;
-    } else {
-      m_isPixel = false;
-    }
-    m_phi = std::atan(SP->x() / SP->y());
+  float m_phi;
+  float m_r;
+  float m_ClusterWidth;
+  GbtsSP(const space_point_t *sp, int id, int combined_id, float ClusterWidth)
+      : SP(sp),
+        gbtsID(id),
+        combined_ID{combined_id},
+        m_ClusterWidth(ClusterWidth) {
+    m_phi = std::atan2(SP->y(), SP->x());
+    m_r = std::sqrt((SP->x() * SP->x()) + (SP->y() * SP->y()));
   };
-  bool isPixel() const { return m_isPixel; }
-  bool isSCT() const { return !m_isPixel; }
   float phi() const { return m_phi; }
-  bool m_isPixel;
-  float m_phi;
+  float r() const { return m_r; }
+  float ClusterWidth() const { return m_ClusterWidth; }
 };
 
 template <typename space_point_t>
@@ -152,7 +152,7 @@ class GbtsDataStorage {
       return -1;
     }
 
-    int binIndex = pL->getEtaBin(sp.SP->z(), sp.SP->r());
+    int binIndex = pL->getEtaBin(sp.SP->z(), sp.r());
 
     if (binIndex == -1) {
       return -2;
@@ -165,7 +165,7 @@ class GbtsDataStorage {
       float max_tau = 100.0;
       // can't do this bit yet as dont have cluster width
       if (useClusterWidth) {
-        float cluster_width = 1;  // temporary while cluster width not available
+        float cluster_width = sp.ClusterWidth();
         min_tau = 6.7 * (cluster_width - 0.2);
         max_tau =
             1.6 + 0.15 / (cluster_width + 0.2) + 6.1 * (cluster_width - 0.2);
@@ -176,7 +176,7 @@ class GbtsDataStorage {
               sp, min_tau, max_tau));  // adding ftf member to nodes
     } else {
       if (useClusterWidth) {
-        float cluster_width = 1;  // temporary while cluster width not available
+        float cluster_width = sp.ClusterWidth();
         if (cluster_width > 0.2) {
           return -3;
         }

Core/include/Acts/Seeding/GbtsTrackingFilter.hpp

@@ -52,7 +52,7 @@ struct GbtsEdgeState {
     //  x' =  x*m_c + y*m_s
     //  y' = -x*m_s + y*m_c
 
-    m_refY = pS->m_n2->m_spGbts.SP->r();
+    m_refY = pS->m_n2->m_spGbts.r();
     m_refX =
         pS->m_n2->m_spGbts.SP->x() * m_c + pS->m_n2->m_spGbts.SP->y() * m_s;
 
@@ -67,7 +67,7 @@ struct GbtsEdgeState {
 
     m_Y[0] = pS->m_n2->m_spGbts.SP->z();
     m_Y[1] = (pS->m_n1->m_spGbts.SP->z() - pS->m_n2->m_spGbts.SP->z()) /
-             (pS->m_n1->m_spGbts.SP->r() - pS->m_n2->m_spGbts.SP->r());
+             (pS->m_n1->m_spGbts.r() - pS->m_n2->m_spGbts.r());
 
     memset(&m_Cx[0][0], 0, sizeof(m_Cx));
     memset(&m_Cy[0][0], 0, sizeof(m_Cy));
@@ -276,7 +276,7 @@ class GbtsTrackingFilter {
     x = pS->m_n1->m_spGbts.SP->x();
     y = pS->m_n1->m_spGbts.SP->y();
     z = pS->m_n1->m_spGbts.SP->z();
-    r = pS->m_n1->m_spGbts.SP->r();
+    r = pS->m_n1->m_spGbts.r();
 
     refX = x * ts.m_c + y * ts.m_s;
     mx = -x * ts.m_s + y * ts.m_c;  // measured X[0]

Core/include/Acts/Seeding/SeedFinderGbts.ipp

@@ -10,7 +10,6 @@
 // TODO: update to C++17 style
 
 #include "Acts/Geometry/Extent.hpp"
-#include "Acts/Seeding/SeedFilter.hpp"
 #include "Acts/Seeding/SeedFinder.hpp"
 #include "Acts/Seeding/SeedFinderGbtsConfig.hpp"
 #include "Acts/Seeding/SeedFinderUtils.hpp"
@@ -47,10 +46,6 @@ void SeedFinderGbts<external_spacepoint_t>::loadSpacePoints(
     const std::vector<GbtsSP<external_spacepoint_t>>& gbtsSPvect) {
   ACTS_VERBOSE("Loading space points");
   for (const auto& gbtssp : gbtsSPvect) {
-    bool is_Pixel = gbtssp.isPixel();
-    if (!is_Pixel) {
-      continue;
-    }
     m_storage->addSpacePoint(gbtssp, (m_config.m_useClusterWidth > 0));
   }
 
@@ -166,11 +161,11 @@ void SeedFinderGbts<external_spacepoint_t>::runGbts_TrackFinder(
                 continue;
               }
 
-              float r1 = n1->m_spGbts.SP->r();
+              float r1 = n1->m_spGbts.r();
               float x1 = n1->m_spGbts.SP->x();
               float y1 = n1->m_spGbts.SP->y();
               float z1 = n1->m_spGbts.SP->z();
-              float phi1 = std::atan(x1 / y1);
+              float phi1 = n1->m_spGbts.phi();
 
               float minPhi = phi1 - deltaPhi;
               float maxPhi = phi1 + deltaPhi;
@@ -199,7 +194,7 @@ void SeedFinderGbts<external_spacepoint_t>::runGbts_TrackFinder(
                   continue;
                 }
 
-                float r2 = n2->m_spGbts.SP->r();
+                float r2 = n2->m_spGbts.r();
 
                 float dr = r2 - r1;
 
@@ -545,7 +540,7 @@ void SeedFinderGbts<external_spacepoint_t>::runGbts_TrackFinder(
 
     for (unsigned int idx_m = 1; idx_m < vSP.size() - 1; idx_m++) {
       const GbtsSP<external_spacepoint_t>& spM = *vSP.at(idx_m);
-      const double pS_r = spM.SP->r();
+      const double pS_r = spM.r();
       const double pS_x = spM.SP->x();
       const double pS_y = spM.SP->y();
       const double cosA = pS_x / pS_r;
@@ -656,31 +651,32 @@ void SeedFinderGbts<external_spacepoint_t>::createSeeds(
     return;
   }
 
-  m_triplets.clear();  // member of class , saying not declared, maybe public?
+  m_triplets.clear();
 
   for (auto& track : vTracks) {
     for (auto& seed : track.m_seeds) {  // access member of GbtsTrigTracklet
-
-      float newQ = seed.Q();  // function of TrigInDetTriplet
-      if (m_config.m_LRTmode) {
-        // In LRT mode penalize pixels in Triplets
-        if (seed.s1().isPixel()) {
-          newQ += 1000;  // functions of TrigSiSpacePointBase
-        }
-        if (seed.s2().isPixel()) {
-          newQ += 1000;
-        }
-        if (seed.s3().isPixel()) {
-          newQ += 1000;
-        }
-      } else {
-        // In normal (non LRT) mode penalise SSS by 1000, PSS (if enabled) and
-        // PPS by 10000
-        if (seed.s3().isSCT()) {
-          newQ += seed.s1().isSCT() ? 1000.0 : 10000.0;
-        }
-      }
-      seed.Q(newQ);
+      // Currently not used, but leaving in to use concept in future development
+      //  float newQ = seed.Q();  // function of TrigInDetTriplet
+      //  if (m_config.m_LRTmode) {
+      //    // In LRT mode penalize pixels in Triplets
+      //    if (seed.s1().isPixel()) {
+      //      newQ += 1000;  // functions of TrigSiSpacePointBase
+      //    }
+      //    if (seed.s2().isPixel()) {
+      //      newQ += 1000;
+      //    }
+      //    if (seed.s3().isPixel()) {
+      //      newQ += 1000;
+      //    }
+      //  } else {
+      //    // In normal (non LRT) mode penalise SSS by 1000, PSS (if enabled)
+      //    and
+      //    // PPS by 10000
+      //    if (seed.s3().isSCT()) {
+      //      newQ += seed.s1().isSCT() ? 1000.0 : 10000.0;
+      //    }
+      //  }
+      //  seed.Q(newQ);
       m_triplets.emplace_back(seed);
     }
   }

Core/include/Acts/Seeding/SeedFinderGbtsConfig.hpp

@@ -20,9 +20,6 @@
 // core algorithm so in acts namespace
 namespace Acts {
 
-template <typename T>
-class SeedFilter;
-
 template <typename SpacePoint>
 struct SeedFinderGbtsConfig {
   // // how many sigmas of scattering angle should be considered?
@@ -31,16 +28,12 @@ struct SeedFinderGbtsConfig {
   // Seed cut
   float minPt = 400. * Acts::UnitConstants::MeV;
 
-  ///////////some declared not filled in by reco: //////
-  std::shared_ptr<Acts::SeedFilter<SpacePoint>> seedFilter;
-
   //   //detector ROI
   //   // derived values, set on SeedFinder construction
   float highland = 0;
   float maxScatteringAngle2 = 0;
   //   bool isInInternalUnits = false;
   /// for load space points
-  unsigned int maxSeedsPerSpM = 5;
 
   // Parameter which can loosen the tolerance of the track seed to form a
   // helix. This is useful for e.g. misaligned seeding.
@@ -50,11 +43,11 @@ struct SeedFinderGbtsConfig {
   float m_nMaxPhiSlice = 53;  // used to calculate phi slices
   bool m_useClusterWidth =
       false;  // bool for use of cluster width in loadSpacePoints function
-  std::string connector_input_file;  // input file for connector object
+  std::string ConnectorInputFile;  // Path to the connector configuration file
+                                   // that defines the layer connections
   std::vector<TrigInDetSiLayer> m_layerGeometry;
 
   // for runGbts_TrackFinder
-  bool m_LRTmode = true;  // eventually want to set from full chain
   bool m_useEtaBinning =
       true;  // bool to use eta binning from geometry structure
   bool m_doubletFilterRZ = true;  // bool applies new Z cuts on doublets

Examples/Algorithms/TrackFinding/include/ActsExamples/TrackFinding/GbtsSeedingAlgorithm.hpp

@@ -10,7 +10,6 @@
 
 #pragma once
 
-#include "Acts/Seeding/SeedFilterConfig.hpp"
 #include "Acts/Seeding/SeedFinderGbts.hpp"
 #include "Acts/Seeding/SeedFinderGbtsConfig.hpp"
 #include "ActsExamples/EventData/Cluster.hpp"
@@ -36,7 +35,6 @@ class GbtsSeedingAlgorithm final : public IAlgorithm {
 
     std::string outputSeeds;
 
-    Acts::SeedFilterConfig seedFilterConfig;
     Acts::SeedFinderGbtsConfig<SimSpacePoint> seedFinderConfig;
     Acts::SeedFinderOptions seedFinderOptions;
 

Examples/Algorithms/TrackFinding/src/GbtsSeedingAlgorithm.cpp

@@ -9,7 +9,6 @@
 #include "ActsExamples/TrackFinding/GbtsSeedingAlgorithm.hpp"
 
 #include "Acts/Geometry/GeometryIdentifier.hpp"
-#include "Acts/Seeding/SeedFilter.hpp"
 #include "ActsExamples/EventData/IndexSourceLink.hpp"
 #include "ActsExamples/EventData/Measurement.hpp"
 #include "ActsExamples/EventData/ProtoTrack.hpp"
@@ -39,8 +38,6 @@ ActsExamples::GbtsSeedingAlgorithm::GbtsSeedingAlgorithm(
   // fill config struct
   m_cfg.layerMappingFile = m_cfg.layerMappingFile;
 
-  m_cfg.seedFilterConfig = m_cfg.seedFilterConfig.toInternalUnits();
-
   m_cfg.seedFinderConfig =
       m_cfg.seedFinderConfig.toInternalUnits().calculateDerivedQuantities();
 
@@ -70,7 +67,7 @@ ActsExamples::GbtsSeedingAlgorithm::GbtsSeedingAlgorithm(
   m_cfg.seedFinderConfig.m_layerGeometry = LayerNumbering();
 
   std::ifstream input_ifstream(
-      m_cfg.seedFinderConfig.connector_input_file.c_str(), std::ifstream::in);
+      m_cfg.seedFinderConfig.ConnectorInputFile.c_str(), std::ifstream::in);
 
   // connector
   std::unique_ptr<Acts::GbtsConnector> inputConnector =
@@ -89,16 +86,14 @@ ActsExamples::ProcessCode ActsExamples::GbtsSeedingAlgorithm::execute(
       MakeGbtsSpacePoints(ctx, m_cfg.ActsGbtsMap);
 
   for (auto sp : GbtsSpacePoints) {
-    ACTS_DEBUG("Gbts space points:  Gbts_id: "
-               << sp.gbtsID << " z: " << sp.SP->z() << " r: " << sp.SP->r()
-               << " ACTS volume:  "
-               << sp.SP->sourceLinks()
-                      .front()
-                      .get<IndexSourceLink>()
-                      .geometryId()
-                      .volume());
+    const auto &links = sp.SP->sourceLinks();
+    if (!links.empty()) {
+      ACTS_DEBUG("Gbts space points:  Gbts_id: "
+                 << sp.gbtsID << " z: " << sp.SP->z() << " r: " << sp.SP->r()
+                 << " ACTS volume:  "
+                 << links.front().get<IndexSourceLink>().geometryId().volume());
+    }
   }
-
   // this is now calling on a core algorithm
   Acts::SeedFinderGbts<SimSpacePoint> finder(
       m_cfg.seedFinderConfig, *m_gbtsGeo,
@@ -222,8 +217,11 @@ ActsExamples::GbtsSeedingAlgorithm::MakeGbtsSpacePoints(
       int eta_mod = Find->second.second;
       int combined_id = Gbts_id * 1000 + eta_mod;
 
+      float ClusterWidth =
+          0;  // false input as this is not available in examples
       // fill Gbts vector with current sapce point and ID
-      gbtsSpacePoints.emplace_back(&spacePoint, Gbts_id, combined_id);
+      gbtsSpacePoints.emplace_back(&spacePoint, Gbts_id, combined_id,
+                                   ClusterWidth);  // make new GbtsSP here !
     }
   }
   ACTS_VERBOSE("Space points successfully assigned Gbts ID");

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

@@ -278,7 +278,7 @@ def addSeeding(
     field: acts.MagneticFieldProvider,
     geoSelectionConfigFile: Optional[Union[Path, str]] = None,
     layerMappingConfigFile: Optional[Union[Path, str]] = None,
-    connector_inputConfigFile: Optional[Union[Path, str]] = None,
+    ConnectorInputConfigFile: Optional[Union[Path, str]] = None,
     seedingAlgorithm: SeedingAlgorithm = SeedingAlgorithm.Default,
     trackSmearingSigmas: TrackSmearingSigmas = TrackSmearingSigmas(),
     initialSigmas: Optional[list] = None,
@@ -430,12 +430,11 @@ def addSeeding(
                 spacePoints,
                 seedFinderConfigArg,
                 seedFinderOptionsArg,
-                seedFilterConfigArg,
                 trackingGeometry,
                 logLevel,
                 layerMappingConfigFile,
                 geoSelectionConfigFile,
-                connector_inputConfigFile,
+                ConnectorInputConfigFile,
             )
         elif seedingAlgorithm == SeedingAlgorithm.Hashing:
             logger.info("Using Hashing seeding")
@@ -1088,24 +1087,22 @@ def addGbtsSeeding(
     spacePoints: str,
     seedFinderConfigArg: SeedFinderConfigArg,
     seedFinderOptionsArg: SeedFinderOptionsArg,
-    seedFilterConfigArg: SeedFilterConfigArg,
     trackingGeometry: acts.TrackingGeometry,
     logLevel: acts.logging.Level = None,
     layerMappingConfigFile: Union[Path, str] = None,
     geoSelectionConfigFile: Union[Path, str] = None,
-    connector_inputConfigFile: Union[Path, str] = None,
+    ConnectorInputConfigFile: Union[Path, str] = None,
 ):
     """Gbts seeding"""
 
     logLevel = acts.examples.defaultLogging(sequence, logLevel)()
     layerMappingFile = str(layerMappingConfigFile)  # turn path into string
-    connector_inputFile = str(connector_inputConfigFile)
+    ConnectorInputFileStr = str(ConnectorInputConfigFile)
     seedFinderConfig = acts.SeedFinderGbtsConfig(
         **acts.examples.defaultKWArgs(
             sigmaScattering=seedFinderConfigArg.sigmaScattering,
-            maxSeedsPerSpM=seedFinderConfigArg.maxSeedsPerSpM,
             minPt=seedFinderConfigArg.minPt,
-            connector_input_file=connector_inputFile,
+            ConnectorInputFile=ConnectorInputFileStr,
             m_useClusterWidth=False,
         ),
     )
@@ -1121,33 +1118,11 @@ def addGbtsSeeding(
             bFieldInZ=seedFinderOptionsArg.bFieldInZ,
         )
     )
-    seedFilterConfig = acts.SeedFilterConfig(
-        **acts.examples.defaultKWArgs(
-            maxSeedsPerSpM=seedFinderConfig.maxSeedsPerSpM,
-            deltaRMin=(
-                seedFinderConfigArg.deltaR[0]
-                if seedFilterConfigArg.deltaRMin is None
-                else seedFilterConfigArg.deltaRMin
-            ),
-            impactWeightFactor=seedFilterConfigArg.impactWeightFactor,
-            zOriginWeightFactor=seedFilterConfigArg.zOriginWeightFactor,
-            compatSeedWeight=seedFilterConfigArg.compatSeedWeight,
-            compatSeedLimit=seedFilterConfigArg.compatSeedLimit,
-            numSeedIncrement=seedFilterConfigArg.numSeedIncrement,
-            seedWeightIncrement=seedFilterConfigArg.seedWeightIncrement,
-            seedConfirmation=seedFilterConfigArg.seedConfirmation,
-            # curvatureSortingInFilter=seedFilterConfigArg.curvatureSortingInFilter,
-            maxSeedsPerSpMConf=seedFilterConfigArg.maxSeedsPerSpMConf,
-            maxQualitySeedsPerSpMConf=seedFilterConfigArg.maxQualitySeedsPerSpMConf,
-            useDeltaRorTopRadius=seedFilterConfigArg.useDeltaRorTopRadius,
-        )
-    )
 
     seedingAlg = acts.examples.GbtsSeedingAlgorithm(
         level=logLevel,
         inputSpacePoints=[spacePoints],
         outputSeeds="seeds",
-        seedFilterConfig=seedFilterConfig,
         seedFinderConfig=seedFinderConfig,
         seedFinderOptions=seedFinderOptions,
         layerMappingFile=layerMappingFile,