Merge branch 'main' into usingnamespacestd

.github/workflows/analysis.yml

@@ -63,9 +63,9 @@ jobs:
           && du -sh build
       - name: Coverage
         run: >
-          pip3 install gcovr==5.0
+          pip3 install gcovr==6.0
           && cd build
-          && /usr/bin/python3 ../CI/test_coverage
+          && /usr/bin/python3 ../CI/test_coverage.py
       - name: Upload coverage
         uses: codecov/codecov-action@v3
         with:

CI/test_coverage.py

@@ -0,0 +1,86 @@
+#!/usr/bin/env python
+import sys
+import os
+import subprocess
+import argparse
+import multiprocessing as mp
+import re
+
+
+if not os.path.exists("CMakeCache.txt"):
+    print("Not in CMake build dir. Not executing")
+    sys.exit(1)
+
+
+def check_output(*args, **kwargs):
+    p = subprocess.Popen(
+        *args, stdout=subprocess.PIPE, stderr=subprocess.STDOUT, **kwargs
+    )
+    p.wait()
+    stdout, stderr = p.communicate()
+    stdout = stdout.decode("utf-8")
+    return (p.returncode, stdout.strip())
+
+
+# call helper function
+def call(cmd):
+    print(" ".join(cmd))
+    try:
+        subprocess.check_call(cmd)
+    except subprocess.CalledProcessError as e:
+        print("Failed, output: ", e.output)
+        raise e
+
+
+p = argparse.ArgumentParser()
+p.add_argument("--gcov", default=check_output(["which", "gcov"])[1])
+args = p.parse_args()
+
+ret, gcovr_exe = check_output(["which", "gcovr"])
+assert ret == 0, "gcovr not installed. Use 'pip install gcovr'."
+
+ret, gcovr_version_text = check_output(["gcovr", "--version"])
+gcovr_version = tuple(
+    map(int, re.match("gcovr (\d+\.\d+)", gcovr_version_text).group(1).split("."))
+)
+
+extra_flags = []
+
+print(f"Found gcovr version {gcovr_version[0]}.{gcovr_version[1]}")
+if gcovr_version < (5,):
+    print("Consider upgrading to a newer gcovr version.")
+elif gcovr_version == (5, 1):
+    assert False and "Version 5.1 does not support parallel processing of gcov data"
+elif gcovr_version >= (6,):
+    extra_flags += ["--exclude-noncode-lines"]
+
+gcovr = [gcovr_exe]
+
+script_dir = os.path.dirname(__file__)
+source_dir = os.path.abspath(os.path.join(script_dir, ".."))
+coverage_dir = os.path.abspath("coverage")
+
+if not os.path.exists(coverage_dir):
+    os.makedirs(coverage_dir)
+
+excludes = ["-e", "../Tests/", "-e", ".*json\.hpp"]
+
+# create the html report
+call(
+    gcovr
+    + ["-r", source_dir]
+    + ["--gcov-executable", args.gcov]
+    + ["-j", str(mp.cpu_count())]
+    + excludes
+    + extra_flags
+    + ["--xml", "-o", "coverage/cov.xml"]
+)
+
+call(
+    gcovr
+    + ["-r", source_dir]
+    + ["-j", str(mp.cpu_count())]
+    + ["--gcov-executable", args.gcov]
+    + excludes
+    + extra_flags
+)

Core/include/Acts/TrackFitting/GlobalChiSquareFitter.hpp

