macro/loc_lib/JetIndicesMatcher_bkup.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
0015     eta_reco.clear();
0016     phi_reco.clear();
0017     pt_reco.clear();
0018
0019     indices_fake.clear();
0020     indices_miss.clear();
0021     indices_matched.clear();
0022 }
0023
0024 /* void JetIndicesMatcher::add_truth(float eta, float phi, float pt) { */
0025 /*     /1* if (fabs(eta)>0.6) return; *1/ */
0026 /*     eta_truth.push_back(eta); */
0027 /*     phi_truth.push_back(phi); */
0028 /*     pt_truth.push_back(pt); */
0029 /* } */
0030
0031 /* void JetIndicesMatcher::add_reco(float eta, float phi, float pt) { */
0032 /*     /1* if (fabs(eta)>0.6) return; *1/ */
0033 /*     eta_reco.push_back(eta); */
0034 /*     phi_reco.push_back(phi); */
0035 /*     pt_reco.push_back(pt); */
0036 /* } */
0037
0038 void JetIndicesMatcher::add_truth(vector<float>& eta, vector<float>& phi, vector<float>& pt) {
0039     vector<array<float,3>> jets;
0040     for (unsigned int i=0; i<eta.size(); ++i) {
0041         jets.push_back({pt[i],eta[i],phi[i]});
0042     }
0043     std::sort(jets.rbegin(), jets.rend());
0044     for (auto jet : jets) {
0045       pt_truth  .push_back(jet[0]);
0046       eta_truth .push_back(jet[1]);
0047       phi_truth .push_back(jet[2]);
0048     }
0049 }
0050 void JetIndicesMatcher::add_reco(vector<float>& eta, vector<float>& phi, vector<float>& pt) {
0051     vector<array<float,3>> jets;
0052     for (unsigned int i=0; i<eta.size(); ++i) {
0053         jets.push_back({pt[i],eta[i],phi[i]});
0054     }
0055     /* cout << " First: " << endl; */
0056     /* for (auto jet : jets) cout << jet[0] <<":"<<jet[1]<<":"<<jet[2] << endl; */
0057     std::sort(jets.rbegin(), jets.rend());
0058     /* cout << " second: " << endl; */
0059     for (auto jet : jets) {
0060       /* cout << jet[0] <<":"<<jet[1]<<":"<<jet[2] << endl; */
0061       pt_reco  .push_back(jet[0]);
0062       eta_reco .push_back(jet[1]);
0063       phi_reco .push_back(jet[2]);
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({T,R});
0102
0103                 /* cout << " ADDED MATCH( " << T << ":"<<R<<")" << endl; */
0104                 break;
0105             }
0106         }
0107         if (!found_match) indices_miss.push_back(T);
0108     }
0109     for (unsigned int R=0;R<eta_reco.size();++R) {
0110         if (pt_reco[R]<min_pT_reco) continue;
0111         if (!is_matched[R]) indices_fake.push_back(R);
0112     }
0113     return {static_cast<unsigned int>(indices_matched.size()),
0114             static_cast<unsigned int>(indices_miss.size()),
0115             static_cast<unsigned int>(indices_fake.size())};
0116 }