File indexing completed on 2023-03-17 11:13:06
0001
0002
0003
0004 #include "TROOT.h"
0005 #include "TStyle.h"
0006 #include "TLatex.h"
0007 #include "TFile.h"
0008 #include "TTree.h"
0009 #include "TChain.h"
0010 #include "TBranch.h"
0011 #include "TLeaf.h"
0012 #include "TCanvas.h"
0013 #include "TLegend.h"
0014 #include "TH1.h"
0015 #include "TH2.h"
0016 #include "TF1.h"
0017 #include "TProfile.h"
0018 #include "TProfile2D.h"
0019 #include "TMath.h"
0020 #include <TError.h>
0021
0022 #include <iostream>
0023 #include <string>
0024 #include <vector>
0025
0026 using namespace std;
0027
0028 void trkMET_turnOn() {
0029 TChain* tree = new TChain("L1TrackNtuple/eventTree");
0030
0031 tree->Add("CheckingJets_CMSSW11_CMS.root");
0032
0033 if (tree->GetEntries() == 0) {
0034 cout << "File doesn't exist or is empty, returning..." << endl;
0035 return;
0036 }
0037
0038
0039 float trueMET = 0;
0040 float trueTkMET = 0;
0041 float trkMET = 0;
0042 vector<float>* pv_L1reco;
0043
0044 vector<float>* gen_pt;
0045 vector<float>* gen_phi;
0046 vector<int>* gen_pdgid;
0047
0048
0049 vector<float>* tp_pt;
0050 vector<float>* tp_eta;
0051 vector<float>* tp_phi;
0052 vector<float>* tp_dxy;
0053 vector<float>* tp_z0;
0054 vector<float>* tp_d0;
0055 vector<int>* tp_pdgid;
0056 vector<int>* tp_nmatch;
0057 vector<int>* tp_nstub;
0058
0059 vector<int>* tp_eventid;
0060
0061 vector<int>* tp_charge;
0062
0063
0064 vector<float>* matchtrk_pt;
0065 vector<float>* matchtrk_eta;
0066 vector<float>* matchtrk_phi;
0067 vector<float>* matchtrk_d0;
0068 vector<float>* matchtrk_z0;
0069 vector<float>* matchtrk_chi2;
0070 vector<float>* matchtrk_MVA1;
0071 vector<float>* matchtrk_bendchi2;
0072 vector<int>* matchtrk_nstub;
0073 vector<int>* matchtrk_seed;
0074
0075
0076 vector<float>* trk_pt;
0077 vector<float>* trk_eta;
0078 vector<float>* trk_phi;
0079 vector<float>* trk_z0;
0080 vector<float>* trk_chi2;
0081 vector<float>* trk_MVA1;
0082 vector<float>* trk_bendchi2;
0083 vector<int>* trk_nstub;
0084 vector<int>* trk_seed;
0085 vector<int>* trk_fake;
0086
0087 TBranch* b_gen_pt;
0088 TBranch* b_gen_phi;
0089 TBranch* b_gen_pdgid;
0090
0091 TBranch* b_tp_pt;
0092 TBranch* b_tp_eta;
0093 TBranch* b_tp_phi;
0094 TBranch* b_tp_dxy;
0095 TBranch* b_tp_z0;
0096 TBranch* b_tp_d0;
0097 TBranch* b_tp_pdgid;
0098 TBranch* b_tp_nmatch;
0099
0100 TBranch* b_tp_nstub;
0101 TBranch* b_tp_eventid;
0102
0103 TBranch* b_tp_charge;
0104
0105 TBranch* b_matchtrk_pt;
0106 TBranch* b_matchtrk_eta;
0107 TBranch* b_matchtrk_phi;
0108 TBranch* b_matchtrk_d0;
0109 TBranch* b_matchtrk_z0;
0110 TBranch* b_matchtrk_chi2;
0111 TBranch* b_matchtrk_MVA1;
0112 TBranch* b_matchtrk_bendchi2;
0113 TBranch* b_matchtrk_nstub;
0114 TBranch* b_matchtrk_seed;
0115
0116 TBranch* b_trk_pt;
0117 TBranch* b_trk_eta;
0118 TBranch* b_trk_phi;
0119 TBranch* b_trk_z0;
0120 TBranch* b_trk_chi2;
0121 TBranch* b_trk_MVA1;
0122 TBranch* b_trk_bendchi2;
0123 TBranch* b_trk_nstub;
0124 TBranch* b_trk_seed;
0125 TBranch* b_trk_fake;
0126
0127 TBranch* b_pv_L1reco;
0128 TBranch* b_trueMET;
0129 TBranch* b_trueTkMET;
0130 TBranch* b_trkMET;
0131
0132 trueMET = 0;
0133 trueTkMET = 0;
0134 trkMET = 0;
0135
0136 gen_pt = 0;
0137 gen_phi = 0;
0138 gen_pdgid = 0;
0139
0140 tp_pt = 0;
0141 tp_eta = 0;
0142 tp_phi = 0;
0143 tp_dxy = 0;
0144 tp_z0 = 0;
0145 tp_d0 = 0;
0146 tp_pdgid = 0;
0147 tp_nmatch = 0;
0148
0149 tp_nstub = 0;
0150 tp_eventid = 0;
0151
0152 tp_charge = 0;
0153
0154 matchtrk_pt = 0;
0155 matchtrk_eta = 0;
0156 matchtrk_phi = 0;
0157 matchtrk_d0 = 0;
0158 matchtrk_z0 = 0;
0159 matchtrk_chi2 = 0;
0160 matchtrk_MVA1 = 0;
0161 matchtrk_bendchi2 = 0;
0162 matchtrk_nstub = 0;
0163 matchtrk_seed = 0;
0164
0165 trk_pt = 0;
0166 trk_eta = 0;
0167 trk_phi = 0;
0168 trk_z0 = 0;
0169 trk_chi2 = 0;
0170 trk_MVA1 = 0;
0171 trk_bendchi2 = 0;
0172 trk_nstub = 0;
0173 trk_seed = 0;
0174 trk_fake = 0;
0175
0176 pv_L1reco = 0;
0177
0178 tree->SetBranchAddress("pv_L1reco", &pv_L1reco, &b_pv_L1reco);
0179 tree->SetBranchAddress("trueMET", &trueMET, &b_trueMET);
0180 tree->SetBranchAddress("trueTkMET", &trueTkMET, &b_trueTkMET);
0181 tree->SetBranchAddress("trkMET", &trkMET, &b_trkMET);
0182 tree->SetBranchAddress("gen_pt", &gen_pt, &b_gen_pt);
0183 tree->SetBranchAddress("gen_phi", &gen_phi, &b_gen_phi);
0184 tree->SetBranchAddress("gen_pdgid", &gen_pdgid, &b_gen_pdgid);
0185
0186 tree->SetBranchAddress("tp_pt", &tp_pt, &b_tp_pt);
0187 tree->SetBranchAddress("tp_eta", &tp_eta, &b_tp_eta);
0188 tree->SetBranchAddress("tp_phi", &tp_phi, &b_tp_phi);
0189 tree->SetBranchAddress("tp_dxy", &tp_dxy, &b_tp_dxy);
0190 tree->SetBranchAddress("tp_z0", &tp_z0, &b_tp_z0);
0191 tree->SetBranchAddress("tp_d0", &tp_d0, &b_tp_d0);
0192 tree->SetBranchAddress("tp_pdgid", &tp_pdgid, &b_tp_pdgid);
0193 tree->SetBranchAddress("tp_nmatch", &tp_nmatch, &b_tp_nmatch);
0194
0195 tree->SetBranchAddress("tp_nstub", &tp_nstub, &b_tp_nstub);
0196 tree->SetBranchAddress("tp_eventid", &tp_eventid, &b_tp_eventid);
0197
0198 tree->SetBranchAddress("tp_charge", &tp_charge, &b_tp_charge);
0199
0200 tree->SetBranchAddress("matchtrk_pt", &matchtrk_pt, &b_matchtrk_pt);
0201 tree->SetBranchAddress("matchtrk_eta", &matchtrk_eta, &b_matchtrk_eta);
0202 tree->SetBranchAddress("matchtrk_phi", &matchtrk_phi, &b_matchtrk_phi);
0203 tree->SetBranchAddress("matchtrk_d0", &matchtrk_d0, &b_matchtrk_d0);
0204 tree->SetBranchAddress("matchtrk_z0", &matchtrk_z0, &b_matchtrk_z0);
0205 tree->SetBranchAddress("matchtrk_chi2", &matchtrk_chi2, &b_matchtrk_chi2);
0206 tree->SetBranchAddress("matchtrk_MVA1", &matchtrk_MVA1, &b_matchtrk_MVA1);
0207 tree->SetBranchAddress("matchtrk_bendchi2", &matchtrk_bendchi2, &b_matchtrk_bendchi2);
0208 tree->SetBranchAddress("matchtrk_nstub", &matchtrk_nstub, &b_matchtrk_nstub);
0209
0210 tree->SetBranchAddress("trk_pt", &trk_pt, &b_trk_pt);
0211 tree->SetBranchAddress("trk_eta", &trk_eta, &b_trk_eta);
0212 tree->SetBranchAddress("trk_phi", &trk_phi, &b_trk_phi);
0213 tree->SetBranchAddress("trk_z0", &trk_z0, &b_trk_z0);
0214 tree->SetBranchAddress("trk_chi2", &trk_chi2, &b_trk_chi2);
0215 tree->SetBranchAddress("trk_MVA1", &trk_MVA1, &b_trk_MVA1);
0216 tree->SetBranchAddress("trk_bendchi2", &trk_bendchi2, &b_trk_bendchi2);
0217 tree->SetBranchAddress("trk_nstub", &trk_nstub, &b_trk_nstub);
0218 tree->SetBranchAddress("trk_fake", &trk_fake, &b_trk_fake);
0219
0220
0221
0222 float numBins = 100.0;
0223 float cutoff = 300.0 - 0.1;
0224 float TP_minPt = 2.0;
0225 float TP_maxEta = 2.4;
0226 float chi2dof_cut = 10.0;
0227 float bendchi2_cut = 2.2;
0228
0229 TH1F* h_trueMET = new TH1F("trueMET", ";trueMET [GeV]; Events", numBins, 0, 500.0);
0230 TH1F* h_trueTkMET = new TH1F("trueTkMET", ";trueTkMET [GeV]; Events", numBins, 0, 500.0);
0231 TH1F* h_recoTkMET = new TH1F("recoTkMET", ";recoTkMET [GeV]; Events", numBins, 0, 500.0);
0232
0233 TH1F* h_trueTkMET_num = new TH1F("trueTkMET_turnon", ";trueMET [GeV]; Events", numBins, 0, 500.0);
0234 TH1F* h_recoTkMET_num = new TH1F("recoTkMET_turnon", ";trueMET [GeV]; Events", numBins, 0, 500.0);
0235
0236
0237 float trueTkMET_rate = 49.0;
0238 float recoTkMET_rate = 49.0;
0239
0240 bool rerun_trueMET = true;
0241 bool rerun_trueTkMET = true;
0242 bool rerun_recoTkMET = true;
0243
0244 int nevt = tree->GetEntries();
0245 cout << "number of events = " << nevt << endl;
0246
0247
0248 for (int i = 0; i < nevt; i++) {
0249 tree->GetEntry(i);
0250 if (i % 10000 == 0)
0251 std::cout << i << "/" << nevt << std::endl;
0252
0253 float trueMET_ntuple = trueMET;
0254 float trueTkMET_ntuple = trueTkMET;
0255 float recoTkMET_ntuple = trkMET;
0256
0257 if (rerun_trueMET) {
0258 float trueMETx = 0.0;
0259 float trueMETy = 0.0;
0260 float trueMET_calc = 0.0;
0261 for (size_t i = 0; i < gen_pt->size(); ++i) {
0262 int id = gen_pdgid->at(i);
0263 float pt = gen_pt->at(i);
0264 float phi = gen_phi->at(i);
0265 bool isNeutrino = false;
0266 if ((fabs(id) == 12 || fabs(id) == 14 || fabs(id) == 16))
0267 isNeutrino = true;
0268 if ((isNeutrino || id == 1000022)) {
0269 trueMETx += pt * cos(phi);
0270 trueMETy += pt * sin(phi);
0271 }
0272 }
0273 trueMET_calc = sqrt(trueMETx * trueMETx + trueMETy * trueMETy);
0274 trueMET_ntuple = trueMET_calc;
0275 }
0276
0277 if (rerun_trueTkMET) {
0278 float trueTkMET_calc = 0;
0279 float trueTkMETx = 0;
0280 float trueTkMETy = 0;
0281
0282 for (int it = 0; it < (int)tp_pt->size(); it++) {
0283 float this_tp_pt = tp_pt->at(it);
0284 float this_tp_eta = tp_eta->at(it);
0285 float this_tp_phi = tp_phi->at(it);
0286 int this_tp_eventID = tp_eventid->at(it);
0287
0288
0289 if (tp_pt->at(it) < TP_minPt)
0290 continue;
0291 if (fabs(this_tp_eta) > TP_maxEta)
0292 continue;
0293 if (tp_nstub->at(it) < 4)
0294 continue;
0295 if (fabs(tp_z0->at(it)) > 15)
0296 continue;
0297 if (this_tp_eventID != 0)
0298 continue;
0299
0300 if (tp_charge->at(it) == 0)
0301 continue;
0302
0303 trueTkMETx += this_tp_pt * cos(this_tp_phi);
0304 trueTkMETy += this_tp_pt * sin(this_tp_phi);
0305 }
0306 trueTkMET_calc = sqrt(trueTkMETx * trueTkMETx + trueTkMETy * trueTkMETy);
0307 trueTkMET_ntuple = trueTkMET_calc;
0308 }
0309
0310 if (rerun_recoTkMET) {
0311 float recoTkMET_calc = 0;
0312 float recoTkMETx = 0;
0313 float recoTkMETy = 0;
0314
0315 for (int it = 0; it < (int)trk_pt->size(); it++) {
0316 float thisTrk_pt = trk_pt->at(it);
0317 float thisTrk_eta = trk_eta->at(it);
0318 int thisTrk_nstub = trk_nstub->at(it);
0319 float thisTrk_chi2dof = trk_chi2->at(it) / (2 * thisTrk_nstub - 4);
0320 float thisTrk_bendchi2 = trk_bendchi2->at(it);
0321 float thisTrk_phi = trk_phi->at(it);
0322 float thisTrk_MVA = trk_MVA1->at(it);
0323 int thisTrk_fake = trk_fake->at(it);
0324 float thisTrk_z0 = trk_z0->at(it);
0325 float deltaZ = thisTrk_z0 - pv_L1reco->at(0);
0326
0327 float deltaZ_cut = 3.0;
0328 if (fabs(thisTrk_eta) >= 0 && fabs(thisTrk_eta) < 0.7)
0329 deltaZ_cut = 0.4;
0330 else if (fabs(thisTrk_eta) >= 0.7 && fabs(thisTrk_eta) < 1.0)
0331 deltaZ_cut = 0.6;
0332 else if (fabs(thisTrk_eta) >= 1.0 && fabs(thisTrk_eta) < 1.2)
0333 deltaZ_cut = 0.76;
0334 else if (fabs(thisTrk_eta) >= 1.2 && fabs(thisTrk_eta) < 1.6)
0335 deltaZ_cut = 1.0;
0336 else if (fabs(thisTrk_eta) >= 1.6 && fabs(thisTrk_eta) < 2.0)
0337 deltaZ_cut = 1.7;
0338 else if (fabs(thisTrk_eta) >= 2.0 && fabs(thisTrk_eta) <= 2.4)
0339 deltaZ_cut = 2.20;
0340
0341 if (thisTrk_pt < TP_minPt || fabs(thisTrk_eta) > TP_maxEta || thisTrk_nstub < 4 || fabs(thisTrk_z0) > 15)
0342 continue;
0343 if (fabs(deltaZ) > deltaZ_cut)
0344 continue;
0345
0346
0347
0348 if (thisTrk_chi2dof < chi2dof_cut && thisTrk_bendchi2 < bendchi2_cut) {
0349 recoTkMETx += thisTrk_pt * cos(thisTrk_phi);
0350 recoTkMETy += thisTrk_pt * sin(thisTrk_phi);
0351 }
0352 }
0353 recoTkMET_calc = sqrt(recoTkMETx * recoTkMETx + recoTkMETy * recoTkMETy);
0354 if (recoTkMET_calc != recoTkMET_ntuple)
0355 std::cout << "Ugh: " << recoTkMET_calc << ", " << recoTkMET_ntuple << std::endl;
0356
0357 recoTkMET_ntuple = recoTkMET_calc;
0358 }
0359
0360 h_trueMET->Fill(trueMET_ntuple);
0361 h_trueTkMET->Fill(trueTkMET_ntuple);
0362 h_recoTkMET->Fill(recoTkMET_ntuple);
0363
0364
0365
0366 if (trueTkMET_ntuple > 49)
0367 h_trueTkMET_num->Fill(trueMET_ntuple);
0368 if (recoTkMET_ntuple > 49)
0369 h_recoTkMET_num->Fill(trueMET_ntuple);
0370 }
0371
0372
0373 h_trueMET->SetLineColor(kBlack);
0374 h_trueMET->SetLineStyle(2);
0375 h_trueTkMET->SetLineColor(kBlack);
0376 h_recoTkMET->SetLineColor(kRed);
0377
0378 TLegend* leg = new TLegend(0.55, 0.15, 0.85, 0.35);
0379 leg->SetBorderSize(0);
0380 leg->SetFillStyle(0);
0381 leg->SetTextSize(0.04);
0382 leg->SetTextFont(42);
0383 leg->AddEntry(h_trueTkMET, "trueTkMET", "l");
0384 leg->AddEntry(h_recoTkMET, "recoTkMET", "l");
0385
0386
0387 h_trueMET->Rebin(2);
0388 h_trueTkMET_num->Rebin(2);
0389 h_recoTkMET_num->Rebin(2);
0390 h_trueMET->Sumw2();
0391 h_trueTkMET->Sumw2();
0392 h_trueTkMET_num->Sumw2();
0393
0394 TH1F* h_TkMET_turnon = (TH1F*)h_trueTkMET_num->Clone();
0395 h_TkMET_turnon->SetName("trueTkMET_turnon");
0396 h_TkMET_turnon->SetStats(0);
0397 h_TkMET_turnon->SetTitle("; Gen-level MET [GeV]; Efficiency");
0398 h_TkMET_turnon->Divide(h_trueTkMET_num, h_trueMET, 1.0, 1.0, "B");
0399 h_TkMET_turnon->SetMarkerStyle(20);
0400 h_TkMET_turnon->SetMarkerSize(0.7);
0401 h_TkMET_turnon->SetMarkerColor(kBlack);
0402 h_TkMET_turnon->SetLineColor(kBlack);
0403
0404 h_recoTkMET_num->Sumw2();
0405 TH1F* h_recoTkMET_turnon = (TH1F*)h_recoTkMET_num->Clone();
0406 h_recoTkMET_turnon->SetName("recoTkMET_turnon");
0407 h_recoTkMET_turnon->SetStats(0);
0408 h_recoTkMET_turnon->SetTitle("; Gen-level MET [GeV]; Efficiency");
0409 h_recoTkMET_turnon->Divide(h_recoTkMET_num, h_trueMET, 1.0, 1.0, "B");
0410 h_recoTkMET_turnon->SetMarkerStyle(20);
0411 h_recoTkMET_turnon->SetMarkerSize(0.7);
0412 h_recoTkMET_turnon->SetMarkerColor(kRed);
0413 h_recoTkMET_turnon->SetLineColor(kRed);
0414
0415 TCanvas* can = new TCanvas("can", "can", 800, 600);
0416 can->SetGridx();
0417 can->SetGridy();
0418 h_TkMET_turnon->SetMinimum(0);
0419 h_TkMET_turnon->SetMaximum(1.1);
0420 h_TkMET_turnon->Draw();
0421 h_recoTkMET_turnon->Draw("same");
0422 leg->Draw("same");
0423
0424 float trueTkMET_95eff = h_TkMET_turnon->GetBinLowEdge(h_TkMET_turnon->FindFirstBinAbove(0.95));
0425 float recoTkMET_95eff = h_recoTkMET_turnon->GetBinLowEdge(h_recoTkMET_turnon->FindFirstBinAbove(0.95));
0426
0427 std::cout << "Online (true) Track MET at " << trueTkMET_rate << " GeV = " << trueTkMET_95eff << " GeV (offline)"
0428 << std::endl;
0429 std::cout << "Online (reco) Track MET at " << recoTkMET_rate << " GeV = " << recoTkMET_95eff << " GeV (offline)"
0430 << std::endl;
0431
0432 can->SaveAs("trkMET_ttbarPU200_TurnOnPlot.root");
0433 }