Skip to content
Snippets Groups Projects
Commit cd7224e6 authored by Xiaocong Ai's avatar Xiaocong Ai Committed by FASER Reco
Browse files

update FaserActsKalmanFilterAlg

parent cc363d92
No related branches found
No related tags found
1 merge request!418alma9-dev merge
Hide whitespace changes
Inline Side-by-side
Tracking/Acts/FaserActsKalmanFilter/src/FaserActsKalmanFilterAlg.cxx
+ 80
− 77
View file @ cd7224e6
......@@ -3,6 +3,7 @@
*/
#include "FaserActsKalmanFilterAlg.h"
#include "FaserActsGeometry/FASERMagneticFieldWrapper.h"
// ATHENA
#include "GaudiKernel/EventContext.h"
......@@ -62,14 +63,13 @@
#include <fstream>
#include <cmath>
using TrajectoryContainer = std::vector<FaserActsRecMultiTrajectory>;
using namespace Acts::UnitLiterals;
using Acts::VectorHelpers::eta;
using Acts::VectorHelpers::perp;
using Acts::VectorHelpers::phi;
using Acts::VectorHelpers::theta;
using ThisMeasurement = Acts::Measurement<IndexSourceLink, Acts::BoundIndices, 2>;
using ThisMeasurement = Acts::Measurement<Acts::BoundIndices, 2>;
using IdentifierMap = std::map<Identifier, Acts::GeometryIdentifier>;
FaserActsKalmanFilterAlg::FaserActsKalmanFilterAlg(const std::string& name, ISvcLocator* pSvcLocator) :
......@@ -80,12 +80,12 @@ StatusCode FaserActsKalmanFilterAlg::initialize() {
ATH_CHECK(m_fieldCondObjInputKey.initialize());
ATH_CHECK(m_trackingGeometryTool.retrieve());
ATH_CHECK(m_trackFinderTool.retrieve());
ATH_CHECK(m_trajectoryWriterTool.retrieve());
ATH_CHECK(m_trajectoryStatesWriterTool.retrieve());
//todo ATH_CHECK(m_trajectoryWriterTool.retrieve());
//todo ATH_CHECK(m_trajectoryStatesWriterTool.retrieve());
// ATH_CHECK(m_protoTrackWriterTool.retrieve());
ATH_CHECK(m_trackCollection.initialize());
ATH_CHECK(detStore()->retrieve(m_idHelper,"FaserSCT_ID"));
m_fit = makeTrackFitterFunction(m_trackingGeometryTool->trackingGeometry());
if (m_actsLogging == "VERBOSE") {
m_logger = Acts::getDefaultLogger("KalmanFitter", Acts::Logging::VERBOSE);
} else if (m_actsLogging == "DEBUG") {
......@@ -93,6 +93,12 @@ StatusCode FaserActsKalmanFilterAlg::initialize() {
} else {
m_logger = Acts::getDefaultLogger("KalmanFitter", Acts::Logging::INFO);
}
auto magneticField = std::make_shared<FASERMagneticFieldWrapper>();
double reverseFilteringMomThreshold = 0.1; //@todo: needs to be validated
//@todo: the multiple scattering and energy loss are originallly turned off
m_fit = makeTrackFitterFunction(m_trackingGeometryTool->trackingGeometry(), magneticField, false, false, reverseFilteringMomThreshold, Acts::FreeToBoundCorrection(false), *m_logger);
return StatusCode::SUCCESS;
}
......@@ -133,62 +139,58 @@ StatusCode FaserActsKalmanFilterAlg::execute() {
int n_trackSeeds = initialTrackParameters->size();
TrajectoryContainer trajectories;
trajectories.reserve(1);
for (int i = 0; i < n_trackSeeds; ++i) {
auto actsTrackContainer = std::make_shared<Acts::VectorTrackContainer>();
auto actsTrackStateContainer = std::make_shared<Acts::VectorMultiTrajectory>();
TrackContainer tracks(actsTrackContainer, actsTrackStateContainer);
Acts::KalmanFitterOptions<MeasurementCalibrator, Acts::VoidOutlierFinder, Acts::VoidReverseFilteringLogic> kfOptions(
geoctx,
magctx,
calctx,
MeasurementCalibrator(measurements->at(i)),
Acts::VoidOutlierFinder(),
Acts::VoidReverseFilteringLogic(),
Acts::LoggerWrapper{*m_logger},
Acts::PropagatorPlainOptions(),
&(*initialSurface)
);
kfOptions.multipleScattering = false;
kfOptions.energyLoss = false;
//@todo: the initialSurface should be targetSurface
FaserActsKalmanFilterAlg::GeneralFitterOptions options{
geoctx, magctx, calctx, &(*initialSurface),
Acts::PropagatorPlainOptions()};
for (int i = 0; i < n_trackSeeds; ++i) {
ATH_MSG_DEBUG("Invoke fitter");
auto result = (*m_fit)(sourceLinks->at(i), initialTrackParameters->at(i), kfOptions);
auto sls = sourceLinks->at(i);
std::vector<Acts::SourceLink> actsSls;
for(const auto& sl: sls){
actsSls.push_back(Acts::SourceLink{sl});
}
auto result = (*m_fit)(actsSls, initialTrackParameters->at(i), options, MeasurementCalibratorAdapter(MeasurementCalibrator(), measurements->at(i)), tracks);
int itrack = 0;
if (result.ok()) {
// Get the fit output object
const auto& fitOutput = result.value();
//const auto& fitOutput = result.value();
std::unique_ptr<Trk::Track> track = makeTrack(geoctx, result, clusters->at(i));
if (track) {
outputTracks->push_back(std::move(track));
}
// The track entry indices container. One element here.
std::vector<size_t> trackTips;
trackTips.reserve(1);
trackTips.emplace_back(fitOutput.lastMeasurementIndex);
// The fitted parameters container. One element (at most) here.
IndexedParams indexedParams;
if (fitOutput.fittedParameters) {
const auto& params = fitOutput.fittedParameters.value();
ATH_MSG_VERBOSE("Fitted paramemeters for track " << itrack);
ATH_MSG_VERBOSE(" parameters: " << params);
ATH_MSG_VERBOSE(" position: " << params.position(geoctx).transpose());
ATH_MSG_VERBOSE(" momentum: " << params.momentum().transpose());
// Push the fitted parameters to the container
indexedParams.emplace(fitOutput.lastMeasurementIndex, std::move(params));
} else {
ATH_MSG_DEBUG("No fitted paramemeters for track " << itrack);
}
// Create a SimMultiTrajectory
trajectories.emplace_back(std::move(fitOutput.fittedStates),
std::move(trackTips), std::move(indexedParams));
// // The track entry indices container. One element here.
// std::vector<size_t> trackTips;
// trackTips.reserve(1);
// trackTips.emplace_back(fitOutput.lastMeasurementIndex);
// // The fitted parameters container. One element (at most) here.
// IndexedParams indexedParams;
// if (fitOutput.fittedParameters) {
// const auto& params = fitOutput.fittedParameters.value();
// ATH_MSG_VERBOSE("Fitted paramemeters for track " << itrack);
// ATH_MSG_VERBOSE(" parameters: " << params);
// ATH_MSG_VERBOSE(" position: " << params.position(geoctx).transpose());
// ATH_MSG_VERBOSE(" momentum: " << params.momentum().transpose());
// // Push the fitted parameters to the container
// indexedParams.emplace(fitOutput.lastMeasurementIndex, std::move(params));
// } else {
// ATH_MSG_DEBUG("No fitted paramemeters for track " << itrack);
// }
// // Create a SimMultiTrajectory
// trajectories.emplace_back(std::move(fitOutput.fittedStates),
// std::move(trackTips), std::move(indexedParams));
} else {
ATH_MSG_WARNING("Fit failed for track " << itrack << " with error" << result.error());
// Fit failed, but still create a empty truth fit track
trajectories.push_back(FaserActsRecMultiTrajectory());
// trajectories.push_back(FaserActsRecMultiTrajectory());
}
}
......@@ -221,86 +223,88 @@ Acts::MagneticFieldContext FaserActsKalmanFilterAlg::getMagneticFieldContext(con
std::unique_ptr<Trk::Track>
FaserActsKalmanFilterAlg::makeTrack(Acts::GeometryContext& geoCtx, TrackFitterResult& fitResult, std::vector<const Tracker::FaserSCT_Cluster*> clusters) const {
using ConstTrackStateProxy =
Acts::detail_lt::TrackStateProxy<IndexSourceLink, 6, true>;
std::unique_ptr<Trk::Track> newtrack = nullptr;
//Get the fit output object
const auto& fitOutput = fitResult.value();
if (fitOutput.fittedParameters) {
DataVector<const Trk::TrackStateOnSurface>* finalTrajectory = new DataVector<const Trk::TrackStateOnSurface>{};
const auto& track = fitResult.value();
if (track.hasReferenceSurface()) {
std::unique_ptr<DataVector<const Trk::TrackStateOnSurface>> finalTrajectory = std::make_unique<DataVector<const Trk::TrackStateOnSurface>>();
std::vector<std::unique_ptr<const Acts::BoundTrackParameters>> actsSmoothedParam;
// Loop over all the output state to create track state
fitOutput.fittedStates.visitBackwards(fitOutput.lastMeasurementIndex, [&](const ConstTrackStateProxy& state) {
for (const auto& state : track.trackStatesReversed()) {
auto flag = state.typeFlags();
if (state.referenceSurface().associatedDetectorElement() != nullptr) {
// We need to determine the type of state
std::bitset<Trk::TrackStateOnSurface::NumberOfTrackStateOnSurfaceTypes> typePattern;
const Trk::TrackParameters *parm;
std::unique_ptr<Trk::TrackParameters> parm;
// State is a hole (no associated measurement), use predicted para meters
if (flag[Acts::TrackStateFlag::HoleFlag] == true) {
if (flag.test(Acts::TrackStateFlag::HoleFlag) == true) {
//@todo: ParticleHypothesis?
const Acts::BoundTrackParameters actsParam(state.referenceSurface().getSharedPtr(),
state.predicted(),
state.predictedCovariance());
parm = ConvertActsTrackParameterToATLAS(actsParam, geoCtx);
state.predictedCovariance(), Acts::ParticleHypothesis::pion());
parm = std::move(ConvertActsTrackParameterToATLAS(actsParam, geoCtx));
// auto boundaryCheck = m_boundaryCheckTool->boundaryCheck(*p arm);
typePattern.set(Trk::TrackStateOnSurface::Hole);
}
// The state was tagged as an outlier, use filtered parameters
else if (flag[Acts::TrackStateFlag::OutlierFlag] == true) {
else if (flag.test(Acts::TrackStateFlag::OutlierFlag) == true) {
//@todo: ParticleHypothesis?
const Acts::BoundTrackParameters actsParam(state.referenceSurface().getSharedPtr(),
state.filtered(), state.filteredCovariance());
parm = ConvertActsTrackParameterToATLAS(actsParam, geoCtx);
state.filtered(), state.filteredCovariance(), Acts::ParticleHypothesis::pion());
parm = std::move(ConvertActsTrackParameterToATLAS(actsParam, geoCtx));
typePattern.set(Trk::TrackStateOnSurface::Outlier);
}
// The state is a measurement state, use smoothed parameters
else {
//@todo: ParticleHypothesis?
const Acts::BoundTrackParameters actsParam(state.referenceSurface().getSharedPtr(),
state.smoothed(), state.smoothedCovariance());
state.smoothed(), state.smoothedCovariance(), Acts::ParticleHypothesis::pion());
actsSmoothedParam.push_back(std::make_unique<const Acts::BoundTrackParameters>(Acts::BoundTrackParameters(actsParam)));
// const auto& psurface=actsParam.referenceSurface();
Acts::Vector2 local(actsParam.parameters()[Acts::eBoundLoc0], actsParam.parameters()[Acts::eBoundLoc1]);
// const Acts::Vector3 dir = Acts::makeDirectionUnitFromPhiTheta(actsParam.parameters()[Acts::eBoundPhi], actsParam.parameters()[Acts::eBoundTheta]);
// auto pos=actsParam.position(tgContext);
parm = ConvertActsTrackParameterToATLAS(actsParam, geoCtx);
parm = std::move(ConvertActsTrackParameterToATLAS(actsParam, geoCtx));
typePattern.set(Trk::TrackStateOnSurface::Measurement);
}
Tracker::FaserSCT_ClusterOnTrack* measState = nullptr;
if (state.hasUncalibrated()) {
const Tracker::FaserSCT_Cluster* fitCluster = clusters.at(state.uncalibrated().index());
std::unique_ptr<Tracker::FaserSCT_ClusterOnTrack> measState = nullptr;
if (state.hasUncalibratedSourceLink()) {
auto sl = state.getUncalibratedSourceLink().template get<IndexSourceLink>();
const Tracker::FaserSCT_Cluster* fitCluster = clusters.at(sl.index());
if (fitCluster->detectorElement() != nullptr) {
measState = new Tracker::FaserSCT_ClusterOnTrack{
measState = std::make_unique<Tracker::FaserSCT_ClusterOnTrack>(
fitCluster,
Trk::LocalParameters{
Trk::DefinedParameter{fitCluster->localPosition()[0], Trk::loc1},
Trk::DefinedParameter{fitCluster->localPosition()[1], Trk::loc2}
},
fitCluster->localCovariance(),
Amg::MatrixX(fitCluster->localCovariance()),
m_idHelper->wafer_hash(fitCluster->detectorElement()->identify())
};
);
}
}
double nDoF = state.calibratedSize();
const Trk::FitQualityOnSurface *quality = new Trk::FitQualityOnSurface(state.chi2(), nDoF);
const Trk::TrackStateOnSurface *perState = new Trk::TrackStateOnSurface(measState, parm, quality, nullptr, typePattern);
const Trk::TrackStateOnSurface *perState = new Trk::TrackStateOnSurface(*quality, std::move(measState), std::move(parm), nullptr, typePattern);
// If a state was succesfully created add it to the trajectory
if (perState) {
finalTrajectory->insert(finalTrajectory->begin(), perState);
(*finalTrajectory).insert((*finalTrajectory).begin(), perState);
}
}
return;
});
} // state has referenceSurface
} //end loop for all states
// Create the track using the states
const Trk::TrackInfo newInfo(Trk::TrackInfo::TrackFitter::KalmanFitter, Trk::ParticleHypothesis::muon);
// Trk::FitQuality* q = nullptr;
// newInfo.setTrackFitter(Trk::TrackInfo::TrackFitter::KalmanFitter ); //Mark the fitter as KalmanFitter
newtrack = std::make_unique<Trk::Track>(newInfo, std::move(*finalTrajectory), nullptr);
}
newtrack = std::make_unique<Trk::Track>(newInfo, std::move(finalTrajectory), nullptr);
} // hasFittedParameters
return newtrack;
}
const Trk::TrackParameters*
std::unique_ptr<Trk::TrackParameters>
FaserActsKalmanFilterAlg ::ConvertActsTrackParameterToATLAS(const Acts::BoundTrackParameters &actsParameter, const Acts::GeometryContext& gctx) const {
using namespace Acts::UnitLiterals;
std::optional<AmgSymMatrix(5)> cov = std::nullopt;
......@@ -321,7 +325,6 @@ FaserActsKalmanFilterAlg ::ConvertActsTrackParameterToATLAS(const Acts::BoundTra
double tqOverP=actsParameter.get<Acts::eBoundQOverP>()*1_MeV;
double p = std::abs(1. / tqOverP);
Amg::Vector3D tmom(p * std::cos(tphi) * std::sin(ttheta), p * std::sin(tphi) * std::sin(ttheta), p * std::cos(ttheta));
const Trk::CurvilinearParameters * curv = new Trk::CurvilinearParameters(pos,tmom,tqOverP>0, cov);
return curv;
return std::make_unique<Trk::CurvilinearParameters>(pos, tmom, tqOverP>0, cov);
}
Tracking/Acts/FaserActsKalmanFilter/src/FaserActsKalmanFilterAlg.h
+ 46
− 20
View file @ cd7224e6
......@@ -24,32 +24,36 @@
#include "GeneratorObjects/McEventCollection.h"
#include "TrackerSimData/TrackerSimDataCollection.h"
#include "TrkTrack/TrackCollection.h"
#include "TrajectoryWriterTool.h"
//todo#include "TrajectoryWriterTool.h"
// ACTS
#include "Acts/MagneticField/ConstantBField.hpp"
#include "Acts/MagneticField/InterpolatedBFieldMap.hpp"
#include "Acts/MagneticField/SharedBField.hpp"
#include "Acts/Propagator/EigenStepper.hpp"
#include "Acts/Propagator/Propagator.hpp"
#include "Acts/Propagator/detail/SteppingLogger.hpp"
#include "Acts/TrackFitting/KalmanFitter.hpp"
#include "Acts/Geometry/TrackingGeometry.hpp"
#include "Acts/EventData/TrackParameters.hpp"
#include "Acts/EventData/TrackContainer.hpp"
#include "Acts/EventData/TrackProxy.hpp"
#include "Acts/EventData/VectorTrackContainer.hpp"
#include "Acts/EventData/VectorMultiTrajectory.hpp"
#include "Acts/Geometry/GeometryIdentifier.hpp"
#include "Acts/Utilities/Helpers.hpp"
#include "Acts/Definitions/Common.hpp"
#include "Acts/Utilities/Result.hpp"
// PACKAGE
#include "FaserActsGeometry/FASERMagneticFieldWrapper.h"
#include "FaserActsGeometryInterfaces/IFaserActsTrackingGeometryTool.h"
#include "FaserActsGeometryInterfaces/IFaserActsExtrapolationTool.h"
#include "FaserActsRecMultiTrajectory.h"
#include "FaserActsKalmanFilter/IndexSourceLink.h"
#include "FaserActsKalmanFilter/Measurement.h"
#include "FaserActsKalmanFilter/ITrackFinderTool.h"
//#include "ProtoTrackWriterTool.h"
#include "RootTrajectoryStatesWriterTool.h"
//todo #include "RootTrajectoryStatesWriterTool.h"
// STL
#include <memory>
......@@ -70,8 +74,9 @@ namespace TrackerDD {
class SCT_DetectorManager;
}
using TrajectoryContainer = std::vector<FaserActsRecMultiTrajectory>;
using BField_t = FASERMagneticFieldWrapper;
//@toberemoved
//using TrajectoryContainer = std::vector<FaserActsRecMultiTrajectory>;
//using BField_t = FASERMagneticFieldWrapper;
//class FaserActsKalmanFilterAlg : public AthReentrantAlgorithm {
class FaserActsKalmanFilterAlg : public AthAlgorithm {
......@@ -83,26 +88,47 @@ public:
// StatusCode execute(const EventContext& ctx) const override;
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>;
using TrackContainer =
Acts::TrackContainer<Acts::VectorTrackContainer,
Acts::VectorMultiTrajectory, std::shared_ptr>;
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 IndexedParams = std::unordered_map<size_t, Acts::BoundTrackParameters>;
using TrackFitterOptions =
Acts::KalmanFitterOptions<MeasurementCalibrator, Acts::VoidOutlierFinder,
Acts::VoidReverseFilteringLogic>;
using TrackFitterResult =
Acts::Result<Acts::KalmanFitterResult<IndexSourceLink>>;
using TrackFitterResult = Acts::Result<TrackContainer::TrackProxy>;
using TrackParameters = Acts::CurvilinearTrackParameters;
using TrackParameters = Acts::BoundTrackParameters;
class TrackFitterFunction {
public:
virtual ~TrackFitterFunction() = default;
virtual TrackFitterResult operator()(const std::vector<IndexSourceLink>&,
const TrackParameters&,
const TrackFitterOptions&) const = 0;
virtual TrackFitterResult operator()(
const std::vector<Acts::SourceLink> &sourceLinks,
const TrackParameters &initialParameters,
const GeneralFitterOptions& options,
const MeasurementCalibratorAdapter& calibrator,
TrackContainer& tracks
) const = 0;
};
static std::shared_ptr<TrackFitterFunction> makeTrackFitterFunction(
std::shared_ptr<const Acts::TrackingGeometry> trackingGeometry);
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;
......@@ -113,15 +139,15 @@ private:
Gaudi::Property<std::string> m_actsLogging {this, "ActsLogging", "VERBOSE"};
std::unique_ptr<Trk::Track> makeTrack(Acts::GeometryContext& tgContext, TrackFitterResult& fitResult, std::vector<const Tracker::FaserSCT_Cluster*> clusters) const;
const Trk::TrackParameters* ConvertActsTrackParameterToATLAS(const Acts::BoundTrackParameters &actsParameter, const Acts::GeometryContext& gctx) const;
std::unique_ptr<Trk::TrackParameters> ConvertActsTrackParameterToATLAS(const Acts::BoundTrackParameters &actsParameter, const Acts::GeometryContext& gctx) const;
// Read handle for conditions object to get the field cache
SG::ReadCondHandleKey<FaserFieldCacheCondObj> m_fieldCondObjInputKey {this, "FaserFieldCacheCondObj", "fieldCondObj", "Name of the Magnetic Field conditions object key"};
ToolHandle<ITrackFinderTool> m_trackFinderTool {this, "TrackFinderTool", "TruthTrackFinderTool"};
ToolHandle<IFaserActsTrackingGeometryTool> m_trackingGeometryTool {this, "TrackingGeometryTool", "FaserActsTrackingGeometryTool"};
ToolHandle<TrajectoryWriterTool> m_trajectoryWriterTool {this, "OutputTool", "TrajectoryWriterTool"};
ToolHandle<RootTrajectoryStatesWriterTool> m_trajectoryStatesWriterTool {this, "RootTrajectoryStatesWriterTool", "RootTrajectoryStatesWriterTool"};
//todo ToolHandle<TrajectoryWriterTool> m_trajectoryWriterTool {this, "OutputTool", "TrajectoryWriterTool"};
//todo ToolHandle<RootTrajectoryStatesWriterTool> m_trajectoryStatesWriterTool {this, "RootTrajectoryStatesWriterTool", "RootTrajectoryStatesWriterTool"};
// ToolHandle<ProtoTrackWriterTool> m_protoTrackWriterTool {this, "ProtoTrackWriterTool", "ProtoTrackWriterTool"};
SG::ReadHandleKey<McEventCollection> m_mcEventKey { this, "McEventCollection", "BeamTruthEvent" };
......
Tracking/Acts/FaserActsKalmanFilter/src/TrackFittingFunction.cxx
+ 107
− 19
View file @ cd7224e6
#include "FaserActsKalmanFilterAlg.h"
#include "FaserActsGeometry/FASERMagneticFieldWrapper.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 Updater = Acts::GainMatrixUpdater;
using Smoother = Acts::GainMatrixSmoother;
using Stepper = Acts::EigenStepper<>;
using Propagator = Acts::Propagator<Stepper, Acts::Navigator>;
using Fitter = Acts::KalmanFitter<Propagator, Updater, Smoother>;
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;
TrackFitterFunctionImpl(Fitter &&f) : trackFitter(std::move(f)) {}
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<IndexSourceLink> &sourceLinks,
const std::vector<Acts::SourceLink> &sourceLinks,
const FaserActsKalmanFilterAlg::TrackParameters &initialParameters,
const FaserActsKalmanFilterAlg::TrackFitterOptions &options)
const override {
return trackFitter.fit(sourceLinks, initialParameters, options);
};
const FaserActsKalmanFilterAlg::GeneralFitterOptions& options,
const MeasurementCalibratorAdapter& calibrator,
FaserActsKalmanFilterAlg::TrackContainer& tracks) const override {
const auto kfOptions = makeKfOptions(options, calibrator);
return trackFitter.fit(sourceLinks.begin(), sourceLinks.end(), initialParameters,
kfOptions, tracks);
}
};
} // namespace
......@@ -36,19 +104,39 @@ struct TrackFitterFunctionImpl
std::shared_ptr<FaserActsKalmanFilterAlg::TrackFitterFunction>
FaserActsKalmanFilterAlg::makeTrackFitterFunction(
std::shared_ptr<const Acts::TrackingGeometry> trackingGeometry) {
auto magneticField = std::make_shared<FASERMagneticFieldWrapper>();
auto stepper = Stepper(std::move(magneticField));
Acts::Navigator::Config cfg{trackingGeometry};
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);
Propagator propagator(std::move(stepper), std::move(navigator));
Fitter trackFitter(std::move(propagator));
return std::make_shared<TrackFitterFunctionImpl>(std::move(trackFitter));
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 {
......
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment