From 6cc821513f6c766986879248b9835a95a32b87e6 Mon Sep 17 00:00:00 2001 From: andiwand Date: Wed, 22 Nov 2023 19:50:27 +0100 Subject: [PATCH 1/4] hacky two way ckf --- .../Algorithms/TrackFinding/CMakeLists.txt | 1 + .../TrackFinding/MyTrackFindingAlgorithm.hpp | 104 +++++++ .../src/MyTrackFindingAlgorithm.cpp | 258 ++++++++++++++++++ Examples/Python/src/TrackFinding.cpp | 25 ++ Examples/Scripts/Python/full_chain_odd.py | 82 +++--- 5 files changed, 439 insertions(+), 31 deletions(-) create mode 100644 Examples/Algorithms/TrackFinding/include/ActsExamples/TrackFinding/MyTrackFindingAlgorithm.hpp create mode 100644 Examples/Algorithms/TrackFinding/src/MyTrackFindingAlgorithm.cpp diff --git a/Examples/Algorithms/TrackFinding/CMakeLists.txt b/Examples/Algorithms/TrackFinding/CMakeLists.txt index ab7b025c8a2..b68372f556d 100644 --- a/Examples/Algorithms/TrackFinding/CMakeLists.txt +++ b/Examples/Algorithms/TrackFinding/CMakeLists.txt @@ -8,6 +8,7 @@ add_library( src/HoughTransformSeeder.cpp src/TrackParamsEstimationAlgorithm.cpp src/SeedingFTFAlgorithm.cpp + src/MyTrackFindingAlgorithm.cpp ) target_include_directories( diff --git a/Examples/Algorithms/TrackFinding/include/ActsExamples/TrackFinding/MyTrackFindingAlgorithm.hpp b/Examples/Algorithms/TrackFinding/include/ActsExamples/TrackFinding/MyTrackFindingAlgorithm.hpp new file mode 100644 index 00000000000..bc1061da20f --- /dev/null +++ b/Examples/Algorithms/TrackFinding/include/ActsExamples/TrackFinding/MyTrackFindingAlgorithm.hpp @@ -0,0 +1,104 @@ +// 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 + +#include "Acts/EventData/MultiTrajectory.hpp" +#include "Acts/EventData/SourceLink.hpp" +#include "Acts/EventData/TrackContainer.hpp" +#include "Acts/EventData/TrackProxy.hpp" +#include "Acts/EventData/VectorMultiTrajectory.hpp" +#include "Acts/EventData/VectorTrackContainer.hpp" +#include "Acts/Geometry/TrackingGeometry.hpp" +#include "Acts/TrackFinding/CombinatorialKalmanFilter.hpp" +#include "Acts/TrackFinding/MeasurementSelector.hpp" +#include "Acts/TrackFinding/SourceLinkAccessorConcept.hpp" +#include "Acts/TrackFinding/TrackSelector.hpp" +#include "Acts/Utilities/Logger.hpp" +#include "Acts/Utilities/Result.hpp" +#include "ActsExamples/EventData/IndexSourceLink.hpp" +#include "ActsExamples/EventData/Measurement.hpp" +#include "ActsExamples/EventData/Track.hpp" +#include "ActsExamples/Framework/DataHandle.hpp" +#include "ActsExamples/Framework/IAlgorithm.hpp" +#include "ActsExamples/Framework/ProcessCode.hpp" +#include "ActsExamples/MagneticField/MagneticField.hpp" + +#include +#include +#include +#include +#include +#include +#include + +namespace Acts { +class MagneticFieldProvider; +class TrackingGeometry; +} // namespace Acts + +namespace ActsExamples { +struct AlgorithmContext; + +class MyTrackFindingAlgorithm final : public IAlgorithm { + public: + struct Config { + /// Input measurements collection. + std::string inputMeasurements; + /// Input source links collection. + std::string inputSourceLinks; + /// Input initial track parameter estimates for for each proto track. + std::string inputInitialTrackParameters; + /// Output find trajectories collection. + std::string outputTracks; + + /// The tracking geometry that should be used. + std::shared_ptr trackingGeometry; + /// The magnetic field that should be used. + std::shared_ptr magneticField; + + /// CKF measurement selector config + Acts::MeasurementSelector::Config measurementSelectorCfg; + /// Track selector config + std::optional trackSelectorCfg; + /// Maximum number of propagation steps + unsigned int maxSteps = 100000; + }; + + /// Constructor of the track finding algorithm + /// + /// @param config is the config struct to configure the algorithm + /// @param level is the logging level + MyTrackFindingAlgorithm(Config config, Acts::Logging::Level level); + + /// Framework execute method of the track finding algorithm + /// + /// @param ctx is the algorithm context that holds event-wise information + /// @return a process code to steer the algorithm flow + ActsExamples::ProcessCode execute( + const ActsExamples::AlgorithmContext& ctx) const final; + + /// Get readonly access to the config parameters + const Config& config() const { return m_cfg; } + + private: + Config m_cfg; + std::optional m_trackSelector; + + ReadDataHandle m_inputMeasurements{this, + "InputMeasurements"}; + ReadDataHandle m_inputSourceLinks{ + this, "InputSourceLinks"}; + + ReadDataHandle m_inputInitialTrackParameters{ + this, "InputInitialTrackParameters"}; + + WriteDataHandle m_outputTracks{this, "OutputTracks"}; +}; + +} // namespace ActsExamples diff --git a/Examples/Algorithms/TrackFinding/src/MyTrackFindingAlgorithm.cpp b/Examples/Algorithms/TrackFinding/src/MyTrackFindingAlgorithm.cpp new file mode 100644 index 00000000000..24b72fcd6be --- /dev/null +++ b/Examples/Algorithms/TrackFinding/src/MyTrackFindingAlgorithm.cpp @@ -0,0 +1,258 @@ +// 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 "ActsExamples/TrackFinding/MyTrackFindingAlgorithm.hpp" + +#include "Acts/Definitions/Algebra.hpp" +#include "Acts/Definitions/Direction.hpp" +#include "Acts/EventData/MultiTrajectory.hpp" +#include "Acts/EventData/TrackContainer.hpp" +#include "Acts/EventData/TrackParameters.hpp" +#include "Acts/EventData/TrackStateType.hpp" +#include "Acts/EventData/VectorMultiTrajectory.hpp" +#include "Acts/EventData/VectorTrackContainer.hpp" +#include "Acts/Propagator/EigenStepper.hpp" +#include "Acts/Propagator/Propagator.hpp" +#include "Acts/Surfaces/PerigeeSurface.hpp" +#include "Acts/Surfaces/Surface.hpp" +#include "Acts/TrackFinding/CombinatorialKalmanFilter.hpp" +#include "Acts/TrackFitting/GainMatrixSmoother.hpp" +#include "Acts/TrackFitting/GainMatrixUpdater.hpp" +#include "Acts/TrackFitting/KalmanFitter.hpp" +#include "Acts/Utilities/Delegate.hpp" +#include "ActsExamples/EventData/Measurement.hpp" +#include "ActsExamples/EventData/MeasurementCalibration.hpp" +#include "ActsExamples/EventData/Track.hpp" +#include "ActsExamples/Framework/AlgorithmContext.hpp" +#include "ActsExamples/Framework/ProcessCode.hpp" + +#include +#include +#include +#include +#include +#include +#include +#include + +ActsExamples::MyTrackFindingAlgorithm::MyTrackFindingAlgorithm( + Config config, Acts::Logging::Level level) + : ActsExamples::IAlgorithm("MyTrackFindingAlgorithm", level), + m_cfg(std::move(config)), + m_trackSelector([this]() -> std::optional { + if (m_cfg.trackSelectorCfg.has_value()) { + return {m_cfg.trackSelectorCfg.value()}; + } else { + return std::nullopt; + } + }()) { + if (m_cfg.inputMeasurements.empty()) { + throw std::invalid_argument("Missing measurements input collection"); + } + if (m_cfg.inputSourceLinks.empty()) { + throw std::invalid_argument("Missing source links input collection"); + } + if (m_cfg.inputInitialTrackParameters.empty()) { + throw std::invalid_argument( + "Missing initial track parameters input collection"); + } + if (m_cfg.outputTracks.empty()) { + throw std::invalid_argument("Missing tracks output collection"); + } + + m_inputMeasurements.initialize(m_cfg.inputMeasurements); + m_inputSourceLinks.initialize(m_cfg.inputSourceLinks); + m_inputInitialTrackParameters.initialize(m_cfg.inputInitialTrackParameters); + m_outputTracks.initialize(m_cfg.outputTracks); +} + +ActsExamples::ProcessCode ActsExamples::MyTrackFindingAlgorithm::execute( + const ActsExamples::AlgorithmContext& ctx) const { + // Read input data + const auto& measurements = m_inputMeasurements(ctx); + const auto& sourceLinks = m_inputSourceLinks(ctx); + const auto& initialParameters = m_inputInitialTrackParameters(ctx); + + // Construct a perigee surface as the target surface + auto pSurface = Acts::Surface::makeShared( + Acts::Vector3{0., 0., 0.}); + + PassThroughCalibrator pcalibrator; + MeasurementCalibratorAdapter calibrator(pcalibrator, measurements); + Acts::GainMatrixUpdater kfUpdater; + Acts::GainMatrixSmoother kfSmoother; + Acts::MeasurementSelector measSel{m_cfg.measurementSelectorCfg}; + + Acts::CombinatorialKalmanFilterExtensions + extensions; + extensions.calibrator.connect<&MeasurementCalibratorAdapter::calibrate>( + &calibrator); + extensions.updater.connect< + &Acts::GainMatrixUpdater::operator()>( + &kfUpdater); + extensions.smoother.connect< + &Acts::GainMatrixSmoother::operator()>( + &kfSmoother); + extensions.measurementSelector + .connect<&Acts::MeasurementSelector::select>( + &measSel); + + IndexSourceLinkAccessor slAccessor; + slAccessor.container = &sourceLinks; + Acts::SourceLinkAccessorDelegate + slAccessorDelegate; + slAccessorDelegate.connect<&IndexSourceLinkAccessor::range>(&slAccessor); + + using Propagator = Acts::Propagator, Acts::Navigator>; + using CKF = + Acts::CombinatorialKalmanFilter; + using Options = + Acts::CombinatorialKalmanFilterOptions; + + Acts::EigenStepper<> stepper(m_cfg.magneticField); + Acts::Navigator navigator({m_cfg.trackingGeometry}, + logger().cloneWithSuffix("Navigator")); + Propagator propagator(std::move(stepper), std::move(navigator), + logger().cloneWithSuffix("Propagator")); + CKF ckf(std::move(propagator), logger().cloneWithSuffix("Finder")); + + Acts::PropagatorPlainOptions fwdPropOptions; + fwdPropOptions.maxSteps = m_cfg.maxSteps; + fwdPropOptions.direction = Acts::Direction::Forward; + + Options fwdOptions(ctx.geoContext, ctx.magFieldContext, ctx.calibContext, + slAccessorDelegate, extensions, fwdPropOptions, + pSurface.get()); + fwdOptions.filterTargetSurface = nullptr; + fwdOptions.smoothing = true; + fwdOptions.smoothingTargetSurface = pSurface.get(); + fwdOptions.smoothingTargetSurfaceStrategy = + Acts::CombinatorialKalmanFilterTargetSurfaceStrategy::first; + + Acts::PropagatorPlainOptions bwdPropOptions; + bwdPropOptions.maxSteps = m_cfg.maxSteps; + bwdPropOptions.direction = Acts::Direction::Backward; + + Options bwdOptions(ctx.geoContext, ctx.magFieldContext, ctx.calibContext, + slAccessorDelegate, extensions, bwdPropOptions, + pSurface.get()); + bwdOptions.filterTargetSurface = pSurface.get(); + bwdOptions.smoothing = true; // TODO this should not be necessary + bwdOptions.smoothingTargetSurface = pSurface.get(); + bwdOptions.smoothingTargetSurfaceStrategy = + Acts::CombinatorialKalmanFilterTargetSurfaceStrategy::last; + + // Perform the track finding for all initial parameters + ACTS_DEBUG("Invoke track finding with " << initialParameters.size() + << " seeds."); + + auto trackContainer = std::make_shared(); + auto trackStateContainer = std::make_shared(); + + auto trackContainerTemp = std::make_shared(); + auto trackStateContainerTemp = + std::make_shared(); + + TrackContainer tracks(trackContainer, trackStateContainer); + TrackContainer tracksTemp(trackContainerTemp, trackStateContainerTemp); + + tracks.addColumn("trackGroup"); + tracksTemp.addColumn("trackGroup"); + Acts::TrackAccessor seedNumber("trackGroup"); + + unsigned int nSeed = 0; + + for (std::size_t iseed = 0; iseed < initialParameters.size(); ++iseed) { + // Clear trackContainerTemp and trackStateContainerTemp + tracksTemp.clear(); + + auto fwdResult = + ckf.findTracks(initialParameters.at(iseed), fwdOptions, tracksTemp); + nSeed++; + + if (!fwdResult.ok()) { + ACTS_WARNING("Forward track finding failed for seed " + << iseed << " with error" << fwdResult.error()); + continue; + } + + auto& fwdTracksForSeed = fwdResult.value(); + for (auto& fwdTrack : fwdTracksForSeed) { + std::optional firstState; + for (auto st : fwdTrack.trackStatesReversed()) { + bool isMeasurement = + st.typeFlags().test(Acts::TrackStateFlag::MeasurementFlag); + bool isOutlier = st.typeFlags().test(Acts::TrackStateFlag::OutlierFlag); + // We are excluding non measurement states and outlier here. Those can + // decrease resolution because only the smoothing corrected the very + // first prediction as filtering is not possible. + if (isMeasurement && !isOutlier) { + firstState = st; + } + } + + Acts::BoundTrackParameters bwdInitialParameters( + firstState->referenceSurface().getSharedPtr(), + firstState->parameters(), firstState->covariance(), + initialParameters.at(iseed).particleHypothesis()); + + auto bwdResult = + ckf.findTracks(bwdInitialParameters, bwdOptions, tracksTemp); + + if (!bwdResult.ok()) { + ACTS_WARNING("Backward track finding failed for seed " + << iseed << " with error" << bwdResult.error()); + continue; + } + + auto& bwdTracks = bwdResult.value(); + for (auto& bwdTrack : bwdTracks) { + if (bwdTrack.nTrackStates() < 2) { + continue; + } + + bwdTrack.reverseTrackStates(true); + seedNumber(bwdTrack) = nSeed; + + auto innermostFwdState = fwdTrack.trackStatesReversed().begin(); + for (auto i = innermostFwdState; + ++i != fwdTrack.trackStatesReversed().end(); + i = ++innermostFwdState) { + } + + (*innermostFwdState).previous() = + (*std::next(bwdTrack.trackStatesReversed().begin())).index(); + bwdTrack.tipIndex() = fwdTrack.tipIndex(); + + Acts::calculateTrackQuantities(bwdTrack); + + if (!m_trackSelector.has_value() || + m_trackSelector->isValidTrack(bwdTrack)) { + auto destProxy = tracks.getTrack(tracks.addTrack()); + destProxy.copyFrom(bwdTrack, true); + } + } + } + } + + ACTS_DEBUG("Track finding with " << tracks.size() << " track candidates."); + + auto constTrackStateContainer = + std::make_shared( + std::move(*trackStateContainer)); + + auto constTrackContainer = std::make_shared( + std::move(*trackContainer)); + + ConstTrackContainer constTracks{constTrackContainer, + constTrackStateContainer}; + + m_outputTracks(ctx, std::move(constTracks)); + return ActsExamples::ProcessCode::SUCCESS; +} diff --git a/Examples/Python/src/TrackFinding.cpp b/Examples/Python/src/TrackFinding.cpp index 43968c37fa2..d6fb7793e85 100644 --- a/Examples/Python/src/TrackFinding.cpp +++ b/Examples/Python/src/TrackFinding.cpp @@ -20,6 +20,7 @@ #include "Acts/Utilities/TypeTraits.hpp" #include "ActsExamples/EventData/Track.hpp" #include "ActsExamples/TrackFinding/HoughTransformSeeder.hpp" +#include "ActsExamples/TrackFinding/MyTrackFindingAlgorithm.hpp" #include "ActsExamples/TrackFinding/SeedingAlgorithm.hpp" #include "ActsExamples/TrackFinding/SeedingFTFAlgorithm.hpp" #include "ActsExamples/TrackFinding/SeedingOrthogonalAlgorithm.hpp" @@ -331,6 +332,30 @@ void addTrackFinding(Context& ctx) { ACTS_PYTHON_STRUCT_END(); } + { + using Alg = ActsExamples::MyTrackFindingAlgorithm; + using Config = Alg::Config; + + auto alg = py::class_>( + mex, "MyTrackFindingAlgorithm") + .def(py::init(), + py::arg("config"), py::arg("level")) + .def_property_readonly("config", &Alg::config); + + auto c = py::class_(alg, "Config").def(py::init<>()); + ACTS_PYTHON_STRUCT_BEGIN(c, Config); + ACTS_PYTHON_MEMBER(inputMeasurements); + ACTS_PYTHON_MEMBER(inputSourceLinks); + ACTS_PYTHON_MEMBER(inputInitialTrackParameters); + ACTS_PYTHON_MEMBER(outputTracks); + ACTS_PYTHON_MEMBER(magneticField); + ACTS_PYTHON_MEMBER(trackingGeometry); + ACTS_PYTHON_MEMBER(measurementSelectorCfg); + ACTS_PYTHON_MEMBER(trackSelectorCfg); + ACTS_PYTHON_MEMBER(maxSteps); + ACTS_PYTHON_STRUCT_END(); + } + ACTS_PYTHON_DECLARE_ALGORITHM(ActsExamples::TrajectoriesToPrototracks, mex, "TrajectoriesToPrototracks", inputTrajectories, outputProtoTracks); diff --git a/Examples/Scripts/Python/full_chain_odd.py b/Examples/Scripts/Python/full_chain_odd.py index 26c951249eb..4d2e9d28428 100755 --- a/Examples/Scripts/Python/full_chain_odd.py +++ b/Examples/Scripts/Python/full_chain_odd.py @@ -77,12 +77,12 @@ MomentumConfig(1.0 * u.GeV, 10.0 * u.GeV, transverse=True), EtaConfig(-3.0, 3.0), PhiConfig(0.0, 360.0 * u.degree), - ParticleConfig(4, acts.PdgParticle.eMuon, randomizeCharge=True), + ParticleConfig(1, acts.PdgParticle.eMuon, randomizeCharge=True), vtxGen=acts.examples.GaussianVertexGenerator( mean=acts.Vector4(0, 0, 0, 0), stddev=acts.Vector4(0.0125 * u.mm, 0.0125 * u.mm, 55.5 * u.mm, 1.0 * u.ns), ), - multiplicity=200, + multiplicity=1, rnd=rnd, ) else: @@ -164,6 +164,55 @@ outputDirRoot=outputDir, ) +s.addAlgorithm( + acts.examples.MyTrackFindingAlgorithm( + level=acts.logging.INFO, + + magneticField=field, + trackingGeometry=trackingGeometry, + + inputMeasurements="measurements", + inputSourceLinks="sourcelinks", + inputInitialTrackParameters="estimatedparameters", + outputTracks="mytracks", + + measurementSelectorCfg=acts.MeasurementSelector.Config( + [ + ( + acts.GeometryIdentifier(), + ( + [], + [30], + [10], + ), + ) + ] + ), + ) +) + +trackStatesWriter = acts.examples.RootTrackStatesWriter( + level=acts.logging.INFO, + inputTracks="mytracks", + inputParticles="particles_selected", + inputSimHits="simhits", + inputMeasurementParticlesMap="measurement_particles_map", + inputMeasurementSimHitsMap="measurement_simhits_map", + filePath=str(outputDir / f"trackstates_hi.root"), + treeName="trackstates", +) +s.addWriter(trackStatesWriter) + +trackSummaryWriter = acts.examples.RootTrackSummaryWriter( + level=acts.logging.INFO, + inputTracks="mytracks", + inputParticles="particles_selected", + inputMeasurementParticlesMap="measurement_particles_map", + filePath=str(outputDir / f"tracksummary_hi.root"), + treeName="tracksummary", +) +s.addWriter(trackSummaryWriter) + addCKFTracks( s, trackingGeometry, @@ -179,33 +228,4 @@ # outputDirCsv=outputDir, ) -if ambiguity_MLSolver: - addAmbiguityResolutionML( - s, - AmbiguityResolutionMLConfig( - maximumSharedHits=3, maximumIterations=1000000, nMeasurementsMin=7 - ), - outputDirRoot=outputDir, - # outputDirCsv=outputDir, - onnxModelFile=os.path.dirname(__file__) - + "/MLAmbiguityResolution/duplicateClassifier.onnx", - ) -else: - addAmbiguityResolution( - s, - AmbiguityResolutionConfig( - maximumSharedHits=3, maximumIterations=1000000, nMeasurementsMin=7 - ), - outputDirRoot=outputDir, - writeCovMat=True, - # outputDirCsv=outputDir, - ) - -addVertexFitting( - s, - field, - vertexFinder=VertexFinder.Iterative, - outputDirRoot=outputDir, -) - s.run() From d9a3ae7bdcb0994903f2e303dcfee63447ef82e8 Mon Sep 17 00:00:00 2001 From: andiwand Date: Thu, 23 Nov 2023 13:17:06 +0100 Subject: [PATCH 2/4] port back to normal track finding algo --- .../TrackFinding/TrackFindingAlgorithm.hpp | 11 +- .../src/TrackFindingAlgorithm.cpp | 84 +++++-- Examples/Python/src/TrackFinding.cpp | 2 +- Examples/Scripts/Python/full_chain_odd.py | 82 +++---- Examples/Scripts/Python/full_chain_odd_hi.py | 231 ++++++++++++++++++ 5 files changed, 337 insertions(+), 73 deletions(-) create mode 100755 Examples/Scripts/Python/full_chain_odd_hi.py diff --git a/Examples/Algorithms/TrackFinding/include/ActsExamples/TrackFinding/TrackFindingAlgorithm.hpp b/Examples/Algorithms/TrackFinding/include/ActsExamples/TrackFinding/TrackFindingAlgorithm.hpp index c15199d1aff..cb40c1cb6d6 100644 --- a/Examples/Algorithms/TrackFinding/include/ActsExamples/TrackFinding/TrackFindingAlgorithm.hpp +++ b/Examples/Algorithms/TrackFinding/include/ActsExamples/TrackFinding/TrackFindingAlgorithm.hpp @@ -91,14 +91,17 @@ class TrackFindingAlgorithm final : public IAlgorithm { std::shared_ptr findTracks; /// CKF measurement selector config Acts::MeasurementSelector::Config measurementSelectorCfg; - /// Compute shared hit information - bool computeSharedHits = false; /// Track selector config std::optional trackSelectorCfg = std::nullopt; - /// Run backward finding - bool backward = false; + + /// Run finding in two directions + bool twoWay = true; + /// Maximum number of propagation steps unsigned int maxSteps = 100000; + + /// Compute shared hit information + bool computeSharedHits = false; }; /// Constructor of the track finding algorithm diff --git a/Examples/Algorithms/TrackFinding/src/TrackFindingAlgorithm.cpp b/Examples/Algorithms/TrackFinding/src/TrackFindingAlgorithm.cpp index bebf8711e20..f0aa3470cd2 100644 --- a/Examples/Algorithms/TrackFinding/src/TrackFindingAlgorithm.cpp +++ b/Examples/Algorithms/TrackFinding/src/TrackFindingAlgorithm.cpp @@ -81,11 +81,6 @@ ActsExamples::ProcessCode ActsExamples::TrackFindingAlgorithm::execute( auto pSurface = Acts::Surface::makeShared( Acts::Vector3{0., 0., 0.}); - Acts::PropagatorPlainOptions pOptions; - pOptions.maxSteps = m_cfg.maxSteps; - pOptions.direction = - m_cfg.backward ? Acts::Direction::Backward : Acts::Direction::Forward; - PassThroughCalibrator pcalibrator; MeasurementCalibratorAdapter calibrator(pcalibrator, measurements); Acts::GainMatrixUpdater kfUpdater; @@ -112,13 +107,30 @@ ActsExamples::ProcessCode ActsExamples::TrackFindingAlgorithm::execute( slAccessorDelegate; slAccessorDelegate.connect<&IndexSourceLinkAccessor::range>(&slAccessor); + Acts::PropagatorPlainOptions firstPropOptions; + firstPropOptions.maxSteps = m_cfg.maxSteps; + firstPropOptions.direction = Acts::Direction::Forward; + + Acts::PropagatorPlainOptions secondPropOptions; + secondPropOptions.maxSteps = m_cfg.maxSteps; + secondPropOptions.direction = firstPropOptions.direction.invert(); + // Set the CombinatorialKalmanFilter options - ActsExamples::TrackFindingAlgorithm::TrackFinderOptions options( + ActsExamples::TrackFindingAlgorithm::TrackFinderOptions firstOptions( ctx.geoContext, ctx.magFieldContext, ctx.calibContext, slAccessorDelegate, - extensions, pOptions, pSurface.get()); - options.smoothingTargetSurfaceStrategy = + extensions, firstPropOptions, pSurface.get()); + firstOptions.smoothing = true; + firstOptions.smoothingTargetSurfaceStrategy = Acts::CombinatorialKalmanFilterTargetSurfaceStrategy::first; + ActsExamples::TrackFindingAlgorithm::TrackFinderOptions secondOptions( + ctx.geoContext, ctx.magFieldContext, ctx.calibContext, slAccessorDelegate, + extensions, firstPropOptions, pSurface.get()); + secondOptions.filterTargetSurface = pSurface.get(); + secondOptions.smoothing = true; + secondOptions.smoothingTargetSurfaceStrategy = + Acts::CombinatorialKalmanFilterTargetSurfaceStrategy::last; + // Perform the track finding for all initial parameters ACTS_DEBUG("Invoke track finding with " << initialParameters.size() << " seeds."); @@ -143,25 +155,63 @@ ActsExamples::ProcessCode ActsExamples::TrackFindingAlgorithm::execute( // Clear trackContainerTemp and trackStateContainerTemp tracksTemp.clear(); - auto result = - (*m_cfg.findTracks)(initialParameters.at(iseed), options, tracksTemp); + auto firstResult = (*m_cfg.findTracks)(initialParameters.at(iseed), + firstOptions, tracksTemp); m_nTotalSeeds++; nSeed++; - if (!result.ok()) { + if (!firstResult.ok()) { m_nFailedSeeds++; ACTS_WARNING("Track finding failed for seed " << iseed << " with error" - << result.error()); + << firstResult.error()); continue; } - auto& tracksForSeed = result.value(); - for (auto& track : tracksForSeed) { - seedNumber(track) = nSeed; + auto& firstTracksForSeed = firstResult.value(); + for (auto& firstTrack : firstTracksForSeed) { + seedNumber(firstTrack) = nSeed; + if (!m_trackSelector.has_value() || - m_trackSelector->isValidTrack(track)) { + m_trackSelector->isValidTrack(firstTrack)) { auto destProxy = tracks.getTrack(tracks.addTrack()); - destProxy.copyFrom(track, true); // make sure we copy track states! + destProxy.copyFrom(firstTrack, true); + } + + if (m_cfg.twoWay) { + auto secondResult = (*m_cfg.findTracks)(initialParameters.at(iseed), + secondOptions, tracksTemp); + + if (!secondResult.ok()) { + ACTS_WARNING("Second track finding failed for seed " + << iseed << " with error" << secondResult.error()); + continue; + } + + auto firstFirstState = + std::next(firstTrack.trackStatesReversed().begin(), + firstTrack.nTrackStates() - 1); + + auto& secondTracksForSeed = secondResult.value(); + for (auto& secondTrack : secondTracksForSeed) { + if (secondTrack.nTrackStates() < 2) { + continue; + } + + secondTrack.reverseTrackStates(true); + seedNumber(secondTrack) = nSeed; + + (*firstFirstState).previous() = + (*std::next(secondTrack.trackStatesReversed().begin())).index(); + secondTrack.tipIndex() = firstTrack.tipIndex(); + + Acts::calculateTrackQuantities(secondTrack); + + if (!m_trackSelector.has_value() || + m_trackSelector->isValidTrack(secondTrack)) { + auto destProxy = tracks.getTrack(tracks.addTrack()); + destProxy.copyFrom(secondTrack, true); + } + } } } } diff --git a/Examples/Python/src/TrackFinding.cpp b/Examples/Python/src/TrackFinding.cpp index d6fb7793e85..418f16b41ef 100644 --- a/Examples/Python/src/TrackFinding.cpp +++ b/Examples/Python/src/TrackFinding.cpp @@ -327,7 +327,7 @@ void addTrackFinding(Context& ctx) { ACTS_PYTHON_MEMBER(findTracks); ACTS_PYTHON_MEMBER(measurementSelectorCfg); ACTS_PYTHON_MEMBER(trackSelectorCfg); - ACTS_PYTHON_MEMBER(backward); + ACTS_PYTHON_MEMBER(twoWay); ACTS_PYTHON_MEMBER(maxSteps); ACTS_PYTHON_STRUCT_END(); } diff --git a/Examples/Scripts/Python/full_chain_odd.py b/Examples/Scripts/Python/full_chain_odd.py index 4d2e9d28428..26c951249eb 100755 --- a/Examples/Scripts/Python/full_chain_odd.py +++ b/Examples/Scripts/Python/full_chain_odd.py @@ -77,12 +77,12 @@ MomentumConfig(1.0 * u.GeV, 10.0 * u.GeV, transverse=True), EtaConfig(-3.0, 3.0), PhiConfig(0.0, 360.0 * u.degree), - ParticleConfig(1, acts.PdgParticle.eMuon, randomizeCharge=True), + ParticleConfig(4, acts.PdgParticle.eMuon, randomizeCharge=True), vtxGen=acts.examples.GaussianVertexGenerator( mean=acts.Vector4(0, 0, 0, 0), stddev=acts.Vector4(0.0125 * u.mm, 0.0125 * u.mm, 55.5 * u.mm, 1.0 * u.ns), ), - multiplicity=1, + multiplicity=200, rnd=rnd, ) else: @@ -164,55 +164,6 @@ outputDirRoot=outputDir, ) -s.addAlgorithm( - acts.examples.MyTrackFindingAlgorithm( - level=acts.logging.INFO, - - magneticField=field, - trackingGeometry=trackingGeometry, - - inputMeasurements="measurements", - inputSourceLinks="sourcelinks", - inputInitialTrackParameters="estimatedparameters", - outputTracks="mytracks", - - measurementSelectorCfg=acts.MeasurementSelector.Config( - [ - ( - acts.GeometryIdentifier(), - ( - [], - [30], - [10], - ), - ) - ] - ), - ) -) - -trackStatesWriter = acts.examples.RootTrackStatesWriter( - level=acts.logging.INFO, - inputTracks="mytracks", - inputParticles="particles_selected", - inputSimHits="simhits", - inputMeasurementParticlesMap="measurement_particles_map", - inputMeasurementSimHitsMap="measurement_simhits_map", - filePath=str(outputDir / f"trackstates_hi.root"), - treeName="trackstates", -) -s.addWriter(trackStatesWriter) - -trackSummaryWriter = acts.examples.RootTrackSummaryWriter( - level=acts.logging.INFO, - inputTracks="mytracks", - inputParticles="particles_selected", - inputMeasurementParticlesMap="measurement_particles_map", - filePath=str(outputDir / f"tracksummary_hi.root"), - treeName="tracksummary", -) -s.addWriter(trackSummaryWriter) - addCKFTracks( s, trackingGeometry, @@ -228,4 +179,33 @@ # outputDirCsv=outputDir, ) +if ambiguity_MLSolver: + addAmbiguityResolutionML( + s, + AmbiguityResolutionMLConfig( + maximumSharedHits=3, maximumIterations=1000000, nMeasurementsMin=7 + ), + outputDirRoot=outputDir, + # outputDirCsv=outputDir, + onnxModelFile=os.path.dirname(__file__) + + "/MLAmbiguityResolution/duplicateClassifier.onnx", + ) +else: + addAmbiguityResolution( + s, + AmbiguityResolutionConfig( + maximumSharedHits=3, maximumIterations=1000000, nMeasurementsMin=7 + ), + outputDirRoot=outputDir, + writeCovMat=True, + # outputDirCsv=outputDir, + ) + +addVertexFitting( + s, + field, + vertexFinder=VertexFinder.Iterative, + outputDirRoot=outputDir, +) + s.run() diff --git a/Examples/Scripts/Python/full_chain_odd_hi.py b/Examples/Scripts/Python/full_chain_odd_hi.py new file mode 100755 index 00000000000..4d2e9d28428 --- /dev/null +++ b/Examples/Scripts/Python/full_chain_odd_hi.py @@ -0,0 +1,231 @@ +#!/usr/bin/env python3 +import os, argparse, pathlib, acts, acts.examples +from acts.examples.simulation import ( + addParticleGun, + MomentumConfig, + EtaConfig, + PhiConfig, + ParticleConfig, + addPythia8, + addFatras, + addGeant4, + ParticleSelectorConfig, + addDigitization, +) +from acts.examples.reconstruction import ( + addSeeding, + TruthSeedRanges, + addCKFTracks, + TrackSelectorConfig, + addAmbiguityResolution, + AmbiguityResolutionConfig, + addAmbiguityResolutionML, + AmbiguityResolutionMLConfig, + addVertexFitting, + VertexFinder, +) +from common import getOpenDataDetectorDirectory +from acts.examples.odd import getOpenDataDetector + +parser = argparse.ArgumentParser(description="Full chain with the OpenDataDetector") + +parser.add_argument("--events", "-n", help="Number of events", type=int, default=100) +parser.add_argument( + "--geant4", help="Use Geant4 instead of fatras", action="store_true" +) +parser.add_argument( + "--ttbar", + help="Use Pythia8 (ttbar, pile-up 200) instead of particle gun", + action="store_true", +) +parser.add_argument( + "--MLSolver", + help="Use the Ml Ambiguity Solver instead of the classical one", + action="store_true", +) + +args = vars(parser.parse_args()) + +ttbar = args["ttbar"] +g4_simulation = args["geant4"] +ambiguity_MLSolver = args["MLSolver"] +u = acts.UnitConstants +geoDir = getOpenDataDetectorDirectory() +outputDir = pathlib.Path.cwd() / "odd_output" +# acts.examples.dump_args_calls(locals()) # show python binding calls + +oddMaterialMap = geoDir / "data/odd-material-maps.root" +oddDigiConfig = geoDir / "config/odd-digi-smearing-config.json" +oddSeedingSel = geoDir / "config/odd-seeding-config.json" +oddMaterialDeco = acts.IMaterialDecorator.fromFile(oddMaterialMap) + +detector, trackingGeometry, decorators = getOpenDataDetector( + geoDir, mdecorator=oddMaterialDeco +) +field = acts.ConstantBField(acts.Vector3(0.0, 0.0, 2.0 * u.T)) +rnd = acts.examples.RandomNumbers(seed=42) + +s = acts.examples.Sequencer( + events=args["events"], + numThreads=1, + outputDir=str(outputDir), +) + +if not ttbar: + addParticleGun( + s, + MomentumConfig(1.0 * u.GeV, 10.0 * u.GeV, transverse=True), + EtaConfig(-3.0, 3.0), + PhiConfig(0.0, 360.0 * u.degree), + ParticleConfig(1, acts.PdgParticle.eMuon, randomizeCharge=True), + vtxGen=acts.examples.GaussianVertexGenerator( + mean=acts.Vector4(0, 0, 0, 0), + stddev=acts.Vector4(0.0125 * u.mm, 0.0125 * u.mm, 55.5 * u.mm, 1.0 * u.ns), + ), + multiplicity=1, + rnd=rnd, + ) +else: + addPythia8( + s, + hardProcess=["Top:qqbar2ttbar=on"], + npileup=50, + vtxGen=acts.examples.GaussianVertexGenerator( + mean=acts.Vector4(0, 0, 0, 0), + stddev=acts.Vector4(0.0125 * u.mm, 0.0125 * u.mm, 55.5 * u.mm, 5.0 * u.ns), + ), + rnd=rnd, + outputDirRoot=outputDir, + # outputDirCsv=outputDir, + ) +if g4_simulation: + if s.config.numThreads != 1: + raise ValueError("Geant 4 simulation does not support multi-threading") + + # Pythia can sometime simulate particles outside the world volume, a cut on the Z of the track help mitigate this effect + # Older version of G4 might not work, this as has been tested on version `geant4-11-00-patch-03` + # For more detail see issue #1578 + addGeant4( + s, + detector, + trackingGeometry, + field, + preSelectParticles=ParticleSelectorConfig( + rho=(0.0, 24 * u.mm), + absZ=(0.0, 1.0 * u.m), + eta=(-3.0, 3.0), + pt=(150 * u.MeV, None), + removeNeutral=True, + ), + outputDirRoot=outputDir, + # outputDirCsv=outputDir, + rnd=rnd, + killVolume=trackingGeometry.worldVolume, + killAfterTime=25 * u.ns, + ) +else: + addFatras( + s, + trackingGeometry, + field, + preSelectParticles=ParticleSelectorConfig( + rho=(0.0, 24 * u.mm), + absZ=(0.0, 1.0 * u.m), + eta=(-3.0, 3.0), + pt=(150 * u.MeV, None), + removeNeutral=True, + ) + if ttbar + else ParticleSelectorConfig(), + enableInteractions=True, + outputDirRoot=outputDir, + # outputDirCsv=outputDir, + rnd=rnd, + ) + +addDigitization( + s, + trackingGeometry, + field, + digiConfigFile=oddDigiConfig, + outputDirRoot=outputDir, + # outputDirCsv=outputDir, + rnd=rnd, +) + +addSeeding( + s, + trackingGeometry, + field, + TruthSeedRanges(pt=(1.0 * u.GeV, None), eta=(-3.0, 3.0), nHits=(9, None)) + if ttbar + else TruthSeedRanges(), + geoSelectionConfigFile=oddSeedingSel, + outputDirRoot=outputDir, +) + +s.addAlgorithm( + acts.examples.MyTrackFindingAlgorithm( + level=acts.logging.INFO, + + magneticField=field, + trackingGeometry=trackingGeometry, + + inputMeasurements="measurements", + inputSourceLinks="sourcelinks", + inputInitialTrackParameters="estimatedparameters", + outputTracks="mytracks", + + measurementSelectorCfg=acts.MeasurementSelector.Config( + [ + ( + acts.GeometryIdentifier(), + ( + [], + [30], + [10], + ), + ) + ] + ), + ) +) + +trackStatesWriter = acts.examples.RootTrackStatesWriter( + level=acts.logging.INFO, + inputTracks="mytracks", + inputParticles="particles_selected", + inputSimHits="simhits", + inputMeasurementParticlesMap="measurement_particles_map", + inputMeasurementSimHitsMap="measurement_simhits_map", + filePath=str(outputDir / f"trackstates_hi.root"), + treeName="trackstates", +) +s.addWriter(trackStatesWriter) + +trackSummaryWriter = acts.examples.RootTrackSummaryWriter( + level=acts.logging.INFO, + inputTracks="mytracks", + inputParticles="particles_selected", + inputMeasurementParticlesMap="measurement_particles_map", + filePath=str(outputDir / f"tracksummary_hi.root"), + treeName="tracksummary", +) +s.addWriter(trackSummaryWriter) + +addCKFTracks( + s, + trackingGeometry, + field, + TrackSelectorConfig( + pt=(1.0 * u.GeV if ttbar else 0.0, None), + absEta=(None, 3.0), + loc0=(-4.0 * u.mm, 4.0 * u.mm), + nMeasurementsMin=7, + ), + outputDirRoot=outputDir, + writeCovMat=True, + # outputDirCsv=outputDir, +) + +s.run() From 15485a7cdefc1cf3fce8a143149d2e66c467d19d Mon Sep 17 00:00:00 2001 From: andiwand Date: Thu, 23 Nov 2023 15:23:28 +0100 Subject: [PATCH 3/4] fix --- .../src/TrackFindingAlgorithm.cpp | 23 +++++++++++++++++-- .../python/acts/examples/reconstruction.py | 2 ++ 2 files changed, 23 insertions(+), 2 deletions(-) diff --git a/Examples/Algorithms/TrackFinding/src/TrackFindingAlgorithm.cpp b/Examples/Algorithms/TrackFinding/src/TrackFindingAlgorithm.cpp index f0aa3470cd2..9b0234633ef 100644 --- a/Examples/Algorithms/TrackFinding/src/TrackFindingAlgorithm.cpp +++ b/Examples/Algorithms/TrackFinding/src/TrackFindingAlgorithm.cpp @@ -125,7 +125,7 @@ ActsExamples::ProcessCode ActsExamples::TrackFindingAlgorithm::execute( ActsExamples::TrackFindingAlgorithm::TrackFinderOptions secondOptions( ctx.geoContext, ctx.magFieldContext, ctx.calibContext, slAccessorDelegate, - extensions, firstPropOptions, pSurface.get()); + extensions, secondPropOptions, pSurface.get()); secondOptions.filterTargetSurface = pSurface.get(); secondOptions.smoothing = true; secondOptions.smoothingTargetSurfaceStrategy = @@ -178,7 +178,26 @@ ActsExamples::ProcessCode ActsExamples::TrackFindingAlgorithm::execute( } if (m_cfg.twoWay) { - auto secondResult = (*m_cfg.findTracks)(initialParameters.at(iseed), + std::optional firstState; + for (auto st : firstTrack.trackStatesReversed()) { + bool isMeasurement = + st.typeFlags().test(Acts::TrackStateFlag::MeasurementFlag); + bool isOutlier = + st.typeFlags().test(Acts::TrackStateFlag::OutlierFlag); + // We are excluding non measurement states and outlier here. Those can + // decrease resolution because only the smoothing corrected the very + // first prediction as filtering is not possible. + if (isMeasurement && !isOutlier) { + firstState = st; + } + } + + Acts::BoundTrackParameters secondInitialParameters( + firstState->referenceSurface().getSharedPtr(), + firstState->parameters(), firstState->covariance(), + initialParameters.at(iseed).particleHypothesis()); + + auto secondResult = (*m_cfg.findTracks)(secondInitialParameters, secondOptions, tracksTemp); if (!secondResult.ok()) { diff --git a/Examples/Python/python/acts/examples/reconstruction.py b/Examples/Python/python/acts/examples/reconstruction.py index cfc8de476ff..8dc6290dfcc 100644 --- a/Examples/Python/python/acts/examples/reconstruction.py +++ b/Examples/Python/python/acts/examples/reconstruction.py @@ -1012,6 +1012,7 @@ def addCKFTracks( field: acts.MagneticFieldProvider, trackSelectorConfig: Optional[TrackSelectorConfig] = None, ckfConfig: CkfConfig = CkfConfig(), + twoWay: bool = True, outputDirCsv: Optional[Union[Path, str]] = None, outputDirRoot: Optional[Union[Path, str]] = None, writeTrajectories: bool = True, @@ -1086,6 +1087,7 @@ def addCKFTracks( ), **acts.examples.defaultKWArgs( maxSteps=ckfConfig.maxSteps, + twoWay=twoWay, ), ) s.addAlgorithm(trackFinder) From 82ed722d3191695a1843b11f5b194aa07962e5f9 Mon Sep 17 00:00:00 2001 From: andiwand Date: Thu, 23 Nov 2023 16:32:21 +0100 Subject: [PATCH 4/4] deduplicate a bit --- .../src/TrackFindingAlgorithm.cpp | 20 ++++++++++++------- 1 file changed, 13 insertions(+), 7 deletions(-) diff --git a/Examples/Algorithms/TrackFinding/src/TrackFindingAlgorithm.cpp b/Examples/Algorithms/TrackFinding/src/TrackFindingAlgorithm.cpp index 9b0234633ef..44f945ab15d 100644 --- a/Examples/Algorithms/TrackFinding/src/TrackFindingAlgorithm.cpp +++ b/Examples/Algorithms/TrackFinding/src/TrackFindingAlgorithm.cpp @@ -171,13 +171,13 @@ ActsExamples::ProcessCode ActsExamples::TrackFindingAlgorithm::execute( for (auto& firstTrack : firstTracksForSeed) { seedNumber(firstTrack) = nSeed; - if (!m_trackSelector.has_value() || - m_trackSelector->isValidTrack(firstTrack)) { - auto destProxy = tracks.getTrack(tracks.addTrack()); - destProxy.copyFrom(firstTrack, true); - } - - if (m_cfg.twoWay) { + if (!m_cfg.twoWay) { + if (!m_trackSelector.has_value() || + m_trackSelector->isValidTrack(firstTrack)) { + auto destProxy = tracks.getTrack(tracks.addTrack()); + destProxy.copyFrom(firstTrack, true); + } + } else { std::optional firstState; for (auto st : firstTrack.trackStatesReversed()) { bool isMeasurement = @@ -213,6 +213,12 @@ ActsExamples::ProcessCode ActsExamples::TrackFindingAlgorithm::execute( auto& secondTracksForSeed = secondResult.value(); for (auto& secondTrack : secondTracksForSeed) { if (secondTrack.nTrackStates() < 2) { + if (!m_trackSelector.has_value() || + m_trackSelector->isValidTrack(firstTrack)) { + auto destProxy = tracks.getTrack(tracks.addTrack()); + destProxy.copyFrom(firstTrack, true); + } + continue; }