JetBkgdSub/offline/MatchJets.C

0001 #include <iostream>
0002 #include <vector>
0003
0004 using namespace std;
0005
0006
0007 Int_t MatchJets(TString infile="../macro/output.root", TString outfile="matched.root", Float_t JetParameter=0.4, Float_t pTHard=30.0)
0008 {
0009
0010     // Set event selection cuts
0011     double dRMax = 0.75*JetParameter;
0012     float pTmin_truth = 10;
0013     float pTmin_reco = 5;
0014     float pThard_min, pThard_max;
0015     float weight;
0016     if(pTHard == 10){
0017         pThard_min = 10;
0018         pThard_max = 30;
0019         weight = 3.210e-6;
0020     }
0021     else{
0022         pThard_min = 30;
0023         pThard_max = 999;
0024         weight = 2.178e-9;
0025     }
0026
0027
0028     //=======================================================
0029     //================ Input File ===========================
0030     //=======================================================
0031     std::cout << "Input file: " << infile.Data() << std::endl;
0032
0033     // open input file
0034     TFile fin(infile.Data(), "READ");
0035     if(!fin.IsOpen())
0036     {
0037         std::cout << "Error: could not open input file" << std::endl;
0038         exit(-1);
0039     }
0040
0041     // get tree
0042     TTree *tree = (TTree*)fin.Get("T");
0043     if(!tree)
0044     {
0045         std::cout << "Error: could not find tree" << std::endl;
0046         exit(-1);
0047     }
0048
0049     int event = 0;
0050     int centrality = 0;
0051     float event_leading_truth_pt = 0;
0052     double rho_area = 0;
0053     double rho_area_sigma = 0;
0054     double rho_mult = 0;
0055
0056     // reco jet variables
0057     std::vector<float> * iter_eta = 0;
0058     std::vector<float> * iter_phi = 0;
0059     std::vector<float> * iter_pt = 0;
0060     std::vector<float> * iter_pt_unsub = 0;
0061
0062     // mult jet variables
0063     std::vector<int> * mult_ncomponent = 0;
0064     std::vector<float> * mult_eta = 0;
0065     std::vector<float> * mult_phi = 0;
0066     std::vector<float> * mult_pt = 0;
0067     std::vector<float> * mult_pt_unsub = 0;
0068     std::vector<int> * mult_nsignal = 0;
0069
0070     // truth jet variables
0071     std::vector<int> * truth_ncomponent = 0;
0072     std::vector<float> * truth_eta = 0;
0073     std::vector<float> * truth_phi = 0;
0074     std::vector<float> * truth_pt = 0;
0075
0076     //rhoA jet variables
0077     std::vector<float> * rhoA_eta = 0;
0078     std::vector<float> * rhoA_phi = 0;
0079     std::vector<float> * rhoA_pt = 0;
0080     std::vector<float> * rhoA_pt_unsub = 0;
0081     std::vector<float> * rhoA_area = 0;
0082
0083     tree->SetBranchAddress("event", &event);
0084     tree->SetBranchAddress("event_leading_truth_pt", &event_leading_truth_pt);
0085     tree->SetBranchAddress("centrality", &centrality);
0086     tree->SetBranchAddress("truth_ncomponent", &truth_ncomponent);
0087     tree->SetBranchAddress("truth_eta", &truth_eta);
0088     tree->SetBranchAddress("truth_phi", &truth_phi);
0089     tree->SetBranchAddress("truth_pt", &truth_pt);
0090     tree->SetBranchAddress("rho_area", &rho_area);
0091     tree->SetBranchAddress("rho_area_sigma", &rho_area_sigma);
0092     tree->SetBranchAddress("rhoA_eta", &rhoA_eta);
0093     tree->SetBranchAddress("rhoA_phi", &rhoA_phi);
0094     tree->SetBranchAddress("rhoA_pt", &rhoA_pt);
0095     tree->SetBranchAddress("rhoA_area", &rhoA_area);
0096     tree->SetBranchAddress("rhoA_pt_unsub", &rhoA_pt_unsub);
0097     tree->SetBranchAddress("rho_mult", &rho_mult);
0098     tree->SetBranchAddress("mult_ncomponent", &mult_ncomponent);
0099     tree->SetBranchAddress("mult_nsignal", &mult_nsignal);
0100     tree->SetBranchAddress("mult_eta", &mult_eta);
0101     tree->SetBranchAddress("mult_phi", &mult_phi);
0102     tree->SetBranchAddress("mult_pt", &mult_pt);
0103     tree->SetBranchAddress("mult_pt_unsub", &mult_pt_unsub);
0104     tree->SetBranchAddress("iter_eta", &iter_eta);
0105     tree->SetBranchAddress("iter_phi", &iter_phi);
0106     tree->SetBranchAddress("iter_pt", &iter_pt);
0107     tree->SetBranchAddress("iter_pt_unsub", &iter_pt_unsub);
0108
0109
0110     //=======================================================
0111
0112     // output file
0113     std::cout << "Output file: " << outfile.Data() << std::endl;
0114     TFile *fout = new TFile(outfile.Data(), "RECREATE");
0115     TTree *out_tree = new TTree("T", "MatchedJets");
0116
0117     // output variables
0118     int out_centrality = 0;
0119     double out_rho_area = 0;
0120     double out_rho_mult = 0;
0121
0122     std::vector<float> matched_pt_iter, matched_pt_area, matched_pt_mult;
0123     std::vector<float> matched_pt_iter_unsub, matched_pt_area_unsub, matched_pt_mult_unsub;
0124     std::vector<float> matched_truth_pt_iter, matched_truth_pt_area, matched_truth_pt_mult;
0125     std::vector<float> unmatched_pt_iter, unmatched_pt_area, unmatched_pt_mult;
0126     std::vector<float> unmatched_pt_iter_unsub, unmatched_pt_area_unsub, unmatched_pt_mult_unsub;
0127     std::vector<float> missed_truth_pt_iter, missed_truth_pt_area, missed_truth_pt_mult;
0128
0129     // set branch addresses
0130     out_tree->Branch("weight", &weight);
0131     out_tree->Branch("centrality", &out_centrality);
0132     out_tree->Branch("rho_area", &out_rho_area);
0133     out_tree->Branch("rho_mult", &out_rho_mult);
0134     out_tree->Branch("matched_pt_iter", &matched_pt_iter);
0135     out_tree->Branch("matched_pt_area", &matched_pt_area);
0136     out_tree->Branch("matched_pt_mult", &matched_pt_mult);
0137     out_tree->Branch("matched_pt_iter_unsub", &matched_pt_iter_unsub);
0138     out_tree->Branch("matched_pt_area_unsub", &matched_pt_area_unsub);
0139     out_tree->Branch("matched_pt_mult_unsub", &matched_pt_mult_unsub);
0140     out_tree->Branch("matched_truth_pt_iter", &matched_truth_pt_iter);
0141     out_tree->Branch("matched_truth_pt_area", &matched_truth_pt_area);
0142     out_tree->Branch("matched_truth_pt_mult", &matched_truth_pt_mult);
0143     out_tree->Branch("unmatched_pt_iter", &unmatched_pt_iter);
0144     out_tree->Branch("unmatched_pt_area", &unmatched_pt_area);
0145     out_tree->Branch("unmatched_pt_mult", &unmatched_pt_mult);
0146     out_tree->Branch("unmatched_pt_iter_unsub", &unmatched_pt_iter_unsub);
0147     out_tree->Branch("unmatched_pt_area_unsub", &unmatched_pt_area_unsub);
0148     out_tree->Branch("unmatched_pt_mult_unsub", &unmatched_pt_mult_unsub);
0149     out_tree->Branch("missed_truth_pt_iter", &missed_truth_pt_iter);
0150     out_tree->Branch("missed_truth_pt_area", &missed_truth_pt_area);
0151     out_tree->Branch("missed_truth_pt_mult", &missed_truth_pt_mult);
0152
0153     // get number of entries
0154     int nEntries = tree->GetEntries();
0155     for (int iEvent = 0; iEvent < nEntries; iEvent++){
0156         // get entry
0157         tree->GetEntry(iEvent);
0158
0159         // event selection
0160         if(event_leading_truth_pt < pThard_min || event_leading_truth_pt > pThard_max) continue;
0161
0162         // clear output variables
0163         out_centrality = centrality;
0164         out_rho_area = rho_area;
0165         out_rho_mult = rho_mult;
0166         // clear vectors
0167         matched_pt_iter.clear();
0168         matched_pt_area.clear();
0169         matched_pt_mult.clear();
0170         matched_pt_iter_unsub.clear();
0171         matched_pt_area_unsub.clear();
0172         matched_pt_mult_unsub.clear();
0173         matched_truth_pt_iter.clear();
0174         matched_truth_pt_area.clear();
0175         matched_truth_pt_mult.clear();
0176         unmatched_pt_iter.clear();
0177         unmatched_pt_area.clear();
0178         unmatched_pt_mult.clear();
0179         unmatched_pt_iter_unsub.clear();
0180         unmatched_pt_area_unsub.clear();
0181         unmatched_pt_mult_unsub.clear();
0182         missed_truth_pt_iter.clear();
0183         missed_truth_pt_area.clear();
0184         missed_truth_pt_mult.clear();
0185
0186         // vector of matched reco jets for each method
0187         std::vector<unsigned int> matched_area_jets, matched_mult_jets, matched_iter_jets;
0188
0189         // loop over truth jets
0190         for(unsigned int ijet =0; ijet < truth_pt->size(); ijet++)
0191         {
0192             // apply truth cut
0193             if(truth_pt->at(ijet) < pTmin_truth) continue;
0194
0195             float dr = 999;
0196             float matched_area_pt, matched_area_unsubpt;
0197             float matched_mult_pt, matched_mult_unsubpt;
0198             float matched_iter_pt, matched_iter_unsubpt;
0199             unsigned int matched_idx_area, matched_idx_mult, matched_idx_iter;
0200
0201             // match area jets to truth jets
0202             for(unsigned int jjet = 0; jjet < rhoA_pt->size(); jjet++)
0203             {
0204                 // apply reco jet cut
0205                 if(rhoA_pt->at(jjet) < pTmin_reco) continue;
0206
0207                 // dr calculation
0208                 float dphi = rhoA_phi->at(jjet) - truth_phi->at(ijet);
0209                 float deta = rhoA_eta->at(jjet) - truth_eta->at(ijet);
0210                 if(dphi > TMath::Pi()) dphi -= 2*TMath::Pi();
0211                 if(dphi < -TMath::Pi()) dphi += 2*TMath::Pi();
0212                 float dr_tmp = TMath::Sqrt(dphi*dphi + deta*deta);
0213                 if(dr_tmp < dr && dr_tmp < dRMax)
0214                 {
0215                     dr = dr_tmp;
0216                     matched_area_pt = rhoA_pt->at(jjet);
0217                     matched_area_unsubpt = rhoA_pt_unsub->at(jjet);
0218                     matched_idx_area = jjet;
0219                 }
0220             }
0221             if(dr > dRMax){
0222                 missed_truth_pt_area.push_back(truth_pt->at(ijet));
0223             }
0224             else{
0225                 matched_truth_pt_area.push_back(truth_pt->at(ijet));
0226                 matched_pt_area.push_back(matched_area_pt);
0227                 matched_pt_area_unsub.push_back(matched_area_unsubpt);
0228                 matched_area_jets.push_back(matched_idx_area);
0229             }
0230
0231             // reset dr
0232             dr = 999;
0233             // match mult jets to truth jets
0234             for(unsigned int jjet = 0; jjet < mult_pt->size(); jjet++)
0235             {
0236                 // apply reco jet cut
0237                 if(mult_pt->at(jjet) < pTmin_reco) continue;
0238
0239                 // dr calculation
0240                 float dphi = mult_phi->at(jjet) - truth_phi->at(ijet);
0241                 float deta = mult_eta->at(jjet) - truth_eta->at(ijet);
0242                 if(dphi > TMath::Pi()) dphi -= 2*TMath::Pi();
0243                 if(dphi < -TMath::Pi()) dphi += 2*TMath::Pi();
0244                 float dr_tmp = TMath::Sqrt(dphi*dphi + deta*deta);
0245                 if(dr_tmp < dr && dr_tmp < dRMax)
0246                 {
0247                     dr = dr_tmp;
0248                     matched_mult_pt = mult_pt->at(jjet);
0249                     matched_mult_unsubpt = mult_pt_unsub->at(jjet);
0250                     matched_idx_mult = jjet;
0251                 }
0252             }
0253             if (dr > dRMax){
0254                 missed_truth_pt_mult.push_back(truth_pt->at(ijet));
0255             }
0256             else{
0257                 matched_truth_pt_mult.push_back(truth_pt->at(ijet));
0258                 matched_pt_mult.push_back(matched_mult_pt);
0259                 matched_pt_mult_unsub.push_back(matched_mult_unsubpt);
0260                 matched_mult_jets.push_back(matched_idx_mult);
0261             }
0262
0263             // // reset dr
0264             dr = 999;
0265             // match iter jets to truth jets
0266             for (unsigned int jjet = 0; jjet < iter_pt->size(); jjet++)
0267             {
0268                 // apply reco jet cut
0269                 if(iter_pt->at(jjet) < pTmin_reco) continue;
0270
0271                 // dr calculation
0272                 float dphi = iter_phi->at(jjet) - truth_phi->at(ijet);
0273                 float deta = iter_eta->at(jjet) - truth_eta->at(ijet);
0274                 if(dphi > TMath::Pi()) dphi -= 2*TMath::Pi();
0275                 if(dphi < -TMath::Pi()) dphi += 2*TMath::Pi();
0276                 float dr_tmp = TMath::Sqrt(dphi*dphi + deta*deta);
0277                 if(dr_tmp < dr && dr_tmp < dRMax)
0278                 {
0279                     dr = dr_tmp;
0280                     matched_iter_pt = iter_pt->at(jjet);
0281                     matched_iter_unsubpt = iter_pt_unsub->at(jjet);
0282                     matched_idx_iter = jjet;
0283                 }
0284             }
0285             if (dr > dRMax){
0286                 missed_truth_pt_iter.push_back(truth_pt->at(ijet));
0287             }
0288             else{
0289                 matched_truth_pt_iter.push_back(truth_pt->at(ijet));
0290                 matched_pt_iter.push_back(matched_iter_pt);
0291                 matched_pt_iter_unsub.push_back(matched_iter_unsubpt);
0292                 matched_iter_jets.push_back(matched_idx_iter);
0293             }
0294
0295
0296         }
0297         // fill unmatched  area jets
0298         for(unsigned int jjet = 0; jjet < rhoA_pt->size(); jjet++)
0299         {
0300             if(std::find(matched_area_jets.begin(), matched_area_jets.end(), jjet) != matched_area_jets.end()) continue;
0301             if(rhoA_pt->at(jjet) < pTmin_reco) continue;
0302             unmatched_pt_area.push_back(rhoA_pt->at(jjet));
0303             unmatched_pt_area_unsub.push_back(rhoA_pt_unsub->at(jjet));
0304         }
0305         // fill unmatched reco mult jets
0306         for(unsigned int jjet = 0; jjet < mult_pt->size(); jjet++)
0307         {
0308             if(std::find(matched_mult_jets.begin(), matched_mult_jets.end(), jjet) != matched_mult_jets.end()) continue;
0309             if(mult_pt->at(jjet) < pTmin_reco) continue;
0310             unmatched_pt_mult.push_back(mult_pt->at(jjet));
0311             unmatched_pt_mult_unsub.push_back(mult_pt_unsub->at(jjet));
0312         }
0313         // fill unmatched reco iter jets
0314         for(unsigned int jjet = 0; jjet < iter_pt->size(); jjet++)
0315         {
0316             if(std::find(matched_iter_jets.begin(), matched_iter_jets.end(), jjet) != matched_iter_jets.end()) continue;
0317             if(iter_pt->at(jjet) < pTmin_reco) continue;
0318             unmatched_pt_iter.push_back(iter_pt->at(jjet));
0319             unmatched_pt_iter_unsub.push_back(iter_pt_unsub->at(jjet));
0320         }
0321
0322         // fill tree
0323         out_tree->Fill();
0324     }
0325
0326     // write output file
0327     fout->Write();
0328     fout->Close();
0329     // close input file
0330     fin.Close();
0331     return 0;
0332 }
0333