diff --git a/CI/physmon/phys_perf_mon.sh b/CI/physmon/phys_perf_mon.sh index c2944f0dd06..8482f92c4b6 100755 --- a/CI/physmon/phys_perf_mon.sh +++ b/CI/physmon/phys_perf_mon.sh @@ -265,6 +265,15 @@ function trackfinding() { $path/performance_finding_ckf_ambi.html \ $path/performance_finding_ckf_ambi fi + + if [ -f $refdir/$path/performance_finding_ckf_ml_solver.root ]; then + run_histcmp \ + $outdir/data/$path/performance_finding_ckf_ml_solver.root \ + $refdir/$path/performance_finding_ckf_ml_solver.root \ + "ML Ambisolver | ${name}" \ + $path/performance_finding_ckf_ml_solver.html \ + $path/performance_finding_ckf_ml_solver + fi } function vertexing() { diff --git a/CI/physmon/reference/trackfinding_ttbar_pu200/performance_finding_ckf_ml_solver.root b/CI/physmon/reference/trackfinding_ttbar_pu200/performance_finding_ckf_ml_solver.root new file mode 100644 index 00000000000..db22a916d9e Binary files /dev/null and b/CI/physmon/reference/trackfinding_ttbar_pu200/performance_finding_ckf_ml_solver.root differ diff --git a/CI/physmon/reference/trackfinding_ttbar_pu200/performance_fitting_ckf_ml_solver.root b/CI/physmon/reference/trackfinding_ttbar_pu200/performance_fitting_ckf_ml_solver.root new file mode 100644 index 00000000000..0a40506d1bb Binary files /dev/null and b/CI/physmon/reference/trackfinding_ttbar_pu200/performance_fitting_ckf_ml_solver.root differ diff --git a/CI/physmon/workflows/physmon_trackfinding_ttbar_pu200.py b/CI/physmon/workflows/physmon_trackfinding_ttbar_pu200.py index 99f12d0170d..4b4c75e1f36 100755 --- a/CI/physmon/workflows/physmon_trackfinding_ttbar_pu200.py +++ b/CI/physmon/workflows/physmon_trackfinding_ttbar_pu200.py @@ -19,7 +19,9 @@ CkfConfig, addCKFTracks, addAmbiguityResolution, + addAmbiguityResolutionML, AmbiguityResolutionConfig, + AmbiguityResolutionMLConfig, addVertexFitting, VertexFinder, TrackSelectorConfig, @@ -134,6 +136,17 @@ outputDirRoot=tp, ) + addAmbiguityResolutionML( + s, + AmbiguityResolutionMLConfig( + maximumSharedHits=3, maximumIterations=1000000, nMeasurementsMin=6 + ), + tracks="ckf_tracks", + outputDirRoot=tp, + onnxModelFile=Path(__file__).resolve().parent.parent.parent.parent + / "thirdparty/OpenDataDetector/data/duplicateClassifier.onnx", + ) + addAmbiguityResolution( s, AmbiguityResolutionConfig( @@ -141,6 +154,7 @@ maximumIterations=100000, nMeasurementsMin=6, ), + tracks="ckf_tracks", outputDirRoot=tp, ) @@ -187,6 +201,17 @@ tp / "performance_fitting_ambi.root", tp / "performance_fitting_ckf_ambi.root", ) + + shutil.move( + tp / "performance_finding_ambiML.root", + tp / "performance_finding_ckf_ml_solver.root", + ) + + shutil.move( + tp / "performance_fitting_ambiML.root", + tp / "performance_fitting_ckf_ml_solver.root", + ) + for vertexing in ["amvf_gauss_notime", "amvf_grid_time"]: shutil.move( tp / f"{vertexing}/performance_vertexing.root", @@ -200,6 +225,8 @@ "performance_fitting_ckf.root", "performance_finding_ckf_ambi.root", "performance_fitting_ckf_ambi.root", + "performance_finding_ckf_ml_solver.root", + "performance_fitting_ckf_ml_solver.root", "performance_vertexing_amvf_gauss_notime.root", "performance_vertexing_amvf_grid_time.root", ]: diff --git a/Core/include/Acts/AmbiguityResolution/AmbiguityNetworkConcept.hpp b/Core/include/Acts/AmbiguityResolution/AmbiguityNetworkConcept.hpp new file mode 100644 index 00000000000..a69e1fee9fe --- /dev/null +++ b/Core/include/Acts/AmbiguityResolution/AmbiguityNetworkConcept.hpp @@ -0,0 +1,51 @@ +// This file is part of the ACTS project. +// +// Copyright (C) 2016 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 https://mozilla.org/MPL/2.0/. + +#pragma once + +#include "Acts/EventData/TrackContainer.hpp" +#include "Acts/EventData/TrackContainerFrontendConcept.hpp" +#include "Acts/EventData/VectorMultiTrajectory.hpp" +#include "Acts/EventData/VectorTrackContainer.hpp" +#include "Acts/Utilities/Concepts.hpp" + +namespace Acts { + +/// @brief Concept for the ambiguity network used in the ambiguity resolution +/// +/// The ambiguity network correspond to the AmbiguityTrackClassifier found in +/// the Onnx plugin. It is used to score the tracks and select the best ones. +/// +/// The constructor of the Ambiguity Solver network should take string as input +/// corresponding to the path of the ONNX model. +/// The implementation of the Ambiguity Solver network should have two methods: +/// - inferScores: takes clusters (a list of track ID associated with a cluster +/// ID) and the track container and return an outputTensor (list of scores for +/// each track in the clusters). +/// - trackSelection: Takes clusters and the output tensor from the inferScores +/// method and return the list of track ID to keep. +/// +/// @tparam N the type of the network +template +concept AmbiguityNetworkConcept = requires( + TrackContainer &tracks, + std::unordered_map> &clusters, + std::vector> &outputTensor, const char *modelPath, + network_t &n) { + { network_t(modelPath) } -> std::same_as; + + { + n.inferScores(clusters, tracks) + } -> std::same_as>>; + { + n.trackSelection(clusters, outputTensor) + } -> std::same_as>; +}; + +} // namespace Acts diff --git a/Core/include/Acts/AmbiguityResolution/AmbiguityResolutionML.hpp b/Core/include/Acts/AmbiguityResolution/AmbiguityResolutionML.hpp new file mode 100644 index 00000000000..66717e3f8ee --- /dev/null +++ b/Core/include/Acts/AmbiguityResolution/AmbiguityResolutionML.hpp @@ -0,0 +1,136 @@ +// This file is part of the ACTS project. +// +// Copyright (C) 2016 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 https://mozilla.org/MPL/2.0/. + +#pragma once + +#include "Acts/AmbiguityResolution/AmbiguityNetworkConcept.hpp" +#include "Acts/Definitions/Units.hpp" +#include "Acts/EventData/TrackContainer.hpp" +#include "Acts/Utilities/Delegate.hpp" +#include "Acts/Utilities/Logger.hpp" + +#include +#include +#include +#include +#include +#include + +namespace Acts { + +/// Generic implementation of the machine learning ambiguity resolution +/// Contains method for data preparations +template +class AmbiguityResolutionML { + public: + struct Config { + /// Path to the model file for the duplicate neural network + std::string inputDuplicateNN = ""; + /// Minimum number of measurement to form a track. + std::size_t nMeasurementsMin = 7; + }; + /// Construct the ambiguity resolution algorithm. + /// + /// @param cfg is the algorithm configuration + /// @param logger is the logging instance + AmbiguityResolutionML(const Config& cfg, + std::unique_ptr logger = getDefaultLogger( + "AmbiguityResolutionML", Logging::INFO)) + : m_cfg{cfg}, + m_duplicateClassifier(m_cfg.inputDuplicateNN.c_str()), + m_logger{std::move(logger)} {} + + /// Associate the hits to the tracks + /// + /// This algorithm performs the mapping of hits ID to track ID. Our final goal + /// is too loop over all the tracks (and their associated hits) by order of + /// decreasing number hits for this we use a multimap where the key is the + /// number of hits as this will automatically perform the sorting. + /// + /// @param tracks is the input track container + /// @param sourceLinkHash is the hash function for the source link, will be used to associate to tracks + /// @param sourceLinkEquality is the equality function for the source link used used to associated hits to tracks + /// @return an ordered list containing pairs of track ID and associated measurement ID + template + std::multimap>> + mapTrackHits(const track_container_t& tracks, + const source_link_hash_t& sourceLinkHash, + const source_link_equality_t& sourceLinkEquality) const { + // A map to store (and generate) the measurement index for each source link + auto measurementIndexMap = + std::unordered_map(0, sourceLinkHash, + sourceLinkEquality); + + // A map to store the track Id and their associated measurements ID, a + // multimap is used to automatically sort the tracks by the number of + // measurements + std::multimap>> + trackMap; + std::size_t trackIndex = 0; + std::vector measurements; + // Loop over all the trajectories in the events + for (const auto& track : tracks) { + // Kick out tracks that do not fulfill our initial requirements + if (track.nMeasurements() < m_cfg.nMeasurementsMin) { + continue; + } + measurements.clear(); + for (auto ts : track.trackStatesReversed()) { + if (ts.typeFlags().test(Acts::TrackStateFlag::MeasurementFlag)) { + SourceLink sourceLink = ts.getUncalibratedSourceLink(); + // assign a new measurement index if the source link was not seen yet + auto emplace = measurementIndexMap.try_emplace( + sourceLink, measurementIndexMap.size()); + measurements.push_back(emplace.first->second); + } + } + trackMap.emplace(track.nMeasurements(), + std::make_pair(trackIndex, measurements)); + ++trackIndex; + } + return trackMap; + } + + /// Select the track associated with each cluster + /// + /// In this algorithm the call the neural network to score the tracks and then + /// select the track with the highest score in each cluster + /// + /// @param clusters is a map of clusters, each cluster correspond to a vector of track ID + /// @param tracks is the input track container + /// @return a vector of trackID corresponding tho the good tracks + template + std::vector solveAmbiguity( + std::unordered_map>& clusters, + const track_container_t& tracks) const { + std::vector> outputTensor = + m_duplicateClassifier.inferScores(clusters, tracks); + std::vector goodTracks = + m_duplicateClassifier.trackSelection(clusters, outputTensor); + + return goodTracks; + } + + private: + // Configuration + Config m_cfg; + + // The neural network for duplicate classification, the network + // implementation is chosen with the AmbiguityNetwork template parameter + AmbiguityNetwork m_duplicateClassifier; + + /// Logging instance + std::unique_ptr m_logger = nullptr; + + /// Private access to logging instance + const Logger& logger() const { return *m_logger; } +}; + +} // namespace Acts diff --git a/Core/include/Acts/EventData/MultiComponentTrackParameters.hpp b/Core/include/Acts/EventData/MultiComponentTrackParameters.hpp index cf932fc0470..5e079af0df7 100644 --- a/Core/include/Acts/EventData/MultiComponentTrackParameters.hpp +++ b/Core/include/Acts/EventData/MultiComponentTrackParameters.hpp @@ -129,11 +129,11 @@ class MultiComponentBoundTrackParameters { /// Get the weight and a GenericBoundTrackParameters object for one component std::pair operator[](std::size_t i) const { - return std::make_pair( + return { std::get(m_components[i]), Parameters(m_surface, std::get(m_components[i]), std::get>(m_components[i]), - m_particleHypothesis)); + m_particleHypothesis)}; } /// Parameters vector. diff --git a/Core/include/Acts/MagneticField/MultiRangeBField.hpp b/Core/include/Acts/MagneticField/MultiRangeBField.hpp new file mode 100644 index 00000000000..60b90097161 --- /dev/null +++ b/Core/include/Acts/MagneticField/MultiRangeBField.hpp @@ -0,0 +1,67 @@ +// This file is part of the ACTS project. +// +// Copyright (C) 2016 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 https://mozilla.org/MPL/2.0/. + +#pragma once +#include "Acts/Definitions/Algebra.hpp" +#include "Acts/MagneticField/MagneticFieldContext.hpp" +#include "Acts/MagneticField/MagneticFieldError.hpp" +#include "Acts/MagneticField/MagneticFieldProvider.hpp" +#include "Acts/Utilities/RangeXD.hpp" + +namespace Acts { + +/// @ingroup MagneticField +/// +/// @brief Magnetic field provider modelling a magnetic field consisting of +/// several (potentially overlapping) regions of constant values. +class MultiRangeBField final : public MagneticFieldProvider { + private: + struct Cache { + explicit Cache(const MagneticFieldContext& /*unused*/); + + std::optional index = {}; + }; + + using BFieldRange = std::pair, Vector3>; + + // The different ranges and their corresponding field vectors. Note that + // regions positioned _later_ in this vector take priority over earlier + // regions. + std::vector fieldRanges; + + public: + /// @brief Construct a magnetic field from a vector of ranges. + /// + /// @warning These ranges are listed in increasing order of precedence, + /// i.e. ranges further along the vector have higher priority. + explicit MultiRangeBField(const std::vector& ranges); + + explicit MultiRangeBField(std::vector&& ranges); + + /// @brief Construct a cache object. + MagneticFieldProvider::Cache makeCache( + const MagneticFieldContext& mctx) const override; + + /// @brief Request the value of the magnetic field at a given position. + /// + /// @param [in] position Global 3D position for the lookup. + /// @param [in, out] cache Cache object. + /// @returns A successful value containing a field vector if the given + /// location is contained inside any of the regions, or a failure value + /// otherwise. + Result getField(const Vector3& position, + MagneticFieldProvider::Cache& cache) const override; + + /// @brief Get the field gradient at a given position. + /// + /// @warning This is not currently implemented. + Result getFieldGradient( + const Vector3& position, ActsMatrix<3, 3>& /*unused*/, + MagneticFieldProvider::Cache& cache) const override; +}; +} // namespace Acts diff --git a/Core/include/Acts/Seeding/BinnedGroupIterator.ipp b/Core/include/Acts/Seeding/BinnedGroupIterator.ipp index 70cd691f5db..8bcd0ac04b8 100644 --- a/Core/include/Acts/Seeding/BinnedGroupIterator.ipp +++ b/Core/include/Acts/Seeding/BinnedGroupIterator.ipp @@ -64,7 +64,7 @@ Acts::BinnedGroupIterator::operator*() const { #pragma GCC diagnostic push #pragma GCC diagnostic ignored "-Wstringop-overread" #endif - return std::make_tuple(std::move(bottoms), global_index, std::move(tops)); + return {std::move(bottoms), global_index, std::move(tops)}; #if defined(__GNUC__) && __GNUC__ >= 12 && !defined(__clang__) #pragma GCC diagnostic pop #endif diff --git a/Core/include/Acts/Seeding/SeedFinder.ipp b/Core/include/Acts/Seeding/SeedFinder.ipp index e6ccb9f589f..b61947bb634 100644 --- a/Core/include/Acts/Seeding/SeedFinder.ipp +++ b/Core/include/Acts/Seeding/SeedFinder.ipp @@ -839,7 +839,7 @@ std::pair SeedFinder:: const external_spacepoint_t& spM, const Acts::Range1D& rMiddleSPRange) const { if (m_config.useVariableMiddleSPRange) { - return std::make_pair(rMiddleSPRange.min(), rMiddleSPRange.max()); + return {rMiddleSPRange.min(), rMiddleSPRange.max()}; } if (!m_config.rRangeMiddleSP.empty()) { /// get zBin position of the middle SP @@ -848,10 +848,9 @@ std::pair SeedFinder:: int zBin = std::distance(m_config.zBinEdges.begin(), pVal); /// protects against zM at the limit of zBinEdges zBin == 0 ? zBin : --zBin; - return std::make_pair(m_config.rRangeMiddleSP[zBin][0], - m_config.rRangeMiddleSP[zBin][1]); + return {m_config.rRangeMiddleSP[zBin][0], m_config.rRangeMiddleSP[zBin][1]}; } - return std::make_pair(m_config.rMinMiddle, m_config.rMaxMiddle); + return {m_config.rMinMiddle, m_config.rMaxMiddle}; } } // namespace Acts diff --git a/Core/include/Acts/TrackFinding/detail/AmbiguityTrackClustering.hpp b/Core/include/Acts/TrackFinding/detail/AmbiguityTrackClustering.hpp index 6bc565f80eb..644a0fa34fd 100644 --- a/Core/include/Acts/TrackFinding/detail/AmbiguityTrackClustering.hpp +++ b/Core/include/Acts/TrackFinding/detail/AmbiguityTrackClustering.hpp @@ -15,7 +15,15 @@ namespace Acts::detail { -/// Clusterise tracks based on shared hits +/// Cluster tracks based on shared hits. +/// +/// In this algorithm we will loop through all the tracks by decreasing number +/// of measurements. Cluster are created when a new track is encountered that +/// doesn't share hits with the leading track of a previous cluster (with the +/// leading track defined as the track that lead to the cluster creation). If a +/// track shares hits with the leading track of a cluster, it is added to that +/// cluster. If a track shares hits with multiple clusters, it is associated to +/// the cluster with the leading track with the most hits. /// /// @param trackMap : Multimap storing pair of track ID and vector of measurement ID. The keys are the number of measurement and are just there to facilitate the ordering. /// @return an unordered map representing the clusters, the keys the ID of the primary track of each cluster and the store a vector of track IDs. diff --git a/Core/include/Acts/TrackFitting/GlobalChiSquareFitter.hpp b/Core/include/Acts/TrackFitting/GlobalChiSquareFitter.hpp index dc0bc221e9a..5e46bc0dba3 100644 --- a/Core/include/Acts/TrackFitting/GlobalChiSquareFitter.hpp +++ b/Core/include/Acts/TrackFitting/GlobalChiSquareFitter.hpp @@ -361,6 +361,29 @@ struct Gx2fSystem { std::size_t m_ndf = 0u; }; +/// @brief Adds a measurement to the GX2F equation system in a modular backend function. +/// +/// This function processes measurement data and integrates it into the GX2F +/// system. +/// +/// @param extendedSystem All parameters of the current equation system to update. +/// @param jacobianFromStart The Jacobian matrix from the start to the current state. +/// @param covarianceMeasurement The covariance matrix of the measurement. +/// @param predicted The predicted state vector based on the track state. +/// @param measurement The measurement vector. +/// @param projector The projection matrix. +/// @param logger A logger instance. +/// +/// @note The dynamic Eigen matrices are suboptimal. We could think of +/// templating again in the future on kMeasDims. We currently use dynamic +/// matrices to reduce the memory during compile time. +void addMeasurementToGx2fSumsBackend( + Gx2fSystem& extendedSystem, + const std::vector& jacobianFromStart, + const Eigen::MatrixXd& covarianceMeasurement, const BoundVector& predicted, + const Eigen::VectorXd& measurement, const Eigen::MatrixXd& projector, + const Logger& logger); + /// @brief Process measurements and fill the aMatrix and bVector /// /// The function processes each measurement for the GX2F Actor fitting process. @@ -370,7 +393,7 @@ struct Gx2fSystem { /// @tparam kMeasDim Number of dimensions of the measurement /// @tparam track_state_t The type of the track state /// -/// @param extendedSystem All parameters of the current equation system +/// @param extendedSystem All parameters of the current equation system to update /// @param jacobianFromStart The Jacobian matrix from start to the current state /// @param trackState The track state to analyse /// @param logger A logger instance @@ -379,44 +402,9 @@ void addMeasurementToGx2fSums(Gx2fSystem& extendedSystem, const std::vector& jacobianFromStart, const track_state_t& trackState, const Logger& logger) { - // First we get back the covariance and try to invert it. If the inversion - // fails, we can already abort. const ActsSquareMatrix covarianceMeasurement = trackState.template calibratedCovariance(); - const auto safeInvCovMeasurement = safeInverse(covarianceMeasurement); - if (!safeInvCovMeasurement) { - ACTS_WARNING("addMeasurementToGx2fSums: safeInvCovMeasurement failed."); - ACTS_VERBOSE(" covarianceMeasurement:\n" << covarianceMeasurement); - return; - } - - // Create an extended Jacobian. This one contains only eBoundSize rows, - // because the rest is irrelevant. We fill it in the next steps. - // TODO make dimsExtendedParams template with unrolling - Eigen::MatrixXd extendedJacobian = - Eigen::MatrixXd::Zero(eBoundSize, extendedSystem.nDims()); - - // This part of the Jacobian comes from the material-less propagation - extendedJacobian.topLeftCorner() = - jacobianFromStart[0]; - - // If we have material, loop here over all Jacobians. We add extra columns for - // their phi-theta projections. These parts account for the propagation of the - // scattering angles. - for (std::size_t matSurface = 1; matSurface < jacobianFromStart.size(); - matSurface++) { - const BoundMatrix jac = jacobianFromStart[matSurface]; - - const ActsMatrix jacPhiTheta = - jac * Gx2fConstants::phiThetaProjector; - - // The position, where we need to insert the values in the extended Jacobian - const std::size_t deltaPosition = eBoundSize + 2 * (matSurface - 1); - - extendedJacobian.block(0, deltaPosition) = jacPhiTheta; - } - const BoundVector predicted = trackState.smoothed(); const ActsVector measurement = @@ -425,54 +413,9 @@ void addMeasurementToGx2fSums(Gx2fSystem& extendedSystem, const ActsMatrix projector = trackState.template projectorSubspaceHelper().projector(); - const Eigen::MatrixXd projJacobian = projector * extendedJacobian; - - const ActsMatrix projPredicted = projector * predicted; - - const ActsVector residual = measurement - projPredicted; - - // Finally contribute to chi2sum, aMatrix, and bVector - extendedSystem.chi2() += - (residual.transpose() * (*safeInvCovMeasurement) * residual)(0, 0); - - extendedSystem.aMatrix() += - (projJacobian.transpose() * (*safeInvCovMeasurement) * projJacobian) - .eval(); - - extendedSystem.bVector() += - (residual.transpose() * (*safeInvCovMeasurement) * projJacobian) - .eval() - .transpose(); - - ACTS_VERBOSE( - "Contributions in addMeasurementToGx2fSums:\n" - << " kMeasDim: " << kMeasDim << "\n" - << " predicted: " << predicted.transpose() << "\n" - << " measurement: " << measurement.transpose() << "\n" - << " covarianceMeasurement:\n" - << covarianceMeasurement << "\n" - << " projector:\n" - << projector.eval() << "\n" - << " projJacobian:\n" - << projJacobian.eval() << "\n" - << " projPredicted: " << (projPredicted.transpose()).eval() << "\n" - << " residual: " << (residual.transpose()).eval() << "\n" - << " extendedJacobian:\n" - << extendedJacobian << "\n" - << " aMatrix contribution:\n" - << (projJacobian.transpose() * (*safeInvCovMeasurement) * projJacobian) - .eval() - << "\n" - << " bVector contribution: " - << (residual.transpose() * (*safeInvCovMeasurement) * projJacobian).eval() - << "\n" - << " chi2sum contribution: " - << (residual.transpose() * (*safeInvCovMeasurement) * residual)(0, 0) - << "\n" - << " safeInvCovMeasurement:\n" - << (*safeInvCovMeasurement)); - - return; + addMeasurementToGx2fSumsBackend(extendedSystem, jacobianFromStart, + covarianceMeasurement, predicted, measurement, + projector, logger); } /// @brief Process material and fill the aMatrix and bVector @@ -694,6 +637,15 @@ std::size_t countMaterialStates( return nMaterialSurfaces; } +/// @brief Solve the gx2f system to get the delta parameters for the update +/// +/// This function computes the delta parameters for the GX2F Actor fitting +/// process by solving the linear equation system [a] * delta = b. It uses the +/// column-pivoting Householder QR decomposition for numerical stability. +/// +/// @param extendedSystem All parameters of the current equation system +Eigen::VectorXd computeGx2fDeltaParams(const Gx2fSystem& extendedSystem); + /// @brief Update parameters (and scattering angles if applicable) /// /// @param params Parameters to be updated @@ -1391,10 +1343,8 @@ class Gx2Fitter { return Experimental::GlobalChiSquareFitterError::NotEnoughMeasurements; } - // calculate delta params [a] * delta = b Eigen::VectorXd deltaParamsExtended = - extendedSystem.aMatrix().colPivHouseholderQr().solve( - extendedSystem.bVector()); + computeGx2fDeltaParams(extendedSystem); ACTS_VERBOSE("aMatrix:\n" << extendedSystem.aMatrix() << "\n" @@ -1558,10 +1508,8 @@ class Gx2Fitter { return Experimental::GlobalChiSquareFitterError::NotEnoughMeasurements; } - // calculate delta params [a] * delta = b Eigen::VectorXd deltaParamsExtended = - extendedSystem.aMatrix().colPivHouseholderQr().solve( - extendedSystem.bVector()); + computeGx2fDeltaParams(extendedSystem); ACTS_VERBOSE("aMatrix:\n" << extendedSystem.aMatrix() << "\n" diff --git a/Core/include/Acts/TrackFitting/detail/KalmanGlobalCovariance.hpp b/Core/include/Acts/TrackFitting/detail/KalmanGlobalCovariance.hpp index cf9a3f36b01..66dfe439fd8 100644 --- a/Core/include/Acts/TrackFitting/detail/KalmanGlobalCovariance.hpp +++ b/Core/include/Acts/TrackFitting/detail/KalmanGlobalCovariance.hpp @@ -95,7 +95,7 @@ globalTrackParametersCovariance(const traj_t& multiTraj, prev_ts = ts; }); - return std::make_pair(fullGlobalTrackParamsCov, stateRowIndices); + return {fullGlobalTrackParamsCov, stateRowIndices}; } } // namespace Acts::detail diff --git a/Core/include/Acts/Utilities/AlgebraHelpers.hpp b/Core/include/Acts/Utilities/AlgebraHelpers.hpp index 06379ec3559..aa522a7dbc2 100644 --- a/Core/include/Acts/Utilities/AlgebraHelpers.hpp +++ b/Core/include/Acts/Utilities/AlgebraHelpers.hpp @@ -178,7 +178,7 @@ inline ActsMatrix blockedMult( /// Calculate the inverse of an Eigen matrix after checking if it can be /// numerically inverted. This allows to catch potential FPEs before they occur. /// For matrices up to 4x4, the inverse is computed directly. For larger -/// matrices, the FullPivLU is used. +/// matrices, and dynamic matrices the FullPivLU is used. /// /// @tparam Derived Eigen derived concrete type /// @tparam Result Eigen result type defaulted to input type @@ -192,12 +192,11 @@ std::optional safeInverse(const MatrixType& m) noexcept { constexpr int cols = MatrixType::ColsAtCompileTime; static_assert(rows == cols); - static_assert(rows != -1); ResultType result; bool invertible = false; - if constexpr (rows > 4) { + if constexpr (rows > 4 || rows == -1) { Eigen::FullPivLU mFullPivLU(m); if (mFullPivLU.isInvertible()) { invertible = true; diff --git a/Core/include/Acts/Utilities/GridBinFinder.ipp b/Core/include/Acts/Utilities/GridBinFinder.ipp index 8b06748ffec..37db6082bb0 100644 --- a/Core/include/Acts/Utilities/GridBinFinder.ipp +++ b/Core/include/Acts/Utilities/GridBinFinder.ipp @@ -45,9 +45,9 @@ std::array, DIM> Acts::GridBinFinder::getSizePerAxis( using value_t = typename std::decay_t; if constexpr (std::is_same_v) { assert(val >= 0); - return std::make_pair(-val, val); + return {-val, val}; } else if constexpr (std::is_same_v, value_t>) { - return std::make_pair(-val.first, val.second); + return {-val.first, val.second}; } else { assert(locPosition.size() > i); assert(locPosition[i] > 0ul); diff --git a/Core/include/Acts/Utilities/TrackHelpers.hpp b/Core/include/Acts/Utilities/TrackHelpers.hpp index 36d37fab714..dc3765be7f6 100644 --- a/Core/include/Acts/Utilities/TrackHelpers.hpp +++ b/Core/include/Acts/Utilities/TrackHelpers.hpp @@ -209,7 +209,7 @@ findTrackStateForExtrapolation( } ACTS_VERBOSE("found intersection at " << intersection.pathLength()); - return std::make_pair(*first, intersection.pathLength()); + return std::pair(*first, intersection.pathLength()); } case TrackExtrapolationStrategy::last: { @@ -229,7 +229,7 @@ findTrackStateForExtrapolation( } ACTS_VERBOSE("found intersection at " << intersection.pathLength()); - return std::make_pair(*last, intersection.pathLength()); + return std::pair(*last, intersection.pathLength()); } case TrackExtrapolationStrategy::firstOrLast: { @@ -256,13 +256,13 @@ findTrackStateForExtrapolation( if (intersectionFirst.isValid() && absDistanceFirst <= absDistanceLast) { ACTS_VERBOSE("using first track state with intersection at " << intersectionFirst.pathLength()); - return std::make_pair(*first, intersectionFirst.pathLength()); + return std::pair(*first, intersectionFirst.pathLength()); } if (intersectionLast.isValid() && absDistanceLast <= absDistanceFirst) { ACTS_VERBOSE("using last track state with intersection at " << intersectionLast.pathLength()); - return std::make_pair(*last, intersectionLast.pathLength()); + return std::pair(*last, intersectionLast.pathLength()); } ACTS_ERROR("no intersection found"); @@ -531,7 +531,7 @@ calculatePredictedResidual(track_state_proxy_t trackState) { MeasurementMatrix residualCovariance = measurementCovariance + predictedCovariance; - return std::pair(residual, residualCovariance); + return {residual, residualCovariance}; } /// Helper function to calculate the filtered residual and its covariance @@ -568,7 +568,7 @@ calculateFilteredResidual(track_state_proxy_t trackState) { MeasurementMatrix residualCovariance = measurementCovariance + filteredCovariance; - return std::pair(residual, residualCovariance); + return {residual, residualCovariance}; } /// Helper function to calculate the smoothed residual and its covariance @@ -605,7 +605,7 @@ calculateSmoothedResidual(track_state_proxy_t trackState) { MeasurementMatrix residualCovariance = measurementCovariance + smoothedCovariance; - return std::pair(residual, residualCovariance); + return {residual, residualCovariance}; } /// Helper function to calculate the predicted chi2 diff --git a/Core/include/Acts/Visualization/Interpolation3D.hpp b/Core/include/Acts/Visualization/Interpolation3D.hpp new file mode 100644 index 00000000000..ea8e649869e --- /dev/null +++ b/Core/include/Acts/Visualization/Interpolation3D.hpp @@ -0,0 +1,96 @@ +// This file is part of the ACTS project. +// +// Copyright (C) 2016 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 https://mozilla.org/MPL/2.0/. + +#pragma once + +#include "Acts/Definitions/Algebra.hpp" + +#include + +namespace Acts::Interpolation3D { + +/// @brief Helper function to interpolate points using a spline +/// from Eigen +/// +/// The only requirement is that the input trajectory type has +/// a method empty() and size() and that the elements can be +/// accessed with operator[] and have themselves a operator[] to +/// access the coordinates. +/// +/// @tparam input_trajectory_type input trajectory type +/// +/// @param inputsRaw input vector points +/// @param nPoints number of interpolation points +/// @param keepOriginalHits keep the original hits in the trajectory +/// +/// @return std::vector interpolated points +template +trajectory_type spline(const trajectory_type& inputsRaw, std::size_t nPoints, + bool keepOriginalHits = false) { + trajectory_type output; + if (inputsRaw.empty()) { + return output; + } + + using InputVectorType = typename trajectory_type::value_type; + + std::vector inputs; + // If input type is a vector of Vector3 we can use it directly + if constexpr (std::is_same_v>) { + inputs = inputsRaw; + } else { + inputs.reserve(inputsRaw.size()); + for (const auto& input : inputsRaw) { + inputs.push_back(Vector3(input[0], input[1], input[2])); + } + } + + // Don't do anything if we have less than 3 points or less interpolation + // points than input points + if (inputsRaw.size() < 3 || nPoints <= inputsRaw.size()) { + return inputsRaw; + } else { + Eigen::MatrixXd points(3, inputs.size()); + for (std::size_t i = 0; i < inputs.size(); ++i) { + points.col(i) = inputs[i].transpose(); + } + Eigen::Spline spline3D = + Eigen::SplineFitting>::Interpolate(points, 2); + + double step = 1. / (nPoints - 1); + for (std::size_t i = 0; i < nPoints; ++i) { + double t = i * step; + InputVectorType point; + point[0] = spline3D(t)[0]; + point[1] = spline3D(t)[1]; + point[2] = spline3D(t)[2]; + output.push_back(point); + } + } + // If we want to keep the original hits, we add them to the output + // (first and last are there anyway) + if (keepOriginalHits) { + output.insert(output.begin(), inputsRaw.begin() + 1, inputsRaw.end() - 1); + // We need to sort the output in distance to first + std::sort(output.begin(), output.end(), + [&inputs](const auto& a, const auto& b) { + const auto ifront = inputs.front(); + double da2 = (a[0] - ifront[0]) * (a[0] - ifront[0]) + + (a[1] - ifront[1]) * (a[1] - ifront[1]) + + (a[2] - ifront[2]) * (a[2] - ifront[2]); + double db2 = (b[0] - ifront[0]) * (b[0] - ifront[0]) + + (b[1] - ifront[1]) * (b[1] - ifront[1]) + + (b[2] - ifront[2]) * (b[2] - ifront[2]); + return da2 < db2; + }); + } + + return output; +} + +} // namespace Acts::Interpolation3D diff --git a/Core/src/Detector/detail/IndexedGridFiller.cpp b/Core/src/Detector/detail/IndexedGridFiller.cpp index 6e2f8cfd114..88e6dd608b0 100644 --- a/Core/src/Detector/detail/IndexedGridFiller.cpp +++ b/Core/src/Detector/detail/IndexedGridFiller.cpp @@ -47,14 +47,10 @@ std::vector Acts::Experimental::detail::binSequence( rBins.reserve(bmax - bmin + 1u + 2 * expand); // handle bmin:/max expand it down (for bound, don't fill underflow) if (type == Acts::AxisBoundaryType::Bound) { - bmin = (static_cast(bmin) - static_cast(expand) > 0) - ? bmin - expand - : 1u; + bmin = bmin > expand ? bmin - expand : 1u; bmax = (bmax + expand <= nBins) ? bmax + expand : nBins; } else if (type == Acts::AxisBoundaryType::Open) { - bmin = (static_cast(bmin) - static_cast(expand) >= 0) - ? bmin - expand - : 0u; + bmin = bmin >= expand ? bmin - expand : 0u; bmax = (bmax + expand <= nBins + 1u) ? bmax + expand : nBins + 1u; } fill_linear(bmin, bmax); @@ -62,14 +58,11 @@ std::vector Acts::Experimental::detail::binSequence( // Close case std::size_t span = bmax - bmin + 1u + 2 * expand; // Safe with respect to the closure point, treat as bound - if (2 * span < nBins && (bmax + expand <= nBins) && - (static_cast(bmin) - static_cast(expand) > 0)) { + if (2 * span < nBins && (bmax + expand <= nBins) && (bmin > expand)) { return binSequence({bmin, bmax}, expand, nBins, Acts::AxisBoundaryType::Bound); } else if (2 * span < nBins) { - bmin = static_cast(bmin) - static_cast(expand) > 0 - ? bmin - expand - : 1u; + bmin = bmin > expand ? bmin - expand : 1u; bmax = bmax + expand <= nBins ? bmax + expand : nBins; fill_linear(bmin, bmax); // deal with expansions over the phi boundary @@ -77,9 +70,8 @@ std::vector Acts::Experimental::detail::binSequence( std::size_t overstep = (bmax + expand - nBins); fill_linear(1u, overstep); } - if (static_cast(bmin) - static_cast(expand) < 1) { - std::size_t understep = - abs(static_cast(bmin) - static_cast(expand)); + if (bmin <= expand) { + std::size_t understep = expand - bmin; fill_linear(nBins - understep, nBins); } std::ranges::sort(rBins); diff --git a/Core/src/MagneticField/CMakeLists.txt b/Core/src/MagneticField/CMakeLists.txt index 4b9ebb0be0f..3b164d39505 100644 --- a/Core/src/MagneticField/CMakeLists.txt +++ b/Core/src/MagneticField/CMakeLists.txt @@ -1,4 +1,8 @@ target_sources( ActsCore - PRIVATE BFieldMapUtils.cpp SolenoidBField.cpp MagneticFieldError.cpp + PRIVATE + BFieldMapUtils.cpp + SolenoidBField.cpp + MagneticFieldError.cpp + MultiRangeBField.cpp ) diff --git a/Core/src/MagneticField/MultiRangeBField.cpp b/Core/src/MagneticField/MultiRangeBField.cpp new file mode 100644 index 00000000000..8899b50d802 --- /dev/null +++ b/Core/src/MagneticField/MultiRangeBField.cpp @@ -0,0 +1,78 @@ +// This file is part of the ACTS project. +// +// Copyright (C) 2016 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 https://mozilla.org/MPL/2.0/. + +#include "Acts/MagneticField/MultiRangeBField.hpp" + +namespace Acts { + +MultiRangeBField::Cache::Cache(const MagneticFieldContext& /*unused*/) {} + +MultiRangeBField::MultiRangeBField(const std::vector& ranges) + : fieldRanges(ranges) {} + +MultiRangeBField::MultiRangeBField( + std::vector&& ranges) + : fieldRanges(std::move(ranges)) {} + +MagneticFieldProvider::Cache MultiRangeBField::makeCache( + const MagneticFieldContext& mctx) const { + return MagneticFieldProvider::Cache(std::in_place_type, mctx); +} + +Result MultiRangeBField::getField( + const Vector3& position, MagneticFieldProvider::Cache& cache) const { + // Because we assume that the regions are in increasing order of + // precedence, we can iterate over the array, remembering the _last_ + // region that contained the given point. At the end of the loop, this + // region will then be the one we query for its field vector. + std::optional foundRange = {}; + + // The cache object for this type contains an optional integer; if the + // integer is set, it indicates the index of the region in which the last + // access succeeded. This acts as a cache because if the current access is + // still within that region, it is guaranteed that none of the regions + // with lower priority -- i.e. that are earlier in the region vector -- + // can be relevant to the current access. Thus, we request the cache index + // if it exists and perform a membership check on it; if that succeeds, we + // remember the corresponding region as a candidate. + if (Cache& lCache = cache.as(); + lCache.index.has_value() && + std::get<0>(fieldRanges[*lCache.index]) + .contains({position[0], position[1], position[2]})) { + foundRange = *lCache.index; + } + + // Now, we iterate over the ranges. If we already have a range candidate, + // we start iteration just after the existing candidate; otherwise we + // start from the beginning. + for (std::size_t i = (foundRange.has_value() ? (*foundRange) + 1 : 0); + i < fieldRanges.size(); ++i) { + if (std::get<0>(fieldRanges[i]) + .contains({position[0], position[1], position[2]})) { + foundRange = i; + } + } + + // Update the cache using the result of this access. + cache.as().index = foundRange; + + // If we found a valid range, return the corresponding vector; otherwise + // return an out-of-bounds error. + if (foundRange.has_value()) { + return Result::success(std::get<1>(fieldRanges[*foundRange])); + } else { + return Result::failure(MagneticFieldError::OutOfBounds); + } +} + +Result MultiRangeBField::getFieldGradient( + const Vector3& position, ActsMatrix<3, 3>& /*unused*/, + MagneticFieldProvider::Cache& cache) const { + return getField(position, cache); +} +} // namespace Acts diff --git a/Core/src/Propagator/detail/CovarianceEngine.cpp b/Core/src/Propagator/detail/CovarianceEngine.cpp index a89016bce3c..ca2f57b88cf 100644 --- a/Core/src/Propagator/detail/CovarianceEngine.cpp +++ b/Core/src/Propagator/detail/CovarianceEngine.cpp @@ -97,8 +97,7 @@ CurvilinearState detail::curvilinearState( pos4, direction, freeParameters[eFreeQOverP], std::move(cov), particleHypothesis); // Create the curvilinear state - return std::make_tuple(std::move(curvilinearParams), fullTransportJacobian, - accumulatedPath); + return {std::move(curvilinearParams), fullTransportJacobian, accumulatedPath}; } void detail::transportCovarianceToBound( diff --git a/Core/src/Propagator/detail/SympyCovarianceEngine.cpp b/Core/src/Propagator/detail/SympyCovarianceEngine.cpp index f10171cc5ec..488a2d69ec9 100644 --- a/Core/src/Propagator/detail/SympyCovarianceEngine.cpp +++ b/Core/src/Propagator/detail/SympyCovarianceEngine.cpp @@ -86,8 +86,7 @@ CurvilinearState sympy::curvilinearState( pos4, direction, freeParameters[eFreeQOverP], std::move(cov), particleHypothesis); // Create the curvilinear state - return std::make_tuple(std::move(curvilinearParams), fullTransportJacobian, - accumulatedPath); + return {std::move(curvilinearParams), fullTransportJacobian, accumulatedPath}; } void sympy::transportCovarianceToBound( diff --git a/Core/src/Surfaces/detail/AlignmentHelper.cpp b/Core/src/Surfaces/detail/AlignmentHelper.cpp index 92c1fb7d2a1..b5308d60c84 100644 --- a/Core/src/Surfaces/detail/AlignmentHelper.cpp +++ b/Core/src/Surfaces/detail/AlignmentHelper.cpp @@ -87,6 +87,6 @@ Acts::detail::RotationToAxes Acts::detail::rotationToLocalAxesDerivative( rotToCompositeLocalYAxis * relRotation(1, 2) + rotToCompositeLocalZAxis * relRotation(2, 2); - return std::make_tuple(std::move(rotToLocalXAxis), std::move(rotToLocalYAxis), - std::move(rotToLocalZAxis)); + return {std::move(rotToLocalXAxis), std::move(rotToLocalYAxis), + std::move(rotToLocalZAxis)}; } diff --git a/Core/src/Surfaces/detail/AnnulusBoundsHelper.cpp b/Core/src/Surfaces/detail/AnnulusBoundsHelper.cpp index b46bfaf5b8a..4ceab6d594b 100644 --- a/Core/src/Surfaces/detail/AnnulusBoundsHelper.cpp +++ b/Core/src/Surfaces/detail/AnnulusBoundsHelper.cpp @@ -63,5 +63,5 @@ Acts::detail::AnnulusBoundsHelper::create(const Transform3& transform, auto annulusBounds = std::make_shared( rMin, rMax, phiMin, phiMax, originShift, phiShift); - return std::make_tuple(annulusBounds, boundsTransform); + return {annulusBounds, boundsTransform}; } diff --git a/Core/src/TrackFinding/AmbiguityTrackClustering.cpp b/Core/src/TrackFinding/AmbiguityTrackClustering.cpp index ff0e95cbb3d..8586f571aec 100644 --- a/Core/src/TrackFinding/AmbiguityTrackClustering.cpp +++ b/Core/src/TrackFinding/AmbiguityTrackClustering.cpp @@ -21,8 +21,9 @@ Acts::detail::clusterDuplicateTracks( // different clusters. std::unordered_map hitToTrack; - // Loop over all the tracks - for (const auto& [_, trackValue] : trackMap) { + // Loop backward over all the tracks + for (auto track = trackMap.rbegin(); track != trackMap.rend(); ++track) { + const auto& trackValue = track->second; std::vector hits = trackValue.second; auto matchedTrack = hitToTrack.end(); // Loop over all the hits in the track diff --git a/Core/src/TrackFitting/GlobalChiSquareFitter.cpp b/Core/src/TrackFitting/GlobalChiSquareFitter.cpp index 50b8b4cab8c..4062124496e 100644 --- a/Core/src/TrackFitting/GlobalChiSquareFitter.cpp +++ b/Core/src/TrackFitting/GlobalChiSquareFitter.cpp @@ -50,3 +50,98 @@ void Acts::Experimental::updateGx2fCovarianceParams( return; } + +void Acts::Experimental::addMeasurementToGx2fSumsBackend( + Gx2fSystem& extendedSystem, + const std::vector& jacobianFromStart, + const Eigen::MatrixXd& covarianceMeasurement, const BoundVector& predicted, + const Eigen::VectorXd& measurement, const Eigen::MatrixXd& projector, + const Logger& logger) { + // First, w try to invert the covariance matrix. If the inversion fails, we + // can already abort. + const auto safeInvCovMeasurement = safeInverse(covarianceMeasurement); + if (!safeInvCovMeasurement) { + ACTS_WARNING("addMeasurementToGx2fSums: safeInvCovMeasurement failed."); + ACTS_VERBOSE(" covarianceMeasurement:\n" << covarianceMeasurement); + return; + } + + // Create an extended Jacobian. This one contains only eBoundSize rows, + // because the rest is irrelevant. We fill it in the next steps. + // TODO make dimsExtendedParams template with unrolling + Eigen::MatrixXd extendedJacobian = + Eigen::MatrixXd::Zero(eBoundSize, extendedSystem.nDims()); + + // This part of the Jacobian comes from the material-less propagation + extendedJacobian.topLeftCorner() = + jacobianFromStart[0]; + + // If we have material, loop here over all Jacobians. We add extra columns for + // their phi-theta projections. These parts account for the propagation of the + // scattering angles. + for (std::size_t matSurface = 1; matSurface < jacobianFromStart.size(); + matSurface++) { + const BoundMatrix jac = jacobianFromStart[matSurface]; + + const ActsMatrix jacPhiTheta = + jac * Gx2fConstants::phiThetaProjector; + + // The position, where we need to insert the values in the extended Jacobian + const std::size_t deltaPosition = eBoundSize + 2 * (matSurface - 1); + + extendedJacobian.template block(0, deltaPosition) = + jacPhiTheta; + } + + const Eigen::MatrixXd projJacobian = projector * extendedJacobian; + + const Eigen::VectorXd projPredicted = projector * predicted; + + const Eigen::VectorXd residual = measurement - projPredicted; + + // Finally contribute to chi2sum, aMatrix, and bVector + extendedSystem.chi2() += + (residual.transpose() * (*safeInvCovMeasurement) * residual)(0, 0); + + extendedSystem.aMatrix() += + (projJacobian.transpose() * (*safeInvCovMeasurement) * projJacobian) + .eval(); + + extendedSystem.bVector() += + (residual.transpose() * (*safeInvCovMeasurement) * projJacobian) + .eval() + .transpose(); + + ACTS_VERBOSE( + "Contributions in addMeasurementToGx2fSums:\n" + << " predicted: " << predicted.transpose() << "\n" + << " measurement: " << measurement.transpose() << "\n" + << " covarianceMeasurement:\n" + << covarianceMeasurement << "\n" + << " projector:\n" + << projector.eval() << "\n" + << " projJacobian:\n" + << projJacobian.eval() << "\n" + << " projPredicted: " << (projPredicted.transpose()).eval() << "\n" + << " residual: " << (residual.transpose()).eval() << "\n" + << " extendedJacobian:\n" + << extendedJacobian << "\n" + << " aMatrix contribution:\n" + << (projJacobian.transpose() * (*safeInvCovMeasurement) * projJacobian) + .eval() + << "\n" + << " bVector contribution: " + << (residual.transpose() * (*safeInvCovMeasurement) * projJacobian).eval() + << "\n" + << " chi2sum contribution: " + << (residual.transpose() * (*safeInvCovMeasurement) * residual)(0, 0) + << "\n" + << " safeInvCovMeasurement:\n" + << (*safeInvCovMeasurement)); +} + +Eigen::VectorXd Acts::Experimental::computeGx2fDeltaParams( + const Acts::Experimental::Gx2fSystem& extendedSystem) { + return extendedSystem.aMatrix().colPivHouseholderQr().solve( + extendedSystem.bVector()); +} diff --git a/Core/src/Utilities/SpacePointUtility.cpp b/Core/src/Utilities/SpacePointUtility.cpp index 32470e1ebb9..6f662fe20bd 100644 --- a/Core/src/Utilities/SpacePointUtility.cpp +++ b/Core/src/Utilities/SpacePointUtility.cpp @@ -99,7 +99,7 @@ SpacePointUtility::globalCoords( tcov = std::nullopt; } - return std::make_tuple(globalPos, globalTime, gcov, tcov); + return {globalPos, globalTime, gcov, tcov}; } Vector2 SpacePointUtility::calcRhoZVars( diff --git a/Core/src/Vertexing/AdaptiveGridTrackDensity.cpp b/Core/src/Vertexing/AdaptiveGridTrackDensity.cpp index d08b9db4913..31daae4b039 100644 --- a/Core/src/Vertexing/AdaptiveGridTrackDensity.cpp +++ b/Core/src/Vertexing/AdaptiveGridTrackDensity.cpp @@ -156,7 +156,7 @@ AdaptiveGridTrackDensity::getMaxZTPosition(DensityMap& densityMap) const { double maxZ = getSpatialBinCenter(bin.first); double maxT = getTemporalBinCenter(bin.second); - return std::make_pair(maxZ, maxT); + return std::pair(maxZ, maxT); } Result diff --git a/Core/src/Vertexing/ImpactPointEstimator.cpp b/Core/src/Vertexing/ImpactPointEstimator.cpp index df9f75fd9fa..51f0ebe24db 100644 --- a/Core/src/Vertexing/ImpactPointEstimator.cpp +++ b/Core/src/Vertexing/ImpactPointEstimator.cpp @@ -217,7 +217,7 @@ Result> getDistanceAndMomentumImpl( Vector4 deltaRStraightTrack{Vector4::Zero()}; deltaRStraightTrack.head() = pcaStraightTrack - vtxPos; - return std::make_pair(deltaRStraightTrack, momDirStraightTrack); + return std::pair(deltaRStraightTrack, momDirStraightTrack); } // Charged particles in a constant B field follow a helical trajectory. In @@ -291,7 +291,7 @@ Result> getDistanceAndMomentumImpl( Vector4 deltaR{Vector4::Zero()}; deltaR.head() = pca - vtxPos; - return std::make_pair(deltaR, momDir); + return std::pair(deltaR, momDir); } } // namespace diff --git a/Core/src/Vertexing/Vertex.cpp b/Core/src/Vertexing/Vertex.cpp index 7c2c54cc7ac..4368f3b2faa 100644 --- a/Core/src/Vertexing/Vertex.cpp +++ b/Core/src/Vertexing/Vertex.cpp @@ -74,7 +74,7 @@ const std::vector& Vertex::tracks() const { } std::pair Vertex::fitQuality() const { - return std::pair(m_chiSquared, m_numberDoF); + return {m_chiSquared, m_numberDoF}; } void Vertex::setPosition(const Vector3& position) { diff --git a/Examples/Algorithms/Digitization/include/ActsExamples/Digitization/DigitizationAlgorithm.hpp b/Examples/Algorithms/Digitization/include/ActsExamples/Digitization/DigitizationAlgorithm.hpp index 331408faa9a..a7e9b31bb7d 100644 --- a/Examples/Algorithms/Digitization/include/ActsExamples/Digitization/DigitizationAlgorithm.hpp +++ b/Examples/Algorithms/Digitization/include/ActsExamples/Digitization/DigitizationAlgorithm.hpp @@ -48,6 +48,10 @@ class DigitizationAlgorithm final : public IAlgorithm { std::string outputMeasurementParticlesMap = "measurement_particles_map"; /// Output collection to map measured hits to simulated hits. std::string outputMeasurementSimHitsMap = "measurement_simhits_map"; + /// Output collection to map particles to measurements. + std::string outputParticleMeasurementsMap = "particle_measurements_map"; + /// Output collection to map particles to simulated hits. + std::string outputSimHitMeasurementsMap = "simhit_measurements_map"; /// Map of surface by identifier to allow local - to global std::unordered_map @@ -122,7 +126,7 @@ class DigitizationAlgorithm final : public IAlgorithm { Config m_cfg; /// Digitizers within geometry hierarchy Acts::GeometryHierarchyMap m_digitizers; - /// Geometric digtizer + /// Geometric digitizer ActsFatras::Channelizer m_channelizer; using CellsMap = @@ -140,6 +144,11 @@ class DigitizationAlgorithm final : public IAlgorithm { WriteDataHandle> m_outputMeasurementSimHitsMap{ this, "OutputMeasurementSimHitsMap"}; + WriteDataHandle> m_outputParticleMeasurementsMap{ + this, "OutputParticleMeasurementsMap"}; + WriteDataHandle> m_outputSimHitMeasurementsMap{ + this, "OutputSimHitMeasurementsMap"}; + /// Construct a fixed-size smearer from a configuration. /// /// It's templated on the smearing dimension given by @tparam kSmearDIM @@ -153,7 +162,7 @@ class DigitizationAlgorithm final : public IAlgorithm { // Copy the geometric configuration impl.geometric = cfg.geometricDigiConfig; // Prepare the smearing configuration - for (int i = 0; i < static_cast(kSmearDIM); ++i) { + for (std::size_t i = 0; i < kSmearDIM; ++i) { impl.smearing.indices[i] = cfg.smearingDigiConfig.at(i).index; impl.smearing.smearFunctions[i] = cfg.smearingDigiConfig.at(i).smearFunction; diff --git a/Examples/Algorithms/Digitization/src/DigitizationAlgorithm.cpp b/Examples/Algorithms/Digitization/src/DigitizationAlgorithm.cpp index d7b4754cb3a..e55ac39d09b 100644 --- a/Examples/Algorithms/Digitization/src/DigitizationAlgorithm.cpp +++ b/Examples/Algorithms/Digitization/src/DigitizationAlgorithm.cpp @@ -16,12 +16,9 @@ #include "ActsExamples/EventData/GeometryContainers.hpp" #include "ActsExamples/EventData/Index.hpp" #include "ActsExamples/Framework/AlgorithmContext.hpp" -#include "ActsExamples/Utilities/Range.hpp" -#include "ActsFatras/EventData/Barcode.hpp" #include #include -#include #include #include #include @@ -61,12 +58,23 @@ DigitizationAlgorithm::DigitizationAlgorithm(Config config, throw std::invalid_argument( "Missing hit-to-simulated-hits map output collection"); } + if (m_cfg.outputParticleMeasurementsMap.empty()) { + throw std::invalid_argument( + "Missing particle-to-measurements map output collection"); + } + if (m_cfg.outputSimHitMeasurementsMap.empty()) { + throw std::invalid_argument( + "Missing particle-to-simulated-hits map output collection"); + } m_outputMeasurements.initialize(m_cfg.outputMeasurements); m_outputClusters.initialize(m_cfg.outputClusters); m_outputMeasurementParticlesMap.initialize( m_cfg.outputMeasurementParticlesMap); m_outputMeasurementSimHitsMap.initialize(m_cfg.outputMeasurementSimHitsMap); + m_outputParticleMeasurementsMap.initialize( + m_cfg.outputParticleMeasurementsMap); + m_outputSimHitMeasurementsMap.initialize(m_cfg.outputSimHitMeasurementsMap); } if (m_cfg.doOutputCells) { @@ -141,7 +149,7 @@ ProcessCode DigitizationAlgorithm::execute(const AlgorithmContext& ctx) const { MeasurementContainer measurements; ClusterContainer clusters; - IndexMultimap measurementParticlesMap; + IndexMultimap measurementParticlesMap; IndexMultimap measurementSimHitsMap; measurements.reserve(simHits.size()); measurementParticlesMap.reserve(simHits.size()); @@ -302,6 +310,12 @@ ProcessCode DigitizationAlgorithm::execute(const AlgorithmContext& ctx) const { m_outputMeasurements(ctx, std::move(measurements)); m_outputClusters(ctx, std::move(clusters)); + // invert them before they are moved + m_outputParticleMeasurementsMap( + ctx, invertIndexMultimap(measurementParticlesMap)); + m_outputSimHitMeasurementsMap(ctx, + invertIndexMultimap(measurementSimHitsMap)); + m_outputMeasurementParticlesMap(ctx, std::move(measurementParticlesMap)); m_outputMeasurementSimHitsMap(ctx, std::move(measurementSimHitsMap)); } diff --git a/Examples/Algorithms/TrackFindingML/CMakeLists.txt b/Examples/Algorithms/TrackFindingML/CMakeLists.txt index d826b224f2c..80f55f60579 100644 --- a/Examples/Algorithms/TrackFindingML/CMakeLists.txt +++ b/Examples/Algorithms/TrackFindingML/CMakeLists.txt @@ -1,7 +1,5 @@ set(SOURCES - src/AmbiguityResolutionML.cpp src/AmbiguityResolutionMLAlgorithm.cpp - src/AmbiguityResolutionMLDBScanAlgorithm.cpp src/SeedFilterMLAlgorithm.cpp ) diff --git a/Examples/Algorithms/TrackFindingML/include/ActsExamples/TrackFindingML/AmbiguityDBScanClustering.hpp b/Examples/Algorithms/TrackFindingML/include/ActsExamples/TrackFindingML/AmbiguityDBScanClustering.hpp deleted file mode 100644 index 45fe4cb5480..00000000000 --- a/Examples/Algorithms/TrackFindingML/include/ActsExamples/TrackFindingML/AmbiguityDBScanClustering.hpp +++ /dev/null @@ -1,79 +0,0 @@ -// This file is part of the ACTS project. -// -// Copyright (C) 2016 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 https://mozilla.org/MPL/2.0/. - -#pragma once - -#include "Acts/EventData/TrackContainer.hpp" -#include "Acts/EventData/TrackContainerFrontendConcept.hpp" -#include "Acts/TrackFinding/detail/AmbiguityTrackClustering.hpp" -#include "Acts/Utilities/DBScan.hpp" - -#include -#include -#include - -namespace Acts { - -/// Clusterise tracks based on shared hits -/// -/// @param trackMap Multimap storing pair of track ID and vector of measurement ID. The keys are the number of measurement and are just there to facilitate the ordering. -/// @param tracks Track container with all the track to be clustered -/// @param epsilon Maximum distance between 2 tracks to be clustered -/// @param minPoints Minimum number of tracks to create a cluster -/// @return an unordered map representing the clusters, the keys the ID of the primary track of each cluster and the store a vector of track IDs. -template -std::unordered_map> dbscanTrackClustering( - std::multimap>>& - trackMap, - const track_container_t& tracks, float epsilon = 0.07, int minPoints = 2) { - // Unordered map associating a vector with all the track ID of a cluster to - // the ID of the first track of the cluster - std::unordered_map> cluster; - // Unordered map associating hits to the ID of the first track of the - // different clusters. - std::unordered_map hitToTrack; - - // Initialize a DBScan of dimension 4 (phi, eta, z, Pt) - using DBSCAN = Acts::DBScan<4, double, 4>; - DBSCAN dbscan(epsilon, minPoints, true); - - std::vector> data; - std::size_t trackID = 0; - std::vector clusterAssignments; - - // Get the input feature of the network for all the tracks - for (const auto& [key, val] : trackMap) { - auto traj = tracks.getTrack(val.first); - data.push_back({Acts::VectorHelpers::eta(traj.momentum()), - Acts::VectorHelpers::phi(traj.momentum())}); - } - std::size_t clusterNb = dbscan.cluster(data, clusterAssignments); - - // Cluster track with DBScan - std::vector< - std::multimap>>> - dbscanClusters(clusterNb); - for (const auto& [key, val] : trackMap) { - std::size_t clusterID = clusterAssignments[trackID]; - dbscanClusters[clusterID].emplace(key, val); - trackID++; - } - - // Perform a subClustering of the DBScan cluster using the measurement ID - // clustering - for (const auto& dbscanCluster : dbscanClusters) { - auto subCluster = Acts::detail::clusterDuplicateTracks(dbscanCluster); - cluster.merge(subCluster); - if (!subCluster.empty()) { - std::cout << "Overlapping track ID, there must be an error" << std::endl; - } - } - return cluster; -} - -} // namespace Acts diff --git a/Examples/Algorithms/TrackFindingML/include/ActsExamples/TrackFindingML/AmbiguityResolutionML.hpp b/Examples/Algorithms/TrackFindingML/include/ActsExamples/TrackFindingML/AmbiguityResolutionML.hpp deleted file mode 100644 index ea967a23c22..00000000000 --- a/Examples/Algorithms/TrackFindingML/include/ActsExamples/TrackFindingML/AmbiguityResolutionML.hpp +++ /dev/null @@ -1,48 +0,0 @@ -// This file is part of the ACTS project. -// -// Copyright (C) 2016 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 https://mozilla.org/MPL/2.0/. - -#pragma once - -#include "ActsExamples/EventData/Track.hpp" -#include "ActsExamples/Framework/IAlgorithm.hpp" - -#include -#include -#include - -namespace ActsExamples { - -/// Generic implementation of the machine learning ambiguity resolution -/// Contains method for data preparations -class AmbiguityResolutionML : public IAlgorithm { - public: - /// Construct the ambiguity resolution algorithm. - /// - /// @param name name of the algorithm - /// @param lvl is the logging level - AmbiguityResolutionML(std::string name, Acts::Logging::Level lvl); - - protected: - /// Associated measurements ID to Tracks ID - /// - /// @param tracks is the input track container - /// @param nMeasurementsMin minimum number of measurement per track - /// @return an ordered list containing pairs of track ID and associated measurement ID - std::multimap>> - mapTrackHits(const ConstTrackContainer& tracks, int nMeasurementsMin) const; - - /// Prepare the output track container to be written - /// - /// @param tracks is the input track container - /// @param goodTracks is list of the IDs of all the tracks we want to keep - ConstTrackContainer prepareOutputTrack( - const ConstTrackContainer& tracks, - std::vector& goodTracks) const; -}; - -} // namespace ActsExamples diff --git a/Examples/Algorithms/TrackFindingML/include/ActsExamples/TrackFindingML/AmbiguityResolutionMLAlgorithm.hpp b/Examples/Algorithms/TrackFindingML/include/ActsExamples/TrackFindingML/AmbiguityResolutionMLAlgorithm.hpp index 2e5d22ac01d..e591b2e201a 100644 --- a/Examples/Algorithms/TrackFindingML/include/ActsExamples/TrackFindingML/AmbiguityResolutionMLAlgorithm.hpp +++ b/Examples/Algorithms/TrackFindingML/include/ActsExamples/TrackFindingML/AmbiguityResolutionMLAlgorithm.hpp @@ -8,10 +8,11 @@ #pragma once +#include "Acts/AmbiguityResolution/AmbiguityResolutionML.hpp" #include "Acts/Plugins/Onnx/AmbiguityTrackClassifier.hpp" #include "ActsExamples/EventData/Track.hpp" #include "ActsExamples/Framework/DataHandle.hpp" -#include "ActsExamples/TrackFindingML/AmbiguityResolutionML.hpp" +#include "ActsExamples/Framework/IAlgorithm.hpp" #include @@ -23,7 +24,10 @@ namespace ActsExamples { /// 1) Cluster together nearby tracks using shared hits /// 2) For each track use a neural network to compute a score /// 3) In each cluster keep the track with the highest score -class AmbiguityResolutionMLAlgorithm final : public AmbiguityResolutionML { +class AmbiguityResolutionMLAlgorithm final : public IAlgorithm { + using AmbiguityResolution = + Acts::AmbiguityResolutionML; + public: struct Config { /// Input track collection. @@ -33,7 +37,11 @@ class AmbiguityResolutionMLAlgorithm final : public AmbiguityResolutionML { /// Output track collection. std::string outputTracks; /// Minimum number of measurement to form a track. - int nMeasurementsMin = 7; + std::size_t nMeasurementsMin = 7; + /// Construct the ML ambiguity resolution configuration. + AmbiguityResolution::Config toAmbiguityResolutionMLConfig() const { + return {inputDuplicateNN, nMeasurementsMin}; + } }; /// Construct the ambiguity resolution algorithm. @@ -53,8 +61,7 @@ class AmbiguityResolutionMLAlgorithm final : public AmbiguityResolutionML { private: Config m_cfg; - // ONNX model for track selection - Acts::AmbiguityTrackClassifier m_duplicateClassifier; + AmbiguityResolution m_ambiML; ReadDataHandle m_inputTracks{this, "InputTracks"}; WriteDataHandle m_outputTracks{this, "OutputTracks"}; }; diff --git a/Examples/Algorithms/TrackFindingML/include/ActsExamples/TrackFindingML/AmbiguityResolutionMLDBScanAlgorithm.hpp b/Examples/Algorithms/TrackFindingML/include/ActsExamples/TrackFindingML/AmbiguityResolutionMLDBScanAlgorithm.hpp deleted file mode 100644 index 6ad387ffdc4..00000000000 --- a/Examples/Algorithms/TrackFindingML/include/ActsExamples/TrackFindingML/AmbiguityResolutionMLDBScanAlgorithm.hpp +++ /dev/null @@ -1,68 +0,0 @@ -// This file is part of the ACTS project. -// -// Copyright (C) 2016 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 https://mozilla.org/MPL/2.0/. - -#pragma once - -#include "Acts/Plugins/Onnx/AmbiguityTrackClassifier.hpp" -#include "ActsExamples/EventData/Track.hpp" -#include "ActsExamples/Framework/DataHandle.hpp" -#include "ActsExamples/TrackFindingML/AmbiguityResolutionML.hpp" - -#include - -namespace ActsExamples { - -/// Evicts tracks that seem to be duplicated and fake. -/// -/// The implementation works as follows: -/// 1) Cluster together nearby tracks using a DBScan -/// 2) Create subcluster based on tracks with shared hits -/// 3) For each track use a neural network to compute a score -/// 4) In each cluster keep the track with the highest score -class AmbiguityResolutionMLDBScanAlgorithm final - : public AmbiguityResolutionML { - public: - struct Config { - /// Input trajectories collection. - std::string inputTracks; - /// Path to the ONNX model for the duplicate neural network - std::string inputDuplicateNN; - /// Output trajectories collection. - std::string outputTracks; - /// Minimum number of measurement to form a track. - int nMeasurementsMin = 7; - /// Maximum distance between 2 tracks to be clustered in the DBScan - float epsilonDBScan = 0.07; - /// Minimum number of tracks to create a cluster in the DBScan - int minPointsDBScan = 2; - }; - - /// Construct the ambiguity resolution algorithm. - /// - /// @param cfg is the algorithm configuration - /// @param lvl is the logging level - AmbiguityResolutionMLDBScanAlgorithm(Config cfg, Acts::Logging::Level lvl); - - /// Run the ambiguity resolution algorithm. - /// - /// @param cxt is the algorithm context with event information - /// @return a process code indication success or failure - ProcessCode execute(const AlgorithmContext& ctx) const final; - - /// Const access to the config - const Config& config() const { return m_cfg; } - - private: - Config m_cfg; - // ONNX model for track selection - Acts::AmbiguityTrackClassifier m_duplicateClassifier; - ReadDataHandle m_inputTracks{this, "InputTracks"}; - WriteDataHandle m_outputTracks{this, "OutputTracks"}; -}; - -} // namespace ActsExamples diff --git a/Examples/Algorithms/TrackFindingML/include/ActsExamples/TrackFindingML/SeedFilterMLAlgorithm.hpp b/Examples/Algorithms/TrackFindingML/include/ActsExamples/TrackFindingML/SeedFilterMLAlgorithm.hpp index 6472a6e8d2b..53659a33e4f 100644 --- a/Examples/Algorithms/TrackFindingML/include/ActsExamples/TrackFindingML/SeedFilterMLAlgorithm.hpp +++ b/Examples/Algorithms/TrackFindingML/include/ActsExamples/TrackFindingML/SeedFilterMLAlgorithm.hpp @@ -12,7 +12,7 @@ #include "ActsExamples/EventData/SimSeed.hpp" #include "ActsExamples/EventData/Track.hpp" #include "ActsExamples/Framework/DataHandle.hpp" -#include "ActsExamples/TrackFindingML/AmbiguityResolutionML.hpp" +#include "ActsExamples/Framework/IAlgorithm.hpp" #include diff --git a/Examples/Algorithms/TrackFindingML/src/AmbiguityResolutionML.cpp b/Examples/Algorithms/TrackFindingML/src/AmbiguityResolutionML.cpp deleted file mode 100644 index f3282a15d94..00000000000 --- a/Examples/Algorithms/TrackFindingML/src/AmbiguityResolutionML.cpp +++ /dev/null @@ -1,76 +0,0 @@ -// This file is part of the ACTS project. -// -// Copyright (C) 2016 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 https://mozilla.org/MPL/2.0/. - -#include "ActsExamples/TrackFindingML/AmbiguityResolutionML.hpp" - -#include "ActsExamples/EventData/IndexSourceLink.hpp" -#include "ActsExamples/EventData/Measurement.hpp" - -ActsExamples::AmbiguityResolutionML::AmbiguityResolutionML( - std::string name, Acts::Logging::Level lvl) - : ActsExamples::IAlgorithm(name, lvl) {} - -std::multimap>> -ActsExamples::AmbiguityResolutionML::mapTrackHits( - const ActsExamples::ConstTrackContainer& tracks, - int nMeasurementsMin) const { - std::multimap>> trackMap; - // Loop over all the trajectories in the events - for (const auto& track : tracks) { - std::vector hits; - int nbMeasurements = 0; - // Store the hits id for the trajectory and compute the number of - // measurement - tracks.trackStateContainer().visitBackwards( - track.tipIndex(), [&](const auto& state) { - if (state.typeFlags().test(Acts::TrackStateFlag::MeasurementFlag)) { - std::size_t indexHit = - state.getUncalibratedSourceLink() - .template get() - .index(); - hits.emplace_back(indexHit); - ++nbMeasurements; - } - }); - if (nbMeasurements < nMeasurementsMin) { - continue; - } - trackMap.emplace(nbMeasurements, std::make_pair(track.index(), hits)); - } - return trackMap; -} - -ActsExamples::ConstTrackContainer -ActsExamples::AmbiguityResolutionML::prepareOutputTrack( - const ActsExamples::ConstTrackContainer& tracks, - std::vector& goodTracks) const { - std::shared_ptr trackStateContainer = - tracks.trackStateContainerHolder(); - auto trackContainer = std::make_shared(); - trackContainer->reserve(goodTracks.size()); - // Temporary empty track state container: we don't change the original one, - // but we need one for filtering - auto tempTrackStateContainer = - std::make_shared(); - - TrackContainer solvedTracks{trackContainer, tempTrackStateContainer}; - solvedTracks.ensureDynamicColumns(tracks); - - for (auto&& iTrack : goodTracks) { - auto destProxy = solvedTracks.makeTrack(); - auto srcProxy = tracks.getTrack(iTrack); - destProxy.copyFrom(srcProxy, false); - destProxy.tipIndex() = srcProxy.tipIndex(); - } - - ConstTrackContainer outputTracks{ - std::make_shared( - std::move(*trackContainer)), - trackStateContainer}; - return outputTracks; -} diff --git a/Examples/Algorithms/TrackFindingML/src/AmbiguityResolutionMLAlgorithm.cpp b/Examples/Algorithms/TrackFindingML/src/AmbiguityResolutionMLAlgorithm.cpp index ba2e9abeb98..805913e60ae 100644 --- a/Examples/Algorithms/TrackFindingML/src/AmbiguityResolutionMLAlgorithm.cpp +++ b/Examples/Algorithms/TrackFindingML/src/AmbiguityResolutionMLAlgorithm.cpp @@ -8,18 +8,30 @@ #include "ActsExamples/TrackFindingML/AmbiguityResolutionMLAlgorithm.hpp" +#include "ActsExamples/EventData/IndexSourceLink.hpp" +#include "ActsExamples/EventData/Measurement.hpp" #include "ActsExamples/Framework/ProcessCode.hpp" #include #include +static std::size_t sourceLinkHash(const Acts::SourceLink& a) { + return static_cast( + a.get().index()); +} + +static bool sourceLinkEquality(const Acts::SourceLink& a, + const Acts::SourceLink& b) { + return a.get().index() == + b.get().index(); +} + ActsExamples::AmbiguityResolutionMLAlgorithm::AmbiguityResolutionMLAlgorithm( ActsExamples::AmbiguityResolutionMLAlgorithm::Config cfg, Acts::Logging::Level lvl) - : ActsExamples::AmbiguityResolutionML("AmbiguityResolutionMLAlgorithm", - lvl), + : ActsExamples::IAlgorithm("AmbiguityResolutionMLAlgorithm", lvl), m_cfg(std::move(cfg)), - m_duplicateClassifier(m_cfg.inputDuplicateNN.c_str()) { + m_ambiML(m_cfg.toAmbiguityResolutionMLConfig(), logger().clone()) { if (m_cfg.inputTracks.empty()) { throw std::invalid_argument("Missing trajectories input collection"); } @@ -34,15 +46,32 @@ ActsExamples::ProcessCode ActsExamples::AmbiguityResolutionMLAlgorithm::execute( const AlgorithmContext& ctx) const { // Read input data const auto& tracks = m_inputTracks(ctx); - // Associate measurement to their respective tracks + // Associate measurement to their respective tracks to prepare the track + // shared hits based clustering std::multimap>> - trackMap = mapTrackHits(tracks, m_cfg.nMeasurementsMin); + trackMap = + m_ambiML.mapTrackHits(tracks, &sourceLinkHash, &sourceLinkEquality); + // Cluster the tracks based on the shared hits auto cluster = Acts::detail::clusterDuplicateTracks(trackMap); // Select the ID of the track we want to keep std::vector goodTracks = - m_duplicateClassifier.solveAmbiguity(cluster, tracks); + m_ambiML.solveAmbiguity(cluster, tracks); // Prepare the output track collection from the IDs - auto outputTracks = prepareOutputTrack(tracks, goodTracks); + TrackContainer solvedTracks{std::make_shared(), + std::make_shared()}; + solvedTracks.ensureDynamicColumns(tracks); + for (auto iTrack : goodTracks) { + auto destProxy = solvedTracks.makeTrack(); + auto srcProxy = tracks.getTrack(iTrack); + destProxy.copyFrom(srcProxy, false); + destProxy.tipIndex() = srcProxy.tipIndex(); + } + + ActsExamples::ConstTrackContainer outputTracks{ + std::make_shared( + std::move(solvedTracks.container())), + tracks.trackStateContainerHolder()}; + m_outputTracks(ctx, std::move(outputTracks)); return ActsExamples::ProcessCode::SUCCESS; diff --git a/Examples/Algorithms/TrackFindingML/src/AmbiguityResolutionMLDBScanAlgorithm.cpp b/Examples/Algorithms/TrackFindingML/src/AmbiguityResolutionMLDBScanAlgorithm.cpp deleted file mode 100644 index 6f0c7529208..00000000000 --- a/Examples/Algorithms/TrackFindingML/src/AmbiguityResolutionMLDBScanAlgorithm.cpp +++ /dev/null @@ -1,55 +0,0 @@ -// This file is part of the ACTS project. -// -// Copyright (C) 2016 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 https://mozilla.org/MPL/2.0/. - -#include "ActsExamples/TrackFindingML/AmbiguityResolutionMLDBScanAlgorithm.hpp" - -#include "ActsExamples/Framework/ProcessCode.hpp" -#include "ActsExamples/Framework/WhiteBoard.hpp" -#include "ActsExamples/TrackFindingML/AmbiguityDBScanClustering.hpp" - -#include -#include - -ActsExamples::AmbiguityResolutionMLDBScanAlgorithm:: - AmbiguityResolutionMLDBScanAlgorithm( - ActsExamples::AmbiguityResolutionMLDBScanAlgorithm::Config cfg, - Acts::Logging::Level lvl) - : ActsExamples::AmbiguityResolutionML( - "AmbiguityResolutionMLDBScanAlgorithm", lvl), - m_cfg(std::move(cfg)), - m_duplicateClassifier(m_cfg.inputDuplicateNN.c_str()) { - if (m_cfg.inputTracks.empty()) { - throw std::invalid_argument("Missing trajectories input collection"); - } - if (m_cfg.outputTracks.empty()) { - throw std::invalid_argument("Missing trajectories output collection"); - } - m_inputTracks.initialize(m_cfg.inputTracks); - m_outputTracks.initialize(m_cfg.outputTracks); -} - -ActsExamples::ProcessCode -ActsExamples::AmbiguityResolutionMLDBScanAlgorithm::execute( - const AlgorithmContext& ctx) const { - // Read input data - const auto& tracks = m_inputTracks(ctx); - // Associate measurement to their respective tracks - std::multimap>> - trackMap = mapTrackHits(tracks, m_cfg.nMeasurementsMin); - // Cluster the tracks using DBscan - auto cluster = Acts::dbscanTrackClustering( - trackMap, tracks, m_cfg.epsilonDBScan, m_cfg.minPointsDBScan); - // Select the ID of the track we want to keep - std::vector goodTracks = - m_duplicateClassifier.solveAmbiguity(cluster, tracks); - // Prepare the output track collection from the IDs - auto outputTracks = prepareOutputTrack(tracks, goodTracks); - m_outputTracks(ctx, std::move(outputTracks)); - - return ActsExamples::ProcessCode::SUCCESS; -} diff --git a/Examples/Algorithms/TrackFitting/include/ActsExamples/TrackFitting/SurfaceSortingAlgorithm.hpp b/Examples/Algorithms/TrackFitting/include/ActsExamples/TrackFitting/SurfaceSortingAlgorithm.hpp index 02f85bff2e1..6315e192e67 100644 --- a/Examples/Algorithms/TrackFitting/include/ActsExamples/TrackFitting/SurfaceSortingAlgorithm.hpp +++ b/Examples/Algorithms/TrackFitting/include/ActsExamples/TrackFitting/SurfaceSortingAlgorithm.hpp @@ -9,25 +9,16 @@ #pragma once #include "Acts/Utilities/Logger.hpp" -#include "ActsExamples/EventData/Index.hpp" -#include "ActsExamples/EventData/IndexSourceLink.hpp" #include "ActsExamples/EventData/Measurement.hpp" #include "ActsExamples/EventData/ProtoTrack.hpp" #include "ActsExamples/EventData/SimHit.hpp" -#include "ActsExamples/EventData/Track.hpp" #include "ActsExamples/Framework/DataHandle.hpp" #include "ActsExamples/Framework/IAlgorithm.hpp" #include "ActsExamples/Framework/ProcessCode.hpp" -#include -#include #include -#include namespace ActsExamples { -struct AlgorithmContext; - -using TrackHitList = std::map; class SurfaceSortingAlgorithm final : public IAlgorithm { public: @@ -55,7 +46,7 @@ class SurfaceSortingAlgorithm final : public IAlgorithm { ReadDataHandle m_inputProtoTracks{this, "InputProtoTracks"}; ReadDataHandle m_inputSimHits{this, "InputSimHits"}; - ReadDataHandle m_inputMeasurementSimHitsMap{ + ReadDataHandle m_inputMeasurementSimHitsMap{ this, "InputMeasurementSimHitsMap"}; WriteDataHandle m_outputProtoTracks{this, "OutputProtoTracks"}; diff --git a/Examples/Algorithms/TrackFitting/src/RefittingAlgorithm.cpp b/Examples/Algorithms/TrackFitting/src/RefittingAlgorithm.cpp index 03aa74fa58f..79eb367589b 100644 --- a/Examples/Algorithms/TrackFitting/src/RefittingAlgorithm.cpp +++ b/Examples/Algorithms/TrackFitting/src/RefittingAlgorithm.cpp @@ -64,7 +64,8 @@ ActsExamples::ProcessCode ActsExamples::RefittingAlgorithm::execute( auto itrack = 0ul; for (const auto& track : inputTracks) { // Check if you are not in picking mode - if (m_cfg.pickTrack > -1 && m_cfg.pickTrack != static_cast(itrack++)) { + if (m_cfg.pickTrack > -1 && + static_cast(m_cfg.pickTrack) != itrack++) { continue; } diff --git a/Examples/Algorithms/TrackFitting/src/SurfaceSortingAlgorithm.cpp b/Examples/Algorithms/TrackFitting/src/SurfaceSortingAlgorithm.cpp index 040634efbe9..57a79ab2db0 100644 --- a/Examples/Algorithms/TrackFitting/src/SurfaceSortingAlgorithm.cpp +++ b/Examples/Algorithms/TrackFitting/src/SurfaceSortingAlgorithm.cpp @@ -9,23 +9,18 @@ #include "ActsExamples/TrackFitting/SurfaceSortingAlgorithm.hpp" #include "ActsExamples/EventData/ProtoTrack.hpp" -#include "ActsExamples/EventData/SimHit.hpp" #include "ActsFatras/EventData/Hit.hpp" #include -#include #include #include #include namespace ActsExamples { -struct AlgorithmContext; -} // namespace ActsExamples -ActsExamples::SurfaceSortingAlgorithm::SurfaceSortingAlgorithm( - Config cfg, Acts::Logging::Level level) - : ActsExamples::IAlgorithm("SurfaceSortingAlgorithm", level), - m_cfg(std::move(cfg)) { +SurfaceSortingAlgorithm::SurfaceSortingAlgorithm(Config cfg, + Acts::Logging::Level level) + : IAlgorithm("SurfaceSortingAlgorithm", level), m_cfg(std::move(cfg)) { if (m_cfg.inputProtoTracks.empty()) { throw std::invalid_argument("Missing input proto track collection"); } @@ -45,15 +40,15 @@ ActsExamples::SurfaceSortingAlgorithm::SurfaceSortingAlgorithm( m_outputProtoTracks.initialize(m_cfg.outputProtoTracks); } -ActsExamples::ProcessCode ActsExamples::SurfaceSortingAlgorithm::execute( - const ActsExamples::AlgorithmContext& ctx) const { +ProcessCode SurfaceSortingAlgorithm::execute( + const AlgorithmContext& ctx) const { const auto& protoTracks = m_inputProtoTracks(ctx); const auto& simHits = m_inputSimHits(ctx); const auto& simHitsMap = m_inputMeasurementSimHitsMap(ctx); ProtoTrackContainer sortedTracks; sortedTracks.reserve(protoTracks.size()); - TrackHitList trackHitList; + std::map trackHitList; for (std::size_t itrack = 0; itrack < protoTracks.size(); ++itrack) { const auto& protoTrack = protoTracks[itrack]; @@ -83,5 +78,7 @@ ActsExamples::ProcessCode ActsExamples::SurfaceSortingAlgorithm::execute( m_outputProtoTracks(ctx, std::move(sortedTracks)); - return ActsExamples::ProcessCode::SUCCESS; + return ProcessCode::SUCCESS; } + +} // namespace ActsExamples diff --git a/Examples/Algorithms/TrackFitting/src/TrackFittingAlgorithm.cpp b/Examples/Algorithms/TrackFitting/src/TrackFittingAlgorithm.cpp index 50360e3d3b4..d79a864bc73 100644 --- a/Examples/Algorithms/TrackFitting/src/TrackFittingAlgorithm.cpp +++ b/Examples/Algorithms/TrackFitting/src/TrackFittingAlgorithm.cpp @@ -103,7 +103,8 @@ ActsExamples::ProcessCode ActsExamples::TrackFittingAlgorithm::execute( std::vector trackSourceLinks; for (std::size_t itrack = 0; itrack < protoTracks.size(); ++itrack) { // Check if you are not in picking mode - if (m_cfg.pickTrack > -1 && m_cfg.pickTrack != static_cast(itrack)) { + if (m_cfg.pickTrack > -1 && + static_cast(m_cfg.pickTrack) != itrack) { continue; } diff --git a/Examples/Algorithms/TruthTracking/ActsExamples/TruthTracking/TrackTruthMatcher.hpp b/Examples/Algorithms/TruthTracking/ActsExamples/TruthTracking/TrackTruthMatcher.hpp index 5f8d89a98b0..6c11d4e0c43 100644 --- a/Examples/Algorithms/TruthTracking/ActsExamples/TruthTracking/TrackTruthMatcher.hpp +++ b/Examples/Algorithms/TruthTracking/ActsExamples/TruthTracking/TrackTruthMatcher.hpp @@ -9,7 +9,7 @@ #pragma once #include "Acts/Utilities/Logger.hpp" -#include "ActsExamples/EventData/SimHit.hpp" +#include "ActsExamples/EventData/Measurement.hpp" #include "ActsExamples/EventData/SimParticle.hpp" #include "ActsExamples/EventData/Track.hpp" #include "ActsExamples/EventData/TruthMatching.hpp" @@ -21,8 +21,6 @@ namespace ActsExamples { -struct AlgorithmContext; - /// Matches tracks to truth particles and vice versa class TrackTruthMatcher final : public IAlgorithm { public: @@ -56,7 +54,7 @@ class TrackTruthMatcher final : public IAlgorithm { ReadDataHandle m_inputTracks{this, "InputTracks"}; ReadDataHandle m_inputParticles{this, "InputParticles"}; - ReadDataHandle m_inputMeasurementParticlesMap{ + ReadDataHandle m_inputMeasurementParticlesMap{ this, "InputMeasurementParticlesMap"}; WriteDataHandle m_outputTrackParticleMatching{ this, "OutputTrackParticleMatching"}; diff --git a/Examples/Algorithms/TruthTracking/ActsExamples/TruthTracking/TruthSeedingAlgorithm.cpp b/Examples/Algorithms/TruthTracking/ActsExamples/TruthTracking/TruthSeedingAlgorithm.cpp index 6ea0225e200..2e5d1c4cd15 100644 --- a/Examples/Algorithms/TruthTracking/ActsExamples/TruthTracking/TruthSeedingAlgorithm.cpp +++ b/Examples/Algorithms/TruthTracking/ActsExamples/TruthTracking/TruthSeedingAlgorithm.cpp @@ -12,7 +12,6 @@ #include "ActsExamples/EventData/IndexSourceLink.hpp" #include "ActsExamples/EventData/SimParticle.hpp" #include "ActsExamples/Utilities/Range.hpp" -#include "ActsFatras/EventData/Particle.hpp" #include #include @@ -25,18 +24,16 @@ #include namespace ActsExamples { -struct AlgorithmContext; -} // namespace ActsExamples -ActsExamples::TruthSeedingAlgorithm::TruthSeedingAlgorithm( - ActsExamples::TruthSeedingAlgorithm::Config cfg, Acts::Logging::Level lvl) - : ActsExamples::IAlgorithm("TruthSeedingAlgorithm", lvl), - m_cfg(std::move(cfg)) { +TruthSeedingAlgorithm::TruthSeedingAlgorithm(Config cfg, + Acts::Logging::Level lvl) + : IAlgorithm("TruthSeedingAlgorithm", lvl), m_cfg(std::move(cfg)) { if (m_cfg.inputParticles.empty()) { throw std::invalid_argument("Missing input truth particles collection"); } - if (m_cfg.inputMeasurementParticlesMap.empty()) { - throw std::invalid_argument("Missing input hit-particles map collection"); + if (m_cfg.inputParticleMeasurementsMap.empty()) { + throw std::invalid_argument( + "Missing input particle-measurements map collection"); } if (m_cfg.inputSpacePoints.empty()) { throw std::invalid_argument("Missing seeds or space point collection"); @@ -65,20 +62,16 @@ ActsExamples::TruthSeedingAlgorithm::TruthSeedingAlgorithm( } m_inputParticles.initialize(m_cfg.inputParticles); - m_inputMeasurementParticlesMap.initialize(m_cfg.inputMeasurementParticlesMap); + m_inputParticleMeasurementsMap.initialize(m_cfg.inputParticleMeasurementsMap); m_outputParticles.initialize(m_cfg.outputParticles); m_outputProtoTracks.initialize(m_cfg.outputProtoTracks); m_outputSeeds.initialize(m_cfg.outputSeeds); } -ActsExamples::ProcessCode ActsExamples::TruthSeedingAlgorithm::execute( - const ActsExamples::AlgorithmContext& ctx) const { +ProcessCode TruthSeedingAlgorithm::execute(const AlgorithmContext& ctx) const { // prepare input collections const auto& particles = m_inputParticles(ctx); - const auto& hitParticlesMap = m_inputMeasurementParticlesMap(ctx); - // compute particle_id -> {hit_id...} map from the - // hit_id -> {particle_id...} map on the fly. - const auto& particleHitsMap = invertIndexMultimap(hitParticlesMap); + const auto& particleMeasurementsMap = m_inputParticleMeasurementsMap(ctx); // construct the combined input container of space point pointers from all // configured input sources. @@ -120,27 +113,28 @@ ActsExamples::ProcessCode ActsExamples::TruthSeedingAlgorithm::execute( } for (const auto& particle : particles) { - // find the corresponding hits for this particle - const auto& hits = - makeRange(particleHitsMap.equal_range(particle.particleId())); - // fill hit indices to create the proto track + // find the corresponding measurements for this particle + const auto& measurements = + makeRange(particleMeasurementsMap.equal_range(particle.particleId())); + // fill measurement indices to create the proto track ProtoTrack track; - track.reserve(hits.size()); - for (const auto& hit : hits) { - track.push_back(hit.second); + track.reserve(measurements.size()); + for (const auto& measurement : measurements) { + track.push_back(measurement.second); } - // The list of hits and the initial start parameters + // The list of measurements and the initial start parameters if (track.size() < 3) { - ACTS_WARNING("Particle " << particle << " has less than 3 hits"); + ACTS_WARNING("Particle " << particle << " has less than 3 measurements"); continue; } // Space points on the proto track std::vector spacePointsOnTrack; spacePointsOnTrack.reserve(track.size()); - // Loop over the hit index on the proto track to find the space points - for (const auto& hitIndex : track) { - auto it = spMap.find(hitIndex); + // Loop over the measurement index on the proto track to find the space + // points + for (const auto& measurementIndex : track) { + auto it = spMap.find(measurementIndex); if (it != spMap.end()) { spacePointsOnTrack.push_back(it->second); } @@ -198,5 +192,7 @@ ActsExamples::ProcessCode ActsExamples::TruthSeedingAlgorithm::execute( m_outputProtoTracks(ctx, std::move(tracks)); m_outputSeeds(ctx, std::move(seeds)); - return ActsExamples::ProcessCode::SUCCESS; + return ProcessCode::SUCCESS; } + +} // namespace ActsExamples diff --git a/Examples/Algorithms/TruthTracking/ActsExamples/TruthTracking/TruthSeedingAlgorithm.hpp b/Examples/Algorithms/TruthTracking/ActsExamples/TruthTracking/TruthSeedingAlgorithm.hpp index a6cb0416160..5807ef54cd0 100644 --- a/Examples/Algorithms/TruthTracking/ActsExamples/TruthTracking/TruthSeedingAlgorithm.hpp +++ b/Examples/Algorithms/TruthTracking/ActsExamples/TruthTracking/TruthSeedingAlgorithm.hpp @@ -11,7 +11,6 @@ #include "Acts/Definitions/Units.hpp" #include "Acts/Utilities/Logger.hpp" #include "ActsExamples/EventData/ProtoTrack.hpp" -#include "ActsExamples/EventData/SimHit.hpp" #include "ActsExamples/EventData/SimParticle.hpp" #include "ActsExamples/EventData/SimSeed.hpp" #include "ActsExamples/EventData/SimSpacePoint.hpp" @@ -23,16 +22,7 @@ #include #include -namespace ActsFatras { -class Barcode; -} // namespace ActsFatras - -namespace Acts { -class TrackingGeometry; -} - namespace ActsExamples { -struct AlgorithmContext; /// Construct track seeds from particles. class TruthSeedingAlgorithm final : public IAlgorithm { @@ -40,8 +30,8 @@ class TruthSeedingAlgorithm final : public IAlgorithm { struct Config { /// The input truth particles that should be used for truth seeding. std::string inputParticles; - /// The input hit-particles map collection. - std::string inputMeasurementParticlesMap; + /// The input particle-measurements map collection. + std::string inputParticleMeasurementsMap; /// Input space point collections. /// /// We allow multiple space point collections to allow different parts of @@ -80,8 +70,8 @@ class TruthSeedingAlgorithm final : public IAlgorithm { Config m_cfg; ReadDataHandle m_inputParticles{this, "InputParticles"}; - ReadDataHandle m_inputMeasurementParticlesMap{ - this, "InputMeasurementParticlesMaps"}; + ReadDataHandle> m_inputParticleMeasurementsMap{ + this, "InputParticleMeasurementsMap"}; std::vector>> m_inputSpacePoints{}; diff --git a/Examples/Algorithms/TruthTracking/ActsExamples/TruthTracking/TruthTrackFinder.cpp b/Examples/Algorithms/TruthTracking/ActsExamples/TruthTracking/TruthTrackFinder.cpp index 54953ce7c2b..84b08693635 100644 --- a/Examples/Algorithms/TruthTracking/ActsExamples/TruthTracking/TruthTrackFinder.cpp +++ b/Examples/Algorithms/TruthTracking/ActsExamples/TruthTracking/TruthTrackFinder.cpp @@ -8,23 +8,15 @@ #include "ActsExamples/TruthTracking/TruthTrackFinder.hpp" -#include "Acts/Utilities/MultiIndex.hpp" -#include "ActsExamples/EventData/Index.hpp" #include "ActsExamples/EventData/ProtoTrack.hpp" #include "ActsExamples/EventData/SimParticle.hpp" #include "ActsExamples/Utilities/Range.hpp" -#include "ActsFatras/EventData/Particle.hpp" -#include #include #include #include namespace ActsExamples { -struct AlgorithmContext; -} // namespace ActsExamples - -using namespace ActsExamples; TruthTrackFinder::TruthTrackFinder(const Config& config, Acts::Logging::Level level) @@ -32,7 +24,7 @@ TruthTrackFinder::TruthTrackFinder(const Config& config, if (m_cfg.inputParticles.empty()) { throw std::invalid_argument("Missing input truth particles collection"); } - if (m_cfg.inputMeasurementParticlesMap.empty()) { + if (m_cfg.inputParticleMeasurementsMap.empty()) { throw std::invalid_argument("Missing input hit-particles map collection"); } if (m_cfg.outputProtoTracks.empty()) { @@ -40,17 +32,14 @@ TruthTrackFinder::TruthTrackFinder(const Config& config, } m_inputParticles.initialize(m_cfg.inputParticles); - m_inputMeasurementParticlesMap.initialize(m_cfg.inputMeasurementParticlesMap); + m_inputParticleMeasurementsMap.initialize(m_cfg.inputParticleMeasurementsMap); m_outputProtoTracks.initialize(m_cfg.outputProtoTracks); } ProcessCode TruthTrackFinder::execute(const AlgorithmContext& ctx) const { // prepare input collections const auto& particles = m_inputParticles(ctx); - const auto& hitParticlesMap = m_inputMeasurementParticlesMap(ctx); - // compute particle_id -> {hit_id...} map from the - // hit_id -> {particle_id...} map on the fly. - const auto& particleHitsMap = invertIndexMultimap(hitParticlesMap); + const auto& particleMeasurementsMap = m_inputParticleMeasurementsMap(ctx); // prepare output collection ProtoTrackContainer tracks; @@ -59,14 +48,15 @@ ProcessCode TruthTrackFinder::execute(const AlgorithmContext& ctx) const { ACTS_VERBOSE("Create prototracks for " << particles.size() << " particles"); for (const auto& particle : particles) { // find the corresponding hits for this particle - const auto& hits = - makeRange(particleHitsMap.equal_range(particle.particleId())); - ACTS_VERBOSE(" - Prototrack from " << hits.size() << " hits"); + const auto& measurements = + makeRange(particleMeasurementsMap.equal_range(particle.particleId())); + ACTS_VERBOSE(" - Prototrack from " << measurements.size() + << " measurements"); // fill hit indices to create the proto track ProtoTrack track; - track.reserve(hits.size()); - for (const auto& hit : hits) { - track.emplace_back(hit.second); + track.reserve(measurements.size()); + for (const auto& measurement : measurements) { + track.emplace_back(measurement.second); } // add proto track to the output collection tracks.emplace_back(std::move(track)); @@ -75,3 +65,5 @@ ProcessCode TruthTrackFinder::execute(const AlgorithmContext& ctx) const { m_outputProtoTracks(ctx, std::move(tracks)); return ProcessCode::SUCCESS; } + +} // namespace ActsExamples diff --git a/Examples/Algorithms/TruthTracking/ActsExamples/TruthTracking/TruthTrackFinder.hpp b/Examples/Algorithms/TruthTracking/ActsExamples/TruthTracking/TruthTrackFinder.hpp index 8f7affb79c3..39c42f651b3 100644 --- a/Examples/Algorithms/TruthTracking/ActsExamples/TruthTracking/TruthTrackFinder.hpp +++ b/Examples/Algorithms/TruthTracking/ActsExamples/TruthTracking/TruthTrackFinder.hpp @@ -10,7 +10,6 @@ #include "Acts/Utilities/Logger.hpp" #include "ActsExamples/EventData/ProtoTrack.hpp" -#include "ActsExamples/EventData/SimHit.hpp" #include "ActsExamples/EventData/SimParticle.hpp" #include "ActsExamples/Framework/DataHandle.hpp" #include "ActsExamples/Framework/IAlgorithm.hpp" @@ -18,12 +17,7 @@ #include -namespace ActsFatras { -class Barcode; -} // namespace ActsFatras - namespace ActsExamples { -struct AlgorithmContext; /// Convert true particle tracks into "reconstructed" proto tracks. /// @@ -37,8 +31,8 @@ class TruthTrackFinder final : public IAlgorithm { struct Config { /// The input truth particles that should be used to create proto tracks. std::string inputParticles; - /// The input hit-particles map collection. - std::string inputMeasurementParticlesMap; + /// The input particle-measurements map collection. + std::string inputParticleMeasurementsMap; /// The output proto tracks collection. std::string outputProtoTracks; }; @@ -55,8 +49,8 @@ class TruthTrackFinder final : public IAlgorithm { ReadDataHandle m_inputParticles{this, "InputParticles"}; - ReadDataHandle m_inputMeasurementParticlesMap{ - this, "InputMeasurementParticlesMap"}; + ReadDataHandle> m_inputParticleMeasurementsMap{ + this, "InputParticleMeasurementsMap"}; WriteDataHandle m_outputProtoTracks{this, "OutputProtoTracks"}; diff --git a/Examples/Algorithms/TruthTracking/ActsExamples/TruthTracking/TruthVertexFinder.hpp b/Examples/Algorithms/TruthTracking/ActsExamples/TruthTracking/TruthVertexFinder.hpp index 2ea7d401866..bbbefa2dd96 100644 --- a/Examples/Algorithms/TruthTracking/ActsExamples/TruthTracking/TruthVertexFinder.hpp +++ b/Examples/Algorithms/TruthTracking/ActsExamples/TruthTracking/TruthVertexFinder.hpp @@ -9,6 +9,7 @@ #pragma once #include "Acts/Utilities/Logger.hpp" +#include "ActsExamples/EventData/Measurement.hpp" #include "ActsExamples/EventData/ProtoVertex.hpp" #include "ActsExamples/EventData/SimParticle.hpp" #include "ActsExamples/EventData/Track.hpp" @@ -19,13 +20,10 @@ #include namespace ActsExamples { -struct AlgorithmContext; /// Group particles into proto vertices using truth information. class TruthVertexFinder final : public IAlgorithm { public: - using HitParticlesMap = ActsExamples::IndexMultimap; - struct Config { /// The input tracks that should be used to create proto vertices. std::string inputTracks; @@ -55,7 +53,7 @@ class TruthVertexFinder final : public IAlgorithm { ReadDataHandle m_inputTracks{this, "InputTracks"}; ReadDataHandle m_inputParticles{this, "InputParticles"}; - ReadDataHandle m_inputMeasurementParticlesMap{ + ReadDataHandle m_inputMeasurementParticlesMap{ this, "InputMeasurementParticlesMap"}; WriteDataHandle m_outputProtoVertices{ this, "OutputProtoVertices"}; diff --git a/Examples/Detectors/ContextualDetector/src/AlignedDetector.cpp b/Examples/Detectors/ContextualDetector/src/AlignedDetector.cpp index 84e050c8470..9a9c5901543 100644 --- a/Examples/Detectors/ContextualDetector/src/AlignedDetector.cpp +++ b/Examples/Detectors/ContextualDetector/src/AlignedDetector.cpp @@ -106,8 +106,7 @@ auto AlignedDetector::finalize( } // return the pair of geometry and the alignment decorator(s) - return std::make_pair( - std::move(aTrackingGeometry), std::move(aContextDecorators)); + return {std::move(aTrackingGeometry), std::move(aContextDecorators)}; } } // namespace ActsExamples diff --git a/Examples/Detectors/DD4hepDetector/src/DD4hepDetector.cpp b/Examples/Detectors/DD4hepDetector/src/DD4hepDetector.cpp index 4ecae308134..a2c12dfaaa7 100644 --- a/Examples/Detectors/DD4hepDetector/src/DD4hepDetector.cpp +++ b/Examples/Detectors/DD4hepDetector/src/DD4hepDetector.cpp @@ -42,8 +42,7 @@ auto DD4hepDetector::finalize( } ContextDecorators dd4ContextDecorators = {}; // return the pair of geometry and empty decorators - return std::make_pair( - std::move(dd4tGeometry), std::move(dd4ContextDecorators)); + return {std::move(dd4tGeometry), std::move(dd4ContextDecorators)}; } auto DD4hepDetector::finalize( diff --git a/Examples/Detectors/GenericDetector/src/GenericDetector.cpp b/Examples/Detectors/GenericDetector/src/GenericDetector.cpp index bf033fc1e3e..b80b769708d 100644 --- a/Examples/Detectors/GenericDetector/src/GenericDetector.cpp +++ b/Examples/Detectors/GenericDetector/src/GenericDetector.cpp @@ -27,8 +27,7 @@ auto GenericDetector::finalize( cfg.volumeLogLevel); ContextDecorators gContextDecorators = {}; // return the pair of geometry and empty decorators - return std::make_pair( - std::move(gGeometry), std::move(gContextDecorators)); + return {std::move(gGeometry), std::move(gContextDecorators)}; } } // namespace ActsExamples diff --git a/Examples/Detectors/TGeoDetector/src/TGeoDetector.cpp b/Examples/Detectors/TGeoDetector/src/TGeoDetector.cpp index 0ece428a438..f594b681d49 100644 --- a/Examples/Detectors/TGeoDetector/src/TGeoDetector.cpp +++ b/Examples/Detectors/TGeoDetector/src/TGeoDetector.cpp @@ -375,8 +375,7 @@ auto TGeoDetector::finalize( ContextDecorators tgeoContextDecorators = {}; // Return the pair of geometry and empty decorators - return std::make_pair( - std::move(tgeoTrackingGeometry), std::move(tgeoContextDecorators)); + return {std::move(tgeoTrackingGeometry), std::move(tgeoContextDecorators)}; } void TGeoDetector::Config::readJson(const std::string& jsonFile) { diff --git a/Examples/Detectors/TelescopeDetector/src/TelescopeDetector.cpp b/Examples/Detectors/TelescopeDetector/src/TelescopeDetector.cpp index bbf1aff4106..ccba63a1c34 100644 --- a/Examples/Detectors/TelescopeDetector/src/TelescopeDetector.cpp +++ b/Examples/Detectors/TelescopeDetector/src/TelescopeDetector.cpp @@ -61,8 +61,7 @@ auto TelescopeDetector::finalize( static_cast(cfg.binValue)); ContextDecorators gContextDecorators = {}; // return the pair of geometry and empty decorators - return std::make_pair( - std::move(gGeometry), std::move(gContextDecorators)); + return {std::move(gGeometry), std::move(gContextDecorators)}; } } // namespace ActsExamples diff --git a/Examples/Framework/include/ActsExamples/EventData/Index.hpp b/Examples/Framework/include/ActsExamples/EventData/Index.hpp index 46e431ea76a..5750737b05f 100644 --- a/Examples/Framework/include/ActsExamples/EventData/Index.hpp +++ b/Examples/Framework/include/ActsExamples/EventData/Index.hpp @@ -30,18 +30,21 @@ using Index = std::uint32_t; template using IndexMultimap = boost::container::flat_multimap; -/// Invert the multimap, i.e. from a -> {b...} to b -> {a...}. +/// Store the inverse of an index multimap, i.e. from a -> {b...} to b -> +/// {a...}. /// /// @note This assumes that the value in the initial multimap is itself a /// sortable index-like object, as would be the case when mapping e.g. /// hit ids to particle ids/ barcodes. template -inline boost::container::flat_multimap invertIndexMultimap( - const IndexMultimap& multimap) { - using InverseMultimap = boost::container::flat_multimap; +using InverseMultimap = boost::container::flat_multimap; +/// Invert the multimap, i.e. from a -> {b...} to b -> {a...} +template +inline InverseMultimap invertIndexMultimap( + const IndexMultimap& multimap) { // switch key-value without enforcing the new ordering (linear copy) - typename InverseMultimap::sequence_type unordered; + typename InverseMultimap::sequence_type unordered; unordered.reserve(multimap.size()); for (auto&& [index, value] : multimap) { // value is now the key and the index is now the value @@ -49,7 +52,7 @@ inline boost::container::flat_multimap invertIndexMultimap( } // adopting the unordered sequence will reestablish the correct order - InverseMultimap inverse; + InverseMultimap inverse; #if BOOST_VERSION < 107800 for (const auto& i : unordered) { inverse.insert(i); diff --git a/Examples/Framework/include/ActsExamples/EventData/Measurement.hpp b/Examples/Framework/include/ActsExamples/EventData/Measurement.hpp index 251d7bc293a..86a8be27aba 100644 --- a/Examples/Framework/include/ActsExamples/EventData/Measurement.hpp +++ b/Examples/Framework/include/ActsExamples/EventData/Measurement.hpp @@ -17,6 +17,7 @@ #include "ActsExamples/EventData/GeometryContainers.hpp" #include "ActsExamples/EventData/IndexSourceLink.hpp" #include "ActsExamples/EventData/MeasurementConcept.hpp" +#include "ActsExamples/EventData/SimParticle.hpp" #include #include @@ -531,4 +532,10 @@ static_assert( std::random_access_iterator && std::random_access_iterator); +using MeasurementSimHitsMap = IndexMultimap; +using MeasurementParticlesMap = IndexMultimap; + +using SimHitMeasurementsMap = InverseMultimap; +using ParticleMeasurementsMap = InverseMultimap; + } // namespace ActsExamples diff --git a/Examples/Framework/include/ActsExamples/EventData/SimHit.hpp b/Examples/Framework/include/ActsExamples/EventData/SimHit.hpp index be10ce2639c..8c7854082fd 100644 --- a/Examples/Framework/include/ActsExamples/EventData/SimHit.hpp +++ b/Examples/Framework/include/ActsExamples/EventData/SimHit.hpp @@ -9,8 +9,6 @@ #pragma once #include "ActsExamples/EventData/GeometryContainers.hpp" -#include "ActsExamples/EventData/Index.hpp" -#include "ActsExamples/EventData/SimParticle.hpp" #include "ActsFatras/EventData/Hit.hpp" namespace ActsExamples { @@ -19,8 +17,4 @@ using SimHit = ::ActsFatras::Hit; /// Store hits ordered by geometry identifier. using SimHitContainer = GeometryIdMultiset; -using HitParticlesMap = IndexMultimap; - -using HitSimHitsMap = IndexMultimap; - } // namespace ActsExamples diff --git a/Examples/Framework/src/EventData/ScalingCalibrator.cpp b/Examples/Framework/src/EventData/ScalingCalibrator.cpp index 50e693e8b15..068f6221a25 100644 --- a/Examples/Framework/src/EventData/ScalingCalibrator.cpp +++ b/Examples/Framework/src/EventData/ScalingCalibrator.cpp @@ -52,22 +52,22 @@ std::pair parseMapKey( std::regex reg("^map_([0-9]+)-([0-9]+)-([0-9]+)_([xy]_.*)$"); std::smatch matches; - if (std::regex_search(mapkey, matches, reg) && matches.size() == 5) { - std::size_t vol = std::stoull(matches[1].str()); - std::size_t lyr = std::stoull(matches[2].str()); - std::size_t mod = std::stoull(matches[3].str()); + if (!std::regex_search(mapkey, matches, reg) || matches.size() != 5) { + throw std::runtime_error("Invalid map key: " + mapkey); + } - Acts::GeometryIdentifier geoId; - geoId.setVolume(vol); - geoId.setLayer(lyr); - geoId.setSensitive(mod); + std::size_t vol = std::stoull(matches[1].str()); + std::size_t lyr = std::stoull(matches[2].str()); + std::size_t mod = std::stoull(matches[3].str()); - std::string var(matches[4].str()); + Acts::GeometryIdentifier geoId; + geoId.setVolume(vol); + geoId.setLayer(lyr); + geoId.setSensitive(mod); - return std::make_pair(geoId, var); - } else { - throw std::runtime_error("Invalid map key: " + mapkey); - } + std::string var(matches[4].str()); + + return {geoId, var}; } std::map diff --git a/Examples/Framework/src/Utilities/Paths.cpp b/Examples/Framework/src/Utilities/Paths.cpp index 02bf5889c68..4de064d6eb0 100644 --- a/Examples/Framework/src/Utilities/Paths.cpp +++ b/Examples/Framework/src/Utilities/Paths.cpp @@ -110,7 +110,7 @@ std::pair ActsExamples::determineEventFilesRange( // should only occur if no files matched and the initial values persisted. if (eventMax < eventMin) { - return std::make_pair(0u, 0u); + return {0u, 0u}; } - return std::make_pair(eventMin, eventMax + 1); + return {eventMin, eventMax + 1}; } diff --git a/Examples/Io/Csv/include/ActsExamples/Io/Csv/CsvSeedWriter.hpp b/Examples/Io/Csv/include/ActsExamples/Io/Csv/CsvSeedWriter.hpp index 670183f7c29..d893c5aae42 100644 --- a/Examples/Io/Csv/include/ActsExamples/Io/Csv/CsvSeedWriter.hpp +++ b/Examples/Io/Csv/include/ActsExamples/Io/Csv/CsvSeedWriter.hpp @@ -8,6 +8,7 @@ #pragma once +#include "ActsExamples/EventData/Measurement.hpp" #include "ActsExamples/EventData/ProtoTrack.hpp" #include "ActsExamples/EventData/SimHit.hpp" #include "ActsExamples/EventData/SimParticle.hpp" @@ -76,9 +77,9 @@ class CsvSeedWriter : public WriterT { ReadDataHandle m_inputParticles{this, "InputParticles"}; ReadDataHandle m_inputSimSeeds{this, "InputSimSeeds"}; ReadDataHandle m_inputSimHits{this, "InputSimHits"}; - ReadDataHandle m_inputMeasurementParticlesMap{ + ReadDataHandle m_inputMeasurementParticlesMap{ this, "InputMeasurementParticlesMap"}; - ReadDataHandle m_inputMeasurementSimHitsMap{ + ReadDataHandle m_inputMeasurementSimHitsMap{ this, "InputMeasurementSimHitsMap"}; /// @brief Struct for brief seed summary info @@ -95,7 +96,7 @@ class CsvSeedWriter : public WriterT { float truthDistance = -1; std::string seedType = "unknown"; ProtoTrack measurementsID; - }; // trackInfo struct + }; }; } // namespace ActsExamples diff --git a/Examples/Io/Csv/include/ActsExamples/Io/Csv/CsvSpacepointWriter.hpp b/Examples/Io/Csv/include/ActsExamples/Io/Csv/CsvSpacePointWriter.hpp similarity index 93% rename from Examples/Io/Csv/include/ActsExamples/Io/Csv/CsvSpacepointWriter.hpp rename to Examples/Io/Csv/include/ActsExamples/Io/Csv/CsvSpacePointWriter.hpp index 1e24eb4bc97..ea7315d4927 100644 --- a/Examples/Io/Csv/include/ActsExamples/Io/Csv/CsvSpacepointWriter.hpp +++ b/Examples/Io/Csv/include/ActsExamples/Io/Csv/CsvSpacePointWriter.hpp @@ -34,7 +34,7 @@ namespace ActsExamples { /// ... /// /// Intrinsically thread-safe as one file per event. -class CsvSpacepointWriter final : public WriterT { +class CsvSpacePointWriter final : public WriterT { public: struct Config { /// Which measurement collection to write. @@ -48,10 +48,10 @@ class CsvSpacepointWriter final : public WriterT { /// Constructor with /// @param config configuration struct /// @param level logging level - CsvSpacepointWriter(const Config& config, Acts::Logging::Level level); + CsvSpacePointWriter(const Config& config, Acts::Logging::Level level); /// Virtual destructor - ~CsvSpacepointWriter() override; + ~CsvSpacePointWriter() override; /// End-of-run hook ProcessCode finalize() override; diff --git a/Examples/Io/Csv/include/ActsExamples/Io/Csv/CsvTrackWriter.hpp b/Examples/Io/Csv/include/ActsExamples/Io/Csv/CsvTrackWriter.hpp index 06a449e2615..d3bbfe989f7 100644 --- a/Examples/Io/Csv/include/ActsExamples/Io/Csv/CsvTrackWriter.hpp +++ b/Examples/Io/Csv/include/ActsExamples/Io/Csv/CsvTrackWriter.hpp @@ -11,7 +11,7 @@ #include "Acts/Definitions/Units.hpp" #include "Acts/EventData/MultiTrajectoryHelpers.hpp" #include "Acts/Utilities/Logger.hpp" -#include "ActsExamples/EventData/SimHit.hpp" +#include "ActsExamples/EventData/Measurement.hpp" #include "ActsExamples/EventData/Track.hpp" #include "ActsExamples/Framework/ProcessCode.hpp" #include "ActsExamples/Framework/WriterT.hpp" @@ -22,12 +22,6 @@ #include #include -namespace ActsExamples { -struct AlgorithmContext; -} // namespace ActsExamples - -using namespace Acts::UnitLiterals; - namespace ActsExamples { /// @class CsvTrackWriter @@ -46,16 +40,24 @@ namespace ActsExamples { class CsvTrackWriter : public WriterT { public: struct Config { - std::string inputTracks; ///< Input track collection - std::string outputDir; ///< where to place output files - std::string fileName = "CKFtracks.csv"; ///< name of the output files - std::string - inputMeasurementParticlesMap; ///< Input hit-particles map collection - std::size_t outputPrecision = 6; ///< floating point precision - std::size_t nMeasurementsMin = 7; ///< Min number of measurements - bool onlyTruthMatched = false; ///< Only write truth matched tracks - double truthMatchProbMin = 0.5; ///< Probability threshold for fake tracks - double ptMin = 1_GeV; ///< Min pt of tracks + /// Input track collection + std::string inputTracks; + /// where to place output files + std::string outputDir; + /// name of the output files + std::string fileName = "CKFtracks.csv"; + /// Input hit-particles map collection + std::string inputMeasurementParticlesMap; + /// floating point precision + std::size_t outputPrecision = 6; + /// Min number of measurements + std::size_t nMeasurementsMin = 7; + /// Only write truth matched tracks + bool onlyTruthMatched = false; + /// Probability threshold for fake tracks + double truthMatchProbMin = 0.5; + /// Min pt of tracks + double ptMin = 1 * Acts::UnitConstants::GeV; }; /// constructor @@ -75,9 +77,10 @@ class CsvTrackWriter : public WriterT { const ConstTrackContainer& tracks) override; private: - Config m_cfg; //!< Nested configuration struct + /// Nested configuration struct + Config m_cfg; - ReadDataHandle m_inputMeasurementParticlesMap{ + ReadDataHandle m_inputMeasurementParticlesMap{ this, "InputMeasurementParticlesMap"}; /// @brief Struct for brief trajectory summary info @@ -91,7 +94,7 @@ class CsvTrackWriter : public WriterT { double truthMatchProb = 0; std::optional fittedParameters; std::vector measurementsID; - }; // TrackInfo struct + }; }; } // namespace ActsExamples diff --git a/Examples/Io/Csv/src/CsvSpacePointWriter.cpp b/Examples/Io/Csv/src/CsvSpacePointWriter.cpp index 9eb440f7692..4e71020a390 100644 --- a/Examples/Io/Csv/src/CsvSpacePointWriter.cpp +++ b/Examples/Io/Csv/src/CsvSpacePointWriter.cpp @@ -6,7 +6,7 @@ // License, v. 2.0. If a copy of the MPL was not distributed with this // file, You can obtain one at https://mozilla.org/MPL/2.0/. -#include "ActsExamples/Io/Csv/CsvSpacepointWriter.hpp" +#include "ActsExamples/Io/Csv/CsvSpacePointWriter.hpp" #include "Acts/EventData/SourceLink.hpp" #include "Acts/Geometry/GeometryIdentifier.hpp" @@ -24,20 +24,20 @@ #include "CsvOutputData.hpp" -ActsExamples::CsvSpacepointWriter::CsvSpacepointWriter( - const ActsExamples::CsvSpacepointWriter::Config& config, +ActsExamples::CsvSpacePointWriter::CsvSpacePointWriter( + const ActsExamples::CsvSpacePointWriter::Config& config, Acts::Logging::Level level) - : WriterT(config.inputSpacepoints, "CsvSpacepointWriter", level), + : WriterT(config.inputSpacepoints, "CsvSpacePointWriter", level), m_cfg(config) {} -ActsExamples::CsvSpacepointWriter::~CsvSpacepointWriter() = default; +ActsExamples::CsvSpacePointWriter::~CsvSpacePointWriter() = default; -ActsExamples::ProcessCode ActsExamples::CsvSpacepointWriter::finalize() { +ActsExamples::ProcessCode ActsExamples::CsvSpacePointWriter::finalize() { // Write the tree return ProcessCode::SUCCESS; } -ActsExamples::ProcessCode ActsExamples::CsvSpacepointWriter::writeT( +ActsExamples::ProcessCode ActsExamples::CsvSpacePointWriter::writeT( const AlgorithmContext& ctx, const SimSpacePointContainer& spacepoints) { // Open per-event file for all components std::string pathSP = diff --git a/Examples/Io/Obj/CMakeLists.txt b/Examples/Io/Obj/CMakeLists.txt index 9ca888f726a..bce7bf19e6b 100644 --- a/Examples/Io/Obj/CMakeLists.txt +++ b/Examples/Io/Obj/CMakeLists.txt @@ -2,6 +2,7 @@ add_library( ActsExamplesIoObj SHARED src/ObjTrackingGeometryWriter.cpp + src/ObjSimHitWriter.cpp src/ObjPropagationStepsWriter.cpp ) diff --git a/Examples/Io/Obj/include/ActsExamples/Plugins/Obj/ObjPropagationStepsWriter.hpp b/Examples/Io/Obj/include/ActsExamples/Io/Obj/ObjPropagationStepsWriter.hpp similarity index 100% rename from Examples/Io/Obj/include/ActsExamples/Plugins/Obj/ObjPropagationStepsWriter.hpp rename to Examples/Io/Obj/include/ActsExamples/Io/Obj/ObjPropagationStepsWriter.hpp diff --git a/Examples/Io/Obj/include/ActsExamples/Io/Obj/ObjSimHitWriter.hpp b/Examples/Io/Obj/include/ActsExamples/Io/Obj/ObjSimHitWriter.hpp new file mode 100644 index 00000000000..b510bd3c035 --- /dev/null +++ b/Examples/Io/Obj/include/ActsExamples/Io/Obj/ObjSimHitWriter.hpp @@ -0,0 +1,90 @@ +// This file is part of the ACTS project. +// +// Copyright (C) 2016 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 https://mozilla.org/MPL/2.0/. + +#pragma once + +#include "Acts/Definitions/Units.hpp" +#include "Acts/Utilities/Logger.hpp" +#include "ActsExamples/EventData/SimHit.hpp" +#include "ActsExamples/Framework/ProcessCode.hpp" +#include "ActsExamples/Framework/WriterT.hpp" + +#include +#include +#include + +namespace ActsExamples { +struct AlgorithmContext; + +/// Write out a simhit collection before detector digitization as wavefront obj +/// file(s per event). +/// +/// This writes two files per event, one for the hits one for the trajectory. +/// The latter can be smoothed using a spline interpolation. +/// +/// event000000001-.obj +/// event000000001-_trajectory.obj +/// event000000002-.obj +/// event000000002-_trajectory.obj +/// +/// +/// The trajectory can be smoothed using a spline interpolation, where +/// nInterpolatedPoints points are added between each hit. +class ObjSimHitWriter : public WriterT { + public: + struct Config { + /// Input sim hit collection to write. + std::string inputSimHits; + /// Where to place output files + std::string outputDir; + /// Output filename stem. + std::string outputStem = "simhits"; + /// Number of decimal digits for floating point precision in output. + std::size_t outputPrecision = std::numeric_limits::max_digits10; + /// Draw line connections between hits + bool drawConnections = true; + /// Momentum threshold for hits + double momentumThreshold = 0.05 * Acts::UnitConstants::GeV; + /// Momentum threshold for trajectories + double momentumThresholdTraj = 0.05 * Acts::UnitConstants::GeV; + /// Number of points to interpolated between hits to smooth the + /// trajectory view in the obj file. + std::size_t nInterpolatedPoints = 4; + /// Keep the original hits in the trajectory file + bool keepOriginalHits = false; + }; + + /// Construct the particle writer. + /// + /// @param config is the configuration object + /// @param level is the logging level + ObjSimHitWriter(const Config& config, Acts::Logging::Level level); + + /// Ensure underlying file is closed. + ~ObjSimHitWriter() override = default; + + /// End-of-run hook + ProcessCode finalize() override; + + /// Get readonly access to the config parameters + const Config& config() const { return m_cfg; } + + protected: + /// Type-specific write implementation. + /// + /// @param[in] ctx is the algorithm context + /// @param[in] hits are the hits to be written + ProcessCode writeT(const AlgorithmContext& ctx, + const SimHitContainer& hits) override; + + private: + Config m_cfg; + std::mutex m_writeMutex; +}; + +} // namespace ActsExamples diff --git a/Examples/Io/Obj/include/ActsExamples/Plugins/Obj/ObjTrackingGeometryWriter.hpp b/Examples/Io/Obj/include/ActsExamples/Io/Obj/ObjTrackingGeometryWriter.hpp similarity index 100% rename from Examples/Io/Obj/include/ActsExamples/Plugins/Obj/ObjTrackingGeometryWriter.hpp rename to Examples/Io/Obj/include/ActsExamples/Io/Obj/ObjTrackingGeometryWriter.hpp diff --git a/Examples/Io/Obj/include/ActsExamples/Plugins/Obj/ObjWriterOptions.hpp b/Examples/Io/Obj/include/ActsExamples/Io/Obj/ObjWriterOptions.hpp similarity index 100% rename from Examples/Io/Obj/include/ActsExamples/Plugins/Obj/ObjWriterOptions.hpp rename to Examples/Io/Obj/include/ActsExamples/Io/Obj/ObjWriterOptions.hpp diff --git a/Examples/Io/Obj/src/ObjPropagationStepsWriter.cpp b/Examples/Io/Obj/src/ObjPropagationStepsWriter.cpp index 6d0388984cc..0f196ce162c 100644 --- a/Examples/Io/Obj/src/ObjPropagationStepsWriter.cpp +++ b/Examples/Io/Obj/src/ObjPropagationStepsWriter.cpp @@ -6,7 +6,7 @@ // License, v. 2.0. If a copy of the MPL was not distributed with this // file, You can obtain one at https://mozilla.org/MPL/2.0/. -#include "ActsExamples/Plugins/Obj/ObjPropagationStepsWriter.hpp" +#include "ActsExamples/Io/Obj/ObjPropagationStepsWriter.hpp" namespace ActsExamples { diff --git a/Examples/Io/Obj/src/ObjSimHitWriter.cpp b/Examples/Io/Obj/src/ObjSimHitWriter.cpp new file mode 100644 index 00000000000..e5552b627e2 --- /dev/null +++ b/Examples/Io/Obj/src/ObjSimHitWriter.cpp @@ -0,0 +1,152 @@ +// This file is part of the ACTS project. +// +// Copyright (C) 2016 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 https://mozilla.org/MPL/2.0/. + +#include "ActsExamples/Io/Obj/ObjSimHitWriter.hpp" + +#include "Acts/Definitions/Algebra.hpp" +#include "Acts/Definitions/Common.hpp" +#include "Acts/Geometry/GeometryIdentifier.hpp" +#include "Acts/Visualization/Interpolation3D.hpp" +#include "ActsExamples/EventData/SimHit.hpp" +#include "ActsExamples/Framework/AlgorithmContext.hpp" +#include "ActsExamples/Utilities/Paths.hpp" +#include "ActsFatras/EventData/Barcode.hpp" +#include "ActsFatras/EventData/Hit.hpp" + +#include +#include +#include +#include + +ActsExamples::ObjSimHitWriter::ObjSimHitWriter( + const ActsExamples::ObjSimHitWriter::Config& config, + Acts::Logging::Level level) + : WriterT(config.inputSimHits, "ObjSimHitWriter", level), m_cfg(config) { + // inputSimHits is already checked by base constructor + if (m_cfg.outputStem.empty()) { + throw std::invalid_argument("Missing output filename stem"); + } +} + +ActsExamples::ProcessCode ActsExamples::ObjSimHitWriter::writeT( + const AlgorithmContext& ctx, const ActsExamples::SimHitContainer& simHits) { + // ensure exclusive access to tree/file while writing + std::scoped_lock lock(m_writeMutex); + + // open per-event file for all simhit components + std::string pathSimHit = perEventFilepath( + m_cfg.outputDir, m_cfg.outputStem + ".obj", ctx.eventNumber); + std::string pathSimTrajectory = perEventFilepath( + m_cfg.outputDir, m_cfg.outputStem + "_trajectory.obj", ctx.eventNumber); + + std::ofstream osHits(pathSimHit, std::ofstream::out | std::ofstream::trunc); + std::ofstream osTrajectory(pathSimTrajectory, + std::ofstream::out | std::ofstream::trunc); + + if (!osHits || !osTrajectory) { + throw std::ios_base::failure("Could not open '" + pathSimHit + "' or '" + + pathSimTrajectory + "' to write"); + } + + // Only hit plotting + if (!m_cfg.drawConnections) { + // Write data from internal immplementation + for (const auto& simHit : simHits) { + // local simhit information in global coord. + const Acts::Vector4& globalPos4 = simHit.fourPosition(); + + osHits << "v " << globalPos4[Acts::ePos0] / Acts::UnitConstants::mm << " " + << globalPos4[Acts::ePos1] / Acts::UnitConstants::mm << " " + << globalPos4[Acts::ePos2] / Acts::UnitConstants::mm << std::endl; + + } // end simHit loop + } else { + // We need to associate first + std::unordered_map> particleHits; + // Pre-loop over hits ... write those below threshold + for (const auto& simHit : simHits) { + double momentum = simHit.momentum4Before().head<3>().norm(); + if (momentum < m_cfg.momentumThreshold) { + ACTS_VERBOSE("Skipping: Hit below threshold: " << momentum); + continue; + } else if (momentum < m_cfg.momentumThresholdTraj) { + ACTS_VERBOSE( + "Skipping (trajectory): Hit below threshold: " << momentum); + osHits << "v " + << simHit.fourPosition()[Acts::ePos0] / Acts::UnitConstants::mm + << " " + << simHit.fourPosition()[Acts::ePos1] / Acts::UnitConstants::mm + << " " + << simHit.fourPosition()[Acts::ePos2] / Acts::UnitConstants::mm + << std::endl; + continue; + } + ACTS_VERBOSE("Accepting: Hit above threshold: " << momentum); + + if (particleHits.find(simHit.particleId().value()) == + particleHits.end()) { + particleHits[simHit.particleId().value()] = {}; + } + particleHits[simHit.particleId().value()].push_back( + simHit.fourPosition()); + } + // Draw loop + std::size_t lOffset = 1; + for (auto& [pId, pHits] : particleHits) { + // Draw the particle hits + std::sort(pHits.begin(), pHits.end(), + [](const Acts::Vector4& a, const Acts::Vector4& b) { + return a[Acts::eTime] < b[Acts::eTime]; + }); + + osHits << "o particle_" << pId << std::endl; + for (const auto& hit : pHits) { + osHits << "v " << hit[Acts::ePos0] / Acts::UnitConstants::mm << " " + << hit[Acts::ePos1] / Acts::UnitConstants::mm << " " + << hit[Acts::ePos2] / Acts::UnitConstants::mm << std::endl; + } + osHits << '\n'; + + // Interpolate the points, a minimum number of 3 hits is necessary for + // that + std::vector trajectory; + if (pHits.size() < 3 || m_cfg.nInterpolatedPoints == 0) { + trajectory = pHits; + } else { + // The total number of points is the number of hits times the number of + // interpolated points plus the number of hits + trajectory = Acts::Interpolation3D::spline( + pHits, pHits.size() * (m_cfg.nInterpolatedPoints + 1) - 1, + m_cfg.keepOriginalHits); + } + + osTrajectory << "o particle_trajectory_" << pId << std::endl; + for (const auto& hit : trajectory) { + osTrajectory << "v " << hit[Acts::ePos0] / Acts::UnitConstants::mm + << " " << hit[Acts::ePos1] / Acts::UnitConstants::mm << " " + << hit[Acts::ePos2] / Acts::UnitConstants::mm << std::endl; + } + // Draw the line + for (std::size_t iv = lOffset + 1; iv < lOffset + trajectory.size(); + ++iv) { + osTrajectory << "l " << iv - 1 << " " << iv << '\n'; + } + osTrajectory << '\n'; + // Increase the offset count + lOffset += trajectory.size(); + } + } + osHits.close(); + osTrajectory.close(); + + return ActsExamples::ProcessCode::SUCCESS; +} + +ActsExamples::ProcessCode ActsExamples::ObjSimHitWriter::finalize() { + return ActsExamples::ProcessCode::SUCCESS; +} diff --git a/Examples/Io/Obj/src/ObjTrackingGeometryWriter.cpp b/Examples/Io/Obj/src/ObjTrackingGeometryWriter.cpp index 4937757c85f..d2dbe52e745 100644 --- a/Examples/Io/Obj/src/ObjTrackingGeometryWriter.cpp +++ b/Examples/Io/Obj/src/ObjTrackingGeometryWriter.cpp @@ -6,7 +6,7 @@ // License, v. 2.0. If a copy of the MPL was not distributed with this // file, You can obtain one at https://mozilla.org/MPL/2.0/. -#include "ActsExamples/Plugins/Obj/ObjTrackingGeometryWriter.hpp" +#include "ActsExamples/Io/Obj/ObjTrackingGeometryWriter.hpp" #include "Acts/Utilities/Logger.hpp" #include "ActsExamples/Framework/AlgorithmContext.hpp" diff --git a/Examples/Io/Root/include/ActsExamples/Io/Root/RootMaterialDecorator.hpp b/Examples/Io/Root/include/ActsExamples/Io/Root/RootMaterialDecorator.hpp index f5bab09fb70..0cc10482c9e 100644 --- a/Examples/Io/Root/include/ActsExamples/Io/Root/RootMaterialDecorator.hpp +++ b/Examples/Io/Root/include/ActsExamples/Io/Root/RootMaterialDecorator.hpp @@ -128,7 +128,7 @@ class RootMaterialDecorator : public Acts::IMaterialDecorator { /// Return the maps const Acts::DetectorMaterialMaps materialMaps() const { - return std::make_pair(m_surfaceMaterialMap, m_volumeMaterialMap); + return {m_surfaceMaterialMap, m_volumeMaterialMap}; } /// Get readonly access to the config parameters diff --git a/Examples/Io/Root/include/ActsExamples/Io/Root/RootTrackParameterWriter.hpp b/Examples/Io/Root/include/ActsExamples/Io/Root/RootTrackParameterWriter.hpp index 9f5b2fa2c42..22540f90399 100644 --- a/Examples/Io/Root/include/ActsExamples/Io/Root/RootTrackParameterWriter.hpp +++ b/Examples/Io/Root/include/ActsExamples/Io/Root/RootTrackParameterWriter.hpp @@ -9,7 +9,7 @@ #pragma once #include "Acts/Utilities/Logger.hpp" -#include "ActsExamples/EventData/Index.hpp" +#include "ActsExamples/EventData/Measurement.hpp" #include "ActsExamples/EventData/ProtoTrack.hpp" #include "ActsExamples/EventData/SimHit.hpp" #include "ActsExamples/EventData/SimParticle.hpp" @@ -25,16 +25,14 @@ class TFile; class TTree; namespace ActsExamples { -struct AlgorithmContext; - -using TrackParameterWriter = WriterT; /// Write out the track parameters from both simulation and those estimated from /// reconstructed seeds into a TTree /// /// Each entry in the TTree corresponds to one seed for optimum writing /// speed. The event number is part of the written data. -class RootTrackParameterWriter final : public TrackParameterWriter { +class RootTrackParameterWriter final + : public WriterT { public: struct Config { /// Input estimated track parameters collection. @@ -87,9 +85,9 @@ class RootTrackParameterWriter final : public TrackParameterWriter { "InputProtoTracks"}; ReadDataHandle m_inputParticles{this, "InputParticles"}; ReadDataHandle m_inputSimHits{this, "InputSimHits"}; - ReadDataHandle m_inputMeasurementParticlesMap{ + ReadDataHandle m_inputMeasurementParticlesMap{ this, "InputMeasurementParticlesMap"}; - ReadDataHandle m_inputMeasurementSimHitsMap{ + ReadDataHandle m_inputMeasurementSimHitsMap{ this, "InputMeasurementSimHitsMap"}; /// Mutex used to protect multi-threaded writes diff --git a/Examples/Io/Root/include/ActsExamples/Io/Root/RootTrackStatesWriter.hpp b/Examples/Io/Root/include/ActsExamples/Io/Root/RootTrackStatesWriter.hpp index bacb794113e..73da1c56e48 100644 --- a/Examples/Io/Root/include/ActsExamples/Io/Root/RootTrackStatesWriter.hpp +++ b/Examples/Io/Root/include/ActsExamples/Io/Root/RootTrackStatesWriter.hpp @@ -9,7 +9,7 @@ #pragma once #include "Acts/Utilities/Logger.hpp" -#include "ActsExamples/EventData/Index.hpp" +#include "ActsExamples/EventData/Measurement.hpp" #include "ActsExamples/EventData/SimHit.hpp" #include "ActsExamples/EventData/SimParticle.hpp" #include "ActsExamples/EventData/Track.hpp" @@ -26,12 +26,8 @@ class TFile; class TTree; -namespace ActsFatras { -class Barcode; -} // namespace ActsFatras namespace ActsExamples { -struct AlgorithmContext; /// @class RootTrackStatesWriter /// @@ -108,7 +104,7 @@ class RootTrackStatesWriter final : public WriterT { ReadDataHandle m_inputTrackParticleMatching{ this, "InputTrackParticleMatching"}; ReadDataHandle m_inputSimHits{this, "InputSimHits"}; - ReadDataHandle m_inputMeasurementSimHitsMap{ + ReadDataHandle m_inputMeasurementSimHitsMap{ this, "InputMeasurementSimHitsMap"}; /// Mutex used to protect multi-threaded writes diff --git a/Examples/Io/Root/include/ActsExamples/Io/Root/SeedingPerformanceWriter.hpp b/Examples/Io/Root/include/ActsExamples/Io/Root/SeedingPerformanceWriter.hpp index 49006631b92..b85ab90a0ca 100644 --- a/Examples/Io/Root/include/ActsExamples/Io/Root/SeedingPerformanceWriter.hpp +++ b/Examples/Io/Root/include/ActsExamples/Io/Root/SeedingPerformanceWriter.hpp @@ -9,7 +9,7 @@ #pragma once #include "Acts/Utilities/Logger.hpp" -#include "ActsExamples/EventData/SimHit.hpp" +#include "ActsExamples/EventData/Measurement.hpp" #include "ActsExamples/EventData/SimParticle.hpp" #include "ActsExamples/EventData/SimSeed.hpp" #include "ActsExamples/Framework/DataHandle.hpp" @@ -24,12 +24,8 @@ class TFile; class TTree; -namespace ActsFatras { -class Barcode; -} // namespace ActsFatras namespace ActsExamples { -struct AlgorithmContext; class SeedingPerformanceWriter final : public WriterT { public: @@ -84,7 +80,7 @@ class SeedingPerformanceWriter final : public WriterT { std::size_t m_nTotalDuplicatedParticles = 0; ReadDataHandle m_inputParticles{this, "InputParticles"}; - ReadDataHandle m_inputMeasurementParticlesMap{ + ReadDataHandle m_inputMeasurementParticlesMap{ this, "InputMeasurementParticlesMaps"}; }; diff --git a/Examples/Io/Root/include/ActsExamples/Io/Root/TrackFinderNTupleWriter.hpp b/Examples/Io/Root/include/ActsExamples/Io/Root/TrackFinderNTupleWriter.hpp index 799471be784..fc090773b98 100644 --- a/Examples/Io/Root/include/ActsExamples/Io/Root/TrackFinderNTupleWriter.hpp +++ b/Examples/Io/Root/include/ActsExamples/Io/Root/TrackFinderNTupleWriter.hpp @@ -18,7 +18,6 @@ #include namespace ActsExamples { -struct AlgorithmContext; /// Write track finder performance measures. /// @@ -31,8 +30,8 @@ class TrackFinderNTupleWriter final : public WriterT { std::string inputTracks; /// Input particles collection. std::string inputParticles; - /// Input hit-particles map collection. - std::string inputMeasurementParticlesMap; + /// Input particle-measurements map collection. + std::string inputParticleMeasurementsMap; /// Input proto track-particle matching. std::string inputTrackParticleMatching; /// Output filename. diff --git a/Examples/Io/Root/src/RootNuclearInteractionParametersWriter.cpp b/Examples/Io/Root/src/RootNuclearInteractionParametersWriter.cpp index 995f7bdf57d..265fde166b9 100644 --- a/Examples/Io/Root/src/RootNuclearInteractionParametersWriter.cpp +++ b/Examples/Io/Root/src/RootNuclearInteractionParametersWriter.cpp @@ -121,7 +121,7 @@ buildNotNormalisedMap(TH1F const* hist) { if (integral == 0.) { histoBorders.clear(); temp_HistoContents.clear(); - return std::make_tuple(histoBorders, temp_HistoContents, integral); + return {histoBorders, temp_HistoContents, integral}; } // Set the bin borders @@ -130,7 +130,7 @@ buildNotNormalisedMap(TH1F const* hist) { } histoBorders[nBins] = hist->GetXaxis()->GetXmax(); - return std::make_tuple(histoBorders, temp_HistoContents, integral); + return {histoBorders, temp_HistoContents, integral}; } /// @brief This function combines neighbouring bins with the same value @@ -169,7 +169,7 @@ std::pair, std::vector> buildMap( // Fast exit if the histogram is empty if (histoContents.empty()) { - return std::make_pair(std::get<0>(map), std::vector()); + return {std::get<0>(map), std::vector()}; } // Set the bin content @@ -183,7 +183,7 @@ std::pair, std::vector> buildMap( auto histoBorders = std::get<0>(map); reduceMap(histoBorders, normalisedHistoContents); - return std::make_pair(histoBorders, normalisedHistoContents); + return {histoBorders, normalisedHistoContents}; } /// @brief This method transforms a probability distribution into components @@ -208,7 +208,7 @@ std::pair, std::vector> buildMap( // Fast exit if the histogram is empty if (histoContents.empty()) { - return std::make_pair(std::get<0>(map), std::vector()); + return {std::get<0>(map), std::vector()}; } // Set the bin content @@ -223,7 +223,7 @@ std::pair, std::vector> buildMap( std::vector histoBorders = std::get<0>(map); reduceMap(histoBorders, normalisedHistoContents); - return std::make_pair(histoBorders, normalisedHistoContents); + return {histoBorders, normalisedHistoContents}; } /// @brief This method builds decomposed cumulative probability distributions diff --git a/Examples/Io/Root/src/RootTrackParameterWriter.cpp b/Examples/Io/Root/src/RootTrackParameterWriter.cpp index a4e74a44e2c..520a875f1a8 100644 --- a/Examples/Io/Root/src/RootTrackParameterWriter.cpp +++ b/Examples/Io/Root/src/RootTrackParameterWriter.cpp @@ -18,16 +18,13 @@ #include "ActsExamples/Validation/TrackClassification.hpp" #include "ActsFatras/EventData/Barcode.hpp" #include "ActsFatras/EventData/Hit.hpp" -#include "ActsFatras/EventData/Particle.hpp" #include #include #include #include -#include #include #include -#include #include #include @@ -43,8 +40,8 @@ namespace ActsExamples { RootTrackParameterWriter::RootTrackParameterWriter( const RootTrackParameterWriter::Config& config, Acts::Logging::Level level) - : TrackParameterWriter(config.inputTrackParameters, - "RootTrackParameterWriter", level), + : WriterT(config.inputTrackParameters, + "RootTrackParameterWriter", level), m_cfg(config) { if (m_cfg.inputProtoTracks.empty()) { throw std::invalid_argument("Missing proto tracks input collection"); diff --git a/Examples/Io/Root/src/RootTrackStatesWriter.cpp b/Examples/Io/Root/src/RootTrackStatesWriter.cpp index 37d7cdf2fc9..ee199e48ff3 100644 --- a/Examples/Io/Root/src/RootTrackStatesWriter.cpp +++ b/Examples/Io/Root/src/RootTrackStatesWriter.cpp @@ -471,13 +471,13 @@ ProcessCode RootTrackStatesWriter::writeT(const AlgorithmContext& ctx, auto getTrackParams = [&](unsigned int ipar) -> std::optional> { if (ipar == ePredicted && state.hasPredicted()) { - return std::make_pair(state.predicted(), state.predictedCovariance()); + return std::pair(state.predicted(), state.predictedCovariance()); } if (ipar == eFiltered && state.hasFiltered()) { - return std::make_pair(state.filtered(), state.filteredCovariance()); + return std::pair(state.filtered(), state.filteredCovariance()); } if (ipar == eSmoothed && state.hasSmoothed()) { - return std::make_pair(state.smoothed(), state.smoothedCovariance()); + return std::pair(state.smoothed(), state.smoothedCovariance()); } if (ipar == eUnbiased && state.hasSmoothed() && state.hasProjector() && state.hasCalibrated()) { diff --git a/Examples/Io/Root/src/TrackFinderNTupleWriter.cpp b/Examples/Io/Root/src/TrackFinderNTupleWriter.cpp index 3ddfc2b6f84..4b9aa94b4d0 100644 --- a/Examples/Io/Root/src/TrackFinderNTupleWriter.cpp +++ b/Examples/Io/Root/src/TrackFinderNTupleWriter.cpp @@ -9,8 +9,7 @@ #include "ActsExamples/Io/Root/TrackFinderNTupleWriter.hpp" #include "Acts/Definitions/Units.hpp" -#include "ActsExamples/EventData/Index.hpp" -#include "ActsExamples/EventData/SimHit.hpp" +#include "ActsExamples/EventData/Measurement.hpp" #include "ActsExamples/EventData/SimParticle.hpp" #include "ActsExamples/EventData/TruthMatching.hpp" #include "ActsExamples/Framework/AlgorithmContext.hpp" @@ -18,7 +17,6 @@ #include "ActsExamples/Utilities/Range.hpp" #include "ActsExamples/Validation/TrackClassification.hpp" #include "ActsFatras/EventData/Barcode.hpp" -#include "ActsFatras/EventData/Particle.hpp" #include #include @@ -32,11 +30,13 @@ #include #include -struct ActsExamples::TrackFinderNTupleWriter::Impl { +namespace ActsExamples { + +struct TrackFinderNTupleWriter::Impl { Config cfg; ReadDataHandle inputParticles; - ReadDataHandle inputMeasurementParticlesMap; + ReadDataHandle inputParticleMeasurementsMap; ReadDataHandle inputTrackParticleMatching; TFile* file = nullptr; @@ -78,7 +78,7 @@ struct ActsExamples::TrackFinderNTupleWriter::Impl { // particle charge in e float prtQ = 0; // particle reconstruction - UShort_t prtNumHits = 0; // number of hits for this particle + UShort_t prtNumMeasurements = 0; // number of hits for this particle UShort_t prtNumTracks = 0; // number of tracks this particle was reconstructed in UShort_t prtNumTracksMajority = @@ -89,7 +89,7 @@ struct ActsExamples::TrackFinderNTupleWriter::Impl { Impl(TrackFinderNTupleWriter* parent, Config&& c, const Acts::Logger& l) : cfg(std::move(c)), inputParticles{parent, "InputParticles"}, - inputMeasurementParticlesMap{parent, "InputMeasurementParticlesMap"}, + inputParticleMeasurementsMap{parent, "InputParticleMeasurementsMap"}, inputTrackParticleMatching{parent, "InputTrackParticleMatching"}, _logger(l) { if (cfg.inputTracks.empty()) { @@ -98,8 +98,9 @@ struct ActsExamples::TrackFinderNTupleWriter::Impl { if (cfg.inputParticles.empty()) { throw std::invalid_argument("Missing particles input collection"); } - if (cfg.inputMeasurementParticlesMap.empty()) { - throw std::invalid_argument("Missing hit-particles map input collection"); + if (cfg.inputParticleMeasurementsMap.empty()) { + throw std::invalid_argument( + "Missing particle-measurements map input collection"); } if (cfg.inputTrackParticleMatching.empty()) { throw std::invalid_argument( @@ -110,7 +111,7 @@ struct ActsExamples::TrackFinderNTupleWriter::Impl { } inputParticles.initialize(cfg.inputParticles); - inputMeasurementParticlesMap.initialize(cfg.inputMeasurementParticlesMap); + inputParticleMeasurementsMap.initialize(cfg.inputParticleMeasurementsMap); inputTrackParticleMatching.initialize(cfg.inputTrackParticleMatching); // the output file can not be given externally since TFile accesses to the @@ -146,7 +147,7 @@ struct ActsExamples::TrackFinderNTupleWriter::Impl { prtTree->Branch("pz", &prtPz); prtTree->Branch("m", &prtM); prtTree->Branch("q", &prtQ); - prtTree->Branch("nhits", &prtNumHits); + prtTree->Branch("nhits", &prtNumMeasurements); prtTree->Branch("ntracks", &prtNumTracks); prtTree->Branch("ntracks_majority", &prtNumTracksMajority); } @@ -155,10 +156,8 @@ struct ActsExamples::TrackFinderNTupleWriter::Impl { void write(std::uint64_t eventId, const ConstTrackContainer& tracks, const SimParticleContainer& particles, - const HitParticlesMap& hitParticlesMap, + const ParticleMeasurementsMap& particleMeasurementsMap, const TrackParticleMatching& trackParticleMatching) { - const auto& particleHitsMap = invertIndexMultimap(hitParticlesMap); - // How often a particle was reconstructed. std::unordered_map reconCount; reconCount.reserve(particles.size()); @@ -209,8 +208,8 @@ struct ActsExamples::TrackFinderNTupleWriter::Impl { for (const auto& phc : particleMatch.contributingParticles) { trkParticleId.push_back(phc.particleId.value()); // count total number of hits for this particle - auto trueParticleHits = - makeRange(particleHitsMap.equal_range(phc.particleId.value())); + auto trueParticleHits = makeRange( + particleMeasurementsMap.equal_range(phc.particleId.value())); trkParticleNumHitsTotal.push_back(trueParticleHits.size()); trkParticleNumHitsOnTrack.push_back(phc.hitCount); } @@ -223,9 +222,9 @@ struct ActsExamples::TrackFinderNTupleWriter::Impl { { std::lock_guard guardPrt(trkMutex); for (const auto& particle : particles) { - // find all hits for this particle - auto hits = - makeRange(particleHitsMap.equal_range(particle.particleId())); + // find all measurements for this particle + auto measurements = makeRange( + particleMeasurementsMap.equal_range(particle.particleId())); // identification prtEventId = eventId; @@ -243,7 +242,7 @@ struct ActsExamples::TrackFinderNTupleWriter::Impl { prtM = particle.mass() / Acts::UnitConstants::GeV; prtQ = particle.charge() / Acts::UnitConstants::e; // reconstruction - prtNumHits = hits.size(); + prtNumMeasurements = measurements.size(); auto nt = reconCount.find(particle.particleId()); prtNumTracks = (nt != reconCount.end()) ? nt->second : 0u; auto nm = majorityCount.find(particle.particleId()); @@ -265,31 +264,31 @@ struct ActsExamples::TrackFinderNTupleWriter::Impl { } }; -ActsExamples::TrackFinderNTupleWriter::TrackFinderNTupleWriter( - ActsExamples::TrackFinderNTupleWriter::Config config, - Acts::Logging::Level level) +TrackFinderNTupleWriter::TrackFinderNTupleWriter( + TrackFinderNTupleWriter::Config config, Acts::Logging::Level level) : WriterT(config.inputTracks, "TrackFinderNTupleWriter", level), m_impl(std::make_unique(this, std::move(config), logger())) {} -ActsExamples::TrackFinderNTupleWriter::~TrackFinderNTupleWriter() = default; +TrackFinderNTupleWriter::~TrackFinderNTupleWriter() = default; -ActsExamples::ProcessCode ActsExamples::TrackFinderNTupleWriter::writeT( - const ActsExamples::AlgorithmContext& ctx, - const ActsExamples::ConstTrackContainer& tracks) { +ProcessCode TrackFinderNTupleWriter::writeT(const AlgorithmContext& ctx, + const ConstTrackContainer& tracks) { const auto& particles = m_impl->inputParticles(ctx); - const auto& hitParticlesMap = m_impl->inputMeasurementParticlesMap(ctx); + const auto& particleMeasurementsMap = + m_impl->inputParticleMeasurementsMap(ctx); const auto& trackParticleMatching = m_impl->inputTrackParticleMatching(ctx); - m_impl->write(ctx.eventNumber, tracks, particles, hitParticlesMap, + m_impl->write(ctx.eventNumber, tracks, particles, particleMeasurementsMap, trackParticleMatching); return ProcessCode::SUCCESS; } -ActsExamples::ProcessCode ActsExamples::TrackFinderNTupleWriter::finalize() { +ProcessCode TrackFinderNTupleWriter::finalize() { m_impl->close(); return ProcessCode::SUCCESS; } -const ActsExamples::TrackFinderNTupleWriter::Config& -ActsExamples::TrackFinderNTupleWriter::config() const { +const TrackFinderNTupleWriter::Config& TrackFinderNTupleWriter::config() const { return m_impl->cfg; } + +} // namespace ActsExamples diff --git a/Examples/Io/Root/src/detail/NuclearInteractionParametrisation.cpp b/Examples/Io/Root/src/detail/NuclearInteractionParametrisation.cpp index 244c30ca039..ee89b828590 100644 --- a/Examples/Io/Root/src/detail/NuclearInteractionParametrisation.cpp +++ b/Examples/Io/Root/src/detail/NuclearInteractionParametrisation.cpp @@ -87,7 +87,7 @@ std::pair calculateMeanAndCovariance( } covariance /= events.size(); - return std::make_pair(mean, covariance); + return {mean, covariance}; } EigenspaceComponents calculateEigenspace(const Vector& mean, @@ -99,7 +99,7 @@ EigenspaceComponents calculateEigenspace(const Vector& mean, // Transform the mean vector into eigenspace Vector meanEigenspace = eigenvectors * mean; - return std::make_tuple(eigenvalues, eigenvectors, meanEigenspace); + return {eigenvalues, eigenvectors, meanEigenspace}; } Parametrisation buildMomentumParameters(const EventCollection& events, @@ -122,7 +122,7 @@ Parametrisation buildMomentumParameters(const EventCollection& events, EigenspaceComponents eigenspaceElements = calculateEigenspace(meanAndCovariance.first, meanAndCovariance.second); // Calculate the cumulative distributions - return std::make_pair(eigenspaceElements, histos); + return {eigenspaceElements, histos}; } EventProperties prepareMomenta(const EventCollection& events, @@ -255,7 +255,7 @@ Parametrisation buildInvariantMassParameters(const EventCollection& events, EigenspaceComponents eigenspaceElements = calculateEigenspace(meanAndCovariance.first, meanAndCovariance.second); // Calculate the cumulative distributions - return std::make_pair(eigenspaceElements, histos); + return {eigenspaceElements, histos}; } std::unordered_map> @@ -341,7 +341,7 @@ cumulativeMultiplicityProbability(const EventCollection& events, } } - return std::make_pair(softHisto, hardHisto); + return {softHisto, hardHisto}; } TVectorF softProbability(const EventCollection& events) { diff --git a/Examples/Python/python/acts/examples/reconstruction.py b/Examples/Python/python/acts/examples/reconstruction.py index a365572cabc..727bf9e9785 100644 --- a/Examples/Python/python/acts/examples/reconstruction.py +++ b/Examples/Python/python/acts/examples/reconstruction.py @@ -234,12 +234,6 @@ def trackSelectorDefaultKWArgs(c): defaults=[None] * 3, ) -AmbiguityResolutionMLDBScanConfig = namedtuple( - "AmbiguityResolutionMLDBScanConfig", - ["nMeasurementsMin", "epsilonDBScan", "minPointsDBScan"], - defaults=[None] * 3, -) - SeedFilterMLDBScanConfig = namedtuple( "SeedFilterMLDBScanConfig", ["epsilonDBScan", "minPointsDBScan", "minSeedScore"], @@ -580,7 +574,7 @@ def addTruthSmearedSeeding( truthTrkFndAlg = acts.examples.TruthTrackFinder( level=logLevel, inputParticles=selectedParticles, - inputMeasurementParticlesMap="measurement_particles_map", + inputParticleMeasurementsMap="particle_measurements_map", outputProtoTracks="truth_particle_tracks", ) s.addAlgorithm(truthTrkFndAlg) @@ -601,7 +595,7 @@ def addTruthEstimatedSeeding( truthSeeding = acts.examples.TruthSeedingAlgorithm( level=logLevel, inputParticles=inputParticles, - inputMeasurementParticlesMap="measurement_particles_map", + inputParticleMeasurementsMap="particle_measurements_map", inputSpacePoints=[spacePoints], outputParticles="truth_seeded_particles", outputProtoTracks="truth_particle_tracks", @@ -1857,7 +1851,7 @@ def addExaTrkX( inputProtoTracks=findingAlg.config.outputProtoTracks, # the original selected particles after digitization inputParticles="particles_initial", - inputMeasurementParticlesMap="measurement_particles_map", + inputParticleMeasurementsMap="particle_measurements_map", inputTrackParticleMatching=matchAlg.config.outputTrackParticleMatching, filePath=str(Path(outputDirRoot) / "performance_track_finding.root"), ) @@ -1972,14 +1966,31 @@ def addScoreBasedAmbiguityResolution( s.addAlgorithm(algScoreBased) s.addWhiteboardAlias("tracks", algScoreBased.config.outputTracks) + matchAlg = acts.examples.TrackTruthMatcher( + level=customLogLevel(), + inputTracks=algScoreBased.config.outputTracks, + inputParticles="particles", + inputMeasurementParticlesMap="measurement_particles_map", + outputTrackParticleMatching="ambi_scorebased_track_particle_matching", + outputParticleTrackMatching="ambi_scorebased_particle_track_matching", + doubleMatching=True, + ) + s.addAlgorithm(matchAlg) + s.addWhiteboardAlias( + "track_particle_matching", matchAlg.config.outputTrackParticleMatching + ) + s.addWhiteboardAlias( + "particle_track_matching", matchAlg.config.outputParticleTrackMatching + ) + addTrackWriters( s, name="ambi_scorebased", tracks=algScoreBased.config.outputTracks, outputDirCsv=outputDirCsv, outputDirRoot=outputDirRoot, - writeSummary=writeTrackStates, - writeStates=writeTrackSummary, + writeSummary=writeTrackSummary, + writeStates=writeTrackStates, writeFitterPerformance=writePerformance, writeFinderPerformance=writePerformance, writeCovMat=writeCovMat, @@ -1995,6 +2006,7 @@ def addScoreBasedAmbiguityResolution( def addAmbiguityResolutionML( s, config: AmbiguityResolutionMLConfig = AmbiguityResolutionMLConfig(), + tracks: str = "tracks", onnxModelFile: Optional[Union[Path, str]] = None, outputDirCsv: Optional[Union[Path, str]] = None, outputDirRoot: Optional[Union[Path, str]] = None, @@ -2008,10 +2020,9 @@ def addAmbiguityResolutionML( from acts.examples import GreedyAmbiguityResolutionAlgorithm customLogLevel = acts.examples.defaultLogging(s, logLevel) - algML = AmbiguityResolutionMLAlgorithm( level=customLogLevel(), - inputTracks="tracks", + inputTracks=tracks, inputDuplicateNN=onnxModelFile, outputTracks="ambiTracksML", **acts.examples.defaultKWArgs( @@ -2033,63 +2044,33 @@ def addAmbiguityResolutionML( s.addAlgorithm(algML) s.addAlgorithm(algGreedy) - addTrackWriters( - s, - name="ambiML", - tracks=algGreedy.config.outputTracks, - outputDirCsv=outputDirCsv, - outputDirRoot=outputDirRoot, - writeSummary=writeTrackStates, - writeStates=writeTrackSummary, - writeFitterPerformance=writePerformance, - writeFinderPerformance=writePerformance, - writeCovMat=writeCovMat, - logLevel=logLevel, - ) - - return s - + s.addWhiteboardAlias("tracks", algGreedy.config.outputTracks) -@acts.examples.NamedTypeArgs( - config=AmbiguityResolutionMLDBScanConfig, -) -def addAmbiguityResolutionMLDBScan( - s, - config: AmbiguityResolutionMLDBScanConfig = AmbiguityResolutionMLDBScanConfig(), - onnxModelFile: Optional[Union[Path, str]] = None, - outputDirCsv: Optional[Union[Path, str]] = None, - outputDirRoot: Optional[Union[Path, str]] = None, - writeTrackSummary: bool = True, - writeTrackStates: bool = False, - writePerformance: bool = True, - writeCovMat=False, - logLevel: Optional[acts.logging.Level] = None, -) -> None: - from acts.examples import AmbiguityResolutionMLDBScanAlgorithm - - customLogLevel = acts.examples.defaultLogging(s, logLevel) - - alg = AmbiguityResolutionMLDBScanAlgorithm( + matchAlg = acts.examples.TrackTruthMatcher( level=customLogLevel(), - inputTracks="tracks", - inputDuplicateNN=onnxModelFile, - outputTracks="ambiTracksMLDBScan", - **acts.examples.defaultKWArgs( - nMeasurementsMin=config.nMeasurementsMin, - epsilonDBScan=config.epsilonDBScan, - minPointsDBScan=config.minPointsDBScan, - ), + inputTracks=algGreedy.config.outputTracks, + inputParticles="particles", + inputMeasurementParticlesMap="measurement_particles_map", + outputTrackParticleMatching="ambiML_track_particle_matching", + outputParticleTrackMatching="ambiML_particle_track_matching", + doubleMatching=True, + ) + s.addAlgorithm(matchAlg) + s.addWhiteboardAlias( + "track_particle_matching", matchAlg.config.outputTrackParticleMatching + ) + s.addWhiteboardAlias( + "particle_track_matching", matchAlg.config.outputParticleTrackMatching ) - s.addAlgorithm(alg) addTrackWriters( s, - name="ambiMLDBScan", - trajectories=alg.config.outputTracks, - outputDirRoot=outputDirRoot, + name="ambiML", + tracks=algGreedy.config.outputTracks, outputDirCsv=outputDirCsv, - writeSummary=writeTrackStates, - writeStates=writeTrackSummary, + outputDirRoot=outputDirRoot, + writeSummary=writeTrackSummary, + writeStates=writeTrackStates, writeFitterPerformance=writePerformance, writeFinderPerformance=writePerformance, writeCovMat=writeCovMat, diff --git a/Examples/Python/python/acts/examples/simulation.py b/Examples/Python/python/acts/examples/simulation.py index b60b05466a0..411f4410361 100644 --- a/Examples/Python/python/acts/examples/simulation.py +++ b/Examples/Python/python/acts/examples/simulation.py @@ -420,6 +420,7 @@ def addFatras( outputSimHits: str = "simhits", outputDirCsv: Optional[Union[Path, str]] = None, outputDirRoot: Optional[Union[Path, str]] = None, + outputDirObj: Optional[Union[Path, str]] = None, logLevel: Optional[acts.logging.Level] = None, ) -> None: """This function steers the detector simulation using Fatras @@ -444,6 +445,8 @@ def addFatras( the output folder for the Csv output, None triggers no output outputDirRoot : Path|str, path, None the output folder for the Root output, None triggers no output + outputDirObj : Path|str, path, None + the output folder for the Obj output, None triggers no output """ customLogLevel = acts.examples.defaultLogging(s, logLevel) @@ -507,6 +510,7 @@ def addFatras( particlesPostSelected, outputDirCsv, outputDirRoot, + outputDirObj, logLevel, ) @@ -519,6 +523,7 @@ def addSimWriters( particlesSimulated: str = "particles_simulated", outputDirCsv: Optional[Union[Path, str]] = None, outputDirRoot: Optional[Union[Path, str]] = None, + outputDirObj: Optional[Union[Path, str]] = None, logLevel: Optional[acts.logging.Level] = None, ) -> None: customLogLevel = acts.examples.defaultLogging(s, logLevel) @@ -563,6 +568,19 @@ def addSimWriters( ) ) + if outputDirObj is not None: + outputDirObj = Path(outputDirObj) + if not outputDirObj.exists(): + outputDirObj.mkdir() + s.addWriter( + acts.examples.ObjSimHitWriter( + level=customLogLevel(), + inputSimHits=simHits, + outputDir=str(outputDirObj), + outputStem="hits", + ) + ) + def getG4DetectorConstructionFactory( detector: Any, @@ -620,6 +638,7 @@ def addGeant4( keepParticlesWithoutHits=True, outputDirCsv: Optional[Union[Path, str]] = None, outputDirRoot: Optional[Union[Path, str]] = None, + outputDirObj: Optional[Union[Path, str]] = None, logLevel: Optional[acts.logging.Level] = None, killVolume: Optional[acts.Volume] = None, killAfterTime: float = float("inf"), @@ -647,6 +666,8 @@ def addGeant4( the output folder for the Csv output, None triggers no output outputDirRoot : Path|str, path, None the output folder for the Root output, None triggers no output + outputDirObj : Path|str, path, None + the output folder for the Obj output, None triggers no output killVolume: acts.Volume, None if given, particles are killed when going outside this volume. killAfterTime: float @@ -737,7 +758,8 @@ def addGeant4( particlesPostSelected, outputDirCsv, outputDirRoot, - logLevel, + outputDirObj, + logLevel=logLevel, ) return s @@ -789,6 +811,8 @@ def addDigitization( outputMeasurements="measurements", outputMeasurementParticlesMap="measurement_particles_map", outputMeasurementSimHitsMap="measurement_simhits_map", + outputParticleMeasurementsMap="particle_measurements_map", + outputSimHitMeasurementsMap="simhit_measurements_map", **acts.examples.defaultKWArgs( doMerge=doMerge, ), diff --git a/Examples/Python/src/Digitization.cpp b/Examples/Python/src/Digitization.cpp index 5bcce4b74ca..cfafb7bf627 100644 --- a/Examples/Python/src/Digitization.cpp +++ b/Examples/Python/src/Digitization.cpp @@ -55,6 +55,8 @@ void addDigitization(Context& ctx) { ACTS_PYTHON_MEMBER(outputClusters); ACTS_PYTHON_MEMBER(outputMeasurementParticlesMap); ACTS_PYTHON_MEMBER(outputMeasurementSimHitsMap); + ACTS_PYTHON_MEMBER(outputParticleMeasurementsMap); + ACTS_PYTHON_MEMBER(outputSimHitMeasurementsMap); ACTS_PYTHON_MEMBER(surfaceByIdentifier); ACTS_PYTHON_MEMBER(randomNumbers); ACTS_PYTHON_MEMBER(doOutputCells); diff --git a/Examples/Python/src/Geometry.cpp b/Examples/Python/src/Geometry.cpp index d69b35a86ba..4c15e0dcf94 100644 --- a/Examples/Python/src/Geometry.cpp +++ b/Examples/Python/src/Geometry.cpp @@ -492,6 +492,21 @@ void addExperimentalGeometry(Context& ctx) { .def(py::init, const Extent&>()); } + { + using RangeXDDim3 = Acts::RangeXD<3u, double>; + + py::class_(m, "RangeXDDim3") + .def(py::init([](const std::array& range0, + const std::array& range1, + const std::array& range2) { + RangeXDDim3 range; + range[0].shrink(range0[0], range0[1]); + range[1].shrink(range1[0], range1[1]); + range[2].shrink(range2[0], range2[1]); + return range; + })); + } + { // The external volume structure builder py::class_>(m, "NullBField") .def(py::init<>()); + py::class_>(m, "MultiRangeBField") + .def(py::init< + std::vector, Acts::Vector3>>>()); + { using Config = Acts::SolenoidBField::Config; diff --git a/Examples/Python/src/Onnx.cpp b/Examples/Python/src/Onnx.cpp index c22e595bb43..8bd75a3047a 100644 --- a/Examples/Python/src/Onnx.cpp +++ b/Examples/Python/src/Onnx.cpp @@ -8,7 +8,6 @@ #include "Acts/Plugins/Python/Utilities.hpp" #include "ActsExamples/TrackFindingML/AmbiguityResolutionMLAlgorithm.hpp" -#include "ActsExamples/TrackFindingML/AmbiguityResolutionMLDBScanAlgorithm.hpp" #include "ActsExamples/TrackFindingML/SeedFilterMLAlgorithm.hpp" #include @@ -31,11 +30,6 @@ void addOnnx(Context& ctx) { inputTracks, inputDuplicateNN, outputTracks, nMeasurementsMin); - ACTS_PYTHON_DECLARE_ALGORITHM( - ActsExamples::AmbiguityResolutionMLDBScanAlgorithm, onnx, - "AmbiguityResolutionMLDBScanAlgorithm", inputTracks, inputDuplicateNN, - outputTracks, nMeasurementsMin, epsilonDBScan, minPointsDBScan); - ACTS_PYTHON_DECLARE_ALGORITHM(ActsExamples::SeedFilterMLAlgorithm, onnx, "SeedFilterMLAlgorithm", inputTrackParameters, inputSimSeeds, inputSeedFilterNN, diff --git a/Examples/Python/src/Output.cpp b/Examples/Python/src/Output.cpp index 4281f2db26e..838649816f6 100644 --- a/Examples/Python/src/Output.cpp +++ b/Examples/Python/src/Output.cpp @@ -6,13 +6,10 @@ // License, v. 2.0. If a copy of the MPL was not distributed with this // file, You can obtain one at https://mozilla.org/MPL/2.0/. -#include "Acts/Definitions/TrackParametrization.hpp" -#include "Acts/Geometry/GeometryHierarchyMap.hpp" #include "Acts/Plugins/Python/Utilities.hpp" #include "Acts/Utilities/Logger.hpp" #include "Acts/Visualization/IVisualization3D.hpp" #include "Acts/Visualization/ViewConfig.hpp" -#include "ActsExamples/Digitization/DigitizationAlgorithm.hpp" #include "ActsExamples/Io/Csv/CsvBFieldWriter.hpp" #include "ActsExamples/Io/Csv/CsvExaTrkXGraphWriter.hpp" #include "ActsExamples/Io/Csv/CsvMeasurementWriter.hpp" @@ -20,11 +17,14 @@ #include "ActsExamples/Io/Csv/CsvProtoTrackWriter.hpp" #include "ActsExamples/Io/Csv/CsvSeedWriter.hpp" #include "ActsExamples/Io/Csv/CsvSimHitWriter.hpp" +#include "ActsExamples/Io/Csv/CsvSpacePointWriter.hpp" #include "ActsExamples/Io/Csv/CsvSpacePointsBucketWriter.hpp" -#include "ActsExamples/Io/Csv/CsvSpacepointWriter.hpp" #include "ActsExamples/Io/Csv/CsvTrackParameterWriter.hpp" #include "ActsExamples/Io/Csv/CsvTrackWriter.hpp" #include "ActsExamples/Io/Csv/CsvTrackingGeometryWriter.hpp" +#include "ActsExamples/Io/Obj/ObjPropagationStepsWriter.hpp" +#include "ActsExamples/Io/Obj/ObjSimHitWriter.hpp" +#include "ActsExamples/Io/Obj/ObjTrackingGeometryWriter.hpp" #include "ActsExamples/Io/Root/RootBFieldWriter.hpp" #include "ActsExamples/Io/Root/RootMaterialTrackWriter.hpp" #include "ActsExamples/Io/Root/RootMaterialWriter.hpp" @@ -46,8 +46,6 @@ #include "ActsExamples/Io/Root/TrackFitterPerformanceWriter.hpp" #include "ActsExamples/Io/Root/VertexNTupleWriter.hpp" #include "ActsExamples/MaterialMapping/IMaterialWriter.hpp" -#include "ActsExamples/Plugins/Obj/ObjPropagationStepsWriter.hpp" -#include "ActsExamples/Plugins/Obj/ObjTrackingGeometryWriter.hpp" #include "ActsExamples/TrackFinding/ITrackParamsLookupReader.hpp" #include "ActsExamples/TrackFinding/ITrackParamsLookupWriter.hpp" @@ -111,6 +109,12 @@ void addOutput(Context& ctx) { "ObjPropagationStepsWriter", collection, outputDir, outputScalor, outputPrecision); + ACTS_PYTHON_DECLARE_WRITER(ActsExamples::ObjSimHitWriter, mex, + "ObjSimHitWriter", inputSimHits, outputDir, + outputStem, outputPrecision, drawConnections, + momentumThreshold, momentumThresholdTraj, + nInterpolatedPoints, keepOriginalHits); + { auto c = py::class_(m, "ViewConfig").def(py::init<>()); @@ -201,7 +205,7 @@ void addOutput(Context& ctx) { ACTS_PYTHON_DECLARE_WRITER(ActsExamples::TrackFinderNTupleWriter, mex, "TrackFinderNTupleWriter", inputTracks, - inputParticles, inputMeasurementParticlesMap, + inputParticles, inputParticleMeasurementsMap, inputTrackParticleMatching, filePath, fileMode, treeNameTracks, treeNameParticles); @@ -371,8 +375,8 @@ void addOutput(Context& ctx) { "CsvSimHitWriter", inputSimHits, outputDir, outputStem, outputPrecision); - ACTS_PYTHON_DECLARE_WRITER(ActsExamples::CsvSpacepointWriter, mex, - "CsvSpacepointWriter", inputSpacepoints, outputDir, + ACTS_PYTHON_DECLARE_WRITER(ActsExamples::CsvSpacePointWriter, mex, + "CsvSpacePointWriter", inputSpacepoints, outputDir, outputPrecision); ACTS_PYTHON_DECLARE_WRITER(ActsExamples::CsvSpacePointsBucketWriter, mex, diff --git a/Examples/Python/src/TruthTracking.cpp b/Examples/Python/src/TruthTracking.cpp index 319ee24e2ee..4799d93dbdb 100644 --- a/Examples/Python/src/TruthTracking.cpp +++ b/Examples/Python/src/TruthTracking.cpp @@ -40,7 +40,7 @@ void addTruthTracking(Context& ctx) { ACTS_PYTHON_DECLARE_ALGORITHM( ActsExamples::TruthTrackFinder, mex, "TruthTrackFinder", inputParticles, - inputMeasurementParticlesMap, outputProtoTracks); + inputParticleMeasurementsMap, outputProtoTracks); ACTS_PYTHON_DECLARE_ALGORITHM(ActsExamples::ParticleTrackParamExtractor, mex, "ParticleTrackParamExtractor", inputParticles, @@ -159,7 +159,7 @@ void addTruthTracking(Context& ctx) { ACTS_PYTHON_DECLARE_ALGORITHM( ActsExamples::TruthSeedingAlgorithm, mex, "TruthSeedingAlgorithm", - inputParticles, inputMeasurementParticlesMap, inputSpacePoints, + inputParticles, inputParticleMeasurementsMap, inputSpacePoints, outputParticles, outputSeeds, outputProtoTracks, deltaRMin, deltaRMax); ACTS_PYTHON_DECLARE_ALGORITHM(ActsExamples::HitSelector, mex, "HitSelector", diff --git a/Examples/Python/tests/test_magnetic_field.py b/Examples/Python/tests/test_magnetic_field.py index a51ebd55ab6..d1336e3ca8a 100644 --- a/Examples/Python/tests/test_magnetic_field.py +++ b/Examples/Python/tests/test_magnetic_field.py @@ -72,3 +72,41 @@ def test_solenoid(conf_const): ) assert isinstance(field, acts.examples.InterpolatedMagneticField2) + + +def test_multiregion_bfield(): + with pytest.raises(TypeError): + acts.MultiRangeBField() + + rs = [ + (acts.RangeXDDim3((0, 3), (0, 3), (0, 3)), acts.Vector3(0.0, 0.0, 2.0)), + (acts.RangeXDDim3((1, 2), (1, 2), (1, 10)), acts.Vector3(2.0, 0.0, 0.0)), + ] + f = acts.MultiRangeBField(rs) + assert f + + ctx = acts.MagneticFieldContext() + assert ctx + + fc = f.makeCache(ctx) + assert fc + + rv = f.getField(acts.Vector3(0.5, 0.5, 0.5), fc) + assert rv[0] == pytest.approx(0.0) + assert rv[1] == pytest.approx(0.0) + assert rv[2] == pytest.approx(2.0) + + rv = f.getField(acts.Vector3(1.5, 1.5, 5.0), fc) + assert rv[0] == pytest.approx(2.0) + assert rv[1] == pytest.approx(0.0) + assert rv[2] == pytest.approx(0.0) + + rv = f.getField(acts.Vector3(2.5, 2.5, 2.5), fc) + assert rv[0] == pytest.approx(0.0) + assert rv[1] == pytest.approx(0.0) + assert rv[2] == pytest.approx(2.0) + + rv = f.getField(acts.Vector3(1.5, 1.5, 1.5), fc) + assert rv[0] == pytest.approx(2.0) + assert rv[1] == pytest.approx(0.0) + assert rv[2] == pytest.approx(0.0) diff --git a/Examples/Scripts/Python/MLAmbiguityResolution/ambiguity_solver_network.py b/Examples/Scripts/Python/MLAmbiguityResolution/ambiguity_solver_network.py index 15f596414f7..5e241fdf717 100644 --- a/Examples/Scripts/Python/MLAmbiguityResolution/ambiguity_solver_network.py +++ b/Examples/Scripts/Python/MLAmbiguityResolution/ambiguity_solver_network.py @@ -16,7 +16,7 @@ def prepareDataSet(data: pd.DataFrame) -> pd.DataFrame: data = data # Remove tracks with less than 7 measurements data = data[data["nMeasurements"] > 6] - # data = data.sort_values("good/duplicate/fake", ascending=False) + data = data.sort_values("good/duplicate/fake", ascending=False) # Remove pure duplicate (tracks purely identical) keep the ones good one if among them. data = data.drop_duplicates( subset=[ diff --git a/Examples/Scripts/Python/full_chain_odd.py b/Examples/Scripts/Python/full_chain_odd.py index 8b86402b820..7ddaa02054f 100755 --- a/Examples/Scripts/Python/full_chain_odd.py +++ b/Examples/Scripts/Python/full_chain_odd.py @@ -136,6 +136,12 @@ default=True, action=argparse.BooleanOptionalAction, ) +parser.add_argument( + "--output-obj", + help="Switch obj output on/off", + default=True, + action=argparse.BooleanOptionalAction, +) args = parser.parse_args() @@ -277,6 +283,7 @@ ), outputDirRoot=outputDir if args.output_root else None, outputDirCsv=outputDir if args.output_csv else None, + outputDirObj=outputDir if args.output_obj else None, rnd=rnd, killVolume=trackingGeometry.highestTrackingVolume, killAfterTime=25 * u.ns, @@ -305,6 +312,7 @@ enableInteractions=True, outputDirRoot=outputDir if args.output_root else None, outputDirCsv=outputDir if args.output_csv else None, + outputDirObj=outputDir if args.output_obj else None, rnd=rnd, ) diff --git a/Examples/Scripts/Python/geant4.py b/Examples/Scripts/Python/geant4.py index 42cb5db95bb..65e154c2704 100755 --- a/Examples/Scripts/Python/geant4.py +++ b/Examples/Scripts/Python/geant4.py @@ -36,6 +36,7 @@ def runGeant4( field, outputDirCsv=outputDir / "csv", outputDirRoot=outputDir, + outputDirObj=outputDir / "obj", rnd=rnd, materialMappings=materialMappings, volumeMappings=volumeMappings, diff --git a/Examples/Scripts/compareRootFiles.hpp b/Examples/Scripts/compareRootFiles.hpp index 464230373be..a7043b15e96 100644 --- a/Examples/Scripts/compareRootFiles.hpp +++ b/Examples/Scripts/compareRootFiles.hpp @@ -34,8 +34,7 @@ class AnyVector { static std::pair*> create(Args&&... args) { std::vector* vector = new std::vector(std::forward(args)...); std::function deleter = [vector] { delete vector; }; - return std::make_pair( - AnyVector{static_cast(vector), std::move(deleter)}, vector); + return {AnyVector{static_cast(vector), std::move(deleter)}, vector}; } // Default-construct a null type-erased vector diff --git a/Fatras/include/ActsFatras/Digitization/UncorrelatedHitSmearer.hpp b/Fatras/include/ActsFatras/Digitization/UncorrelatedHitSmearer.hpp index caf5efa8c61..a4825b2a99b 100644 --- a/Fatras/include/ActsFatras/Digitization/UncorrelatedHitSmearer.hpp +++ b/Fatras/include/ActsFatras/Digitization/UncorrelatedHitSmearer.hpp @@ -86,7 +86,7 @@ struct BoundParametersSmearer { ParametersVector par = ParametersVector::Zero(); CovarianceMatrix cov = CovarianceMatrix::Zero(); - for (int i = 0; i < static_cast(kSize); ++i) { + for (std::size_t i = 0; i < kSize; ++i) { auto res = smearFunctions[i](boundParams[indices[i]], rng); if (!res.ok()) { return Result::failure(res.error()); diff --git a/Fatras/include/ActsFatras/Physics/ElectroMagnetic/PhotonConversion.hpp b/Fatras/include/ActsFatras/Physics/ElectroMagnetic/PhotonConversion.hpp index af45481258e..582cc073bc4 100644 --- a/Fatras/include/ActsFatras/Physics/ElectroMagnetic/PhotonConversion.hpp +++ b/Fatras/include/ActsFatras/Physics/ElectroMagnetic/PhotonConversion.hpp @@ -129,8 +129,8 @@ std::pair PhotonConversion::generatePathLimits( // Fast exit if not a photon or the energy is too low if (particle.pdg() != Acts::PdgParticle::eGamma || particle.absoluteMomentum() < (2 * kElectronMass)) { - return std::make_pair(std::numeric_limits::infinity(), - std::numeric_limits::infinity()); + return {std::numeric_limits::infinity(), + std::numeric_limits::infinity()}; } // Use for the moment only Al data - Yung Tsai - Rev.Mod.Particle Physics Vol. @@ -155,11 +155,11 @@ std::pair PhotonConversion::generatePathLimits( std::uniform_real_distribution uniformDistribution{0., 1.}; // This is a transformation of eq. 3.75 - return std::make_pair(-9. / 7. * - std::log(conversionProbScaleFactor * - (1 - uniformDistribution(generator))) / - (1. - xi), - std::numeric_limits::infinity()); + return {-9. / 7. * + std::log(conversionProbScaleFactor * + (1 - uniformDistribution(generator))) / + (1. - xi), + std::numeric_limits::infinity()}; } template diff --git a/Fatras/include/ActsFatras/Physics/NuclearInteraction/NuclearInteraction.hpp b/Fatras/include/ActsFatras/Physics/NuclearInteraction/NuclearInteraction.hpp index d372cf8ffc9..d0a0d84c36c 100644 --- a/Fatras/include/ActsFatras/Physics/NuclearInteraction/NuclearInteraction.hpp +++ b/Fatras/include/ActsFatras/Physics/NuclearInteraction/NuclearInteraction.hpp @@ -58,8 +58,8 @@ struct NuclearInteraction { const Particle& particle) const { // Fast exit: No parameterisation provided if (multiParticleParameterisation.empty()) { - return std::make_pair(std::numeric_limits::infinity(), - std::numeric_limits::infinity()); + return {std::numeric_limits::infinity(), + std::numeric_limits::infinity()}; } // Find the parametrisation that corresponds to the particle type for (const auto& particleParametrisation : multiParticleParameterisation) { @@ -77,15 +77,13 @@ struct NuclearInteraction { // Set the L0 limit if not done already const auto& distribution = parametrisation.nuclearInteractionProbability; - auto limits = - std::make_pair(std::numeric_limits::infinity(), - sampleContinuousValues( - uniformDistribution(generator), distribution)); - return limits; + return {std::numeric_limits::infinity(), + sampleContinuousValues(uniformDistribution(generator), + distribution)}; } } - return std::make_pair(std::numeric_limits::infinity(), - std::numeric_limits::infinity()); + return {std::numeric_limits::infinity(), + std::numeric_limits::infinity()}; } /// This method performs a nuclear interaction. @@ -483,14 +481,14 @@ NuclearInteraction::sampleKinematics( if (trials == nMatchingTrialsTotal) { return std::nullopt; } - // Re-sampole invariant masses if no fitting momenta were found + // Re-sample invariant masses if no fitting momenta were found if (trials++ % nMatchingTrials == 0) { invariantMasses = sampleInvariantMasses(generator, parameters); } else { momenta = sampleMomenta(generator, parameters, momentum); } } - return std::make_pair(momenta, invariantMasses); + return std::pair(momenta, invariantMasses); } template diff --git a/Fatras/src/Physics/NuclearInteraction/NuclearInteraction.cpp b/Fatras/src/Physics/NuclearInteraction/NuclearInteraction.cpp index 84ab93921fe..ad1e0bbb6a6 100644 --- a/Fatras/src/Physics/NuclearInteraction/NuclearInteraction.cpp +++ b/Fatras/src/Physics/NuclearInteraction/NuclearInteraction.cpp @@ -141,7 +141,7 @@ std::pair NuclearInteraction::globalAngle(double phi1, const float theta = std::acos(vectorSum.z() / vectorSum.norm()); const float phi = std::atan2(vectorSum.y(), vectorSum.x()); - return std::make_pair(phi, theta); + return {phi, theta}; } bool NuclearInteraction::match(const Acts::ActsDynamicVector& momenta, diff --git a/Plugins/DD4hep/src/DD4hepBlueprintFactory.cpp b/Plugins/DD4hep/src/DD4hepBlueprintFactory.cpp index bf653fff801..06ff97153dd 100644 --- a/Plugins/DD4hep/src/DD4hepBlueprintFactory.cpp +++ b/Plugins/DD4hep/src/DD4hepBlueprintFactory.cpp @@ -207,7 +207,7 @@ Acts::Experimental::DD4hepBlueprintFactory::extractExternals( aux += "vol. binning : " + binningString; } // Return the tuple - return std::make_tuple(transform, bValueType, bValues, bBinning, aux); + return {transform, bValueType, bValues, bBinning, aux}; } std::tuple, @@ -309,6 +309,5 @@ Acts::Experimental::DD4hepBlueprintFactory::extractInternals( std::make_shared(geoIdCfg); } - return std::make_tuple(internalsBuilder, rootsFinderBuilder, geoIdGenerator, - aux, ext); + return {internalsBuilder, rootsFinderBuilder, geoIdGenerator, aux, ext}; } diff --git a/Plugins/Geant4/src/Geant4Converters.cpp b/Plugins/Geant4/src/Geant4Converters.cpp index 1608f45bb10..6db950f570a 100644 --- a/Plugins/Geant4/src/Geant4Converters.cpp +++ b/Plugins/Geant4/src/Geant4Converters.cpp @@ -110,7 +110,7 @@ Acts::Geant4ShapeConverter::cylinderBounds(const G4Tubs& g4Tubs) { } double thickness = g4Tubs.GetOuterRadius() - g4Tubs.GetInnerRadius(); auto cBounds = std::make_shared(tArray); - return std::make_tuple(std::move(cBounds), thickness); + return {std::move(cBounds), thickness}; } std::tuple, double> @@ -133,7 +133,7 @@ Acts::Geant4ShapeConverter::radialBounds(const G4Tubs& g4Tubs) { } double thickness = g4Tubs.GetZHalfLength() * 2; auto rBounds = std::make_shared(tArray); - return std::make_tuple(std::move(rBounds), thickness); + return {std::move(rBounds), thickness}; } std::shared_ptr Acts::Geant4ShapeConverter::lineBounds( @@ -173,7 +173,7 @@ Acts::Geant4ShapeConverter::rectangleBounds(const G4Box& g4Box) { } auto rBounds = std::make_shared(hG4XYZ[std::abs(rAxes[0u])], hG4XYZ[std::abs(rAxes[1u])]); - return std::make_tuple(std::move(rBounds), rAxes, thickness); + return {std::move(rBounds), rAxes, thickness}; } std::tuple, std::array, double> @@ -226,7 +226,7 @@ Acts::Geant4ShapeConverter::trapezoidBounds(const G4Trd& g4Trd) { auto tBounds = std::make_shared( halfLengthXminY, halfLengthXmaxY, halfLengthY); - return std::make_tuple(std::move(tBounds), rAxes, thickness); + return {std::move(tBounds), rAxes, thickness}; } std::tuple, std::array, double> @@ -234,19 +234,19 @@ Acts::Geant4ShapeConverter::planarBounds(const G4VSolid& g4Solid) { const G4Box* box = dynamic_cast(&g4Solid); if (box != nullptr) { auto [rBounds, axes, thickness] = rectangleBounds(*box); - return std::make_tuple(std::move(rBounds), axes, thickness); + return {std::move(rBounds), axes, thickness}; } const G4Trd* trd = dynamic_cast(&g4Solid); if (trd != nullptr) { auto [tBounds, axes, thickness] = trapezoidBounds(*trd); - return std::make_tuple(std::move(tBounds), axes, thickness); + return {std::move(tBounds), axes, thickness}; } std::shared_ptr pBounds = nullptr; std::array rAxes = {}; double rThickness = 0.; - return std::make_tuple(std::move(pBounds), rAxes, rThickness); + return {std::move(pBounds), rAxes, rThickness}; } namespace { diff --git a/Plugins/GeoModel/src/GeoModelBlueprintCreater.cpp b/Plugins/GeoModel/src/GeoModelBlueprintCreater.cpp index 429321828d6..081834d7069 100644 --- a/Plugins/GeoModel/src/GeoModelBlueprintCreater.cpp +++ b/Plugins/GeoModel/src/GeoModelBlueprintCreater.cpp @@ -320,7 +320,7 @@ Acts::GeoModelBlueprintCreater::createInternalStructureBuilder( const std::vector& internalConstraints) const { // Check if the internals entry is empty if (entry.internals.empty()) { - return std::make_tuple(nullptr, Extent()); + return {nullptr, Extent()}; } // Build a layer structure @@ -393,10 +393,10 @@ Acts::GeoModelBlueprintCreater::createInternalStructureBuilder( lsbCfg.nMinimalSurfaces = surfaces.size() + 1u; } - return std::make_tuple( + return { std::make_shared( lsbCfg, m_logger->clone(entry.name + "_LayerStructureBuilder")), - internalExtent); + internalExtent}; } else { throw std::invalid_argument( @@ -404,7 +404,7 @@ Acts::GeoModelBlueprintCreater::createInternalStructureBuilder( entry.internals[1u] + "' / or now kdt surfaces provided."); } } - return std::make_tuple(nullptr, Extent()); + return {nullptr, Extent()}; } std::tuple, @@ -434,5 +434,5 @@ Acts::GeoModelBlueprintCreater::parseBounds( "supported for the moment."); } - return std::make_tuple(boundsType, extent, boundValues, translation); + return {boundsType, extent, boundValues, translation}; } diff --git a/Plugins/GeoModel/src/detail/GeoModelExtentHelper.cpp b/Plugins/GeoModel/src/detail/GeoModelExtentHelper.cpp index daab4fd10e6..3734aefb8fe 100644 --- a/Plugins/GeoModel/src/detail/GeoModelExtentHelper.cpp +++ b/Plugins/GeoModel/src/detail/GeoModelExtentHelper.cpp @@ -111,7 +111,7 @@ Acts::detail::GeoModelExentHelper::extentFromTable( BinningValue bValue = bvCyl.at(iv); double val = std::numeric_limits::max(); bool isMin = (iv % 2 == 0); - // Case "e" : exxternal extent + // Case "e" : external extent if (value == "e") { // External parameters do not constrain it if (!externalExtent.constrains(bValue)) { @@ -157,5 +157,5 @@ Acts::detail::GeoModelExentHelper::extentFromTable( } } - return std::make_tuple(boundsType, extent); + return {boundsType, extent}; } diff --git a/Plugins/Onnx/include/Acts/Plugins/Onnx/AmbiguityTrackClassifier.hpp b/Plugins/Onnx/include/Acts/Plugins/Onnx/AmbiguityTrackClassifier.hpp index 3e1450f2b85..b0ddbd09bc7 100644 --- a/Plugins/Onnx/include/Acts/Plugins/Onnx/AmbiguityTrackClassifier.hpp +++ b/Plugins/Onnx/include/Acts/Plugins/Onnx/AmbiguityTrackClassifier.hpp @@ -102,23 +102,6 @@ class AmbiguityTrackClassifier { return goodTracks; } - /// Select the track associated with each cluster - /// - /// @param clusters is a map of clusters, each cluster correspond to a vector of track ID - /// @param tracks is the input track container - /// @return a vector of trackID corresponding tho the good tracks - template - std::vector solveAmbiguity( - std::unordered_map>& clusters, - const track_container_t& tracks) const { - std::vector> outputTensor = - inferScores(clusters, tracks); - std::vector goodTracks = - trackSelection(clusters, outputTensor); - - return goodTracks; - } - private: // ONNX environment Ort::Env m_env; diff --git a/Tests/UnitTests/Core/Geometry/LayerCreatorTests.cpp b/Tests/UnitTests/Core/Geometry/LayerCreatorTests.cpp index 13aec3ff545..20f43a3b68e 100644 --- a/Tests/UnitTests/Core/Geometry/LayerCreatorTests.cpp +++ b/Tests/UnitTests/Core/Geometry/LayerCreatorTests.cpp @@ -230,7 +230,7 @@ struct LayerCreatorFixture { } } - return std::make_pair(res, pairs); + return {res, pairs}; } }; diff --git a/Tests/UnitTests/Core/Geometry/SurfaceArrayCreatorTests.cpp b/Tests/UnitTests/Core/Geometry/SurfaceArrayCreatorTests.cpp index 00b6509ef25..6d69c1667da 100644 --- a/Tests/UnitTests/Core/Geometry/SurfaceArrayCreatorTests.cpp +++ b/Tests/UnitTests/Core/Geometry/SurfaceArrayCreatorTests.cpp @@ -228,7 +228,7 @@ struct SurfaceArrayCreatorFixture { } } - return std::make_pair(res, pairs); + return {res, pairs}; } }; diff --git a/Tests/UnitTests/Core/MagneticField/CMakeLists.txt b/Tests/UnitTests/Core/MagneticField/CMakeLists.txt index 60a829bf548..4600cc8352e 100644 --- a/Tests/UnitTests/Core/MagneticField/CMakeLists.txt +++ b/Tests/UnitTests/Core/MagneticField/CMakeLists.txt @@ -1,4 +1,5 @@ add_unittest(ConstantBField ConstantBFieldTests.cpp) add_unittest(InterpolatedBFieldMap InterpolatedBFieldMapTests.cpp) add_unittest(SolenoidBField SolenoidBFieldTests.cpp) +add_unittest(MultiRangeBField MultiRangeBFieldTests.cpp) add_unittest(MagneticFieldProvider MagneticFieldProviderTests.cpp) diff --git a/Tests/UnitTests/Core/MagneticField/MultiRangeBFieldTests.cpp b/Tests/UnitTests/Core/MagneticField/MultiRangeBFieldTests.cpp new file mode 100644 index 00000000000..3644db058b4 --- /dev/null +++ b/Tests/UnitTests/Core/MagneticField/MultiRangeBFieldTests.cpp @@ -0,0 +1,74 @@ +// This file is part of the ACTS project. +// +// Copyright (C) 2016 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 https://mozilla.org/MPL/2.0/. + +#include +#include + +#include "Acts/MagneticField/MagneticFieldContext.hpp" +#include "Acts/MagneticField/MultiRangeBField.hpp" +#include "Acts/Utilities/Result.hpp" + +namespace Acts::Test { + +MagneticFieldContext mfContext = MagneticFieldContext(); + +BOOST_AUTO_TEST_CASE(TestMultiRangeBField) { + std::vector, Vector3>> inputs; + + inputs.emplace_back(RangeXD<3, double>{{0., 0., 0.}, {3., 3., 3.}}, + Vector3{0., 0., 2.}); + inputs.emplace_back(RangeXD<3, double>{{1., 1., 1.}, {2., 2., 10.}}, + Vector3{2., 0., 0.}); + + const MultiRangeBField bfield(std::move(inputs)); + + auto bcache = bfield.makeCache(mfContext); + + // Test a point inside the first volume. + { + Result r = bfield.getField({0.5, 0.5, 0.5}, bcache); + BOOST_CHECK(r.ok()); + BOOST_CHECK_CLOSE((*r)[0], 0., 0.05); + BOOST_CHECK_CLOSE((*r)[1], 0., 0.05); + BOOST_CHECK_CLOSE((*r)[2], 2., 0.05); + } + + // Test a point inside the second volume, not overlapping the first. + { + Result r = bfield.getField({1.5, 1.5, 5.}, bcache); + BOOST_CHECK(r.ok()); + BOOST_CHECK_CLOSE((*r)[0], 2., 0.05); + BOOST_CHECK_CLOSE((*r)[1], 0., 0.05); + BOOST_CHECK_CLOSE((*r)[2], 0., 0.05); + } + + // Test a point inside the first volume again. + { + Result r = bfield.getField({2.5, 2.5, 2.5}, bcache); + BOOST_CHECK(r.ok()); + BOOST_CHECK_CLOSE((*r)[0], 0., 0.05); + BOOST_CHECK_CLOSE((*r)[1], 0., 0.05); + BOOST_CHECK_CLOSE((*r)[2], 2., 0.05); + } + + // Test a point inside the second volume in the overlap region. + { + Result r = bfield.getField({1.5, 1.5, 1.5}, bcache); + BOOST_CHECK(r.ok()); + BOOST_CHECK_CLOSE((*r)[0], 2., 0.05); + BOOST_CHECK_CLOSE((*r)[1], 0., 0.05); + BOOST_CHECK_CLOSE((*r)[2], 0., 0.05); + } + + // Test a point outside all volumes. + { + Result r = bfield.getField({-1., -1., -1}, bcache); + BOOST_CHECK(!r.ok()); + } +} +} // namespace Acts::Test diff --git a/Tests/UnitTests/Core/SpacePointFormation/SpacePointBuilderTests.cpp b/Tests/UnitTests/Core/SpacePointFormation/SpacePointBuilderTests.cpp index c20239bb078..6d043178320 100644 --- a/Tests/UnitTests/Core/SpacePointFormation/SpacePointBuilderTests.cpp +++ b/Tests/UnitTests/Core/SpacePointFormation/SpacePointBuilderTests.cpp @@ -93,7 +93,7 @@ std::pair stripEnds( auto gPos1 = surface->localToGlobal(gctx, lpos1, globalFakeMom); auto gPos2 = surface->localToGlobal(gctx, lpos2, globalFakeMom); - return std::make_pair(gPos1, gPos2); + return {gPos1, gPos2}; } // Create a test context diff --git a/Tests/UnitTests/Core/Utilities/AlgebraHelpersTests.cpp b/Tests/UnitTests/Core/Utilities/AlgebraHelpersTests.cpp index 86165ddb2ba..87e47d3e3a9 100644 --- a/Tests/UnitTests/Core/Utilities/AlgebraHelpersTests.cpp +++ b/Tests/UnitTests/Core/Utilities/AlgebraHelpersTests.cpp @@ -51,6 +51,15 @@ BOOST_AUTO_TEST_CASE(safeInverseLargeMatrix) { BOOST_CHECK_EQUAL(*identityInv, identity); } +BOOST_AUTO_TEST_CASE(safeInverseDynamicMatrix) { + Eigen::MatrixXd identity{Eigen::MatrixXd::Identity(2, 2)}; + + auto identityInv = Acts::safeInverse(identity); + + BOOST_CHECK(identityInv); + BOOST_CHECK_EQUAL(*identityInv, identity); +} + BOOST_AUTO_TEST_CASE(SafeInverseBadSmallMatrix) { Eigen::Matrix m; m << 1, 1, 2, 2; @@ -111,13 +120,6 @@ BOOST_AUTO_TEST_CASE(SafeInverseFPELargeMatrix) { // auto mInv = Acts::safeInverse(m); // } -/// This test should not compile -// BOOST_AUTO_TEST_CASE(SafeInverseDynamicMatrix) { -// Eigen::MatrixXd m{Eigen::MatrixXd::Identity(2, 2)}; -// -// auto mInv = Acts::safeInverse(m); -// } - BOOST_AUTO_TEST_SUITE_END() BOOST_AUTO_TEST_SUITE(SafeExp) diff --git a/Tests/UnitTests/Core/Vertexing/VertexingDataHelper.hpp b/Tests/UnitTests/Core/Vertexing/VertexingDataHelper.hpp index d984b4c8c44..383b5b84d2c 100644 --- a/Tests/UnitTests/Core/Vertexing/VertexingDataHelper.hpp +++ b/Tests/UnitTests/Core/Vertexing/VertexingDataHelper.hpp @@ -143,7 +143,7 @@ readTracksAndVertexCSV(const std::string& toolString, vertices.push_back(vertexInfo); } - return std::make_tuple(beamspotConstraint, vertices, tracks); + return {beamspotConstraint, vertices, tracks}; } } // namespace Acts::Test diff --git a/Tests/UnitTests/Core/Visualization/CMakeLists.txt b/Tests/UnitTests/Core/Visualization/CMakeLists.txt index ca6768e876b..97d76217917 100644 --- a/Tests/UnitTests/Core/Visualization/CMakeLists.txt +++ b/Tests/UnitTests/Core/Visualization/CMakeLists.txt @@ -1,5 +1,6 @@ add_unittest(Visualization3D Visualization3DTests.cpp) add_unittest(EventDataView3D EventDataView3DTests.cpp) +add_unittest(Interpolation3D Interpolation3DTests.cpp) add_unittest(PrimitivesView3D PrimitivesView3DTests.cpp) add_unittest(SurfaceView3D SurfaceView3DTests.cpp) add_unittest(TrackingGeometryView3D TrackingGeometryView3DTests.cpp) diff --git a/Tests/UnitTests/Core/Visualization/Interpolation3DTests.cpp b/Tests/UnitTests/Core/Visualization/Interpolation3DTests.cpp new file mode 100644 index 00000000000..ba19c55b03f --- /dev/null +++ b/Tests/UnitTests/Core/Visualization/Interpolation3DTests.cpp @@ -0,0 +1,92 @@ +// This file is part of the ACTS project. +// +// Copyright (C) 2016 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 https://mozilla.org/MPL/2.0/. + +#include + +#include "Acts/Tests/CommonHelpers/FloatComparisons.hpp" +#include "Acts/Visualization/Interpolation3D.hpp" + +#include + +namespace Acts::Test { + +BOOST_AUTO_TEST_SUITE(Visualization) + +BOOST_AUTO_TEST_CASE(SplineInterpolationEigen) { + /// Define the input vector + double R = 10.; + std::vector inputs; + + // Interpolate the points options + std::vector trajectory; + + // Empty in empty out check + trajectory = Acts::Interpolation3D::spline(inputs, 10); + BOOST_CHECK(trajectory.empty()); + + for (double phi = 0; phi < 2 * std::numbers::pi; + phi += std::numbers::pi / 4) { + inputs.push_back(Acts::Vector3(R * cos(phi), R * sin(phi), 0.)); + } + + // (0) - nothing happens + trajectory = Acts::Interpolation3D::spline(inputs, 1); + // Check input and output size are the same + BOOST_CHECK_EQUAL(trajectory.size(), inputs.size()); + + // (1) - interpolate between the points with 12 points in total + trajectory = Acts::Interpolation3D::spline(inputs, 12); + // Check the output size is correct + BOOST_CHECK_EQUAL(trajectory.size(), 12); + + for (const auto& point : trajectory) { + // Check the interpolated points are on the circle + // with a tolerance of course + CHECK_CLOSE_ABS(point.norm(), R, 0.1); + // Verify points remain in the XY plane + CHECK_CLOSE_ABS(point.z(), 0., 0.1); + } +} + +BOOST_AUTO_TEST_CASE(SplineInterpolationArray) { + /// Define the input vector + std::vector> inputs; + + for (double x = 0; x < 10; x += 1) { + inputs.push_back({x, x * x, 0.}); + } + + // This time we keep the original hits + auto trajectory = Acts::Interpolation3D::spline(inputs, 100, true); + + // Check the outpu type is correct + constexpr bool isOutput = + std::is_same_v; + BOOST_CHECK(isOutput); + + // Check the output size is correct + BOOST_CHECK_EQUAL(trajectory.size(), 108); +} + +BOOST_AUTO_TEST_CASE(SplineInterpolationErrors) { + std::vector> inputs; + + // Test with single point + inputs.push_back({0., 0., 0.}); + auto result = Acts::Interpolation3D::spline(inputs, 10); + BOOST_CHECK_EQUAL(result.size(), 1); + + // Test with two points + inputs.push_back({1., 1., 1.}); + result = Acts::Interpolation3D::spline(inputs, 10); + BOOST_CHECK_EQUAL(result.size(), 2); +} + +BOOST_AUTO_TEST_SUITE_END() + +} // namespace Acts::Test diff --git a/Tests/UnitTests/Plugins/Cuda/Seeding/SeedFinderCudaTest.cpp b/Tests/UnitTests/Plugins/Cuda/Seeding/SeedFinderCudaTest.cpp index 1b31f083fd8..87d61b1df71 100644 --- a/Tests/UnitTests/Plugins/Cuda/Seeding/SeedFinderCudaTest.cpp +++ b/Tests/UnitTests/Plugins/Cuda/Seeding/SeedFinderCudaTest.cpp @@ -231,7 +231,7 @@ int main(int argc, char** argv) { -> std::tuple> { Acts::Vector3 position(sp.x(), sp.y(), sp.z()); Acts::Vector2 variance(sp.varianceR, sp.varianceZ); - return std::make_tuple(position, variance, std::nullopt); + return {position, variance, std::nullopt}; }; // setup spacepoint grid config diff --git a/Tests/UnitTests/Plugins/Cuda/Seeding2/main.cpp b/Tests/UnitTests/Plugins/Cuda/Seeding2/main.cpp index b63d549e37f..3b134aca9a0 100644 --- a/Tests/UnitTests/Plugins/Cuda/Seeding2/main.cpp +++ b/Tests/UnitTests/Plugins/Cuda/Seeding2/main.cpp @@ -113,7 +113,7 @@ int main(int argc, char* argv[]) { -> std::tuple> { Acts::Vector3 position(sp.x(), sp.y(), sp.z()); Acts::Vector2 covariance(sp.m_varianceR, sp.m_varianceZ); - return std::make_tuple(position, covariance, std::nullopt); + return {position, covariance, std::nullopt}; }; // extent used to store r range for middle spacepoint diff --git a/Tests/UnitTests/Plugins/Geant4/Geant4SurfaceProviderTests.cpp b/Tests/UnitTests/Plugins/Geant4/Geant4SurfaceProviderTests.cpp index 49521ae786d..7508fe9a148 100644 --- a/Tests/UnitTests/Plugins/Geant4/Geant4SurfaceProviderTests.cpp +++ b/Tests/UnitTests/Plugins/Geant4/Geant4SurfaceProviderTests.cpp @@ -123,7 +123,7 @@ ConstructGeant4World() { parser.SetOutputFileOverwrite(true); parser.Write(gdmlPath.string(), physWorld); - return std::make_tuple(physWorld, names); + return {physWorld, names}; } auto gctx = Acts::GeometryContext(); diff --git a/docs/core/magnetic_field.md b/docs/core/magnetic_field.md index 1924234bde0..488d5890203 100644 --- a/docs/core/magnetic_field.md +++ b/docs/core/magnetic_field.md @@ -235,6 +235,28 @@ analytical implementation and is much faster to lookup: :::{doxygenfunction} Acts::solenoidFieldMap ::: +### Multi-range constant field + +The multi-range constant field allows modelling cases where a magnetic field +can be described as multiple (potentially overlapping) regions, each of which +has its own constant magnetic field. This provides more flexibility than the +{class}`Acts::ConstantBField` while providing higher performance than +{class}`Acts::InterpolatedBFieldMap`. + +This magnetic field provider is configured using a list of pairs, where each +pair defines a region in three-dimensional space as well as a field vector. +Magnetic field lookup then proceeds by finding the _last_ region in the +user-provided list that contains the requested coordinate and returning the +corresponding field vector. + +The implementation uses a simple caching mechanism to store the last matched +region, providing improved performance for consecutive lookups within the same +region. This is thread-safe when each thread uses its own cache instance. The +field configuration itself is immutable after construction. + +:::{doxygenclass} Acts::MultiRangeBField +::: + ## Full provider interface :::{doxygenclass} Acts::MagneticFieldProvider diff --git a/docs/tracking.md b/docs/tracking.md index 54e582d1e0f..442a3793adc 100644 --- a/docs/tracking.md +++ b/docs/tracking.md @@ -97,16 +97,17 @@ position $l$ and the momentum vector $\vec p$ are shown. Aside from the nominal quantities captured in $\vec x$, the related uncertainties and correlations need to be taken into account as well. They -can be expressed as a $5\times 5$ covariance matrix like +can be expressed as a $6\times 6$ covariance matrix like \begin{equation*} C = \begin{bmatrix} - \sigma^2(l_0)& \text{cov}(l_0,l_1) & \text{cov}(l_0, \phi) & \text{cov}(l_0, \theta) & \text{cov}(l_0, q/p) \\ - . & \sigma^2(l_1) & \text{cov}(l_1, \phi) & \text{cov}(l_1, \theta) & \text{cov}(l_1, q/p) \\ - . & . & \sigma^2(\phi) & \text{cov}(\phi,\theta) & \text{cov}(\phi, q/p) \\ - . & . & . & \sigma^2(\theta) & \text{cov}(\theta, q/p) \\ - . & . & . & . & \sigma^2(q/p) + \sigma^2(l_0)& \text{cov}(l_0,l_1) & \text{cov}(l_0, \phi) & \text{cov}(l_0, \theta) & \text{cov}(l_0, q/p) & \text{cov}(l_0, t) \\ + . & \sigma^2(l_1) & \text{cov}(l_1, \phi) & \text{cov}(l_1, \theta) & \text{cov}(l_1, q/p) & \text{cov}(l_1, t) \\ + . & . & \sigma^2(\phi) & \text{cov}(\phi,\theta) & \text{cov}(\phi, q/p) & \text{cov}(\phi, t) \\ + . & . & . & \sigma^2(\theta) & \text{cov}(\theta, q/p) & \text{cov}(\theta, t) \\ + . & . & . & . & \sigma^2(q/p) & \text{cov}(q/p, t) \\ + . & . & . & . & . & \sigma^2(t) \end{bmatrix} \end{equation*}