feat: Support also 1d-measurements in GX2F

Core/include/Acts/TrackFitting/GlobalChiSquareFitter.hpp

@@ -202,16 +202,70 @@ struct Gx2FitterResult {
   Result<void> result{Result<void>::success()};
 
   // collectors
-  std::vector<ActsVector<2>> collectorResiduals;
-  std::vector<ActsSquareMatrix<2>> collectorCovariance;
-  std::vector<BoundMatrix> collectorJacobians;
+  std::vector<ActsScalar> collectorResiduals;
+  std::vector<ActsScalar> collectorCovariances;
+  std::vector<BoundVector> collectorProjectedJacobians;
 
   BoundMatrix jacobianFromStart = BoundMatrix::Identity();
 
   // Count how many surfaces have been hit
   size_t surfaceCount = 0;
 };
 
+/// Collector for the GX2F Actor
+/// The collector prepares each measurement for the actual fitting process. Each
+/// n-dimensional measurement is split into n 1-dimensional linearly independent
+/// measurements. Then the collector saves the following information:
+/// - Residual: Calculated from measurement and prediction
+/// - Covariance: The covariance of the measurement
+/// - Projected Jacobian: This implicitly contains the measurement type
+/// It also checks if the covariance is above a threshold, to detect and avoid
+/// too small covariances for a stable fit.
+///
+/// @tparam measDim Number of dimensions of the measurement
+/// @tparam traj_t The trajectory type
+///
+/// @param trackStateProxy is the current track state
+/// @param result is the mutable result/cache object
+/// @param logger a logger instance
+template <size_t measDim, typename traj_t>
+void collector(typename traj_t::TrackStateProxy& trackStateProxy,
+               Gx2FitterResult<traj_t>& result, const Logger& logger) {
+  auto predicted = trackStateProxy.predicted();
+  auto measurement = trackStateProxy.template calibrated<measDim>();
+  auto covarianceMeasurement =
+      trackStateProxy.template calibratedCovariance<measDim>();
+  // calculate residuals and return with covariances and jacobians
+  auto projJacobian =
+      (trackStateProxy.effectiveProjector() * result.jacobianFromStart).eval();
+  auto projPredicted =
+      (trackStateProxy.effectiveProjector() * predicted).eval();
+
+  ACTS_VERBOSE("Processing and collecting measurements in Actor:\n"
+               << "\tMeasurement:\t" << measurement.transpose()
+               << "\n\tPredicted:\t" << predicted.transpose()
+               << "\n\tProjector:\t" << trackStateProxy.effectiveProjector()
+               << "\n\tProjected Jacobian:\t" << projJacobian
+               << "\n\tCovariance Measurements:\t" << covarianceMeasurement);
+
+  for (size_t i = 0; i < measDim; i++) {
+    if (covarianceMeasurement(i, i) < 1e-10) {
+      ACTS_WARNING("Invalid covariance of measurement: cov(" << i << "," << i
+                                                             << ") ~ 0")
+      continue;
+    }
+
+    result.collectorResiduals.push_back(measurement[i] - projPredicted[i]);
+    result.collectorCovariances.push_back(covarianceMeasurement(i, i));
+    result.collectorProjectedJacobians.push_back(projJacobian.row(i));
+
+    ACTS_VERBOSE("\tSplitting the measurement:\n"
+                 << "\t\tResidual:\t" << measurement[i] - projPredicted[i]
+                 << "\n\t\tCovariance:\t" << covarianceMeasurement(i, i)
+                 << "\n\t\tProjected Jacobian:\t" << projJacobian.row(i));
+  }
+}
+
 /// Global Chi Square fitter (GX2F) implementation.
 ///
 /// @tparam propagator_t Type of the propagation class
