macro/plot/JES_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    JetIndicesMatcher jetmatcher{0.4, 30., 5.};
0015
0016    TFile* fout = new TFile(voi_stem(str_inp,".root").c_str(),"recreate");
0017
0018    array<TH1D*,10> harr_sub1, harr_rhoA;
0019    for (int i=0; i<10; ++i) {
0020      int i0 = i*10;
0021      int i1 = i*10+10;
0022      harr_sub1[i] = new TH1D(Form("JES_sub1_%i",i),
0023          Form("JES %i-%i%% MBD Centrality;p_{T}_{truth};#frac{p_{T}^{jet,calo}_{SUB1}-p_{T}^{truth}}{p_{T}^{truth}}",
0024            i0,i1), 300, -3., 3.);
0025      harr_rhoA[i] = new TH1D(Form("JES_rhoA_%i",i),
0026          Form("JES %i-%i%% MBD Centrality;p_{T}_{truth};#frac{(p_{T}^{jet,calo}-#rho^{md.bkgd}A_{jet})-p_{T}^{truth}}{p_{T}^{truth}}",
0027            i0,i1), 300, -3., 3.);
0028    }
0029
0030
0031    Long64_t nbytes = 0, nb = 0;
0032    for (Long64_t jentry=0; jentry<nentries;jentry++) {
0033       Long64_t ientry = LoadTree(jentry);
0034       if (ientry < 0) break;
0035       nb = fChain->GetEntry(jentry);   nbytes += nb;
0036       // if (Cut(ientry) < 0) continue;
0037
0038       int k = ((int) cent_mdb)/10;
0039       if (k < 0 || k > 9) {
0040         cout << " WARNING centrality in un expected bin" << endl;
0041         continue;
0042       }
0043
0044       // match the sub1 jets
0045       jetmatcher.reset();
0046       jetmatcher.add_truth (*TruthJetEta, *TruthJetPhi, *TruthJetPt);
0047       jetmatcher.add_reco  (*sub1JetEta,  *sub1JetPhi, *sub1JetPt);
0048       jetmatcher.do_matching();
0049
0050       for (auto pp : jetmatcher.match) {
0051         auto T = (*TruthJetPt)[pp.first];
0052         if (T<30) continue;
0053         auto M = (*sub1JetPt)[pp.second];
0054         harr_sub1[k]->Fill((M-T)/T);
0055       }
0056
0057       // match rhoA jets
0058       jetmatcher.reset();
0059       jetmatcher.add_truth (*TruthJetEta, *TruthJetPhi, *TruthJetPt);
0060       jetmatcher.add_reco  (*rhoAJetEta,  *rhoAJetPhi, *rhoAJetPtLessRhoA);
0061       jetmatcher.do_matching();
0062
0063       for (auto pp : jetmatcher.match) {
0064         auto T = (*TruthJetPt)[pp.first];
0065         if (T<30) continue;
0066         auto M = (*rhoAJetPtLessRhoA)[pp.second];
0067         harr_rhoA[k]->Fill((M-T)/T);
0068       }
0069
0070       // if (Cut(ientry) < 0) continue;
0071       if (false) if (jentry>10) {
0072         cout << " Breaking on jentry(" << jentry << ")" << endl;
0073         break;
0074       }
0075    }
0076
0077    for (int k=0;k<10;++k) {
0078      harr_sub1[k]->Write();
0079      harr_rhoA[k]->Write();
0080    }
0081    fout->Close();
0082 }