PPG09_JetSpectrum_Analysis/offline/analysis_data.C

0001 #include <TChain.h>
0002 #include <TFile.h>
0003 #include <TH1.h>
0004 #include <TH2.h>
0005 #include <TF1.h>
0006 #include <TMath.h>
0007 #include <iostream>
0008 #include <vector>
0009 #include <string>
0010 #include <cmath>
0011 #include "unfold_Def.h"
0012
0013 void get_leading_subleading_jet(int& leadingjet_index, int& subleadingjet_index, std::vector<float>* jet_et);
0014 bool match_leading_subleading_jet(float leadingjet_phi, float subleadingjet_phi);
0015 void get_jet_filter(std::vector<bool>& jet_filter, std::vector<float>* jet_e, std::vector<float>* jet_pt, std::vector<float>* jet_eta, float zvertex, float jet_radius);
0016 int get_zvertex_index(float zvertex);
0017 void get_calibjet(std::vector<float>& calibjet_pt, std::vector<float>& calibjet_eta, std::vector<float>& calibjet_phi, std::vector<bool>& jet_filter, std::vector<float>* jet_pt, std::vector<float>* jet_eta, std::vector<float>* jet_phi, const std::vector<std::vector<TF1*>>& f_corr, int zvertex_index, double jet_radius);
0018
0019 ////////////////////////////////////////// Main Function //////////////////////////////////////////
0020 void analysis_data(int runnumber, int nseg, int iseg, double jet_radius) {
0021   /////////////// General Set up ///////////////
0022   int jet_radius_index = (int)(10 * jet_radius);
0023   TFile *f_out = new TFile(Form("output_data/output_r0%d_%d_%d_%d.root", jet_radius_index, runnumber, iseg, iseg+nseg),"RECREATE");
0024
0025   /////////////// Read Files ///////////////
0026   const char* baseDirJet = "/sphenix/tg/tg01/jets/hanpuj/CaloDataAna_ppg09";
0027   TChain chain("ttree");
0028   for (int i = iseg; i < iseg + nseg; ++i) {
0029     chain.Add(Form("%s/Run%d/OutDir%d/output_data.root", baseDirJet, runnumber, i));
0030   }
0031   chain.SetBranchStatus("*", 0);
0032   // Read necessary branches for analysis.
0033   float zvertex; chain.SetBranchStatus("z_vertex", 1); chain.SetBranchAddress("z_vertex", &zvertex);
0034   int gl1_trigger_vector_scaled[64]; chain.SetBranchStatus("gl1_trigger_vector_scaled", 1); chain.SetBranchAddress("gl1_trigger_vector_scaled", gl1_trigger_vector_scaled);
0035   std::vector<float>* unsubjet_e = nullptr; chain.SetBranchStatus(Form("unsubjet0%d_e", jet_radius_index), 1); chain.SetBranchAddress(Form("unsubjet0%d_e", jet_radius_index), &unsubjet_e);
0036   std::vector<float>* unsubjet_pt = nullptr; chain.SetBranchStatus(Form("unsubjet0%d_pt", jet_radius_index), 1); chain.SetBranchAddress(Form("unsubjet0%d_pt", jet_radius_index), &unsubjet_pt);
0037   std::vector<float>* unsubjet_eta = nullptr; chain.SetBranchStatus(Form("unsubjet0%d_eta", jet_radius_index), 1); chain.SetBranchAddress(Form("unsubjet0%d_eta", jet_radius_index), &unsubjet_eta);
0038   std::vector<float>* unsubjet_phi = nullptr; chain.SetBranchStatus(Form("unsubjet0%d_phi", jet_radius_index), 1); chain.SetBranchAddress(Form("unsubjet0%d_phi", jet_radius_index), &unsubjet_phi);
0039   std::vector<float>* unsubjet_time = nullptr; chain.SetBranchStatus(Form("unsubjet0%d_time", jet_radius_index), 1); chain.SetBranchAddress(Form("unsubjet0%d_time", jet_radius_index), &unsubjet_time);
0040   // Read R=0.4 jet branches for event selection.
0041   std::vector<float>* unsubjet04_e_reader = nullptr, *unsubjet04_pt_reader = nullptr, *unsubjet04_eta_reader = nullptr, *unsubjet04_phi_reader = nullptr, *unsubjet04_time_reader = nullptr;
0042   if (jet_radius_index != 4) {
0043     chain.SetBranchStatus("unsubjet04_e", 1); chain.SetBranchAddress("unsubjet04_e", &unsubjet04_e_reader);
0044     chain.SetBranchStatus("unsubjet04_pt", 1); chain.SetBranchAddress("unsubjet04_pt", &unsubjet04_pt_reader);
0045     chain.SetBranchStatus("unsubjet04_eta", 1); chain.SetBranchAddress("unsubjet04_eta", &unsubjet04_eta_reader);
0046     chain.SetBranchStatus("unsubjet04_phi", 1); chain.SetBranchAddress("unsubjet04_phi", &unsubjet04_phi_reader);
0047     chain.SetBranchStatus("unsubjet04_time", 1); chain.SetBranchAddress("unsubjet04_time", &unsubjet04_time_reader);
0048   }
0049   std::vector<float>* &unsubjet04_e = (jet_radius_index != 4) ? unsubjet04_e_reader : unsubjet_e;
0050   std::vector<float>* &unsubjet04_pt = (jet_radius_index != 4) ? unsubjet04_pt_reader : unsubjet_pt;
0051   std::vector<float>* &unsubjet04_eta = (jet_radius_index != 4) ? unsubjet04_eta_reader : unsubjet_eta;
0052   std::vector<float>* &unsubjet04_phi = (jet_radius_index != 4) ? unsubjet04_phi_reader : unsubjet_phi;
0053   std::vector<float>* &unsubjet04_time = (jet_radius_index != 4) ? unsubjet04_time_reader : unsubjet_time;
0054
0055   /////////////// JES func ///////////////
0056   TFile *corrFile = new TFile("input_jescalib.root", "READ");
0057   if (!corrFile || corrFile->IsZombie()) {
0058     std::cout << "Error: cannot open JES_Calib_Default.root" << std::endl;
0059     return;
0060   }
0061   std::vector<std::vector<TF1 *>> f_corr;
0062   int nZvrtxBins = 5, nEtaBins = 4;
0063   f_corr.resize(nZvrtxBins);
0064   for (auto &row : f_corr) {
0065     row.resize(nEtaBins);
0066   }
0067   for (int iz = 0; iz < nZvrtxBins; ++iz) {
0068     if (iz <= 2 && jet_radius_index <= 4) {
0069       for (int ieta = 0; ieta < nEtaBins; ++ieta) {
0070         f_corr[iz][ieta] = (TF1*)corrFile->Get(Form("JES_Calib_Func_R0%d_Z%d_Eta%d", jet_radius_index, iz, ieta));
0071         if (!f_corr[iz][ieta]) {
0072           std::cout << "Error: cannot open f_corr for zbin " << iz << " and etabin " << ieta << std::endl;
0073           return;
0074         }
0075       }
0076     } else if (iz <= 2 && jet_radius_index > 4) {
0077       for (int ieta = 0; ieta < nEtaBins; ++ieta) {
0078         f_corr[iz][ieta] = (TF1*)corrFile->Get(Form("JES_Calib_Func_R0%d_Default", jet_radius_index));
0079         if (!f_corr[iz][ieta]) {
0080           std::cout << "Error: cannot open f_corr for radius " << jet_radius_index << ", zbin " << iz << " and etabin " << ieta << std::endl;
0081           return;
0082         }
0083       }
0084     } else if (iz > 2) {
0085       f_corr[iz][0] = (TF1*)corrFile->Get(Form("JES_Calib_Func_R0%d_Z%d", jet_radius_index, iz));
0086       if (!f_corr[iz][0]) {
0087         std::cout << "Error: cannot open f_corr for radius " << jet_radius_index << " and zbin " << iz << std::endl;
0088         return;
0089       }
0090     } else {
0091       std::cout << "Error: invalid z-vertex bin: iz = " << iz << " with jet radius = " << jet_radius_index << std::endl;
0092       return;
0093     }
0094   }
0095
0096   /////////////// Jet Trigger Efficiency ///////////////
0097   TFile *f_jettrigger = new TFile("input_jetefficiency.root", "READ");
0098   if (!f_jettrigger) {
0099     std::cout << "Error: cannot open input_jetefficiency.root" << std::endl;
0100     return;
0101   }
0102   TF1* f_jettrigger_nominal = (TF1*)f_jettrigger->Get(Form("jettrig_0%d_pt_nominal", jet_radius_index));
0103   TF1* f_jettrigger_up = (TF1*)f_jettrigger->Get(Form("jettrig_0%d_pt_up", jet_radius_index));
0104   TF1* f_jettrigger_down = (TF1*)f_jettrigger->Get(Form("jettrig_0%d_pt_down", jet_radius_index));
0105   if (!f_jettrigger_nominal || !f_jettrigger_up || !f_jettrigger_down) {
0106     std::cout << "Error: cannot open jet trigger functions" << std::endl;
0107     return;
0108   }
0109
0110   /////////////// Histograms ///////////////
0111   TH1D *h_recojet_pt_record = new TH1D("h_recojet_pt_record", ";p_{T} [GeV]", 1000, 0, 100);
0112   TH1D *h_calibjet_pt_record = new TH1D("h_calibjet_pt_record", ";p_{T} [GeV]", 1000, 0, 100);
0113   TH1D *h_calibjet_pt_scaled_nominal = new TH1D("h_calibjet_pt_scaled_nominal", ";p_{T} [GeV]", calibnpt, calibptbins);
0114   TH1D *h_calibjet_pt_scaled_jettrigup = new TH1D("h_calibjet_pt_scaled_jettrigup", ";p_{T} [GeV]", calibnpt, calibptbins);
0115   TH1D *h_calibjet_pt_scaled_jettrigdown = new TH1D("h_calibjet_pt_scaled_jettrigdown", ";p_{T} [GeV]", calibnpt, calibptbins);
0116   TH1D *h_calibjet_pt_scaled_jettimingup = new TH1D("h_calibjet_pt_scaled_jettimingup", ";p_{T} [GeV]", calibnpt, calibptbins);
0117   TH1D *h_calibjet_pt_scaled_jettimingdown = new TH1D("h_calibjet_pt_scaled_jettimingdown", ";p_{T} [GeV]", calibnpt, calibptbins);
0118
0119   /////////////// Event Loop ///////////////
0120   std::cout << "Data analysis started." << std::endl;
0121   Long64_t n_events = chain.GetEntries();
0122   std::cout << "Total number of events: " << n_events << std::endl;
0123   // Event variables setup.
0124   std::vector<bool> jet_filter;
0125   std::vector<float> calibjet_pt, calibjet_eta, calibjet_phi;
0126   for (Long64_t ie = 0; ie < n_events; ++ie) { // event loop start
0127     // Load event.
0128     if (ie % 10000 == 0) {
0129       std::cout << "Processing event " << ie << "..." << std::endl;
0130     }
0131     chain.GetEntry(ie);
0132
0133     // Trigger check.
0134     if (gl1_trigger_vector_scaled[22] != 1) continue;
0135
0136     // Get z-vertex index and re-assign no-zvertex events.
0137     int zvertex_index = get_zvertex_index(zvertex);
0138     if (zvertex < -990) zvertex = 0;
0139
0140     // Event selection.
0141     bool pass_dijet = false;
0142     bool pass_timing = false;
0143     bool pass_njet = false;
0144     // Leading and subleading jet with energy and R = 0.4 for dijet and timing requirement.
0145     int leadingunsubjet04_index = -1;
0146     int subleadingunsubjet04_index = -1;
0147     get_leading_subleading_jet(leadingunsubjet04_index, subleadingunsubjet04_index, unsubjet04_e);
0148     if (leadingunsubjet04_index < 0) continue;
0149     // Dijet requirement.
0150     if (subleadingunsubjet04_index >= 0 && unsubjet04_e->at(subleadingunsubjet04_index) / (float) unsubjet04_e->at(leadingunsubjet04_index) > 0.3) {
0151       pass_dijet = match_leading_subleading_jet(unsubjet04_phi->at(leadingunsubjet04_index), unsubjet04_phi->at(subleadingunsubjet04_index));
0152     }
0153     // Timing requirement.
0154     if (subleadingunsubjet04_index >= 0 && fabs(unsubjet04_time->at(leadingunsubjet04_index)*17.6+2) < 6 && fabs(unsubjet04_time->at(leadingunsubjet04_index)*17.6 - unsubjet04_time->at(subleadingunsubjet04_index)*17.6) < 3) {
0155       pass_timing = true;
0156     }
0157     // Njet requirement.
0158     int n_jets = 0;
0159     for (int ij = 0; ij < unsubjet04_e->size(); ++ij) {
0160       if (unsubjet04_e->at(ij) > 5.0) {
0161         n_jets++;
0162       }
0163     }
0164     if (n_jets <= 9) {
0165       pass_njet = true;
0166     }
0167     if (!pass_dijet || !pass_timing || !pass_njet) continue;
0168
0169     // Leading jet with pT for trigger efficiency.
0170     int leadingunsubjet_index = -1;
0171     int subleadingunsubjet_index = -1;
0172     get_leading_subleading_jet(leadingunsubjet_index, subleadingunsubjet_index, unsubjet_pt);
0173     if (leadingunsubjet_index < 0) continue;
0174     double leadingunsubjet_pt = unsubjet_pt->at(leadingunsubjet_index);
0175
0176     // Get jet trigger scale factors.
0177     double jettrigger_nominal = f_jettrigger_nominal->Eval(leadingunsubjet_pt);
0178     if (jettrigger_nominal < 0.01) jettrigger_nominal = 0.01;
0179     double scale_jettrig_nominal = 1.0 / jettrigger_nominal;
0180
0181     double jettrigger_up = f_jettrigger_up->Eval(leadingunsubjet_pt);
0182     if (jettrigger_up < 0.01) jettrigger_up = 0.01;
0183     double scale_jettrig_up = 1.0 / jettrigger_up;
0184
0185     double jettrigger_down = f_jettrigger_down->Eval(leadingunsubjet_pt);
0186     if (jettrigger_down < 0.01) jettrigger_down = 0.01;
0187     double scale_jettrig_down = 1.0 / jettrigger_down;
0188
0189     // Get jet timing scale factors.
0190     double jettiming_nominal = 0.92;
0191     double scale_jettiming_nominal = 1.0 / jettiming_nominal;
0192
0193     double jettiming_up = 0.95;
0194     double scale_jettiming_up = 1.0 / jettiming_up;
0195
0196     double jettiming_down = 0.89;
0197     double scale_jettiming_down = 1.0 / jettiming_down;
0198
0199     // Fill histograms.
0200     get_jet_filter(jet_filter, unsubjet_e, unsubjet_pt, unsubjet_eta, zvertex, jet_radius);
0201     for (int ij = 0; ij < jet_filter.size(); ++ij) {
0202       if (jet_filter.at(ij)) continue;
0203       h_recojet_pt_record->Fill(unsubjet_pt->at(ij));
0204     }
0205
0206     get_calibjet(calibjet_pt, calibjet_eta, calibjet_phi, jet_filter, unsubjet_pt, unsubjet_eta, unsubjet_phi, f_corr, zvertex_index, jet_radius);
0207     for (int ij = 0; ij < calibjet_pt.size(); ++ij) {
0208       h_calibjet_pt_record->Fill(calibjet_pt[ij]);
0209       h_calibjet_pt_scaled_nominal->Fill(calibjet_pt[ij], scale_jettrig_nominal * scale_jettiming_nominal);
0210       h_calibjet_pt_scaled_jettrigup->Fill(calibjet_pt[ij], scale_jettrig_up * scale_jettiming_nominal);
0211       h_calibjet_pt_scaled_jettrigdown->Fill(calibjet_pt[ij], scale_jettrig_down * scale_jettiming_nominal);
0212       h_calibjet_pt_scaled_jettimingup->Fill(calibjet_pt[ij], scale_jettrig_nominal * scale_jettiming_up);
0213       h_calibjet_pt_scaled_jettimingdown->Fill(calibjet_pt[ij], scale_jettrig_nominal * scale_jettiming_down);
0214     }
0215   } // event loop end
0216
0217   // Write histograms.
0218   std::cout << "Writing histograms..." << std::endl;
0219   f_out->cd();
0220   h_recojet_pt_record->Write();
0221   h_calibjet_pt_record->Write();
0222   h_calibjet_pt_scaled_nominal->Write();
0223   h_calibjet_pt_scaled_jettrigup->Write();
0224   h_calibjet_pt_scaled_jettrigdown->Write();
0225   h_calibjet_pt_scaled_jettimingup->Write();
0226   h_calibjet_pt_scaled_jettimingdown->Write();
0227   f_out->Close();
0228   std::cout << "All done!" << std::endl;
0229 }
0230
0231 ////////////////////////////////////////// Helper Functions //////////////////////////////////////////
0232 float get_deta(float eta1, float eta2) {
0233   return eta1 - eta2;
0234 }
0235
0236 float get_dphi(float phi1, float phi2) {
0237   float dphi1 = phi1 - phi2;
0238   float dphi2 = phi1 - phi2 + 2*TMath::Pi();
0239   float dphi3 = phi1 - phi2 - 2*TMath::Pi();
0240   if (fabs(dphi1) > fabs(dphi2)) {
0241     dphi1 = dphi2;
0242   }
0243   if (fabs(dphi1) > fabs(dphi3)) {
0244     dphi1 = dphi3;
0245   }
0246   return dphi1;
0247 }
0248
0249 float get_dR(float eta1, float phi1, float eta2, float phi2) {
0250   float deta = get_deta(eta1, eta2);
0251   float dphi = get_dphi(phi1, phi2);
0252   return sqrt(deta*deta + dphi*dphi);
0253 }
0254
0255 float get_emcal_mineta_zcorrected(float zvertex) {
0256   float minz_EM = -130.23;
0257   float radius_EM = 93.5;
0258   float z = minz_EM - zvertex;
0259   float eta_zcorrected = asinh(z / (float)radius_EM);
0260   return eta_zcorrected;
0261 }
0262
0263 float get_emcal_maxeta_zcorrected(float zvertex) {
0264   float maxz_EM = 130.23;
0265   float radius_EM = 93.5;
0266   float z = maxz_EM - zvertex;
0267   float eta_zcorrected = asinh(z / (float)radius_EM);
0268   return eta_zcorrected;
0269 }
0270
0271 float get_ihcal_mineta_zcorrected(float zvertex) {
0272   float minz_IH = -170.299;
0273   float radius_IH = 127.503;
0274   float z = minz_IH - zvertex;
0275   float eta_zcorrected = asinh(z / (float)radius_IH);
0276   return eta_zcorrected;
0277 }
0278
0279 float get_ihcal_maxeta_zcorrected(float zvertex) {
0280   float maxz_IH = 170.299;
0281   float radius_IH = 127.503;
0282   float z = maxz_IH - zvertex;
0283   float eta_zcorrected = asinh(z / (float)radius_IH);
0284   return eta_zcorrected;
0285 }
0286
0287 float get_ohcal_mineta_zcorrected(float zvertex) {
0288   float minz_OH = -301.683;
0289   float radius_OH = 225.87;
0290   float z = minz_OH - zvertex;
0291   float eta_zcorrected = asinh(z / (float)radius_OH);
0292   return eta_zcorrected;
0293 }
0294
0295 float get_ohcal_maxeta_zcorrected(float zvertex) {
0296   float maxz_OH = 301.683;
0297   float radius_OH = 225.87;
0298   float z = maxz_OH - zvertex;
0299   float eta_zcorrected = asinh(z / (float)radius_OH);
0300   return eta_zcorrected;
0301 }
0302
0303 bool check_bad_jet_eta(float jet_eta, float zvertex, float jet_radius) {
0304   float emcal_mineta = get_emcal_mineta_zcorrected(zvertex);
0305   float emcal_maxeta = get_emcal_maxeta_zcorrected(zvertex);
0306   float ihcal_mineta = get_ihcal_mineta_zcorrected(zvertex);
0307   float ihcal_maxeta = get_ihcal_maxeta_zcorrected(zvertex);
0308   float ohcal_mineta = get_ohcal_mineta_zcorrected(zvertex);
0309   float ohcal_maxeta = get_ohcal_maxeta_zcorrected(zvertex);
0310   float minlimit = emcal_mineta;
0311   if (ihcal_mineta > minlimit) minlimit = ihcal_mineta;
0312   if (ohcal_mineta > minlimit) minlimit = ohcal_mineta;
0313   float maxlimit = emcal_maxeta;
0314   if (ihcal_maxeta < maxlimit) maxlimit = ihcal_maxeta;
0315   if (ohcal_maxeta < maxlimit) maxlimit = ohcal_maxeta;
0316   minlimit += jet_radius;
0317   maxlimit -= jet_radius;
0318   return jet_eta < minlimit || jet_eta > maxlimit;
0319 }
0320
0321 int get_zvertex_index(float zvertex) {
0322   if (zvertex >= -60.0 && zvertex < -30.0) return 1;
0323   else if (zvertex >= -30.0 && zvertex < 30.0) return 0;
0324   else if (zvertex >= 30.0 && zvertex < 60) return 2;
0325   else if ((zvertex < -60.0 && zvertex >= -900) || zvertex >= 60.0) return 3;
0326   else if (zvertex < -900) return 4;
0327   return -1; // should not reach here
0328 }
0329
0330 int get_eta_index(float jet_eta, int zvertex_index, double jet_radius) {
0331   float eta_low = 0;
0332   float eta_high = 0;
0333   if (zvertex_index == 0) {// -30 < zvrtx < 30
0334     eta_low = -1.2 + jet_radius;
0335     eta_high = 1.2 - jet_radius;
0336   } else if (zvertex_index == 1) {// -60 < zvrtx < -30
0337     eta_low = -0.95 + jet_radius;
0338     eta_high = 1.25 - jet_radius;
0339   } else if (zvertex_index == 2) {// 30 < zvrtx < inf
0340     eta_low = -1.25 + jet_radius;
0341     eta_high = 0.95 - jet_radius;
0342   } else if (zvertex_index == 3 || zvertex_index == 4) {
0343     return 0;
0344   }
0345   float threshold1 = eta_low + ((eta_high - eta_low) / 4.0);
0346   float threshold2 = eta_low + ((eta_high - eta_low) / 2.0);
0347   float threshold3 = eta_low + (3.0 * (eta_high - eta_low) / 4.0);
0348   if (jet_eta < threshold1) return 0;  // -inf to threshold1
0349   if (jet_eta < threshold2) return 1;  // threshold1 to threshold2
0350   if (jet_eta < threshold3) return 2;  // threshold2 to threshold3
0351   return 3;  // threshold3 to inf
0352 }
0353
0354 ////////////////////////////////////////// Functions //////////////////////////////////////////
0355 void get_leading_subleading_jet(int& leadingjet_index, int& subleadingjet_index, std::vector<float>* jet_et) {
0356   leadingjet_index = -1;
0357   subleadingjet_index = -1;
0358   float leadingjet_et = -9999;
0359   float subleadingjet_et = -9999;
0360   for (int ij = 0; ij < jet_et->size(); ++ij) {
0361     float jetet = jet_et->at(ij);
0362     if (jetet > leadingjet_et) {
0363       subleadingjet_et = leadingjet_et;
0364       subleadingjet_index = leadingjet_index;
0365       leadingjet_et = jetet;
0366       leadingjet_index = ij;
0367     } else if (jetet > subleadingjet_et) {
0368       subleadingjet_et = jetet;
0369       subleadingjet_index = ij;
0370     }
0371   }
0372 }
0373
0374 bool match_leading_subleading_jet(float leadingjet_phi, float subleadingjet_phi) {
0375   float dijet_min_phi = 3*TMath::Pi()/4.;
0376   float dphi = fabs(get_dphi(leadingjet_phi, subleadingjet_phi));
0377   return dphi > dijet_min_phi;
0378 }
0379
0380 void get_jet_filter(std::vector<bool>& jet_filter, std::vector<float>* jet_e, std::vector<float>* jet_pt, std::vector<float>* jet_eta, float zvertex, float jet_radius) {
0381   jet_filter.clear();
0382   int njet = jet_e->size();
0383   for (int ij = 0; ij < njet; ++ij) {
0384     jet_filter.push_back(jet_e->at(ij) < 0 || check_bad_jet_eta(jet_eta->at(ij), zvertex, jet_radius) || jet_pt->at(ij) < 1);
0385   }
0386 }
0387
0388 void get_calibjet(std::vector<float>& calibjet_pt, std::vector<float>& calibjet_eta, std::vector<float>& calibjet_phi, std::vector<bool>& jet_filter, std::vector<float>* jet_pt, std::vector<float>* jet_eta, std::vector<float>* jet_phi, const std::vector<std::vector<TF1*>>& f_corr, int zvertex_index, double jet_radius) {
0389   calibjet_pt.clear();
0390   calibjet_eta.clear();
0391   calibjet_phi.clear();
0392   for (int ij = 0; ij < jet_filter.size(); ++ij) {
0393     if (jet_filter.at(ij)) continue;
0394     int eta_index;
0395     if (zvertex_index == 0 || zvertex_index == 1 || zvertex_index == 2) {
0396       eta_index = get_eta_index(jet_eta->at(ij), zvertex_index, jet_radius);
0397     } else {
0398       eta_index = 0;
0399     }
0400     double calib_pt = f_corr[zvertex_index][eta_index]->Eval(jet_pt->at(ij));
0401     if (calib_pt < calibptbins[0] || calib_pt > calibptbins[calibnpt]) continue;
0402     calibjet_pt.push_back(calib_pt);
0403     calibjet_eta.push_back(jet_eta->at(ij));
0404     calibjet_phi.push_back(jet_phi->at(ij));
0405   }
0406 }