@@ -363,41 +417,32 @@ class Gx2Fitter {
 
           // Fill the track state
           trackStateProxy.predicted() = std::move(boundParams.parameters());
-          auto predicted = trackStateProxy.predicted();
 
           // We have predicted parameters, so calibrate the uncalibrated input
           // measurement
           extensions.calibrator(state.geoContext, *calibrationContext,
                                 sourcelink_it->second, trackStateProxy);
 
-          const size_t measdimPlaceholder = 2;
-          auto measurement =
-              trackStateProxy.template calibrated<measdimPlaceholder>();
-          auto covarianceMeasurement =
-              trackStateProxy
-                  .template calibratedCovariance<measdimPlaceholder>();
-          // calculate residuals and return with covariances and jacobians
-          ActsVector<2> residual;
-          for (long i = 0; i < measurement.size(); i++) {
-            residual[i] = measurement[i] - predicted[i];
+          if (trackStateProxy.calibratedSize() == 1) {
+            collector<1>(trackStateProxy, result, *actorLogger);
+          } else if (trackStateProxy.calibratedSize() == 2) {
+            collector<2>(trackStateProxy, result, *actorLogger);
+          } else {
+            ACTS_WARNING(
+                "Only measurements of 1 and 2 dimensions are implemented yet.");
           }
-          ACTS_VERBOSE("Measurement in Actor:\n" << measurement);
-          result.collectorResiduals.push_back(residual);
-          result.collectorCovariance.push_back(covarianceMeasurement);
 
           if (boundParams.covariance().has_value()) {
             trackStateProxy.predictedCovariance() =
                 std::move(*boundParams.covariance());
           }
 
-          result.collectorJacobians.push_back(result.jacobianFromStart);
-
           trackStateProxy.jacobian() = std::move(jacobian);
           trackStateProxy.pathLength() = std::move(pathLength);
         }
       }
 
