Merge remote-tracking branch 'upstream/main' into SeedWriting

.zenodo.json

@@ -42,7 +42,7 @@
       "name": "Xiaocong Ai"
     },
     {
-      "affiliation": "Universitaet Freiburg",
+      "affiliation": "Universitaet Regensburg",
       "name": "Benjamin Huth"
     },
     {
@@ -113,6 +113,11 @@
       "name": "Stephen Nicholas Swatman",
       "orcid": "0000-0002-3747-3229"
     }
+    {
+      "affiliation": "CERN / TU Wien",
+      "name": "Felix Russo",
+      "orcid": "0009-0005-8975-2245"
+    }
   ],
   "access_right": "open",
   "related_identifiers": [

AUTHORS

@@ -15,6 +15,7 @@ The following people have contributed to the project (in alphabetical order):
 - Fabian Klimpel, CERN
 - Robert Langenberg, UMass
 - Alexander J. Pfleger, CERN, University of Graz
+- Felix Russo, CERN, TU Wien
 - Andreas Salzburger, CERN
 - Bastian Schlag, CERN, JGU Mainz
 - Andreas Stefl, CERN, TU Wien

CITATION.cff

@@ -86,6 +86,10 @@ authors:
   family-names: Swatman
   affiliation: CERN / University of Amsterdam
   orcid: https://orcid.org/0000-0002-3747-3229
+- given-names: Felix
+  family-names: Russo
+  affiliation: CERN / TU Wien
+  orcid: https://orcid.org/0009-0005-8975-2245
 version: 10.0.0
 date-released: 2021-07-28
 repository-code: https://github.com/acts-project/acts

Core/include/Acts/Detector/LayerStructureBuilder.hpp

@@ -86,6 +86,9 @@ class LayerStructureBuilder : public IInternalStructureBuilder {
     std::vector<ProtoSupport> supports = {};
     /// Definition of Binnings
     std::vector<ProtoBinning> binnings = {};
+    /// Minimum number of surfaces to build an internal structure
+    /// - otherwise the tryAll options is used
+    unsigned int nMinimalSurfaces = 4u;
     /// Polyhedron approximations
     unsigned int nSegments = 1u;
     /// Extra information, mainly for screen output

Core/include/Acts/Seeding/EstimateTrackParamsFromSeed.hpp

@@ -214,24 +214,32 @@ std::optional<BoundVector> estimateTrackParamsFromSeed(
   Vector3 local1 = transform.inverse() * spGlobalPositions[1];
   Vector3 local2 = transform.inverse() * spGlobalPositions[2];
 
-  // Lambda to transform the coordinates to the (u, v) space
-  auto uvTransform = [](const Vector3& local) -> Vector2 {
-    Vector2 uv;
-    ActsScalar denominator = local.x() * local.x() + local.y() * local.y();
-    uv.x() = local.x() / denominator;
-    uv.y() = local.y() / denominator;
-    return uv;
-  };
-  // The uv1.y() should be zero
-  Vector2 uv1 = uvTransform(local1);
-  Vector2 uv2 = uvTransform(local2);
+  // In the new frame the bottom sp is at the origin, while the middle
+  // sp in along the x axis. As such, the x-coordinate of the circle is
+  // at: x-middle / 2.
+  // The y coordinate can be found by using the straight line passing
+  // between the mid point between the middle and top sp and perpendicular to
+  // the line connecting them
+  Vector2 circleCenter;
+  circleCenter(0) = 0.5 * local1(0);
+
+  ActsScalar deltaX21 = local2(0) - local1(0);
+  ActsScalar sumX21 = local2(0) + local1(0);
+  // straight line connecting the two points
+  // y = a * x + c (we don't care about c right now)
+  // we simply need the slope
+  ActsScalar a = local2(1) / deltaX21;
+  // Perpendicular line is then y = -1/a *x + b
+  // we can evaluate b given we know a already by imposing
+  // the line passes through P = (0.5 * (x2 + x1), 0.5 * y2)
+  ActsScalar b = 0.5 * (local2(1) + 1. / a * sumX21);
+  circleCenter(1) = -1. / a * circleCenter(0) + b;
+  // Radius is distance between circleCenter and first sp, which is at (0, 0) in
+  // the new frame
+  // Sign depends on the slope a (positive vs negative)
+  int sign = a > 0 ? -1 : 1;
+  ActsScalar rho = sign / circleCenter.norm();
 
-  // A,B are slope and intercept of the straight line in the u,v plane
-  // connecting the three points
-  ActsScalar A = (uv2.y() - uv1.y()) / (uv2.x() - uv1.x());
-  ActsScalar B = uv2.y() - A * uv2.x();
-  // Curvature (with a sign) estimate
-  ActsScalar rho = -2.0 * B / std::hypot(1., A);
   // The projection of the top space point on the transverse plane of the new
   // frame
   ActsScalar rn = local2.x() * local2.x() + local2.y() * local2.y();
@@ -241,6 +249,7 @@ std::optional<BoundVector> estimateTrackParamsFromSeed(
       local2.z() * std::sqrt(1. / rn) / (1. + G * rho * rho * rn);
   // The momentum direction in the new frame (the center of the circle has the
   // coordinate (-1.*A/(2*B), 1./(2*B)))
+  ActsScalar A = -circleCenter(0) / circleCenter(1);
   Vector3 transDirection(1., A, std::hypot(1, A) * invTanTheta);
   // Transform it back to the original frame
   Vector3 direction = rotation * transDirection.normalized();

Core/include/Acts/Seeding/SeedConfirmationRangeConfig.hpp

@@ -15,6 +15,11 @@
 
 namespace Acts {
 
+// defaults experimental cuts to no operation in both seeding algorithms
+inline bool noopExperimentCuts(float /*bottomRadius*/, float /*cotTheta*/) {
+  return true;
+}
+
 /// @brief contains parameters for  seed confirmation
 struct SeedConfirmationRangeConfig {
   // z minimum and maximum of middle component of the seed used to define the

Core/include/Acts/Seeding/SeedFinder.ipp

@@ -206,7 +206,11 @@ SeedFinder<external_spacepoint_t, platform_t>::getCompatibleDoublets(
     const float deltaRMinSP, const float deltaRMaxSP, const float uIP,
     const float uIP2, const float cosPhiM, const float sinPhiM) const {
   float impactMax = m_config.impactMax;
-  if constexpr (candidateType == Acts::SpacePointCandidateType::eBottom) {
+
+  constexpr bool isBottomCandidate =
+      candidateType == Acts::SpacePointCandidateType::eBottom;
+
+  if constexpr (isBottomCandidate) {
     impactMax = -impactMax;
   }
 
@@ -252,7 +256,7 @@ SeedFinder<external_spacepoint_t, platform_t>::getCompatibleDoublets(
     // the iterator so we don't need to look at the other SPs again
     for (; min_itr != otherSPs.end(); ++min_itr) {
       const auto& otherSP = *min_itr;
-      if constexpr (candidateType == Acts::SpacePointCandidateType::eBottom) {
+      if constexpr (isBottomCandidate) {
         // if r-distance is too big, try next SP in bin
         if ((rM - otherSP->radius()) <= deltaRMaxSP) {
           break;
@@ -272,7 +276,7 @@ SeedFinder<external_spacepoint_t, platform_t>::getCompatibleDoublets(
     for (; min_itr != otherSPs.end(); ++min_itr) {
       const auto& otherSP = *min_itr;
 
-      if constexpr (candidateType == Acts::SpacePointCandidateType::eBottom) {
+      if constexpr (isBottomCandidate) {
         deltaR = (rM - otherSP->radius());
 
         // if r-distance is too small, try next SP in bin
@@ -288,7 +292,7 @@ SeedFinder<external_spacepoint_t, platform_t>::getCompatibleDoublets(
         }
       }
 
-      if constexpr (candidateType == Acts::SpacePointCandidateType::eBottom) {
+      if constexpr (isBottomCandidate) {
         deltaZ = (zM - otherSP->z());
       } else {
         deltaZ = (otherSP->z() - zM);
@@ -379,6 +383,14 @@ SeedFinder<external_spacepoint_t, platform_t>::getCompatibleDoublets(
         const float iDeltaR = std::sqrt(iDeltaR2);
         const float cotTheta = deltaZ * iDeltaR;
 
+        // discard bottom-middle dublets in a certain (r, eta) region according
+        // to detector specific cuts
+        if constexpr (isBottomCandidate) {
+          if (!m_config.experimentCuts(otherSP->radius(), cotTheta)) {
+            continue;
+          }
+        }
+
         const float Er =
             ((varianceZM + otherSP->varianceZ()) +
              (cotTheta * cotTheta) * (varianceRM + otherSP->varianceR())) *
@@ -420,6 +432,14 @@ SeedFinder<external_spacepoint_t, platform_t>::getCompatibleDoublets(
       const float iDeltaR = std::sqrt(iDeltaR2);
       const float cotTheta = deltaZ * iDeltaR;
 
+      // discard bottom-middle dublets in a certain (r, eta) region according
+      // to detector specific cuts
+      if constexpr (isBottomCandidate) {
+        if (!m_config.experimentCuts(otherSP->radius(), cotTheta)) {
+          continue;
+        }
+      }
+
       const float Er =
           ((varianceZM + otherSP->varianceZ()) +
            (cotTheta * cotTheta) * (varianceRM + otherSP->varianceR())) *

Core/include/Acts/Seeding/SeedFinderConfig.hpp

@@ -13,7 +13,9 @@
 #include "Acts/Seeding/SeedConfirmationRangeConfig.hpp"
 #include "Acts/Utilities/Delegate.hpp"
 
+#include <limits>
 #include <memory>
+#include <vector>
 
 namespace Acts {
 
@@ -171,6 +173,10 @@ struct SeedFinderConfig {
   // Returns position of the center of the top strip.
   Delegate<Acts::Vector3(const SpacePoint&)> getTopStripCenterPosition;
 
+  // Delegate to apply experiment specific cuts
+  Delegate<bool(float /*bottomRadius*/, float /*cotTheta*/)> experimentCuts{
+      DelegateFuncTag<&noopExperimentCuts>{}};
+
   bool isInInternalUnits = false;
 
   SeedFinderConfig toInternalUnits() const {

Core/include/Acts/Seeding/SeedFinderOrthogonal.ipp

@@ -226,6 +226,15 @@ bool SeedFinderOrthogonal<external_spacepoint_t>::validTuple(
     return false;
   }
 
+  /*
+   * Cut: Ensure that inner-middle dublet is in a certain (r, eta) region of the
+   * detector according to detector specific cuts.
+   */
+  const float rInner = (isMiddleInverted) ? rH : rL;
+  if (!m_config.experimentCuts(rInner, cotTheta)) {
+    return false;
+  }
+
   return true;
 }
 
@@ -453,7 +462,7 @@ void SeedFinderOrthogonal<external_spacepoint_t>::filterCandidates(
     }
 
     // continue if number of top SPs is smaller than minimum required for filter
-    if (top.size() < minCompatibleTopSPs) {
+    if (top_valid.size() < minCompatibleTopSPs) {
       continue;
     }
 

Core/include/Acts/Seeding/SeedFinderOrthogonalConfig.hpp

@@ -11,6 +11,7 @@
 #include "Acts/Definitions/Algebra.hpp"
 #include "Acts/Definitions/Units.hpp"
 #include "Acts/Seeding/SeedConfirmationRangeConfig.hpp"
+#include "Acts/Utilities/Delegate.hpp"
 
 #include <memory>
 
@@ -109,6 +110,10 @@ struct SeedFinderOrthogonalConfig {
   float highland = 0;
   float maxScatteringAngle2 = 0;
 
+  // Delegate to apply experiment specific cuts
+  Delegate<bool(float /*bottomRadius*/, float /*cotTheta*/)> experimentCuts{
+      DelegateFuncTag<&noopExperimentCuts>{}};
+
   bool isInInternalUnits = false;
 
   SeedFinderOrthogonalConfig calculateDerivedQuantities() const {

Core/include/Acts/TrackFinding/CombinatorialKalmanFilter.hpp

@@ -510,7 +510,7 @@ class CombinatorialKalmanFilter {
         // Return error if filtering finds no tracks
         if (result.lastTrackIndices.empty()) {
           // @TODO: Tracks like this should not be in the final output!
-          ACTS_WARNING("No tracks found");
+          ACTS_DEBUG("No tracks found");
           result.finished = true;
         } else {
           if (!smoothing) {

Core/include/Acts/TrackFitting/GaussianSumFitter.hpp

@@ -305,7 +305,7 @@ struct GaussianSumFitter {
       return return_error_or_abort(fwdResult.error());
     }
 
-    auto& fwdGsfResult =
+    const auto& fwdGsfResult =
         fwdResult->template get<typename GsfActor::result_type>();
 
     if (!fwdGsfResult.result.ok()) {
@@ -321,8 +321,8 @@ struct GaussianSumFitter {
     ACTS_VERBOSE("- processed states: " << fwdGsfResult.processedStates);
     ACTS_VERBOSE("- measurement states: " << fwdGsfResult.measurementStates);
 
-    std::size_t nInvalidBetheHeitler = fwdGsfResult.nInvalidBetheHeitler;
-    double maxPathXOverX0 = fwdGsfResult.maxPathXOverX0;
+    std::size_t nInvalidBetheHeitler = fwdGsfResult.nInvalidBetheHeitler.val();
+    double maxPathXOverX0 = fwdGsfResult.maxPathXOverX0.val();
 
     //////////////////
     // Backward pass
@@ -408,8 +408,14 @@ struct GaussianSumFitter {
           GsfError::NoMeasurementStatesCreatedBackward);
     }
 
-    nInvalidBetheHeitler += bwdGsfResult.nInvalidBetheHeitler;
-    maxPathXOverX0 = std::max(maxPathXOverX0, bwdGsfResult.maxPathXOverX0);
+    // For the backward pass we want the counters at in end (= at the
+    // interaction point) and not at the last measurement surface
+    bwdGsfResult.nInvalidBetheHeitler.update();
+    bwdGsfResult.maxPathXOverX0.update();
+    bwdGsfResult.sumPathXOverX0.update();
+    nInvalidBetheHeitler += bwdGsfResult.nInvalidBetheHeitler.val();
+    maxPathXOverX0 =
+        std::max(maxPathXOverX0, bwdGsfResult.maxPathXOverX0.val());
 
     if (nInvalidBetheHeitler > 0) {
       ACTS_WARNING("Encountered " << nInvalidBetheHeitler
@@ -479,12 +485,23 @@ struct GaussianSumFitter {
 
       if (trackContainer.hasColumn(
               hashString(GsfConstants::kFinalMultiComponentStateColumn))) {
-        ACTS_DEBUG("Add final multi-component state to track")
+        ACTS_DEBUG("Add final multi-component state to track");
         track.template component<GsfConstants::FinalMultiComponentState>(
             GsfConstants::kFinalMultiComponentStateColumn) = std::move(params);
       }
     }
 
+    if (trackContainer.hasColumn(
+            hashString(GsfConstants::kFwdMaxMaterialXOverX0))) {
+      track.template component<double>(GsfConstants::kFwdMaxMaterialXOverX0) =
+          fwdGsfResult.maxPathXOverX0.val();
+    }
+    if (trackContainer.hasColumn(
+            hashString(GsfConstants::kFwdSumMaterialXOverX0))) {
+      track.template component<double>(GsfConstants::kFwdSumMaterialXOverX0) =
+          fwdGsfResult.sumPathXOverX0.val();
+    }
+
     calculateTrackQuantities(track);
 
     return track;

Core/include/Acts/TrackFitting/GlobalChiSquareFitter.hpp

@@ -46,6 +46,10 @@
 namespace Acts {
 namespace Experimental {
 
+namespace Gx2fConstants {
+constexpr std::string_view gx2fnUpdateColumn = "Gx2fnUpdateColumn";
+}  // namespace Gx2fConstants
+
 /// Extension struct which holds delegates to customize the KF behavior
 template <typename traj_t>
 struct Gx2FitterExtensions {
@@ -789,7 +793,8 @@ class Gx2Fitter {
 
     ACTS_VERBOSE("final covariance:\n" << fullCovariancePredicted);
 
-    if (!trackContainer.hasColumn(Acts::hashString("Gx2fnUpdateColumn"))) {
+    if (!trackContainer.hasColumn(
+            Acts::hashString(Gx2fConstants::gx2fnUpdateColumn))) {
       trackContainer.template addColumn<std::size_t>("Gx2fnUpdateColumn");
     }
 
@@ -800,7 +805,8 @@ class Gx2Fitter {
     track.covariance() = fullCovariancePredicted;
     track.setReferenceSurface(params.referenceSurface().getSharedPtr());
 
-    if (trackContainer.hasColumn(Acts::hashString("Gx2fnUpdateColumn"))) {
+    if (trackContainer.hasColumn(
+            Acts::hashString(Gx2fConstants::gx2fnUpdateColumn))) {
       ACTS_DEBUG("Add nUpdate to track")
       track.template component<std::size_t>("Gx2fnUpdateColumn") = nUpdate;
     }

Core/include/Acts/TrackFitting/GsfOptions.hpp

@@ -39,6 +39,10 @@ constexpr std::string_view kFinalMultiComponentStateColumn =
     "gsf-final-multi-component-state";
 using FinalMultiComponentState =
     std::optional<Acts::MultiComponentBoundTrackParameters>;
+constexpr std::string_view kFwdSumMaterialXOverX0 =
+    "gsf-fwd-sum-material-x-over-x0";
+constexpr std::string_view kFwdMaxMaterialXOverX0 =
+    "gsf-fwd-max-material-x-over-x0";
 }  // namespace GsfConstants
 
 /// The extensions needed for the GSF