From c3ceb12f9cd17f9da3f27198230b87a826ea99ba Mon Sep 17 00:00:00 2001 From: "Alexander J. Pfleger" <70842573+AJPfleger@users.noreply.github.com> Date: Fri, 10 Nov 2023 17:41:22 +0100 Subject: [PATCH] feat: GX2F: abort condition `relChi2changeCutOff` (#2586) A new abort condition for the Global Chi Square Fitter: `relChi2changeCutOff`. Before we could only steer with `nUpdateMax` when to stop to iterate. This relative cut-off checks if the `chi2sum` already converged. TODO: detect in the unitTest, that the system stopped before reaching `nUpdateMax`. blocked by: - https://github.com/acts-project/acts/pull/2574 --- .../TrackFitting/GlobalChiSquareFitter.hpp | 52 ++++++-- .../UnitTests/Core/TrackFitting/Gx2fTests.cpp | 121 +++++++++++++++++- 2 files changed, 156 insertions(+), 17 deletions(-) diff --git a/Core/include/Acts/TrackFitting/GlobalChiSquareFitter.hpp b/Core/include/Acts/TrackFitting/GlobalChiSquareFitter.hpp index f079418493f..22f71daa4c1 100644 --- a/Core/include/Acts/TrackFitting/GlobalChiSquareFitter.hpp +++ b/Core/include/Acts/TrackFitting/GlobalChiSquareFitter.hpp @@ -114,7 +114,8 @@ struct Gx2FitterOptions { const FreeToBoundCorrection& freeToBoundCorrection_ = FreeToBoundCorrection(false), const std::size_t nUpdateMax_ = 5, - const bool zeroField_ = false) + const bool zeroField_ = false, + double relChi2changeCutOff_ = 1e-5) : geoContext(gctx), magFieldContext(mctx), calibrationContext(cctx), @@ -125,7 +126,8 @@ struct Gx2FitterOptions { energyLoss(eLoss), freeToBoundCorrection(freeToBoundCorrection_), nUpdateMax(nUpdateMax_), - zeroField(zeroField_) {} + zeroField(zeroField_), + relChi2changeCutOff(relChi2changeCutOff_) {} /// Contexts are required and the options must not be default-constructible. Gx2FitterOptions() = delete; @@ -155,11 +157,14 @@ struct Gx2FitterOptions { /// transformation FreeToBoundCorrection freeToBoundCorrection; - /// Max number of iterations during the fit + /// Max number of iterations during the fit (abort condition) std::size_t nUpdateMax = 5; /// Disables the QoP fit in case of missing B-field bool zeroField = false; + + /// Check for convergence (abort condition). Set to 0 to skip. + double relChi2changeCutOff = 1e-7; }; template @@ -623,6 +628,7 @@ class Gx2Fitter { start_parameters_t params = sParameters; BoundVector deltaParams = BoundVector::Zero(); double chi2sum = 0; + double oldChi2sum = std::numeric_limits::max(); BoundMatrix aMatrix = BoundMatrix::Zero(); BoundVector bVector = BoundVector::Zero(); @@ -638,7 +644,9 @@ class Gx2Fitter { // Iterate the fit and improve result. Abort after n steps or after // convergence - for (std::size_t nUpdate = 0; nUpdate < gx2fOptions.nUpdateMax; nUpdate++) { + // nUpdate is initialized outside to save its state for the track + size_t nUpdate = 0; + for (nUpdate = 0; nUpdate < gx2fOptions.nUpdateMax; nUpdate++) { ACTS_VERBOSE("nUpdate = " << nUpdate + 1 << "/" << gx2fOptions.nUpdateMax); @@ -719,18 +727,27 @@ class Gx2Fitter { .solve(bVector.topLeftCorner()); } - ACTS_VERBOSE("chi2sum = " << chi2sum); - ACTS_VERBOSE("aMatrix:\n" << aMatrix); - ACTS_VERBOSE("bVector:\n" << bVector); - ACTS_VERBOSE("deltaParams:\n" << deltaParams); + ACTS_VERBOSE("aMatrix:\n" + << aMatrix << "\n" + << "bVector:\n" + << bVector << "\n" + << "deltaParams:\n" + << deltaParams << "\n" + << "oldChi2sum = " << oldChi2sum << "\n" + << "chi2sum = " << chi2sum); tipIndex = gx2fResult.lastMeasurementIndex; - // TODO check delta params and abort - // similar to: - // if (sum(delta_params) < 1e-3) { - // break; - // } + if ((gx2fOptions.relChi2changeCutOff != 0) && (nUpdate > 0) && + (std::abs(chi2sum / oldChi2sum - 1) < + gx2fOptions.relChi2changeCutOff)) { + ACTS_VERBOSE("Abort with relChi2changeCutOff after " + << nUpdate + 1 << "/" << gx2fOptions.nUpdateMax + << " iterations."); + break; + } + + oldChi2sum = chi2sum; } ACTS_DEBUG("Finished to iterate"); ACTS_VERBOSE("final params:\n" << params); @@ -770,6 +787,10 @@ class Gx2Fitter { ACTS_VERBOSE("final covariance:\n" << fullCovariancePredicted); + if (!trackContainer.hasColumn(Acts::hashString("Gx2fnUpdateColumn"))) { + trackContainer.template addColumn("Gx2fnUpdateColumn"); + } + // Prepare track for return auto track = trackContainer.getTrack(trackContainer.addTrack()); track.tipIndex() = tipIndex; @@ -777,6 +798,11 @@ class Gx2Fitter { track.covariance() = fullCovariancePredicted; track.setReferenceSurface(params.referenceSurface().getSharedPtr()); + if (trackContainer.hasColumn(Acts::hashString("Gx2fnUpdateColumn"))) { + ACTS_DEBUG("Add nUpdate to track") + track.template component("Gx2fnUpdateColumn") = nUpdate; + } + // TODO write test for calculateTrackQuantities calculateTrackQuantities(track); diff --git a/Tests/UnitTests/Core/TrackFitting/Gx2fTests.cpp b/Tests/UnitTests/Core/TrackFitting/Gx2fTests.cpp index d0181c160ed..c9b9d5a05dd 100644 --- a/Tests/UnitTests/Core/TrackFitting/Gx2fTests.cpp +++ b/Tests/UnitTests/Core/TrackFitting/Gx2fTests.cpp @@ -214,7 +214,7 @@ BOOST_AUTO_TEST_CASE(NoFit) { Experimental::Gx2FitterOptions gx2fOptions( geoCtx, magCtx, calCtx, extensions, PropagatorPlainOptions(), rSurface, - false, false, FreeToBoundCorrection(false), 0, true); + false, false, FreeToBoundCorrection(false), 0, true, 0); Acts::TrackContainer tracks{Acts::VectorTrackContainer{}, Acts::VectorMultiTrajectory{}}; @@ -232,6 +232,10 @@ BOOST_AUTO_TEST_CASE(NoFit) { BOOST_CHECK_EQUAL(track.nHoles(), 0u); BOOST_CHECK_EQUAL(track.parameters(), startParametersFit.parameters()); BOOST_CHECK_EQUAL(track.covariance(), BoundMatrix::Identity()); + BOOST_CHECK_EQUAL( + (track + .template component()), + 0); ACTS_INFO("*** Test: NoFit -- Finish"); } @@ -301,7 +305,7 @@ BOOST_AUTO_TEST_CASE(Fit5Iterations) { const Experimental::Gx2FitterOptions gx2fOptions( tgContext, mfContext, calContext, extensions, PropagatorPlainOptions(), - rSurface, false, false, FreeToBoundCorrection(false), 5, true); + rSurface, false, false, FreeToBoundCorrection(false), 5, true, 0); Acts::TrackContainer tracks{Acts::VectorTrackContainer{}, Acts::VectorMultiTrajectory{}}; @@ -326,6 +330,10 @@ BOOST_AUTO_TEST_CASE(Fit5Iterations) { BOOST_CHECK_EQUAL(track.parameters()[eBoundQOverP], 1); BOOST_CHECK_CLOSE(track.parameters()[eBoundTime], 12591.2832360000, 1e-6); BOOST_CHECK_CLOSE(track.covariance().determinant(), 1e-27, 4e0); + BOOST_CHECK_EQUAL( + (track + .template component()), + 5); ACTS_INFO("*** Test: Fit5Iterations -- Finish"); } @@ -404,7 +412,7 @@ BOOST_AUTO_TEST_CASE(MixedDetector) { const Experimental::Gx2FitterOptions gx2fOptions( tgContext, mfContext, calContext, extensions, PropagatorPlainOptions(), - rSurface, false, false, FreeToBoundCorrection(false), 5, true); + rSurface, false, false, FreeToBoundCorrection(false), 5, true, 0); Acts::TrackContainer tracks{Acts::VectorTrackContainer{}, Acts::VectorMultiTrajectory{}}; @@ -429,6 +437,10 @@ BOOST_AUTO_TEST_CASE(MixedDetector) { 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); + BOOST_CHECK_EQUAL( + (track + .template component()), + 5); ACTS_INFO("*** Test: MixedDetector -- Finish"); } @@ -497,7 +509,7 @@ BOOST_AUTO_TEST_CASE(FitWithBfield) { const Experimental::Gx2FitterOptions gx2fOptions( tgContext, mfContext, calContext, extensions, PropagatorPlainOptions(), - rSurface, false, false, FreeToBoundCorrection(false), 5, false); + rSurface, false, false, FreeToBoundCorrection(false), 5, false, 0); Acts::TrackContainer tracks{Acts::VectorTrackContainer{}, Acts::VectorMultiTrajectory{}}; @@ -523,10 +535,111 @@ BOOST_AUTO_TEST_CASE(FitWithBfield) { BOOST_CHECK_CLOSE(track.parameters()[eBoundQOverP], 0.5, 2e-1); BOOST_CHECK_CLOSE(track.parameters()[eBoundTime], 12591.2832360000, 1e-6); BOOST_CHECK_CLOSE(track.covariance().determinant(), 8e-35, 4e0); + BOOST_CHECK_EQUAL( + (track + .template component()), + 5); ACTS_INFO("*** Test: FitWithBfield -- Finish"); } +BOOST_AUTO_TEST_CASE(relChi2changeCutOff) { + ACTS_LOCAL_LOGGER(Acts::getDefaultLogger("Gx2fTests", logLevel)); + ACTS_INFO("*** Test: relChi2changeCutOff -- Start"); + + // Create a test context + GeometryContext tgContext = GeometryContext(); + + Detector detector; + const size_t nSurfaces = 5; + detector.geometry = makeToyDetector(tgContext, nSurfaces); + + ACTS_DEBUG("Go to propagator"); + + auto parametersMeasurements = makeParameters(); + auto startParametersFit = makeParameters(7_mm, 11_mm, 15_mm, 42_ns, 10_degree, + 80_degree, 1_GeV, 1_e); + + // 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}}; + MeasurementResolutionMap resolutions = { + {Acts::GeometryIdentifier().setVolume(0), resPixel}}; + + // simulation propagator + using SimPropagator = + Acts::Propagator; + 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 = ¶metersMeasurements.referenceSurface(); + + using RecoStepper = EigenStepper<>; + const auto recoPropagator = + makeConstantFieldPropagator(detector.geometry, 0_T); + + using RecoPropagator = decltype(recoPropagator); + using Gx2Fitter = + Experimental::Gx2Fitter; + Gx2Fitter fitter(recoPropagator, gx2fLogger->clone()); + + Experimental::Gx2FitterExtensions extensions; + extensions.calibrator + .connect<&testSourceLinkCalibrator>(); + 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), 500, true, 1e-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(), nSurfaces - 1); + BOOST_CHECK(track.hasReferenceSurface()); + BOOST_CHECK_EQUAL(track.nMeasurements(), nSurfaces); + 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], -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_EQUAL(track.parameters()[eBoundQOverP], 1); + BOOST_CHECK_CLOSE(track.parameters()[eBoundTime], 12591.2832360000, 1e-6); + BOOST_CHECK_CLOSE(track.covariance().determinant(), 1e-27, 4e0); + BOOST_CHECK_LT( + (track + .template component()), + 10); + + ACTS_INFO("*** Test: relChi2changeCutOff -- Finish"); +} BOOST_AUTO_TEST_SUITE_END() } // namespace Test } // namespace Acts