g4simulation/g4eval/TruthRecoTrackMatching.h

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