PPG09_JetSpectrum_Analysis/offline/analysis_sim.C

0001 #include <TChain.h>
0002 #include <TFile.h>
0003 #include <TH1F.h>
0004 #include <TH2F.h>
0005 #include <TF1.h>
0006 #include <TMath.h>
0007 #include <iostream>
0008 #include <vector>
0009 #include <string>
0010 #include <cmath>
0011 #include "TRandom.h"
0012 #include "TRandom3.h"
0013 #include "unfold_Def.h"
0014
0015 void get_leading_subleading_jet(int& leadingjet_index, int& subleadingjet_index, std::vector<float>* jet_et);
0016 bool match_leading_subleading_jet(float leadingjet_phi, float subleadingjet_phi);
0017 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);
0018 int get_zvertex_index(float zvertex);
0019 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, float jes_para, float jer_para, bool jet_background);
0020 void get_truthjet(std::vector<float>& goodtruthjet_pt, std::vector<float>& goodtruthjet_eta, std::vector<float>& goodtruthjet_phi, std::vector<bool>& truthjet_filter, std::vector<float>* jet_pt, std::vector<float>* jet_eta, std::vector<float>* jet_phi);
0021 void match_meas_truth(std::vector<float>& meas_eta, std::vector<float>& meas_phi, std::vector<float>& meas_matched, std::vector<float>& truth_eta, std::vector<float>& truth_phi, std::vector<float>& truth_matched, float jet_radius);
0022 void fill_response_matrix(TH1D*& h_truth, TH1D*& h_meas, TH2D*& h_resp, TH1D*& h_fake, TH1D*& h_miss, double scale, std::vector<float>& meas_pt, std::vector<float>& meas_matched, std::vector<float>& truth_pt, std::vector<float>& truth_matched);
0023 TRandom3 randGen(1234);
0024
0025 ////////////////////////////////////////// Main Function //////////////////////////////////////////
0026 void analysis_sim(std::string runtype, int nseg, int iseg, double jet_radius) {
0027   ////////// General Set up //////////
0028   int jet_radius_index = (int)(10 * jet_radius);
0029   TFile *f_out = new TFile(Form("output_sim/output_r0%d_%s_%d_%d.root", jet_radius_index, runtype.c_str(), iseg, iseg+nseg), "RECREATE");
0030
0031   double truthjet_pt_min = 0, truthjet_pt_max = 3000;
0032   if (runtype == "MB") {
0033     if (jet_radius == 0.2) {truthjet_pt_min = 0; truthjet_pt_max = 5;}
0034     else if (jet_radius == 0.3) {truthjet_pt_min = 0; truthjet_pt_max = 6;}
0035     else if (jet_radius == 0.4) {truthjet_pt_min = 0; truthjet_pt_max = 7;}
0036     else if (jet_radius == 0.5) {truthjet_pt_min = 0; truthjet_pt_max = 10;}
0037     else if (jet_radius == 0.6) {truthjet_pt_min = 0; truthjet_pt_max = 11;}
0038   } else if (runtype == "Jet5GeV") {
0039     if (jet_radius == 0.2) {truthjet_pt_min = 5; truthjet_pt_max = 12;}
0040     else if (jet_radius == 0.3) {truthjet_pt_min = 6; truthjet_pt_max = 13;}
0041     else if (jet_radius == 0.4) {truthjet_pt_min = 7; truthjet_pt_max = 14;}
0042     else if (jet_radius == 0.5) {truthjet_pt_min = 10; truthjet_pt_max = 15;}
0043     else if (jet_radius == 0.6) {truthjet_pt_min = 11; truthjet_pt_max = 17;}
0044   } else if (runtype == "Jet10GeV") {
0045     if (jet_radius == 0.2) {truthjet_pt_min = 12; truthjet_pt_max = 15;}
0046     else if (jet_radius == 0.3) {truthjet_pt_min = 13; truthjet_pt_max = 16;}
0047     else if (jet_radius == 0.4) {truthjet_pt_min = 14; truthjet_pt_max = 17;}
0048     else if (jet_radius == 0.5) {truthjet_pt_min = 15; truthjet_pt_max = 24;}
0049     else if (jet_radius == 0.6) {truthjet_pt_min = 17; truthjet_pt_max = 26;}
0050   } else if (runtype == "Jet15GeV") {
0051     if (jet_radius == 0.2) {truthjet_pt_min = 15; truthjet_pt_max = 20;}
0052     else if (jet_radius == 0.3) {truthjet_pt_min = 16; truthjet_pt_max = 21;}
0053     else if (jet_radius == 0.4) {truthjet_pt_min = 17; truthjet_pt_max = 22;}
0054     else if (jet_radius == 0.5) {truthjet_pt_min = 24; truthjet_pt_max = 30;}
0055     else if (jet_radius == 0.6) {truthjet_pt_min = 26; truthjet_pt_max = 35;}
0056   } else if (runtype == "Jet20GeV") {
0057     if (jet_radius == 0.2) {truthjet_pt_min = 20; truthjet_pt_max = 31;}
0058     else if (jet_radius == 0.3) {truthjet_pt_min = 21; truthjet_pt_max = 33;}
0059     else if (jet_radius == 0.4) {truthjet_pt_min = 22; truthjet_pt_max = 35;}
0060     else if (jet_radius == 0.5) {truthjet_pt_min = 30; truthjet_pt_max = 40;}
0061     else if (jet_radius == 0.6) {truthjet_pt_min = 35; truthjet_pt_max = 45;}
0062   } else if (runtype == "Jet30GeV") {
0063     if (jet_radius == 0.2) {truthjet_pt_min = 31; truthjet_pt_max = 50;}
0064     else if (jet_radius == 0.3) {truthjet_pt_min = 33; truthjet_pt_max = 51;}
0065     else if (jet_radius == 0.4) {truthjet_pt_min = 35; truthjet_pt_max = 52;}
0066     else if (jet_radius == 0.5) {truthjet_pt_min = 40; truthjet_pt_max = 60;}
0067     else if (jet_radius == 0.6) {truthjet_pt_min = 45; truthjet_pt_max = 63;}
0068   } else if (runtype == "Jet50GeV") {
0069     if (jet_radius == 0.2) {truthjet_pt_min = 50; truthjet_pt_max = 70;}
0070     else if (jet_radius == 0.3) {truthjet_pt_min = 51; truthjet_pt_max = 70;}
0071     else if (jet_radius == 0.4) {truthjet_pt_min = 52; truthjet_pt_max = 71;}
0072     else if (jet_radius == 0.5) {truthjet_pt_min = 60; truthjet_pt_max = 75;}
0073     else if (jet_radius == 0.6) {truthjet_pt_min = 63; truthjet_pt_max = 79;}
0074   } else if (runtype == "Jet70GeV") {
0075     if (jet_radius == 0.2) {truthjet_pt_min = 70; truthjet_pt_max = 3000;}
0076     else if (jet_radius == 0.3) {truthjet_pt_min = 70; truthjet_pt_max = 3000;}
0077     else if (jet_radius == 0.4) {truthjet_pt_min = 71; truthjet_pt_max = 3000;}
0078     else if (jet_radius == 0.5) {truthjet_pt_min = 75; truthjet_pt_max = 3000;}
0079     else if (jet_radius == 0.6) {truthjet_pt_min = 79; truthjet_pt_max = 3000;}
0080   }
0081
0082   ////////// Files //////////
0083   TChain chain("ttree");
0084   for (int i = iseg; i < iseg + nseg; ++i) {
0085     chain.Add(Form("/sphenix/tg/tg01/jets/hanpuj/JES_MC_run28/%s/OutDir%d/output_sim.root", runtype.c_str(), i));
0086   }
0087   chain.SetBranchStatus("*", 0);
0088   // Read necessary branches for analysis
0089   float zvertex; chain.SetBranchStatus("z_vertex", 1); chain.SetBranchAddress("z_vertex", &zvertex);
0090   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);
0091   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);
0092   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);
0093   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);
0094   std::vector<float>* truthjet_e = nullptr; chain.SetBranchStatus(Form("truthjet0%d_e", jet_radius_index), 1); chain.SetBranchAddress(Form("truthjet0%d_e", jet_radius_index), &truthjet_e);
0095   std::vector<float>* truthjet_pt = nullptr; chain.SetBranchStatus(Form("truthjet0%d_pt", jet_radius_index), 1); chain.SetBranchAddress(Form("truthjet0%d_pt", jet_radius_index), &truthjet_pt);
0096   std::vector<float>* truthjet_eta = nullptr; chain.SetBranchStatus(Form("truthjet0%d_eta", jet_radius_index), 1); chain.SetBranchAddress(Form("truthjet0%d_eta", jet_radius_index), &truthjet_eta);
0097   std::vector<float>* truthjet_phi = nullptr; chain.SetBranchStatus(Form("truthjet0%d_phi", jet_radius_index), 1); chain.SetBranchAddress(Form("truthjet0%d_phi", jet_radius_index), &truthjet_phi);
0098   // Read R=0.4 jet branches for event selection.
0099   std::vector<float>* unsubjet04_e_reader = nullptr, *unsubjet04_pt_reader = nullptr, *unsubjet04_eta_reader = nullptr, *unsubjet04_phi_reader = nullptr;
0100   if (jet_radius_index != 4) {
0101     chain.SetBranchStatus("unsubjet04_e", 1); chain.SetBranchAddress("unsubjet04_e", &unsubjet04_e_reader);
0102     chain.SetBranchStatus("unsubjet04_pt", 1); chain.SetBranchAddress("unsubjet04_pt", &unsubjet04_pt_reader);
0103     chain.SetBranchStatus("unsubjet04_eta", 1); chain.SetBranchAddress("unsubjet04_eta", &unsubjet04_eta_reader);
0104     chain.SetBranchStatus("unsubjet04_phi", 1); chain.SetBranchAddress("unsubjet04_phi", &unsubjet04_phi_reader);
0105   }
0106   std::vector<float>* &unsubjet04_e = (jet_radius_index != 4) ? unsubjet04_e_reader : unsubjet_e;
0107   std::vector<float>* &unsubjet04_pt = (jet_radius_index != 4) ? unsubjet04_pt_reader : unsubjet_pt;
0108   std::vector<float>* &unsubjet04_eta = (jet_radius_index != 4) ? unsubjet04_eta_reader : unsubjet_eta;
0109   std::vector<float>* &unsubjet04_phi = (jet_radius_index != 4) ? unsubjet04_phi_reader : unsubjet_phi;
0110
0111   /////////////// JES func ///////////////
0112   TFile *corrFile = new TFile("input_jescalib.root", "READ");
0113   if (!corrFile || corrFile->IsZombie()) {
0114     std::cout << "Error: cannot open JES_Calib_Default.root" << std::endl;
0115     return;
0116   }
0117   std::vector<std::vector<TF1 *>> f_corr;
0118   int nZvrtxBins = 5, nEtaBins = 4;
0119   f_corr.resize(nZvrtxBins);
0120   for (auto &row : f_corr) {
0121     row.resize(nEtaBins);
0122   }
0123   for (int iz = 0; iz < nZvrtxBins; ++iz) {
0124     if (iz <= 2 && jet_radius_index <= 4) {
0125       for (int ieta = 0; ieta < nEtaBins; ++ieta) {
0126         f_corr[iz][ieta] = (TF1*)corrFile->Get(Form("JES_Calib_Func_R0%d_Z%d_Eta%d", jet_radius_index, iz, ieta));
0127         if (!f_corr[iz][ieta]) {
0128           std::cout << "Error: cannot open f_corr for zbin " << iz << " and etabin " << ieta << std::endl;
0129           return;
0130         }
0131       }
0132     } else if (iz <= 2 && jet_radius_index > 4) {
0133       for (int ieta = 0; ieta < nEtaBins; ++ieta) {
0134         f_corr[iz][ieta] = (TF1*)corrFile->Get(Form("JES_Calib_Func_R0%d_Default", jet_radius_index));
0135         if (!f_corr[iz][ieta]) {
0136           std::cout << "Error: cannot open f_corr for radius " << jet_radius_index << ", zbin " << iz << " and etabin " << ieta << std::endl;
0137           return;
0138         }
0139       }
0140     } else if (iz > 2) {
0141       f_corr[iz][0] = (TF1*)corrFile->Get(Form("JES_Calib_Func_R0%d_Z%d", jet_radius_index, iz));
0142       if (!f_corr[iz][0]) {
0143         std::cout << "Error: cannot open f_corr for radius " << jet_radius_index << " and zbin " << iz << std::endl;
0144         return;
0145       }
0146     } else {
0147       std::cout << "Error: invalid z-vertex bin: iz = " << iz << " with jet radius = " << jet_radius_index << std::endl;
0148       return;
0149     }
0150   }
0151
0152   /////////////// Z-vertex Reweighting ///////////////
0153   TFile *f_zvertexreweight = new TFile("input_zvertexreweight.root", "READ");
0154   if (!f_zvertexreweight || f_zvertexreweight->IsZombie()) {
0155     std::cout << "Error: cannot open input_zvertexreweight.root" << std::endl;
0156     return;
0157   }
0158   TH1D *h_zvertex_reweight = (TH1D*)f_zvertexreweight->Get("h_zvertex_reweight");
0159   if (!h_zvertex_reweight) {
0160     std::cout << "Error: cannot open h_zvertex_reweight" << std::endl;
0161     return;
0162   }
0163
0164   ////////// Histograms //////////
0165   TH1D *h_event = new TH1D("h_event", ";Event Number", 1, 0, 1);
0166   TH1D *h_recojet_pt_record = new TH1D("h_recojet_pt_record", ";p_{T} [GeV]", 1000, 0, 100);
0167   TH1D *h_recojet_pt_passcut_record = new TH1D("h_recojet_pt_passcut_record", ";p_{T} [GeV]", 1000, 0, 100);
0168   TH1D *h_truthjet_pt_record = new TH1D("h_truthjet_pt_record", ";p_{T}^{Truth jet} [GeV]", 1000, 0, 100);
0169   TH1D *h_truthjet_pt_passcut_record = new TH1D("h_truthjet_pt_passcut_record", ";p_{T}^{Truth jet} [GeV]", 1000, 0, 100);
0170
0171   TH1D *h_truth_nominal = new TH1D("h_truth_nominal", ";p_{T}^{Truth jet} [GeV]", truthnpt, truthptbins);
0172   TH1D *h_measure_nominal = new TH1D("h_measure_nominal", ";p_{T}^{Calib jet} [GeV]", calibnpt, calibptbins);
0173   TH2D *h_respmatrix_nominal = new TH2D("h_respmatrix_nominal", ";p_{T}^{Calib jet} [GeV];p_{T}^{Truth jet} [GeV]", calibnpt, calibptbins, truthnpt, truthptbins);
0174   TH1D *h_fake_nominal = new TH1D("h_fake_nominal", ";p_{T}^{Calib jet} [GeV]", calibnpt, calibptbins);
0175   TH1D *h_miss_nominal = new TH1D("h_miss_nominal", ";p_{T}^{Truth jet} [GeV]", truthnpt, truthptbins);
0176
0177   // JES up/down histograms
0178   TH1D *h_truth_jesup = new TH1D("h_truth_jesup", ";p_{T}^{Truth jet} [GeV]", truthnpt, truthptbins);
0179   TH1D *h_measure_jesup = new TH1D("h_measure_jesup", ";p_{T}^{Calib jet} [GeV]", calibnpt, calibptbins);
0180   TH2D *h_respmatrix_jesup = new TH2D("h_respmatrix_jesup", ";p_{T}^{Calib jet} [GeV];p_{T}^{Truth jet} [GeV]", calibnpt, calibptbins, truthnpt, truthptbins);
0181   TH1D *h_fake_jesup = new TH1D("h_fake_jesup", ";p_{T}^{Calib jet} [GeV]", calibnpt, calibptbins);
0182   TH1D *h_miss_jesup = new TH1D("h_miss_jesup", ";p_{T}^{Truth jet} [GeV]", truthnpt, truthptbins);
0183
0184   TH1D *h_truth_jesdown = new TH1D("h_truth_jesdown", ";p_{T}^{Truth jet} [GeV]", truthnpt, truthptbins);
0185   TH1D *h_measure_jesdown = new TH1D("h_measure_jesdown", ";p_{T}^{Calib jet} [GeV]", calibnpt, calibptbins);
0186   TH2D *h_respmatrix_jesdown = new TH2D("h_respmatrix_jesdown", ";p_{T}^{Calib jet} [GeV];p_{T}^{Truth jet} [GeV]", calibnpt, calibptbins, truthnpt, truthptbins);
0187   TH1D *h_fake_jesdown = new TH1D("h_fake_jesdown", ";p_{T}^{Calib jet} [GeV]", calibnpt, calibptbins);
0188   TH1D *h_miss_jesdown = new TH1D("h_miss_jesdown", ";p_{T}^{Truth jet} [GeV]", truthnpt, truthptbins);
0189
0190   // JER up/down histograms
0191   TH1D *h_truth_jerup = new TH1D("h_truth_jerup", ";p_{T}^{Truth jet} [GeV]", truthnpt, truthptbins);
0192   TH1D *h_measure_jerup = new TH1D("h_measure_jerup", ";p_{T}^{Calib jet} [GeV]", calibnpt, calibptbins);
0193   TH2D *h_respmatrix_jerup = new TH2D("h_respmatrix_jerup", ";p_{T}^{Calib jet} [GeV];p_{T}^{Truth jet} [GeV]", calibnpt, calibptbins, truthnpt, truthptbins);
0194   TH1D *h_fake_jerup = new TH1D("h_fake_jerup", ";p_{T}^{Calib jet} [GeV]", calibnpt, calibptbins);
0195   TH1D *h_miss_jerup = new TH1D("h_miss_jerup", ";p_{T}^{Truth jet} [GeV]", truthnpt, truthptbins);
0196
0197   TH1D *h_truth_jerdown = new TH1D("h_truth_jerdown", ";p_{T}^{Truth jet} [GeV]", truthnpt, truthptbins);
0198   TH1D *h_measure_jerdown = new TH1D("h_measure_jerdown", ";p_{T}^{Calib jet} [GeV]", calibnpt, calibptbins);
0199   TH2D *h_respmatrix_jerdown = new TH2D("h_respmatrix_jerdown", ";p_{T}^{Calib jet} [GeV];p_{T}^{Truth jet} [GeV]", calibnpt, calibptbins, truthnpt, truthptbins);
0200   TH1D *h_fake_jerdown = new TH1D("h_fake_jerdown", ";p_{T}^{Calib jet} [GeV]", calibnpt, calibptbins);
0201   TH1D *h_miss_jerdown = new TH1D("h_miss_jerdown", ";p_{T}^{Truth jet} [GeV]", truthnpt, truthptbins);
0202
0203   // Closure test histograms
0204   TH1D *h_fullclosure_measure = new TH1D("h_fullclosure_measure", ";p_{T}^{Calib jet} [GeV]", calibnpt, calibptbins);
0205   TH1D *h_fullclosure_truth = new TH1D("h_fullclosure_truth", ";p_{T}^{Truth jet} [GeV]", truthnpt, truthptbins);
0206   TH2D *h_fullclosure_respmatrix = new TH2D("h_fullclosure_respmatrix", ";p_{T}^{Calib jet} [GeV];p_{T}^{Truth jet} [GeV]", calibnpt, calibptbins, truthnpt, truthptbins);
0207   TH1D *h_fullclosure_fake = new TH1D("h_fullclosure_fake", ";p_{T}^{Calib jet} [GeV]", calibnpt, calibptbins);
0208   TH1D *h_fullclosure_miss = new TH1D("h_fullclosure_miss", ";p_{T}^{Truth jet} [GeV]", truthnpt, truthptbins);
0209
0210   TH1D *h_halfclosure_measure = new TH1D("h_halfclosure_measure", ";p_{T}^{Calib jet} [GeV]", calibnpt, calibptbins);
0211   TH1D *h_halfclosure_truth = new TH1D("h_halfclosure_truth", ";p_{T}^{Truth jet} [GeV]", truthnpt, truthptbins);
0212   TH2D *h_halfclosure_respmatrix = new TH2D("h_halfclosure_respmatrix", ";p_{T}^{Calib jet} [GeV];p_{T}^{Truth jet} [GeV]", calibnpt, calibptbins, truthnpt , truthptbins);
0213   TH1D *h_halfclosure_fake = new TH1D("h_halfclosure_fake", ";p_{T}^{Calib jet} [GeV]", calibnpt, calibptbins);
0214   TH1D *h_halfclosure_miss = new TH1D("h_halfclosure_miss", ";p_{T}^{Truth jet} [GeV]", truthnpt, truthptbins);
0215
0216   TH1D *h_halfclosure_measure_check = new TH1D("h_halfclosure_measure_check", ";p_{T}^{Calib jet} [GeV]", calibnpt, calibptbins);
0217   TH1D *h_halfclosure_truth_check = new TH1D("h_halfclosure_truth_check", ";p_{T}^{Truth jet} [GeV]", truthnpt, truthptbins);
0218   TH2D *h_halfclosure_respmatrix_check = new TH2D("h_halfclosure_respmatrix_check", ";p_{T}^{Calib jet} [GeV];p_{T}^{Truth jet} [GeV]", calibnpt, calibptbins, truthnpt , truthptbins);
0219   TH1D *h_halfclosure_fake_check = new TH1D("h_halfclosure_fake_check", ";p_{T}^{Calib jet} [GeV]", calibnpt, calibptbins);
0220   TH1D *h_halfclosure_miss_check = new TH1D("h_halfclosure_miss_check", ";p_{T}^{Truth jet} [GeV]", truthnpt, truthptbins);
0221
0222   ////////// Event Loop //////////
0223   std::cout << "Data analysis started." << std::endl;
0224   int n_events = chain.GetEntries();
0225   //n_events = 10;
0226   std::cout << "Total number of events: " << n_events << std::endl;
0227   // Event variables setup.
0228   std::vector<bool> jet_filter, truthjet_filter;
0229   std::vector<float> goodtruthjet_pt, goodtruthjet_eta, goodtruthjet_phi, goodtruthjet_matched;
0230   std::vector<float> calibjet_pt_nominal, calibjet_eta_nominal, calibjet_phi_nominal, calibjet_matched_nominal;
0231   std::vector<float> calibjet_pt_jesup, calibjet_eta_jesup, calibjet_phi_jesup, calibjet_matched_jesup;
0232   std::vector<float> calibjet_pt_jesdown, calibjet_eta_jesdown, calibjet_phi_jesdown, calibjet_matched_jesdown;
0233   std::vector<float> calibjet_pt_jerup, calibjet_eta_jerup, calibjet_phi_jerup, calibjet_matched_jerup;
0234   std::vector<float> calibjet_pt_jerdown, calibjet_eta_jerdown, calibjet_phi_jerdown, calibjet_matched_jerdown;
0235   for (int ie = 0; ie < n_events; ++ie) { // event loop start
0236     // Load event.
0237     if (ie % 10000 == 0) {
0238       std::cout << "Processing event " << ie << "..." << std::endl;
0239     }
0240     chain.GetEntry(ie);
0241
0242     // Get z-vertex index, reweighting factor and re-assign no-zvertex events.
0243     int zvertex_index = get_zvertex_index(zvertex);
0244     double scale_zvertexreweight = h_zvertex_reweight->GetBinContent(h_zvertex_reweight->GetXaxis()->FindBin(zvertex));
0245     if (zvertex < -900) zvertex = 0;
0246
0247     // Do data set efficiency check.
0248     int leadingtruthjet_index = -1;
0249     float leadingtruthjet_pt = -9999;
0250     for (int ij = 0; ij < truthjet_pt->size(); ++ij) {
0251       if (truthjet_eta->at(ij) > 1.1 - jet_radius || truthjet_eta->at(ij) < -1.1 + jet_radius) continue;
0252       float jetpt = truthjet_pt->at(ij);
0253       if (jetpt > leadingtruthjet_pt) {
0254         leadingtruthjet_pt = jetpt;
0255         leadingtruthjet_index = ij;
0256       }
0257     }
0258     if (leadingtruthjet_index < 0) continue;
0259     if (leadingtruthjet_pt < truthjet_pt_min || leadingtruthjet_pt > truthjet_pt_max) continue;
0260
0261     // Event selection.
0262     bool pass_dijet = false;
0263     bool pass_njet = false;
0264     // Leading and subleading jet with energy and R = 0.4 for dijet and timing requirement.
0265     int leadingunsubjet04_index = -1;
0266     int subleadingunsubjet04_index = -1;
0267     get_leading_subleading_jet(leadingunsubjet04_index, subleadingunsubjet04_index, unsubjet04_e);
0268     if (leadingunsubjet04_index < 0) continue;
0269     // Dijet requirement.
0270     if (subleadingunsubjet04_index >= 0 && unsubjet04_e->at(subleadingunsubjet04_index) / (float) unsubjet04_e->at(leadingunsubjet04_index) > 0.3) {
0271       pass_dijet = match_leading_subleading_jet(unsubjet04_phi->at(leadingunsubjet04_index), unsubjet04_phi->at(subleadingunsubjet04_index));
0272     }
0273     // Njet requirement.
0274     int n_jets = 0;
0275     for (int ij = 0; ij < unsubjet04_e->size(); ++ij) {
0276       if (unsubjet04_e->at(ij) > 5.0) {
0277         n_jets++;
0278       }
0279     }
0280     if (n_jets <= 9) {
0281       pass_njet = true;
0282     }
0283
0284     // Do jet filter
0285     get_jet_filter(jet_filter, unsubjet_e, unsubjet_pt, unsubjet_eta, zvertex, jet_radius);
0286     get_jet_filter(truthjet_filter, truthjet_e, truthjet_pt, truthjet_eta, zvertex, jet_radius);
0287     for (int ij = 0; ij < unsubjet_pt->size(); ++ij) {
0288       if (jet_filter[ij]) continue;
0289       h_recojet_pt_record->Fill(unsubjet_pt->at(ij));
0290       if (pass_dijet && pass_njet) h_recojet_pt_passcut_record->Fill(unsubjet_pt->at(ij));
0291     }
0292     for (int ij = 0; ij < truthjet_pt->size(); ++ij) {
0293       if (truthjet_filter[ij]) continue;
0294       h_truthjet_pt_record->Fill(truthjet_pt->at(ij));
0295       if (pass_dijet && pass_njet) h_truthjet_pt_passcut_record->Fill(truthjet_pt->at(ij));
0296     }
0297
0298     // Get truth jet and calib jet.
0299     get_truthjet(goodtruthjet_pt, goodtruthjet_eta, goodtruthjet_phi, truthjet_filter, truthjet_pt, truthjet_eta, truthjet_phi);
0300
0301     get_calibjet(calibjet_pt_nominal, calibjet_eta_nominal, calibjet_phi_nominal, jet_filter, unsubjet_pt, unsubjet_eta, unsubjet_phi, f_corr, zvertex_index, jet_radius, 1, 0.1, (!pass_dijet||!pass_njet));
0302     get_calibjet(calibjet_pt_jesup, calibjet_eta_jesup, calibjet_phi_jesup, jet_filter, unsubjet_pt, unsubjet_eta, unsubjet_phi, f_corr, zvertex_index, jet_radius, 1.06, 0.1, (!pass_dijet||!pass_njet));
0303     get_calibjet(calibjet_pt_jesdown, calibjet_eta_jesdown, calibjet_phi_jesdown, jet_filter, unsubjet_pt, unsubjet_eta, unsubjet_phi, f_corr, zvertex_index, jet_radius, 0.94, 0.1, (!pass_dijet||!pass_njet));
0304     get_calibjet(calibjet_pt_jerup, calibjet_eta_jerup, calibjet_phi_jerup, jet_filter, unsubjet_pt, unsubjet_eta, unsubjet_phi, f_corr, zvertex_index, jet_radius, 1, 0.13, (!pass_dijet||!pass_njet));
0305     get_calibjet(calibjet_pt_jerdown, calibjet_eta_jerdown, calibjet_phi_jerdown, jet_filter, unsubjet_pt, unsubjet_eta, unsubjet_phi, f_corr, zvertex_index, jet_radius, 1, 0.07, (!pass_dijet||!pass_njet));
0306
0307     // Match truth jet and calib jet.
0308     match_meas_truth(calibjet_eta_nominal, calibjet_phi_nominal, calibjet_matched_nominal, goodtruthjet_eta, goodtruthjet_phi, goodtruthjet_matched, jet_radius);
0309     fill_response_matrix(h_truth_nominal, h_measure_nominal, h_respmatrix_nominal, h_fake_nominal, h_miss_nominal, scale_zvertexreweight, calibjet_pt_nominal, calibjet_matched_nominal, goodtruthjet_pt, goodtruthjet_matched);
0310
0311     fill_response_matrix(h_fullclosure_truth, h_fullclosure_measure, h_fullclosure_respmatrix, h_fullclosure_fake, h_fullclosure_miss, scale_zvertexreweight, calibjet_pt_nominal, calibjet_matched_nominal, goodtruthjet_pt, goodtruthjet_matched);
0312     if (ie % 2 == 0) {
0313       fill_response_matrix(h_halfclosure_truth, h_halfclosure_measure, h_halfclosure_respmatrix, h_halfclosure_fake, h_halfclosure_miss, scale_zvertexreweight, calibjet_pt_nominal, calibjet_matched_nominal, goodtruthjet_pt, goodtruthjet_matched);
0314     } else {
0315       fill_response_matrix(h_halfclosure_truth_check, h_halfclosure_measure_check, h_halfclosure_respmatrix_check, h_halfclosure_fake_check, h_halfclosure_miss_check, scale_zvertexreweight, calibjet_pt_nominal, calibjet_matched_nominal, goodtruthjet_pt, goodtruthjet_matched);
0316     }
0317
0318     match_meas_truth(calibjet_eta_jesup, calibjet_phi_jesup, calibjet_matched_jesup, goodtruthjet_eta, goodtruthjet_phi, goodtruthjet_matched, jet_radius);
0319     fill_response_matrix(h_truth_jesup, h_measure_jesup, h_respmatrix_jesup, h_fake_jesup, h_miss_jesup, scale_zvertexreweight, calibjet_pt_jesup, calibjet_matched_jesup, goodtruthjet_pt, goodtruthjet_matched);
0320
0321     match_meas_truth(calibjet_eta_jesdown, calibjet_phi_jesdown, calibjet_matched_jesdown, goodtruthjet_eta, goodtruthjet_phi, goodtruthjet_matched, jet_radius);
0322     fill_response_matrix(h_truth_jesdown, h_measure_jesdown, h_respmatrix_jesdown, h_fake_jesdown, h_miss_jesdown, scale_zvertexreweight, calibjet_pt_jesdown, calibjet_matched_jesdown, goodtruthjet_pt, goodtruthjet_matched);
0323
0324     match_meas_truth(calibjet_eta_jerup, calibjet_phi_jerup, calibjet_matched_jerup, goodtruthjet_eta, goodtruthjet_phi, goodtruthjet_matched, jet_radius);
0325     fill_response_matrix(h_truth_jerup, h_measure_jerup, h_respmatrix_jerup, h_fake_jerup, h_miss_jerup, scale_zvertexreweight, calibjet_pt_jerup, calibjet_matched_jerup, goodtruthjet_pt, goodtruthjet_matched);
0326
0327     match_meas_truth(calibjet_eta_jerdown, calibjet_phi_jerdown, calibjet_matched_jerdown, goodtruthjet_eta, goodtruthjet_phi, goodtruthjet_matched, jet_radius);
0328     fill_response_matrix(h_truth_jerdown, h_measure_jerdown, h_respmatrix_jerdown, h_fake_jerdown, h_miss_jerdown, scale_zvertexreweight, calibjet_pt_jerdown, calibjet_matched_jerdown, goodtruthjet_pt, goodtruthjet_matched);
0329   } // event loop end
0330   // Fill event histograms.
0331   h_event->Fill(0.5, n_events);
0332
0333   // Write histograms.
0334   std::cout << "Writing histograms..." << std::endl;
0335   f_out->cd();
0336   h_event->Write();
0337   h_recojet_pt_record->Write();
0338   h_recojet_pt_passcut_record->Write();
0339   h_truthjet_pt_record->Write();
0340   h_truthjet_pt_passcut_record->Write();
0341
0342   h_truth_nominal->Write();
0343   h_measure_nominal->Write();
0344   h_respmatrix_nominal->Write();
0345   h_fake_nominal->Write();
0346   h_miss_nominal->Write();
0347
0348   h_truth_jesup->Write();
0349   h_measure_jesup->Write();
0350   h_respmatrix_jesup->Write();
0351   h_fake_jesup->Write();
0352   h_miss_jesup->Write();
0353
0354   h_truth_jesdown->Write();
0355   h_measure_jesdown->Write();
0356   h_respmatrix_jesdown->Write();
0357   h_fake_jesdown->Write();
0358   h_miss_jesdown->Write();
0359
0360   h_truth_jerup->Write();
0361   h_measure_jerup->Write();
0362   h_respmatrix_jerup->Write();
0363   h_fake_jerup->Write();
0364   h_miss_jerup->Write();
0365
0366   h_truth_jerdown->Write();
0367   h_measure_jerdown->Write();
0368   h_respmatrix_jerdown->Write();
0369   h_fake_jerdown->Write();
0370   h_miss_jerdown->Write();
0371
0372   h_fullclosure_measure->Write();
0373   h_fullclosure_truth->Write();
0374   h_fullclosure_respmatrix->Write();
0375   h_fullclosure_fake->Write();
0376   h_fullclosure_miss->Write();
0377
0378   h_halfclosure_measure->Write();
0379   h_halfclosure_truth->Write();
0380   h_halfclosure_respmatrix->Write();
0381   h_halfclosure_fake->Write();
0382   h_halfclosure_miss->Write();
0383
0384   h_halfclosure_measure_check->Write();
0385   h_halfclosure_truth_check->Write();
0386   h_halfclosure_respmatrix_check->Write();
0387   h_halfclosure_fake_check->Write();
0388   h_halfclosure_miss_check->Write();
0389
0390   f_out->Close();
0391   std::cout << "All done!" << std::endl;
0392 }
0393
0394 ////////////////////////////////////////// Helper Functions //////////////////////////////////////////
0395 float get_deta(float eta1, float eta2) {
0396   return eta1 - eta2;
0397 }
0398
0399 float get_dphi(float phi1, float phi2) {
0400   float dphi1 = phi1 - phi2;
0401   float dphi2 = phi1 - phi2 + 2*TMath::Pi();
0402   float dphi3 = phi1 - phi2 - 2*TMath::Pi();
0403   if (fabs(dphi1) > fabs(dphi2)) {
0404     dphi1 = dphi2;
0405   }
0406   if (fabs(dphi1) > fabs(dphi3)) {
0407     dphi1 = dphi3;
0408   }
0409   return dphi1;
0410 }
0411
0412 float get_dR(float eta1, float phi1, float eta2, float phi2) {
0413   float deta = get_deta(eta1, eta2);
0414   float dphi = get_dphi(phi1, phi2);
0415   return sqrt(deta*deta + dphi*dphi);
0416 }
0417
0418 float get_emcal_mineta_zcorrected(float zvertex) {
0419   float minz_EM = -130.23;
0420   float radius_EM = 93.5;
0421   float z = minz_EM - zvertex;
0422   float eta_zcorrected = asinh(z / (float)radius_EM);
0423   return eta_zcorrected;
0424 }
0425
0426 float get_emcal_maxeta_zcorrected(float zvertex) {
0427   float maxz_EM = 130.23;
0428   float radius_EM = 93.5;
0429   float z = maxz_EM - zvertex;
0430   float eta_zcorrected = asinh(z / (float)radius_EM);
0431   return eta_zcorrected;
0432 }
0433
0434 float get_ihcal_mineta_zcorrected(float zvertex) {
0435   float minz_IH = -170.299;
0436   float radius_IH = 127.503;
0437   float z = minz_IH - zvertex;
0438   float eta_zcorrected = asinh(z / (float)radius_IH);
0439   return eta_zcorrected;
0440 }
0441
0442 float get_ihcal_maxeta_zcorrected(float zvertex) {
0443   float maxz_IH = 170.299;
0444   float radius_IH = 127.503;
0445   float z = maxz_IH - zvertex;
0446   float eta_zcorrected = asinh(z / (float)radius_IH);
0447   return eta_zcorrected;
0448 }
0449
0450 float get_ohcal_mineta_zcorrected(float zvertex) {
0451   float minz_OH = -301.683;
0452   float radius_OH = 225.87;
0453   float z = minz_OH - zvertex;
0454   float eta_zcorrected = asinh(z / (float)radius_OH);
0455   return eta_zcorrected;
0456 }
0457
0458 float get_ohcal_maxeta_zcorrected(float zvertex) {
0459   float maxz_OH = 301.683;
0460   float radius_OH = 225.87;
0461   float z = maxz_OH - zvertex;
0462   float eta_zcorrected = asinh(z / (float)radius_OH);
0463   return eta_zcorrected;
0464 }
0465
0466 bool check_bad_jet_eta(float jet_eta, float zvertex, float jet_radius) {
0467   float emcal_mineta = get_emcal_mineta_zcorrected(zvertex);
0468   float emcal_maxeta = get_emcal_maxeta_zcorrected(zvertex);
0469   float ihcal_mineta = get_ihcal_mineta_zcorrected(zvertex);
0470   float ihcal_maxeta = get_ihcal_maxeta_zcorrected(zvertex);
0471   float ohcal_mineta = get_ohcal_mineta_zcorrected(zvertex);
0472   float ohcal_maxeta = get_ohcal_maxeta_zcorrected(zvertex);
0473   float minlimit = emcal_mineta;
0474   if (ihcal_mineta > minlimit) minlimit = ihcal_mineta;
0475   if (ohcal_mineta > minlimit) minlimit = ohcal_mineta;
0476   float maxlimit = emcal_maxeta;
0477   if (ihcal_maxeta < maxlimit) maxlimit = ihcal_maxeta;
0478   if (ohcal_maxeta < maxlimit) maxlimit = ohcal_maxeta;
0479   minlimit += jet_radius;
0480   maxlimit -= jet_radius;
0481   return jet_eta < minlimit || jet_eta > maxlimit;
0482 }
0483
0484 int get_zvertex_index(float zvertex) {
0485   if (zvertex >= -60.0 && zvertex < -30.0) return 1;
0486   else if (zvertex >= -30.0 && zvertex < 30.0) return 0;
0487   else if (zvertex >= 30.0 && zvertex < 60) return 2;
0488   else if ((zvertex < -60.0 && zvertex >= -900) || zvertex >= 60.0) return 3;
0489   else if (zvertex < -900) return 4;
0490   return -1; // should not reach here
0491 }
0492
0493 int get_eta_index(float jet_eta, int zvertex_index, double jet_radius) {
0494   float eta_low = 0;
0495   float eta_high = 0;
0496   if (zvertex_index == 0) {// -30 < zvrtx < 30
0497     eta_low = -1.2 + jet_radius;
0498     eta_high = 1.2 - jet_radius;
0499   } else if (zvertex_index == 1) {// -60 < zvrtx < -30
0500     eta_low = -0.95 + jet_radius;
0501     eta_high = 1.25 - jet_radius;
0502   } else if (zvertex_index == 2) {// 30 < zvrtx < inf
0503     eta_low = -1.25 + jet_radius;
0504     eta_high = 0.95 - jet_radius;
0505   } else if (zvertex_index == 3 || zvertex_index == 4) {
0506     return 0;
0507   }
0508   float threshold1 = eta_low + ((eta_high - eta_low) / 4.0);
0509   float threshold2 = eta_low + ((eta_high - eta_low) / 2.0);
0510   float threshold3 = eta_low + (3.0 * (eta_high - eta_low) / 4.0);
0511   if (jet_eta < threshold1) return 0;  // -inf to threshold1
0512   if (jet_eta < threshold2) return 1;  // threshold1 to threshold2
0513   if (jet_eta < threshold3) return 2;  // threshold2 to threshold3
0514   return 3;  // threshold3 to inf
0515 }
0516
0517 ////////////////////////////////////////// Functions //////////////////////////////////////////
0518 void get_leading_subleading_jet(int& leadingjet_index, int& subleadingjet_index, std::vector<float>* jet_et) {
0519   leadingjet_index = -1;
0520   subleadingjet_index = -1;
0521   float leadingjet_et = -9999;
0522   float subleadingjet_et = -9999;
0523   for (int ij = 0; ij < jet_et->size(); ++ij) {
0524     float jetet = jet_et->at(ij);
0525     if (jetet > leadingjet_et) {
0526       subleadingjet_et = leadingjet_et;
0527       subleadingjet_index = leadingjet_index;
0528       leadingjet_et = jetet;
0529       leadingjet_index = ij;
0530     } else if (jetet > subleadingjet_et) {
0531       subleadingjet_et = jetet;
0532       subleadingjet_index = ij;
0533     }
0534   }
0535 }
0536
0537 bool match_leading_subleading_jet(float leadingjet_phi, float subleadingjet_phi) {
0538   float dijet_min_phi = 3*TMath::Pi()/4.;
0539   float dphi = fabs(get_dphi(leadingjet_phi, subleadingjet_phi));
0540   return dphi > dijet_min_phi;
0541 }
0542
0543 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) {
0544   jet_filter.clear();
0545   int njet = jet_e->size();
0546   for (int ij = 0; ij < njet; ++ij) {
0547     jet_filter.push_back(jet_e->at(ij) < 0 || check_bad_jet_eta(jet_eta->at(ij), zvertex, jet_radius));
0548   }
0549 }
0550
0551 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, float jes_para, float jer_para, bool jet_background) {
0552   calibjet_pt.clear();
0553   calibjet_eta.clear();
0554   calibjet_phi.clear();
0555   for (int ij = 0; ij < jet_filter.size(); ++ij) {
0556     if (jet_filter.at(ij) || jet_background) continue;
0557     int eta_index;
0558     if (zvertex_index == 0 || zvertex_index == 1 || zvertex_index == 2) {
0559       eta_index = get_eta_index(jet_eta->at(ij), zvertex_index, jet_radius);
0560     } else {
0561       eta_index = 0;
0562     }
0563     double calib_pt = f_corr[zvertex_index][eta_index]->Eval(jet_pt->at(ij)) * (1 + randGen.Gaus(0.0, jer_para)) * jes_para;
0564     //std::cout << "Debug: jet_pt = " << jet_pt->at(ij) << ", calib_pt = " << calib_pt << " with zvertex = " << zvertex_index << ", eta = " << jet_eta->at(ij) << std::endl;
0565     if (calib_pt < calibptbins[0] || calib_pt > calibptbins[calibnpt]) continue;
0566     calibjet_pt.push_back(calib_pt);
0567     calibjet_eta.push_back(jet_eta->at(ij));
0568     calibjet_phi.push_back(jet_phi->at(ij));
0569   }
0570 }
0571
0572 void get_truthjet(std::vector<float>& goodtruthjet_pt, std::vector<float>& goodtruthjet_eta, std::vector<float>& goodtruthjet_phi, std::vector<bool>& truthjet_filter, std::vector<float>* jet_pt, std::vector<float>* jet_eta, std::vector<float>* jet_phi) {
0573   goodtruthjet_pt.clear();
0574   goodtruthjet_eta.clear();
0575   goodtruthjet_phi.clear();
0576   for (int ij = 0; ij < truthjet_filter.size(); ++ij) {
0577     if (truthjet_filter.at(ij)) continue;
0578     if (jet_pt->at(ij) < truthptbins[0] || jet_pt->at(ij) > truthptbins[truthnpt]) continue;
0579     goodtruthjet_pt.push_back(jet_pt->at(ij));
0580     goodtruthjet_eta.push_back(jet_eta->at(ij));
0581     goodtruthjet_phi.push_back(jet_phi->at(ij));
0582   }
0583 }
0584
0585 void match_meas_truth(std::vector<float>& meas_eta, std::vector<float>& meas_phi, std::vector<float>& meas_matched, std::vector<float>& truth_eta, std::vector<float>& truth_phi, std::vector<float>& truth_matched, float jet_radius) {
0586   meas_matched.assign(meas_eta.size(), -1);
0587   truth_matched.assign(truth_eta.size(), -1);
0588   float max_match_dR = jet_radius * 0.75;
0589   for (int im = 0; im < meas_eta.size(); ++im) {
0590     float min_dR = 100;
0591     int match_index = -9999;
0592     for (int it = 0; it < truth_eta.size(); ++it) {
0593       float dR = get_dR(meas_eta[im], meas_phi[im], truth_eta[it], truth_phi[it]);
0594       if (dR < min_dR) {
0595         match_index = it;
0596         min_dR = dR;
0597       }
0598     }
0599     if (min_dR < max_match_dR) {
0600       if (truth_matched[match_index] == -1) {
0601         meas_matched[im] = match_index;
0602         truth_matched[match_index] = im;
0603       } else {
0604         float dR1 = get_dR(meas_eta[truth_matched[match_index]], meas_phi[truth_matched[match_index]], truth_eta[match_index], truth_phi[match_index]);
0605         if (min_dR < dR1) {
0606           meas_matched[truth_matched[match_index]] = -1;
0607           meas_matched[im] = match_index;
0608           truth_matched[match_index] = im;
0609         }
0610       }
0611     }
0612   }
0613 }
0614
0615 void fill_response_matrix(TH1D*& h_truth, TH1D*& h_meas, TH2D*& h_resp, TH1D*& h_fake, TH1D*& h_miss, double scale, std::vector<float>& meas_pt, std::vector<float>& meas_matched, std::vector<float>& truth_pt, std::vector<float>& truth_matched) {
0616   for (int im = 0; im < meas_pt.size(); ++im) {
0617     h_meas->Fill(meas_pt[im], scale);
0618     if (meas_matched[im] < 0) {
0619       h_fake->Fill(meas_pt[im], scale);
0620     } else {
0621       h_resp->Fill(meas_pt[im], truth_pt[meas_matched[im]], scale);
0622     }
0623   }
0624   for (int it = 0; it < truth_pt.size(); ++it) {
0625     h_truth->Fill(truth_pt[it], scale);
0626     if (truth_matched[it] < 0) {
0627       h_miss->Fill(truth_pt[it], scale);
0628     }
0629   }
0630 }