From 60ff9805719ea948b9ff8ef35b1870d39eae2864 Mon Sep 17 00:00:00 2001 From: Andreas Salzburger Date: Wed, 22 Nov 2023 10:17:20 +0100 Subject: [PATCH 1/2] chore: update to newer podio versions (#2606) This PR has a small modification to access the index of `podio::ObjectID` directly, and thus makes ACTS compatible with newer `podio` versions. Co-authored-by: Paul Gessinger <1058585+paulgessinger@users.noreply.github.com> --- .../ActsExamples/Io/EDM4hep/EDM4hepUtil.hpp | 14 ++++++++++++++ Examples/Io/EDM4hep/src/EDM4hepParticleReader.cpp | 3 +-- Examples/Io/EDM4hep/src/EDM4hepSimHitReader.cpp | 7 +++---- Examples/Io/EDM4hep/src/EDM4hepUtil.cpp | 3 ++- 4 files changed, 20 insertions(+), 7 deletions(-) diff --git a/Examples/Io/EDM4hep/include/ActsExamples/Io/EDM4hep/EDM4hepUtil.hpp b/Examples/Io/EDM4hep/include/ActsExamples/Io/EDM4hep/EDM4hepUtil.hpp index 4574b0c0407..11b731125bc 100644 --- a/Examples/Io/EDM4hep/include/ActsExamples/Io/EDM4hep/EDM4hepUtil.hpp +++ b/Examples/Io/EDM4hep/include/ActsExamples/Io/EDM4hep/EDM4hepUtil.hpp @@ -118,5 +118,19 @@ void writeTrajectory(const Acts::GeometryContext& gctx, double Bz, const Acts::ParticleHypothesis& particleHypothesis, const IndexMultimap& hitParticlesMap); +/// Helper function to either return an id as is, or unpack an index from it +/// if it is a podio::ObjectID. +/// @tparam T The type of the id. +/// @param o The id to convert. +/// @return The id as an unsigned integer. +template +uint64_t podioObjectIDToInteger(T&& o) { + if constexpr (!std::is_same_v) { + return o; + } else { + return o.index; + } +} + } // namespace EDM4hepUtil } // namespace ActsExamples diff --git a/Examples/Io/EDM4hep/src/EDM4hepParticleReader.cpp b/Examples/Io/EDM4hep/src/EDM4hepParticleReader.cpp index 0cbf3ee9a1a..96c6d5572b2 100644 --- a/Examples/Io/EDM4hep/src/EDM4hepParticleReader.cpp +++ b/Examples/Io/EDM4hep/src/EDM4hepParticleReader.cpp @@ -55,8 +55,7 @@ ProcessCode EDM4hepParticleReader::read(const AlgorithmContext& ctx) { auto particle = EDM4hepUtil::readParticle(mcParticle, [](const edm4hep::MCParticle& p) { ActsFatras::Barcode result; - // TODO dont use podio internal id - result.setParticle(p.id()); + result.setParticle(EDM4hepUtil::podioObjectIDToInteger(p.id())); return result; }); unordered.push_back(particle); diff --git a/Examples/Io/EDM4hep/src/EDM4hepSimHitReader.cpp b/Examples/Io/EDM4hep/src/EDM4hepSimHitReader.cpp index 00b5de279d9..aa036811f3f 100644 --- a/Examples/Io/EDM4hep/src/EDM4hepSimHitReader.cpp +++ b/Examples/Io/EDM4hep/src/EDM4hepSimHitReader.cpp @@ -54,8 +54,7 @@ ProcessCode EDM4hepSimHitReader::read(const AlgorithmContext& ctx) { auto particle = EDM4hepUtil::readParticle( mcParticle, [](const edm4hep::MCParticle& p) { ActsFatras::Barcode result; - // TODO dont use podio internal id - result.setParticle(p.id()); + result.setParticle(EDM4hepUtil::podioObjectIDToInteger(p.id())); return result; }); unordered.push_back(particle); @@ -78,8 +77,8 @@ ProcessCode EDM4hepSimHitReader::read(const AlgorithmContext& ctx) { simTrackerHit, [](const edm4hep::MCParticle& particle) { ActsFatras::Barcode result; - // TODO dont use podio internal id - result.setParticle(particle.id()); + result.setParticle( + EDM4hepUtil::podioObjectIDToInteger(particle.id())); return result; }, [&](std::uint64_t cellId) { diff --git a/Examples/Io/EDM4hep/src/EDM4hepUtil.cpp b/Examples/Io/EDM4hep/src/EDM4hepUtil.cpp index 7fe55d41bfe..aef0d6cb191 100644 --- a/Examples/Io/EDM4hep/src/EDM4hepUtil.cpp +++ b/Examples/Io/EDM4hep/src/EDM4hepUtil.cpp @@ -152,7 +152,8 @@ Measurement EDM4hepUtil::readMeasurement( // no need for digitization as we only want to identify the sensor Acts::GeometryIdentifier geometryId = geometryMapper(from.getCellID()); - IndexSourceLink sourceLink{geometryId, from.id()}; + IndexSourceLink sourceLink{ + geometryId, static_cast(podioObjectIDToInteger(from.id()))}; auto pos = from.getPosition(); auto cov = from.getCovMatrix(); From 3fd1e0f99f98e15c97b7f0d4c469b12aad2cb09f Mon Sep 17 00:00:00 2001 From: Benjamin Huth <37871400+benjaminhuth@users.noreply.github.com> Date: Wed, 22 Nov 2023 11:53:09 +0100 Subject: [PATCH 2/2] feat: Collect and write material statistics in GSF (#2695) This refactors the collection of material statistics in the GSF, so it only accumulates values until the last measurement in the forward pass (not collecting material in the solenoid etc...) Also writes the sum and the maximum material of the forward pass to the tracksummary. --- .../Acts/TrackFitting/GaussianSumFitter.hpp | 29 +++++++++--- Core/include/Acts/TrackFitting/GsfOptions.hpp | 4 ++ .../Acts/TrackFitting/detail/GsfActor.hpp | 45 +++++++++++++------ .../Acts/TrackFitting/detail/GsfUtils.hpp | 17 +++++++ .../TrackFitting/src/GsfFitterFunction.cpp | 9 ++++ .../Io/Root/RootTrackSummaryWriter.hpp | 5 +++ .../Io/Root/src/RootTrackSummaryWriter.cpp | 28 ++++++++++++ Examples/Python/src/Output.cpp | 8 ++-- Examples/Python/tests/root_file_hashes.txt | 4 +- Examples/Scripts/Python/truth_tracking_gsf.py | 1 + 10 files changed, 125 insertions(+), 25 deletions(-) diff --git a/Core/include/Acts/TrackFitting/GaussianSumFitter.hpp b/Core/include/Acts/TrackFitting/GaussianSumFitter.hpp index adc918b6eac..b2c6c46589c 100644 --- a/Core/include/Acts/TrackFitting/GaussianSumFitter.hpp +++ b/Core/include/Acts/TrackFitting/GaussianSumFitter.hpp @@ -305,7 +305,7 @@ struct GaussianSumFitter { return return_error_or_abort(fwdResult.error()); } - auto& fwdGsfResult = + const auto& fwdGsfResult = fwdResult->template get(); if (!fwdGsfResult.result.ok()) { @@ -321,8 +321,8 @@ struct GaussianSumFitter { ACTS_VERBOSE("- processed states: " << fwdGsfResult.processedStates); ACTS_VERBOSE("- measurement states: " << fwdGsfResult.measurementStates); - std::size_t nInvalidBetheHeitler = fwdGsfResult.nInvalidBetheHeitler; - double maxPathXOverX0 = fwdGsfResult.maxPathXOverX0; + std::size_t nInvalidBetheHeitler = fwdGsfResult.nInvalidBetheHeitler.val(); + double maxPathXOverX0 = fwdGsfResult.maxPathXOverX0.val(); ////////////////// // Backward pass @@ -408,8 +408,14 @@ struct GaussianSumFitter { GsfError::NoMeasurementStatesCreatedBackward); } - nInvalidBetheHeitler += bwdGsfResult.nInvalidBetheHeitler; - maxPathXOverX0 = std::max(maxPathXOverX0, bwdGsfResult.maxPathXOverX0); + // For the backward pass we want the counters at in end (= at the + // interaction point) and not at the last measurement surface + bwdGsfResult.nInvalidBetheHeitler.update(); + bwdGsfResult.maxPathXOverX0.update(); + bwdGsfResult.sumPathXOverX0.update(); + nInvalidBetheHeitler += bwdGsfResult.nInvalidBetheHeitler.val(); + maxPathXOverX0 = + std::max(maxPathXOverX0, bwdGsfResult.maxPathXOverX0.val()); if (nInvalidBetheHeitler > 0) { ACTS_WARNING("Encountered " << nInvalidBetheHeitler @@ -479,12 +485,23 @@ struct GaussianSumFitter { if (trackContainer.hasColumn( hashString(GsfConstants::kFinalMultiComponentStateColumn))) { - ACTS_DEBUG("Add final multi-component state to track") + ACTS_DEBUG("Add final multi-component state to track"); track.template component( GsfConstants::kFinalMultiComponentStateColumn) = std::move(params); } } + if (trackContainer.hasColumn( + hashString(GsfConstants::kFwdMaxMaterialXOverX0))) { + track.template component(GsfConstants::kFwdMaxMaterialXOverX0) = + fwdGsfResult.maxPathXOverX0.val(); + } + if (trackContainer.hasColumn( + hashString(GsfConstants::kFwdSumMaterialXOverX0))) { + track.template component(GsfConstants::kFwdSumMaterialXOverX0) = + fwdGsfResult.sumPathXOverX0.val(); + } + calculateTrackQuantities(track); return track; diff --git a/Core/include/Acts/TrackFitting/GsfOptions.hpp b/Core/include/Acts/TrackFitting/GsfOptions.hpp index 488d07db12f..09bbc4e940f 100644 --- a/Core/include/Acts/TrackFitting/GsfOptions.hpp +++ b/Core/include/Acts/TrackFitting/GsfOptions.hpp @@ -39,6 +39,10 @@ constexpr std::string_view kFinalMultiComponentStateColumn = "gsf-final-multi-component-state"; using FinalMultiComponentState = std::optional; +constexpr std::string_view kFwdSumMaterialXOverX0 = + "gsf-fwd-sum-material-x-over-x0"; +constexpr std::string_view kFwdMaxMaterialXOverX0 = + "gsf-fwd-max-material-x-over-x0"; } // namespace GsfConstants /// The extensions needed for the GSF diff --git a/Core/include/Acts/TrackFitting/detail/GsfActor.hpp b/Core/include/Acts/TrackFitting/detail/GsfActor.hpp index f348451b660..ca022bc8bc0 100644 --- a/Core/include/Acts/TrackFitting/detail/GsfActor.hpp +++ b/Core/include/Acts/TrackFitting/detail/GsfActor.hpp @@ -56,9 +56,10 @@ struct GsfResult { std::vector visitedSurfaces; std::vector surfacesVisitedBwdAgain; - /// Statistics about encounterings of invalid bethe heitler approximations - std::size_t nInvalidBetheHeitler = 0; - double maxPathXOverX0 = 0.0; + /// Statistics about material encounterings + Updatable nInvalidBetheHeitler; + Updatable maxPathXOverX0; + Updatable sumPathXOverX0; // Propagate potential errors to the outside Result result{Result::success()}; @@ -232,6 +233,15 @@ struct GsfActor { return; } + // Update the counters. Note that this should be done before potential + // material interactions, because if this is our last measurement this would + // not influence the fit anymore. + if (haveMeasurement) { + result.maxPathXOverX0.update(); + result.sumPathXOverX0.update(); + result.nInvalidBetheHeitler.update(); + } + for (auto cmp : stepper.componentIterable(state.stepping)) { auto singleState = cmp.singleState(state); cmp.singleStepper(stepper).transportCovarianceToBound( @@ -332,6 +342,7 @@ struct GsfActor { std::vector& componentCache, result_type& result) const { auto cmps = stepper.componentIterable(state.stepping); + double pathXOverX0 = 0.0; for (auto [idx, cmp] : zip(tmpStates.tips, cmps)) { auto proxy = tmpStates.traj.getTrackState(idx); @@ -339,18 +350,23 @@ struct GsfActor { proxy.filtered(), proxy.filteredCovariance(), stepper.particleHypothesis(state.stepping)); - applyBetheHeitler(state, navigator, bound, tmpStates.weights.at(idx), - componentCache, result); + pathXOverX0 += + applyBetheHeitler(state, navigator, bound, tmpStates.weights.at(idx), + componentCache, result); } + + // Store average material seen by the components + // Should not be too broadly distributed + result.sumPathXOverX0.tmp() += pathXOverX0 / tmpStates.tips.size(); } template - void applyBetheHeitler(const propagator_state_t& state, - const navigator_t& navigator, - const BoundTrackParameters& old_bound, - const double old_weight, - std::vector& componentCaches, - result_type& result) const { + double applyBetheHeitler(const propagator_state_t& state, + const navigator_t& navigator, + const BoundTrackParameters& old_bound, + const double old_weight, + std::vector& componentCaches, + result_type& result) const { const auto& surface = *navigator.currentSurface(state.navigation); const auto p_prev = old_bound.absoluteMomentum(); @@ -365,11 +381,12 @@ struct GsfActor { slab.scaleThickness(pathCorrection); const double pathXOverX0 = slab.thicknessInX0(); + result.maxPathXOverX0.tmp() = + std::max(result.maxPathXOverX0.tmp(), pathXOverX0); // Emit a warning if the approximation is not valid for this x/x0 if (!m_cfg.bethe_heitler_approx->validXOverX0(pathXOverX0)) { - ++result.nInvalidBetheHeitler; - result.maxPathXOverX0 = std::max(result.maxPathXOverX0, pathXOverX0); + ++result.nInvalidBetheHeitler.tmp(); ACTS_DEBUG( "Bethe-Heitler approximation encountered invalid value for x/x0=" << pathXOverX0 << " at surface " << surface.geometryId()); @@ -428,6 +445,8 @@ struct GsfActor { // Set the remaining things and push to vector componentCaches.push_back({new_weight, new_pars, new_cov}); } + + return pathXOverX0; } /// Remove components with low weights and renormalize from the component diff --git a/Core/include/Acts/TrackFitting/detail/GsfUtils.hpp b/Core/include/Acts/TrackFitting/detail/GsfUtils.hpp index 41e671a6fa8..49f61eed3f3 100644 --- a/Core/include/Acts/TrackFitting/detail/GsfUtils.hpp +++ b/Core/include/Acts/TrackFitting/detail/GsfUtils.hpp @@ -239,5 +239,22 @@ struct MultiTrajectoryProjector { } }; +/// Small Helper class that allows to carry a temporary value until we decide to +/// update the actual value. The temporary value is deliberately only accessible +/// with a mutable reference +template +class Updatable { + T m_tmp{}; + T m_val{}; + + public: + Updatable() : m_tmp(0), m_val(0) {} + + T &tmp() { return m_tmp; } + void update() { m_val = m_tmp; } + + const T &val() const { return m_val; } +}; + } // namespace detail } // namespace Acts diff --git a/Examples/Algorithms/TrackFitting/src/GsfFitterFunction.cpp b/Examples/Algorithms/TrackFitting/src/GsfFitterFunction.cpp index 636a6b74c3d..536c24c147a 100644 --- a/Examples/Algorithms/TrackFitting/src/GsfFitterFunction.cpp +++ b/Examples/Algorithms/TrackFitting/src/GsfFitterFunction.cpp @@ -78,6 +78,7 @@ struct GsfFitterFunctionImpl final : public ActsExamples::TrackFitterFunction { std::size_t maxComponents = 0; double weightCutoff = 0; + const double momentumCutoff = 0; // 500_MeV; bool abortOnError = false; bool disableAllMaterialHandling = false; MixtureReductionAlgorithm reductionAlg = @@ -146,6 +147,14 @@ struct GsfFitterFunctionImpl final : public ActsExamples::TrackFitterFunction { tracks.template addColumn(key); } + if (!tracks.hasColumn(Acts::hashString(kFwdMaxMaterialXOverX0))) { + tracks.template addColumn(std::string(kFwdMaxMaterialXOverX0)); + } + + if (!tracks.hasColumn(Acts::hashString(kFwdSumMaterialXOverX0))) { + tracks.template addColumn(std::string(kFwdSumMaterialXOverX0)); + } + return fitter.fit(sourceLinks.begin(), sourceLinks.end(), initialParameters, gsfOptions, tracks); } diff --git a/Examples/Io/Root/include/ActsExamples/Io/Root/RootTrackSummaryWriter.hpp b/Examples/Io/Root/include/ActsExamples/Io/Root/RootTrackSummaryWriter.hpp index d2dd29e15be..ab1a13e2457 100644 --- a/Examples/Io/Root/include/ActsExamples/Io/Root/RootTrackSummaryWriter.hpp +++ b/Examples/Io/Root/include/ActsExamples/Io/Root/RootTrackSummaryWriter.hpp @@ -67,6 +67,8 @@ class RootTrackSummaryWriter final : public WriterT { std::string fileMode = "RECREATE"; /// Switch for adding full covariance matrix to output file. bool writeCovMat = false; + /// Write GSF specific things (for now only some material statistics) + bool writeGsfSpecific = false; }; /// Constructor @@ -226,6 +228,9 @@ class RootTrackSummaryWriter final : public WriterT { std::vector m_cov_eT_eTHETA; std::vector m_cov_eT_eQOP; std::vector m_cov_eT_eT; + + std::vector m_gsf_max_material_fwd; + std::vector m_gsf_sum_material_fwd; }; } // namespace ActsExamples diff --git a/Examples/Io/Root/src/RootTrackSummaryWriter.cpp b/Examples/Io/Root/src/RootTrackSummaryWriter.cpp index 47fe85b20fd..356a5dfbecd 100644 --- a/Examples/Io/Root/src/RootTrackSummaryWriter.cpp +++ b/Examples/Io/Root/src/RootTrackSummaryWriter.cpp @@ -13,6 +13,7 @@ #include "Acts/EventData/MultiTrajectoryHelpers.hpp" #include "Acts/EventData/VectorMultiTrajectory.hpp" #include "Acts/Surfaces/Surface.hpp" +#include "Acts/TrackFitting/GsfOptions.hpp" #include "Acts/Utilities/Intersection.hpp" #include "Acts/Utilities/MultiIndex.hpp" #include "Acts/Utilities/Result.hpp" @@ -137,6 +138,12 @@ ActsExamples::RootTrackSummaryWriter::RootTrackSummaryWriter( m_outputTree->Branch("pull_eTHETA_fit", &m_pull_eTHETA_fit); m_outputTree->Branch("pull_eQOP_fit", &m_pull_eQOP_fit); m_outputTree->Branch("pull_eT_fit", &m_pull_eT_fit); + + if (m_cfg.writeGsfSpecific) { + m_outputTree->Branch("max_material_fwd", &m_gsf_max_material_fwd); + m_outputTree->Branch("sum_material_fwd", &m_gsf_sum_material_fwd); + } + if (m_cfg.writeCovMat == true) { // create one branch for every entry of covariance matrix // one block for every row of the matrix, every entry gets own branch @@ -440,6 +447,24 @@ ActsExamples::ProcessCode ActsExamples::RootTrackSummaryWriter::writeT( m_pull_eT_fit.push_back(pull[Acts::eBoundTime]); m_hasFittedParams.push_back(hasFittedParams); + + if (m_cfg.writeGsfSpecific) { + using namespace Acts::GsfConstants; + if (tracks.hasColumn(Acts::hashString(kFwdMaxMaterialXOverX0))) { + m_gsf_max_material_fwd.push_back( + track.template component(kFwdMaxMaterialXOverX0)); + } else { + m_gsf_max_material_fwd.push_back(NaNfloat); + } + + if (tracks.hasColumn(Acts::hashString(kFwdSumMaterialXOverX0))) { + m_gsf_sum_material_fwd.push_back( + track.template component(kFwdSumMaterialXOverX0)); + } else { + m_gsf_sum_material_fwd.push_back(NaNfloat); + } + } + if (m_cfg.writeCovMat) { // write all entries of covariance matrix to output file // one branch for every entry of the matrix. @@ -549,6 +574,9 @@ ActsExamples::ProcessCode ActsExamples::RootTrackSummaryWriter::writeT( m_pull_eQOP_fit.clear(); m_pull_eT_fit.clear(); + m_gsf_max_material_fwd.clear(); + m_gsf_sum_material_fwd.clear(); + if (m_cfg.writeCovMat) { m_cov_eLOC0_eLOC0.clear(); m_cov_eLOC0_eLOC1.clear(); diff --git a/Examples/Python/src/Output.cpp b/Examples/Python/src/Output.cpp index e5dde8abcea..90b2915540e 100644 --- a/Examples/Python/src/Output.cpp +++ b/Examples/Python/src/Output.cpp @@ -323,10 +323,10 @@ void addOutput(Context& ctx) { inputTracks, inputParticles, inputSimHits, inputMeasurementParticlesMap, inputMeasurementSimHitsMap, filePath, treeName, fileMode); - ACTS_PYTHON_DECLARE_WRITER(ActsExamples::RootTrackSummaryWriter, mex, - "RootTrackSummaryWriter", inputTracks, - inputParticles, inputMeasurementParticlesMap, - filePath, treeName, fileMode, writeCovMat); + ACTS_PYTHON_DECLARE_WRITER( + ActsExamples::RootTrackSummaryWriter, mex, "RootTrackSummaryWriter", + inputTracks, inputParticles, inputMeasurementParticlesMap, filePath, + treeName, fileMode, writeCovMat, writeGsfSpecific); ACTS_PYTHON_DECLARE_WRITER( ActsExamples::VertexPerformanceWriter, mex, "VertexPerformanceWriter", diff --git a/Examples/Python/tests/root_file_hashes.txt b/Examples/Python/tests/root_file_hashes.txt index ade2c97a26f..5ca47d7a7d0 100644 --- a/Examples/Python/tests/root_file_hashes.txt +++ b/Examples/Python/tests/root_file_hashes.txt @@ -35,9 +35,9 @@ test_truth_tracking_kalman[odd-1000.0]__trackstates_fitter.root: f616b2234db239b test_truth_tracking_kalman[odd-1000.0]__tracksummary_fitter.root: a53253b509b5779a5d856f8c09d76a499b217b55ba4b0e52d9076ffad726463f test_truth_tracking_kalman[odd-1000.0]__performance_track_finder.root: 39aec6316cceb90e314e16b02947faa691c18f57c3a851a25e547a8fc05a4593 test_truth_tracking_gsf[generic]__trackstates_gsf.root: 8422ad83843d69aac7d1e37d6d3e4019b1abfaf854d9246c1fb8a4fc5d6c3f13 -test_truth_tracking_gsf[generic]__tracksummary_gsf.root: fca19c92a21bec029325189d7a404e5da4851b6cca886170e14486863b86c29e +test_truth_tracking_gsf[generic]__tracksummary_gsf.root: fe5bce60f13f6221290a68340772beefe5bb5389f424b23c14787a944aa5d84f test_truth_tracking_gsf[odd]__trackstates_gsf.root: 9e1c0acaf1aed2347992c0a4c076ca27ce71d770521b33680f30f9734dd8494e -test_truth_tracking_gsf[odd]__tracksummary_gsf.root: ca78cba367449062e985ee5aacd4db88ef721cbe56a772b818f1586770ec05ba +test_truth_tracking_gsf[odd]__tracksummary_gsf.root: 09e06b254b18e18cb3086038e2f8b7ac96a33d980ac14395bc9a02227c4433bd test_particle_gun__particles.root: 8549ba6e20338004ab8ba299fc65e1ee5071985b46df8f77f887cb6fef56a8ec test_material_mapping__material-map_tracks.root: 2deac7a48ff1185ba3889cfd4cc0bd0a6de57293048dfd6766f3c3907232e45e test_material_mapping__propagation-material.root: 84b04ebd5721550867f0f193b5eb4e9f3f205df542cb46090d320be6bff565c5 diff --git a/Examples/Scripts/Python/truth_tracking_gsf.py b/Examples/Scripts/Python/truth_tracking_gsf.py index 75e54d2557a..f7738585506 100755 --- a/Examples/Scripts/Python/truth_tracking_gsf.py +++ b/Examples/Scripts/Python/truth_tracking_gsf.py @@ -139,6 +139,7 @@ def runTruthTrackingGsf( inputParticles="truth_seeds_selected", inputMeasurementParticlesMap="measurement_particles_map", filePath=str(outputDir / "tracksummary_gsf.root"), + writeGsfSpecific=True, ) )