refactor(gx2f): make material work for |eta| > 2

Core/include/Acts/TrackFitting/GlobalChiSquareFitter.hpp

@@ -392,10 +392,8 @@ void addMeasurementToGx2fSums(Gx2fSystem& extendedSystem,
   }
 
   // Create an extended Jacobian. This one contains only eBoundSize rows,
-  // because the rest is irrelevant. We fill it in the next steps
+  // because the rest is irrelevant. We fill it in the next steps.
   // TODO make dimsExtendedParams template with unrolling
-
-  // We create an empty Jacobian and fill it in the next steps
   Eigen::MatrixXd extendedJacobian =
       Eigen::MatrixXd::Zero(eBoundSize, extendedSystem.nDims());
 
@@ -1246,7 +1244,7 @@ class Gx2Fitter {
 
       auto& gx2fActor = propagatorOptions.actorList.template get<GX2FActor>();
       gx2fActor.inputMeasurements = &inputMeasurements;
-      gx2fActor.multipleScattering = multipleScattering;
+      gx2fActor.multipleScattering = false;
       gx2fActor.extensions = gx2fOptions.extensions;
       gx2fActor.calibrationContext = &gx2fOptions.calibrationContext.get();
       gx2fActor.actorLogger = m_actorLogger.get();
@@ -1301,38 +1299,8 @@ class Gx2Fitter {
       track.tipIndex() = tipIndex;
       track.linkForward();
 
-      // Count the material surfaces, to set up the system. In the multiple
-      // scattering case, we need to extend our system.
-      std::size_t nMaterialSurfaces = 0;
-      if (multipleScattering) {
-        ACTS_DEBUG("Count the valid material surfaces.");
-        for (const auto& trackState : track.trackStates()) {
-          const auto typeFlags = trackState.typeFlags();
-          const bool stateHasMaterial =
-              typeFlags.test(TrackStateFlag::MaterialFlag);
-
-          if (!stateHasMaterial) {
-            continue;
-          }
-
-          // Get and store geoId for the current material surface
-          const GeometryIdentifier geoId =
-              trackState.referenceSurface().geometryId();
-
-          const auto scatteringMapId = scatteringMap.find(geoId);
-          assert(scatteringMapId != scatteringMap.end() &&
-                 "No scattering angles found for material surface.");
-          if (!scatteringMapId->second.materialIsValid()) {
-            continue;
-          }
-
-          nMaterialSurfaces++;
-        }
-      }
-
-      // We need 6 dimensions for the bound parameters and 2 * nMaterialSurfaces
-      // dimensions for the scattering angles.
-      const std::size_t dimsExtendedParams = eBoundSize + 2 * nMaterialSurfaces;
+      // We need 6 dimensions for the bound parameters.
+      const std::size_t dimsExtendedParams = eBoundSize;
 
       // System that we fill with the information gathered by the actor and
       // evaluate later
@@ -1344,8 +1312,8 @@ class Gx2Fitter {
       // all stored material in each propagation.
       std::vector<GeometryIdentifier> geoIdVector;
 
-      fillGx2fSystem(track, extendedSystem, multipleScattering, scatteringMap,
-                     geoIdVector, *m_addToSumLogger);
+      fillGx2fSystem(track, extendedSystem, false, scatteringMap, geoIdVector,
+                     *m_addToSumLogger);
 
       chi2sum = extendedSystem.chi2();
 
@@ -1394,10 +1362,7 @@ class Gx2Fitter {
       }
 
       if (extendedSystem.chi2() > oldChi2sum + 1e-5) {
-        ACTS_DEBUG("chi2 not converging monotonically");
-
-        updateGx2fCovarianceParams(fullCovariancePredicted, extendedSystem);
-        break;
+        ACTS_DEBUG("chi2 not converging monotonically in update " << nUpdate);
       }
 
       // If this is the final iteration, update the covariance and break.
@@ -1418,26 +1383,187 @@ class Gx2Fitter {
         break;
       }
 
-      if (multipleScattering) {
-        // update the scattering angles
-        for (std::size_t matSurface = 0; matSurface < nMaterialSurfaces;
-             matSurface++) {
-          const std::size_t deltaPosition = eBoundSize + 2 * matSurface;
-          const GeometryIdentifier geoId = geoIdVector[matSurface];
-          const auto scatteringMapId = scatteringMap.find(geoId);
-          assert(scatteringMapId != scatteringMap.end() &&
-                 "No scattering angles found for material surface.");
-          scatteringMapId->second.scatteringAngles().block<2, 1>(2, 0) +=
-              deltaParamsExtended.block<2, 1>(deltaPosition, 0).eval();
-        }
-      }
-
       oldChi2sum = extendedSystem.chi2();
     }
     ACTS_DEBUG("Finished to iterate");
     ACTS_VERBOSE("Final parameters: " << params.parameters().transpose());
     /// Finish Fitting /////////////////////////////////////////////////////////
 
+    /// Actual MATERIAL Fitting ////////////////////////////////////////////////
+    ACTS_DEBUG("Start to evaluate material");
+    if (multipleScattering) {
+      // set up propagator and co
+      Acts::GeometryContext geoCtx = gx2fOptions.geoContext;
+      Acts::MagneticFieldContext magCtx = gx2fOptions.magFieldContext;
+      // Set options for propagator
+      PropagatorOptions propagatorOptions(geoCtx, magCtx);
+
+      // Add the measurement surface as external surface to the navigator.
+      // We will try to hit those surface by ignoring boundary checks.
+      for (const auto& [surfaceId, _] : inputMeasurements) {
+        propagatorOptions.navigation.insertExternalSurface(surfaceId);
+      }
+
+      auto& gx2fActor = propagatorOptions.actorList.template get<GX2FActor>();
+      gx2fActor.inputMeasurements = &inputMeasurements;
+      gx2fActor.multipleScattering = true;
+      gx2fActor.extensions = gx2fOptions.extensions;
+      gx2fActor.calibrationContext = &gx2fOptions.calibrationContext.get();
+      gx2fActor.actorLogger = m_actorLogger.get();
+      gx2fActor.scatteringMap = &scatteringMap;
+      gx2fActor.parametersWithHypothesis = &params;
+
+      auto propagatorState = m_propagator.makeState(params, propagatorOptions);
+
+      auto& r = propagatorState.template get<Gx2FitterResult<traj_t>>();
+      r.fittedStates = &trajectoryTempBackend;
+
+      // Clear the track container. It could be more performant to update the
+      // existing states, but this needs some more thinking.
+      trackContainerTemp.clear();
+
+      auto propagationResult = m_propagator.template propagate(propagatorState);
+
+      // Run the fitter
+      auto result = m_propagator.template makeResult(std::move(propagatorState),
+                                                     propagationResult,
+                                                     propagatorOptions, false);
+
+      if (!result.ok()) {
+        ACTS_ERROR("Propagation failed: " << result.error());
+        return result.error();
+      }
+
+      // TODO Improve Propagator + Actor [allocate before loop], rewrite
+      // makeMeasurements
+      auto& propRes = *result;
+      GX2FResult gx2fResult = std::move(propRes.template get<GX2FResult>());
+
+      if (!gx2fResult.result.ok()) {
+        ACTS_INFO("GlobalChiSquareFitter failed in actor: "
+                  << gx2fResult.result.error() << ", "
+                  << gx2fResult.result.error().message());
+        return gx2fResult.result.error();
+      }
+
+      auto track = trackContainerTemp.makeTrack();
+      tipIndex = gx2fResult.lastMeasurementIndex;
+
+      // It could happen, that no measurements were found. Then the track would
+      // be empty and the following operations would be invalid. Usually, this
+      // only happens during the first iteration, due to bad initial parameters.
+      if (tipIndex == Acts::MultiTrajectoryTraits::kInvalid) {
+        ACTS_INFO("Did not find any measurements in material fit.");
+        return Experimental::GlobalChiSquareFitterError::NotEnoughMeasurements;
+      }
+
+      track.tipIndex() = tipIndex;
+      track.linkForward();
+
+      // Count the material surfaces, to set up the system. In the multiple
+      // scattering case, we need to extend our system.
+      std::size_t nMaterialSurfaces = 0;
+      ACTS_DEBUG("Count the valid material surfaces.");
+      for (const auto& trackState : track.trackStates()) {
+        const auto typeFlags = trackState.typeFlags();
+        const bool stateHasMaterial =
+            typeFlags.test(TrackStateFlag::MaterialFlag);
+
+        if (!stateHasMaterial) {
+          continue;
+        }
+
+        // Get and store geoId for the current material surface
+        const GeometryIdentifier geoId =
+            trackState.referenceSurface().geometryId();
+
+        const auto scatteringMapId = scatteringMap.find(geoId);
+        assert(scatteringMapId != scatteringMap.end() &&
+               "No scattering angles found for material surface.");
+        if (!scatteringMapId->second.materialIsValid()) {
+          continue;
+        }
+
+        nMaterialSurfaces++;
+      }
+
+      // We need 6 dimensions for the bound parameters and 2 * nMaterialSurfaces
+      // dimensions for the scattering angles.
+      const std::size_t dimsExtendedParams = eBoundSize + 2 * nMaterialSurfaces;
+
+      // System that we fill with the information gathered by the actor and
+      // evaluate later
+      Gx2fSystem extendedSystem{dimsExtendedParams};
+
+      // This vector stores the IDs for each visited material in order. We use
+      // it later for updating the scattering angles. We cannot use
+      // scatteringMap directly, since we cannot guarantee, that we will visit
+      // all stored material in each propagation.
+      std::vector<GeometryIdentifier> geoIdVector;
+
+      fillGx2fSystem(track, extendedSystem, true, scatteringMap, geoIdVector,
+                     *m_addToSumLogger);
+
+      chi2sum = extendedSystem.chi2();
+
+      // This check takes into account the evaluated dimensions of the
+      // measurements. To fit, we need at least NDF+1 measurements. However, we
+      // count n-dimensional measurements for n measurements, reducing the
+      // effective number of needed measurements. We might encounter the case,
+      // where we cannot use some (parts of a) measurements, maybe if we do not
+      // support that kind of measurement. This is also taken into account here.
+      // We skip the check during the first iteration, since we cannot guarantee
+      // to hit all/enough measurement surfaces with the initial parameter
+      // guess.
+      if ((nUpdate > 0) && !extendedSystem.isWellDefined()) {
+        ACTS_INFO("Not enough measurements. Require "
+                  << extendedSystem.findRequiredNdf() + 1 << ", but only "
+                  << extendedSystem.ndf() << " could be used.");
+        return Experimental::GlobalChiSquareFitterError::NotEnoughMeasurements;
+      }
+
+      // calculate delta params [a] * delta = b
+      Eigen::VectorXd deltaParamsExtended =
+          extendedSystem.aMatrix().colPivHouseholderQr().solve(
+              extendedSystem.bVector());
+
+      deltaParams = deltaParamsExtended.topLeftCorner<eBoundSize, 1>().eval();
+
+      ACTS_VERBOSE("aMatrix:\n"
+                   << extendedSystem.aMatrix() << "\n"
+                   << "bVector:\n"
+                   << extendedSystem.bVector() << "\n"
+                   << "deltaParams:\n"
+                   << deltaParams << "\n"
+                   << "deltaParamsExtended:\n"
+                   << deltaParamsExtended << "\n"
+                   << "oldChi2sum = " << oldChi2sum << "\n"
+                   << "chi2sum = " << extendedSystem.chi2());
+
+      chi2sum = extendedSystem.chi2();
+
+      // update params
+      params.parameters() += deltaParams;
+
+      // update the scattering angles
+      for (std::size_t matSurface = 0; matSurface < nMaterialSurfaces;
+           matSurface++) {
+        const std::size_t deltaPosition = eBoundSize + 2 * matSurface;
+        const GeometryIdentifier geoId = geoIdVector[matSurface];
+        const auto scatteringMapId = scatteringMap.find(geoId);
+        assert(scatteringMapId != scatteringMap.end() &&
+               "No scattering angles found for material surface.");
+        scatteringMapId->second.scatteringAngles().block<2, 1>(2, 0) +=
+            deltaParamsExtended.block<2, 1>(deltaPosition, 0).eval();
+      }
+
+      updateGx2fCovarianceParams(fullCovariancePredicted, extendedSystem);
+    }
+    ACTS_DEBUG("Finished to evaluate material");
+    ACTS_VERBOSE(
+        "Final parameters after material: " << params.parameters().transpose());
+    /// Finish MATERIAL Fitting ////////////////////////////////////////////////
+
     ACTS_VERBOSE("Final scattering angles:");
     for (const auto& [key, value] : scatteringMap) {
       if (!value.materialIsValid()) {

Tests/UnitTests/Core/TrackFitting/Gx2fTests.cpp

@@ -1102,14 +1102,14 @@ BOOST_AUTO_TEST_CASE(Material) {
   // Parameters
   // We need quite coarse checks here, since on different builds
   // the created measurements differ in the randomness
-  BOOST_CHECK_CLOSE(track.parameters()[eBoundLoc0], -11., 7e0);
-  BOOST_CHECK_CLOSE(track.parameters()[eBoundLoc1], -15., 6e0);
-  BOOST_CHECK_CLOSE(track.parameters()[eBoundPhi], 1e-5, 1e3);
-  BOOST_CHECK_CLOSE(track.parameters()[eBoundTheta], M_PI / 2, 1e-3);
+  BOOST_CHECK_CLOSE(track.parameters()[eBoundLoc0], -11., 26e0);
+  BOOST_CHECK_CLOSE(track.parameters()[eBoundLoc1], -15., 15e0);
+  BOOST_CHECK_CLOSE(track.parameters()[eBoundPhi], 1e-5, 1.1e3);
+  BOOST_CHECK_CLOSE(track.parameters()[eBoundTheta], M_PI / 2, 1e-1);
   BOOST_CHECK_EQUAL(track.parameters()[eBoundQOverP], 1);
   BOOST_CHECK_CLOSE(track.parameters()[eBoundTime],
                     startParametersFit.parameters()[eBoundTime], 1e-6);
-  //  BOOST_CHECK_CLOSE(track.covariance().determinant(), 1e-27, 4e0);
+  BOOST_CHECK_CLOSE(track.covariance().determinant(), 3.5e-27, 1e1);
 
   // Convergence
   BOOST_CHECK_EQUAL(