From a22f0fcba2147815ab758aef29b7c64896d8b68f Mon Sep 17 00:00:00 2001 From: Xiaocong Ai <xiaocong.ai@cern.ch> Date: Wed, 26 Jun 2024 17:24:20 +0200 Subject: [PATCH] implement TrackFitterFunction and TrackFinderFunction to avoid duplicate code --- .../FaserActsGeometryContainers.h | 17 +- .../FaserActsKalmanFilter/IndexSourceLink.h | 14 +- .../src/FaserActsKalmanFilterAlg.cxx | 2 +- .../src/FaserActsKalmanFilterAlg.h | 49 +---- .../src/KalmanFitterTool.cxx | 125 +---------- .../src/KalmanFitterTool.h | 57 +---- .../src/TrackFinderFunction.cxx | 77 +++++++ .../src/TrackFinderFunction.h | 69 ++++++ .../src/TrackFitterFunction.cxx | 166 ++++++++++++++ .../src/TrackFitterFunction.h | 104 +++++++++ .../src/TrackFittingFunction.cxx | 207 ------------------ 11 files changed, 444 insertions(+), 443 deletions(-) create mode 100644 Tracking/Acts/FaserActsKalmanFilter/src/TrackFinderFunction.cxx create mode 100644 Tracking/Acts/FaserActsKalmanFilter/src/TrackFinderFunction.h create mode 100644 Tracking/Acts/FaserActsKalmanFilter/src/TrackFitterFunction.cxx create mode 100644 Tracking/Acts/FaserActsKalmanFilter/src/TrackFitterFunction.h delete mode 100644 Tracking/Acts/FaserActsKalmanFilter/src/TrackFittingFunction.cxx diff --git a/Tracking/Acts/FaserActsKalmanFilter/FaserActsKalmanFilter/FaserActsGeometryContainers.h b/Tracking/Acts/FaserActsKalmanFilter/FaserActsKalmanFilter/FaserActsGeometryContainers.h index 2f049c7d6..241c08a00 100644 --- a/Tracking/Acts/FaserActsKalmanFilter/FaserActsKalmanFilter/FaserActsGeometryContainers.h +++ b/Tracking/Acts/FaserActsKalmanFilter/FaserActsKalmanFilter/FaserActsGeometryContainers.h @@ -88,6 +88,7 @@ template <typename T> using GeometryIdMultimap = GeometryIdMultiset<std::pair<Acts::GeometryIdentifier, T>>; + /// The accessor for the GeometryIdMultiset container /// /// It wraps up a few lookup methods to be used in the Combinatorial Kalman @@ -101,20 +102,4 @@ struct GeometryIdMultisetAccessor { // pointer to the container const Container* container = nullptr; - - // count the number of elements with requested geoId - size_t count(const Acts::GeometryIdentifier& geoId) const { - assert(container != nullptr); - return container->count(geoId); - } - - // get the range of elements with requested geoId - std::pair<Iterator, Iterator> range( - const Acts::GeometryIdentifier& geoId) const { - assert(container != nullptr); - return container->equal_range(geoId); - } - - // get the element using the iterator - const Value& at(const Iterator& it) const { return *it; } }; diff --git a/Tracking/Acts/FaserActsKalmanFilter/FaserActsKalmanFilter/IndexSourceLink.h b/Tracking/Acts/FaserActsKalmanFilter/FaserActsKalmanFilter/IndexSourceLink.h index 68097833c..a9b59c5f5 100644 --- a/Tracking/Acts/FaserActsKalmanFilter/FaserActsKalmanFilter/IndexSourceLink.h +++ b/Tracking/Acts/FaserActsKalmanFilter/FaserActsKalmanFilter/IndexSourceLink.h @@ -97,4 +97,16 @@ using IndexSourceLinkContainer = GeometryIdMultiset<IndexSourceLink>; /// /// It wraps up a few lookup methods to be used in the Combinatorial Kalman /// Filter -using IndexSourceLinkAccessor = GeometryIdMultisetAccessor<IndexSourceLink>; +struct IndexSourceLinkAccessor : GeometryIdMultisetAccessor<IndexSourceLink> { + using BaseIterator = GeometryIdMultisetAccessor<IndexSourceLink>::Iterator; + + using Iterator = Acts::SourceLinkAdapterIterator<BaseIterator>; + + // get the range of elements with requested geoId + std::pair<Iterator, Iterator> range(const Acts::Surface& surface) const { + assert(container != nullptr); + auto [begin, end] = container->equal_range(surface.geometryId()); + return {Iterator{begin}, Iterator{end}}; + } +}; + diff --git a/Tracking/Acts/FaserActsKalmanFilter/src/FaserActsKalmanFilterAlg.cxx b/Tracking/Acts/FaserActsKalmanFilter/src/FaserActsKalmanFilterAlg.cxx index c72c34c51..7a8d541d8 100755 --- a/Tracking/Acts/FaserActsKalmanFilter/src/FaserActsKalmanFilterAlg.cxx +++ b/Tracking/Acts/FaserActsKalmanFilter/src/FaserActsKalmanFilterAlg.cxx @@ -144,7 +144,7 @@ StatusCode FaserActsKalmanFilterAlg::execute() { FaserActsTrackContainer tracks(actsTrackContainer, actsTrackStateContainer); //@todo: the initialSurface should be targetSurface - FaserActsKalmanFilterAlg::GeneralFitterOptions options{ + GeneralFitterOptions options{ geoctx, magctx, calctx, &(*initialSurface), Acts::PropagatorPlainOptions()}; diff --git a/Tracking/Acts/FaserActsKalmanFilter/src/FaserActsKalmanFilterAlg.h b/Tracking/Acts/FaserActsKalmanFilter/src/FaserActsKalmanFilterAlg.h index 5bc6cfd7f..9ef73935d 100755 --- a/Tracking/Acts/FaserActsKalmanFilter/src/FaserActsKalmanFilterAlg.h +++ b/Tracking/Acts/FaserActsKalmanFilter/src/FaserActsKalmanFilterAlg.h @@ -26,6 +26,7 @@ #include "TrkTrack/TrackCollection.h" //todo#include "TrajectoryWriterTool.h" + // ACTS #include "Acts/MagneticField/ConstantBField.hpp" #include "Acts/MagneticField/InterpolatedBFieldMap.hpp" @@ -50,9 +51,10 @@ #include "FaserActsGeometryInterfaces/IFaserActsTrackingGeometryTool.h" #include "FaserActsGeometryInterfaces/IFaserActsExtrapolationTool.h" #include "FaserActsKalmanFilter/IndexSourceLink.h" -#include "FaserActsKalmanFilter/Measurement.h" #include "FaserActsKalmanFilter/ITrackFinderTool.h" #include "FaserActsTrack.h" +#include "TrackFitterFunction.h" + //#include "ProtoTrackWriterTool.h" //todo #include "RootTrajectoryStatesWriterTool.h" @@ -66,6 +68,8 @@ class EventContext; class IAthRNGSvc; class FaserSCT_ID; class TTree; +class TrackFitterFunction; + namespace ActsExtrapolationDetail { class VariantPropagator; @@ -75,9 +79,8 @@ namespace TrackerDD { class SCT_DetectorManager; } -//@toberemoved -//using TrajectoryContainer = std::vector<FaserActsRecMultiTrajectory>; -//using BField_t = FASERMagneticFieldWrapper; +using TrackFitterResult = Acts::Result<FaserActsTrackContainer::TrackProxy>; + //class FaserActsKalmanFilterAlg : public AthReentrantAlgorithm { class FaserActsKalmanFilterAlg : public AthAlgorithm { @@ -90,44 +93,6 @@ public: StatusCode execute() override; StatusCode finalize() override; -//@toberemoved -// using IndexedParams = std::unordered_map<size_t, Acts::BoundTrackParameters>; -// using TrackFitterOptions = -// Acts::KalmanFitterOptions<MeasurementCalibrator, Acts::VoidOutlierFinder, -// Acts::VoidReverseFilteringLogic>; - - struct GeneralFitterOptions { - std::reference_wrapper<const Acts::GeometryContext> geoContext; - std::reference_wrapper<const Acts::MagneticFieldContext> magFieldContext; - std::reference_wrapper<const Acts::CalibrationContext> calibrationContext; - const Acts::Surface* referenceSurface = nullptr; - Acts::PropagatorPlainOptions propOptions; - }; - - using TrackFitterResult = Acts::Result<FaserActsTrackContainer::TrackProxy>; - - using TrackParameters = Acts::BoundTrackParameters; - - class TrackFitterFunction { - public: - virtual ~TrackFitterFunction() = default; - virtual TrackFitterResult operator()( - const std::vector<Acts::SourceLink> &sourceLinks, - const TrackParameters &initialParameters, - const GeneralFitterOptions& options, - const MeasurementCalibratorAdapter& calibrator, - FaserActsTrackContainer& tracks - ) const = 0; - }; - - static std::shared_ptr<TrackFitterFunction> makeTrackFitterFunction( - std::shared_ptr<const Acts::TrackingGeometry> trackingGeometry, - std::shared_ptr<const Acts::MagneticFieldProvider> magneticField, - bool multipleScattering, bool energyLoss, - double reverseFilteringMomThreshold, - Acts::FreeToBoundCorrection freeToBoundCorrection, - const Acts::Logger& logger); - virtual Acts::MagneticFieldContext getMagneticFieldContext(const EventContext& ctx) const; private: diff --git a/Tracking/Acts/FaserActsKalmanFilter/src/KalmanFitterTool.cxx b/Tracking/Acts/FaserActsKalmanFilter/src/KalmanFitterTool.cxx index b5f87b010..f5f64c970 100644 --- a/Tracking/Acts/FaserActsKalmanFilter/src/KalmanFitterTool.cxx +++ b/Tracking/Acts/FaserActsKalmanFilter/src/KalmanFitterTool.cxx @@ -84,7 +84,7 @@ KalmanFitterTool::getUnbiasedResidual(const EventContext &ctx, const Acts::Geome FaserActsTrackContainer tracks(actsTrackContainer, actsTrackStateContainer); //@todo: the initialSurface should be targetSurface - KalmanFitterTool::GeneralFitterOptions options{ + GeneralFitterOptions options{ gctx, mfContext, calibContext, &(*pSurface), Acts::PropagatorPlainOptions()}; @@ -180,7 +180,7 @@ KalmanFitterTool::getUnbiasedResidual(const EventContext &ctx, const Acts::Geome FaserActsTrackContainer tracks(actsTrackContainer, actsTrackStateContainer); //@todo: the initialSurface should be targetSurface - KalmanFitterTool::GeneralFitterOptions options{ + GeneralFitterOptions options{ gctx, mfContext, calibContext, &(*pSurface), Acts::PropagatorPlainOptions()}; @@ -279,7 +279,7 @@ KalmanFitterTool::getUnbiasedResidual(const EventContext &ctx, const Acts::Geome FaserActsTrackContainer tracks(actsTrackContainer, actsTrackStateContainer); //@todo: the initialSurface should be targetSurface - KalmanFitterTool::GeneralFitterOptions options{ + GeneralFitterOptions options{ gctx, mfContext, calibContext, &(*pSurface), Acts::PropagatorPlainOptions()}; @@ -377,7 +377,7 @@ KalmanFitterTool::fit(const EventContext &ctx, const Acts::GeometryContext &gctx FaserActsTrackContainer tracks(actsTrackContainer, actsTrackStateContainer); //@todo: the initialSurface should be targetSurface - KalmanFitterTool::GeneralFitterOptions options{ + GeneralFitterOptions options{ gctx, mfContext, calibContext, &(*pSurface), Acts::PropagatorPlainOptions()}; @@ -429,123 +429,6 @@ KalmanFitterTool::fit(const EventContext &ctx, const Acts::GeometryContext &gctx } -namespace { - -using Stepper = Acts::EigenStepper<>; -using Propagator = Acts::Propagator<Stepper, Acts::Navigator>; -using Fitter = Acts::KalmanFitter<Propagator, Acts::VectorMultiTrajectory>; - - -struct SimpleReverseFilteringLogic { - double momentumThreshold = 0; - - bool doBackwardFiltering( - Acts::VectorMultiTrajectory::ConstTrackStateProxy trackState) const { - auto momentum = fabs(1 / trackState.filtered()[Acts::eBoundQOverP]); - return (momentum <= momentumThreshold); - } -}; - - -struct TrackFitterFunctionImpl - : public KalmanFitterTool::TrackFitterFunction { - Fitter trackFitter; - - Acts::GainMatrixUpdater kfUpdater; - Acts::GainMatrixSmoother kfSmoother; - SimpleReverseFilteringLogic reverseFilteringLogic; - - bool multipleScattering = false; - bool energyLoss = false; - Acts::FreeToBoundCorrection freeToBoundCorrection; - - IndexSourceLink::SurfaceAccessor slSurfaceAccessor; - - TrackFitterFunctionImpl(Fitter &&f, const Acts::TrackingGeometry& trkGeo) : trackFitter(std::move(f)),slSurfaceAccessor{trkGeo} {} - - template <typename calibrator_t> - auto makeKfOptions(const KalmanFitterTool::GeneralFitterOptions& options, - const calibrator_t& calibrator, const FaserActsOutlierFinder& outlierFinder) const { - Acts::KalmanFitterExtensions<Acts::VectorMultiTrajectory> extensions; - extensions.updater.connect< - &Acts::GainMatrixUpdater::operator()<Acts::VectorMultiTrajectory>>( - &kfUpdater); - extensions.smoother.connect< - &Acts::GainMatrixSmoother::operator()<Acts::VectorMultiTrajectory>>( - &kfSmoother); - extensions.reverseFilteringLogic - .connect<&SimpleReverseFilteringLogic::doBackwardFiltering>( - &reverseFilteringLogic); - extensions.outlierFinder.connect<&FaserActsOutlierFinder::operator()<Acts::VectorMultiTrajectory>>( - &outlierFinder); - - Acts::KalmanFitterOptions<Acts::VectorMultiTrajectory> kfOptions( - options.geoContext, options.magFieldContext, options.calibrationContext, - extensions, options.propOptions, &(*options.referenceSurface)); - - kfOptions.referenceSurfaceStrategy = - Acts::KalmanFitterTargetSurfaceStrategy::first; - kfOptions.multipleScattering = multipleScattering; - kfOptions.energyLoss = energyLoss; - kfOptions.freeToBoundCorrection = freeToBoundCorrection; - kfOptions.extensions.calibrator.connect<&calibrator_t::calibrate>( - &calibrator); - kfOptions.extensions.surfaceAccessor - .connect<&IndexSourceLink::SurfaceAccessor::operator()>( - &slSurfaceAccessor); - - return kfOptions; - } - - KalmanFitterTool::TrackFitterResult operator()( - const std::vector<Acts::SourceLink> &sourceLinks, - const KalmanFitterTool::TrackParameters &initialParameters, - const KalmanFitterTool::GeneralFitterOptions& options, - const MeasurementCalibratorAdapter& calibrator, - const FaserActsOutlierFinder& outlierFinder, - FaserActsTrackContainer& tracks) const override { - const auto kfOptions = makeKfOptions(options, calibrator, outlierFinder); - return trackFitter.fit(sourceLinks.begin(), sourceLinks.end(), initialParameters, - kfOptions, tracks); - } -}; - -} // namespace - - -std::shared_ptr<KalmanFitterTool::TrackFitterFunction> -KalmanFitterTool::makeTrackFitterFunction( - std::shared_ptr<const Acts::TrackingGeometry> trackingGeometry, - std::shared_ptr<const Acts::MagneticFieldProvider> magneticField, - bool multipleScattering, bool energyLoss, - double reverseFilteringMomThreshold, - Acts::FreeToBoundCorrection freeToBoundCorrection, - const Acts::Logger& logger) { - // Stepper should be copied into the fitters - const Stepper stepper(std::move(magneticField)); - - // Standard fitter - const auto& geo = *trackingGeometry; - Acts::Navigator::Config cfg{std::move(trackingGeometry)}; - cfg.resolvePassive = false; - cfg.resolveMaterial = true; - cfg.resolveSensitive = true; - Acts::Navigator navigator(cfg, logger.cloneWithSuffix("Navigator")); - Propagator propagator(stepper, std::move(navigator), - logger.cloneWithSuffix("Propagator")); - Fitter trackFitter(std::move(propagator), logger.cloneWithSuffix("Fitter")); - - // build the fitter function. owns the fitter object. - auto fitterFunction = std::make_shared<TrackFitterFunctionImpl>( - std::move(trackFitter), geo); - fitterFunction->multipleScattering = multipleScattering; - fitterFunction->energyLoss = energyLoss; - fitterFunction->reverseFilteringLogic.momentumThreshold = - reverseFilteringMomThreshold; - fitterFunction->freeToBoundCorrection = freeToBoundCorrection; - - return fitterFunction; -} Acts::MagneticFieldContext KalmanFitterTool::getMagneticFieldContext(const EventContext& ctx) const { SG::ReadCondHandle<FaserFieldCacheCondObj> readHandle{m_fieldCondObjInputKey, ctx}; diff --git a/Tracking/Acts/FaserActsKalmanFilter/src/KalmanFitterTool.h b/Tracking/Acts/FaserActsKalmanFilter/src/KalmanFitterTool.h index b58a31269..8484d37d8 100644 --- a/Tracking/Acts/FaserActsKalmanFilter/src/KalmanFitterTool.h +++ b/Tracking/Acts/FaserActsKalmanFilter/src/KalmanFitterTool.h @@ -21,6 +21,8 @@ #include "FaserActsTrack.h" #include "FaserActsGeometry/FASERMagneticFieldWrapper.h" +#include "TrackFitterFunction.h" + struct TSOS4Residual{ double fit_local_x; double fit_local_y; @@ -36,29 +38,6 @@ struct TSOS4Residual{ class FaserSCT_ID; -//set the cluster to be removed as outlier in order to get the unbiased residual -/// Outlier finder using a Chi2 cut. -struct FaserActsOutlierFinder { - double StateChiSquaredPerNumberDoFCut = 10000.; - double cluster_z = -10000.; - - template <typename traj_t> - bool operator()(typename traj_t::ConstTrackStateProxy state) const { - //remove the whole IFT - if(cluster_z<-10000){ - IndexSourceLink sl = state.getUncalibratedSourceLink().template get<IndexSourceLink>(); - if(sl.hit()->globalPosition().z()<-100)return true; - if(abs(sl.hit()->globalPosition().z()-cluster_z)<3)return true; - } - - if (not state.hasCalibrated() or not state.hasPredicted()) { - return false; - } - - return bool(state.chi2() > StateChiSquaredPerNumberDoFCut * state.calibratedSize()); - } -}; - class KalmanFitterTool : virtual public AthAlgTool { public: @@ -67,38 +46,6 @@ public: virtual StatusCode initialize() override; virtual StatusCode finalize() override; - struct GeneralFitterOptions { - std::reference_wrapper<const Acts::GeometryContext> geoContext; - std::reference_wrapper<const Acts::MagneticFieldContext> magFieldContext; - std::reference_wrapper<const Acts::CalibrationContext> calibrationContext; - const Acts::Surface* referenceSurface = nullptr; - Acts::PropagatorPlainOptions propOptions; - }; - using TrackFitterResult = Acts::Result<FaserActsTrackContainer::TrackProxy>; - using TrackParameters = Acts::BoundTrackParameters; - //using IndexedParams = std::unordered_map<size_t, TrackParameters>; - - class TrackFitterFunction { - public: - virtual ~TrackFitterFunction() = default; - virtual TrackFitterResult operator()( - const std::vector<Acts::SourceLink> &sourceLinks, - const TrackParameters &initialParameters, - const GeneralFitterOptions & kfOptions, - const MeasurementCalibratorAdapter& calibrator, - const FaserActsOutlierFinder& outlierFinder, - FaserActsTrackContainer& tracks - ) const = 0; - }; - - static std::shared_ptr<TrackFitterFunction> makeTrackFitterFunction( - std::shared_ptr<const Acts::TrackingGeometry> trackingGeometry, - std::shared_ptr<const Acts::MagneticFieldProvider> magneticField, - bool multipleScattering, bool energyLoss, - double reverseFilteringMomThreshold, - Acts::FreeToBoundCorrection freeToBoundCorrection, - const Acts::Logger& logger); - virtual Acts::MagneticFieldContext getMagneticFieldContext(const EventContext& ctx) const; std::unique_ptr<Trk::Track> fit(const EventContext &ctx, const Acts::GeometryContext &gctx, diff --git a/Tracking/Acts/FaserActsKalmanFilter/src/TrackFinderFunction.cxx b/Tracking/Acts/FaserActsKalmanFilter/src/TrackFinderFunction.cxx new file mode 100644 index 000000000..cbdfef9c9 --- /dev/null +++ b/Tracking/Acts/FaserActsKalmanFilter/src/TrackFinderFunction.cxx @@ -0,0 +1,77 @@ +#include "TrackFinderFunction.h" + + +struct TrackFinderFunctionImpl + final : public TrackFinderFunction { + Acts::GainMatrixUpdater kfUpdater; + Acts::GainMatrixSmoother kfSmoother; + + CKF trackFinder; + + TrackFinderFunctionImpl(CKF &&f) : trackFinder(std::move(f)) {} + + template <typename calibrator_t> + auto makeCkfOptions(const GeneralFitterOptions& options, + const calibrator_t& calibrator, const IndexSourceLinkAccessor& slAccessor, const Acts::MeasurementSelector& measSel) const { + + Acts::CombinatorialKalmanFilterExtensions<Acts::VectorMultiTrajectory> + extensions; + extensions.calibrator.connect<&MeasurementCalibratorAdapter::calibrate>( + &calibrator); + extensions.updater.connect< + &Acts::GainMatrixUpdater::operator()<Acts::VectorMultiTrajectory>>( + &kfUpdater); + extensions.smoother.connect< + &Acts::GainMatrixSmoother::operator()<Acts::VectorMultiTrajectory>>( + &kfSmoother); + extensions.measurementSelector + .connect<&Acts::MeasurementSelector::select<Acts::VectorMultiTrajectory>>( + &measSel); + + Acts::SourceLinkAccessorDelegate<IndexSourceLinkAccessor::Iterator> + slAccessorDelegate; + slAccessorDelegate.connect<&IndexSourceLinkAccessor::range>(&slAccessor); + + TrackFinderOptions ckfOptions(options.geoContext, options.magFieldContext, options.calibrationContext, slAccessorDelegate, + extensions, options.propOptions, &(*options.referenceSurface)); + + //todo: make this configurable + ckfOptions.smoothingTargetSurfaceStrategy = + Acts::CombinatorialKalmanFilterTargetSurfaceStrategy::first; + + return ckfOptions; + } + + TrackFinderResult operator()( + const Acts::BoundTrackParameters& initialParameters, + const GeneralFitterOptions& options, + const MeasurementCalibratorAdapter& calibrator, + const IndexSourceLinkAccessor& slAccessor, + const Acts::MeasurementSelector& measSel, + FaserActsTrackContainer& tracks) const override { + auto ckfOptions = makeCkfOptions(options, calibrator, slAccessor, measSel); + return trackFinder.findTracks(initialParameters, ckfOptions, tracks); + } +}; + + + +std::shared_ptr<TrackFinderFunction> makeTrackFinderFunction( + std::shared_ptr<const Acts::TrackingGeometry> trackingGeometry, + std::shared_ptr<const Acts::MagneticFieldProvider> magneticField, + bool resolvePassive, bool resolveMaterial, bool resolveSensitive, + const Acts::Logger& logger) { + Stepper stepper(std::move(magneticField)); + Navigator::Config cfg{std::move(trackingGeometry)}; + cfg.resolvePassive = resolvePassive; + cfg.resolveMaterial = resolveMaterial; + cfg.resolveSensitive = resolveSensitive; + Navigator navigator(cfg, logger.cloneWithSuffix("Navigator")); + Propagator propagator(std::move(stepper), std::move(navigator), + logger.cloneWithSuffix("Propagator")); + CKF trackFinder(std::move(propagator), logger.cloneWithSuffix("Finder")); + + // build the track finder functions. owns the track finder object. + return std::make_shared<TrackFinderFunctionImpl>(std::move(trackFinder)); +} + diff --git a/Tracking/Acts/FaserActsKalmanFilter/src/TrackFinderFunction.h b/Tracking/Acts/FaserActsKalmanFilter/src/TrackFinderFunction.h new file mode 100644 index 000000000..8d879d159 --- /dev/null +++ b/Tracking/Acts/FaserActsKalmanFilter/src/TrackFinderFunction.h @@ -0,0 +1,69 @@ +#ifndef FASERACTSKALMANFILTER_TRACKFINDERFUNCTION_H +#define FASERACTSKALMANFILTER_TRACKFINDERFUNCTION_H + + +#include "FaserActsTrack.h" +#include "TrackFitterFunction.h" +#include "FaserActsKalmanFilter/IndexSourceLink.h" +#include "FaserActsKalmanFilter/Measurement.h" + + +#include "Acts/Definitions/Direction.hpp" +#include "Acts/EventData/MultiTrajectory.hpp" +#include "Acts/EventData/TrackContainer.hpp" +#include "Acts/EventData/TrackStatePropMask.hpp" +#include "Acts/EventData/VectorMultiTrajectory.hpp" +#include "Acts/EventData/VectorTrackContainer.hpp" +#include "Acts/Propagator/EigenStepper.hpp" +#include "Acts/Propagator/Navigator.hpp" +#include "Acts/Propagator/Propagator.hpp" +#include "Acts/TrackFinding/CombinatorialKalmanFilter.hpp" +#include "Acts/TrackFitting/GainMatrixSmoother.hpp" +#include "Acts/TrackFitting/GainMatrixUpdater.hpp" +#include "Acts/TrackFinding/MeasurementSelector.hpp" +#include "Acts/Utilities/Logger.hpp" + +//namespace { + +using Updater = Acts::GainMatrixUpdater; +using Smoother = Acts::GainMatrixSmoother; + +using Stepper = Acts::EigenStepper<>; +using Navigator = Acts::Navigator; +using Propagator = Acts::Propagator<Stepper, Navigator>; +using CKF = + Acts::CombinatorialKalmanFilter<Propagator, Acts::VectorMultiTrajectory>; + + +using TrackFinderOptions = + Acts::CombinatorialKalmanFilterOptions<IndexSourceLinkAccessor::Iterator, + Acts::VectorMultiTrajectory>; +using TrackFinderResult = + Acts::Result<std::vector<FaserActsTrackContainer::TrackProxy>>; + + +class TrackFinderFunction { + public: + + virtual ~TrackFinderFunction() = default; + virtual TrackFinderResult operator()( + const Acts::BoundTrackParameters& initialParameters, + const GeneralFitterOptions& options, + const MeasurementCalibratorAdapter& calibrator, + const IndexSourceLinkAccessor& slAccessor, + const Acts::MeasurementSelector& measSel, + FaserActsTrackContainer& tracks) const = 0; +}; + + + std::shared_ptr<TrackFinderFunction> makeTrackFinderFunction( + std::shared_ptr<const Acts::TrackingGeometry> trackingGeometry, + std::shared_ptr<const Acts::MagneticFieldProvider> magneticField, + bool resolvePassive, bool resolveMaterial, bool resolveSensitive, + const Acts::Logger& logger); + + +//} // namespace + +#endif // FASERACTSKALMANFILTER_TRACKFINDERFUNCTION_H + diff --git a/Tracking/Acts/FaserActsKalmanFilter/src/TrackFitterFunction.cxx b/Tracking/Acts/FaserActsKalmanFilter/src/TrackFitterFunction.cxx new file mode 100644 index 000000000..ad4d40e53 --- /dev/null +++ b/Tracking/Acts/FaserActsKalmanFilter/src/TrackFitterFunction.cxx @@ -0,0 +1,166 @@ +#include "TrackFitterFunction.h" +#include "TrackerPrepRawData/FaserSCT_Cluster.h" + + +template <typename traj_t> +bool FaserActsOutlierFinder::operator()(typename traj_t::ConstTrackStateProxy state) const { + //remove the whole IFT + if(cluster_z<-10000){ + IndexSourceLink sl = state.getUncalibratedSourceLink().template get<IndexSourceLink>(); + if(sl.hit()->globalPosition().z()<-100)return true; + if(abs(sl.hit()->globalPosition().z()-cluster_z)<3)return true; + } + + if (not state.hasCalibrated() or not state.hasPredicted()) { + return false; + } + + return bool(state.chi2() > StateChiSquaredPerNumberDoFCut * state.calibratedSize()); + } + + +struct TrackFitterFunctionImpl final : public TrackFitterFunction { + Fitter trackFitter; + + Acts::GainMatrixUpdater kfUpdater; + Acts::GainMatrixSmoother kfSmoother; + SimpleReverseFilteringLogic reverseFilteringLogic; + + bool multipleScattering = false; + bool energyLoss = false; + Acts::FreeToBoundCorrection freeToBoundCorrection; + + IndexSourceLink::SurfaceAccessor slSurfaceAccessor; + + TrackFitterFunctionImpl(Fitter &&f, const Acts::TrackingGeometry& trkGeo) : trackFitter(std::move(f)),slSurfaceAccessor{trkGeo} {} + + template <typename calibrator_t> + auto makeKfOptions(const GeneralFitterOptions& options, + const calibrator_t& calibrator) const { + Acts::KalmanFitterExtensions<Acts::VectorMultiTrajectory> extensions; + extensions.updater.connect< + &Acts::GainMatrixUpdater::operator()<Acts::VectorMultiTrajectory>>( + &kfUpdater); + extensions.smoother.connect< + &Acts::GainMatrixSmoother::operator()<Acts::VectorMultiTrajectory>>( + &kfSmoother); + extensions.reverseFilteringLogic + .connect<&SimpleReverseFilteringLogic::doBackwardFiltering>( + &reverseFilteringLogic); + + Acts::KalmanFitterOptions<Acts::VectorMultiTrajectory> kfOptions( + options.geoContext, options.magFieldContext, options.calibrationContext, + extensions, options.propOptions, &(*options.referenceSurface)); + + kfOptions.referenceSurfaceStrategy = + Acts::KalmanFitterTargetSurfaceStrategy::first; + kfOptions.multipleScattering = multipleScattering; + kfOptions.energyLoss = energyLoss; + kfOptions.freeToBoundCorrection = freeToBoundCorrection; + kfOptions.extensions.calibrator.connect<&calibrator_t::calibrate>( + &calibrator); + kfOptions.extensions.surfaceAccessor + .connect<&IndexSourceLink::SurfaceAccessor::operator()>( + &slSurfaceAccessor); + + return kfOptions; + } + + + template <typename calibrator_t> + auto makeKfOptions(const GeneralFitterOptions& options, + const calibrator_t& calibrator, const FaserActsOutlierFinder& outlierFinder) const { + Acts::KalmanFitterExtensions<Acts::VectorMultiTrajectory> extensions; + extensions.updater.connect< + &Acts::GainMatrixUpdater::operator()<Acts::VectorMultiTrajectory>>( + &kfUpdater); + extensions.smoother.connect< + &Acts::GainMatrixSmoother::operator()<Acts::VectorMultiTrajectory>>( + &kfSmoother); + extensions.reverseFilteringLogic + .connect<&SimpleReverseFilteringLogic::doBackwardFiltering>( + &reverseFilteringLogic); + extensions.outlierFinder.connect<&FaserActsOutlierFinder::operator()<Acts::VectorMultiTrajectory>>( + &outlierFinder); + + Acts::KalmanFitterOptions<Acts::VectorMultiTrajectory> kfOptions( + options.geoContext, options.magFieldContext, options.calibrationContext, + extensions, options.propOptions, &(*options.referenceSurface)); + + kfOptions.referenceSurfaceStrategy = + Acts::KalmanFitterTargetSurfaceStrategy::first; + kfOptions.multipleScattering = multipleScattering; + kfOptions.energyLoss = energyLoss; + kfOptions.freeToBoundCorrection = freeToBoundCorrection; + kfOptions.extensions.calibrator.connect<&calibrator_t::calibrate>( + &calibrator); + kfOptions.extensions.surfaceAccessor + .connect<&IndexSourceLink::SurfaceAccessor::operator()>( + &slSurfaceAccessor); + + return kfOptions; + } + + + //fit with default outlier finder, i.e. not finding outlier + TrackFitterResult operator()( + const std::vector<Acts::SourceLink> &sourceLinks, + const Acts::BoundTrackParameters &initialParameters, + const GeneralFitterOptions& options, + const MeasurementCalibratorAdapter& calibrator, + FaserActsTrackContainer& tracks) const override { + const auto kfOptions = makeKfOptions(options, calibrator); + return trackFitter.fit(sourceLinks.begin(), sourceLinks.end(), initialParameters, + kfOptions, tracks); + } + + + //fit with faser outlier finder + TrackFitterResult operator()( + const std::vector<Acts::SourceLink> &sourceLinks, + const Acts::BoundTrackParameters &initialParameters, + const GeneralFitterOptions& options, + const MeasurementCalibratorAdapter& calibrator, + const FaserActsOutlierFinder& outlierFinder, + FaserActsTrackContainer& tracks) const override { + const auto kfOptions = makeKfOptions(options, calibrator, outlierFinder); + return trackFitter.fit(sourceLinks.begin(), sourceLinks.end(), initialParameters, + kfOptions, tracks); + } + +}; + + +std::shared_ptr<TrackFitterFunction> makeTrackFitterFunction( + std::shared_ptr<const Acts::TrackingGeometry> trackingGeometry, + std::shared_ptr<const Acts::MagneticFieldProvider> magneticField, + bool multipleScattering, bool energyLoss, + double reverseFilteringMomThreshold, + Acts::FreeToBoundCorrection freeToBoundCorrection, + const Acts::Logger& logger) { + // Stepper should be copied into the fitters + const Stepper stepper(std::move(magneticField)); + + // Standard fitter + const auto& geo = *trackingGeometry; + Acts::Navigator::Config cfg{std::move(trackingGeometry)}; + cfg.resolvePassive = false; + cfg.resolveMaterial = true; + cfg.resolveSensitive = true; + Acts::Navigator navigator(cfg, logger.cloneWithSuffix("Navigator")); + Propagator propagator(stepper, std::move(navigator), + logger.cloneWithSuffix("Propagator")); + Fitter trackFitter(std::move(propagator), logger.cloneWithSuffix("Fitter")); + + // build the fitter function. owns the fitter object. + auto fitterFunction = std::make_shared<TrackFitterFunctionImpl>( + std::move(trackFitter), geo); + fitterFunction->multipleScattering = multipleScattering; + fitterFunction->energyLoss = energyLoss; + fitterFunction->reverseFilteringLogic.momentumThreshold = + reverseFilteringMomThreshold; + fitterFunction->freeToBoundCorrection = freeToBoundCorrection; + + return fitterFunction; +} + diff --git a/Tracking/Acts/FaserActsKalmanFilter/src/TrackFitterFunction.h b/Tracking/Acts/FaserActsKalmanFilter/src/TrackFitterFunction.h new file mode 100644 index 000000000..855611fcd --- /dev/null +++ b/Tracking/Acts/FaserActsKalmanFilter/src/TrackFitterFunction.h @@ -0,0 +1,104 @@ +#ifndef FASERACTSKALMANFILTER_TRACKFITTERFUNCTION_H +#define FASERACTSKALMANFILTER_TRACKFITTERFUNCTION_H + + +#include "FaserActsTrack.h" +#include "FaserActsKalmanFilter/IndexSourceLink.h" +#include "FaserActsKalmanFilter/Measurement.h" + +#include "Acts/Definitions/Direction.hpp" +#include "Acts/Definitions/TrackParametrization.hpp" +#include "Acts/EventData/MultiTrajectory.hpp" +#include "Acts/EventData/TrackContainer.hpp" +#include "Acts/EventData/TrackStatePropMask.hpp" +#include "Acts/EventData/VectorMultiTrajectory.hpp" +#include "Acts/EventData/VectorTrackContainer.hpp" +#include "Acts/EventData/detail/CorrectedTransformationFreeToBound.hpp" +#include "Acts/Geometry/GeometryIdentifier.hpp" +#include "Acts/Propagator/DirectNavigator.hpp" +#include "Acts/Propagator/EigenStepper.hpp" +#include "Acts/Propagator/Navigator.hpp" +#include "Acts/Propagator/Propagator.hpp" +#include "Acts/TrackFitting/GainMatrixSmoother.hpp" +#include "Acts/TrackFitting/GainMatrixUpdater.hpp" +#include "Acts/TrackFitting/KalmanFitter.hpp" +#include "Acts/Utilities/Delegate.hpp" +#include "Acts/Utilities/Logger.hpp" + +//#include "TrackerPrepRawData/FaserSCT_Cluster.h" + +//namespace { + +using Stepper = Acts::EigenStepper<>; +using Propagator = Acts::Propagator<Stepper, Acts::Navigator>; +using Fitter = Acts::KalmanFitter<Propagator, Acts::VectorMultiTrajectory>; + + +struct GeneralFitterOptions { + std::reference_wrapper<const Acts::GeometryContext> geoContext; + std::reference_wrapper<const Acts::MagneticFieldContext> magFieldContext; + std::reference_wrapper<const Acts::CalibrationContext> calibrationContext; + const Acts::Surface* referenceSurface = nullptr; + Acts::PropagatorPlainOptions propOptions; +}; + + +struct SimpleReverseFilteringLogic { + double momentumThreshold = 0; + + bool doBackwardFiltering( + Acts::VectorMultiTrajectory::ConstTrackStateProxy trackState) const { + auto momentum = fabs(1 / trackState.filtered()[Acts::eBoundQOverP]); + return (momentum <= momentumThreshold); + } +}; + + +//set the cluster to be removed as outlier in order to get the unbiased residual +/// Outlier finder using a Chi2 cut. +struct FaserActsOutlierFinder { + double StateChiSquaredPerNumberDoFCut = 10000.; + double cluster_z = -10000.; + + template <typename traj_t> + bool operator()(typename traj_t::ConstTrackStateProxy state) const; +}; + + +/// Fit function that takes the above parameters and runs a fit +/// @note This is separated into a virtual interface to keep compilation units +/// small. +class TrackFitterFunction { + public: + using TrackFitterResult = Acts::Result<FaserActsTrackContainer::TrackProxy>; + + virtual ~TrackFitterFunction() = default; + + virtual TrackFitterResult operator()(const std::vector<Acts::SourceLink>&, + const Acts::BoundTrackParameters&, + const GeneralFitterOptions&, + const MeasurementCalibratorAdapter&, + FaserActsTrackContainer&) const = 0; + + virtual TrackFitterResult operator()(const std::vector<Acts::SourceLink>&, + const Acts::BoundTrackParameters&, + const GeneralFitterOptions&, + const MeasurementCalibratorAdapter&, + const FaserActsOutlierFinder& outlierFinder, + FaserActsTrackContainer&) const = 0; + +}; + + +std::shared_ptr<TrackFitterFunction> makeTrackFitterFunction( +std::shared_ptr<const Acts::TrackingGeometry> trackingGeometry, +std::shared_ptr<const Acts::MagneticFieldProvider> magneticField, +bool multipleScattering, bool energyLoss, +double reverseFilteringMomThreshold, +Acts::FreeToBoundCorrection freeToBoundCorrection, +const Acts::Logger& logger); + + +//} // namespace + +#endif // FASERACTSKALMANFILTER_TRACKFITTERFUNCTION_H diff --git a/Tracking/Acts/FaserActsKalmanFilter/src/TrackFittingFunction.cxx b/Tracking/Acts/FaserActsKalmanFilter/src/TrackFittingFunction.cxx deleted file mode 100644 index 7a1e31dfc..000000000 --- a/Tracking/Acts/FaserActsKalmanFilter/src/TrackFittingFunction.cxx +++ /dev/null @@ -1,207 +0,0 @@ -#include "FaserActsKalmanFilterAlg.h" -#include "FaserActsTrack.h" -//#include "FaserActsGeometry/FASERMagneticFieldWrapper.h" -//#include "FaserActsKalmanFilter/IndexSourceLink.h" - -#include "Acts/Definitions/Direction.hpp" -#include "Acts/Definitions/TrackParametrization.hpp" -#include "Acts/EventData/MultiTrajectory.hpp" -#include "Acts/EventData/TrackContainer.hpp" -#include "Acts/EventData/TrackStatePropMask.hpp" -#include "Acts/EventData/VectorMultiTrajectory.hpp" -#include "Acts/EventData/VectorTrackContainer.hpp" -#include "Acts/EventData/detail/CorrectedTransformationFreeToBound.hpp" -#include "Acts/Geometry/GeometryIdentifier.hpp" -#include "Acts/Propagator/DirectNavigator.hpp" -#include "Acts/Propagator/EigenStepper.hpp" -#include "Acts/Propagator/Navigator.hpp" -#include "Acts/Propagator/Propagator.hpp" -#include "Acts/TrackFitting/GainMatrixSmoother.hpp" -#include "Acts/TrackFitting/GainMatrixUpdater.hpp" -#include "Acts/TrackFitting/KalmanFitter.hpp" -#include "Acts/Utilities/Delegate.hpp" -#include "Acts/Utilities/Logger.hpp" - - -namespace { - -using Stepper = Acts::EigenStepper<>; -using Propagator = Acts::Propagator<Stepper, Acts::Navigator>; -using Fitter = Acts::KalmanFitter<Propagator, Acts::VectorMultiTrajectory>; - - -struct SimpleReverseFilteringLogic { - double momentumThreshold = 0; - - bool doBackwardFiltering( - Acts::VectorMultiTrajectory::ConstTrackStateProxy trackState) const { - auto momentum = fabs(1 / trackState.filtered()[Acts::eBoundQOverP]); - return (momentum <= momentumThreshold); - } -}; - -struct TrackFitterFunctionImpl - : public FaserActsKalmanFilterAlg::TrackFitterFunction { - Fitter trackFitter; - - Acts::GainMatrixUpdater kfUpdater; - Acts::GainMatrixSmoother kfSmoother; - SimpleReverseFilteringLogic reverseFilteringLogic; - - bool multipleScattering = false; - bool energyLoss = false; - Acts::FreeToBoundCorrection freeToBoundCorrection; - - IndexSourceLink::SurfaceAccessor slSurfaceAccessor; - - TrackFitterFunctionImpl(Fitter &&f, const Acts::TrackingGeometry& trkGeo) : trackFitter(std::move(f)),slSurfaceAccessor{trkGeo} {} - - template <typename calibrator_t> - auto makeKfOptions(const FaserActsKalmanFilterAlg::GeneralFitterOptions& options, - const calibrator_t& calibrator) const { - Acts::KalmanFitterExtensions<Acts::VectorMultiTrajectory> extensions; - extensions.updater.connect< - &Acts::GainMatrixUpdater::operator()<Acts::VectorMultiTrajectory>>( - &kfUpdater); - extensions.smoother.connect< - &Acts::GainMatrixSmoother::operator()<Acts::VectorMultiTrajectory>>( - &kfSmoother); - extensions.reverseFilteringLogic - .connect<&SimpleReverseFilteringLogic::doBackwardFiltering>( - &reverseFilteringLogic); - - Acts::KalmanFitterOptions<Acts::VectorMultiTrajectory> kfOptions( - options.geoContext, options.magFieldContext, options.calibrationContext, - extensions, options.propOptions, &(*options.referenceSurface)); - - kfOptions.referenceSurfaceStrategy = - Acts::KalmanFitterTargetSurfaceStrategy::first; - kfOptions.multipleScattering = multipleScattering; - kfOptions.energyLoss = energyLoss; - kfOptions.freeToBoundCorrection = freeToBoundCorrection; - kfOptions.extensions.calibrator.connect<&calibrator_t::calibrate>( - &calibrator); - kfOptions.extensions.surfaceAccessor - .connect<&IndexSourceLink::SurfaceAccessor::operator()>( - &slSurfaceAccessor); - - return kfOptions; - } - - FaserActsKalmanFilterAlg::TrackFitterResult operator()( - const std::vector<Acts::SourceLink> &sourceLinks, - const FaserActsKalmanFilterAlg::TrackParameters &initialParameters, - const FaserActsKalmanFilterAlg::GeneralFitterOptions& options, - const MeasurementCalibratorAdapter& calibrator, - FaserActsTrackContainer& tracks) const override { - const auto kfOptions = makeKfOptions(options, calibrator); - return trackFitter.fit(sourceLinks.begin(), sourceLinks.end(), initialParameters, - kfOptions, tracks); - } -}; - -} // namespace - - -std::shared_ptr<FaserActsKalmanFilterAlg::TrackFitterFunction> -FaserActsKalmanFilterAlg::makeTrackFitterFunction( - std::shared_ptr<const Acts::TrackingGeometry> trackingGeometry, - std::shared_ptr<const Acts::MagneticFieldProvider> magneticField, - bool multipleScattering, bool energyLoss, - double reverseFilteringMomThreshold, - Acts::FreeToBoundCorrection freeToBoundCorrection, - const Acts::Logger& logger) { - // Stepper should be copied into the fitters - const Stepper stepper(std::move(magneticField)); - - // Standard fitter - const auto& geo = *trackingGeometry; - Acts::Navigator::Config cfg{std::move(trackingGeometry)}; - cfg.resolvePassive = false; - cfg.resolveMaterial = true; - cfg.resolveSensitive = true; - Acts::Navigator navigator(cfg, logger.cloneWithSuffix("Navigator")); - Propagator propagator(stepper, std::move(navigator), - logger.cloneWithSuffix("Propagator")); - Fitter trackFitter(std::move(propagator), logger.cloneWithSuffix("Fitter")); - - // build the fitter function. owns the fitter object. - auto fitterFunction = std::make_shared<TrackFitterFunctionImpl>( - std::move(trackFitter), geo); - fitterFunction->multipleScattering = multipleScattering; - fitterFunction->energyLoss = energyLoss; - fitterFunction->reverseFilteringLogic.momentumThreshold = - reverseFilteringMomThreshold; - fitterFunction->freeToBoundCorrection = freeToBoundCorrection; - - return fitterFunction; -} - -// The following are obsolete -/* - -namespace ActsExtrapolationDetail { -using VariantPropagatorBase = boost::variant< - Acts::Propagator<Acts::EigenStepper<>, Acts::Navigator>, - Acts::Propagator<Acts::EigenStepper<>, Acts::Navigator>>; - -class VariantPropagator : public VariantPropagatorBase { -public: - using VariantPropagatorBase::VariantPropagatorBase; -}; -} // namespace ActsExtrapolationDetail - - -namespace { -template <typename Fitter> -struct FitterFunctionImpl { - Fitter fitter; - FitterFunctionImpl(Fitter&& f) : fitter(std::move(f)) {} - FaserActsKalmanFilterAlg::TrackFitterResult operator()( - const std::vector<IndexSourceLink>& sourceLinks, - const Acts::CurvilinearTrackParameters& initialParameters, - const Acts::KalmanFitterOptions<MeasurementCalibrator, Acts::VoidOutlierFinder, Acts::VoidReverseFilteringLogic>& options) const { - return fitter.fit(sourceLinks, initialParameters, options); - }; -}; -} // namespace - -std::shared_ptr<FaserActsKalmanFilterAlg::TrackFitterFunction> -FaserActsKalmanFilterAlg::makeTrackFitterFunction(std::shared_ptr<const Acts::TrackingGeometry> trackingGeometry) { - - const std::string fieldMode = "FASER"; - const std::vector<double> constantFieldVector = {0., 0., 0.55}; - - Acts::Navigator::Config cfg{trackingGeometry}; - cfg.resolvePassive = false; - cfg.resolveMaterial = true; - cfg.resolveSensitive = true; - Acts::Navigator navigator( cfg ); - - std::unique_ptr<ActsExtrapolationDetail::VariantPropagator> varProp; - - if (fieldMode == "FASER") { - auto bField = std::make_shared<FASERMagneticFieldWrapper>(); - auto stepper = Acts::EigenStepper<>(std::move(bField)); - auto propagator = Acts::Propagator<decltype(stepper), - Acts::Navigator>(std::move(stepper), std::move(navigator)); - varProp = std::make_unique<ActsExtrapolationDetail::VariantPropagator>(propagator); - } else if (fieldMode == "Constant") { - Acts::Vector3 constantFieldVector = Acts::Vector3( - constantFieldVector[0], constantFieldVector[1], constantFieldVector[2]); - auto bField = std::make_shared<Acts::ConstantBField>(constantFieldVector); - auto stepper = Acts::EigenStepper<>(std::move(bField)); - auto propagator = Acts::Propagator<decltype(stepper), - Acts::Navigator>(std::move(stepper), std::move(navigator)); - varProp = std::make_unique<ActsExtrapolationDetail::VariantPropagator>(propagator); - } - - return boost::apply_visitor([&](const auto& propagator) -> TrackFitterFunction { - using Updater = Acts::GainMatrixUpdater; - using Smoother = Acts::GainMatrixSmoother; - using Fitter = Acts::KalmanFitter<typename std::decay_t<decltype(propagator)>, Updater, Smoother>; - Fitter fitter(std::move(propagator)); - return FitterFunctionImpl<Fitter>(std::move(fitter)); - }, *varProp); -} - */ -- GitLab