Skip to content

Commit

Permalink
add chi2 cutoff
Browse files Browse the repository at this point in the history
  • Loading branch information
@AJPfleger
AJPfleger committed Oct 26, 2023
1 parent 4b9e46f commit 8eb88ed
Show file tree
Hide file tree
Showing 2 changed files with 126 additions and 16 deletions.
39 changes: 27 additions & 12 deletions Core/include/Acts/TrackFitting/GlobalChiSquareFitter.hpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -113,7 +113,8 @@ struct Gx2FitterOptions {
bool eLoss = false,
const FreeToBoundCorrection& freeToBoundCorrection_ =
FreeToBoundCorrection(false),
const size_t nUpdateMax_ = 5, const bool zeroField_ = false)
const size_t nUpdateMax_ = 5, const bool zeroField_ = false,
double relChi2changeCutOff_ = 1e-5)
: geoContext(gctx),
magFieldContext(mctx),
calibrationContext(cctx),
Expand All @@ -124,7 +125,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 @@ -154,11 +156,14 @@ struct Gx2FitterOptions {
/// transformation
FreeToBoundCorrection freeToBoundCorrection;

/// Max number of iterations during the fit
/// Max number of iterations during the fit (abort condition)
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 @@ -622,6 +627,7 @@ class Gx2Fitter {
start_parameters_t params = sParameters;
BoundVector deltaParams = BoundVector::Zero();
double chi2sum = 0;
double oldChi2sum = DBL_MAX;
BoundMatrix aMatrix = BoundMatrix::Zero();
BoundVector bVector = BoundVector::Zero();

Expand Down Expand Up @@ -726,18 +732,27 @@ class Gx2Fitter {
}
}

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
103 changes: 99 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 Down Expand Up @@ -301,7 +301,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 Down Expand Up @@ -404,7 +404,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 Down Expand Up @@ -497,7 +497,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 Down Expand Up @@ -527,6 +527,101 @@ BOOST_AUTO_TEST_CASE(FitWithBfield) {
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);
// TODO check for number of updates

ACTS_INFO("*** Test: relChi2changeCutOff -- Finish");
}

BOOST_AUTO_TEST_SUITE_END()
} // namespace Test
} // namespace Acts

0 comments on commit 8eb88ed

Please sign in to comment.