macro/plot/MatchFakeMiss_Sub1rhoA.cc

0001 #define Sub1rhoA_cxx
0002 #include "Sub1rhoA.h"
0003 #include <TH2.h>
0004 #include <TStyle.h>
0005 #include <TCanvas.h>
0006
0007 void Sub1rhoA::Loop(string str_inp)
0008 {
0009    if (fChain == 0) return;
0010
0011    Long64_t nentries = fChain->GetEntriesFast();
0012
0013
0014
0015    TFile* fout = new TFile(voi_stem(str_inp,".root").c_str(),"recreate");
0016
0017    array<TH1D*,10> harr_sub1, harr_rhoA;
0018
0019
0020    array<TH1D*,10> fake_sub1, fake_rhoA, miss_sub1, miss_rhoA;
0021    array<TH2D*,10> match_sub1, match_rhoA;
0022
0023    int   nbins_fake = 40;
0024    float lobin_fake = -20.;
0025    float hibin_fake = 80.;
0026    int   nbins_miss = 40;
0027    float lobin_miss = -20.;
0028    float hibin_miss = 80.;
0029
0030
0031    for (int i=0; i<10; ++i) {
0032      int i0 = i*10;
0033      int i1 = i*10+10;
0034
0035      fake_sub1[i] = new TH1D(Form("fake_sub1_%i",i), Form("cent: %i-%i%% MBC Cent;N_{jets} Fake;p_{T}^{SUB1}",i0, i1), nbins_fake, lobin_fake, hibin_fake);
0036      fake_rhoA[i] = new TH1D(Form("fake_rhoA_%i",i), Form("cent: %i-%i%% MBC Cent;N_{jets} Fake;p_{T}-#rho#timesA",i0, i1), nbins_fake, lobin_fake, hibin_fake);
0037
0038      miss_sub1[i] = new TH1D(Form("miss_sub1_%i",i), Form("cent: %i-%i%% MBC Cent;N_{jets} miss;p_{T}^{truth}",i0, i1), nbins_miss, lobin_miss, hibin_miss);
0039      miss_rhoA[i] = new TH1D(Form("miss_rhoA_%i",i), Form("cent: %i-%i%% MBC Cent;N_{jets} miss;p_{T}^{truth}",i0, i1), nbins_miss, lobin_miss, hibin_miss);
0040
0041      match_sub1[i] = new TH2D(Form("match_sub1_%i",i), Form("cent: %i-%i%% MBC Cent;p_{T}^{SUB1};p_{T}^{truth}",i0, i1),      nbins_fake, lobin_fake, hibin_fake, nbins_miss, lobin_miss, hibin_miss);
0042      match_rhoA[i] = new TH2D(Form("match_rhoA_%i",i), Form("cent: %i-%i%% MBC Cent;p_{T}-#rho#timesA;p_{T}^{truth}",i0, i1), nbins_fake, lobin_fake, hibin_fake, nbins_miss, lobin_miss, hibin_miss);
0043    }
0044
0045    JetIndicesMatcher jetmatcher{0.4, 30., 5.};
0046
0047    Long64_t nbytes = 0, nb = 0;
0048    for (Long64_t jentry=0; jentry<nentries;jentry++) {
0049       Long64_t ientry = LoadTree(jentry);
0050       if (ientry < 0) break;
0051       nb = fChain->GetEntry(jentry);   nbytes += nb;
0052       // if (Cut(ientry) < 0) continue;
0053
0054       int k = ((int) cent_mdb)/10;
0055       if (k < 0 || k > 9) {
0056         cout << " WARNING centrality in un expected bin" << endl;
0057         continue;
0058       }
0059
0060       // match the sub1 jets
0061       jetmatcher.reset();
0062       jetmatcher.add_truth (*TruthJetEta, *TruthJetPhi, *TruthJetPt);
0063       jetmatcher.add_reco  (*sub1JetEta,  *sub1JetPhi, *sub1JetPt);
0064       jetmatcher.do_matching();
0065
0066       for (auto pp : jetmatcher.match) {
0067         auto T = (*TruthJetPt)[pp.first];
0068         if (T<30) {
0069             cout << " WARNING WHY HERE!!! A0" << endl;
0070             continue;
0071         }
0072         auto M = (*sub1JetPt)[pp.second];
0073         match_sub1[k] ->Fill(M,T);
0074       }
0075       for (auto F : jetmatcher.indices_fake) {
0076         auto M = (*sub1JetPt)[F];
0077         fake_sub1[k]->Fill(M);
0078       }
0079       for (auto M : jetmatcher.indices_miss) {
0080         auto T = (*TruthJetPt)[M];
0081         miss_sub1[k]->Fill(T);
0082       }
0083
0084
0085       jetmatcher.reset();
0086       jetmatcher.add_truth (*TruthJetEta, *TruthJetPhi, *TruthJetPt);
0087       jetmatcher.add_reco  (*rhoAJetEta,  *rhoAJetPhi, *rhoAJetPtLessRhoA);
0088       jetmatcher.do_matching();
0089
0090       for (auto pp : jetmatcher.match) {
0091         auto T = (*TruthJetPt)[pp.first];
0092         if (T<30) {
0093             cout << " WARNING WHY HERE!!! A1" << endl;
0094             continue;
0095         }
0096         auto M = (*rhoAJetPtLessRhoA)[pp.second];
0097         match_rhoA[k]->Fill(M,T);
0098       }
0099
0100       for (auto F : jetmatcher.indices_fake) {
0101         auto pt = (*rhoAJetPtLessRhoA)[F];
0102         /* if (pt<6) cout << " pt["<<F<<"]: " << pt << endl; */
0103         fake_rhoA[k]->Fill(pt);
0104       }
0105
0106       for (auto M : jetmatcher.indices_miss) {
0107         auto pt = (*TruthJetPt)[M];
0108         miss_rhoA[k]->Fill(pt);
0109       }
0110
0111       // if (Cut(ientry) < 0) continue;
0112       if (false) if (jentry>10) {
0113         cout << " Breaking on jentry(" << jentry << ")" << endl;
0114         break;
0115       }
0116    }
0117
0118    for (int k=0;k<10;++k) {
0119      fake_sub1[k]->Write();
0120      miss_sub1[k]->Write();
0121      match_sub1[k]->Write();
0122      fake_rhoA[k]->Write();
0123      miss_rhoA[k]->Write();
0124      match_rhoA[k]->Write();
0125    }
0126    fout->Close();
0127 }