-      if (result.surfaceCount > 11) {
+      if (result.surfaceCount > 17) {
         ACTS_INFO("Actor: finish due to limit. Result might be garbage.");
         result.finished = true;
       }
@@ -529,10 +574,10 @@ class Gx2Fitter {
 
       ACTS_VERBOSE("gx2fResult.collectorResiduals.size() = "
                    << gx2fResult.collectorResiduals.size());
-      ACTS_VERBOSE("gx2fResult.collectorCovariance.size() = "
-                   << gx2fResult.collectorCovariance.size());
-      ACTS_VERBOSE("gx2fResult.collectorJacobians.size() = "
-                   << gx2fResult.collectorJacobians.size());
+      ACTS_VERBOSE("gx2fResult.collectorCovariances.size() = "
+                   << gx2fResult.collectorCovariances.size());
+      ACTS_VERBOSE("gx2fResult.collectorProjectedJacobians.size() = "
+                   << gx2fResult.collectorProjectedJacobians.size());
 
       chi2sum = 0;
       aMatrix = BoundMatrix::Zero();
@@ -541,27 +586,15 @@ class Gx2Fitter {
       // TODO generalize for non-2D measurements
       for (size_t iMeas = 0; iMeas < gx2fResult.collectorResiduals.size();
            iMeas++) {
-        ActsMatrix<2, eBoundSize> proj;
-
-        for (size_t i_ = 0; i_ < 2; i_++) {
-          for (size_t j_ = 0; j_ < eBoundSize; j_++) {
-            proj(i_, j_) = 0;
-          }
-        }
-        proj(0, 0) = 1;
-        proj(1, 1) = 1;
-
         const auto ri = gx2fResult.collectorResiduals[iMeas];
-        const auto covi = gx2fResult.collectorCovariance[iMeas];
-        const auto coviInv = covi.inverse();
+        const auto covi = gx2fResult.collectorCovariances[iMeas];
         const auto projectedJacobian =
-            proj * gx2fResult.collectorJacobians[iMeas];
+            gx2fResult.collectorProjectedJacobians[iMeas];
 
-        const double chi2meas = (ri.transpose() * coviInv * ri).eval()(0);
+        const double chi2meas = ri / covi * ri;
         const BoundMatrix aMatrixMeas =
-            projectedJacobian.transpose() * coviInv * projectedJacobian;
-        const BoundVector bVectorMeas =
-            projectedJacobian.transpose() * coviInv * ri;
+            projectedJacobian * projectedJacobian.transpose() / covi;
+        const BoundVector bVectorMeas = projectedJacobian / covi * ri;
 
         chi2sum += chi2meas;
         aMatrix += aMatrixMeas;
@@ -591,6 +624,7 @@ class Gx2Fitter {
       // }
     }
     ACTS_DEBUG("Finished to iterate");
+    ACTS_VERBOSE("final params:\n" << params);
     /// Finish Fitting /////////////////////////////////////////////////////////
 
     // Calculate covariance of the fitted parameters with inverse of [a]

Tests/UnitTests/Core/TrackFitting/Gx2fTests.cpp

@@ -144,8 +144,8 @@ std::shared_ptr<const TrackingGeometry> makeToyDetector(
 
   // Inner Volume - Build volume configuration
   CuboidVolumeBuilder::VolumeConfig volumeConfig;
-  volumeConfig.position = {2.5_m, 0., 0.};
-  volumeConfig.length = {5_m, 1_m, 1_m};
+  volumeConfig.position = {nSurfaces / 2. * 1_m, 0., 0.};
+  volumeConfig.length = {nSurfaces * 1_m, 1_m, 1_m};
   volumeConfig.layerCfg = layerConfig;
   volumeConfig.name = "Test volume";
   volumeConfig.volumeMaterial =
@@ -158,8 +158,8 @@ std::shared_ptr<const TrackingGeometry> makeToyDetector(
 
   // Outer volume - Build TrackingGeometry configuration
   CuboidVolumeBuilder::Config config;
-  config.position = {2.5_m, 0., 0.};
-  config.length = {5_m, 1_m, 1_m};
+  config.position = {nSurfaces / 2. * 1_m, 0., 0.};
+  config.length = {nSurfaces * 1_m, 1_m, 1_m};
   config.volumeCfg = {volumeConfig};
 
   cvb.setConfig(config);
@@ -357,6 +357,116 @@ BOOST_AUTO_TEST_CASE(Fit5Iterations) {
   ACTS_INFO("*** Test: Fit5Iterations -- Finish");
 }
 
+BOOST_AUTO_TEST_CASE(MixedDetector) {
+  ACTS_LOCAL_LOGGER(Acts::getDefaultLogger("Gx2fTests", logLevel));
+  ACTS_INFO("*** Test: MixedDetector -- Start");
+
+  // Create a test context
+  GeometryContext tgContext = GeometryContext();
+
+  Detector detector;
+  const size_t nSurfaces = 7;
+  detector.geometry = makeToyDetector(tgContext, nSurfaces);
+
+  ACTS_DEBUG("Go to propagator");
+
+  auto parametersMeasurements = makeParameters();
+  auto startParametersFit = makeParameters(10_mm, 10_mm, 10_mm, 42_ns,
+                                           10_degree, 80_degree, 1_GeV, 1_e);
+  //  auto startParametersFit = parametersMeasurements;
+  // Context objects
+  Acts::GeometryContext geoCtx;
+  Acts::MagneticFieldContext magCtx;
+  // Acts::CalibrationContext calCtx;
+  std::default_random_engine rng(42);
+
+  MeasurementResolution resPixel = {MeasurementType::eLoc01, {25_um, 50_um}};
+  MeasurementResolution resStrip0 = {MeasurementType::eLoc0, {25_um}};
+  MeasurementResolution resStrip1 = {MeasurementType::eLoc1, {50_um}};
+  MeasurementResolutionMap resolutions = {
+      {Acts::GeometryIdentifier().setVolume(2).setLayer(2), resPixel},
+      {Acts::GeometryIdentifier().setVolume(2).setLayer(4), resStrip0},
+      {Acts::GeometryIdentifier().setVolume(2).setLayer(6), resStrip1},
+      {Acts::GeometryIdentifier().setVolume(2).setLayer(8), resPixel},
+      {Acts::GeometryIdentifier().setVolume(2).setLayer(10), resStrip0},
+      {Acts::GeometryIdentifier().setVolume(2).setLayer(12), resStrip1},
+      {Acts::GeometryIdentifier().setVolume(2).setLayer(14), resPixel},
+  };
+
+  // simulation propagator
+  using SimPropagator =
+      Acts::Propagator<Acts::StraightLineStepper, Acts::Navigator>;
+  SimPropagator simPropagator = makeStraightPropagator(detector.geometry);
+  auto measurements = createMeasurements(
+      simPropagator, geoCtx, magCtx, parametersMeasurements, resolutions, rng);
+  auto sourceLinks = prepareSourceLinks(measurements.sourceLinks);
+  ACTS_VERBOSE("sourceLinks.size() = " << sourceLinks.size());
+
+  BOOST_REQUIRE_EQUAL(sourceLinks.size(), nSurfaces);
+
+  ACTS_DEBUG("Start fitting");
+  ACTS_VERBOSE("startParameter unsmeared:\n" << parametersMeasurements);
+  ACTS_VERBOSE("startParameter fit:\n" << startParametersFit);
+
+  const Surface* rSurface = &parametersMeasurements.referenceSurface();
+
+  Navigator::Config cfg{detector.geometry};
+  cfg.resolvePassive = false;
+  cfg.resolveMaterial = true;
+  cfg.resolveSensitive = true;
+  Navigator rNavigator(cfg);
+  // Configure propagation with deactivated B-field
+  auto bField = std::make_shared<ConstantBField>(Vector3(0., 0., 0.));
+  using RecoStepper = EigenStepper<>;
+  RecoStepper rStepper(bField);
+  using RecoPropagator = Propagator<RecoStepper, Navigator>;
+  RecoPropagator rPropagator(rStepper, rNavigator);
+
+  using Gx2Fitter =
+      Experimental::Gx2Fitter<RecoPropagator, VectorMultiTrajectory>;
+  Gx2Fitter Fitter(rPropagator, gx2fLogger->clone());
+
+  Experimental::Gx2FitterExtensions<VectorMultiTrajectory> extensions;
+  extensions.calibrator
+      .connect<&testSourceLinkCalibrator<VectorMultiTrajectory>>();
+  TestSourceLink::SurfaceAccessor surfaceAccessor{*detector.geometry};
+  extensions.surfaceAccessor
+      .connect<&TestSourceLink::SurfaceAccessor::operator()>(&surfaceAccessor);
+
+  MagneticFieldContext mfContext;
+  CalibrationContext calContext;
+
+  const Experimental::Gx2FitterOptions gx2fOptions(
+      tgContext, mfContext, calContext, extensions, PropagatorPlainOptions(),
+      rSurface, false, false, FreeToBoundCorrection(false), 5);
+
+  Acts::TrackContainer tracks{Acts::VectorTrackContainer{},
+                              Acts::VectorMultiTrajectory{}};
+
+  // Fit the track
+  auto res = Fitter.fit(sourceLinks.begin(), sourceLinks.end(),
+                        startParametersFit, gx2fOptions, tracks);
+
+  BOOST_REQUIRE(res.ok());
+
+  auto& track = *res;
+  BOOST_CHECK_EQUAL(track.tipIndex(), Acts::MultiTrajectoryTraits::kInvalid);
+  BOOST_CHECK(!track.hasReferenceSurface());
+  BOOST_CHECK_EQUAL(track.nMeasurements(), 0u);
+  BOOST_CHECK_EQUAL(track.nHoles(), 0u);
+  // We need quite coarse checks here, since on different builds
+  // the created measurements differ in the randomness
+  BOOST_CHECK_CLOSE(track.parameters()[eBoundLoc0], -10., 6e0);
+  BOOST_CHECK_CLOSE(track.parameters()[eBoundLoc1], -10., 6e0);
+  BOOST_CHECK_CLOSE(track.parameters()[eBoundPhi], 1e-5, 1e3);
+  BOOST_CHECK_CLOSE(track.parameters()[eBoundTheta], M_PI / 2, 1e-3);
+  BOOST_CHECK_EQUAL(track.parameters()[eBoundQOverP], 1);
+  BOOST_CHECK_CLOSE(track.parameters()[eBoundTime], 12591.2832360000, 1e-6);
+  BOOST_CHECK_CLOSE(track.covariance().determinant(), 2e-28, 1e0);
+
+  ACTS_INFO("*** Test: MixedDetector -- Finish");
+}
+
 BOOST_AUTO_TEST_SUITE_END()
 }  // namespace Test
 }  // namespace Acts