Skip to content

Commit

Permalink
Merge branch 'main' into acts-odd-light-part3
Browse files Browse the repository at this point in the history
  • Loading branch information
@asalzburger
asalzburger authored Sep 20, 2023
2 parents 0332a1b + 15eb6b3 commit c3a6257
Show file tree
Hide file tree
Showing 2 changed files with 189 additions and 45 deletions.
116 changes: 75 additions & 41 deletions Core/include/Acts/TrackFitting/GlobalChiSquareFitter.hpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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
Expand Down Expand Up @@ -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;
}
Expand Down Expand Up @@ -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();
Expand All @@ -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;
Expand Down Expand Up @@ -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]
Expand Down
118 changes: 114 additions & 4 deletions Tests/UnitTests/Core/TrackFitting/Gx2fTests.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -145,8 +145,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 =
Expand All @@ -159,8 +159,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);
Expand Down Expand Up @@ -358,6 +358,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

0 comments on commit c3a6257

Please sign in to comment.