diff --git a/JugTrack/CMakeLists.txt b/JugTrack/CMakeLists.txt index ffec289083c8cb193b10bcd493330e0e93cd6bc0..53425a5901895fbd8fa6bdb27152bf1336ea8605 100644 --- a/JugTrack/CMakeLists.txt +++ b/JugTrack/CMakeLists.txt @@ -33,8 +33,8 @@ gaudi_add_module(JugTrackPlugins src/components/TrackerSourcesLinker.cpp src/components/TrackFindingAlgorithm.cpp src/components/TrackFindingAlgorithmFunction.cpp - #src/components/TrackFittingAlgorithm.cpp - #src/components/TrackFittingFunction.cpp + src/components/TrackFittingAlgorithm.cpp + src/components/TrackFittingFunction.cpp src/components/TestACTSLogger.cpp src/components/TrackParamTruthInit.cpp src/components/TrackParamClusterInit.cpp diff --git a/JugTrack/JugTrack/SimMultiTrajectory.hpp b/JugTrack/JugTrack/SimMultiTrajectory.hpp deleted file mode 100644 index d2bf7ab5a6ad65cdc179e0852dbe34db16bcb7dc..0000000000000000000000000000000000000000 --- a/JugTrack/JugTrack/SimMultiTrajectory.hpp +++ /dev/null @@ -1,161 +0,0 @@ -// This file is part of the Acts project. -// -// Copyright (C) 2019 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 "ACTFW/Validation/ProtoTrackClassification.hpp" -#include "Acts/EventData/MultiTrajectory.hpp" -#include "Acts/EventData/TrackParameters.hpp" - -#include "JugTrack/IndexSourceLink.hpp" -#include <unordered_map> -#include <utility> - -namespace Jug { - - /// Associate a particle to its hit count within a proto track. - struct ParticleHitCount { - uint64_t particleId; - std::size_t hitCount; - }; - -using IndexedParams = std::unordered_map<size_t, Acts::BoundTrackParameters>; - -/// @brief Struct for truth track fitting/finding result with -/// Acts::KalmanFitter/Acts::CombinatorialKalmanFilter -/// -/// It contains a MultiTrajectory with a vector of entry indices for individual -/// trajectories, and a map of fitted parameters indexed by the entry index. -/// In case of track fitting, there is at most one trajectory in the -/// MultiTrajectory; In case of track finding, there could be multiple -/// trajectories in the MultiTrajectory. -/// -/// @note The MultiTrajectory is thought to be empty if there is no entry index -struct SimMultiTrajectory { - public: - /// @brief Default constructor - /// - SimMultiTrajectory() = default; - - /// @brief Constructor from multiTrajectory and fitted track parameters - /// - /// @param multiTraj The multiTrajectory - /// @param tTips The entry indices for trajectories in multiTrajectory - /// @param parameters The fitted track parameters indexed by trajectory entry - /// index - SimMultiTrajectory(const Acts::MultiTrajectory<IndexSourceLink>& multiTraj, - const std::vector<size_t>& tTips, - const IndexedParams& parameters) - : m_multiTrajectory(multiTraj), - m_trackTips(tTips), - m_trackParameters(parameters) {} - - /// @brief Copy constructor - /// - /// @param rhs The source SimMultiTrajectory - SimMultiTrajectory(const SimMultiTrajectory& rhs) - : m_multiTrajectory(rhs.m_multiTrajectory), - m_trackTips(rhs.m_trackTips), - m_trackParameters(rhs.m_trackParameters) {} - - /// @brief Copy move constructor - /// - /// @param rhs The source SimMultiTrajectory - SimMultiTrajectory(SimMultiTrajectory&& rhs) - : m_multiTrajectory(std::move(rhs.m_multiTrajectory)), - m_trackTips(std::move(rhs.m_trackTips)), - m_trackParameters(std::move(rhs.m_trackParameters)) {} - - /// @brief Default destructor - /// - ~SimMultiTrajectory() = default; - - /// @brief assignment operator - /// - /// @param rhs The source SimMultiTrajectory - SimMultiTrajectory& operator=(const SimMultiTrajectory& rhs) { - m_multiTrajectory = rhs.m_multiTrajectory; - m_trackTips = rhs.m_trackTips; - m_trackParameters = rhs.m_trackParameters; - return *this; - } - - /// @brief assignment move operator - /// - /// @param rhs The source SimMultiTrajectory - SimMultiTrajectory& operator=(SimMultiTrajectory&& rhs) { - m_multiTrajectory = std::move(rhs.m_multiTrajectory); - m_trackTips = std::move(rhs.m_trackTips); - m_trackParameters = std::move(rhs.m_trackParameters); - return *this; - } - - /// @brief Indicator if a trajectory exists - /// - /// @param entryIndex The trajectory entry index - /// - /// @return Whether there is trajectory with provided entry index - bool hasTrajectory(const size_t& entryIndex) const { - return std::count(m_trackTips.begin(), m_trackTips.end(), entryIndex) > 0; - } - - /// @brief Indicator if there is fitted track parameters for one trajectory - /// - /// @param entryIndex The trajectory entry index - /// - /// @return Whether having fitted track parameters or not - bool hasTrackParameters(const size_t& entryIndex) const { - return m_trackParameters.count(entryIndex) > 0; - } - - /// @brief Getter for multiTrajectory - /// - /// @return The multiTrajectory with trajectory entry indices - /// - /// @note It could return an empty multiTrajectory - std::pair<std::vector<size_t>, Acts::MultiTrajectory<IndexSourceLink>> - trajectory() const { - return std::make_pair(m_trackTips, m_multiTrajectory); - } - - /// @brief Getter of fitted track parameters for one trajectory - /// - /// @param entryIndex The trajectory entry index - /// - /// @return The fitted track parameters of the trajectory - const Acts::BoundTrackParameters& trackParameters(const size_t& entryIndex) const { - auto it = m_trackParameters.find(entryIndex); - if (it != m_trackParameters.end()) { - return it->second; - } else { - throw std::runtime_error( - "No fitted track parameters for trajectory with entry index = " + - std::to_string(entryIndex)); - } - } - - /// @brief Counter of associated truth particles for one trajectory - /// - /// @param entryIndex The trajectory entry index - /// - /// @return The truth particle counts in ascending order - std::vector<ParticleHitCount> identifyMajorityParticle( - const size_t& entryIndex) const; - - private: - // The multiTrajectory - Acts::MultiTrajectory<IndexSourceLink> m_multiTrajectory; - - // The entry indices of trajectories stored in multiTrajectory - std::vector<size_t> m_trackTips = {}; - - // The fitted parameters at the provided surface for individual trajectories - IndexedParams m_trackParameters = {}; -}; - -} // namespace FW diff --git a/JugTrack/JugTrack/Track.hpp b/JugTrack/JugTrack/Track.hpp index 533dc42132b7272089afb680f2699f95a0b42f3d..ec1ee629f7b8e00f14c323c08ff0c67af5f0b9b7 100644 --- a/JugTrack/JugTrack/Track.hpp +++ b/JugTrack/JugTrack/Track.hpp @@ -1,35 +1,22 @@ // This file is part of the Acts project. // -// Copyright (C) 2019 CERN for the benefit of the Acts project +// Copyright (C) 2019-2020 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/. -/// @file -/// @brief All track-related shared types. - #pragma once -#include "JugTrack/SimMultiTrajectory.hpp" -//#include "ACTFW/EventData/SimSourceLink.hpp" -#include "Acts/EventData/MultiTrajectory.hpp" #include "Acts/EventData/TrackParameters.hpp" -#include "JugTrack/IndexSourceLink.hpp" #include <vector> namespace Jug { - /// (Reconstructed) track parameters e.g. close to the vertex. - using TrackParameters = ::Acts::BoundTrackParameters; - /// Container of reconstructed track states for multiple tracks. - using TrackParametersContainer = std::vector<TrackParameters>; - - /// MultiTrajectory definition - using Trajectory = Acts::MultiTrajectory<IndexSourceLink>; - - /// Container for the truth fitting/finding track(s) - using TrajectoryContainer = std::vector<SimMultiTrajectory>; +/// (Reconstructed) track parameters e.g. close to the vertex. +using TrackParameters = ::Acts::BoundTrackParameters; +/// Container of reconstructed track states for multiple tracks. +using TrackParametersContainer = std::vector<TrackParameters>; -} // namespace Jug +} // namespace ActsExamples diff --git a/JugTrack/JugTrack/Trajectories.hpp b/JugTrack/JugTrack/Trajectories.hpp new file mode 100644 index 0000000000000000000000000000000000000000..30d1983dfed683f9c559ed04bded319c79862d9b --- /dev/null +++ b/JugTrack/JugTrack/Trajectories.hpp @@ -0,0 +1,95 @@ +#ifndef JugTrack_Trajectories_HH +#define JugTrack_Trajectories_HH + +#include "Acts/EventData/MultiTrajectory.hpp" +#include "JugTrack/IndexSourceLink.hpp" +#include "JugTrack/Track.hpp" + +#include <algorithm> +#include <unordered_map> +#include <vector> + +namespace Jug { + +/// Store reconstructed trajectories from track finding/fitting. +/// +/// It contains a MultiTrajectory with a vector of entry indices for +/// individual trajectories, and a map of fitted parameters indexed by the +/// entry index. In case of track fitting, there is at most one trajectory +/// in the MultiTrajectory; In case of track finding, there could be +/// multiple trajectories in the MultiTrajectory. +struct Trajectories final { + public: + /// (Reconstructed) trajectory with multiple states. + using MultiTrajectory = ::Acts::MultiTrajectory<IndexSourceLink>; + /// Fitted parameters identified by indices in the multi trajectory. + using IndexedParameters = std::unordered_map<size_t, TrackParameters>; + + /// Default construct an empty object. Required for container compatibility + /// and to signal an error. + Trajectories() = default; + /// Construct from fitted multi trajectory and parameters. + /// + /// @param multiTraj The multi trajectory + /// @param tTips Tip indices that identify valid trajectories + /// @param parameters Fitted track parameters indexed by trajectory index + Trajectories(const MultiTrajectory& multiTraj, + const std::vector<size_t>& tTips, + const IndexedParameters& parameters) + : m_multiTrajectory(multiTraj), + m_trackTips(tTips), + m_trackParameters(parameters) {} + + /// Return true if there exists no valid trajectory. + bool empty() const { return m_trackTips.empty(); } + + /// Access the underlying multi trajectory. + const MultiTrajectory& multiTrajectory() const { return m_multiTrajectory; } + + /// Access the tip indices that identify valid trajectories. + const std::vector<size_t>& tips() const { return m_trackTips; } + + /// Check if a trajectory exists for the given index. + /// + /// @param entryIndex The trajectory entry index + /// @return Whether there is trajectory with provided entry index + bool hasTrajectory(size_t entryIndex) const { + return (0 < std::count(m_trackTips.begin(), m_trackTips.end(), entryIndex)); + } + + /// Check if fitted track parameters exists for the given index. + /// + /// @param entryIndex The trajectory entry index + /// @return Whether having fitted track parameters or not + bool hasTrackParameters(size_t entryIndex) const { + return (0 < m_trackParameters.count(entryIndex)); + } + + /// Access the fitted track parameters for the given index. + /// + /// @param entryIndex The trajectory entry index + /// @return The fitted track parameters of the trajectory + const TrackParameters& trackParameters(size_t entryIndex) const { + auto it = m_trackParameters.find(entryIndex); + if (it == m_trackParameters.end()) { + throw std::runtime_error( + "No fitted track parameters for trajectory with entry index = " + + std::to_string(entryIndex)); + } + return it->second; + } + + private: + // The multiTrajectory + MultiTrajectory m_multiTrajectory; + // The entry indices of trajectories stored in multiTrajectory + std::vector<size_t> m_trackTips = {}; + // The fitted parameters at the provided surface for individual trajectories + IndexedParameters m_trackParameters = {}; +}; + +/// Container for multiple trajectories. +using TrajectoriesContainer = std::vector<Trajectories>; + +} // namespace Jug +#endif diff --git a/JugTrack/src/components/ParticlesFromTrackFit.cpp b/JugTrack/src/components/ParticlesFromTrackFit.cpp index b06fad29db23a033909f344a426c3abfa301750b..a2811bbfd43c151c22e16b91f437e3ebf3f1b75d 100644 --- a/JugTrack/src/components/ParticlesFromTrackFit.cpp +++ b/JugTrack/src/components/ParticlesFromTrackFit.cpp @@ -24,6 +24,7 @@ #include "eicd/TrackParametersCollection.h" #include "JugTrack/IndexSourceLink.hpp" #include "JugTrack/Track.hpp" +#include "JugTrack/Trajectories.hpp" #include "Acts/Utilities/Helpers.hpp" @@ -37,7 +38,7 @@ namespace Jug { class ParticlesFromTrackFit : public GaudiAlgorithm { public: //DataHandle<eic::RawTrackerHitCollection> m_inputHitCollection{"inputHitCollection", Gaudi::DataHandle::Reader, this}; - DataHandle<TrajectoryContainer> m_inputTrajectories{"inputTrajectories", Gaudi::DataHandle::Reader, this}; + DataHandle<TrajectoriesContainer> m_inputTrajectories{"inputTrajectories", Gaudi::DataHandle::Reader, this}; DataHandle<eic::ParticleCollection> m_outputParticles{"outputParticles", Gaudi::DataHandle::Writer, this}; DataHandle<eic::TrackParametersCollection> m_outputTrackParameters{"outputTrackParameters", Gaudi::DataHandle::Writer, this}; @@ -58,92 +59,105 @@ namespace Jug { StatusCode execute() override { // input collection - const TrajectoryContainer* trajectories = m_inputTrajectories.get(); + const TrajectoriesContainer* trajectories = m_inputTrajectories.get(); // create output collections auto rec_parts = m_outputParticles.createAndPut(); auto track_pars = m_outputTrackParameters.createAndPut(); debug() << std::size(*trajectories) << " trajectories " << endmsg; - for(const auto& traj : *trajectories) { - //traj.trajectory().first - const auto& [trackTips, mj] = traj.trajectory(); - if (trackTips.empty()) { - debug() << "Empty multiTrajectory." << endmsg; - continue; - } - - // Get the entry index for the single trajectory - auto& trackTip = trackTips.front(); - - // Collect the trajectory summary info - auto trajState = Acts::MultiTrajectoryHelpers::trajectoryState(mj, trackTip); - int m_nMeasurements = trajState.nMeasurements; - int m_nStates = trajState.nStates; - - // Get the fitted track parameter - bool m_hasFittedParams = false; - if (traj.hasTrackParameters(trackTip)) { - m_hasFittedParams = true; - const auto& boundParam = traj.trackParameters(trackTip); - const auto& parameter = boundParam.parameters(); - const auto& covariance = *boundParam.covariance(); - debug() << "loc 0 = " << parameter[Acts::eBoundLoc0] << endmsg; - debug() << "loc 1 = " << parameter[Acts::eBoundLoc1] << endmsg; - debug() << "phi = " << parameter[Acts::eBoundPhi] << endmsg; - debug() << "theta = " << parameter[Acts::eBoundTheta] << endmsg; - debug() << "q/p = " << parameter[Acts::eBoundQOverP] << endmsg; - debug() << "p = " << 1.0/parameter[Acts::eBoundQOverP] << endmsg; - - debug() << "err phi = " << sqrt(covariance(Acts::eBoundPhi, Acts::eBoundPhi)) << endmsg; - debug() << "err th = " << sqrt(covariance(Acts::eBoundTheta, Acts::eBoundTheta)) << endmsg; - debug() << "err q/p = " << sqrt(covariance(Acts::eBoundQOverP, Acts::eBoundQOverP))<< endmsg; - - debug() << " chi2 = " << trajState.chi2Sum << endmsg; - - eic::TrackParameters pars({ - parameter[Acts::eBoundLoc0], parameter[Acts::eBoundLoc1], parameter[Acts::eBoundPhi], - parameter[Acts::eBoundTheta], parameter[Acts::eBoundQOverP],parameter[Acts::eBoundTime], - sqrt(covariance(Acts::eBoundLoc0, Acts::eBoundLoc0)), sqrt(covariance(Acts::eBoundLoc1, Acts::eBoundLoc1)), - sqrt(covariance(Acts::eBoundPhi, Acts::eBoundPhi)), sqrt(covariance(Acts::eBoundTheta, Acts::eBoundTheta)), - sqrt(covariance(Acts::eBoundQOverP, Acts::eBoundQOverP)),sqrt(covariance(Acts::eBoundTime, Acts::eBoundTime))}); - track_pars->push_back(pars); - - //m_ePHI_fit = parameter[Acts::eBoundPhi]; - //m_eTHETA_fit = parameter[Acts::eBoundTheta]; - //m_eQOP_fit = parameter[Acts::eBoundQOverP]; - //m_eT_fit = parameter[Acts::eBoundTime]; - //m_err_eLOC0_fit = sqrt(covariance(Acts::eBoundLoc0, Acts::eBoundLoc0)); - //m_err_eLOC1_fit = sqrt(covariance(Acts::eBoundLoc1, Acts::eBoundLoc1)); - //m_err_ePHI_fit = sqrt(covariance(Acts::eBoundPhi, Acts::eBoundPhi)); - //m_err_eTHETA_fit = sqrt(covariance(Acts::eBoundTheta, Acts::eBoundTheta)); - //m_err_eQOP_fit = sqrt(covariance(Acts::eBoundQOverP, Acts::eBoundQOverP)); - //m_err_eT_fit = sqrt(covariance(Acts::eBoundTime, Acts::eBoundTime)); - } + // Loop over the trajectories + for (size_t itraj = 0; itraj < trajectories->size(); ++itraj) { + const auto& traj = (*trajectories)[itraj]; - auto tsize = traj.trajectory().first.size(); - debug() << "# fitted parameters : " << tsize << endmsg; - if(tsize == 0 ) continue; - traj.trajectory().second.visitBackwards(tsize-1, [&](auto&& trackstate) { - //debug() << trackstate.hasPredicted() << endmsg; - //debug() << trackstate.predicted() << endmsg; - auto params = trackstate.predicted() ;//<< endmsg; - - double p0 = (1.0 / params[Acts::eBoundQOverP]) / Acts::UnitConstants::GeV; - debug() << "track predicted p = " << p0 << " GeV" << endmsg; - if ( std::abs(p0) > 500) { - debug() << "skipping" << endmsg; - return; + // Get the entry index for the single trajectory + // The trajectory entry indices and the multiTrajectory + const auto& mj = traj.multiTrajectory(); + const auto& trackTips = traj.tips(); + if (trackTips.empty()) { + debug() << "Empty multiTrajectory." << endmsg; + continue; } - eic::Particle p({params[Acts::eBoundPhi], params[Acts::eBoundTheta], 1.0 / std::abs(params[Acts::eBoundQOverP]), 0.000511}, - {0.0, 0.0, 0.0, params[Acts::eBoundTime]}, - (long long)11 * params[Acts::eBoundQOverP] / std::abs(params[Acts::eBoundQOverP]), 0); - //debug() << p << endmsg; - rec_parts->push_back(p); - }); + auto& trackTip = trackTips.front(); + + // Collect the trajectory summary info + auto trajState = Acts::MultiTrajectoryHelpers::trajectoryState(mj, trackTip); + int m_nMeasurements = trajState.nMeasurements; + int m_nStates = trajState.nStates; + + // Get the fitted track parameter + bool m_hasFittedParams = false; + if (traj.hasTrackParameters(trackTip)) { + m_hasFittedParams = true; + const auto& boundParam = traj.trackParameters(trackTip); + const auto& parameter = boundParam.parameters(); + const auto& covariance = *boundParam.covariance(); + debug() << "loc 0 = " << parameter[Acts::eBoundLoc0] << endmsg; + debug() << "loc 1 = " << parameter[Acts::eBoundLoc1] << endmsg; + debug() << "phi = " << parameter[Acts::eBoundPhi] << endmsg; + debug() << "theta = " << parameter[Acts::eBoundTheta] << endmsg; + debug() << "q/p = " << parameter[Acts::eBoundQOverP] << endmsg; + debug() << "p = " << 1.0 / parameter[Acts::eBoundQOverP] << endmsg; + + debug() << "err phi = " << sqrt(covariance(Acts::eBoundPhi, Acts::eBoundPhi)) << endmsg; + debug() << "err th = " << sqrt(covariance(Acts::eBoundTheta, Acts::eBoundTheta)) + << endmsg; + debug() << "err q/p = " << sqrt(covariance(Acts::eBoundQOverP, Acts::eBoundQOverP)) + << endmsg; + + debug() << " chi2 = " << trajState.chi2Sum << endmsg; + + eic::TrackParameters pars({parameter[Acts::eBoundLoc0], parameter[Acts::eBoundLoc1], + parameter[Acts::eBoundPhi], parameter[Acts::eBoundTheta], + parameter[Acts::eBoundQOverP], parameter[Acts::eBoundTime], + sqrt(covariance(Acts::eBoundLoc0, Acts::eBoundLoc0)), + sqrt(covariance(Acts::eBoundLoc1, Acts::eBoundLoc1)), + sqrt(covariance(Acts::eBoundPhi, Acts::eBoundPhi)), + sqrt(covariance(Acts::eBoundTheta, Acts::eBoundTheta)), + sqrt(covariance(Acts::eBoundQOverP, Acts::eBoundQOverP)), + sqrt(covariance(Acts::eBoundTime, Acts::eBoundTime))}); + track_pars->push_back(pars); + + // m_ePHI_fit = parameter[Acts::eBoundPhi]; + // m_eTHETA_fit = parameter[Acts::eBoundTheta]; + // m_eQOP_fit = parameter[Acts::eBoundQOverP]; + // m_eT_fit = parameter[Acts::eBoundTime]; + // m_err_eLOC0_fit = sqrt(covariance(Acts::eBoundLoc0, Acts::eBoundLoc0)); + // m_err_eLOC1_fit = sqrt(covariance(Acts::eBoundLoc1, Acts::eBoundLoc1)); + // m_err_ePHI_fit = sqrt(covariance(Acts::eBoundPhi, Acts::eBoundPhi)); + // m_err_eTHETA_fit = sqrt(covariance(Acts::eBoundTheta, Acts::eBoundTheta)); + // m_err_eQOP_fit = sqrt(covariance(Acts::eBoundQOverP, Acts::eBoundQOverP)); + // m_err_eT_fit = sqrt(covariance(Acts::eBoundTime, Acts::eBoundTime)); + } + auto tsize = trackTips.size(); + debug() << "# fitted parameters : " << tsize << endmsg; + if (tsize == 0) + continue; + + mj.visitBackwards(tsize - 1, [&](auto&& trackstate) { + // debug() << trackstate.hasPredicted() << endmsg; + // debug() << trackstate.predicted() << endmsg; + auto params = trackstate.predicted(); //<< endmsg; + + double p0 = (1.0 / params[Acts::eBoundQOverP]) / Acts::UnitConstants::GeV; + debug() << "track predicted p = " << p0 << " GeV" << endmsg; + if (std::abs(p0) > 500) { + debug() << "skipping" << endmsg; + return; + } + + eic::Particle p({params[Acts::eBoundPhi], params[Acts::eBoundTheta], + 1.0 / std::abs(params[Acts::eBoundQOverP]), 0.000511}, + {0.0, 0.0, 0.0, params[Acts::eBoundTime]}, + (long long)11 * params[Acts::eBoundQOverP] / + std::abs(params[Acts::eBoundQOverP]), + 0); + // debug() << p << endmsg; + rec_parts->push_back(p); + }); } return StatusCode::SUCCESS; } diff --git a/JugTrack/src/components/TrackFindingAlgorithm.cpp b/JugTrack/src/components/TrackFindingAlgorithm.cpp index cd8c1a8f12328bbf91cd7b4d1eefae2b4755e81c..b34807b1a2b66d6d479685151d05694bd8492295 100644 --- a/JugTrack/src/components/TrackFindingAlgorithm.cpp +++ b/JugTrack/src/components/TrackFindingAlgorithm.cpp @@ -47,6 +47,7 @@ #include <random> #include <stdexcept> + namespace Jug::Reco { using namespace Acts::UnitLiterals; diff --git a/JugTrack/src/components/TrackFindingAlgorithm.h b/JugTrack/src/components/TrackFindingAlgorithm.h index ec934e6ce43a14029f2e35e2065a75dd3ddeafbc..d119633e8c9f2387ee983ccd64ffa8499055d1e1 100644 --- a/JugTrack/src/components/TrackFindingAlgorithm.h +++ b/JugTrack/src/components/TrackFindingAlgorithm.h @@ -21,6 +21,7 @@ #include "JugTrack/IndexSourceLink.hpp" #include "JugTrack/Measurement.hpp" #include "JugTrack/Track.hpp" +#include "JugTrack/Trajectories.hpp" #include "eicd/TrackerHitCollection.h" @@ -59,7 +60,7 @@ namespace Jug::Reco { Gaudi::DataHandle::Reader, this}; DataHandle<TrackParametersContainer> m_inputInitialTrackParameters{ "inputInitialTrackParameters", Gaudi::DataHandle::Reader, this}; - DataHandle<TrajectoryContainer> m_outputTrajectories{"outputTrajectories", + DataHandle<TrajectoriesContainer> m_outputTrajectories{"outputTrajectories", Gaudi::DataHandle::Writer, this}; TrackFinderFunction m_trackFinderFunc; SmartIF<IGeoSvc> m_geoSvc; diff --git a/JugTrack/src/components/TrackFittingAlgorithm.cpp b/JugTrack/src/components/TrackFittingAlgorithm.cpp index fe0ae75ed62cf5b9546f543fbc7bbb80c9a520ca..5e12f2f08b08d395cee51ab0ca50277c69a9d2ab 100644 --- a/JugTrack/src/components/TrackFittingAlgorithm.cpp +++ b/JugTrack/src/components/TrackFittingAlgorithm.cpp @@ -30,11 +30,10 @@ #include "JugBase/DataHandle.h" #include "JugBase/IGeoSvc.h" #include "JugTrack/GeometryContainers.hpp" -#include "JugTrack/SourceLinks.h" +#include "JugTrack/IndexSourceLink.hpp" #include "JugTrack/Track.hpp" #include "JugTrack/BField.h" #include "JugTrack/Measurement.hpp" -#include "JugTrack/SourceLinks.h" #include "eicd/TrackerHitCollection.h" @@ -53,6 +52,7 @@ namespace Jug::Reco { { declareProperty("inputSourceLinks", m_inputSourceLinks, ""); declareProperty("initialTrackParameters", m_initialTrackParameters, ""); + declareProperty("inputMeasurements", m_inputMeasurements, ""); declareProperty("foundTracks", m_foundTracks, ""); declareProperty("outputTrajectories", m_outputTrajectories, ""); } @@ -80,8 +80,9 @@ namespace Jug::Reco { StatusCode TrackFittingAlgorithm::execute() { // Read input data - const SourceLinkContainer* src_links = m_inputSourceLinks.get(); + const IndexSourceLinkContainer* src_links = m_inputSourceLinks.get(); const TrackParametersContainer* init_trk_params = m_initialTrackParameters.get(); + const MeasurementContainer* measurements = m_inputMeasurements.get(); // TrajectoryContainer trajectories; auto trajectories = m_outputTrajectories.createAndPut(); @@ -91,39 +92,83 @@ namespace Jug::Reco { auto pSurface = Acts::Surface::makeShared<Acts::PerigeeSurface>(Acts::Vector3{0., 0., 0.}); ACTS_LOCAL_LOGGER(Acts::getDefaultLogger("TrackFittingAlgorithm Logger", Acts::Logging::INFO)); + std::vector<IndexSourceLink> trackSourceLinks; + // Perform the track finding for each starting parameter // @TODO: use seeds from track seeding algorithm as starting parameter + // initial track params and proto tracks might likely have the same size. for (std::size_t iseed = 0; iseed < init_trk_params->size(); ++iseed) { + // this will eventually be per-track. for now we assume all src links are from the same single track! + for (auto& lnk : (*src_links)) { + // auto sourceLink = sourceLinks.nth(hitIndex); + trackSourceLinks.push_back(lnk); + } + + const auto& initialParams = (*init_trk_params)[iseed]; Acts::PropagatorPlainOptions pOptions; pOptions.maxSteps = 10000; - Acts::KalmanFitterOptions<Acts::VoidOutlierFinder> kfOptions( - m_geoctx, m_fieldctx, m_calibctx, - Acts::VoidOutlierFinder(), Acts::LoggerWrapper{logger()}, Acts::PropagatorPlainOptions(), &(*pSurface)); + //Acts::KalmanFitterOptions<Acts::VoidOutlierFinder> kfOptions( + // m_geoctx, m_fieldctx, m_calibctx, Acts::VoidOutlierFinder(), + // Acts::LoggerWrapper{logger()}, Acts::PropagatorPlainOptions(), &(*pSurface)); + // // Set the KalmanFitter options + Acts::KalmanFitterOptions<MeasurementCalibrator, Acts::VoidOutlierFinder> kfOptions( + m_geoctx, m_fieldctx, m_calibctx, + MeasurementCalibrator(*measurements), Acts::VoidOutlierFinder(), + Acts::LoggerWrapper{logger()}, Acts::PropagatorPlainOptions(), &(*pSurface)); debug() << "Invoke track fitting ... " << iseed << endmsg; - auto result = m_trackFittingFunc(*src_links, initialParams, kfOptions); + auto result = fitTrack(trackSourceLinks, initialParams, kfOptions); debug() << "fitting done." << endmsg; - if (result.ok()) { - // Get the track finding output object - const auto& trackFindingOutput = result.value(); - // Create a SimMultiTrajectory - trajectories->emplace_back(std::move(trackFindingOutput.fittedStates), std::move(trackFindingOutput.trackTips), - std::move(trackFindingOutput.fittedParameters)); + // if (result.ok()) { + // // Get the track finding output object + // const auto& trackFindingOutput = result.value(); + // // Create a SimMultiTrajectory + // trajectories->emplace_back(std::move(trackFindingOutput.fittedStates), + // std::move(trackFindingOutput.trackTips), + // std::move(trackFindingOutput.fittedParameters)); + //} else { + // debug() << "Track finding failed for truth seed " << iseed << endmsg; + // ACTS_WARNING("Track finding failed for truth seed " << iseed << " with error" << + // result.error()); + // // Track finding failed, but still create an empty SimMultiTrajectory + // // trajectories->push_back(SimMultiTrajectory()); + //} + if (result.ok()) + { + // Get the fit output object + const auto& fitOutput = result.value(); + // 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. + Trajectories::IndexedParameters indexedParams; + //if (fitOutput.fittedParameters) { + // const auto& params = fitOutput.fittedParameters.value(); + // ACTS_VERBOSE("Fitted paramemeters for track " << itrack); + // ACTS_VERBOSE(" " << params.parameters().transpose()); + // // Push the fitted parameters to the container + // indexedParams.emplace(fitOutput.lastMeasurementIndex, std::move(params)); + //} else { + // ACTS_DEBUG("No fitted paramemeters for track " << itrack); + //} + // store the result + trajectories->emplace_back(std::move(fitOutput.fittedStates), std::move(trackTips), + std::move(indexedParams)); } else { - debug() << "Track finding failed for truth seed " << iseed << endmsg; - ACTS_WARNING("Track finding failed for truth seed " << iseed << " with error" << result.error()); - // Track finding failed, but still create an empty SimMultiTrajectory - // trajectories->push_back(SimMultiTrajectory()); + ACTS_WARNING("Fit failed for track " << iseed << " with error" << result.error()); + // Fit failed. Add an empty result so the output container has + // the same number of entries as the input. + trajectories->push_back(Trajectories()); } } // ctx.eventStore.add(m_cfg.outputTrajectories, std::move(trajectories)); return StatusCode::SUCCESS; - /////////////////////////// // acts example diff --git a/JugTrack/src/components/TrackFittingAlgorithm.h b/JugTrack/src/components/TrackFittingAlgorithm.h index 4322fe5a92606cf060f7c914d5d008ddca72ab0a..567dc4f408bd30631edb877f7a9a4fc886f940bc 100644 --- a/JugTrack/src/components/TrackFittingAlgorithm.h +++ b/JugTrack/src/components/TrackFittingAlgorithm.h @@ -24,15 +24,19 @@ #include <stdexcept> #include <vector> -#include "JugTrack/SourceLinks.h" +#include "JugTrack/IndexSourceLink.hpp" #include "JugTrack/Track.hpp" #include "JugTrack/BField.h" #include "JugTrack/Measurement.hpp" +#include "JugTrack/Trajectories.hpp" #include "eicd/TrackerHitCollection.h" //#include "Acts/Surfaces/PerigeeSurface.hpp" +#include "Acts/Definitions/Common.hpp" +#include "Acts/Geometry/TrackingGeometry.hpp" +#include "Acts/MagneticField/MagneticFieldProvider.hpp" #include "Acts/TrackFitting/KalmanFitter.hpp" #include "Acts/TrackFitting/GainMatrixSmoother.hpp" #include "Acts/TrackFitting/GainMatrixUpdater.hpp" @@ -56,32 +60,26 @@ namespace Jug::Reco { class TrackFittingAlgorithm : public GaudiAlgorithm { public: + /// Track fitter function that takes input measurements, initial trackstate + /// and fitter options and returns some track-fitter-specific result. + using TrackFitterOptions = + Acts::KalmanFitterOptions<MeasurementCalibrator, Acts::VoidOutlierFinder>; + //using TrackFinderResult = Acts::Result<Acts::CombinatorialKalmanFilterResult<SourceLink>>; - using FitterResult = Acts::Result<Acts::KalmanFitterResult<SourceLink>>; - /// Fit function that takes input measurements, initial trackstate and fitter + using FitterResult = Acts::Result<Acts::KalmanFitterResult<IndexSourceLink>>; + /// Fit function that takes input measurements, initial trackstate and fitter using FitterFunction = std::function<FitterResult( - const std::vector<SourceLink>&, const TrackParameters&, - const Acts::KalmanFitterOptions<Acts::VoidOutlierFinder>&)>; - - - /// Track fitter function that takes input measurements, initial trackstate - /// and fitter options and returns some track-fitter-specific result. - using TrackFitterOptions = - Acts::KalmanFitterOptions< Acts::VoidOutlierFinder>; - using TrackFitterResult = - Acts::Result<Acts::KalmanFitterResult<SourceLink>>; - using TrackFitterFunction = std::function<TrackFitterResult( - const std::vector<SourceLink>&, const TrackParameters&, - const TrackFitterOptions&)>; + const std::vector<IndexSourceLink>&, const TrackParameters&, const TrackFitterOptions&)>; public: - DataHandle<SourceLinkContainer> m_inputSourceLinks{"inputSourceLinks", Gaudi::DataHandle::Reader, this}; + DataHandle<IndexSourceLinkContainer> m_inputSourceLinks{"inputSourceLinks", Gaudi::DataHandle::Reader, this}; DataHandle<TrackParametersContainer> m_initialTrackParameters{"initialTrackParameters", Gaudi::DataHandle::Reader, this}; - DataHandle<TrajectoryContainer> m_foundTracks{"foundTracks", Gaudi::DataHandle::Reader, this}; - DataHandle<TrajectoryContainer> m_outputTrajectories{"outputTrajectories", Gaudi::DataHandle::Writer, this}; + DataHandle<MeasurementContainer> m_inputMeasurements{"inputMeasurements", Gaudi::DataHandle::Reader, this}; + DataHandle<TrajectoriesContainer> m_foundTracks{"foundTracks", Gaudi::DataHandle::Reader, this}; + DataHandle<TrajectoriesContainer> m_outputTrajectories{"outputTrajectories", Gaudi::DataHandle::Writer, this}; - FitterFunction m_trackFittingFunc; + FitterFunction m_trackFittingFunc; SmartIF<IGeoSvc> m_geoSvc; std::shared_ptr<Acts::ConstantBField> m_BField = nullptr; Acts::GeometryContext m_geoctx; @@ -96,15 +94,36 @@ namespace Jug::Reco { * The magnetic field is intentionally given by-value since the variant * contains shared_ptr anyways. */ - static FitterFunction makeTrackFittingFunction(std::shared_ptr<const Acts::TrackingGeometry> trackingGeometry, - BFieldVariant magneticField); - + static FitterFunction + makeTrackFittingFunction(std::shared_ptr<const Acts::TrackingGeometry> trackingGeometry, + std::shared_ptr<const Acts::MagneticFieldProvider> magneticField); + // BFieldVariant magneticField); StatusCode initialize() override; StatusCode execute() override; + private: + /// Helper function to call correct FitterFunction + FitterResult fitTrack( + const std::vector<IndexSourceLink>& sourceLinks, + const TrackParameters& initialParameters, + const TrackFitterOptions& options + ) const; + //, const std::vector<const Acts::Surface*>& surfSequence) const; }; + inline TrackFittingAlgorithm::FitterResult TrackFittingAlgorithm::fitTrack( + const std::vector<IndexSourceLink>& sourceLinks, const TrackParameters& initialParameters, + const Acts::KalmanFitterOptions<MeasurementCalibrator, Acts::VoidOutlierFinder>& options) + const + { + // const std::vector<const Acts::Surface*>& surfSequence) const + // if (m_cfg.directNavigation) { + // return m_cfg.dFit(sourceLinks, initialParameters, options, surfSequence); + //} + return m_trackFittingFunc(sourceLinks, initialParameters, options); + } + } // namespace Jug::Reco #endif diff --git a/JugTrack/src/components/TrackFittingFunction.cpp b/JugTrack/src/components/TrackFittingFunction.cpp index 9eafb46ffba71f01639e3a6368322266e9d02bc8..53abe7ac138c90c2a2a2e0303bf411fb7e749eae 100644 --- a/JugTrack/src/components/TrackFittingFunction.cpp +++ b/JugTrack/src/components/TrackFittingFunction.cpp @@ -10,7 +10,9 @@ #include "Acts/TrackFitting/GainMatrixSmoother.hpp" #include "Acts/TrackFitting/GainMatrixUpdater.hpp" #include "Acts/Utilities/Helpers.hpp" -#include "Acts/Utilities/ParameterDefinitions.hpp" +#include "Acts/Definitions/Common.hpp" +#include "Acts/Geometry/TrackingGeometry.hpp" +#include "Acts/MagneticField/MagneticFieldProvider.hpp" //#include "ActsExamples/Plugins/BField/ScalableBField.hpp" #include "TrackFittingAlgorithm.h" @@ -22,11 +24,12 @@ namespace { TrackFitterFunctionImpl(track_fitter_t&& f) : trackFitter(std::move(f)) {} - Jug::Reco::TrackFittingAlgorithm::FitterResult operator()( - const std::vector<Jug::SourceLink>& sourceLinks, - const Jug::TrackParameters& initialParameters, - const Acts::KalmanFitterOptions<Acts::VoidOutlierFinder>& options) const { - return fitter.fit(sourceLinks, initialParameters, options); + Jug::Reco::TrackFittingAlgorithm::FitterResult + operator()(const std::vector<Jug::IndexSourceLink>& sourceLinks, + const Jug::TrackParameters& initialParameters, + const Jug::Reco::TrackFittingAlgorithm::TrackFitterOptions& options) const + { + return trackFitter.fit(sourceLinks, initialParameters, options); } //TrackFittingAlgorithm::TrackFitterResult @@ -41,37 +44,63 @@ namespace { } // namespace namespace Jug::Reco { - TrackFittingAlgorithm::TrackFitterFunction TrackFittingAlgorithm::makeTrackFittingFunction( - std::shared_ptr<const Acts::TrackingGeometry> trackingGeometry, Options::BFieldVariant magneticField) - { - using Updater = Acts::GainMatrixUpdater; - using Smoother = Acts::GainMatrixSmoother; - // unpack the magnetic field variant and instantiate the corresponding fitter. - return std::visit( - [trackingGeometry](auto&& inputField) -> TrackFitterFunction { - // each entry in the variant is already a shared_ptr - // need ::element_type to get the real magnetic field type - using InputMagneticField = typename std::decay_t<decltype(inputField)>::element_type; - using MagneticField = Acts::SharedBField<InputMagneticField>; - using Stepper = Acts::EigenStepper<MagneticField>; - using Navigator = Acts::Navigator; - using Propagator = Acts::Propagator<Stepper, Navigator>; - using Fitter = Acts::KalmanFitter<Propagator, Updater, Smoother>; + //using InputMagneticField = typename std::decay_t<decltype(inputField)>::element_type; + //using MagneticField = Acts::SharedBField<InputMagneticField>; + //using Stepper = Acts::EigenStepper<MagneticField>; + //using Navigator = Acts::Navigator; + //using Propagator = Acts::Propagator<Stepper, Navigator>; + //using Fitter = Acts::KalmanFitter<Propagator, Updater, Smoother>; - // construct all components for the fitter - MagneticField field(std::move(inputField)); - Stepper stepper(std::move(field)); - Navigator navigator(trackingGeometry); - navigator.resolvePassive = false; - navigator.resolveMaterial = true; - navigator.resolveSensitive = true; - Propagator propagator(std::move(stepper), std::move(navigator)); - Fitter trackFitter(std::move(propagator)); + 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 DirectPropagator = Acts::Propagator<Stepper, Acts::DirectNavigator>; + using DirectFitter = Acts::KalmanFitter<DirectPropagator, Updater, Smoother>; - // build the fitter functions. owns the fitter object. - return TrackFitterFunctionImpl<Fitter>(std::move(trackFitter)); - }, - std::move(magneticField)); + TrackFittingAlgorithm::FitterFunction TrackFittingAlgorithm::makeTrackFittingFunction( + std::shared_ptr<const Acts::TrackingGeometry> trackingGeometry, + std::shared_ptr<const Acts::MagneticFieldProvider> magneticField) + + { + Stepper stepper(std::move(magneticField)); + //Acts::Navigator::Config cfg{trackingGeometry}; + //cfg.resolvePassive = false; + //cfg.resolveMaterial = true; + //cfg.resolveSensitive = true; + Acts::Navigator navigator(trackingGeometry); + navigator.resolvePassive = false; + navigator.resolveMaterial = true; + navigator.resolveSensitive = true; + Propagator propagator(std::move(stepper), std::move(navigator)); + Fitter trackFitter(std::move(propagator)); + + // build the fitter functions. owns the fitter object. + return TrackFitterFunctionImpl(std::move(trackFitter)); } + // using Updater = Acts::GainMatrixUpdater; + // using Smoother = Acts::GainMatrixSmoother; + + // // unpack the magnetic field variant and instantiate the corresponding fitter. + // return std::visit( + // [trackingGeometry](auto&& inputField) -> FitterFunction { + // // each entry in the variant is already a shared_ptr + // // need ::element_type to get the real magnetic field type + // // construct all components for the fitter + // MagneticField field(std::move(inputField)); + // Stepper stepper(std::move(field)); + // Navigator navigator(trackingGeometry); + // navigator.resolvePassive = false; + // navigator.resolveMaterial = true; + // navigator.resolveSensitive = true; + // Propagator propagator(std::move(stepper), std::move(navigator)); + // Fitter trackFitter(std::move(propagator)); + + // // build the fitter functions. owns the fitter object. + // return TrackFitterFunctionImpl<Fitter>(std::move(trackFitter)); + // }, + // std::move(magneticField)); + //} } // namespace Jug::Reco