Skip to content

Commit

Permalink
feat: GX2F: abort condition relChi2changeCutOff (#2586)
Browse files Browse the repository at this point in the history 
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:
 - #2574
  • Loading branch information
@AJPfleger
AJPfleger authored Nov 10, 2023
1 parent cf77a53 commit c3ceb12
Show file tree
Hide file tree
Showing 2 changed files with 156 additions and 17 deletions.
52 changes: 39 additions & 13 deletions Core/include/Acts/TrackFitting/GlobalChiSquareFitter.hpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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),
Expand All @@ -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;
Expand Down Expand Up @@ -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 <typename traj_t>
Expand Down Expand Up @@ -623,6 +628,7 @@ class Gx2Fitter {
start_parameters_t params = sParameters;
BoundVector deltaParams = BoundVector::Zero();
double chi2sum = 0;
double oldChi2sum = std::numeric_limits<double>::max();
BoundMatrix aMatrix = BoundMatrix::Zero();
BoundVector bVector = BoundVector::Zero();

Expand All @@ -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);

Expand Down Expand Up @@ -719,18 +727,27 @@ class Gx2Fitter {
.solve(bVector.topLeftCorner<reducedMatrixSize, 1>());
}

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);
Expand Down Expand Up @@ -770,13 +787,22 @@ class Gx2Fitter {

ACTS_VERBOSE("final covariance:\n" << fullCovariancePredicted);

if (!trackContainer.hasColumn(Acts::hashString("Gx2fnUpdateColumn"))) {
trackContainer.template addColumn<size_t>("Gx2fnUpdateColumn");
}

// Prepare track for return
auto track = trackContainer.getTrack(trackContainer.addTrack());
track.tipIndex() = tipIndex;
track.parameters() = params.parameters();
track.covariance() = fullCovariancePredicted;
track.setReferenceSurface(params.referenceSurface().getSharedPtr());

if (trackContainer.hasColumn(Acts::hashString("Gx2fnUpdateColumn"))) {
ACTS_DEBUG("Add nUpdate to track")
track.template component<std::size_t>("Gx2fnUpdateColumn") = nUpdate;
}

// TODO write test for calculateTrackQuantities
calculateTrackQuantities(track);

Expand Down
121 changes: 117 additions & 4 deletions Tests/UnitTests/Core/TrackFitting/Gx2fTests.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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{}};
Expand All @@ -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<std::size_t, hashString("Gx2fnUpdateColumn")>()),
0);

ACTS_INFO("*** Test: NoFit -- Finish");
}
Expand Down Expand Up @@ -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{}};
Expand All @@ -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<std::size_t, hashString("Gx2fnUpdateColumn")>()),
5);

ACTS_INFO("*** Test: Fit5Iterations -- Finish");
}
Expand Down Expand Up @@ -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{}};
Expand All @@ -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<std::size_t, hashString("Gx2fnUpdateColumn")>()),
5);

ACTS_INFO("*** Test: MixedDetector -- Finish");
}
Expand Down Expand Up @@ -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{}};
Expand All @@ -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<std::size_t, hashString("Gx2fnUpdateColumn")>()),
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<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();

using RecoStepper = EigenStepper<>;
const auto recoPropagator =
makeConstantFieldPropagator<RecoStepper>(detector.geometry, 0_T);

using RecoPropagator = decltype(recoPropagator);
using Gx2Fitter =
Experimental::Gx2Fitter<RecoPropagator, VectorMultiTrajectory>;
Gx2Fitter fitter(recoPropagator, 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), 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<std::size_t, hashString("Gx2fnUpdateColumn")>()),
10);

ACTS_INFO("*** Test: relChi2changeCutOff -- Finish");
}
BOOST_AUTO_TEST_SUITE_END()
} // namespace Test
} // namespace Acts

0 comments on commit c3ceb12

Please sign in to comment.