g4simulation/g4eval/TruthRecoTrackMatching.h

0001 #ifndef TRUTHTRKMATCHER__H
0002 #define TRUTHTRKMATCHER__H
0003 #include "TrackClusEvaluator.h"
0004
0005 #include <fun4all/SubsysReco.h>
0006
0007 #include <trackbase/ActsGeometry.h>
0008 #include <trackbase/TrkrDefs.h>
0009 #include <trackbase_historic/TrackSeed.h>
0010
0011 #include <array>
0012 #include <iostream>
0013 #include <map>
0014 #include <set>
0015 #include <string>
0016 #include <tuple>
0017 #include <vector>
0018
0019 class EmbRecoMatchContainer;
0020 class PHCompositeNode;
0021 class PHG4TpcCylinderGeom;
0022 class PHG4TpcCylinderGeomContainer;
0023 class PHG4TruthInfoContainer;
0024 class SvtxTrack;
0025 class SvtxTrackMap;
0026 class TrackSeedContainer;
0027 class TrkrCluster;
0028 class TrkrClusterContainer;
0029 class TrkrTruthTrack;
0030 class TrkrTruthTrackContainer;
0031 class TrkrClusterIsMatcher;
0032 class TTree;
0033 class TFile;
0034
0035 class TruthRecoTrackMatching : public SubsysReco
0036 {
0037   //--------------------------------------------------
0038   // Standard public interface
0039   //--------------------------------------------------
0040  public:
0041   TruthRecoTrackMatching( // Criteria to match a TrkrClusterContainer and track
0042       /*-----------------------------------------------------------------------------------
0043       * Input criteria for Truth Track (with nT clusters) to reco track (with nR clusters) :
0044       *  - nmin_match  : minimum number of clusters in Truth and Reco track that must
0045       *                  match for tracks to match
0046       *  - cutoff_dphi : maximum distance out in |phi_truth-phi_reco| for a matched track
0047       *  - same_dphi   : within this distance, all tracks must be checked
0048       *  - cutoff_deta :   like cutoff_dphi but for eta
0049       *  - same_deta   :   like same_dphi but for eta
0050       *  - cluster_nzwidths   : z-distance allowed for the truth center to be from
0051       *                         reco-cluster center:
0052       *                         |z_true-z_reco| must be <= dz_reco * cluster_nzwidths
0053       *  - cluster_nphiwidths : like cluster_nzwidths for phi
0054       *--------------------------------------------------------*/
0055         TrkrClusterIsMatcher* _ismatcher
0056       , const unsigned short _nmin_match = 4
0057       , const float _nmin_ratio          = 0.
0058       , const float _cutoff_dphi         = 0.3
0059       , const float _same_dphi           = 0.05
0060       , const float _cutoff_deta         = 0.3
0061       , const float _same_deta           = 0.05
0062       , const unsigned short _max_nreco_per_truth = 4
0063       , const unsigned short _max_ntruth_per_reco = 4);  // for some output kinematis
0064
0065   ~TruthRecoTrackMatching() override = default;
0066   int Init(PHCompositeNode*) override { return 0; };
0067   int InitRun(PHCompositeNode*) override;        //`
0068   int process_event(PHCompositeNode*) override;  //`
0069   int End(PHCompositeNode*) override;
0070
0071   TrackClusEvaluator    m_TCEval   ;
0072   TrkrClusterIsMatcher* m_ismatcher {nullptr};
0073
0074   int createNodes(PHCompositeNode* topNode);
0075
0076   void set_min_cl_match             (unsigned short val) { m_nmincluster_match = val; };
0077   void set_min_cl_ratio             (float val) { m_nmincluster_ratio  = val  ; };
0078   void set_cutoff_deta              (float val) { m_cutoff_deta        = val; };
0079   void set_cutoff_dphi              (float val) { m_cutoff_dphi        = val; };
0080   void set_nmin_truth_cluster_ratio (float val) { m_nmincluster_ratio  = val; };
0081   void set_smallsearch_deta         (float val) { m_same_deta          = val; };
0082   void set_smallsearch_dphi         (float val) { m_same_dphi          = val; };
0083
0084   void set_max_nreco_per_truth(unsigned short val) { m_max_nreco_per_truth = val; };
0085   void set_max_ntruth_per_reco(unsigned short val) { m_max_ntruth_per_reco = val; };
0086
0087  private:
0088   //--------------------------------------------------
0089   // Internal functions
0090   //--------------------------------------------------
0091
0092   //--------------------------------------------------
0093   // Constant parameters for track matching
0094   //--------------------------------------------------
0095
0096   unsigned short m_nmincluster_match;  // minimum of matched clustered to keep a truth to emb match
0097   float          m_nmincluster_ratio;           // minimum ratio of truth clustered that must be matched in reconstructed track
0098
0099   double m_cutoff_dphi;  // how far in |phi_truth-phi_reco| to match
0100   double m_same_dphi;    // |phi_truth-phi_reco| to auto-evaluate (is < _m_cutoff_dphi)
0101   double m_cutoff_deta;  //  how far in |eta_truth-eta_reco| to match
0102   double m_same_deta;    // |eta_truth-eta_reco| to auto-evaluate (is < m_cutoff_deta)
0103
0104   /* double m_cluster_nzwidths;   // cutoff in *getPhiSize() in cluster for |cluster_phi_truth-cluster_phi_reco| to match */
0105   /* double m_cluster_nphiwidths; // same for eta */
0106
0107   unsigned short m_max_nreco_per_truth;
0108   unsigned short m_max_ntruth_per_reco;
0109
0110   /* std::array<double, 55> m_phistep {0.}; // the phistep squared */
0111   /* double m_zstep {0.}; */
0112
0113   std::map<unsigned short, unsigned short> m_nmatched_index_true;
0114
0115   std::map<unsigned short, unsigned short>* m_nmatched_id_reco = nullptr;
0116   std::map<unsigned short, unsigned short>* m_nmatched_id_true = nullptr;
0117
0118   //--------------------------------------------------
0119   // Data from input nodes
0120   //--------------------------------------------------
0121   PHG4TruthInfoContainer* m_PHG4TruthInfoContainer = nullptr;  // Get the truth track ids
0122   SvtxTrackMap* m_SvtxTrackMap = nullptr;
0123   TrkrClusterContainer* m_TruthClusterContainer = nullptr;
0124   TrkrClusterContainer* m_RecoClusterContainer = nullptr;
0125   TrkrTruthTrackContainer* m_TrkrTruthTrackContainer = nullptr;
0126   PHG4TpcCylinderGeomContainer* m_PHG4TpcCylinderGeomContainer = nullptr;
0127
0128   ActsGeometry* m_ActsGeometry = nullptr;
0129
0130   // Output data node:
0131   EmbRecoMatchContainer* m_EmbRecoMatchContainer = nullptr;
0132
0133   //--------------------------------------------------
0134   //    RECO data for a "table" of reconstructed tracks
0135   //--------------------------------------------------
0136   using RECOentry = std::tuple<float, float, float, unsigned short>;
0137   static constexpr int RECOphi = 0;
0138   static constexpr int RECOeta = 1;
0139   static constexpr int RECOpt = 2;
0140   static constexpr int RECOid = 3;
0141
0142  public:
0143   using RECOvec = std::vector<RECOentry>;
0144   using RECOiter = RECOvec::iterator;
0145   using RECO_pair_iter = std::pair<RECOiter, RECOiter>;
0146
0147  private:
0148   struct CompRECOtoPhi
0149   {
0150     bool operator()(const RECOentry& lhs, const double& rhs) { return std::get<RECOphi>(lhs) < rhs; }
0151     bool operator()(const double& lhs, const RECOentry& rhs) { return lhs < std::get<RECOphi>(rhs); }
0152     bool operator()(const RECOentry& lhs, const RECOentry& rhs) { return std::get<RECOphi>(lhs) < std::get<RECOphi>(rhs); }
0153   };
0154   struct CompRECOtoEta
0155   {
0156     bool operator()(const RECOentry& lhs, const double& rhs) { return std::get<RECOeta>(lhs) < rhs; }
0157     bool operator()(const double& lhs, const RECOentry& rhs) { return lhs < std::get<RECOeta>(rhs); }
0158     bool operator()(const RECOentry& lhs, const RECOentry& rhs) { return std::get<RECOeta>(lhs) < std::get<RECOeta>(rhs); }
0159   };
0160   struct CompRECOtoPt
0161   {
0162     bool operator()(const RECOentry& lhs, const double& rhs) { return std::get<RECOpt>(lhs) < rhs; }
0163     bool operator()(const double& lhs, const RECOentry& rhs) { return lhs < std::get<RECOpt>(rhs); }
0164     bool operator()(const RECOentry& lhs, const RECOentry& rhs) { return std::get<RECOpt>(lhs) < std::get<RECOpt>(rhs); }
0165   };
0166
0167   // sorting tuple is by: phi, eta, pT,  index, is_matched
0168   //--------------------------------------------------
0169   //    PossibleMatches (just array<unsinged short, 5>)
0170   //--------------------------------------------------
0171  public:
0172   using PossibleMatch = std::array<unsigned short, 5>;
0173
0174  private:
0175   static constexpr int PM_nmatch = 0;
0176   static constexpr int PM_ntrue = 1;
0177   static constexpr int PM_nreco = 2;
0178   static constexpr int PM_idtrue = 3;
0179   static constexpr int PM_idreco = 4;
0180   struct SortPossibleMatch
0181   {
0182     // Sort by most matched clusters first, then smallest number of truth clusters, then smallest number of reco clusters
0183     bool operator()(const PossibleMatch& lhs, const PossibleMatch& rhs)
0184     {
0185       if (lhs[PM_nmatch] != rhs[PM_nmatch]) return lhs[PM_nmatch] > rhs[PM_nmatch];
0186       if (lhs[PM_ntrue] != rhs[PM_ntrue]) return lhs[PM_ntrue] < rhs[PM_ntrue];
0187       if (lhs[PM_nreco] != rhs[PM_nreco]) return lhs[PM_nreco] < rhs[PM_nreco];
0188       return false;
0189     }
0190   };
0191
0192   //--------------------------------------------------
0193   //    PossibleMatches (just array<unsinged short, 5>)
0194   //--------------------------------------------------
0195  public:
0196   //--------------------------------------------------
0197   // non-node member data
0198   //--------------------------------------------------
0199   RECOvec recoData{};  // "sv" = Sort Vector
0200
0201   //--------------------------------------------------
0202   // Member functions
0203   //--------------------------------------------------
0204
0205   //    Main functions
0206   // -------------------------------------------------------------------
0207   std::pair<std::vector<unsigned short>, std::vector<unsigned short>>
0208   find_box_matches(float truth_phi, float truth_eta, float truth_pt);                         // will populate to truth_to_reco_map and
0209   void match_tracks_in_box(std::vector<std::pair<unsigned short, unsigned short>>& indices);  // pairs of {id_true, id_reco}
0210
0211   //    Helper functions
0212   // -------------------------------------------------------------------
0213   float delta_outer_pt(float) const;
0214   float delta_inner_pt(float) const;
0215   float abs_dphi(float phi0, float phi1);
0216   float sigma_CompMatchClusters(PossibleMatch&);
0217   bool skip_match(PossibleMatch& match);
0218   bool at_nmax_index_true(unsigned short);  // test if the truth track already has maximum matches matches
0219   bool at_nmax_id_reco(unsigned short);     // "           reco                                          "
0220
0221   std::pair<bool, float> compare_cluster_pair(TrkrDefs::cluskey key_T,
0222                                               TrkrDefs::cluskey key_R, TrkrDefs::hitsetkey key, bool calc_sigma = false);
0223
0224   // ------------------------------------------------------------------
0225   // output for diagnositics:
0226   //  If there is output, put all the clusters into the x, y, z
0227   // ------------------------------------------------------------------
0228   TTree* m_diag_tree = nullptr;
0229   TFile* m_diag_file = nullptr;
0230   bool m_write_diag = false;
0231   void set_diagnostic_file(const std::string& file_name);
0232
0233   std::vector<int> m_trkid_reco_matched{};
0234   std::vector<int> m_cnt_reco_matched{};
0235   std::vector<unsigned int> m_i0_reco_matched{};
0236   std::vector<unsigned int> m_i1_reco_matched{};
0237   std::vector<unsigned int> m_layer_reco_matched{};
0238   std::vector<float> m_x_reco_matched{};
0239   std::vector<float> m_y_reco_matched{};
0240   std::vector<float> m_z_reco_matched{};
0241
0242   std::vector<int> m_trkid_reco_notmatched{};
0243   std::vector<int> m_cnt_reco_notmatched{};
0244   std::vector<unsigned int> m_i0_reco_notmatched{};
0245   std::vector<unsigned int> m_i1_reco_notmatched{};
0246   std::vector<unsigned int> m_layer_reco_notmatched{};
0247   std::vector<float> m_x_reco_notmatched{};
0248   std::vector<float> m_y_reco_notmatched{};
0249   std::vector<float> m_z_reco_notmatched{};
0250
0251   std::vector<int> m_trkid_true_matched{};
0252   std::vector<int> m_cnt_true_matched{};
0253   std::vector<unsigned int> m_i0_true_matched{};
0254   std::vector<unsigned int> m_i1_true_matched{};
0255   std::vector<unsigned int> m_layer_true_matched{};
0256   std::vector<float> m_x_true_matched{};
0257   std::vector<float> m_y_true_matched{};
0258   std::vector<float> m_z_true_matched{};
0259
0260   std::vector<int> m_trkid_true_notmatched{};
0261   std::vector<int> m_cnt_true_notmatched{};
0262   std::vector<unsigned int> m_i0_true_notmatched{};
0263   std::vector<unsigned int> m_i1_true_notmatched{};
0264   std::vector<unsigned int> m_layer_true_notmatched{};
0265   std::vector<float> m_x_true_notmatched{};
0266   std::vector<float> m_y_true_notmatched{};
0267   std::vector<float> m_z_true_notmatched{};
0268
0269   int m_event = 0;
0270
0271   void clear_branch_vectors();
0272   void fill_tree();
0273 };
0274
0275 #endif