L1TTrackMatch/test/trkMET_rate.C

0001 //To find the rate from a (very very large) min bias sample,
0002 //Runs off of the track object ntuplizer
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   // tree->Add("/uscms_data/d3/csavard/working_dir/for_Emily/CMSSW_11_1_3/src/L1Trigger/TrackFindingTracklet/test/crab_output/crab_nu_gun/results/NuGun_PU200_TDR.root");
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   // define leafs & branches
0040   float trueMET = 0;
0041   float trueTkMET = 0;
0042   float trkMET = 0;
0043   vector<float>* pv_L1reco;
0044   //gen particles
0045   vector<float>* gen_pt;
0046   vector<float>* gen_phi;
0047   vector<int>* gen_pdgid;
0048
0049   // tracking particles
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   // vector<int>*   tp_nLayers;
0060   vector<int>* tp_eventid;
0061   // vector<int>*   tp_signal;
0062   vector<int>* tp_charge;
0063
0064   // *L1 track* properties, for tracking particles matched to a L1 track
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   // all L1 tracks
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   // TBranch* b_tp_nLayers;
0101   TBranch* b_tp_nstub;
0102   TBranch* b_tp_eventid;
0103   // TBranch* b_tp_signal;
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   // tp_nLayers = 0;
0150   tp_nstub = 0;
0151   tp_eventid = 0;
0152   // tp_signal = 0;
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   // tree->SetBranchAddress("tp_nLayers", &tp_nLayers, &b_tp_nLayers);
0196   tree->SetBranchAddress("tp_nstub", &tp_nstub, &b_tp_nstub);
0197   tree->SetBranchAddress("tp_eventid", &tp_eventid, &b_tp_eventid);
0198   // tree->SetBranchAddress("tp_signal",    &tp_signal,    &b_tp_signal);
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   //Need trkMET from tracking particles, and trkMET from tracks
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);  //using tracking particles
0235   TH1D* recoTkMET_thresh = new TH1D(
0236       "recoTkMET_thresh", "recoTkMET_thresh;recoTkMET threshold (GeV);Occurrences", numBins, 0, numBins);  //using tracks
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   // event loop
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;  //grab these from ntuple
0251     float recoTkMET_ntuple = 0;  //grab these from ntuple
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       // tracking particle loop
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         // kinematic cuts
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         // if (this_tp_signal!=0) continue;
0280         // if (tp_dxy->at(it) > 1) continue;
0281         if (tp_charge->at(it) == 0)
0282           continue;
0283
0284         float deltaZ_cut = 3.0;  // cuts out PU
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       }  // end of tracking particle loop
0304
0305       trueTkMET_calc = sqrt(trueTkMETx * trueTkMETx + trueTkMETy * trueTkMETy);
0306       trueTkMET_ntuple = trueTkMET_calc;
0307     }  //re-run trueTkMET
0308
0309     if (rerun_recoTkMET) {  //also useful for checking different track quality cuts
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;  // cuts out PU
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       }  //end loop over tracks
0351
0352       recoTkMET_calc = sqrt(recoTkMETx * recoTkMETx + recoTkMETy * recoTkMETy);
0353       recoTkMET_ntuple = recoTkMET_calc;
0354     }  //re-run reco tkMET
0355
0356     trueTkMET_thresh->Fill(trueTkMET_ntuple);
0357     recoTkMET_thresh->Fill(recoTkMET_ntuple);
0358   }  // end event loop
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   // fCumulative_true->Rebin(2); fCumulative_reco->Rebin(2);
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 }