File indexing completed on 2024-04-06 12:21:20
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 TH1D* GetCumulative(TH1D* plot);
0028
0029 void trkMET_rate() {
0030 TChain* tree = new TChain("L1TrackNtuple/eventTree");
0031
0032 tree->Add("CheckingJets_CMSSW11_CMS.root");
0033
0034 if (tree->GetEntries() == 0) {
0035 cout << "File doesn't exist or is empty, returning..." << endl;
0036 return;
0037 }
0038
0039
0040 float trueMET = 0;
0041 float trueTkMET = 0;
0042 float trkMET = 0;
0043 vector<float>* pv_L1reco;
0044
0045 vector<float>* gen_pt;
0046 vector<float>* gen_phi;
0047 vector<int>* gen_pdgid;
0048
0049
0050 vector<float>* tp_pt;
0051 vector<float>* tp_eta;
0052 vector<float>* tp_phi;
0053 vector<float>* tp_dxy;
0054 vector<float>* tp_z0;
0055 vector<float>* tp_d0;
0056 vector<int>* tp_pdgid;
0057 vector<int>* tp_nmatch;
0058 vector<int>* tp_nstub;
0059
0060 vector<int>* tp_eventid;
0061
0062 vector<int>* tp_charge;
0063
0064
0065 vector<float>* matchtrk_pt;
0066 vector<float>* matchtrk_eta;
0067 vector<float>* matchtrk_phi;
0068 vector<float>* matchtrk_d0;
0069 vector<float>* matchtrk_z0;
0070 vector<float>* matchtrk_chi2;
0071 vector<float>* matchtrk_MVA1;
0072 vector<float>* matchtrk_bendchi2;
0073 vector<int>* matchtrk_nstub;
0074 vector<int>* matchtrk_seed;
0075
0076
0077 vector<float>* trk_pt;
0078 vector<float>* trk_eta;
0079 vector<float>* trk_phi;
0080 vector<float>* trk_z0;
0081 vector<float>* trk_chi2;
0082 vector<float>* trk_MVA1;
0083 vector<float>* trk_bendchi2;
0084 vector<int>* trk_nstub;
0085 vector<int>* trk_seed;
0086 vector<int>* trk_fake;
0087
0088 TBranch* b_gen_pt;
0089 TBranch* b_gen_phi;
0090 TBranch* b_gen_pdgid;
0091
0092 TBranch* b_tp_pt;
0093 TBranch* b_tp_eta;
0094 TBranch* b_tp_phi;
0095 TBranch* b_tp_dxy;
0096 TBranch* b_tp_z0;
0097 TBranch* b_tp_d0;
0098 TBranch* b_tp_pdgid;
0099 TBranch* b_tp_nmatch;
0100
0101 TBranch* b_tp_nstub;
0102 TBranch* b_tp_eventid;
0103
0104 TBranch* b_tp_charge;
0105
0106 TBranch* b_matchtrk_pt;
0107 TBranch* b_matchtrk_eta;
0108 TBranch* b_matchtrk_phi;
0109 TBranch* b_matchtrk_d0;
0110 TBranch* b_matchtrk_z0;
0111 TBranch* b_matchtrk_chi2;
0112 TBranch* b_matchtrk_MVA1;
0113 TBranch* b_matchtrk_bendchi2;
0114 TBranch* b_matchtrk_nstub;
0115 TBranch* b_matchtrk_seed;
0116
0117 TBranch* b_trk_pt;
0118 TBranch* b_trk_eta;
0119 TBranch* b_trk_phi;
0120 TBranch* b_trk_z0;
0121 TBranch* b_trk_chi2;
0122 TBranch* b_trk_MVA1;
0123 TBranch* b_trk_bendchi2;
0124 TBranch* b_trk_nstub;
0125 TBranch* b_trk_seed;
0126 TBranch* b_trk_fake;
0127
0128 TBranch* b_pv_L1reco;
0129 TBranch* b_trueMET;
0130 TBranch* b_trueTkMET;
0131 TBranch* b_trkMET;
0132
0133 trueMET = 0;
0134 trueTkMET = 0;
0135 trkMET = 0;
0136
0137 gen_pt = 0;
0138 gen_phi = 0;
0139 gen_pdgid = 0;
0140
0141 tp_pt = 0;
0142 tp_eta = 0;
0143 tp_phi = 0;
0144 tp_dxy = 0;
0145 tp_z0 = 0;
0146 tp_d0 = 0;
0147 tp_pdgid = 0;
0148 tp_nmatch = 0;
0149
0150 tp_nstub = 0;
0151 tp_eventid = 0;
0152
0153 tp_charge = 0;
0154
0155 matchtrk_pt = 0;
0156 matchtrk_eta = 0;
0157 matchtrk_phi = 0;
0158 matchtrk_d0 = 0;
0159 matchtrk_z0 = 0;
0160 matchtrk_chi2 = 0;
0161 matchtrk_MVA1 = 0;
0162 matchtrk_bendchi2 = 0;
0163 matchtrk_nstub = 0;
0164 matchtrk_seed = 0;
0165
0166 trk_pt = 0;
0167 trk_eta = 0;
0168 trk_phi = 0;
0169 trk_z0 = 0;
0170 trk_chi2 = 0;
0171 trk_MVA1 = 0;
0172 trk_bendchi2 = 0;
0173 trk_nstub = 0;
0174 trk_seed = 0;
0175 trk_fake = 0;
0176
0177 pv_L1reco = 0;
0178
0179 tree->SetBranchAddress("pv_L1reco", &pv_L1reco, &b_pv_L1reco);
0180 tree->SetBranchAddress("trueMET", &trueMET, &b_trueMET);
0181 tree->SetBranchAddress("trueTkMET", &trueTkMET, &b_trueTkMET);
0182 tree->SetBranchAddress("trkMET", &trkMET, &b_trkMET);
0183 tree->SetBranchAddress("gen_pt", &gen_pt, &b_gen_pt);
0184 tree->SetBranchAddress("gen_phi", &gen_phi, &b_gen_phi);
0185 tree->SetBranchAddress("gen_pdgid", &gen_pdgid, &b_gen_pdgid);
0186
0187 tree->SetBranchAddress("tp_pt", &tp_pt, &b_tp_pt);
0188 tree->SetBranchAddress("tp_eta", &tp_eta, &b_tp_eta);
0189 tree->SetBranchAddress("tp_phi", &tp_phi, &b_tp_phi);
0190 tree->SetBranchAddress("tp_dxy", &tp_dxy, &b_tp_dxy);
0191 tree->SetBranchAddress("tp_z0", &tp_z0, &b_tp_z0);
0192 tree->SetBranchAddress("tp_d0", &tp_d0, &b_tp_d0);
0193 tree->SetBranchAddress("tp_pdgid", &tp_pdgid, &b_tp_pdgid);
0194 tree->SetBranchAddress("tp_nmatch", &tp_nmatch, &b_tp_nmatch);
0195
0196 tree->SetBranchAddress("tp_nstub", &tp_nstub, &b_tp_nstub);
0197 tree->SetBranchAddress("tp_eventid", &tp_eventid, &b_tp_eventid);
0198
0199 tree->SetBranchAddress("tp_charge", &tp_charge, &b_tp_charge);
0200
0201 tree->SetBranchAddress("matchtrk_pt", &matchtrk_pt, &b_matchtrk_pt);
0202 tree->SetBranchAddress("matchtrk_eta", &matchtrk_eta, &b_matchtrk_eta);
0203 tree->SetBranchAddress("matchtrk_phi", &matchtrk_phi, &b_matchtrk_phi);
0204 tree->SetBranchAddress("matchtrk_d0", &matchtrk_d0, &b_matchtrk_d0);
0205 tree->SetBranchAddress("matchtrk_z0", &matchtrk_z0, &b_matchtrk_z0);
0206 tree->SetBranchAddress("matchtrk_chi2", &matchtrk_chi2, &b_matchtrk_chi2);
0207 tree->SetBranchAddress("matchtrk_MVA1", &matchtrk_MVA1, &b_matchtrk_MVA1);
0208 tree->SetBranchAddress("matchtrk_bendchi2", &matchtrk_bendchi2, &b_matchtrk_bendchi2);
0209 tree->SetBranchAddress("matchtrk_nstub", &matchtrk_nstub, &b_matchtrk_nstub);
0210
0211 tree->SetBranchAddress("trk_pt", &trk_pt, &b_trk_pt);
0212 tree->SetBranchAddress("trk_eta", &trk_eta, &b_trk_eta);
0213 tree->SetBranchAddress("trk_phi", &trk_phi, &b_trk_phi);
0214 tree->SetBranchAddress("trk_z0", &trk_z0, &b_trk_z0);
0215 tree->SetBranchAddress("trk_chi2", &trk_chi2, &b_trk_chi2);
0216 tree->SetBranchAddress("trk_MVA1", &trk_MVA1, &b_trk_MVA1);
0217 tree->SetBranchAddress("trk_bendchi2", &trk_bendchi2, &b_trk_bendchi2);
0218 tree->SetBranchAddress("trk_nstub", &trk_nstub, &b_trk_nstub);
0219 tree->SetBranchAddress("trk_fake", &trk_fake, &b_trk_fake);
0220
0221
0222
0223 float numBins = 100.0;
0224 float cutoff = 200.0 - 0.1;
0225 float TP_minPt = 2.0;
0226 float TP_maxEta = 2.4;
0227 float chi2dof_cut = 10.0;
0228 float bendchi2_cut = 2.2;
0229
0230 TH1D* trueTkMET_thresh = new TH1D("trueTkMET_thresh",
0231 "trueTkMET_thresh;trueTkMET threshold (GeV);Occurrences",
0232 numBins,
0233 0,
0234 numBins);
0235 TH1D* recoTkMET_thresh = new TH1D(
0236 "recoTkMET_thresh", "recoTkMET_thresh;recoTkMET threshold (GeV);Occurrences", numBins, 0, numBins);
0237
0238 bool rerun_trueTkMET = false;
0239 bool rerun_recoTkMET = false;
0240
0241 int nevt = tree->GetEntries();
0242 cout << "number of events = " << nevt << endl;
0243
0244
0245 for (int i = 0; i < nevt; i++) {
0246 tree->GetEntry(i);
0247 if (i % 10000 == 0)
0248 std::cout << i << "/" << nevt << std::endl;
0249
0250 float trueTkMET_ntuple = 0;
0251 float recoTkMET_ntuple = 0;
0252
0253 if (!rerun_trueTkMET)
0254 trueTkMET_ntuple = trueTkMET;
0255 if (!rerun_recoTkMET)
0256 recoTkMET_ntuple = trkMET;
0257
0258 if (rerun_trueTkMET) {
0259 float trueTkMET_calc = 0;
0260 float trueTkMETx = 0;
0261 float trueTkMETy = 0;
0262
0263 for (int it = 0; it < (int)tp_pt->size(); it++) {
0264 float this_tp_pt = tp_pt->at(it);
0265 float this_tp_eta = tp_eta->at(it);
0266 float this_tp_phi = tp_phi->at(it);
0267 int this_tp_signal = tp_eventid->at(it);
0268 float deltaZ = std::abs(tp_z0->at(it) - pv_L1reco->at(0));
0269
0270
0271 if (tp_pt->at(it) < TP_minPt)
0272 continue;
0273 if (std::abs(this_tp_eta) > TP_maxEta)
0274 continue;
0275 if (tp_nstub->at(it) < 4)
0276 continue;
0277 if (std::abs(tp_z0->at(it)) > 15)
0278 continue;
0279
0280
0281 if (tp_charge->at(it) == 0)
0282 continue;
0283
0284 float deltaZ_cut = 3.0;
0285 if (std::abs(this_tp_eta) >= 0 && std::abs(this_tp_eta) < 0.7)
0286 deltaZ_cut = 0.4;
0287 else if (std::abs(this_tp_eta) >= 0.7 && std::abs(this_tp_eta) < 1.0)
0288 deltaZ_cut = 0.6;
0289 else if (std::abs(this_tp_eta) >= 1.0 && std::abs(this_tp_eta) < 1.2)
0290 deltaZ_cut = 0.76;
0291 else if (std::abs(this_tp_eta) >= 1.2 && std::abs(this_tp_eta) < 1.6)
0292 deltaZ_cut = 1.0;
0293 else if (std::abs(this_tp_eta) >= 1.6 && std::abs(this_tp_eta) < 2.0)
0294 deltaZ_cut = 1.7;
0295 else if (std::abs(this_tp_eta) >= 2.0 && std::abs(this_tp_eta) <= 2.4)
0296 deltaZ_cut = 2.20;
0297
0298 if (deltaZ > deltaZ_cut)
0299 continue;
0300
0301 trueTkMETx += this_tp_pt * cos(this_tp_phi);
0302 trueTkMETy += this_tp_pt * sin(this_tp_phi);
0303 }
0304
0305 trueTkMET_calc = sqrt(trueTkMETx * trueTkMETx + trueTkMETy * trueTkMETy);
0306 trueTkMET_ntuple = trueTkMET_calc;
0307 }
0308
0309 if (rerun_recoTkMET) {
0310 float recoTkMET_calc = 0;
0311 float recoTkMETx = 0;
0312 float recoTkMETy = 0;
0313
0314 for (int it = 0; it < (int)trk_pt->size(); it++) {
0315 float thisTrk_pt = trk_pt->at(it);
0316 float thisTrk_eta = trk_eta->at(it);
0317 int thisTrk_nstub = trk_nstub->at(it);
0318 float thisTrk_chi2dof = trk_chi2->at(it) / (2 * thisTrk_nstub - 4);
0319 float thisTrk_bendchi2 = trk_bendchi2->at(it);
0320 float thisTrk_phi = trk_phi->at(it);
0321 float thisTrk_MVA = trk_MVA1->at(it);
0322 int thisTrk_fake = trk_fake->at(it);
0323 float thisTrk_z0 = trk_z0->at(it);
0324 float deltaZ = thisTrk_z0 - pv_L1reco->at(0);
0325
0326 float deltaZ_cut = 3.0;
0327 if (std::abs(thisTrk_eta) >= 0 && std::abs(thisTrk_eta) < 0.7)
0328 deltaZ_cut = 0.4;
0329 else if (std::abs(thisTrk_eta) >= 0.7 && std::abs(thisTrk_eta) < 1.0)
0330 deltaZ_cut = 0.6;
0331 else if (std::abs(thisTrk_eta) >= 1.0 && std::abs(thisTrk_eta) < 1.2)
0332 deltaZ_cut = 0.76;
0333 else if (std::abs(thisTrk_eta) >= 1.2 && std::abs(thisTrk_eta) < 1.6)
0334 deltaZ_cut = 1.0;
0335 else if (std::abs(thisTrk_eta) >= 1.6 && std::abs(thisTrk_eta) < 2.0)
0336 deltaZ_cut = 1.7;
0337 else if (std::abs(thisTrk_eta) >= 2.0 && std::abs(thisTrk_eta) <= 2.4)
0338 deltaZ_cut = 2.20;
0339
0340 if (thisTrk_pt < TP_minPt || std::abs(thisTrk_eta) > TP_maxEta || thisTrk_nstub < 4 ||
0341 std::abs(thisTrk_z0) > 15)
0342 continue;
0343 if (std::abs(deltaZ) > deltaZ_cut)
0344 continue;
0345
0346 if (thisTrk_chi2dof < chi2dof_cut && thisTrk_bendchi2 < bendchi2_cut) {
0347 recoTkMETx += thisTrk_pt * cos(thisTrk_phi);
0348 recoTkMETy += thisTrk_pt * sin(thisTrk_phi);
0349 }
0350 }
0351
0352 recoTkMET_calc = sqrt(recoTkMETx * recoTkMETx + recoTkMETy * recoTkMETy);
0353 recoTkMET_ntuple = recoTkMET_calc;
0354 }
0355
0356 trueTkMET_thresh->Fill(trueTkMET_ntuple);
0357 recoTkMET_thresh->Fill(recoTkMET_ntuple);
0358 }
0359
0360
0361 trueTkMET_thresh->SetStats(0);
0362 recoTkMET_thresh->SetStats(0);
0363 trueTkMET_thresh->SetLineColor(kBlack);
0364 recoTkMET_thresh->SetLineColor(kRed);
0365 trueTkMET_thresh->SetTitle("");
0366 recoTkMET_thresh->SetTitle("");
0367
0368 TLegend* leg = new TLegend(0.48, 0.68, 0.88, 0.88);
0369 leg->SetBorderSize(0);
0370 leg->SetFillStyle(0);
0371 leg->SetTextSize(0.04);
0372 leg->SetTextFont(42);
0373 leg->AddEntry(trueTkMET_thresh, "trueTkMET", "l");
0374 leg->AddEntry(recoTkMET_thresh, "recoTkMET", "l");
0375
0376 TH1D* fCumulative_true = GetCumulative(trueTkMET_thresh);
0377 TH1D* fCumulative_reco = GetCumulative(recoTkMET_thresh);
0378 fCumulative_true->Scale(4.0e3 / fCumulative_true->GetBinContent(1));
0379 fCumulative_reco->Scale(4.0e3 / fCumulative_reco->GetBinContent(1));
0380
0381 fCumulative_true->SetTitle("MinBias Events PU 200; Track MET [GeV]; L1 Rate [kHz]");
0382 fCumulative_true->SetMarkerSize(0.7);
0383 fCumulative_true->SetMarkerStyle(20);
0384 fCumulative_true->SetMarkerColor(fCumulative_true->GetLineColor());
0385
0386 fCumulative_reco->SetMarkerSize(0.7);
0387 fCumulative_reco->SetMarkerStyle(20);
0388 fCumulative_reco->SetMarkerColor(fCumulative_reco->GetLineColor());
0389
0390
0391
0392 TCanvas* can = new TCanvas("can", "can");
0393 can->SetGridx();
0394 can->SetGridy();
0395 can->SetLogy();
0396
0397 fCumulative_true->GetYaxis()->SetRangeUser(1E-01, 1E05);
0398
0399 fCumulative_true->Draw();
0400 fCumulative_reco->Draw("same");
0401 leg->Draw("same");
0402
0403 TLine* line_rate = new TLine(0, 35, numBins, 35);
0404 line_rate->SetLineColor(kBlack);
0405 line_rate->SetLineWidth(2);
0406 line_rate->Draw("SAME");
0407
0408 std::cout << "Track MET Threshold at 35kHz "
0409 << fCumulative_reco->GetBinLowEdge(fCumulative_reco->FindLastBinAbove(35)) << std::endl;
0410
0411 can->SaveAs("trkMET_minBias_RatePlot.root");
0412 }
0413
0414 TH1D* GetCumulative(TH1D* plot) {
0415 std::string newName = Form("cumulative_%s", plot->GetName());
0416 TH1D* temp = (TH1D*)plot->Clone(newName.c_str());
0417 temp->SetDirectory(0);
0418 for (int i = 0; i < plot->GetNbinsX() + 1; i++) {
0419 double content(0.0), error2(0.0);
0420 for (int j = i; j < plot->GetNbinsX() + 1; j++) {
0421 content += plot->GetBinContent(j);
0422 error2 += plot->GetBinError(j) * plot->GetBinError(j);
0423 }
0424 temp->SetBinContent(i, content);
0425 temp->SetBinError(i, TMath::Sqrt(error2));
0426 }
0427 return temp;
0428 }