JetSpectrum/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_jet(int& leadingjet_index, std::vector<float>* jet_pt);
0016 void get_leading_subleading_jet(int& leadingjet_index, int& subleadingjet_index, std::vector<float>* jet_et);
0017 bool match_leading_subleading_jet(float leadingjet_phi, float subleadingjet_phi);
0018 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);
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, std::vector<bool>* jet_background, TF1* f_corr, float jes_para, float jer_para);
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, std::vector<float>& meas_pt, std::vector<float>& meas_matched, std::vector<float>& truth_pt, std::vector<float>& truth_matched, float weight_scale);
0023
0024 TF1* f_jer = new TF1("f_jer", "sqrt( 0.095077098*0.095077098 + (0.63134847*0.63134847/x) + (2.1664610*2.1664610/(x*x)) )", 0, 100);
0025 TRandom3 randGen(1234);
0026
0027 ////////////////////////////////////////// Main Function //////////////////////////////////////////
0028 void analysis_sim(std::string runtype, int nseg, int iseg)  {
0029   ////////// General Set up //////////
0030   double weight_scale = 1.0, truthjet_pt_min = 0, truthjet_pt_max = 3000;
0031   if (runtype == "MB") {
0032     weight_scale = MB_scale;
0033     truthjet_pt_min = 0;
0034     truthjet_pt_max = 14;
0035   } else if (runtype == "Jet10GeV") {
0036     weight_scale = Jet10GeV_scale;
0037     truthjet_pt_min = 14;
0038     truthjet_pt_max = 30;
0039   } else if (runtype == "Jet30GeV") {
0040     weight_scale = Jet30GeV_scale;
0041     truthjet_pt_min = 30;
0042     truthjet_pt_max = 3000;
0043   } else {
0044     std::cout << "Unknown runtype" << std::endl;
0045     return;
0046   }
0047
0048   ////////// Files //////////
0049   TFile *f_out = new TFile(Form("output_sim/output_%s_%d_%d.root", runtype.c_str(), iseg, iseg+nseg), "RECREATE");
0050   TChain chain("ttree");
0051   for (int i = iseg; i < iseg + nseg; ++i) {
0052     chain.Add(Form("/sphenix/tg/tg01/jets/hanpuj/JES_MC_run21/%s/OutDir%d/output_sim.root", runtype.c_str(), i));
0053   }
0054   chain.SetBranchStatus("*", 0);
0055
0056   float zvertex; chain.SetBranchStatus("z_vertex", 1); chain.SetBranchAddress("z_vertex", &zvertex);
0057   std::vector<float>* unsubjet_e = nullptr; chain.SetBranchStatus("unsubjet04_e", 1); chain.SetBranchAddress("unsubjet04_e", &unsubjet_e);
0058   std::vector<float>* unsubjet_pt = nullptr; chain.SetBranchStatus("unsubjet04_pt", 1); chain.SetBranchAddress("unsubjet04_pt", &unsubjet_pt);
0059   std::vector<float>* unsubjet_et = nullptr; chain.SetBranchStatus("unsubjet04_et", 1); chain.SetBranchAddress("unsubjet04_et", &unsubjet_et);
0060   std::vector<float>* unsubjet_eta = nullptr; chain.SetBranchStatus("unsubjet04_eta", 1); chain.SetBranchAddress("unsubjet04_eta", &unsubjet_eta);
0061   std::vector<float>* unsubjet_phi = nullptr; chain.SetBranchStatus("unsubjet04_phi", 1); chain.SetBranchAddress("unsubjet04_phi", &unsubjet_phi);
0062   std::vector<float>* unsubjet_emcal_calo_e = nullptr; chain.SetBranchStatus("unsubjet04_emcal_calo_e", 1); chain.SetBranchAddress("unsubjet04_emcal_calo_e", &unsubjet_emcal_calo_e);
0063   std::vector<float>* unsubjet_ihcal_calo_e = nullptr; chain.SetBranchStatus("unsubjet04_ihcal_calo_e", 1); chain.SetBranchAddress("unsubjet04_ihcal_calo_e", &unsubjet_ihcal_calo_e);
0064   std::vector<float>* unsubjet_ohcal_calo_e = nullptr; chain.SetBranchStatus("unsubjet04_ohcal_calo_e", 1); chain.SetBranchAddress("unsubjet04_ohcal_calo_e", &unsubjet_ohcal_calo_e);
0065
0066   std::vector<float>* truthjet_e = nullptr; chain.SetBranchStatus("truthjet04_e", 1); chain.SetBranchAddress("truthjet04_e", &truthjet_e);
0067   std::vector<float>* truthjet_pt = nullptr; chain.SetBranchStatus("truthjet04_pt", 1); chain.SetBranchAddress("truthjet04_pt", &truthjet_pt);
0068   std::vector<float>* truthjet_eta = nullptr; chain.SetBranchStatus("truthjet04_eta", 1); chain.SetBranchAddress("truthjet04_eta", &truthjet_eta);
0069   std::vector<float>* truthjet_phi = nullptr; chain.SetBranchStatus("truthjet04_phi", 1); chain.SetBranchAddress("truthjet04_phi", &truthjet_phi);
0070
0071   ////////// JES func //////////
0072   TFile *corrFile = new TFile("/sphenix/user/hanpuj/JES_MC_Calibration/offline/JES_Calib_Default.root", "READ");
0073   if (!corrFile) {
0074     std::cout << "Error: cannot open JES_Calib_Default.root" << std::endl;
0075     return;
0076   }
0077   TF1 *f_corr = (TF1*)corrFile->Get("JES_Calib_Default_Func");
0078   if (!f_corr) {
0079     std::cout << "Error: cannot open f_corr" << std::endl;
0080     return;
0081   }
0082
0083   ////////// Histograms //////////
0084   TH1D* h_zvertex = new TH1D("h_zvertex", ";Z-vertex [cm]", 60, -30, 30);
0085   TH1D *h_recojet_pt_record = new TH1D("h_recojet_pt_record", ";p_{T} [GeV]", 1000, 0, 100);
0086   TH1D *h_recojet_pt_record_dijet = new TH1D("h_recojet_pt_record_dijet", ";p_{T} [GeV]", 1000, 0, 100);
0087   TH1D *h_recojet_pt_record_frac = new TH1D("h_recojet_pt_record_frac", ";p_{T} [GeV]", 1000, 0, 100);
0088
0089   TH1D* h_truth_calib_dijet = new TH1D("h_truth_calib_dijet", ";p_{T}^{Truth jet} [GeV]", truthnpt, truthptbins);
0090   TH1D* h_measure_calib_dijet = new TH1D("h_measure_calib_dijet", ";p_{T}^{Calib jet} [GeV]", calibnpt, calibptbins);
0091   TH2D* h_respmatrix_calib_dijet = new TH2D("h_respmatrix_calib_dijet", ";p_{T}^{Calib jet} [GeV];p_{T}^{Truth jet} [GeV]", calibnpt, calibptbins, truthnpt, truthptbins);
0092   TH1D* h_fake_calib_dijet = new TH1D("h_fake_calib_dijet", ";p_{T}^{Calib jet} [GeV]", calibnpt, calibptbins);
0093   TH1D* h_miss_calib_dijet = new TH1D("h_miss_calib_dijet", ";p_{T}^{Truth jet} [GeV]", truthnpt, truthptbins);
0094
0095   TH1D* h_truth_calib_dijet_jesdown = new TH1D("h_truth_calib_dijet_jesdown", ";p_{T}^{Truth jet} [GeV]", truthnpt, truthptbins);
0096   TH1D* h_measure_calib_dijet_jesdown = new TH1D("h_measure_calib_dijet_jesdown", ";p_{T}^{Calib jet} [GeV]", calibnpt, calibptbins);
0097   TH2D* h_respmatrix_calib_dijet_jesdown = new TH2D("h_respmatrix_calib_dijet_jesdown", ";p_{T}^{Calib jet} [GeV];p_{T}^{Truth jet} [GeV]", calibnpt, calibptbins, truthnpt, truthptbins);
0098   TH1D* h_fake_calib_dijet_jesdown = new TH1D("h_fake_calib_dijet_jesdown", ";p_{T}^{Calib jet} [GeV]", calibnpt, calibptbins);
0099   TH1D* h_miss_calib_dijet_jesdown = new TH1D("h_miss_calib_dijet_jesdown", ";p_{T}^{Truth jet} [GeV]", truthnpt, truthptbins);
0100
0101   TH1D* h_truth_calib_dijet_jesup = new TH1D("h_truth_calib_dijet_jesup", ";p_{T}^{Truth jet} [GeV]", truthnpt, truthptbins);
0102   TH1D* h_measure_calib_dijet_jesup = new TH1D("h_measure_calib_dijet_jesup", ";p_{T}^{Calib jet} [GeV]", calibnpt, calibptbins);
0103   TH2D* h_respmatrix_calib_dijet_jesup = new TH2D("h_respmatrix_calib_dijet_jesup", ";p_{T}^{Calib jet} [GeV];p_{T}^{Truth jet} [GeV]", calibnpt, calibptbins, truthnpt, truthptbins);
0104   TH1D* h_fake_calib_dijet_jesup = new TH1D("h_fake_calib_dijet_jesup", ";p_{T}^{Calib jet} [GeV]", calibnpt, calibptbins);
0105   TH1D* h_miss_calib_dijet_jesup = new TH1D("h_miss_calib_dijet_jesup", ";p_{T}^{Truth jet} [GeV]", truthnpt, truthptbins);
0106
0107   TH1D* h_truth_calib_dijet_jerdown = new TH1D("h_truth_calib_dijet_jerdown", ";p_{T}^{Truth jet} [GeV]", truthnpt, truthptbins);
0108   TH1D* h_measure_calib_dijet_jerdown = new TH1D("h_measure_calib_dijet_jerdown", ";p_{T}^{Calib jet} [GeV]", calibnpt, calibptbins);
0109   TH2D* h_respmatrix_calib_dijet_jerdown = new TH2D("h_respmatrix_calib_dijet_jerdown", ";p_{T}^{Calib jet} [GeV];p_{T}^{Truth jet} [GeV]", calibnpt, calibptbins, truthnpt, truthptbins);
0110   TH1D* h_fake_calib_dijet_jerdown = new TH1D("h_fake_calib_dijet_jerdown", ";p_{T}^{Calib jet} [GeV]", calibnpt, calibptbins);
0111   TH1D* h_miss_calib_dijet_jerdown = new TH1D("h_miss_calib_dijet_jerdown", ";p_{T}^{Truth jet} [GeV]", truthnpt, truthptbins);
0112
0113   TH1D* h_truth_calib_dijet_jerup = new TH1D("h_truth_calib_dijet_jerup", ";p_{T}^{Truth jet} [GeV]", truthnpt, truthptbins);
0114   TH1D* h_measure_calib_dijet_jerup = new TH1D("h_measure_calib_dijet_jerup", ";p_{T}^{Calib jet} [GeV]", calibnpt, calibptbins);
0115   TH2D* h_respmatrix_calib_dijet_jerup = new TH2D("h_respmatrix_calib_dijet_jerup", ";p_{T}^{Calib jet} [GeV];p_{T}^{Truth jet} [GeV]", calibnpt, calibptbins, truthnpt, truthptbins);
0116   TH1D* h_fake_calib_dijet_jerup = new TH1D("h_fake_calib_dijet_jerup", ";p_{T}^{Calib jet} [GeV]", calibnpt, calibptbins);
0117   TH1D* h_miss_calib_dijet_jerup = new TH1D("h_miss_calib_dijet_jerup", ";p_{T}^{Truth jet} [GeV]", truthnpt, truthptbins);
0118
0119   TH1D* h_truth_calib_frac = new TH1D("h_truth_calib_frac", ";p_{T}^{Truth jet} [GeV]", truthnpt, truthptbins);
0120   TH1D* h_measure_calib_frac = new TH1D("h_measure_calib_frac", ";p_{T}^{Calib jet} [GeV]", calibnpt, calibptbins);
0121   TH2D* h_respmatrix_calib_frac = new TH2D("h_respmatrix_calib_frac", ";p_{T}^{Calib jet} [GeV];p_{T}^{Truth jet} [GeV]", calibnpt, calibptbins, truthnpt, truthptbins);
0122   TH1D* h_fake_calib_frac = new TH1D("h_fake_calib_frac", ";p_{T}^{Calib jet} [GeV]", calibnpt, calibptbins);
0123   TH1D* h_miss_calib_frac = new TH1D("h_miss_calib_frac", ";p_{T}^{Truth jet} [GeV]", truthnpt, truthptbins);
0124
0125   TH1D* h_truth_calib_frac_jesdown = new TH1D("h_truth_calib_frac_jesdown", ";p_{T}^{Truth jet} [GeV]", truthnpt, truthptbins);
0126   TH1D* h_measure_calib_frac_jesdown = new TH1D("h_measure_calib_frac_jesdown", ";p_{T}^{Calib jet} [GeV]", calibnpt, calibptbins);
0127   TH2D* h_respmatrix_calib_frac_jesdown = new TH2D("h_respmatrix_calib_frac_jesdown", ";p_{T}^{Calib jet} [GeV];p_{T}^{Truth jet} [GeV]", calibnpt, calibptbins, truthnpt, truthptbins);
0128   TH1D* h_fake_calib_frac_jesdown = new TH1D("h_fake_calib_frac_jesdown", ";p_{T}^{Calib jet} [GeV]", calibnpt, calibptbins);
0129   TH1D* h_miss_calib_frac_jesdown = new TH1D("h_miss_calib_frac_jesdown", ";p_{T}^{Truth jet} [GeV]", truthnpt, truthptbins);
0130
0131   TH1D* h_truth_calib_frac_jesup = new TH1D("h_truth_calib_frac_jesup", ";p_{T}^{Truth jet} [GeV]", truthnpt, truthptbins);
0132   TH1D* h_measure_calib_frac_jesup = new TH1D("h_measure_calib_frac_jesup", ";p_{T}^{Calib jet} [GeV]", calibnpt, calibptbins);
0133   TH2D* h_respmatrix_calib_frac_jesup = new TH2D("h_respmatrix_calib_frac_jesup", ";p_{T}^{Calib jet} [GeV];p_{T}^{Truth jet} [GeV]", calibnpt, calibptbins, truthnpt, truthptbins);
0134   TH1D* h_fake_calib_frac_jesup = new TH1D("h_fake_calib_frac_jesup", ";p_{T}^{Calib jet} [GeV]", calibnpt, calibptbins);
0135   TH1D* h_miss_calib_frac_jesup = new TH1D("h_miss_calib_frac_jesup", ";p_{T}^{Truth jet} [GeV]", truthnpt, truthptbins);
0136
0137   TH1D* h_truth_calib_frac_jerdown = new TH1D("h_truth_calib_frac_jerdown", ";p_{T}^{Truth jet} [GeV]", truthnpt, truthptbins);
0138   TH1D* h_measure_calib_frac_jerdown = new TH1D("h_measure_calib_frac_jerdown", ";p_{T}^{Calib jet} [GeV]", calibnpt, calibptbins);
0139   TH2D* h_respmatrix_calib_frac_jerdown = new TH2D("h_respmatrix_calib_frac_jerdown", ";p_{T}^{Calib jet} [GeV];p_{T}^{Truth jet} [GeV]", calibnpt, calibptbins, truthnpt, truthptbins);
0140   TH1D* h_fake_calib_frac_jerdown = new TH1D("h_fake_calib_frac_jerdown", ";p_{T}^{Calib jet} [GeV]", calibnpt, calibptbins);
0141   TH1D* h_miss_calib_frac_jerdown = new TH1D("h_miss_calib_frac_jerdown", ";p_{T}^{Truth jet} [GeV]", truthnpt, truthptbins);
0142
0143   TH1D* h_truth_calib_frac_jerup = new TH1D("h_truth_calib_frac_jerup", ";p_{T}^{Truth jet} [GeV]", truthnpt, truthptbins);
0144   TH1D* h_measure_calib_frac_jerup = new TH1D("h_measure_calib_frac_jerup", ";p_{T}^{Calib jet} [GeV]", calibnpt, calibptbins);
0145   TH2D* h_respmatrix_calib_frac_jerup = new TH2D("h_respmatrix_calib_frac_jerup", ";p_{T}^{Calib jet} [GeV];p_{T}^{Truth jet} [GeV]", calibnpt, calibptbins, truthnpt, truthptbins);
0146   TH1D* h_fake_calib_frac_jerup = new TH1D("h_fake_calib_frac_jerup", ";p_{T}^{Calib jet} [GeV]", calibnpt, calibptbins);
0147   TH1D* h_miss_calib_frac_jerup = new TH1D("h_miss_calib_frac_jerup", ";p_{T}^{Truth jet} [GeV]", truthnpt, truthptbins);
0148
0149   //TH1D* h_truth_calib_frac_half1 = new TH1D("h_truth_calib_frac_half1", ";p_{T}^{Truth jet} [GeV]", truthnpt, truthptbins);
0150   //TH1D* h_measure_calib_frac_half1 = new TH1D("h_measure_calib_frac_half1", ";p_{T}^{Calib jet} [GeV]", calibnpt, calibptbins);
0151   //TH2D* h_respmatrix_calib_frac_half1 = new TH2D("h_respmatrix_calib_frac_half1", ";p_{T}^{Calib jet} [GeV];p_{T}^{Truth jet} [GeV]", calibnpt, calibptbins, truthnpt, truthptbins);
0152   //TH1D* h_fake_calib_frac_half1 = new TH1D("h_fake_calib_frac_half1", ";p_{T}^{Calib jet} [GeV]", calibnpt, calibptbins);
0153   //TH1D* h_miss_calib_frac_half1 = new TH1D("h_miss_calib_frac_half1", ";p_{T}^{Truth jet} [GeV]", truthnpt, truthptbins);
0154
0155   //TH1D* h_truth_calib_frac_half2 = new TH1D("h_truth_calib_frac_half2", ";p_{T}^{Truth jet} [GeV]", truthnpt, truthptbins);
0156   //TH1D* h_measure_calib_frac_half2 = new TH1D("h_measure_calib_frac_half2", ";p_{T}^{Calib jet} [GeV]", calibnpt, calibptbins);
0157   //TH2D* h_respmatrix_calib_frac_half2 = new TH2D("h_respmatrix_calib_frac_half2", ";p_{T}^{Calib jet} [GeV];p_{T}^{Truth jet} [GeV]", calibnpt, calibptbins, truthnpt, truthptbins);
0158   //TH1D* h_fake_calib_frac_half2 = new TH1D("h_fake_calib_frac_half2", ";p_{T}^{Calib jet} [GeV]", calibnpt, calibptbins);
0159   //TH1D* h_miss_calib_frac_half2 = new TH1D("h_miss_calib_frac_half2", ";p_{T}^{Truth jet} [GeV]", truthnpt, truthptbins);
0160
0161   ////////// Event Loop //////////
0162   std::cout << "Data analysis started." << std::endl;
0163   int n_events = chain.GetEntries();
0164   std::cout << "Total number of events: " << n_events << std::endl;
0165   // Event variables setup.
0166   std::vector<bool>* unsubjet_background_dijet = new std::vector<bool>;
0167   std::vector<bool>* unsubjet_background_frac = new std::vector<bool>;
0168   std::vector<bool> jet_filter, truthjet_filter;
0169   std::vector<float> goodtruthjet_pt, goodtruthjet_eta, goodtruthjet_phi, goodtruthjet_matched;
0170   std::vector<float> calibjet_pt_dijet, calibjet_eta_dijet, calibjet_phi_dijet, calibjet_matched_dijet;
0171   std::vector<float> calibjet_pt_dijet_jesdown, calibjet_eta_dijet_jesdown, calibjet_phi_dijet_jesdown, calibjet_matched_dijet_jesdown;
0172   std::vector<float> calibjet_pt_dijet_jesup, calibjet_eta_dijet_jesup, calibjet_phi_dijet_jesup, calibjet_matched_dijet_jesup;
0173   std::vector<float> calibjet_pt_dijet_jerdown, calibjet_eta_dijet_jerdown, calibjet_phi_dijet_jerdown, calibjet_matched_dijet_jerdown;
0174   std::vector<float> calibjet_pt_dijet_jerup, calibjet_eta_dijet_jerup, calibjet_phi_dijet_jerup, calibjet_matched_dijet_jerup;
0175   std::vector<float> calibjet_pt_frac, calibjet_eta_frac, calibjet_phi_frac, calibjet_matched_frac;
0176   std::vector<float> calibjet_pt_frac_jesdown, calibjet_eta_frac_jesdown, calibjet_phi_frac_jesdown, calibjet_matched_frac_jesdown;
0177   std::vector<float> calibjet_pt_frac_jesup, calibjet_eta_frac_jesup, calibjet_phi_frac_jesup, calibjet_matched_frac_jesup;
0178   std::vector<float> calibjet_pt_frac_jerdown, calibjet_eta_frac_jerdown, calibjet_phi_frac_jerdown, calibjet_matched_frac_jerdown;
0179   std::vector<float> calibjet_pt_frac_jerup, calibjet_eta_frac_jerup, calibjet_phi_frac_jerup, calibjet_matched_frac_jerup;
0180   for (int ie = 0; ie < n_events; ++ie) { // event loop start
0181     // Load event.
0182     if (ie % 1000 == 0) {
0183       std::cout << "Processing event " << ie << "..." << std::endl;
0184     }
0185     chain.GetEntry(ie);
0186
0187     // Z-vertex cut.
0188     if (fabs(zvertex) > 30) continue;
0189
0190     // Get truth leading jet
0191     int leadingtruthjet_index = -9999;
0192     get_leading_jet(leadingtruthjet_index, truthjet_pt);
0193     if (leadingtruthjet_index < 0) continue;
0194     float leadingtruthjet_pt = truthjet_pt->at(leadingtruthjet_index);
0195     if (leadingtruthjet_pt < truthjet_pt_min || leadingtruthjet_pt > truthjet_pt_max) continue;
0196
0197     // Get reco leading jet and subleading jet. Do basic jet cuts.
0198     int leadingunsubjet_index = -1;
0199     int subleadingunsubjet_index = -1;
0200     get_leading_subleading_jet(leadingunsubjet_index, subleadingunsubjet_index, unsubjet_pt);
0201     if (leadingunsubjet_index < 0) continue;
0202
0203     // Fill z-vertex histogram.
0204     h_zvertex->Fill(zvertex);
0205
0206     // Beam background cut. Jet level recording.
0207     unsubjet_background_dijet->clear();
0208     unsubjet_background_frac->clear();
0209     // Dijet cut.
0210     bool match_dijet = false;
0211     if (subleadingunsubjet_index >= 0 && unsubjet_e->at(subleadingunsubjet_index) / (float) unsubjet_e->at(leadingunsubjet_index) > 0.3) {
0212       match_dijet = match_leading_subleading_jet(unsubjet_phi->at(leadingunsubjet_index), unsubjet_phi->at(subleadingunsubjet_index));
0213     }
0214     bool isbeambackground_dijet = true;
0215     if (match_dijet) isbeambackground_dijet = false;
0216     for (int ij = 0; ij < unsubjet_e->size(); ++ij) {
0217       if (isbeambackground_dijet) {
0218         unsubjet_background_dijet->push_back(true);
0219       } else {
0220         unsubjet_background_dijet->push_back(false);
0221       }
0222     }
0223     // Fraction cut.
0224     for (int ij = 0; ij < unsubjet_e->size(); ++ij) {
0225       double emfrac = unsubjet_emcal_calo_e->at(ij) / (double)(unsubjet_e->at(ij));
0226       double ihfrac = unsubjet_ihcal_calo_e->at(ij) / (double)(unsubjet_e->at(ij));
0227       double ohfrac = unsubjet_ohcal_calo_e->at(ij) / (double)(unsubjet_e->at(ij));
0228       if (emfrac < 0.1 || emfrac > 0.9 || ohfrac < 0.1 || ohfrac > 0.9 || ihfrac > 0.9) {
0229         unsubjet_background_frac->push_back(true);
0230       } else {
0231         unsubjet_background_frac->push_back(false);
0232       }
0233     }
0234
0235     // Do jet filter
0236     get_jet_filter(jet_filter, unsubjet_e, unsubjet_pt, unsubjet_eta, zvertex, jet_radius);
0237     get_jet_filter(truthjet_filter, truthjet_e, truthjet_pt, truthjet_eta, zvertex, jet_radius);
0238     for (int ij = 0; ij < unsubjet_pt->size(); ++ij) {
0239       if (jet_filter[ij]) continue;
0240       h_recojet_pt_record->Fill(unsubjet_pt->at(ij), weight_scale);
0241       if (!unsubjet_background_dijet->at(ij)) h_recojet_pt_record_dijet->Fill(unsubjet_pt->at(ij), weight_scale);
0242       if (!unsubjet_background_frac->at(ij)) h_recojet_pt_record_frac->Fill(unsubjet_pt->at(ij), weight_scale);
0243     }
0244
0245     // Get truth jet and calib jet.
0246     get_truthjet(goodtruthjet_pt, goodtruthjet_eta, goodtruthjet_phi, truthjet_filter, truthjet_pt, truthjet_eta, truthjet_phi);
0247     get_calibjet(calibjet_pt_dijet, calibjet_eta_dijet, calibjet_phi_dijet, jet_filter, unsubjet_pt, unsubjet_eta, unsubjet_phi, unsubjet_background_dijet, f_corr, 1, 0.1);
0248     get_calibjet(calibjet_pt_dijet_jesdown, calibjet_eta_dijet_jesdown, calibjet_phi_dijet_jesdown, jet_filter, unsubjet_pt, unsubjet_eta, unsubjet_phi, unsubjet_background_dijet, f_corr, 0.94, 0.1);
0249     get_calibjet(calibjet_pt_dijet_jesup, calibjet_eta_dijet_jesup, calibjet_phi_dijet_jesup, jet_filter, unsubjet_pt, unsubjet_eta, unsubjet_phi, unsubjet_background_dijet, f_corr, 1.06, 0.1);
0250     get_calibjet(calibjet_pt_dijet_jerdown, calibjet_eta_dijet_jerdown, calibjet_phi_dijet_jerdown, jet_filter, unsubjet_pt, unsubjet_eta, unsubjet_phi, unsubjet_background_dijet, f_corr, 1, 0.05);
0251     get_calibjet(calibjet_pt_dijet_jerup, calibjet_eta_dijet_jerup, calibjet_phi_dijet_jerup, jet_filter, unsubjet_pt, unsubjet_eta, unsubjet_phi, unsubjet_background_dijet, f_corr, 1, 0.15);
0252     get_calibjet(calibjet_pt_frac, calibjet_eta_frac, calibjet_phi_frac, jet_filter, unsubjet_pt, unsubjet_eta, unsubjet_phi, unsubjet_background_frac, f_corr, 1, 0.1);
0253     get_calibjet(calibjet_pt_frac_jesdown, calibjet_eta_frac_jesdown, calibjet_phi_frac_jesdown, jet_filter, unsubjet_pt, unsubjet_eta, unsubjet_phi, unsubjet_background_frac, f_corr, 0.94, 0.1);
0254     get_calibjet(calibjet_pt_frac_jesup, calibjet_eta_frac_jesup, calibjet_phi_frac_jesup, jet_filter, unsubjet_pt, unsubjet_eta, unsubjet_phi, unsubjet_background_frac, f_corr, 1.06, 0.1);
0255     get_calibjet(calibjet_pt_frac_jerdown, calibjet_eta_frac_jerdown, calibjet_phi_frac_jerdown, jet_filter, unsubjet_pt, unsubjet_eta, unsubjet_phi, unsubjet_background_frac, f_corr, 1, 0.05);
0256     get_calibjet(calibjet_pt_frac_jerup, calibjet_eta_frac_jerup, calibjet_phi_frac_jerup, jet_filter, unsubjet_pt, unsubjet_eta, unsubjet_phi, unsubjet_background_frac, f_corr, 1, 0.15);
0257
0258     // Match truth jet and calib jet.
0259     match_meas_truth(calibjet_eta_dijet, calibjet_phi_dijet, calibjet_matched_dijet, goodtruthjet_eta, goodtruthjet_phi, goodtruthjet_matched, jet_radius);
0260     fill_response_matrix(h_truth_calib_dijet, h_measure_calib_dijet, h_respmatrix_calib_dijet, h_fake_calib_dijet, h_miss_calib_dijet, calibjet_pt_dijet, calibjet_matched_dijet, goodtruthjet_pt, goodtruthjet_matched, weight_scale);
0261     match_meas_truth(calibjet_eta_dijet_jesdown, calibjet_phi_dijet_jesdown, calibjet_matched_dijet_jesdown, goodtruthjet_eta, goodtruthjet_phi, goodtruthjet_matched, jet_radius);
0262     fill_response_matrix(h_truth_calib_dijet_jesdown, h_measure_calib_dijet_jesdown, h_respmatrix_calib_dijet_jesdown, h_fake_calib_dijet_jesdown, h_miss_calib_dijet_jesdown, calibjet_pt_dijet_jesdown, calibjet_matched_dijet_jesdown, goodtruthjet_pt, goodtruthjet_matched, weight_scale);
0263     match_meas_truth(calibjet_eta_dijet_jesup, calibjet_phi_dijet_jesup, calibjet_matched_dijet_jesup, goodtruthjet_eta, goodtruthjet_phi, goodtruthjet_matched, jet_radius);
0264     fill_response_matrix(h_truth_calib_dijet_jesup, h_measure_calib_dijet_jesup, h_respmatrix_calib_dijet_jesup, h_fake_calib_dijet_jesup, h_miss_calib_dijet_jesup, calibjet_pt_dijet_jesup, calibjet_matched_dijet_jesup, goodtruthjet_pt, goodtruthjet_matched, weight_scale);
0265     match_meas_truth(calibjet_eta_dijet_jerdown, calibjet_phi_dijet_jerdown, calibjet_matched_dijet_jerdown, goodtruthjet_eta, goodtruthjet_phi, goodtruthjet_matched, jet_radius);
0266     fill_response_matrix(h_truth_calib_dijet_jerdown, h_measure_calib_dijet_jerdown, h_respmatrix_calib_dijet_jerdown, h_fake_calib_dijet_jerdown, h_miss_calib_dijet_jerdown, calibjet_pt_dijet_jerdown, calibjet_matched_dijet_jerdown, goodtruthjet_pt, goodtruthjet_matched, weight_scale);
0267     match_meas_truth(calibjet_eta_dijet_jerup, calibjet_phi_dijet_jerup, calibjet_matched_dijet_jerup, goodtruthjet_eta, goodtruthjet_phi, goodtruthjet_matched, jet_radius);
0268     fill_response_matrix(h_truth_calib_dijet_jerup, h_measure_calib_dijet_jerup, h_respmatrix_calib_dijet_jerup, h_fake_calib_dijet_jerup, h_miss_calib_dijet_jerup, calibjet_pt_dijet_jerup, calibjet_matched_dijet_jerup, goodtruthjet_pt, goodtruthjet_matched, weight_scale);
0269     match_meas_truth(calibjet_eta_frac, calibjet_phi_frac, calibjet_matched_frac, goodtruthjet_eta, goodtruthjet_phi, goodtruthjet_matched, jet_radius);
0270     fill_response_matrix(h_truth_calib_frac, h_measure_calib_frac, h_respmatrix_calib_frac, h_fake_calib_frac, h_miss_calib_frac, calibjet_pt_frac, calibjet_matched_frac, goodtruthjet_pt, goodtruthjet_matched, weight_scale);
0271     //if(ie%2 == 0){
0272     //  fill_response_matrix(h_truth_calib_frac_half1, h_measure_calib_frac_half1, h_respmatrix_calib_frac_half1, h_fake_calib_frac_half1, h_miss_calib_frac_half1, calibjet_pt_frac, calibjet_matched_frac, goodtruthjet_pt, goodtruthjet_matched, weight_scale);
0273     //} else {
0274     //  fill_response_matrix(h_truth_calib_frac_half2, h_measure_calib_frac_half2, h_respmatrix_calib_frac_half2, h_fake_calib_frac_half2, h_miss_calib_frac_half2, calibjet_pt_frac, calibjet_matched_frac, goodtruthjet_pt, goodtruthjet_matched, weight_scale);
0275     //}
0276     match_meas_truth(calibjet_eta_frac_jesdown, calibjet_phi_frac_jesdown, calibjet_matched_frac_jesdown, goodtruthjet_eta, goodtruthjet_phi, goodtruthjet_matched, jet_radius);
0277     fill_response_matrix(h_truth_calib_frac_jesdown, h_measure_calib_frac_jesdown, h_respmatrix_calib_frac_jesdown, h_fake_calib_frac_jesdown, h_miss_calib_frac_jesdown, calibjet_pt_frac_jesdown, calibjet_matched_frac_jesdown, goodtruthjet_pt, goodtruthjet_matched, weight_scale);
0278     match_meas_truth(calibjet_eta_frac_jesup, calibjet_phi_frac_jesup, calibjet_matched_frac_jesup, goodtruthjet_eta, goodtruthjet_phi, goodtruthjet_matched, jet_radius);
0279     fill_response_matrix(h_truth_calib_frac_jesup, h_measure_calib_frac_jesup, h_respmatrix_calib_frac_jesup, h_fake_calib_frac_jesup, h_miss_calib_frac_jesup, calibjet_pt_frac_jesup, calibjet_matched_frac_jesup, goodtruthjet_pt, goodtruthjet_matched, weight_scale);
0280     match_meas_truth(calibjet_eta_frac_jerdown, calibjet_phi_frac_jerdown, calibjet_matched_frac_jerdown, goodtruthjet_eta, goodtruthjet_phi, goodtruthjet_matched, jet_radius);
0281     fill_response_matrix(h_truth_calib_frac_jerdown, h_measure_calib_frac_jerdown, h_respmatrix_calib_frac_jerdown, h_fake_calib_frac_jerdown, h_miss_calib_frac_jerdown, calibjet_pt_frac_jerdown, calibjet_matched_frac_jerdown, goodtruthjet_pt, goodtruthjet_matched, weight_scale);
0282     match_meas_truth(calibjet_eta_frac_jerup, calibjet_phi_frac_jerup, calibjet_matched_frac_jerup, goodtruthjet_eta, goodtruthjet_phi, goodtruthjet_matched, jet_radius);
0283     fill_response_matrix(h_truth_calib_frac_jerup, h_measure_calib_frac_jerup, h_respmatrix_calib_frac_jerup, h_fake_calib_frac_jerup, h_miss_calib_frac_jerup, calibjet_pt_frac_jerup, calibjet_matched_frac_jerup, goodtruthjet_pt, goodtruthjet_matched, weight_scale);
0284   } // event loop end
0285
0286   // Write histograms.
0287   std::cout << "Writing histograms..." << std::endl;
0288   f_out->cd();
0289
0290   h_zvertex->Write();
0291   h_recojet_pt_record->Write();
0292   h_recojet_pt_record_dijet->Write();
0293   h_recojet_pt_record_frac->Write();
0294
0295   h_truth_calib_dijet->Write();
0296   h_measure_calib_dijet->Write();
0297   h_respmatrix_calib_dijet->Write();
0298   h_fake_calib_dijet->Write();
0299   h_miss_calib_dijet->Write();
0300
0301   h_truth_calib_dijet_jesdown->Write();
0302   h_measure_calib_dijet_jesdown->Write();
0303   h_respmatrix_calib_dijet_jesdown->Write();
0304   h_fake_calib_dijet_jesdown->Write();
0305   h_miss_calib_dijet_jesdown->Write();
0306
0307   h_truth_calib_dijet_jesup->Write();
0308   h_measure_calib_dijet_jesup->Write();
0309   h_respmatrix_calib_dijet_jesup->Write();
0310   h_fake_calib_dijet_jesup->Write();
0311   h_miss_calib_dijet_jesup->Write();
0312
0313   h_truth_calib_dijet_jerdown->Write();
0314   h_measure_calib_dijet_jerdown->Write();
0315   h_respmatrix_calib_dijet_jerdown->Write();
0316   h_fake_calib_dijet_jerdown->Write();
0317   h_miss_calib_dijet_jerdown->Write();
0318
0319   h_truth_calib_dijet_jerup->Write();
0320   h_measure_calib_dijet_jerup->Write();
0321   h_respmatrix_calib_dijet_jerup->Write();
0322   h_fake_calib_dijet_jerup->Write();
0323   h_miss_calib_dijet_jerup->Write();
0324
0325   h_truth_calib_frac->Write();
0326   h_measure_calib_frac->Write();
0327   h_respmatrix_calib_frac->Write();
0328   h_fake_calib_frac->Write();
0329   h_miss_calib_frac->Write();
0330
0331   h_truth_calib_frac_jesdown->Write();
0332   h_measure_calib_frac_jesdown->Write();
0333   h_respmatrix_calib_frac_jesdown->Write();
0334   h_fake_calib_frac_jesdown->Write();
0335   h_miss_calib_frac_jesdown->Write();
0336
0337   h_truth_calib_frac_jesup->Write();
0338   h_measure_calib_frac_jesup->Write();
0339   h_respmatrix_calib_frac_jesup->Write();
0340   h_fake_calib_frac_jesup->Write();
0341   h_miss_calib_frac_jesup->Write();
0342
0343   h_truth_calib_frac_jerdown->Write();
0344   h_measure_calib_frac_jerdown->Write();
0345   h_respmatrix_calib_frac_jerdown->Write();
0346   h_fake_calib_frac_jerdown->Write();
0347   h_miss_calib_frac_jerdown->Write();
0348
0349   h_truth_calib_frac_jerup->Write();
0350   h_measure_calib_frac_jerup->Write();
0351   h_respmatrix_calib_frac_jerup->Write();
0352   h_fake_calib_frac_jerup->Write();
0353   h_miss_calib_frac_jerup->Write();
0354
0355   //h_truth_calib_frac_half1->Write();
0356   //h_measure_calib_frac_half1->Write();
0357   //h_respmatrix_calib_frac_half1->Write();
0358   //h_fake_calib_frac_half1->Write();
0359   //h_miss_calib_frac_half1->Write();
0360
0361   //h_truth_calib_frac_half2->Write();
0362   //h_measure_calib_frac_half2->Write();
0363   //h_respmatrix_calib_frac_half2->Write();
0364   //h_fake_calib_frac_half2->Write();
0365   //h_miss_calib_frac_half2->Write();
0366
0367   f_out->Close();
0368   std::cout << "All done!" << std::endl;
0369 }
0370
0371 ////////////////////////////////////////// Helper Functions //////////////////////////////////////////
0372 float get_deta(float eta1, float eta2) {
0373   return eta1 - eta2;
0374 }
0375
0376 float get_dphi(float phi1, float phi2) {
0377   float dphi1 = phi1 - phi2;
0378   float dphi2 = phi1 - phi2 + 2*TMath::Pi();
0379   float dphi3 = phi1 - phi2 - 2*TMath::Pi();
0380   if (fabs(dphi1) > fabs(dphi2)) {
0381     dphi1 = dphi2;
0382   }
0383   if (fabs(dphi1) > fabs(dphi3)) {
0384     dphi1 = dphi3;
0385   }
0386   return dphi1;
0387 }
0388
0389 float get_dR(float eta1, float phi1, float eta2, float phi2) {
0390   float deta = get_deta(eta1, eta2);
0391   float dphi = get_dphi(phi1, phi2);
0392   return sqrt(deta*deta + dphi*dphi);
0393 }
0394
0395 float get_emcal_mineta_zcorrected(float zvertex) {
0396   float minz_EM = -130.23;
0397   float radius_EM = 93.5;
0398   float z = minz_EM - zvertex;
0399   float eta_zcorrected = asinh(z / (float)radius_EM);
0400   return eta_zcorrected;
0401 }
0402
0403 float get_emcal_maxeta_zcorrected(float zvertex) {
0404   float maxz_EM = 130.23;
0405   float radius_EM = 93.5;
0406   float z = maxz_EM - zvertex;
0407   float eta_zcorrected = asinh(z / (float)radius_EM);
0408   return eta_zcorrected;
0409 }
0410
0411 float get_ihcal_mineta_zcorrected(float zvertex) {
0412   float minz_IH = -170.299;
0413   float radius_IH = 127.503;
0414   float z = minz_IH - zvertex;
0415   float eta_zcorrected = asinh(z / (float)radius_IH);
0416   return eta_zcorrected;
0417 }
0418
0419 float get_ihcal_maxeta_zcorrected(float zvertex) {
0420   float maxz_IH = 170.299;
0421   float radius_IH = 127.503;
0422   float z = maxz_IH - zvertex;
0423   float eta_zcorrected = asinh(z / (float)radius_IH);
0424   return eta_zcorrected;
0425 }
0426
0427 float get_ohcal_mineta_zcorrected(float zvertex) {
0428   float minz_OH = -301.683;
0429   float radius_OH = 225.87;
0430   float z = minz_OH - zvertex;
0431   float eta_zcorrected = asinh(z / (float)radius_OH);
0432   return eta_zcorrected;
0433 }
0434
0435 float get_ohcal_maxeta_zcorrected(float zvertex) {
0436   float maxz_OH = 301.683;
0437   float radius_OH = 225.87;
0438   float z = maxz_OH - zvertex;
0439   float eta_zcorrected = asinh(z / (float)radius_OH);
0440   return eta_zcorrected;
0441 }
0442
0443 bool check_bad_jet_eta(float jet_eta, float zertex, float jet_radius) {
0444   float emcal_mineta = get_emcal_mineta_zcorrected(zertex);
0445   float emcal_maxeta = get_emcal_maxeta_zcorrected(zertex);
0446   float ihcal_mineta = get_ihcal_mineta_zcorrected(zertex);
0447   float ihcal_maxeta = get_ihcal_maxeta_zcorrected(zertex);
0448   float ohcal_mineta = get_ohcal_mineta_zcorrected(zertex);
0449   float ohcal_maxeta = get_ohcal_maxeta_zcorrected(zertex);
0450   float minlimit = emcal_mineta;
0451   if (ihcal_mineta > minlimit) minlimit = ihcal_mineta;
0452   if (ohcal_mineta > minlimit) minlimit = ohcal_mineta;
0453   float maxlimit = emcal_maxeta;
0454   if (ihcal_maxeta < maxlimit) maxlimit = ihcal_maxeta;
0455   if (ohcal_maxeta < maxlimit) maxlimit = ohcal_maxeta;
0456   minlimit += jet_radius;
0457   maxlimit -= jet_radius;
0458   return jet_eta < minlimit || jet_eta > maxlimit;
0459 }
0460
0461 ////////////////////////////////////////// Functions //////////////////////////////////////////
0462 void get_leading_jet(int& leadingjet_index, std::vector<float>* jet_pt) {
0463   leadingjet_index = -1;
0464   float leadingjet_pt = -9999;
0465   for (int ij = 0; ij < jet_pt->size(); ++ij) {
0466     float jetpt = jet_pt->at(ij);
0467     if (jetpt > leadingjet_pt) {
0468       leadingjet_pt = jetpt;
0469       leadingjet_index = ij;
0470     }
0471   }
0472 }
0473
0474 void get_leading_subleading_jet(int& leadingjet_index, int& subleadingjet_index, std::vector<float>* jet_et) {
0475   leadingjet_index = -1;
0476   subleadingjet_index = -1;
0477   float leadingjet_et = -9999;
0478   float subleadingjet_et = -9999;
0479   for (int ij = 0; ij < jet_et->size(); ++ij) {
0480     float jetet = jet_et->at(ij);
0481     if (jetet > leadingjet_et) {
0482       subleadingjet_et = leadingjet_et;
0483       subleadingjet_index = leadingjet_index;
0484       leadingjet_et = jetet;
0485       leadingjet_index = ij;
0486     } else if (jetet > subleadingjet_et) {
0487       subleadingjet_et = jetet;
0488       subleadingjet_index = ij;
0489     }
0490   }
0491 }
0492
0493 bool match_leading_subleading_jet(float leadingjet_phi, float subleadingjet_phi) {
0494   float dijet_min_phi = 3*TMath::Pi()/4.;
0495   float dphi = get_dphi(leadingjet_phi, subleadingjet_phi);
0496   return dphi > dijet_min_phi;
0497 }
0498
0499 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) {
0500   jet_filter.clear();
0501   int njet = jet_e->size();
0502   for (int ij = 0; ij < njet; ++ij) {
0503     jet_filter.push_back(jet_e->at(ij) < 0 || check_bad_jet_eta(jet_eta->at(ij), zvertex, jet_radius));
0504   }
0505 }
0506
0507 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, std::vector<bool>* jet_background, TF1* f_corr, float jes_para, float jer_para) {
0508   calibjet_pt.clear();
0509   calibjet_eta.clear();
0510   calibjet_phi.clear();
0511   for (int ij = 0; ij < jet_filter.size(); ++ij) {
0512     if (jet_filter.at(ij) || jet_background->at(ij)) continue;
0513     double calib_pt = f_corr->Eval(jet_pt->at(ij)) * (1 + randGen.Gaus(0.0, jer_para)) * jes_para;
0514     if (calib_pt < calibptbins[0] || calib_pt > calibptbins[calibnpt]) continue;
0515     calibjet_pt.push_back(calib_pt);
0516     calibjet_eta.push_back(jet_eta->at(ij));
0517     calibjet_phi.push_back(jet_phi->at(ij));
0518   }
0519 }
0520
0521 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) {
0522   goodtruthjet_pt.clear();
0523   goodtruthjet_eta.clear();
0524   goodtruthjet_phi.clear();
0525   for (int ij = 0; ij < truthjet_filter.size(); ++ij) {
0526     if (truthjet_filter.at(ij)) continue;
0527     if (jet_pt->at(ij) < truthptbins[0] || jet_pt->at(ij) > truthptbins[truthnpt]) continue;
0528     goodtruthjet_pt.push_back(jet_pt->at(ij));
0529     goodtruthjet_eta.push_back(jet_eta->at(ij));
0530     goodtruthjet_phi.push_back(jet_phi->at(ij));
0531   }
0532 }
0533
0534 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) {
0535   meas_matched.assign(meas_eta.size(), -1);
0536   truth_matched.assign(truth_eta.size(), -1);
0537   float max_match_dR = jet_radius * 0.75;
0538   for (int im = 0; im < meas_eta.size(); ++im) {
0539     float min_dR = 100;
0540     int match_index = -9999;
0541     for (int it = 0; it < truth_eta.size(); ++it) {
0542       float dR = get_dR(meas_eta[im], meas_phi[im], truth_eta[it], truth_phi[it]);
0543       if (dR < min_dR) {
0544         match_index = it;
0545         min_dR = dR;
0546       }
0547     }
0548     if (min_dR < max_match_dR) {
0549       if (truth_matched[match_index] == -1) {
0550         meas_matched[im] = match_index;
0551         truth_matched[match_index] = im;
0552       } else {
0553         float dR1 = get_dR(meas_eta[truth_matched[match_index]], meas_phi[truth_matched[match_index]], truth_eta[match_index], truth_phi[match_index]);
0554         if (min_dR < dR1) {
0555           meas_matched[truth_matched[match_index]] = -1;
0556           meas_matched[im] = match_index;
0557           truth_matched[match_index] = im;
0558         }
0559       }
0560     }
0561   }
0562 }
0563
0564 void fill_response_matrix(TH1D*& h_truth, TH1D*& h_meas, TH2D*& h_resp, TH1D*& h_fake, TH1D*& h_miss, std::vector<float>& meas_pt, std::vector<float>& meas_matched, std::vector<float>& truth_pt, std::vector<float>& truth_matched, float weight_scale) {
0565   for (int im = 0; im < meas_pt.size(); ++im) {
0566     h_meas->Fill(meas_pt[im], weight_scale);
0567     if (meas_matched[im] < 0) {
0568       h_fake->Fill(meas_pt[im], weight_scale);
0569     } else {
0570       h_resp->Fill(meas_pt[im], truth_pt[meas_matched[im]], weight_scale);
0571     }
0572   }
0573   for (int it = 0; it < truth_pt.size(); ++it) {
0574     h_truth->Fill(truth_pt[it], weight_scale);
0575     if (truth_matched[it] < 0) {
0576       h_miss->Fill(truth_pt[it], weight_scale);
0577     }
0578   }
0579 }