New algorithm to decorate Acts track particles with truth links.

The original TrackFindingValidationAlg, which already computes truth
matching probability etc.,  is split into a base now used by the
TrackFindingValidationAlg and the new TrackParticleTruthDecorationAlg.
The latter decorates track particles with links to truth particles and
the associated matching probability, hit efficiency and hit purity.

Moreover the base algorithm is modified in such a way that the
statistics gathering is compile time optional. Although, it is
currently enabled for both algorithms.

Tracking/Acts/ActsConfig/python/ActsTruthConfig.py

@@ -97,18 +97,13 @@ def ActsTruthAssociationAlgCfg(flags,
         acc.merge(ActsStripClusterToTruthAssociationAlgCfg(flags, **extractChildKwargs(prefix="StripClusterToTruthAssociationAlg.", **kwargs) ))
         
     return acc
-    
-    
-def ActsTrackFindingValidationAlgCfg(flags,
-                                     name: str = 'ActsTracksValidationAlg',
-                                     **kwargs: dict) -> ComponentAccumulator:
-    acc = ComponentAccumulator()
-    kwargs.setdefault('TruthParticleHitCounts','TruthParticleHitCounts')
-    kwargs.setdefault('TrackToTruthAssociationMap','ActsTracksToTruthParticles')
+
+def setDefaultTruthMatchingArgs(kwargs) :
     kwargs.setdefault('MatchWeights',[0.,               # other
                                       10., 5.,           # ID (pixel, strips)
                                       0.,  0., 0. , 0.,  # MS
                                       0. ])             # HGTD
+    # weights used for hit purity and hit efficiencies
     kwargs.setdefault('CountWeights',[0.,               # other
                                       1.,1.,            # ID (pixel, strips)
                                       0., 0., 0. , 0.,  # MS
@@ -117,6 +112,40 @@ def ActsTrackFindingValidationAlgCfg(flags,
     kwargs.setdefault('ShowDetailedTables',False)
     kwargs.setdefault('PdgIdCategorisation',False)
     kwargs.setdefault('StatisticEtaBins',[eta/10. for eta in range(5, 40, 5)])
+
+def ActsTrackParticleTruthDecorationAlgCfg(flags,
+                                           name: str = 'ActsTrackParticleTruthDecorationAlg',
+                                           **kwargs: dict) -> ComponentAccumulator:
+    acc = ComponentAccumulator()
+    kwargs.setdefault('TrackToTruthAssociationMaps','ActsCombinedTracksToTruthParticleAssociation')
+    kwargs.setdefault('TrackParticleContainerName','ActsCombinedTracksParticlesAlt')
+    kwargs.setdefault('TruthParticleHitCounts','TruthParticleHitCounts')
+    # weights used for computing the matching probability and identifying the best match
+    setDefaultTruthMatchingArgs(kwargs)
+    kwargs.setdefault('ComputeTrackRecoEfficiency',False)
+
+    if 'TruthSelectionTool' not in kwargs:
+        # should be as tight or looser as the TruthSelectionTool when analysing the truth matches
+        from InDetPhysValMonitoring.InDetPhysValMonitoringConfig import InDetRttTruthSelectionToolCfg
+        kwargs.setdefault("TruthSelectionTool", acc.popToolsAndMerge(
+            InDetRttTruthSelectionToolCfg(flags,
+                                          name='RelaxedInDetRttTruthSelectionTool',
+                                          requireOnlyPrimary=False,
+                                          minPt=500.,
+                                          maxEta=4.5
+                                          )))
+
+    acc.addEventAlgo( CompFactory.ActsTrk.TrackParticleTruthDecorationAlg(name=name, **kwargs) )
+    return acc
+
+def ActsTrackFindingValidationAlgCfg(flags,
+                                     name: str = 'ActsTracksValidationAlg',
+                                     **kwargs: dict) -> ComponentAccumulator:
+    acc = ComponentAccumulator()
+    kwargs.setdefault('TruthParticleHitCounts','TruthParticleHitCounts')
+    kwargs.setdefault('TrackToTruthAssociationMap','ActsTracksToTruthParticles')
+    setDefaultTruthMatchingArgs(kwargs)
+    kwargs.setdefault('ComputeTrackRecoEfficiency',True)
     
     if 'TruthSelectionTool' not in kwargs:
         from InDetPhysValMonitoring.InDetPhysValMonitoringConfig import InDetRttTruthSelectionToolCfg

Tracking/Acts/ActsTruth/src/TrackFindingValidationAlg.h

@@ -14,178 +14,24 @@
 
 // Handle Keys
 #include "StoreGate/ReadHandleKey.h"
-#include "StoreGate/WriteHandleKey.h"
 
+#include "TrackTruthMatchingBaseAlg.h"
 #include "ActsEvent/TrackToTruthParticleAssociation.h"
-#include "ActsEvent/TruthParticleHitCounts.h"
-#include "ActsGeometryInterfaces/IActsTrackingGeometryTool.h"
-#include "TrkTruthTrackInterfaces/IAthSelectionTool.h"
-
-#include <mutex>
-#include "ActsInterop/StatUtils.h"
-#include "ActsInterop/TableUtils.h"
-
-#include <string>
-#include <memory>
-#include <array>
-#include <atomic>
-#include <type_traits>
-
-#include <cmath>
-#include <iomanip>
-#include <ostream>
-#include <string>
-#include <sstream>
-#include <vector>
 
 namespace ActsTrk
 {
-  constexpr bool TrackFindingValidationDebugHists = false;
-
-  class TrackFindingValidationAlg : public AthReentrantAlgorithm
+  class TrackFindingValidationAlg : public TrackTruthMatchingBaseAlg
   {
   public:
-    TrackFindingValidationAlg(const std::string &name,
-                               ISvcLocator *pSvcLocator);
+    using TrackTruthMatchingBaseAlg::TrackTruthMatchingBaseAlg;
 
     virtual StatusCode initialize() override;
     virtual StatusCode finalize() override;
     virtual StatusCode execute(const EventContext &ctx) const override;
 
   private:
-     SG::ReadHandleKey<TruthParticleHitCounts>  m_truthHitCounts
-        {this, "TruthParticleHitCounts","", "Map from truth particle to hit counts." };
      SG::ReadHandleKey<TrackToTruthParticleAssociation>  m_trackToTruth
         {this, "TrackToTruthAssociationMap","", "Association map from tracks to generator particles." };
-
-     Gaudi::Property<std::vector<float> > m_weightsForProb
-        {this, "MatchWeights", {}, "Weights applied to the counts per measurement type for weighted sums"
-                                  " which are used compute the match probability." };
-     Gaudi::Property<std::vector<float> > m_weights
-        {this, "CountWeights", {}, "Weights applied to the counts per measurement type for weighted sums"
-                                   " which are used to compute hit efficiencies and purities." };
-
-     Gaudi::Property<std::vector<float>> m_statEtaBins
-        {this, "StatisticEtaBins", {-4, -2.6, -2, 0, 2., 2.6, 4}, "Gather statistics separately for these eta bins."};
-     Gaudi::Property<std::vector<float>> m_statPtBins
-        {this, "StatisticPtBins", {1.e3,2.5e3,10e3, 100e3}, "Gather statistics separately for these pt bins."};
-     Gaudi::Property<bool> m_pdgIdCategorisation
-        {this, "PdgIdCategorisation", false, "Categorise by pdg id."};
-     Gaudi::Property<bool> m_showRawCounts
-        {this, "ShowRawCounts", false, "Show all counters."};
-     Gaudi::Property<bool> m_printDetails
-        {this, "ShowDetailedTables", false, "Show more details; stat. uncert., RMS, entries"};
-
-    ToolHandle<IAthSelectionTool> m_truthSelectionTool{this, "TruthSelectionTool","AthTruthSelectionTool", "Truth selection tool (for efficiencies and resolutions)"};
-
-     // helper struct for compile time optional statistics
-     template <bool IsDebug>
-     struct DebugCounter {
-        struct Empty {
-           template <typename... T_Args>
-           Empty(T_Args... ) {}
-        };
-        mutable typename std::conditional<IsDebug,
-                                          std::mutex,
-                                          Empty>::type m_mutex ATLAS_THREAD_SAFE;
-        mutable typename std::conditional<IsDebug,
-                                          ActsUtils::StatHist,
-                                          Empty>::type m_measPerTruthParticleWithoutCounts ATLAS_THREAD_SAFE {20,-.5,20.-.5};
-        mutable typename std::conditional<IsDebug,
-                                          ActsUtils::StatHist,
-                                          Empty>::type m_bestMatchProb ATLAS_THREAD_SAFE {20,0.,1.};
-        mutable typename std::conditional<IsDebug,
-                                          ActsUtils::StatHist,
-                                          Empty>::type m_nextToBestMatchProb ATLAS_THREAD_SAFE {20,0.,1.};
-
-        template <class T_OutStream>
-        void dumpStatistics(T_OutStream &out) const;
-        void fillMeasForTruthParticleWithoutCount(double weighted_measurement_sum) const;
-        void fillTruthMatchProb(const std::array<float,2> &best_match_prob) const;
-     };
-     DebugCounter<TrackFindingValidationDebugHists> m_debugCounter;
-
-     // s_NMeasurementTypes is equal to the number of UncalibMeasType, but we have to remove the "Other" option, which
-     // corresponds to an unknown type
-     constexpr static unsigned int s_NMeasurementTypes = static_cast<unsigned int>(xAOD::UncalibMeasType::nTypes) - 1u;
-
-     // statistics counter
-     enum ECounter {
-        NTracksTotal,
-        NTruthWithCountsTotal,
-        MissingTruthParticleHitCounts,
-        NoAssociatedTruthParticle,
-        NoSelectedTruthParticle,
-        TruthParticleNoNoiseMismatch,
-        kNCounter
-     };
-     //     mutable std::array< std::atomic<std::size_t>, kNCounter > m_counter ATLAS_THREAD_SAFE {};
-     mutable std::array< std::size_t, kNCounter > m_counter ATLAS_THREAD_SAFE {};
-
-     enum ECategorisedCounter {
-        kNTotalParticles,
-        kNParticleWithAssociatedTrack,
-        kNParticleWithMultipleAssociatedTracks,
-        kNTotalTracks,
-        kNCategorisedCounter
-     };
-     enum ECategorisedStat {
-        kHitEfficiency,
-        kHitPurity,
-        kMatchProbability,
-        kNCategorisedStat
-     };
-
-     constexpr static int s_pdgIdMax = 1000000000;  // categorise all truth particles with this or larger PDG ID as "Other"
-
-     mutable std::mutex m_statMutex ATLAS_THREAD_SAFE;
-
-     mutable std::vector< std::array< std::size_t, kNCategorisedCounter> >  m_counterPerEta ATLAS_THREAD_SAFE;
-     mutable std::vector< int >                                             m_pdgId ATLAS_THREAD_SAFE;
-     mutable std::vector< std::array< std::size_t, kNCategorisedCounter> >  m_counterPerPdgId ATLAS_THREAD_SAFE;
-     mutable std::vector< std::array< ActsUtils::Stat, kNCategorisedStat> > m_statPerEta ATLAS_THREAD_SAFE;
-     mutable std::vector< std::array< ActsUtils::Stat, kNCategorisedStat> > m_statPerPdgId ATLAS_THREAD_SAFE;
-     mutable ActsUtils::StatHist                                            m_truthSelectionCuts ATLAS_THREAD_SAFE;
-
-     bool m_useAbsEtaForStat = false;
-
-     /// @brief check that bins are in increasing order.
-     /// will cause a FATAL error if bins are not in increasing order.
-     void checkBinOrder( const std::vector<float> &bin_edges, const std::string &bin_label) const;
-
-     /// @brief Return the category based on the provided eta value
-     /// @param pt the pt of the truth particle
-     /// @param eta the eta of the truth particle
-     /// @return a bin assigned to the give eta value
-     std::size_t getPtEtaStatCategory(float pt, float eta) const;
-
-     /// @brief Return the category based on the PDG ID
-     /// @param pt the pt of the truth particle
-     /// @param pdg_id the PDG ID
-     /// @return 0 or a slot associated to a single PDG ID (absolut value)
-     std::size_t getPtPdgIdStatCategory(float pt, int pdg_id) const;
-     void initStatTables();
-     void printStatTables() const;
-     void printCategories(const std::vector<std::string> &pt_category_labels,
-                          const std::vector<std::string> &eta_category_labels,
-                          std::vector<std::string> &counter_labels,
-                          std::vector< std::array< ActsUtils::Stat, kNCategorisedStat> > &stat_per_category,
-                          std::vector< std::array< std::size_t, kNCategorisedCounter> >  &counts_per_category,
-                          const std::string &top_left_label,
-                          bool print_sub_categories) const;
-     void printData2D(const std::vector<std::string> &row_category_labels,
-                      const std::vector<std::string> &col_category_labels,
-                      const std::string &top_left_label,
-                      std::vector< std::array< ActsUtils::Stat, kNCategorisedStat> > &stat_per_category,
-                      std::vector< std::array< std::size_t, kNCategorisedCounter> >  &counts_per_category,
-                      bool rotate) const;
-
-     StatusCode checkMatchWeights();
-
-     static  double weightedCountSum(const ActsTrk::HitCounterArray &counts,
-                                     const std::vector<float> &weights);
-     static  double noiseCorrection(const ActsTrk::HitCounterArray &noise_counts,
-                                    const std::vector<float> &weights);
   };
 
 } // namespace

Tracking/Acts/ActsTruth/src/TrackParticleTruthDecorationAlg.cxx

+/*
+  Copyright (C) 2002-2024 CERN for the benefit of the ATLAS collaboration
+*/
+#undef NDEBUG
+#include "TrackParticleTruthDecorationAlg.h"
+#include "ActsEvent/TrackContainer.h"
+#include "xAODTruth/TruthParticle.h"
+#include "xAODTruth/TruthParticleContainer.h"
+#include <unordered_map>
+#include "decoratorUtils.h"
+
+
+namespace ActsTrk
+{
+
+  StatusCode TrackParticleTruthDecorationAlg::initialize()
+  {
+     StatusCode sc = TrackTruthMatchingBaseAlg::initialize();
+     ATH_CHECK( m_trackToTruth.initialize() );
+     ATH_CHECK( m_trkParticleName.initialize() );
+     std::vector<std::string> float_decor_names(kNFloatDecorators);
+     float_decor_names[kMatchingProbability]="truthMatchProbability";
+     float_decor_names[kHitPurity]="truthHitPurity";
+     float_decor_names[kHitEfficiency]="truthHitEfficiency";
+     createDecoratorKeys(*this,m_trkParticleName,"" /*prefix ? */, float_decor_names,m_floatDecor);
+     assert( m_floatDecor.size() == kNFloatDecorators);
+     std::vector<std::string> link_decor_names;
+     link_decor_names.push_back("truthParticleLink");
+     createDecoratorKeys(*this,m_trkParticleName,"" /*prefix ? */, link_decor_names,m_linkDecor);
+     assert( m_linkDecor.size() == 1);
+     return sc;
+  }
+
+  StatusCode TrackParticleTruthDecorationAlg::finalize()
+  {
+     StatusCode sc = TrackTruthMatchingBaseAlg::finalize();
+     return sc;
+  }
+
+  StatusCode TrackParticleTruthDecorationAlg::execute(const EventContext &ctx) const
+  {
+    const TruthParticleHitCounts &truth_particle_hit_counts = getTruthParticleHitCounts(ctx);
+    // @TODO or use simply a vector ?
+    std::unordered_map<const ActsTrk::TrackContainerBase *, const ActsTrk::TrackToTruthParticleAssociation *> truth_association_map;
+    truth_association_map.reserve( m_trackToTruth.size());
+    for (const SG::ReadHandleKey<TrackToTruthParticleAssociation> &truth_association_key : m_trackToTruth) {
+       SG::ReadHandle<TrackToTruthParticleAssociation> track_to_truth_handle = SG::makeHandle(truth_association_key, ctx);
+       if (!track_to_truth_handle.isValid()) {
+          ATH_MSG_ERROR("No track to truth particle association for key " << truth_association_key.key() );
+          return StatusCode::FAILURE;
+       }
+       truth_association_map.insert(std::make_pair( track_to_truth_handle->sourceContainer(), track_to_truth_handle.cptr() ));
+    }
+
+    SG::ReadHandle<xAOD::TrackParticleContainer> track_particle_handle = SG::makeHandle(m_trkParticleName, ctx);
+    if (!track_particle_handle.isValid()) {
+       ATH_MSG_ERROR("No track particle container for key " << track_particle_handle.key() );
+       return StatusCode::FAILURE;
+    }
+    std::vector< SG::WriteDecorHandle<xAOD::TrackParticleContainer,float > >
+       float_decor( createDecorators<xAOD::TrackParticleContainer, float >(m_floatDecor, ctx) );
+    SG::WriteDecorHandle<xAOD::TrackParticleContainer, ElementLink<xAOD::TruthParticleContainer> >
+       link_decor(m_linkDecor.at(0), ctx);
+
+    EventStat event_stat(truthSelectionTool(),
+                         perEtaSize(),
+                         perPdgIdSize(),
+                         track_particle_handle->size());
+
+    static const SG::AuxElement::ConstAccessor<ElementLink<ActsTrk::TrackContainer> > actsTrackLink("actsTrack");
+
+    std::pair<const ActsTrk::TrackContainerBase *, const ActsTrk::TrackToTruthParticleAssociation *>
+       the_track_truth_association{ nullptr, nullptr};
+    ElementLink<xAOD::TruthParticleContainer> ref_truth_link;
+    for(const xAOD::TrackParticle *track_particle : *track_particle_handle) {
+       ElementLink<ActsTrk::TrackContainer> link_to_track = actsTrackLink(*track_particle);
+       static_assert( std::is_same<ElementLink<ActsTrk::TrackContainer>::ElementConstReference,
+                      std::optional<ActsTrk::TrackContainer::ConstTrackProxy> >::value);
+       TruthMatchResult truth_match{} ;
+
+       {
+          std::optional<ActsTrk::TrackContainer::ConstTrackProxy> optional_track = *link_to_track;
+          if (optional_track.has_value()) {
+             const ActsTrk::TrackContainerBase *track_container = &(optional_track.value().container());
+             if (track_container != the_track_truth_association.first && track_container ) {
+                std::unordered_map<const ActsTrk::TrackContainerBase *, const ActsTrk::TrackToTruthParticleAssociation *>::const_iterator
+                   truth_association_map_iter = truth_association_map.find( track_container );
+                if (truth_association_map_iter != truth_association_map.end()) {
+                   the_track_truth_association = *truth_association_map_iter;
+                }
+             }
+             if (the_track_truth_association.second) {
+                truth_match =  analyseTrackTruth(truth_particle_hit_counts,
+                                                 (*the_track_truth_association.second).at(optional_track.value().index()),
+                                                 event_stat);
+
+                const xAOD::TruthParticle *truth_particle = truth_match.m_truthParticle;
+
+                // decorate track particle with link to truth particle, matching probability etc.
+                if (truth_particle) {
+                   if (!ref_truth_link.isValid()) {
+                      const xAOD::TruthParticleContainer *truth_particle_container
+                         = dynamic_cast<const xAOD::TruthParticleContainer *>(truth_particle->container());
+                      if (!truth_particle_container) {
+                         ATH_MSG_ERROR("Valid truth particle not part of a xAOD::TruthParticleContainer");
+                      }
+                      else {
+                         ref_truth_link=  ElementLink<xAOD::TruthParticleContainer>(*truth_particle_container,0u,ctx);
+                      }
+                   }
+                   assert( truth_particle->container() == ref_truth_link.getStorableObjectPointer() );
+                   link_decor(*track_particle) = ElementLink<xAOD::TruthParticleContainer>(ref_truth_link, truth_particle->index());
+                }
+                else {
+                   link_decor(*track_particle) = ElementLink<xAOD::TruthParticleContainer>();
+                }
+             }
+          }
+       }
+       float_decor[kMatchingProbability](*track_particle) = truth_match.m_matchProbability;
+       float_decor[kHitPurity](*track_particle) = truth_match.m_hitPurity;
+       float_decor[kHitEfficiency](*track_particle) = truth_match.m_hitEfficiency;
+    }
+    postProcessEventStat(truth_particle_hit_counts,
+                         track_particle_handle->size(),
+                         event_stat);
+    return StatusCode::SUCCESS;
+  }
+
+}

Tracking/Acts/ActsTruth/src/TrackParticleTruthDecorationAlg.h

+/*
+  Copyright (C) 2002-2024 CERN for the benefit of the ATLAS collaboration
+*/
+
+#ifndef ACTSTRK_TRACKPARTICLETRUTHDECORATIONALG_H
+#define ACTSTRK_TRACKPARTICLETRUTHDECORATIONALG_H 1
+
+#undef NDEBUG
+// Base Class
+#include "AthenaBaseComps/AthReentrantAlgorithm.h"
+
+// Gaudi includes
+#include "Gaudi/Property.h"
+
+// Handle Keys
+#include "StoreGate/ReadHandleKey.h"
+
+#include "TrackTruthMatchingBaseAlg.h"
+#include "ActsEvent/TrackToTruthParticleAssociation.h"
+#include "xAODTracking/TrackParticleContainer.h"
+
+namespace ActsTrk
+{
+  class TrackParticleTruthDecorationAlg : public TrackTruthMatchingBaseAlg
+  {
+  public:
+    using TrackTruthMatchingBaseAlg::TrackTruthMatchingBaseAlg;
+
+    virtual StatusCode initialize() override;
+    virtual StatusCode finalize() override;
+    virtual StatusCode execute(const EventContext &ctx) const override;
+
+  private:
+     SG::ReadHandleKeyArray<TrackToTruthParticleAssociation>  m_trackToTruth
+        {this, "TrackToTruthAssociationMaps",{},
+         "Association maps from tracks to generator particles for all Acts tracks linked from the given track particle." };
+
+     SG::ReadHandleKey<xAOD::TrackParticleContainer>  m_trkParticleName
+        {this,"TrackParticleContainerName", "InDetTrackParticles",""};
+
+     enum FloatDecorations { kMatchingProbability, kHitPurity, kHitEfficiency, kNFloatDecorators};
+     std::vector<SG::WriteDecorHandleKey<xAOD::TrackParticleContainer> >   m_linkDecor;
+     std::vector<SG::WriteDecorHandleKey<xAOD::TrackParticleContainer> >   m_floatDecor;
+  };
+
+} // namespace
+
+#endif

Tracking/Acts/ActsTruth/src/TrackTruthMatchingBaseAlg.h

+/*
+  Copyright (C) 2002-2024 CERN for the benefit of the ATLAS collaboration
+*/
+
+#ifndef ACTSTRK_TRACKTRUTHMATCHINGBASEALG_H
+#define ACTSTRK_TRACKTRUTHMATCHINGBASEALG_H 1
+
+// Base Class
+#include "AthenaBaseComps/AthReentrantAlgorithm.h"
+
+// Gaudi includes
+#include "Gaudi/Property.h"
+
+// Handle Keys
+#include "StoreGate/ReadHandleKey.h"
+#include "StoreGate/WriteHandleKey.h"
+
+#include "ActsEvent/TrackToTruthParticleAssociation.h"
+#include "ActsEvent/TruthParticleHitCounts.h"
+#include "ActsGeometryInterfaces/IActsTrackingGeometryTool.h"
+#include "TrkTruthTrackInterfaces/IAthSelectionTool.h"
+
+#include <mutex>
+#include "ActsInterop/StatUtils.h"
+#include "ActsInterop/TableUtils.h"
+
+#include <string>
+#include <memory>
+#include <array>
+#include <atomic>
+#include <type_traits>
+
+#include <cmath>
+#include <iomanip>
+#include <ostream>
+#include <sstream>
+#include <vector>
+
+namespace ActsTrk
+{
+  constexpr bool TrackFindingValidationDebugHists = false;
+  constexpr bool TrackFindingValidationDetailedStat = true;
+
+  class TrackTruthMatchingBaseAlg : public AthReentrantAlgorithm
+  {
+     template <bool DetailEnabled>
+     struct BaseStat;
+
+     template <bool DetailEnabled>
+     friend struct BaseStat;
+
+  public:
+    TrackTruthMatchingBaseAlg(const std::string &name,
+                               ISvcLocator *pSvcLocator);
+
+    virtual StatusCode initialize() override;
+    virtual StatusCode finalize() override;
+
+  protected:
+     const TruthParticleHitCounts &getTruthParticleHitCounts(const EventContext &ctx) const {
+        SG::ReadHandle<TruthParticleHitCounts> truth_particle_hit_counts_handle = SG::makeHandle(m_truthHitCounts, ctx);
+        if (!truth_particle_hit_counts_handle.isValid()) {
+           ATH_MSG_ERROR("No truth particle hit count map for key " << m_truthHitCounts.key() );
+           std::runtime_error("Failed to get truth particle hit count map");
+        }
+        return *truth_particle_hit_counts_handle;
+     }
+
+     const IAthSelectionTool &truthSelectionTool() const { return *m_truthSelectionTool; }
+
+     std::size_t perEtaSize() const  {
+        return m_detailedStat.perEtaSize();
+     }
+     std::size_t perPdgIdSize() const {
+        return m_detailedStat.perPdgIdSize();
+     }
+
+     template <bool DetailEnabled>
+     struct EventStatBase : public BaseStat<DetailEnabled> {
+        static constexpr bool doDetail = DetailEnabled;
+        
+        EventStatBase(const IAthSelectionTool &truth_selection_tool,
+                  std::size_t per_eta_size,
+                  std::size_t per_pdg_size,
+                  [[maybe_unused]] std::size_t track_to_truth_size)
+           : BaseStat<DetailEnabled>(truth_selection_tool, per_eta_size, per_pdg_size),
+             m_nTruthCuts(truth_selection_tool.nCuts())
+        {
+           if constexpr(DetailEnabled) {
+              m_truthParticlesWithAssociatedTrack.reserve(track_to_truth_size);
+           }
+
+        };
+        void fill([[maybe_unused]] unsigned int eta_category_i,
+                  [[maybe_unused]] unsigned int pdg_id_category_i,
+                  [[maybe_unused]] float hit_efficiency,
+                  [[maybe_unused]] float hit_purity,
+                  [[maybe_unused]] float match_prob,
+                  [[maybe_unused]] const xAOD::TruthParticle *best_match) {
+           BaseStat<DetailEnabled>::fill(eta_category_i,
+                                         pdg_id_category_i,
+                                         hit_efficiency,
+                                         hit_purity,
+                                         match_prob);
+           if constexpr(DetailEnabled) {
+              if (!m_truthParticlesWithAssociatedTrack.insert(best_match).second) {
+                 // truth particle already had a best match
+                 ++this->m_counterPerEta[eta_category_i][kNParticleWithMultipleAssociatedTracks];
+                 ++this->m_counterPerPdgId[pdg_id_category_i][kNParticleWithMultipleAssociatedTracks];
+              }
+              else {
+                 ++this->m_counterPerEta[eta_category_i][kNParticleWithAssociatedTrack];
+                 ++this->m_counterPerPdgId[pdg_id_category_i][kNParticleWithAssociatedTrack];
+              }
+           }
+        }
+        using TruthParticleSet = std::conditional< DetailEnabled,
+                                                   std::unordered_set<const xAOD::TruthParticle *>,
+                                                   typename BaseStat<DetailEnabled>::Empty >::type;
+        TruthParticleSet m_truthParticlesWithAssociatedTrack;
+
+        unsigned int m_nTruthParticleWithoutAssociatedCounts =0u;
+        unsigned int m_nTracksWithoutAssociatedTruthParticle =0u;
+        unsigned int m_nTracksWithoutSelectedTruthParticle =0u;
+        unsigned int m_nTruthParticleNonoiseMismatches=0u;
+
+        unsigned int m_nTruthCuts;
+     };
+     using EventStat = EventStatBase<TrackFindingValidationDetailedStat>;
+
+     /** Match result returned by @ref analyseTrackTruth
+      */
+     struct TruthMatchResult {
+        const xAOD::TruthParticle *m_truthParticle;  ///< best matching truth particle or nullptr
+        float m_matchProbability;                    ///< the matching probability based on weighted hit sums
+        float m_hitPurity;                           ///< fraction of hits originting from best match over total reco hits
+        float m_hitEfficiency;                       ///< fraction of hits originting from best match over total best match hits
+     };
+     /**
+         @param return tuple containing pointer to best matching truth particle, matching probability, hit purity and hit efficiency.
+     */
+     TruthMatchResult
+     analyseTrackTruth(const TruthParticleHitCounts &truth_particle_hit_counts,
+                       const HitCountsPerTrack &track_hit_counts,
+                       EventStat &event_stat) const;
+
+     void postProcessEventStat(const TruthParticleHitCounts &truth_particle_hit_counts,
+                               std::size_t n_tracks,
+                               EventStat &event_stat) const;
+     
+  private:
+     SG::ReadHandleKey<TruthParticleHitCounts>  m_truthHitCounts
+        {this, "TruthParticleHitCounts","", "Map from truth particle to hit counts." };
+
+     Gaudi::Property<std::vector<float> > m_weightsForProb
+        {this, "MatchWeights", {}, "Weights applied to the counts per measurement type for weighted sums"
+                                  " which are used compute the match probability." };
+     Gaudi::Property<std::vector<float> > m_weights
+        {this, "CountWeights", {}, "Weights applied to the counts per measurement type for weighted sums"
+                                   " which are used to compute hit efficiencies and purities." };
+
+     // Empty struct emulating the interface of a Gaudi property to replace optional properties if disabled.
+     template <class Base>
+     struct DummyProperty {
+        template <class OWNER>
+        DummyProperty( OWNER*, std::string, Base&&, std::string ) {}
+        Base value() const { return Base{}; }
+        operator Base() const { return Base{}; }
+        
+        // some dummy implementations for vector properties:
+        std::size_t size() const {return 0u; }
+        bool empty() const {return true; }
+        auto operator[](std::size_t idx) { return value().at(idx);}
+        auto begin() const {return value().end(); }
+        auto end() const {return value().end(); }
+
+        // Delegate operator() to the value
+        template <class... Args>
+        decltype( std::declval<Base>()( std::declval<Args&&>()... ) ) operator()( Args&&... args ) const
+           noexcept( noexcept( std::declval<Base>()( std::declval<Args&&>()... ) ) ) {
+           return value()( std::forward<Args>( args )... );
+        }
+     };
+
+     template <class Base>
+     using Property = std::conditional< TrackFindingValidationDetailedStat, Gaudi::Property<Base>, DummyProperty<Base> >::type;
+     
+     Property<std::vector<float>> m_statEtaBins
+        {this, "StatisticEtaBins", {-4, -2.6, -2, 0, 2., 2.6, 4}, "Gather statistics separately for these eta bins."};
+     Property<std::vector<float>> m_statPtBins
+        {this, "StatisticPtBins", {1.e3,2.5e3,10e3, 100e3}, "Gather statistics separately for these pt bins."};
+     Property<bool> m_pdgIdCategorisation
+        {this, "PdgIdCategorisation", false, "Categorise by pdg id."};
+     Property<bool> m_showRawCounts
+        {this, "ShowRawCounts", false, "Show all counters."};
+     Property<bool> m_printDetails
+        {this, "ShowDetailedTables", false, "Show more details; stat. uncert., RMS, entries"};
+     Property<bool> m_computeTrackRecoEfficiency
+        {this, "ComputeTrackRecoEfficiency", true, "Compute and print track reconstruction efficiency."};
+
+    ToolHandle<IAthSelectionTool> m_truthSelectionTool{this, "TruthSelectionTool","AthTruthSelectionTool", "Truth selection tool (for efficiencies and resolutions)"};
+
+
+     // helper struct for compile time optional statistics
+     template <bool IsDebug>
+     struct DebugCounter {
+        struct Empty {
+           template <typename... T_Args>
+           Empty(T_Args... ) {}
+        };
+        mutable typename std::conditional<IsDebug,
+                                          std::mutex,
+                                          Empty>::type m_mutex ATLAS_THREAD_SAFE;
+        mutable typename std::conditional<IsDebug,
+                                          ActsUtils::StatHist,
+                                          Empty>::type m_measPerTruthParticleWithoutCounts ATLAS_THREAD_SAFE {20,-.5,20.-.5};
+        mutable typename std::conditional<IsDebug,
+                                          ActsUtils::StatHist,
+                                          Empty>::type m_bestMatchProb ATLAS_THREAD_SAFE {20,0.,1.};
+        mutable typename std::conditional<IsDebug,
+                                          ActsUtils::StatHist,
+                                          Empty>::type m_nextToBestMatchProb ATLAS_THREAD_SAFE {20,0.,1.};
+
+        template <class T_OutStream>
+        void dumpStatistics(T_OutStream &out) const;
+        void fillMeasForTruthParticleWithoutCount(double weighted_measurement_sum) const;
+        void fillTruthMatchProb(const std::array<float,2> &best_match_prob) const;
+     };
+     DebugCounter<TrackFindingValidationDebugHists> m_debugCounter;
+
+     // s_NMeasurementTypes is equal to the number of UncalibMeasType, but we have to remove the "Other" option, which
+     // corresponds to an unknown type
+     constexpr static unsigned int s_NMeasurementTypes = static_cast<unsigned int>(xAOD::UncalibMeasType::nTypes) - 1u;
+
+     // statistics counter
+     enum ECounter {
+        NTracksTotal,
+        NTruthWithCountsTotal,
+        MissingTruthParticleHitCounts,
+        NoAssociatedTruthParticle,
+        NoSelectedTruthParticle,
+        TruthParticleNoNoiseMismatch,
+        kNCounter
+     };
+
+     enum ECategorisedCounter {
+        kNTotalParticles,
+        kNParticleWithAssociatedTrack,
+        kNParticleWithMultipleAssociatedTracks,
+        kNTotalTracks,
+        kNCategorisedCounter
+     };
+     enum ECategorisedStat {
+        kHitEfficiency,
+        kHitPurity,
+        kMatchProbability,
+        kNCategorisedStat
+     };
+
+     constexpr static int s_pdgIdMax = 1000000000;  // categorise all truth particles with this or larger PDG ID as "Other"
+
+
+     bool m_useAbsEtaForStat = false;
+
+     template <bool DetailEnabled>
+     struct BaseStat {
+        BaseStat() = default;
+        BaseStat([[maybe_unused]] const IAthSelectionTool &truth_selection_tool,
+                    [[maybe_unused]] std::size_t per_eta_size,
+                    [[maybe_unused]] std::size_t per_pdg_size)
+           : m_truthSelectionCuts(truth_selection_tool.nCuts()+1, -0.5,truth_selection_tool.nCuts()+.5)
+        {
+           if constexpr(DetailEnabled) {
+              m_counterPerEta.resize(per_eta_size);
+              m_counterPerPdgId.resize( per_pdg_size);
+              m_statPerEta.resize( per_eta_size );
+              m_statPerPdgId.resize( per_pdg_size );
+           }
+
+        }
+        void reset(const IAthSelectionTool &truth_selection_tool,
+                   [[maybe_unused]] std::size_t per_eta_size,
+                   [[maybe_unused]] std::size_t per_pdg_size)
+        {
+           m_truthSelectionCuts.setBinning(truth_selection_tool.nCuts()+1, -0.5,truth_selection_tool.nCuts()+.5);
+           if constexpr(DetailEnabled) {
+              m_counterPerEta.clear();
+              m_counterPerPdgId.clear();
+              m_statPerEta.clear();
+              m_statPerPdgId.clear();
+              m_counterPerEta.resize(per_eta_size);
+              m_counterPerPdgId.resize( per_pdg_size);
+              m_statPerEta.resize( per_eta_size );
+              m_statPerPdgId.resize( per_pdg_size );
+           }
+        }
+        void fill([[maybe_unused]] unsigned int eta_category_i,
+                  [[maybe_unused]] unsigned int pdg_id_category_i,
+                  [[maybe_unused]] float hit_efficiency,
+                  [[maybe_unused]] float hit_purity,
+                  [[maybe_unused]] float match_prob) {
+           if (DetailEnabled) {
+              assert( eta_category_i <m_statPerEta.size());
+              m_statPerEta[eta_category_i][kHitEfficiency].add( hit_efficiency);
+              m_statPerEta[eta_category_i][kHitPurity].add( hit_purity);
+              m_statPerEta[eta_category_i][kMatchProbability].add( match_prob);
+              assert( pdg_id_category_i <m_statPerPdgId.size());
+              m_statPerPdgId[pdg_id_category_i][kHitEfficiency].add( hit_efficiency);
+              m_statPerPdgId[pdg_id_category_i][kHitPurity].add( hit_purity);
+              m_statPerPdgId[pdg_id_category_i][kMatchProbability].add( match_prob);
+              assert( eta_category_i < m_counterPerEta.size());
+              assert( pdg_id_category_i <m_counterPerPdgId.size());
+              ++m_counterPerEta[eta_category_i][kNTotalTracks];
+              ++m_counterPerPdgId[pdg_id_category_i][kNTotalTracks];
+           }
+        }
+        void incrementTotal([[maybe_unused]] unsigned int eta_category_i,
+                            [[maybe_unused]] unsigned int pdg_id_category_i) {
+           if constexpr(DetailEnabled) {
+              ++m_counterPerEta[eta_category_i][kNTotalParticles];
+              ++m_counterPerPdgId[pdg_id_category_i][kNTotalParticles];
+           }
+        }
+
+        BaseStat<DetailEnabled> &operator+=(const BaseStat<DetailEnabled> &event_stat);
+        
+        void printStatTables(const TrackTruthMatchingBaseAlg &parent,
+                             const std::vector<float> &statPtBins,
+                             const std::vector<float> &statEtaBins,
+                             std::vector< int > &pdgId,
+                             bool printDetails,
+                             bool pdgIdCategorisation,
+                             bool useAbsEtaForStat);
+        
+
+        std::size_t perEtaSize() const  {
+           if constexpr(DetailEnabled) { return m_counterPerEta.size(); }
+           else { return 0u; }
+        }
+        std::size_t perPdgIdSize() const {
+           if constexpr(DetailEnabled) { return m_counterPerPdgId.size(); }
+           else { return 0u; }
+        }
+        ActsUtils::StatHist                                            m_truthSelectionCuts;
+        // per event statistics
+        struct Empty {};
+        using CounterArrayVec = std::conditional< DetailEnabled,
+                                                  std::vector< std::array< std::size_t, kNCategorisedCounter> >,
+                                                  Empty >::type;
+        using StatArrayVec = std::conditional< DetailEnabled,
+                                               std::vector< std::array<ActsUtils::Stat, kNCategorisedStat> >,
+                                               Empty >::type;
+        CounterArrayVec m_counterPerEta;
+        CounterArrayVec m_counterPerPdgId;
+        StatArrayVec m_statPerEta;
+        StatArrayVec m_statPerPdgId;
+        
+     };
+     
+
+     mutable std::mutex                           m_statMutex ATLAS_THREAD_SAFE;
+     mutable std::array< std::size_t, kNCounter > m_counter ATLAS_THREAD_SAFE {};
+     mutable std::vector< int >                   m_pdgId ATLAS_THREAD_SAFE;
+     mutable BaseStat<TrackFindingValidationDetailedStat> m_detailedStat ATLAS_THREAD_SAFE;
+
+     /// @brief check that bins are in increasing order.
+     /// will cause a FATAL error if bins are not in increasing order.
+     void checkBinOrder( const std::vector<float> &bin_edges, const std::string &bin_label) const;
+
+     /// @brief Return the category based on the provided eta value
+     /// @param pt the pt of the truth particle
+     /// @param eta the eta of the truth particle
+     /// @return a bin assigned to the give eta value
+     std::size_t getPtEtaStatCategory(float pt, float eta) const;
+
+     /// @brief Return the category based on the PDG ID
+     /// @param pt the pt of the truth particle
+     /// @param pdg_id the PDG ID
+     /// @return 0 or a slot associated to a single PDG ID (absolut value)
+     std::size_t getPtPdgIdStatCategory(float pt, int pdg_id) const;
+     void initStatTables();
+     void printStatTables() const;
+     void printCategories(const std::vector<std::string> &pt_category_labels,
+                          const std::vector<std::string> &eta_category_labels,
+                          std::vector<std::string> &counter_labels,
+                          std::vector< std::array< ActsUtils::Stat, kNCategorisedStat> > &stat_per_category,
+                          std::vector< std::array< std::size_t, kNCategorisedCounter> >  &counts_per_category,
+                          const std::string &top_left_label,
+                          bool print_sub_categories) const;
+     void printData2D(const std::vector<std::string> &row_category_labels,
+                      const std::vector<std::string> &col_category_labels,
+                      const std::string &top_left_label,
+                      std::vector< std::array< ActsUtils::Stat, kNCategorisedStat> > &stat_per_category,
+                      std::vector< std::array< std::size_t, kNCategorisedCounter> >  &counts_per_category,
+                      bool rotate) const;
+
+     StatusCode checkMatchWeights();
+
+     static  double weightedCountSum(const ActsTrk::HitCounterArray &counts,
+                                     const std::vector<float> &weights);
+     static  double noiseCorrection(const ActsTrk::HitCounterArray &noise_counts,
+                                    const std::vector<float> &weights);
+
+  };
+
+} // namespace
+
+#endif

Tracking/Acts/ActsTruth/src/components/ActsTruth_entries.cxx

@@ -5,6 +5,7 @@
 #include "../TrackToTruthAssociationAlg.h"
 #include "../TruthParticleHitCountAlg.h"
 #include "../TrackFindingValidationAlg.h"
+#include "../TrackParticleTruthDecorationAlg.h"
 
 // Algorithms
 DECLARE_COMPONENT( ActsTrk::PixelClusterToTruthAssociationAlg )
@@ -12,3 +13,4 @@ DECLARE_COMPONENT( ActsTrk::StripClusterToTruthAssociationAlg )
 DECLARE_COMPONENT( ActsTrk::TrackToTruthAssociationAlg )
 DECLARE_COMPONENT( ActsTrk::TruthParticleHitCountAlg )
 DECLARE_COMPONENT( ActsTrk::TrackFindingValidationAlg )
+DECLARE_COMPONENT( ActsTrk::TrackParticleTruthDecorationAlg )

Tracking/Acts/ActsTruth/src/decoratorUtils.h

+/*
+  Copyright (C) 2002-2024 CERN for the benefit of the ATLAS collaboration
+*/
+
+#ifndef ACTSTRACK_SAFEDECORATOR_H
+#define ACTSTRACK_SAFEDECORATOR_H
+/**
+ * copied from InDetPhysValMonitoring/src/safeDecorator.h
+ */
+
+
+
+#include "AthContainers/AuxElement.h"
+#include "StoreGate/WriteDecorHandleKey.h"
+#include "StoreGate/WriteDecorHandle.h"
+#include "StoreGate/ReadDecorHandleKey.h"
+#include "StoreGate/ReadDecorHandle.h"
+#include "GaudiKernel/EventContext.h"
+
+#include <iostream>
+#include <utility>
+#include <vector>
+
+#include <cstdlib>
+
+namespace ActsTrk {
+
+  // convenience method to create several decorator handles from keys
+  template <class T_Cont, class T>
+  std::vector<SG::WriteDecorHandle<T_Cont,T> >
+  createDecorators(  const std::vector<SG::WriteDecorHandleKey<T_Cont> > &keys,
+                     const EventContext &ctx) {
+    std::vector<SG::WriteDecorHandle<T_Cont,T> > out;
+    out.reserve(keys.size());
+    for( const SG::WriteDecorHandleKey<T_Cont> &a_key : keys) {
+      out.emplace_back(a_key,ctx);
+      if (not out.back().isValid()) {
+        std::stringstream msg;
+        msg << "Failed to create decorator handdle " << a_key.key();
+        throw std::runtime_error( msg.str() );
+      }
+    }
+    return out;
+  }
+
+
+  // convenience method to create several decorator handle keys.
+  template<class T_Parent, class T_Cont>
+  void createDecoratorKeys(T_Parent &parent,
+                           const SG::ReadHandleKey<T_Cont> &container_key,
+                           const std::string &prefix,
+                           const std::vector<std::string> &decor_names,
+                           std::vector<SG::WriteDecorHandleKey<T_Cont> > &decor_out) {
+    decor_out.clear();
+    decor_out.reserve(decor_names.size());
+    for (const std::string &a_decor_name : decor_names) {
+      assert( !a_decor_name.empty() );
+      decor_out.emplace_back(container_key.key()+"."+prefix+a_decor_name);
+      // need to declare handles, otherwise the scheduler would not pick up the data dependencies
+      // introduced by the decorations
+      parent.declare(decor_out.back());
+      decor_out.back().setOwner(&parent);
+      decor_out.back().initialize().ignore();
+    }
+  }
+}
+#endif