TrackFindingTracklet/test/L1TrackQualityPlot.C

0001 // ----------------------------------------------------------------------------------------------------------------
0002 // Basic example ROOT script for making tracking performance plots using the ntuples produced by L1TrackNtupleMaker.cc
0003 //
0004 // e.g. in ROOT do: [0] .L L1TrackQualityPlot.C++
0005 //                  [1] L1TrackQualityPlot("TTbar_PU200_D76")
0006 //
0007 // By Claire Savard, July 2020
0008 // Based off of L1TrackNtuplePlot.C
0009 // ----------------------------------------------------------------------------------------------------------------
0010
0011 #include "TROOT.h"
0012 #include "TStyle.h"
0013 #include "TLatex.h"
0014 #include "TFile.h"
0015 #include "TTree.h"
0016 #include "TChain.h"
0017 #include "TBranch.h"
0018 #include "TLeaf.h"
0019 #include "TCanvas.h"
0020 #include "TLegend.h"
0021 #include "TH1.h"
0022 #include "TH2.h"
0023 #include "TF1.h"
0024 #include "TProfile.h"
0025 #include "TProfile2D.h"
0026 #include "TMath.h"
0027 #include "TGraph.h"
0028 #include "TError.h"
0029 #include "TGraphErrors.h"
0030 #include "TGraphPainter.h"
0031 #include "TSystem.h"
0032
0033 #include <iostream>
0034 #include <string>
0035 #include <vector>
0036 #include <algorithm>
0037 #include <cmath>
0038
0039 using namespace std;
0040
0041 void SetPlotStyle();
0042
0043 // ----------------------------------------------------------------------------------------------------------------
0044 // Main script
0045 // ----------------------------------------------------------------------------------------------------------------
0046
0047 void L1TrackQualityPlot(TString type, TString type_dir = "", TString treeName = "") {
0048   // type:              this is the name of the input file you want to process (minus ".root" extension)
0049   // type_dir:          this is the directory containing the input file you want to process. Note that this must end with a "/", as in "EventSets/"
0050
0051   gROOT->SetBatch();
0052   gErrorIgnoreLevel = kWarning;
0053
0054   SetPlotStyle();
0055   gSystem->mkdir("MVA_plots");
0056
0057   // ----------------------------------------------------------------------------------------------------------------
0058   // define input options
0059
0060   // these are the LOOSE cuts, baseline scenario for efficiency and rate plots ==> configure as appropriate
0061   int L1Tk_minNstub = 4;
0062   float L1Tk_maxChi2 = 999999;
0063   float L1Tk_maxChi2dof = 999999.;
0064
0065   // ----------------------------------------------------------------------------------------------------------------
0066   // read ntuples
0067   TChain* tree = new TChain("L1TrackNtuple" + treeName + "/eventTree");
0068   tree->Add(type_dir + type + ".root");
0069
0070   if (tree->GetEntries() == 0) {
0071     cout << "File doesn't exist or is empty, returning..."
0072          << endl;  //cout's kept in this file as it is an example standalone plotting script, not running in central CMSSW
0073     return;
0074   }
0075
0076   // ----------------------------------------------------------------------------------------------------------------
0077   // define leafs & branches
0078   // all L1 tracks
0079   vector<float>* trk_pt;
0080   vector<float>* trk_eta;
0081   vector<float>* trk_phi;
0082   vector<float>* trk_chi2;
0083   vector<float>* trk_chi2rphi;
0084   vector<float>* trk_chi2rz;
0085   vector<int>* trk_nstub;
0086   vector<int>* trk_lhits;
0087   vector<int>* trk_dhits;
0088   vector<int>* trk_seed;
0089   vector<int>* trk_hitpattern;
0090   vector<unsigned int>* trk_phiSector;
0091   vector<int>* trk_fake;
0092   vector<int>* trk_genuine;
0093   vector<int>* trk_loose;
0094   vector<float>* trk_MVA1;
0095   vector<float>* trk_matchtp_pdgid;
0096
0097   TBranch* b_trk_pt;
0098   TBranch* b_trk_eta;
0099   TBranch* b_trk_phi;
0100   TBranch* b_trk_chi2;
0101   TBranch* b_trk_chi2rphi;
0102   TBranch* b_trk_chi2rz;
0103   TBranch* b_trk_nstub;
0104   TBranch* b_trk_lhits;
0105   TBranch* b_trk_dhits;
0106   TBranch* b_trk_phiSector;
0107   TBranch* b_trk_seed;
0108   TBranch* b_trk_hitpattern;
0109   TBranch* b_trk_fake;
0110   TBranch* b_trk_genuine;
0111   TBranch* b_trk_loose;
0112   TBranch* b_trk_MVA1;
0113   TBranch* b_trk_matchtp_pdgid;
0114
0115   trk_pt = 0;
0116   trk_eta = 0;
0117   trk_phi = 0;
0118   trk_chi2 = 0;
0119   trk_chi2rphi = 0;
0120   trk_chi2rz = 0;
0121   trk_nstub = 0;
0122   trk_lhits = 0;
0123   trk_dhits = 0;
0124   trk_phiSector = 0;
0125   trk_seed = 0;
0126   trk_hitpattern = 0;
0127   trk_fake = 0;
0128   trk_genuine = 0;
0129   trk_loose = 0;
0130   trk_MVA1 = 0;
0131   trk_matchtp_pdgid = 0;
0132
0133   tree->SetBranchAddress("trk_pt", &trk_pt, &b_trk_pt);
0134   tree->SetBranchAddress("trk_eta", &trk_eta, &b_trk_eta);
0135   tree->SetBranchAddress("trk_phi", &trk_phi, &b_trk_phi);
0136   tree->SetBranchAddress("trk_chi2", &trk_chi2, &b_trk_chi2);
0137   tree->SetBranchAddress("trk_chi2rphi", &trk_chi2rphi, &b_trk_chi2rphi);
0138   tree->SetBranchAddress("trk_chi2rz", &trk_chi2rz, &b_trk_chi2rz);
0139   tree->SetBranchAddress("trk_nstub", &trk_nstub, &b_trk_nstub);
0140   tree->SetBranchAddress("trk_lhits", &trk_lhits, &b_trk_lhits);
0141   tree->SetBranchAddress("trk_dhits", &trk_dhits, &b_trk_dhits);
0142   tree->SetBranchAddress("trk_phiSector", &trk_phiSector, &b_trk_phiSector);
0143   tree->SetBranchAddress("trk_seed", &trk_seed, &b_trk_seed);
0144   tree->SetBranchAddress("trk_hitpattern", &trk_hitpattern, &b_trk_hitpattern);
0145   tree->SetBranchAddress("trk_fake", &trk_fake, &b_trk_fake);
0146   tree->SetBranchAddress("trk_genuine", &trk_genuine, &b_trk_genuine);
0147   tree->SetBranchAddress("trk_loose", &trk_loose, &b_trk_loose);
0148   tree->SetBranchAddress("trk_MVA1", &trk_MVA1, &b_trk_MVA1);
0149   tree->SetBranchAddress("trk_matchtp_pdgid", &trk_matchtp_pdgid, &b_trk_matchtp_pdgid);
0150
0151   // ----------------------------------------------------------------------------------------------------------------
0152   // histograms
0153   // ----------------------------------------------------------------------------------------------------------------
0154
0155   TH1F* h_trk_MVA1 = new TH1F("trk_MVA1", "; MVA1; L1 tracks", 50, 0, 1);
0156
0157   TH1F* h_trk_MVA1_real = new TH1F("trk_MVA1_real", ";MVA1; L1 tracks", 50, 0, 1);
0158   h_trk_MVA1_real->SetLineColor(3);
0159   TH1F* h_trk_MVA1_fake = new TH1F("trk_MVA1_fake", ";MVA1; L1 tracks", 50, 0, 1);
0160   h_trk_MVA1_fake->SetLineColor(4);
0161
0162   // ----------------------------------------------------------------------------------------------------------------
0163   //        * * * * *     S T A R T   O F   A C T U A L   R U N N I N G   O N   E V E N T S     * * * * *
0164   // ----------------------------------------------------------------------------------------------------------------
0165
0166   int nevt = tree->GetEntries();
0167   cout << "number of events = " << nevt << endl;
0168
0169   // ----------------------------------------------------------------------------------------------------------------
0170   // event loop
0171   vector<float> MVA1s;
0172   vector<float> fakes;
0173   vector<float> etas;
0174   vector<float> pts;
0175   vector<int> pdgids;
0176   for (int i = 0; i < nevt; i++) {
0177     tree->GetEntry(i, 0);
0178
0179     for (int it = 0; it < (int)trk_pt->size(); it++) {
0180       // ----------------------------------------------------------------------------------------------------------------
0181       // track properties
0182
0183       float MVA1 = trk_MVA1->at(it);
0184       float fake = trk_fake->at(it);
0185       float eta = trk_eta->at(it);
0186       float pt = trk_pt->at(it);
0187       float pdgid = trk_matchtp_pdgid->at(it);
0188
0189       MVA1s.push_back(MVA1);
0190       fakes.push_back(fake);
0191       etas.push_back(eta);
0192       pts.push_back(pt);
0193       pdgids.push_back(pdgid);
0194
0195       h_trk_MVA1->Fill(MVA1);
0196       if (fake == 1.)
0197         h_trk_MVA1_real->Fill(MVA1);
0198       else if (fake == 0.)
0199         h_trk_MVA1_fake->Fill(MVA1);
0200     }
0201   }
0202
0203   // -------------------------------------------------------------------------------------------
0204   // create ROC curve
0205   // ROC = Receiver Operating Characteristic Curve, a plot of True Positive Rate vs False Positive Rate
0206   // TPR = True Positive Rate or Identification efficiency, fraction of real tracks correctly identified as real
0207   // FPR = False Positive Rate or Fake Rate, fraction of fake tracks incorrectly identified as real
0208   // dt = Decision Threshold or cut on the MVA output, below this identify track as fake, above identify as real
0209   // -------------------------------------------------------------------------------------------
0210
0211   vector<float> TPR, TPR_mu, TPR_el, TPR_had;
0212   vector<float> FPR;
0213   vector<float> dec_thresh;
0214   int n = 100;  //num of entries on ROC curve
0215   for (int i = 0; i < n; i++) {
0216     float dt = (float)i / (n - 1);                   //make sure it starts at (0,0) and ends at (1,1)
0217     float TP = 0, TP_mu = 0, TP_el = 0, TP_had = 0;  //True Positives
0218     float FP = 0;                                    //False Positives
0219     float P = 0, P_mu = 0, P_el = 0, P_had = 0;      //Total Positives
0220     float N = 0;                                     //Total Negatives
0221     for (int k = 0; k < (int)MVA1s.size(); k++) {
0222       if (fakes.at(k)) {
0223         P++;
0224         if (MVA1s.at(k) > dt)
0225           TP++;
0226         if (abs(pdgids.at(k)) == 13) {  //muons
0227           P_mu++;
0228           if (MVA1s.at(k) > dt)
0229             TP_mu++;
0230         } else if (abs(pdgids.at(k)) == 11) {  //electrons
0231           P_el++;
0232           if (MVA1s.at(k) > dt)
0233             TP_el++;
0234         } else if (abs(pdgids.at(k)) > 37 && abs(pdgids.at(k)) != 999) {  //hadrons
0235           P_had++;
0236           if (MVA1s.at(k) > dt)
0237             TP_had++;
0238         }
0239       } else {
0240         N++;
0241         if (MVA1s.at(k) > dt)
0242           FP++;
0243       }
0244     }
0245     TPR.push_back((float)TP / P);
0246     TPR_mu.push_back((float)TP_mu / P_mu);
0247     TPR_el.push_back((float)TP_el / P_el);
0248     TPR_had.push_back((float)TP_had / P_had);
0249     FPR.push_back((float)FP / N);
0250     dec_thresh.push_back(dt);
0251   }
0252
0253   // calculate AUC (Area under the ROC curve)
0254   float AUC = 0., AUC_mu = 0., AUC_el = 0., AUC_had = 0.;
0255   for (int i = 0; i < n - 1; i++) {
0256     AUC += (TPR[i] + TPR[i + 1]) / 2 * (FPR[i] - FPR[i + 1]);
0257     AUC_mu += (TPR_mu[i] + TPR_mu[i + 1]) / 2 * (FPR[i] - FPR[i + 1]);
0258     AUC_el += (TPR_el[i] + TPR_el[i + 1]) / 2 * (FPR[i] - FPR[i + 1]);
0259     AUC_had += (TPR_had[i] + TPR_had[i + 1]) / 2 * (FPR[i] - FPR[i + 1]);
0260   }
0261
0262   TGraph* ROC = new TGraph(n, FPR.data(), TPR.data());
0263   ROC->SetName("ROC");
0264   ROC->SetTitle(("ROC curve (AUC = " + to_string(AUC) + "); FPR; TPR").c_str());
0265   ROC->SetLineWidth(4);
0266
0267   TGraph* ROC_mu = new TGraph(n, FPR.data(), TPR_mu.data());
0268   ROC_mu->SetName("ROC_mu");
0269   ROC_mu->SetTitle(("ROC curve (muons, AUC = " + to_string(AUC_mu) + "); FPR; TPR").c_str());
0270   ROC_mu->SetLineWidth(4);
0271
0272   TGraph* ROC_el = new TGraph(n, FPR.data(), TPR_el.data());
0273   ROC_el->SetName("ROC_el");
0274   ROC_el->SetTitle(("ROC curve (electrons, AUC = " + to_string(AUC_el) + "); FPR; TPR").c_str());
0275   ROC_el->SetLineWidth(4);
0276
0277   TGraph* ROC_had = new TGraph(n, FPR.data(), TPR_had.data());
0278   ROC_had->SetName("ROC_had");
0279   ROC_had->SetTitle(("ROC curve (hadrons, AUC = " + to_string(AUC_had) + "); FPR; TPR").c_str());
0280   ROC_had->SetLineWidth(4);
0281
0282   TGraph* TPR_vs_dt = new TGraph(n, dec_thresh.data(), TPR.data());
0283   TPR_vs_dt->SetName("TPR_vs_dt");
0284   TPR_vs_dt->SetTitle("TPR vs decision threshold; decision thresh.; TPR");
0285   TPR_vs_dt->SetLineColor(3);
0286   TPR_vs_dt->SetLineWidth(4);
0287
0288   TGraph* FPR_vs_dt = new TGraph(n, dec_thresh.data(), FPR.data());
0289   FPR_vs_dt->SetName("FPR_vs_dt");
0290   FPR_vs_dt->SetTitle("FPR vs decision threshold; decision thresh.; FPR");
0291   FPR_vs_dt->SetLineColor(4);
0292   FPR_vs_dt->SetLineWidth(4);
0293
0294   TGraph* TPR_vs_dt_mu = new TGraph(n, dec_thresh.data(), TPR_mu.data());
0295   TPR_vs_dt_mu->SetName("TPR_vs_dt_mu");
0296   TPR_vs_dt_mu->SetTitle("TPR vs decision threshold (muons); decision thresh.; TPR");
0297   TPR_vs_dt_mu->SetLineColor(3);
0298   TPR_vs_dt_mu->SetLineWidth(4);
0299
0300   TGraph* TPR_vs_dt_el = new TGraph(n, dec_thresh.data(), TPR_el.data());
0301   TPR_vs_dt_el->SetName("TPR_vs_dt_el");
0302   TPR_vs_dt_el->SetTitle("TPR vs decision threshold (electrons); decision thresh.; TPR");
0303   TPR_vs_dt_el->SetLineColor(3);
0304   TPR_vs_dt_el->SetLineWidth(4);
0305
0306   TGraph* TPR_vs_dt_had = new TGraph(n, dec_thresh.data(), TPR_had.data());
0307   TPR_vs_dt_had->SetName("TPR_vs_dt_had");
0308   TPR_vs_dt_had->SetTitle("TPR vs decision threshold (hadrons); decision thresh.; TPR");
0309   TPR_vs_dt_had->SetLineColor(3);
0310   TPR_vs_dt_had->SetLineWidth(4);
0311
0312   // -------------------------------------------------------------------------------------------
0313   // create TPR vs. eta and FPR vs. eta
0314   // -------------------------------------------------------------------------------------------
0315
0316   vector<float> TPR_eta, TPR_eta_mu, TPR_eta_el, TPR_eta_had;
0317   vector<float> TPR_eta_err, TPR_eta_err_mu, TPR_eta_err_el, TPR_eta_err_had;
0318   vector<float> FPR_eta, FPR_eta_err;
0319   vector<float> eta_range, eta_range_err;
0320   n = 20;
0321   float eta_low = -2.4;
0322   float eta_high = 2.4;
0323   float eta_temp = eta_low;
0324   float eta_step = (eta_high - eta_low) / n;
0325   float dt = .5;
0326   for (int ct = 0; ct < n; ct++) {
0327     float TP = 0, TP_mu = 0, TP_el = 0, TP_had = 0;
0328     float FP = 0;
0329     float P = 0, P_mu = 0, P_el = 0, P_had = 0;
0330     float N = 0;
0331     for (int k = 0; k < (int)etas.size(); k++) {
0332       if (etas.at(k) > eta_temp && etas.at(k) <= (eta_temp + eta_step)) {
0333         if (fakes.at(k)) {
0334           P++;
0335           if (MVA1s.at(k) > dt)
0336             TP++;
0337           if (abs(pdgids.at(k)) == 13) {  //muons
0338             P_mu++;
0339             if (MVA1s.at(k) > dt)
0340               TP_mu++;
0341           } else if (abs(pdgids.at(k)) == 11) {  //electrons
0342             P_el++;
0343             if (MVA1s.at(k) > dt)
0344               TP_el++;
0345           } else if (abs(pdgids.at(k)) > 37 && abs(pdgids.at(k)) != 999) {  //hadrons
0346             P_had++;
0347             if (MVA1s.at(k) > dt)
0348               TP_had++;
0349           }
0350         } else {
0351           N++;
0352           if (MVA1s.at(k) > dt)
0353             FP++;
0354         }
0355       }
0356     }
0357
0358     //use min function to return 0 if no data filled
0359     TPR_eta.push_back(min(TP / P, P));
0360     TPR_eta_mu.push_back(min(TP_mu / P_mu, P_mu));
0361     TPR_eta_el.push_back(min(TP_el / P_el, P_el));
0362     TPR_eta_had.push_back(min(TP_had / P_had, P_had));
0363     TPR_eta_err.push_back(min((float)sqrt(TP * (P - TP) / pow(P, 3)), P));
0364     TPR_eta_err_mu.push_back(min((float)sqrt(TP_mu * (P_mu - TP_mu) / pow(P_mu, 3)), P_mu));
0365     TPR_eta_err_el.push_back(min((float)sqrt(TP_mu * (P_el - TP_el) / pow(P_el, 3)), P_el));
0366     TPR_eta_err_had.push_back(min((float)sqrt(TP_had * (P_had - TP_had) / pow(P_had, 3)), P_had));
0367
0368     FPR_eta.push_back(min(FP / N, N));
0369     FPR_eta_err.push_back(min((float)sqrt(FP * (N - FP) / pow(N, 3)), N));
0370
0371     //fill eta range
0372     eta_range.push_back(eta_temp + eta_step / 2);
0373     eta_range_err.push_back(eta_step / 2);
0374
0375     eta_temp += eta_step;
0376   }
0377
0378   TGraphErrors* TPR_vs_eta =
0379       new TGraphErrors(n, eta_range.data(), TPR_eta.data(), eta_range_err.data(), TPR_eta_err.data());
0380   TPR_vs_eta->SetName("TPR_vs_eta");
0381   TPR_vs_eta->SetTitle("TPR vs. #eta; #eta; TPR");
0382
0383   TGraphErrors* FPR_vs_eta =
0384       new TGraphErrors(n, eta_range.data(), FPR_eta.data(), eta_range_err.data(), FPR_eta_err.data());
0385   FPR_vs_eta->SetName("FPR_vs_eta");
0386   FPR_vs_eta->SetTitle("FPR vs. #eta; #eta; FPR");
0387
0388   TGraphErrors* TPR_vs_eta_mu =
0389       new TGraphErrors(n, eta_range.data(), TPR_eta_mu.data(), eta_range_err.data(), TPR_eta_err_mu.data());
0390   TPR_vs_eta_mu->SetName("TPR_vs_eta_mu");
0391   TPR_vs_eta_mu->SetTitle("TPR vs. #eta (muons); #eta; TPR");
0392
0393   TGraphErrors* TPR_vs_eta_el =
0394       new TGraphErrors(n, eta_range.data(), TPR_eta_el.data(), eta_range_err.data(), TPR_eta_err_el.data());
0395   TPR_vs_eta_el->SetName("TPR_vs_eta_el");
0396   TPR_vs_eta_el->SetTitle("TPR vs. #eta (electrons); #eta; TPR");
0397
0398   TGraphErrors* TPR_vs_eta_had =
0399       new TGraphErrors(n, eta_range.data(), TPR_eta_had.data(), eta_range_err.data(), TPR_eta_err_had.data());
0400   TPR_vs_eta_had->SetName("TPR_vs_eta_had");
0401   TPR_vs_eta_had->SetTitle("TPR vs. #eta (hadrons); #eta; TPR");
0402
0403   // -------------------------------------------------------------------------------------------
0404   // create TPR vs. pt and FPR vs. pt
0405   // -------------------------------------------------------------------------------------------
0406
0407   vector<float> TPR_pt, TPR_pt_mu, TPR_pt_el, TPR_pt_had;
0408   vector<float> TPR_pt_err, TPR_pt_err_mu, TPR_pt_err_el, TPR_pt_err_had;
0409   vector<float> FPR_pt, FPR_pt_err;
0410   vector<float> pt_range, pt_range_err;
0411   n = 10;
0412   float logpt_low = log10(2);     //set low pt in log
0413   float logpt_high = log10(100);  //set high pt in log
0414   float logpt_temp = logpt_low;
0415   float logpt_step = (logpt_high - logpt_low) / n;
0416   dt = .5;
0417   for (int ct = 0; ct < n; ct++) {
0418     float TP = 0, TP_mu = 0, TP_el = 0, TP_had = 0;
0419     float FP = 0;
0420     float P = 0, P_mu = 0, P_el = 0, P_had = 0;
0421     float N = 0;
0422     for (int k = 0; k < (int)pts.size(); k++) {
0423       if (pts.at(k) > pow(10, logpt_temp) && pts.at(k) <= (pow(10, logpt_temp + logpt_step))) {
0424         if (fakes.at(k)) {
0425           P++;
0426           if (MVA1s.at(k) > dt)
0427             TP++;
0428           if (abs(pdgids.at(k)) == 13) {  //muons
0429             P_mu++;
0430             if (MVA1s.at(k) > dt)
0431               TP_mu++;
0432           } else if (abs(pdgids.at(k)) == 11) {  //electrons
0433             P_el++;
0434             if (MVA1s.at(k) > dt)
0435               TP_el++;
0436           } else if (abs(pdgids.at(k)) > 37 && abs(pdgids.at(k)) != 999) {  //hadrons
0437             P_had++;
0438             if (MVA1s.at(k) > dt)
0439               TP_had++;
0440           }
0441         } else {
0442           N++;
0443           if (MVA1s.at(k) > dt)
0444             FP++;
0445         }
0446       }
0447     }
0448
0449     //use min function to return 0 if no data filled
0450     TPR_pt.push_back(min(TP / P, P));
0451     TPR_pt_mu.push_back(min(TP_mu / P_mu, P_mu));
0452     TPR_pt_el.push_back(min(TP_el / P_el, P_el));
0453     TPR_pt_had.push_back(min(TP_had / P_had, P_had));
0454     TPR_pt_err.push_back(min((float)sqrt(TP * (P - TP) / pow(P, 3)), P));
0455     TPR_pt_err_mu.push_back(min((float)sqrt(TP_mu * (P_mu - TP_mu) / pow(P_mu, 3)), P_mu));
0456     TPR_pt_err_el.push_back(min((float)sqrt(TP_el * (P_el - TP_el) / pow(P_el, 3)), P_el));
0457     TPR_pt_err_had.push_back(min((float)sqrt(TP_had * (P_had - TP_had) / pow(P_had, 3)), P_had));
0458
0459     FPR_pt.push_back(min(FP / N, N));
0460     FPR_pt_err.push_back(min((float)sqrt(FP * (N - FP) / pow(N, 3)), N));
0461
0462     //fill pt range
0463     pt_range.push_back((pow(10, logpt_temp) + pow(10, logpt_temp + logpt_step)) / 2);  //halfway in bin
0464     pt_range_err.push_back((pow(10, logpt_temp + logpt_step) - pow(10, logpt_temp)) / 2);
0465
0466     logpt_temp += logpt_step;
0467   }
0468
0469   TGraphErrors* TPR_vs_pt = new TGraphErrors(n, pt_range.data(), TPR_pt.data(), pt_range_err.data(), TPR_pt_err.data());
0470   TPR_vs_pt->SetName("TPR_vs_pt");
0471   TPR_vs_pt->SetTitle("TPR vs. p_{T}; p_{T}; TPR");
0472
0473   TGraphErrors* FPR_vs_pt = new TGraphErrors(n, pt_range.data(), FPR_pt.data(), pt_range_err.data(), FPR_pt_err.data());
0474   FPR_vs_pt->SetName("FPR_vs_pt");
0475   FPR_vs_pt->SetTitle("FPR vs. p_{T}; p_{T}; FPR");
0476
0477   TGraphErrors* TPR_vs_pt_mu =
0478       new TGraphErrors(n, pt_range.data(), TPR_pt_mu.data(), pt_range_err.data(), TPR_pt_err_mu.data());
0479   TPR_vs_pt_mu->SetName("TPR_vs_pt_mu");
0480   TPR_vs_pt_mu->SetTitle("TPR vs. p_{T} (muons); p_{T}; TPR");
0481
0482   TGraphErrors* TPR_vs_pt_el =
0483       new TGraphErrors(n, pt_range.data(), TPR_pt_el.data(), pt_range_err.data(), TPR_pt_err_el.data());
0484   TPR_vs_pt_el->SetName("TPR_vs_pt_el");
0485   TPR_vs_pt_el->SetTitle("TPR vs. p_{T} (electrons); p_{T}; TPR");
0486
0487   TGraphErrors* TPR_vs_pt_had =
0488       new TGraphErrors(n, pt_range.data(), TPR_pt_had.data(), pt_range_err.data(), TPR_pt_err_had.data());
0489   TPR_vs_pt_had->SetName("TPR_vs_pt_had");
0490   TPR_vs_pt_had->SetTitle("TPR vs. p_{T} (hadrons); p_{T}; TPR");
0491
0492   // -------------------------------------------------------------------------------------------
0493   // output file for histograms and graphs
0494   // -------------------------------------------------------------------------------------------
0495
0496   TFile* fout = new TFile(type_dir + "MVAoutput_" + type + treeName + ".root", "recreate");
0497   TCanvas c;
0498
0499   // -------------------------------------------------------------------------------------------
0500   // draw and save plots
0501   // -------------------------------------------------------------------------------------------
0502
0503   h_trk_MVA1->Draw();
0504   h_trk_MVA1->Write();
0505   c.SaveAs("MVA_plots/trk_MVA.pdf");
0506
0507   h_trk_MVA1_real->Draw();
0508   h_trk_MVA1_fake->Draw("same");
0509   h_trk_MVA1_fake->SetTitle("Performance vs. decision threshold; decision thresh.; performance measure");
0510   TLegend* leg1 = new TLegend();
0511   leg1->AddEntry(h_trk_MVA1_real, "real", "l");
0512   leg1->AddEntry(h_trk_MVA1_fake, "fake", "l");
0513   leg1->Draw("same");
0514   c.Write("trk_MVA_rf");
0515   c.SaveAs("MVA_plots/trk_MVA_rf.pdf");
0516
0517   ROC->Draw("AL");
0518   ROC->Write();
0519   c.SaveAs("MVA_plots/ROC.pdf");
0520
0521   ROC_mu->Draw("AL");
0522   ROC_mu->Write();
0523   c.SaveAs("MVA_plots/ROC_mu.pdf");
0524
0525   ROC_el->Draw("AL");
0526   ROC_el->Write();
0527   c.SaveAs("MVA_plots/ROC_el.pdf");
0528
0529   ROC_had->Draw("AL");
0530   ROC_had->Write();
0531   c.SaveAs("MVA_plots/ROC_had.pdf");
0532   c.Clear();
0533
0534   TPR_vs_dt->Draw();
0535   FPR_vs_dt->Draw("same");
0536   TPR_vs_dt->SetTitle("Performance vs. decision threshold; decision thresh.; performance measure");
0537   TLegend* leg2 = new TLegend();
0538   leg2->AddEntry(TPR_vs_dt, "TPR", "l");
0539   leg2->AddEntry(FPR_vs_dt, "FPR", "l");
0540   leg2->Draw("same");
0541   c.Write("TPR_FPR_vs_dt");
0542   c.SaveAs("MVA_plots/TPR_FPR_vs_dt.pdf");
0543   c.Clear();
0544
0545   TPR_vs_dt_mu->Draw();
0546   FPR_vs_dt->Draw("same");
0547   TPR_vs_dt_mu->SetTitle("Performance vs. decision threshold (muons); decision thresh.; performance measure");
0548   TLegend* leg3 = new TLegend();
0549   leg3->AddEntry(TPR_vs_dt_mu, "TPR", "l");
0550   leg3->AddEntry(FPR_vs_dt, "FPR", "l");
0551   leg3->Draw("same");
0552   c.Write("TPR_FPR_vs_dt_mu");
0553   c.SaveAs("MVA_plots/TPR_FPR_vs_dt_mu.pdf");
0554   c.Clear();
0555
0556   TPR_vs_dt_el->Draw();
0557   FPR_vs_dt->Draw("same");
0558   TPR_vs_dt_el->SetTitle("Performance vs. decision threshold (electrons); decision thresh.; performance measure");
0559   TLegend* leg4 = new TLegend();
0560   leg4->AddEntry(TPR_vs_dt_el, "TPR", "l");
0561   leg4->AddEntry(FPR_vs_dt, "FPR", "l");
0562   leg4->Draw("same");
0563   c.Write("TPR_FPR_vs_dt_el");
0564   c.SaveAs("MVA_plots/TPR_FPR_vs_dt_el.pdf");
0565   c.Clear();
0566
0567   TPR_vs_dt_had->Draw();
0568   FPR_vs_dt->Draw("same");
0569   TPR_vs_dt_had->SetTitle("Performance vs. decision threshold (hadrons); decision thresh.; performance measure");
0570   TLegend* leg5 = new TLegend();
0571   leg5->AddEntry(TPR_vs_dt_had, "TPR", "l");
0572   leg5->AddEntry(FPR_vs_dt, "FPR", "l");
0573   leg5->Draw("same");
0574   c.Write("TPR_FPR_vs_dt_had");
0575   c.SaveAs("MVA_plots/TPR_FPR_vs_dt_had.pdf");
0576   c.Clear();
0577
0578   TPR_vs_eta->Draw("ap");
0579   TPR_vs_eta->Write();
0580   c.SaveAs("MVA_plots/TPR_vs_eta.pdf");
0581
0582   TPR_vs_eta_mu->Draw("ap");
0583   TPR_vs_eta_mu->Write();
0584   c.SaveAs("MVA_plots/TPR_vs_eta_mu.pdf");
0585
0586   TPR_vs_eta_el->Draw("ap");
0587   TPR_vs_eta_el->Write();
0588   c.SaveAs("MVA_plots/TPR_vs_eta_el.pdf");
0589
0590   TPR_vs_eta_had->Draw("ap");
0591   TPR_vs_eta_had->Write();
0592   c.SaveAs("MVA_plots/TPR_vs_eta_had.pdf");
0593
0594   FPR_vs_eta->Draw("ap");
0595   FPR_vs_eta->Write();
0596   c.SaveAs("MVA_plots/FPR_vs_eta.pdf");
0597
0598   TPR_vs_pt->Draw("ap");
0599   TPR_vs_pt->Write();
0600   c.SaveAs("MVA_plots/TPR_vs_pt.pdf");
0601
0602   TPR_vs_pt_mu->Draw("ap");
0603   TPR_vs_pt_mu->Write();
0604   c.SaveAs("MVA_plots/TPR_vs_pt_mu.pdf");
0605
0606   TPR_vs_pt_el->Draw("ap");
0607   TPR_vs_pt_el->Write();
0608   c.SaveAs("MVA_plots/TPR_vs_pt_el.pdf");
0609
0610   TPR_vs_pt_had->Draw("ap");
0611   TPR_vs_pt_had->Write();
0612   c.SaveAs("MVA_plots/TPR_vs_pt_had.pdf");
0613
0614   FPR_vs_pt->Draw("ap");
0615   FPR_vs_pt->Write();
0616   c.SaveAs("MVA_plots/FPR_vs_pt.pdf");
0617
0618   fout->Close();
0619 }
0620
0621 void SetPlotStyle() {
0622   // from ATLAS plot style macro
0623
0624   // use plain black on white colors
0625   gStyle->SetFrameBorderMode(0);
0626   gStyle->SetFrameFillColor(0);
0627   gStyle->SetCanvasBorderMode(0);
0628   gStyle->SetCanvasColor(0);
0629   gStyle->SetPadBorderMode(0);
0630   gStyle->SetPadColor(0);
0631   gStyle->SetStatColor(0);
0632   gStyle->SetHistLineColor(1);
0633
0634   gStyle->SetPalette(1);
0635
0636   // set the paper & margin sizes
0637   gStyle->SetPaperSize(20, 26);
0638   gStyle->SetPadTopMargin(0.05);
0639   gStyle->SetPadRightMargin(0.05);
0640   gStyle->SetPadBottomMargin(0.16);
0641   gStyle->SetPadLeftMargin(0.16);
0642
0643   // set title offsets (for axis label)
0644   gStyle->SetTitleXOffset(1.4);
0645   gStyle->SetTitleYOffset(1.4);
0646
0647   // use large fonts
0648   gStyle->SetTextFont(42);
0649   gStyle->SetTextSize(0.05);
0650   gStyle->SetLabelFont(42, "x");
0651   gStyle->SetTitleFont(42, "x");
0652   gStyle->SetLabelFont(42, "y");
0653   gStyle->SetTitleFont(42, "y");
0654   gStyle->SetLabelFont(42, "z");
0655   gStyle->SetTitleFont(42, "z");
0656   gStyle->SetLabelSize(0.05, "x");
0657   gStyle->SetTitleSize(0.05, "x");
0658   gStyle->SetLabelSize(0.05, "y");
0659   gStyle->SetTitleSize(0.05, "y");
0660   gStyle->SetLabelSize(0.05, "z");
0661   gStyle->SetTitleSize(0.05, "z");
0662
0663   // use bold lines and markers
0664   //gStyle->SetMarkerStyle(20);
0665   gStyle->SetMarkerSize(1.2);
0666   gStyle->SetHistLineWidth(4.);
0667   gStyle->SetLineStyleString(4, "[12 12]");
0668
0669   // get rid of error bar caps
0670   gStyle->SetEndErrorSize(0.);
0671
0672   // do not display any of the standard histogram decorations
0673   gStyle->SetOptTitle(0);
0674   gStyle->SetOptStat(0);
0675   gStyle->SetOptFit(0);
0676
0677   // put tick marks on top and RHS of plots
0678   gStyle->SetPadTickX(1);
0679   gStyle->SetPadTickY(1);
0680 }