TrackTruthMatchingBaseAlg.h

Go to the documentation of this file.

/*
  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 "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:
    using AthReentrantAlgorithm::AthReentrantAlgorithm;
 
    virtual StatusCode initialize() override;
    virtual StatusCode finalize() override;
 
  protected:
     const TruthParticleHitCounts &getTruthParticleHitCounts(const EventContext &ctx) const {
        const TruthParticleHitCounts* truth_particle_hit_counts_handle{};
       if (!SG::get(truth_particle_hit_counts_handle, m_truthHitCounts, ctx).isSuccess()) {
           ATH_MSG_ERROR("No truth particle hit count map for key " << m_truthHitCounts.key() );
           throw std::runtime_error("Failed to get truth particle hit count map");
        }
        return *truth_particle_hit_counts_handle;
     }
     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_t< DetailEnabled,
                                                   std::unordered_set<const xAOD::TruthParticle *>,
                                                   typename BaseStat<DetailEnabled>::Empty >;
        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{0u};
     };
     using EventStat = EventStatBase<TrackFindingValidationDetailedStat>;

     struct TruthMatchResult {
        const xAOD::TruthParticle *m_truthParticle{nullptr};  
        float m_matchProbability{0.f};                    
        float m_hitPurity{0.f};                           
        float m_hitEfficiency{0.f};                       
     };

     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;
  protected:
      ToolHandle<IAthSelectionTool> m_truthSelectionTool{this, "TruthSelectionTool","AthTruthSelectionTool", "Truth selection tool (for efficiencies and resolutions)"};
  private:
     SG::ReadHandleKey<TruthParticleHitCounts>  m_truthHitCounts
        {this, "TruthParticleHitCounts","", "Map from truth particle to hit counts." };
 
     Gaudi::Property<std::vector<float> > m_weightsForProb
        {this, "MatchWeights", std::vector<float>(s_NMeasurementTypes,1.f), 
         "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", std::vector<float>(s_NMeasurementTypes,1.f), 
         "Weights applied to the counts per measurement type for weighted sums"
         " which are used to compute hit efficiencies and purities." };
      Gaudi::Property<bool> m_truthLinkedFakeRate
        {this, "TruthLinkedFakeRate", false, "Compute the match probability from the truth-linked measurements only:"
                                             " measurements without truth particle weigh 1 and cannot be the best match, so a"
                                             " track made only of them has probability 0. By default the probability is"
                                             " computed the way the legacy truth matching does (DetailedTrackTruthBuilder,"
                                             " TruthMatchRatio): measurements without truth particle keep their full weight and"
                                             " together compete for the best match as one pseudo particle; if it wins, the"
                                             " probability is their weighted fraction and there is no matching truth particle." };
      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"};
     Gaudi::Property<bool> m_computeTrackRecoEfficiency
        {this, "ComputeTrackRecoEfficiency", true, "Compute and print track reconstruction efficiency."};
 
     // 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
     constexpr static unsigned int s_NMeasurementTypes = Acts::toUnderlying(xAOD::UncalibMeasType::nTypes);
 
     // 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;

     void checkBinOrder( const std::vector<float> &bin_edges, const std::string &bin_label) const;

     std::size_t getPtEtaStatCategory(float pt, float eta) const;

     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