macro/loc_lib/JetIndicesMatcher.cc

0001 #include "JetIndicesMatcher.h"
0002 #include "cmath"
0003 using namespace std;
0004
0005 JetIndicesMatcher::JetIndicesMatcher(float R, float min_T, float min_R)
0006   : R2 {R*R}, min_pT_truth{min_T}, min_pT_reco{min_R}
0007 {}
0008
0009 void JetIndicesMatcher::reset()
0010 {
0011     eta_truth.clear();
0012     phi_truth.clear();
0013     pt_truth.clear();
0014     index_truth.clear();
0015
0016     eta_reco.clear();
0017     phi_reco.clear();
0018     pt_reco.clear();
0019     index_reco.clear();
0020
0021     indices_fake.clear();
0022     indices_miss.clear();
0023     indices_matched.clear();
0024 }
0025
0026 /* void JetIndicesMatcher::add_truth(float eta, float phi, float pt) { */
0027 /*     /1* if (fabs(eta)>0.6) return; *1/ */
0028 /*     eta_truth.push_back(eta); */
0029 /*     phi_truth.push_back(phi); */
0030 /*     pt_truth.push_back(pt); */
0031 /* } */
0032
0033 /* void JetIndicesMatcher::add_reco(float eta, float phi, float pt) { */
0034 /*     /1* if (fabs(eta)>0.6) return; *1/ */
0035 /*     eta_reco.push_back(eta); */
0036 /*     phi_reco.push_back(phi); */
0037 /*     pt_reco.push_back(pt); */
0038 /* } */
0039
0040 void JetIndicesMatcher::add_truth(vector<float>& eta, vector<float>& phi, vector<float>& pt) {
0041     vector<array<float,4>> jets;
0042     for (unsigned int i=0; i<eta.size(); ++i) {
0043         jets.push_back({pt[i],eta[i],phi[i], (float) i});
0044     }
0045     std::sort(jets.rbegin(), jets.rend());
0046     for (auto jet : jets) {
0047       pt_truth  .push_back(jet[0]);
0048       eta_truth .push_back(jet[1]);
0049       phi_truth .push_back(jet[2]);
0050       index_truth .push_back((int)jet[3]);
0051     }
0052 }
0053 void JetIndicesMatcher::add_reco(vector<float>& eta, vector<float>& phi, vector<float>& pt) {
0054     vector<array<float,4>> jets;
0055     for (unsigned int i=0; i<eta.size(); ++i) {
0056         jets.push_back({pt[i],eta[i],phi[i], (float)i});
0057     }
0058     std::sort(jets.rbegin(), jets.rend());
0059     for (auto jet : jets) {
0060       pt_reco  .push_back(jet[0]);
0061       eta_reco .push_back(jet[1]);
0062       phi_reco .push_back(jet[2]);
0063       index_reco .push_back((int)jet[3]);
0064     }
0065 }
0066 array<unsigned int,3> JetIndicesMatcher::do_matching() {
0067     indices_fake.clear();
0068     indices_miss.clear();
0069     indices_matched.clear();
0070
0071     vector<bool> is_matched (eta_reco.size(),false);
0072
0073     //FIXME
0074           /* cout << " Matching Truth Jets " << endl; */
0075           /* for (int i = 0; i< pt_truth.size(); ++i) { */
0076           /*   cout << Form("  Tjet [%2i] pt:phi:eta [%5.2f,%5.2f,%5.2f]", */
0077           /*       i, pt_truth[i],phi_truth[i],eta_truth[i]) <<endl; */
0078           /* } */
0079           /* cout << " In matcher : Reco Jets " << endl; */
0080           /* for (int i = 0; i< pt_reco.size(); ++i) { */
0081           /*   cout << Form("  Mjet [%2i] pt:phi:eta [%5.2f,%5.2f,%5.2f]", */
0082           /*       i, pt_reco[i],phi_reco[i],eta_reco[i]) <<endl; */
0083           /* } */
0084
0085
0086     for (unsigned int T=0;T<eta_truth.size();++T) {
0087       if (pt_truth[T]<min_pT_truth) continue;
0088         /* if (fabs(eta_truth[T])>0.6) continue; */
0089         bool found_match { false };
0090         for (unsigned int R=0;R<eta_reco.size();++R) {
0091           if (pt_reco[R]<min_pT_reco) continue;
0092             /* if (fabs(eta_reco[R])>0.6) continue; */
0093             if (is_matched[R]) continue;
0094             float dphi = fabs(phi_truth[T]-phi_reco[R]);
0095             while (dphi>M_PI) dphi = fabs(dphi - 2*M_PI);
0096             const float deta = eta_truth[T]-eta_reco[R];
0097             const float R2_comp = deta*deta + dphi*dphi;
0098             if (R2_comp<=R2) {
0099                 found_match = true;
0100                 is_matched[R] = true;
0101                 indices_matched.push_back({index_truth[T],index_reco[R]});
0102
0103                 /* cout << " ADDED MATCH( " << T << ":"<<R<<")" << endl; */
0104                 break;
0105             }
0106         }
0107         if (!found_match) indices_miss.push_back(index_truth[T]);
0108     }
0109     for (unsigned int R=0;R<eta_reco.size();++R) {
0110       /* cout << " min_pT_reco: " << min_pT_reco << " and["<<R<<"] " << pt_reco[R] << endl; */
0111         if (pt_reco[R]<min_pT_reco) continue;
0112         if (!is_matched[R]) {
0113           indices_fake.push_back(index_reco[R]);
0114           /* cout << " indices fake: "; for (auto I : indices_fake) cout << " " << I; cout << endl; */
0115         }
0116     }
0117     return {static_cast<unsigned int>(indices_matched.size()),
0118             static_cast<unsigned int>(indices_miss.size()),
0119             static_cast<unsigned int>(indices_fake.size())};
0120 }