Merge branch 'main' into detector-navigation-tests

.github/workflows/builds.yml

@@ -409,7 +409,7 @@ jobs:
           -DACTS_USE_SYSTEM_EIGEN3=OFF
           -DACTS_BUILD_EXAMPLES=ON
       - name: Build Examples
-        run: ${SETUP} && cmake --build build
+        run: ${SETUP} && ${PRELOAD} && cmake --build build
       - name: Install Examples
         run: ${SETUP} && cmake --build build --target install
       - name: Run Examples

Core/include/Acts/Propagator/EigenStepper.ipp

@@ -233,8 +233,41 @@ Acts::Result<double> Acts::EigenStepper<E, A>::step(
       return EigenStepperError::StepInvalid;
     }
 
-    // for moment, only update the transport part
-    state.stepping.jacTransport = D * state.stepping.jacTransport;
+    // See the documentation of Acts::blockedMult for a description of blocked
+    // matrix multiplication. However, we can go one further. Let's assume that
+    // some of these sub-matrices are zero matrices 0₈ and identity matrices
+    // I₈, namely:
+    //
+    // D₁₁ = I₈, J₁₁ = I₈, D₂₁ = 0₈, J₂₁ = 0₈
+    //
+    // Which gives:
+    //
+    // K₁₁ = I₈  * I₈  + D₁₂ * 0₈  = I₈
+    // K₁₂ = I₈  * J₁₂ + D₁₂ * J₂₂ = J₁₂ + D₁₂ * J₂₂
+    // K₂₁ = 0₈  * I₈  + D₂₂ * 0₈  = 0₈
+    // K₂₂ = 0₈  * J₁₂ + D₂₂ * J₂₂ = D₂₂ * J₂₂
+    //
+    // Furthermore, we're constructing K in place of J, and since
+    // K₁₁ = I₈ = J₁₁ and K₂₁ = 0₈ = D₂₁, we don't actually need to touch those
+    // sub-matrices at all!
+    if ((D.topLeftCorner<4, 4>().isIdentity()) &&
+        (D.bottomLeftCorner<4, 4>().isZero()) &&
+        (state.stepping.jacTransport.template topLeftCorner<4, 4>()
+             .isIdentity()) &&
+        (state.stepping.jacTransport.template bottomLeftCorner<4, 4>()
+             .isZero())) {
+      state.stepping.jacTransport.template topRightCorner<4, 4>() +=
+          D.topRightCorner<4, 4>() *
+          state.stepping.jacTransport.template bottomRightCorner<4, 4>();
+      state.stepping.jacTransport.template bottomRightCorner<4, 4>() =
+          (D.bottomRightCorner<4, 4>() *
+           state.stepping.jacTransport.template bottomRightCorner<4, 4>())
+              .eval();
+    } else {
+      // For safety purposes, we provide a full matrix multiplication as a
+      // backup strategy.
+      state.stepping.jacTransport = D * state.stepping.jacTransport;
+    }
   } else {
     if (!state.stepping.extension.finalize(state, *this, navigator, h)) {
       return EigenStepperError::StepInvalid;

Core/include/Acts/Seeding/BinnedGroup.hpp

@@ -40,6 +40,12 @@ class BinnedGroup {
               std::array<std::vector<std::size_t>, DIM> navigation =
                   std::array<std::vector<std::size_t>, DIM>());
 
+  BinnedGroup(grid_t&& grid, std::vector<bool> mask,
+              const Acts::GridBinFinder<DIM>& bottomFinder,
+              const Acts::GridBinFinder<DIM>& topFinder,
+              std::array<std::vector<std::size_t>, DIM> navigation =
+                  std::array<std::vector<std::size_t>, DIM>());
+
   BinnedGroup(grid_t& grid, const Acts::GridBinFinder<DIM>& bottomFinder,
               const Acts::GridBinFinder<DIM>& topFinder,
               std::array<std::vector<std::size_t>, DIM> navigation =
@@ -72,6 +78,11 @@ class BinnedGroup {
   /// @return Mutable reference to the stored grid
   grid_t& grid();
 
+  /// @brief Retrieve the mask
+  /// Only const accessor is supported
+  /// @return The mask
+  const std::vector<bool>& mask() const;
+
   /// @brief Get the begin iterator
   /// @return The iterator
   Acts::BinnedGroupIterator<grid_t> begin() const;
@@ -82,6 +93,9 @@ class BinnedGroup {
  private:
   /// @brief The N-dimentional grid
   grid_t m_grid;
+  /// @brief The mask to be applied to the grid. The size of this vector
+  /// corresponds to the global bins in the grid
+  std::vector<bool> m_mask{};
   /// @brief The Grid Bin Finder for bottom candidates
   const Acts::GridBinFinder<DIM>* m_bottomBinFinder{nullptr};
   /// @brief The Grid Bin Finder for top candidates

Core/include/Acts/Seeding/BinnedGroup.ipp

@@ -16,6 +16,7 @@ BinnedGroup<grid_t>::BinnedGroup(
     std::array<std::vector<std::size_t>, Acts::BinnedGroup<grid_t>::DIM>
         navigation)
     : m_grid(std::move(grid)),
+      m_mask(m_grid.size(true), true),
       m_bottomBinFinder(&bottomFinder),
       m_topBinFinder(&topFinder),
       m_bins(std::move(navigation)) {
@@ -31,6 +32,38 @@ BinnedGroup<grid_t>::BinnedGroup(
   }
 }
 
+template <typename grid_t>
+BinnedGroup<grid_t>::BinnedGroup(
+    grid_t&& grid, std::vector<bool> mask,
+    const Acts::GridBinFinder<Acts::BinnedGroup<grid_t>::DIM>& bottomFinder,
+    const Acts::GridBinFinder<Acts::BinnedGroup<grid_t>::DIM>& topFinder,
+    std::array<std::vector<std::size_t>, Acts::BinnedGroup<grid_t>::DIM>
+        navigation)
+    : m_grid(std::move(grid)),
+      m_mask(std::move(mask)),
+      m_bottomBinFinder(&bottomFinder),
+      m_topBinFinder(&topFinder),
+      m_bins(std::move(navigation)) {
+  // Check the elements in the mask corresponds to all the global bins in the
+  // grid so that we can check if a global bin is masked
+  if (m_mask.size() != m_grid.size(true)) {
+    throw std::invalid_argument(
+        "Provided mask does not match the grid. The number of entries must "
+        "correspond to the number of global bins in the grid.");
+  }
+
+  /// If navigation is not defined for all axes, then we default that to a
+  /// std::iota from 1ul
+  std::array<std::size_t, DIM> numLocBins = m_grid.numLocalBins();
+  for (std::size_t i(0ul); i < DIM; ++i) {
+    if (!m_bins[i].empty()) {
+      continue;
+    }
+    m_bins[i].resize(numLocBins[i]);
+    std::iota(m_bins[i].begin(), m_bins[i].end(), 1ul);
+  }
+}
+
 template <typename grid_t>
 const grid_t& BinnedGroup<grid_t>::grid() const {
   return m_grid;
@@ -41,6 +74,11 @@ grid_t& BinnedGroup<grid_t>::grid() {
   return m_grid;
 }
 
+template <typename grid_t>
+const std::vector<bool>& BinnedGroup<grid_t>::mask() const {
+  return m_mask;
+}
+
 template <typename grid_t>
 Acts::BinnedGroupIterator<grid_t> BinnedGroup<grid_t>::begin() const {
   return Acts::BinnedGroupIterator<grid_t>(

Core/include/Acts/Seeding/BinnedGroupIterator.ipp

@@ -84,7 +84,18 @@ void Acts::BinnedGroupIterator<grid_t>::findNotEmptyBin() {
   /// Iterate on the grid till we find a not-empty bin
   /// We start from the current bin configuration and move forward
   std::size_t dimCollection = (*m_gridItr).size();
-  while (dimCollection == 0ul && ++m_gridItr != m_gridItrEnd) {
+  // Check if the current global bin is masked. This only makes sense if
+  // we have not reached the end of the mask
+  bool passesMask = false;
+  if (m_gridItr != m_gridItrEnd) {
+    passesMask = m_group->mask().at(m_gridItr.globalBinIndex());
+  }
+  // loop and only stop when we find a non-empty bin which is not masked
+  while ((dimCollection == 0ul || !passesMask) && ++m_gridItr != m_gridItrEnd) {
     dimCollection = (*m_gridItr).size();
+    if (dimCollection == 0ul) {
+      continue;
+    }
+    passesMask = m_group->mask().at(m_gridItr.globalBinIndex());
   }
 }

Core/include/Acts/Seeding/SeedFinder.ipp

@@ -65,6 +65,10 @@ void SeedFinder<external_spacepoint_t, grid_t, platform_t>::createSeedsForGroup(
 
   // Get the middle space point candidates
   const auto& middleSPs = grid.at(middleSPsIdx);
+  // Return if somehow there are no middle sp candidates
+  if (middleSPs.size() == 0) {
+    return;
+  }
 
   // neighbours
   // clear previous results
@@ -83,6 +87,11 @@ void SeedFinder<external_spacepoint_t, grid_t, platform_t>::createSeedsForGroup(
         grid, idx, middleSPs.front()->radius() + m_config.deltaRMinTopSP);
   }
 
+  // Return if there are no bottom or top candidates
+  if (state.bottomNeighbours.size() == 0 || state.topNeighbours.size() == 0) {
+    return;
+  }
+
   for (const auto& spM : middleSPs) {
     float rM = spM->radius();
 

Core/include/Acts/Utilities/VectorHelpers.hpp

@@ -139,6 +139,9 @@ inline std::array<ActsScalar, 4> evaluateTrigonomics(const Vector3& direction) {
   // can be turned into cosine/sine
   const ActsScalar cosTheta = z;
   const ActsScalar sinTheta = std::sqrt(1 - z * z);
+  assert(sinTheta != 0 &&
+         "VectorHelpers: Vector is parallel to the z-axis "
+         "which leads to division by zero");
   const ActsScalar invSinTheta = 1. / sinTheta;
   const ActsScalar cosPhi = x * invSinTheta;
   const ActsScalar sinPhi = y * invSinTheta;

Core/include/Acts/Vertexing/AdaptiveMultiVertexFinder.hpp

@@ -112,9 +112,13 @@ class AdaptiveMultiVertexFinder final : public IVertexFinder {
     // Include also single track vertices
     bool addSingleTrackVertices = false;
 
-    // Use 3d information for evaluating the vertex distance significance
-    // for vertex merging/splitting
-    bool do3dSplitting = false;
+    // If doFullSplitting == true, we check the 3D distance (if useTime ==
+    // false) or the 4D distance (if useTime == true) of the vertices to
+    // determine whether they are merged.
+    // If doFullSplitting == false, we check the z distance (if useTime ==
+    // false) or the z-t distance (if useTime == true) of the vertices to
+    // determine whether they are merged.
+    bool doFullSplitting = false;
 
     // Maximum vertex contamination value
     double maximumVertexContamination = 0.5;
@@ -344,18 +348,19 @@ class AdaptiveMultiVertexFinder final : public IVertexFinder {
   /// @param fitterState The vertex fitter state
   ///
   /// @return Keep new vertex
-  bool keepNewVertex(Vertex& vtx, const std::vector<Vertex*>& allVertices,
-                     VertexFitterState& fitterState) const;
+  Result<bool> keepNewVertex(Vertex& vtx,
+                             const std::vector<Vertex*>& allVertices,
+                             VertexFitterState& fitterState) const;
 
   /// @brief Method that evaluates if the new vertex candidate is
   /// merged with one of the previously found vertices
   ///
   /// @param vtx The vertex candidate
-  /// @param allVertices All so far found vertices
+  /// @param allVertices All vertices that were found so far
   ///
-  /// @return Vertex is merged
-  bool isMergedVertex(const Vertex& vtx,
-                      const std::vector<Vertex*>& allVertices) const;
+  /// @return Bool indicating whether the vertex is merged
+  Result<bool> isMergedVertex(const Vertex& vtx,
+                              const std::vector<Vertex*>& allVertices) const;
 
   /// @brief Method that deletes last vertex from list of all vertices
   /// and refits all vertices afterwards

Core/include/Acts/Vertexing/VertexingError.hpp

@@ -21,6 +21,9 @@ enum class VertexingError {
   NotConverged,
   ElementNotFound,
   NoCovariance,
+  SingularMatrix,
+  NonPositiveVariance,
+  MatrixNotPositiveDefinite,
   InvalidInput,
 };
 

Core/src/Vertexing/AdaptiveMultiVertexFinder.cpp

@@ -125,10 +125,12 @@ Acts::Result<std::vector<Acts::Vertex>> AdaptiveMultiVertexFinder::find(
         break;
       }
     }
-    // MARK: fpeMaskBegin(FLTUND, 1, #2590)
-    bool keepVertex = isGoodVertex &&
-                      keepNewVertex(vtxCandidate, allVerticesPtr, fitterState);
-    // MARK: fpeMaskEnd(FLTUND)
+    auto keepNewVertexResult =
+        keepNewVertex(vtxCandidate, allVerticesPtr, fitterState);
+    if (!keepNewVertexResult.ok()) {
+      return keepNewVertexResult.error();
+    }
+    bool keepVertex = isGoodVertex && *keepNewVertexResult;
     ACTS_DEBUG("New vertex will be saved: " << keepVertex);
 
     // Delete vertex from allVertices list again if it's not kept
@@ -464,7 +466,7 @@ bool AdaptiveMultiVertexFinder::removeTrackIfIncompatible(
   return true;
 }
 
-bool AdaptiveMultiVertexFinder::keepNewVertex(
+Result<bool> AdaptiveMultiVertexFinder::keepNewVertex(
     Vertex& vtx, const std::vector<Vertex*>& allVertices,
     VertexFitterState& fitterState) const {
   double contamination = 0.;
@@ -475,23 +477,30 @@ bool AdaptiveMultiVertexFinder::keepNewVertex(
         fitterState.tracksAtVerticesMap.at(std::make_pair(trk, &vtx));
     double trackWeight = trkAtVtx.trackWeight;
     contaminationNum += trackWeight * (1. - trackWeight);
+    // MARK: fpeMaskBegin(FLTUND, 1, #2590)
     contaminationDeNom += trackWeight * trackWeight;
+    // MARK: fpeMaskEnd(FLTUND)
   }
   if (contaminationDeNom != 0) {
     contamination = contaminationNum / contaminationDeNom;
   }
   if (contamination > m_cfg.maximumVertexContamination) {
-    return false;
+    return Result<bool>::success(false);
   }
 
-  if (isMergedVertex(vtx, allVertices)) {
-    return false;
+  auto isMergedVertexResult = isMergedVertex(vtx, allVertices);
+  if (!isMergedVertexResult.ok()) {
+    return Result<bool>::failure(isMergedVertexResult.error());
   }
 
-  return true;
+  if (*isMergedVertexResult) {
+    return Result<bool>::success(false);
+  }
+
+  return Result<bool>::success(true);
 }
 
-bool AdaptiveMultiVertexFinder::isMergedVertex(
+Result<bool> AdaptiveMultiVertexFinder::isMergedVertex(
     const Vertex& vtx, const std::vector<Vertex*>& allVertices) const {
   const Vector4& candidatePos = vtx.fullPosition();
   const SquareMatrix4& candidateCov = vtx.fullCovariance();
@@ -504,24 +513,36 @@ bool AdaptiveMultiVertexFinder::isMergedVertex(
     const SquareMatrix4& otherCov = otherVtx->fullCovariance();
 
     double significance = 0;
-    if (!m_cfg.do3dSplitting) {
-      const double deltaZPos = otherPos[eZ] - candidatePos[eZ];
-      const double sumVarZ = otherCov(eZ, eZ) + candidateCov(eZ, eZ);
-      if (sumVarZ <= 0) {
-        // TODO FIXME this should never happen
-        continue;
+    if (!m_cfg.doFullSplitting) {
+      if (m_cfg.useTime) {
+        const Vector2 deltaZT = otherPos.tail<2>() - candidatePos.tail<2>();
+        SquareMatrix2 sumCovZT = candidateCov.bottomRightCorner<2, 2>() +
+                                 otherCov.bottomRightCorner<2, 2>();
+        auto sumCovZTInverse = safeInverse(sumCovZT);
+        if (!sumCovZTInverse) {
+          ACTS_ERROR("Vertex z-t covariance matrix is singular.");
+          return Result<bool>::failure(VertexingError::SingularMatrix);
+        }
+        significance = std::sqrt(deltaZT.dot(*sumCovZTInverse * deltaZT));
+      } else {
+        const double deltaZPos = otherPos[eZ] - candidatePos[eZ];
+        const double sumVarZ = otherCov(eZ, eZ) + candidateCov(eZ, eZ);
+        if (sumVarZ <= 0) {
+          ACTS_ERROR("Variance of the vertex's z-coordinate is not positive.");
+          return Result<bool>::failure(VertexingError::SingularMatrix);
+        }
+        // Use only z significance
+        significance = std::abs(deltaZPos) / std::sqrt(sumVarZ);
       }
-      // Use only z significance
-      significance = std::abs(deltaZPos) / std::sqrt(sumVarZ);
     } else {
       if (m_cfg.useTime) {
         // Use 4D information for significance
         const Vector4 deltaPos = otherPos - candidatePos;
         SquareMatrix4 sumCov = candidateCov + otherCov;
         auto sumCovInverse = safeInverse(sumCov);
         if (!sumCovInverse) {
-          // TODO FIXME this should never happen
-          continue;
+          ACTS_ERROR("Vertex 4D covariance matrix is singular.");
+          return Result<bool>::failure(VertexingError::SingularMatrix);
         }
         significance = std::sqrt(deltaPos.dot(*sumCovInverse * deltaPos));
       } else {
@@ -531,17 +552,25 @@ bool AdaptiveMultiVertexFinder::isMergedVertex(
             candidateCov.topLeftCorner<3, 3>() + otherCov.topLeftCorner<3, 3>();
         auto sumCovInverse = safeInverse(sumCov);
         if (!sumCovInverse) {
-          // TODO FIXME this should never happen
-          continue;
+          ACTS_ERROR("Vertex 3D covariance matrix is singular.");
+          return Result<bool>::failure(VertexingError::SingularMatrix);
         }
         significance = std::sqrt(deltaPos.dot(*sumCovInverse * deltaPos));
       }
     }
+    if (significance < 0.) {
+      ACTS_ERROR(
+          "Found a negative significance (i.e., a negative chi2) when checking "
+          "if vertices are merged. This should never happen since the vertex "
+          "covariance should be positive definite, and thus its inverse "
+          "should be positive definite as well.");
+      return Result<bool>::failure(VertexingError::MatrixNotPositiveDefinite);
+    }
     if (significance < m_cfg.maxMergeVertexSignificance) {
-      return true;
+      return Result<bool>::success(true);
     }
   }
-  return false;
+  return Result<bool>::success(false);
 }
 
 Acts::Result<void> AdaptiveMultiVertexFinder::deleteLastVertex(