From f90906cab97cddec83fbc1b858b2135b361e4e73 Mon Sep 17 00:00:00 2001 From: "Alexander J. Pfleger" <70842573+AJPfleger@users.noreply.github.com> Date: Thu, 24 Aug 2023 18:37:59 +0700 Subject: [PATCH] feat: GX2F (experimental) + unit test (#2305) A first implementation of the new Global Chi Square Fitter (GX2F). It is probably not performant and lacks many features but it is stable and testable. That's why it should live in the namespace `Experimental`. --- .../TrackFitting/GlobalChiSquareFitter.hpp | 618 ++++++++++++++++++ .../Core/TrackFitting/CMakeLists.txt | 1 + .../UnitTests/Core/TrackFitting/Gx2fTests.cpp | 355 ++++++++++ 3 files changed, 974 insertions(+) create mode 100644 Core/include/Acts/TrackFitting/GlobalChiSquareFitter.hpp create mode 100644 Tests/UnitTests/Core/TrackFitting/Gx2fTests.cpp diff --git a/Core/include/Acts/TrackFitting/GlobalChiSquareFitter.hpp b/Core/include/Acts/TrackFitting/GlobalChiSquareFitter.hpp new file mode 100644 index 00000000000..d1522f99114 --- /dev/null +++ b/Core/include/Acts/TrackFitting/GlobalChiSquareFitter.hpp @@ -0,0 +1,618 @@ +// This file is part of the Acts project. +// +// Copyright (C) 2023 CERN for the benefit of the Acts project +// +// This Source Code Form is subject to the terms of the Mozilla Public +// License, v. 2.0. If a copy of the MPL was not distributed with this +// file, You can obtain one at http://mozilla.org/MPL/2.0/. + +#pragma once + +// Workaround for building on clang+libstdc++ +#include "Acts/Utilities/detail/ReferenceWrapperAnyCompat.hpp" + +#include "Acts/Definitions/Algebra.hpp" +#include "Acts/EventData/Measurement.hpp" +#include "Acts/EventData/MeasurementHelpers.hpp" +#include "Acts/EventData/MultiTrajectory.hpp" +#include "Acts/EventData/MultiTrajectoryHelpers.hpp" +#include "Acts/EventData/SourceLink.hpp" +#include "Acts/EventData/TrackHelpers.hpp" +#include "Acts/EventData/TrackParameters.hpp" +#include "Acts/EventData/VectorMultiTrajectory.hpp" +#include "Acts/Geometry/GeometryContext.hpp" +#include "Acts/MagneticField/MagneticFieldContext.hpp" +#include "Acts/Material/MaterialSlab.hpp" +#include "Acts/Propagator/AbortList.hpp" +#include "Acts/Propagator/ActionList.hpp" +#include "Acts/Propagator/ConstrainedStep.hpp" +#include "Acts/Propagator/DirectNavigator.hpp" +#include "Acts/Propagator/Navigator.hpp" +#include "Acts/Propagator/Propagator.hpp" +#include "Acts/Propagator/StandardAborters.hpp" +#include "Acts/Propagator/StraightLineStepper.hpp" +#include "Acts/Propagator/detail/PointwiseMaterialInteraction.hpp" +#include "Acts/TrackFitting/detail/KalmanUpdateHelpers.hpp" +#include "Acts/TrackFitting/detail/VoidKalmanComponents.hpp" +#include "Acts/Utilities/CalibrationContext.hpp" +#include "Acts/Utilities/Delegate.hpp" +#include "Acts/Utilities/Logger.hpp" +#include "Acts/Utilities/Result.hpp" + +#include +#include +#include + +namespace Acts { +namespace Experimental { + +/// Extension struct which holds delegates to customize the KF behavior +template +struct Gx2FitterExtensions { + using TrackStateProxy = typename MultiTrajectory::TrackStateProxy; + using ConstTrackStateProxy = + typename MultiTrajectory::ConstTrackStateProxy; + using Parameters = typename TrackStateProxy::Parameters; + + using Calibrator = Delegate; + + using Updater = Delegate(const GeometryContext&, TrackStateProxy, + Direction, const Logger&)>; + + using OutlierFinder = Delegate; + + /// The Calibrator is a dedicated calibration algorithm that allows + /// to calibrate measurements using track information, this could be + /// e.g. sagging for wires, module deformations, etc. + Calibrator calibrator; + + /// The updater incorporates measurement information into the track parameters + Updater updater; + + /// Determines whether a measurement is supposed to be considered as an + /// outlier + OutlierFinder outlierFinder; + + // TODO get an own Calibrator and Updater + /// Default constructor which connects the default void components + Gx2FitterExtensions() { + calibrator.template connect<&voidKalmanCalibrator>(); + updater.template connect<&voidKalmanUpdater>(); + outlierFinder.template connect<&voidOutlierFinder>(); + } +}; + +/// Combined options for the Global-Chi-Square fitter. +/// +/// @tparam traj_t The trajectory type +template +struct Gx2FitterOptions { + /// PropagatorOptions with context. + /// + /// @param gctx The geometry context for this fit + /// @param mctx The magnetic context for this fit + /// @param cctx The calibration context for this fit + /// @param extensions_ The KF extensions + /// @param pOptions The plain propagator options + /// @param rSurface The reference surface for the fit to be expressed at + /// @param mScattering Whether to include multiple scattering + /// @param eLoss Whether to include energy loss + /// @param freeToBoundCorrection_ Correction for non-linearity effect during transform from free to bound + /// @param nUpdateMax_ Max number of iterations for updating the parameters + Gx2FitterOptions(const GeometryContext& gctx, + const MagneticFieldContext& mctx, + std::reference_wrapper cctx, + Gx2FitterExtensions extensions_, + const PropagatorPlainOptions& pOptions, + const Surface* rSurface = nullptr, bool mScattering = false, + bool eLoss = false, + const FreeToBoundCorrection& freeToBoundCorrection_ = + FreeToBoundCorrection(false), + const size_t nUpdateMax_ = 5) + : geoContext(gctx), + magFieldContext(mctx), + calibrationContext(cctx), + extensions(extensions_), + propagatorPlainOptions(pOptions), + referenceSurface(rSurface), + multipleScattering(mScattering), + energyLoss(eLoss), + freeToBoundCorrection(freeToBoundCorrection_), + nUpdateMax(nUpdateMax_) {} + + /// Contexts are required and the options must not be default-constructible. + Gx2FitterOptions() = delete; + + /// Context object for the geometry + std::reference_wrapper geoContext; + /// Context object for the magnetic field + std::reference_wrapper magFieldContext; + /// context object for the calibration + std::reference_wrapper calibrationContext; + + Gx2FitterExtensions extensions; + + /// The trivial propagator options + PropagatorPlainOptions propagatorPlainOptions; + + /// The reference Surface + const Surface* referenceSurface = nullptr; + + /// Whether to consider multiple scattering + bool multipleScattering = false; + + /// Whether to consider energy loss + bool energyLoss = false; + + /// Whether to include non-linear correction during global to local + /// transformation + FreeToBoundCorrection freeToBoundCorrection; + + /// Max number of iterations during the fit + size_t nUpdateMax = 5; +}; + +template +struct Gx2FitterResult { + // Fitted states that the actor has handled. + traj_t* fittedStates{nullptr}; + + // This is the index of the 'tip' of the track stored in multitrajectory. + // This corresponds to the last measurement state in the multitrajectory. + // Since this KF only stores one trajectory, it is unambiguous. + // SIZE_MAX is the start of a trajectory. + size_t lastMeasurementIndex = Acts::MultiTrajectoryTraits::kInvalid; + + // This is the index of the 'tip' of the states stored in multitrajectory. + // This corresponds to the last state in the multitrajectory. + // Since this KF only stores one trajectory, it is unambiguous. + // SIZE_MAX is the start of a trajectory. + size_t lastTrackIndex = Acts::MultiTrajectoryTraits::kInvalid; + + // The optional Parameters at the provided surface + std::optional fittedParameters; + + // Counter for states with non-outlier measurements + size_t measurementStates = 0; + + // Counter for measurements holes + // A hole correspond to a surface with an associated detector element with no + // associated measurement. Holes are only taken into account if they are + // between the first and last measurements. + size_t measurementHoles = 0; + + // Counter for handled states + size_t processedStates = 0; + + // Indicator if track fitting has been done + bool finished = false; + + // Measurement surfaces without hits + std::vector missedActiveSurfaces; + + // Measurement surfaces handled in both forward and + // backward filtering + std::vector passedAgainSurfaces; + + Result result{Result::success()}; + + // collectors + std::vector> collectorResiduals; + std::vector> collectorCovariance; + std::vector collectorJacobians; + + BoundMatrix jacobianFromStart = BoundMatrix::Identity(); + + // Count how many surfaces have been hit + size_t surfaceCount = 0; +}; + +/// Global Chi Square fitter (GX2F) implementation. +/// +/// @tparam propagator_t Type of the propagation class +/// +/// TODO Write description +template +class Gx2Fitter { + /// The navigator type + using Gx2fNavigator = typename propagator_t::Navigator; + + /// The navigator has DirectNavigator type or not + static constexpr bool isDirectNavigator = + std::is_same::value; + + public: + Gx2Fitter(propagator_t pPropagator, + std::unique_ptr _logger = + getDefaultLogger("Gx2Fitter", Logging::INFO)) + : m_propagator(std::move(pPropagator)), + m_logger{std::move(_logger)}, + m_actorLogger{m_logger->cloneWithSuffix("Actor")} {} + + private: + /// The propagator for the transport and material update + propagator_t m_propagator; + + /// The logger instance + std::unique_ptr m_logger; + std::unique_ptr m_actorLogger; + + const Logger& logger() const { return *m_logger; } + + /// @brief Propagator Actor plugin for the GX2F + /// + /// @tparam parameters_t The type of parameters used for "local" parameters. + /// @tparam calibrator_t The type of calibrator + /// @tparam outlier_finder_t Type of the outlier finder class + /// + /// The GX2FnActor does not rely on the measurements to be + /// sorted along the track. /// TODO is this true? + template + class Actor { + public: + /// Broadcast the result_type + using result_type = Gx2FitterResult; + + /// The target surface + const Surface* targetSurface = nullptr; + + /// Allows retrieving measurements for a surface + const std::map* inputMeasurements = nullptr; + + /// Whether to consider multiple scattering. + bool multipleScattering = false; /// TODO implement later + + /// Whether to consider energy loss. + bool energyLoss = false; /// TODO implement later + + /// Whether to include non-linear correction during global to local + /// transformation + FreeToBoundCorrection freeToBoundCorrection; + + /// Input MultiTrajectory + std::shared_ptr> outputStates; + + /// The logger instance + const Logger* actorLogger{nullptr}; + + /// Logger helper + const Logger& logger() const { return *actorLogger; } + + Gx2FitterExtensions extensions; + + /// The Surface being + SurfaceReached targetReached; + + /// @brief Gx2f actor operation + /// + /// @tparam propagator_state_t is the type of Propagator state + /// @tparam stepper_t Type of the stepper + /// @tparam navigator_t Type of the navigator + /// + /// @param state is the mutable propagator state object + /// @param stepper The stepper in use + /// @param navigator The navigator in use + /// @param result is the mutable result state object + template + void operator()(propagator_state_t& state, const stepper_t& stepper, + const navigator_t& navigator, result_type& result, + const Logger& /*logger*/) const { + assert(result.fittedStates && "No MultiTrajectory set"); + + if (result.finished) { + return; + } + + // Add the measurement surface as external surface to navigator. + // We will try to hit those surface by ignoring boundary checks. + if constexpr (not isDirectNavigator) { + if (result.processedStates == 0) { + for (auto measurementIt = inputMeasurements->begin(); + measurementIt != inputMeasurements->end(); measurementIt++) { + navigator.insertExternalSurface(state.navigation, + measurementIt->first); + } + } + } + + // Update: + // - Waiting for a current surface + auto surface = navigator.currentSurface(state.navigation); + // std::string direction = state.stepping.navDir.toString(); + if (surface != nullptr) { + ++result.surfaceCount; + ACTS_VERBOSE("Measurement surface " << surface->geometryId() + << " detected."); + + // check if measurement surface + auto sourcelink_it = inputMeasurements->find(surface->geometryId()); + + if (sourcelink_it != inputMeasurements->end()) { + stepper.transportCovarianceToBound(state.stepping, *surface, + freeToBoundCorrection); + auto res = stepper.boundState(state.stepping, *surface, false, + freeToBoundCorrection); + if (!res.ok()) { + std::cout << "dbgActor: res = stepper.boundState res not ok" + << std::endl; + return; + } + auto& [boundParams, jacobian, pathLength] = *res; + result.jacobianFromStart = jacobian * result.jacobianFromStart; + + // add a full TrackState entry multi trajectory + // (this allocates storage for all components, we will set them later) + auto fittedStates = *result.fittedStates; + const auto newTrackIndex = fittedStates.addTrackState( + TrackStatePropMask::All, result.lastTrackIndex); + + // now get track state proxy back + auto trackStateProxy = fittedStates.getTrackState(newTrackIndex); + trackStateProxy.setReferenceSurface(surface->getSharedPtr()); + // assign the source link to the track state + trackStateProxy.setUncalibratedSourceLink(sourcelink_it->second); + + // 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, sourcelink_it->second, + trackStateProxy); + + const size_t measdimPlaceholder = 2; + auto measurement = + trackStateProxy.template calibrated(); + auto covarianceMeasurement = + trackStateProxy + .template calibratedCovariance(); + // 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]; + } + 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) { + ACTS_INFO("Actor: finish due to limit. Result might be garbage."); + result.finished = true; + } + } + }; + + /// Aborter can stay like this probably + template + class Aborter { + public: + /// Broadcast the result_type + using action_type = Actor; + + template + bool operator()(propagator_state_t& /*state*/, const stepper_t& /*stepper*/, + const navigator_t& /*navigator*/, const result_t& result, + const Logger& /*logger*/) const { + if (!result.result.ok() or result.finished) { + return true; + } + return false; + } + }; + + public: + /// Fit implementation + /// + /// @tparam source_link_iterator_t Iterator type used to pass source links + /// @tparam start_parameters_t Type of the initial parameters + /// @tparam parameters_t Type of parameters used for local parameters + /// @tparam track_container_t Type of the track container backend + /// @tparam holder_t Type defining track container backend ownership + /// + /// @param it Begin iterator for the fittable uncalibrated measurements + /// @param end End iterator for the fittable uncalibrated measurements + /// @param sParameters The initial track parameters + /// @param gx2fOptions Gx2FitterOptions steering the fit + /// @param trackContainer Input track container storage to append into + /// @note The input measurements are given in the form of @c SourceLink s. + /// It's the calibrators job to turn them into calibrated measurements used in + /// the fit. + /// + /// @return the output as an output track + template class holder_t, + bool _isdn = isDirectNavigator> + auto fit(source_link_iterator_t it, source_link_iterator_t end, + const start_parameters_t& sParameters, + const Gx2FitterOptions& gx2fOptions, + TrackContainer& trackContainer) + const -> std::enable_if_t< + !_isdn, Result::TrackProxy>> { + // Preprocess Measurements (Sourcelinks -> map) + // To be able to find measurements later, we put them into a map + // We need to copy input SourceLinks anyway, so the map can own them. + ACTS_VERBOSE("Preparing " << std::distance(it, end) + << " input measurements"); + std::map inputMeasurements; + + for (; it != end; ++it) { + SourceLink sl = *it; + auto geoId = sl.geometryId(); + inputMeasurements.emplace(geoId, std::move(sl)); + } + ACTS_VERBOSE("inputMeasurements.size() = " << inputMeasurements.size()); + + /// Fully understand Aborter, Actor, Result later + // Create the ActionList and AbortList + using GX2FAborter = Aborter; + using GX2FActor = Actor; + + using GX2FResult = typename GX2FActor::result_type; + using Actors = Acts::ActionList; + using Aborters = Acts::AbortList; + + using PropagatorOptions = Acts::PropagatorOptions; + + const size_t reducedMatrixSize = 4; + Acts::CurvilinearTrackParameters params = sParameters; + BoundVector deltaParams = BoundVector::Zero(); + double chi2sum = 0; + BoundMatrix aMatrix = BoundMatrix::Zero(); + BoundVector bVector = BoundVector::Zero(); + + ACTS_VERBOSE("params:\n" << params); + + /// Actual Fitting ///////////////////////////////////////////////////////// + ACTS_DEBUG("Start to iterate"); + + // Iterate the fit and improve result. Abort after n steps or after + // convergence + for (size_t nUpdate = 0; nUpdate < gx2fOptions.nUpdateMax; nUpdate++) { + ACTS_VERBOSE("nUpdate = " << nUpdate + 1 << "/" + << gx2fOptions.nUpdateMax); + + // update params + params.parameters() += deltaParams; + ACTS_VERBOSE("updated params:\n" << params); + + // set up propagator and co + Acts::GeometryContext geoCtx = gx2fOptions.geoContext; + Acts::MagneticFieldContext magCtx = gx2fOptions.magFieldContext; + // Set options for propagator + PropagatorOptions propagatorOptions(geoCtx, magCtx); + auto& gx2fActor = propagatorOptions.actionList.template get(); + gx2fActor.inputMeasurements = &inputMeasurements; + gx2fActor.extensions = gx2fOptions.extensions; + gx2fActor.actorLogger = m_actorLogger.get(); + + typename propagator_t::template action_list_t_result_t< + CurvilinearTrackParameters, Actors> + inputResult; + + auto& r = inputResult.template get>(); + + r.fittedStates = &trackContainer.trackStateContainer(); + // propagate with params and return jacobians and residuals + auto result = m_propagator.template propagate(params, propagatorOptions, + std::move(inputResult)); + + // TODO Improve Propagator + Actor [allocate before loop], rewrite + // makeMeasurements + auto& propRes = *result; + auto gx2fResult = std::move(propRes.template get()); + + ACTS_VERBOSE("gx2fResult.collectorResiduals.size() = " + << gx2fResult.collectorResiduals.size()); + ACTS_VERBOSE("gx2fResult.collectorCovariance.size() = " + << gx2fResult.collectorCovariance.size()); + ACTS_VERBOSE("gx2fResult.collectorJacobians.size() = " + << gx2fResult.collectorJacobians.size()); + + chi2sum = 0; + aMatrix = BoundMatrix::Zero(); + bVector = BoundVector::Zero(); + + // 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 projectedJacobian = + proj * gx2fResult.collectorJacobians[iMeas]; + + const double chi2meas = (ri.transpose() * coviInv * ri).eval()(0); + const BoundMatrix aMatrixMeas = + projectedJacobian.transpose() * coviInv * projectedJacobian; + const BoundVector bVectorMeas = + projectedJacobian.transpose() * coviInv * ri; + + chi2sum += chi2meas; + aMatrix += aMatrixMeas; + bVector += bVectorMeas; + } + + // calculate delta params [a] * delta = b + deltaParams = BoundVector::Zero(); + const ActsVector deltaParamsReduced = + aMatrix.topLeftCorner() + .colPivHouseholderQr() + .solve(bVector.topLeftCorner()); + + for (size_t idp = 0; idp < reducedMatrixSize; idp++) { + deltaParams(idp, 0) = deltaParamsReduced(idp, 0); + } + + ACTS_VERBOSE("chi2sum = " << chi2sum); + ACTS_VERBOSE("aMatrix:\n" << aMatrix); + ACTS_VERBOSE("bVector:\n" << bVector); + ACTS_VERBOSE("deltaParams:\n" << deltaParams); + + // TODO check delta params and abort + // similar to: + // if (sum(delta_params) < 1e-3) { + // break; + // } + } + ACTS_DEBUG("Finished to iterate"); + /// Finish Fitting ///////////////////////////////////////////////////////// + + // Calculate covariance of the fitted parameters with inverse of [a] + BoundMatrix fullCovariancePredicted = BoundMatrix::Identity(); + if (aMatrix.topLeftCorner() + .determinant() != 0) { + fullCovariancePredicted + .template topLeftCorner() = + aMatrix.topLeftCorner() + .inverse(); + } else if (gx2fOptions.nUpdateMax > 0) { + // TODO + std::cout << "det(a) == 0. This shouldn't happen. Implement real ERROR" + << std::endl; + } + + // Prepare track for return + auto track = trackContainer.getTrack(trackContainer.addTrack()); + track.parameters() = params.parameters(); + track.covariance() = fullCovariancePredicted; + // TODO track.tipIndex() = gx2fResult.lastMeasurementIndex; + // TODO track.setReferenceSurface(params.referenceSurface().getSharedPtr()); + // TODO track.nMeasurements() = gx2fResult.measurementStates; + // TODO track.nHoles() = gx2fResult.measurementHoles; + calculateTrackQuantities( + track); // TODO write test for calculateTrackQuantities + + // Return the converted Track + return track; + } +}; +} // namespace Experimental +} // namespace Acts diff --git a/Tests/UnitTests/Core/TrackFitting/CMakeLists.txt b/Tests/UnitTests/Core/TrackFitting/CMakeLists.txt index f2b4b97dabe..c415bbb571c 100644 --- a/Tests/UnitTests/Core/TrackFitting/CMakeLists.txt +++ b/Tests/UnitTests/Core/TrackFitting/CMakeLists.txt @@ -4,3 +4,4 @@ add_unittest(KalmanFitter KalmanFitterTests.cpp) add_unittest(Gsf GsfTests.cpp) add_unittest(GsfComponentMerging GsfComponentMergingTests.cpp) add_unittest(GsfMixtureReduction GsfMixtureReductionTests.cpp) +add_unittest(Gx2f Gx2fTests.cpp) diff --git a/Tests/UnitTests/Core/TrackFitting/Gx2fTests.cpp b/Tests/UnitTests/Core/TrackFitting/Gx2fTests.cpp new file mode 100644 index 00000000000..e423a14a7c9 --- /dev/null +++ b/Tests/UnitTests/Core/TrackFitting/Gx2fTests.cpp @@ -0,0 +1,355 @@ +// This file is part of the Acts project. +// +// Copyright (C) 2023 CERN for the benefit of the Acts project +// +// This Source Code Form is subject to the terms of the Mozilla Public +// License, v. 2.0. If a copy of the MPL was not distributed with this +// file, You can obtain one at http://mozilla.org/MPL/2.0/. + +#include + +#include "Acts/Definitions/Algebra.hpp" +#include "Acts/Definitions/Units.hpp" +#include "Acts/EventData/VectorMultiTrajectory.hpp" +#include "Acts/EventData/VectorTrackContainer.hpp" +#include "Acts/Geometry/CuboidVolumeBuilder.hpp" +#include "Acts/Geometry/Layer.hpp" +#include "Acts/Geometry/TrackingGeometry.hpp" +#include "Acts/Geometry/TrackingGeometryBuilder.hpp" +#include "Acts/MagneticField/ConstantBField.hpp" +#include "Acts/Material/HomogeneousSurfaceMaterial.hpp" +#include "Acts/Material/HomogeneousVolumeMaterial.hpp" +#include "Acts/Material/MaterialSlab.hpp" +#include "Acts/Propagator/EigenStepper.hpp" +#include "Acts/Propagator/Navigator.hpp" +#include "Acts/Propagator/Propagator.hpp" +#include "Acts/Propagator/StraightLineStepper.hpp" +#include "Acts/Surfaces/RectangleBounds.hpp" +#include "Acts/Tests/CommonHelpers/DetectorElementStub.hpp" +#include "Acts/Tests/CommonHelpers/MeasurementsCreator.hpp" +#include "Acts/Tests/CommonHelpers/PredefinedMaterials.hpp" +#include "Acts/TrackFitting/GainMatrixUpdater.hpp" +#include "Acts/TrackFitting/GlobalChiSquareFitter.hpp" +#include "Acts/Utilities/Logger.hpp" + +#include + +#include "FitterTestsCommon.hpp" + +using namespace Acts::UnitLiterals; + +Acts::Logging::Level logLevel = Acts::Logging::VERBOSE; +const auto gx2fLogger = Acts::getDefaultLogger("Gx2f", logLevel); + +namespace Acts { +namespace Test { + +//// Construct initial track parameters. +Acts::CurvilinearTrackParameters makeParameters( + const ActsScalar x = 0.0_m, const ActsScalar y = 0.0_m, + const ActsScalar z = 0.0_m, const ActsScalar w = 42_ns, + const ActsScalar phi = 0_degree, const ActsScalar theta = 90_degree, + const ActsScalar p = 1_GeV, const ActsScalar q = 1_e) { + // create covariance matrix from reasonable standard deviations + Acts::BoundVector stddev; + stddev[Acts::eBoundLoc0] = 100_um; + stddev[Acts::eBoundLoc1] = 100_um; + stddev[Acts::eBoundTime] = 25_ns; + stddev[Acts::eBoundPhi] = 2_degree; + stddev[Acts::eBoundTheta] = 2_degree; + stddev[Acts::eBoundQOverP] = 1 / 100_GeV; + Acts::BoundSquareMatrix cov = stddev.cwiseProduct(stddev).asDiagonal(); + // define a track in the transverse plane along x + Acts::Vector4 mPos4(x, y, z, w); + return Acts::CurvilinearTrackParameters(mPos4, phi, theta, p, q, cov); +} + +// Construct a straight-line propagator. +auto makeStraightPropagator(std::shared_ptr geo) { + Acts::Navigator::Config cfg{std::move(geo)}; + cfg.resolvePassive = false; + cfg.resolveMaterial = true; + cfg.resolveSensitive = true; + Acts::Navigator navigator(cfg); + Acts::StraightLineStepper stepper; + return Acts::Propagator( + stepper, std::move(navigator)); +} + +static std::vector prepareSourceLinks( + const std::vector& sourceLinks) { + std::vector result; + std::transform(sourceLinks.begin(), sourceLinks.end(), + std::back_inserter(result), + [](const auto& sl) { return Acts::SourceLink{sl}; }); + return result; +} + +std::shared_ptr makeToyDetector( + const GeometryContext& tgContext, const size_t nSurfaces = 5) { + if (nSurfaces < 1) { + throw std::invalid_argument("At least 1 surfaces needs to be created."); + } + // Construct builder + CuboidVolumeBuilder cvb; + + // Create configurations for surfaces + std::vector surfaceConfig; + for (unsigned int i = 1; i <= nSurfaces; i++) { + // Position of the surfaces + CuboidVolumeBuilder::SurfaceConfig cfg; + cfg.position = {i * UnitConstants::m, 0, 0.}; + + // Rotation of the surfaces + double rotationAngle = M_PI * 0.5; + Vector3 xPos(cos(rotationAngle), 0., sin(rotationAngle)); + Vector3 yPos(0., 1., 0.); + Vector3 zPos(-sin(rotationAngle), 0., cos(rotationAngle)); + cfg.rotation.col(0) = xPos; + cfg.rotation.col(1) = yPos; + cfg.rotation.col(2) = zPos; + /// Shape of the surface + // Boundaries of the surfaces + cfg.rBounds = + std::make_shared(RectangleBounds(0.5_m, 0.5_m)); + + // Material of the surfaces + MaterialSlab matProp(makeBeryllium(), 0.5_mm); + cfg.surMat = std::make_shared(matProp); + + // Thickness of the detector element + cfg.thickness = 1_um; + + cfg.detElementConstructor = + [](const Transform3& trans, + const std::shared_ptr& bounds, + double thickness) { + return new DetectorElementStub(trans, bounds, thickness); + }; + surfaceConfig.push_back(cfg); + } + + // Build layer configurations + std::vector layerConfig; + for (auto& sCfg : surfaceConfig) { + CuboidVolumeBuilder::LayerConfig cfg; + cfg.surfaceCfg = {sCfg}; + layerConfig.push_back(cfg); + } + + for (auto& cfg : layerConfig) { + cfg.surfaces = {}; + } + + // Inner Volume - Build volume configuration + CuboidVolumeBuilder::VolumeConfig volumeConfig; + volumeConfig.position = {2.5_m, 0., 0.}; + volumeConfig.length = {5_m, 1_m, 1_m}; + volumeConfig.layerCfg = layerConfig; + volumeConfig.name = "Test volume"; + volumeConfig.volumeMaterial = + std::make_shared(makeBeryllium()); + + volumeConfig.layers.clear(); + for (auto& lay : volumeConfig.layerCfg) { + lay.active = true; + } + + // Outer volume - Build TrackingGeometry configuration + CuboidVolumeBuilder::Config config; + config.position = {2.5_m, 0., 0.}; + config.length = {5_m, 1_m, 1_m}; + config.volumeCfg = {volumeConfig}; + + cvb.setConfig(config); + + TrackingGeometryBuilder::Config tgbCfg; + + tgbCfg.trackingVolumeBuilders.push_back( + [=](const auto& context, const auto& inner, const auto&) { + return cvb.trackingVolume(context, inner, nullptr); + }); + + TrackingGeometryBuilder tgb(tgbCfg); + + std::unique_ptr detector = + tgb.trackingGeometry(tgContext); + return detector; +} + +struct Detector { + // geometry + std::shared_ptr geometry; +}; + +BOOST_AUTO_TEST_SUITE(Gx2fTest) + +// This test checks if the call to the fitter works and no errors occur in the +// framework, without fitting and updating any parameters +BOOST_AUTO_TEST_CASE(NoFit) { + ACTS_LOCAL_LOGGER(Acts::getDefaultLogger("Gx2fTests", logLevel)); + ACTS_INFO("*** Test: NoFit -- Start"); + + // Context objects + Acts::GeometryContext geoCtx; + Acts::MagneticFieldContext magCtx; + Acts::CalibrationContext calCtx; + std::default_random_engine rng(42); + + Detector detector; + const size_t nSurfaces = 5; + detector.geometry = makeToyDetector(geoCtx, nSurfaces); + + auto parametersMeasurements = makeParameters(); + auto startParametersFit = makeParameters(0.1_m, 0.1_m, 0.1_m, 42_ns, + 10_degree, 80_degree, 1_GeV, 1_e); + + MeasurementResolution resPixel = {MeasurementType::eLoc01, {25_um, 50_um}}; + MeasurementResolutionMap resolutions = { + {Acts::GeometryIdentifier().setVolume(0), resPixel}}; + + // propagator + using SimPropagator = + Acts::Propagator; + SimPropagator simPropagator = makeStraightPropagator(detector.geometry); + auto measurements = createMeasurements( + simPropagator, geoCtx, magCtx, parametersMeasurements, resolutions, rng); + auto sourceLinks = prepareSourceLinks(measurements.sourceLinks); + + using Gx2Fitter = + Experimental::Gx2Fitter; + Gx2Fitter Fitter(simPropagator, gx2fLogger->clone()); + + const Surface* rSurface = ¶metersMeasurements.referenceSurface(); + + Experimental::Gx2FitterExtensions extensions; + extensions.calibrator + .connect<&testSourceLinkCalibrator>(); + + Experimental::Gx2FitterOptions gx2fOptions( + geoCtx, magCtx, calCtx, extensions, PropagatorPlainOptions(), rSurface, + false, false, FreeToBoundCorrection(false), 0); + + 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); + BOOST_CHECK_EQUAL(track.parameters(), startParametersFit.parameters()); + BOOST_CHECK_EQUAL(track.covariance(), BoundMatrix::Identity()); + + ACTS_INFO("*** Test: NoFit -- Finish"); +} + +BOOST_AUTO_TEST_CASE(Fit5Iterations) { + ACTS_LOCAL_LOGGER(Acts::getDefaultLogger("Gx2fTests", logLevel)); + ACTS_INFO("*** Test: Fit5Iterations -- 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(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, {100_um}}; + // MeasurementResolution resStrip1 = {MeasurementType::eLoc1, {150_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(); + + 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(Vector3(0., 0., 0.)); + using RecoStepper = EigenStepper<>; + RecoStepper rStepper(bField); + using RecoPropagator = Propagator; + RecoPropagator rPropagator(rStepper, rNavigator); + + using Gx2Fitter = + Experimental::Gx2Fitter; + Gx2Fitter Fitter(rPropagator, gx2fLogger->clone()); + + Experimental::Gx2FitterExtensions extensions; + extensions.calibrator + .connect<&testSourceLinkCalibrator>(); + + 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(), 1e-27, 4e0); + + ACTS_INFO("*** Test: Fit5Iterations -- Finish"); +} + +BOOST_AUTO_TEST_SUITE_END() +} // namespace Test +} // namespace Acts