L1TTrackMatch/test/trkMET_turnOn.C

0001 //To make basic turn-on curve for trkMET
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
0028 void trkMET_turnOn() {
0029   TChain* tree = new TChain("L1TrackNtuple/eventTree");
0030   // tree->Add("/uscms_data/d3/csavard/working_dir/for_Emily/CMSSW_11_1_3/src/L1Trigger/TrackFindingTracklet/test/crab_output/crab_ttbar_pu200/results/TTbar_PU200_D49.root");
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   // define leafs & branches
0039   float trueMET = 0;
0040   float trueTkMET = 0;
0041   float trkMET = 0;
0042   vector<float>* pv_L1reco;
0043   //gen particles
0044   vector<float>* gen_pt;
0045   vector<float>* gen_phi;
0046   vector<int>* gen_pdgid;
0047
0048   // tracking particles
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   // vector<int>*   tp_nLayers;
0059   vector<int>* tp_eventid;
0060   // vector<int>*   tp_signal;
0061   vector<int>* tp_charge;
0062
0063   // *L1 track* properties, for tracking particles matched to a L1 track
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   // all L1 tracks
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   // TBranch* b_tp_nLayers;
0100   TBranch* b_tp_nstub;
0101   TBranch* b_tp_eventid;
0102   // TBranch* b_tp_signal;
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   // tp_nLayers = 0;
0149   tp_nstub = 0;
0150   tp_eventid = 0;
0151   // tp_signal = 0;
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   // tree->SetBranchAddress("tp_nLayers", &tp_nLayers, &b_tp_nLayers);
0195   tree->SetBranchAddress("tp_nstub", &tp_nstub, &b_tp_nstub);
0196   tree->SetBranchAddress("tp_eventid", &tp_eventid, &b_tp_eventid);
0197   // tree->SetBranchAddress("tp_signal",    &tp_signal,    &b_tp_signal);
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   //Need trkMET from tracking particles, and trkMET from tracks
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   //From rate code
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   // event loop
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;      //grab these from ntuple
0254     float trueTkMET_ntuple = trueTkMET;  //grab these from ntuple
0255     float recoTkMET_ntuple = trkMET;     //grab these from ntuple
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)) {  //only gen parts saved are status==1
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     }  //re-run true MET (gen particles)
0276
0277     if (rerun_trueTkMET) {
0278       float trueTkMET_calc = 0;
0279       float trueTkMETx = 0;
0280       float trueTkMETy = 0;
0281       // tracking particle loop
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         // kinematic cuts
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;  // hard interaction only
0299         // if (tp_dxy->at(it) > 1) continue;
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       }  // end of tracking particle loop
0306       trueTkMET_calc = sqrt(trueTkMETx * trueTkMETx + trueTkMETy * trueTkMETy);
0307       trueTkMET_ntuple = trueTkMET_calc;
0308     }  // re-run true trk MET
0309
0310     if (rerun_recoTkMET) {  //also useful for checking performance of new track quality cuts
0311       float recoTkMET_calc = 0;
0312       float recoTkMETx = 0;
0313       float recoTkMETy = 0;
0314       // loop over all tracks
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;  // cuts out PU
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         // if (thisTrk_pt>200.0) thisTrk_pt = 200.0;
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       }  //end loop over tracks
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     }  //re-run reco trk MET
0359
0360     h_trueMET->Fill(trueMET_ntuple);
0361     h_trueTkMET->Fill(trueTkMET_ntuple);
0362     h_recoTkMET->Fill(recoTkMET_ntuple);
0363
0364     //get rate cuts from Rate Plot script
0365     //if reco MET passes cut, fill with trueMET
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   }  // end event loop
0371
0372   //Draw
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   // calculate the recoTkMET trigger turn on curve
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 }