@@ -103,6 +103,7 @@ struct Gx2FitterOptions {
   /// @param eLoss Whether to include energy loss
   /// @param freeToBoundCorrection_ Correction for non-linearity effect during transform from free to bound
   /// @param nUpdateMax_ Max number of iterations for updating the parameters
+  /// @param zeroField_ Disables the QoP fit in case of missing B-field
   Gx2FitterOptions(const GeometryContext& gctx,
                    const MagneticFieldContext& mctx,
                    std::reference_wrapper<const CalibrationContext> cctx,
@@ -112,7 +113,7 @@ struct Gx2FitterOptions {
                    bool eLoss = false,
                    const FreeToBoundCorrection& freeToBoundCorrection_ =
                        FreeToBoundCorrection(false),
-                   const size_t nUpdateMax_ = 5)
+                   const size_t nUpdateMax_ = 5, const bool zeroField_ = false)
       : geoContext(gctx),
         magFieldContext(mctx),
         calibrationContext(cctx),
@@ -122,7 +123,8 @@ struct Gx2FitterOptions {
         multipleScattering(mScattering),
         energyLoss(eLoss),
         freeToBoundCorrection(freeToBoundCorrection_),
-        nUpdateMax(nUpdateMax_) {}
+        nUpdateMax(nUpdateMax_),
+        zeroField(zeroField_) {}
 
   /// Contexts are required and the options must not be default-constructible.
   Gx2FitterOptions() = delete;
@@ -154,6 +156,9 @@ struct Gx2FitterOptions {
 
   /// Max number of iterations during the fit
   size_t nUpdateMax = 5;
+
+  /// Disables the QoP fit in case of missing B-field
+  bool zeroField = false;
 };
 
 template <typename traj_t>
@@ -614,15 +619,14 @@ class Gx2Fitter {
 
     using PropagatorOptions = Acts::PropagatorOptions<Actors, Aborters>;
 
-    const size_t reducedMatrixSize = 4;
     start_parameters_t params = sParameters;
     BoundVector deltaParams = BoundVector::Zero();
     double chi2sum = 0;
     BoundMatrix aMatrix = BoundMatrix::Zero();
     BoundVector bVector = BoundVector::Zero();
 
-    // Create a index of the 'tip' of the track stored in multitrajectory. It is
-    // needed outside the update loop. It will be updated with each iteration
+    // Create an index of the 'tip' of the track stored in multitrajectory. It
+    // is needed outside the update loop. It will be updated with each iteration
     // and used for the final track
     size_t tipIndex = Acts::MultiTrajectoryTraits::kInvalid;
 
@@ -700,13 +704,18 @@ class Gx2Fitter {
 
       // calculate delta params [a] * delta = b
       deltaParams = BoundVector::Zero();
-      const ActsVector<reducedMatrixSize> deltaParamsReduced =
-          aMatrix.topLeftCorner<reducedMatrixSize, reducedMatrixSize>()
-              .colPivHouseholderQr()
-              .solve(bVector.topLeftCorner<reducedMatrixSize, 1>());
-
-      for (size_t idp = 0; idp < reducedMatrixSize; idp++) {
-        deltaParams(idp, 0) = deltaParamsReduced(idp, 0);
+      if (gx2fOptions.zeroField) {
+        constexpr size_t reducedMatrixSize = 4;
+        deltaParams.topLeftCorner<reducedMatrixSize, 1>() =
+            aMatrix.topLeftCorner<reducedMatrixSize, reducedMatrixSize>()
+                .colPivHouseholderQr()
+                .solve(bVector.topLeftCorner<reducedMatrixSize, 1>());
+      } else {
+        constexpr size_t reducedMatrixSize = 5;
+        deltaParams.topLeftCorner<reducedMatrixSize, 1>() =
+            aMatrix.topLeftCorner<reducedMatrixSize, reducedMatrixSize>()
+                .colPivHouseholderQr()
+                .solve(bVector.topLeftCorner<reducedMatrixSize, 1>());
       }
 
       ACTS_VERBOSE("chi2sum = " << chi2sum);
@@ -728,17 +737,38 @@ class Gx2Fitter {
 
     // Calculate covariance of the fitted parameters with inverse of [a]
     BoundMatrix fullCovariancePredicted = BoundMatrix::Identity();
-    if (aMatrix.topLeftCorner<reducedMatrixSize, reducedMatrixSize>()
-            .determinant() != 0) {
-      fullCovariancePredicted
-          .template topLeftCorner<reducedMatrixSize, reducedMatrixSize>() =
-          aMatrix.topLeftCorner<reducedMatrixSize, reducedMatrixSize>()
-              .inverse();
-    } else if (gx2fOptions.nUpdateMax > 0) {
-      ACTS_ERROR("det(a) == 0. This should not happen ever.");
-      return Experimental::GlobalChiSquareFitterError::DetAIsZero;
+    bool aMatrixIsInvertible = false;
+    if (gx2fOptions.zeroField) {
+      constexpr size_t reducedMatrixSize = 4;
+
+      auto safeReducedCovariance = safeInverse(
+          aMatrix.topLeftCorner<reducedMatrixSize, reducedMatrixSize>().eval());
+      if (safeReducedCovariance) {
+        aMatrixIsInvertible = true;
+        fullCovariancePredicted
+            .topLeftCorner<reducedMatrixSize, reducedMatrixSize>() =
+            *safeReducedCovariance;
+      }
+    } else {
+      constexpr size_t reducedMatrixSize = 5;
+
+      auto safeReducedCovariance = safeInverse(
+          aMatrix.topLeftCorner<reducedMatrixSize, reducedMatrixSize>().eval());
+      if (safeReducedCovariance) {
+        aMatrixIsInvertible = true;
+        fullCovariancePredicted
+            .topLeftCorner<reducedMatrixSize, reducedMatrixSize>() =
+            *safeReducedCovariance;
+      }
     }
 
+    if (!aMatrixIsInvertible && gx2fOptions.nUpdateMax > 0) {
+      ACTS_ERROR("aMatrix is not invertible.");
+      return Experimental::GlobalChiSquareFitterError::AIsNotInvertible;
+    }
+
+    ACTS_VERBOSE("final covariance:\n" << fullCovariancePredicted);
+
     // Prepare track for return
     auto track = trackContainer.getTrack(trackContainer.addTrack());
     track.tipIndex() = tipIndex;

Core/include/Acts/TrackFitting/GlobalChiSquareFitterError.hpp

@@ -16,7 +16,7 @@ namespace Experimental {
 
 enum class GlobalChiSquareFitterError {
   // ensure all values are non-zero
-  DetAIsZero = 1,
+  AIsNotInvertible = 1,
 };
 
 std::error_code make_error_code(

Core/src/TrackFitting/GlobalChiSquareFitterError.cpp

@@ -24,8 +24,8 @@ class GlobalChiSquareFitterErrorCategory : public std::error_category {
     using Acts::Experimental::GlobalChiSquareFitterError;
 
     switch (static_cast<GlobalChiSquareFitterError>(c)) {
-      case GlobalChiSquareFitterError::DetAIsZero:
-        return "Gx2f has det(a) == 0 during update.";
+      case GlobalChiSquareFitterError::AIsNotInvertible:
+        return "Gx2f: aMatrix is not invertible.";
       default:
         return "unknown";
     }