TrackFindingTracklet/test/L1TrackNtuplePlot.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: .L L1TrackNtuplePlot.C++, L1TrackNtuplePlot("TTbar_PU200_hybrid")
0005 //
0006 // By Louise Skinnari, June 2013
0007 // ----------------------------------------------------------------------------------------------------------------
0008
0009 #include "TROOT.h"
0010 #include "TStyle.h"
0011 #include "TLatex.h"
0012 #include "TFile.h"
0013 #include "TTree.h"
0014 #include "TChain.h"
0015 #include "TBranch.h"
0016 #include "TLeaf.h"
0017 #include "TCanvas.h"
0018 #include "TLegend.h"
0019 #include "TH1.h"
0020 #include "TH2.h"
0021 #include "TF1.h"
0022 #include "TProfile.h"
0023 #include "TProfile2D.h"
0024 #include "TMath.h"
0025 #include <TError.h>
0026 #include "TSystem.h"
0027
0028 #include <iostream>
0029 #include <string>
0030 #include <vector>
0031 #include <algorithm>
0032
0033 using namespace std;
0034
0035 void SetPlotStyle();
0036 void mySmallText(Double_t x, Double_t y, Color_t color, char* text);
0037
0038 double getIntervalContainingFractionOfEntries(TH1* histogram, double interval, int minEntries = 5);
0039 void makeResidualIntervalPlot(
0040     TString type, TString dir, TString variable, TH1F* h_68, TH1F* h_90, TH1F* h_99, double minY, double maxY);
0041
0042 // ----------------------------------------------------------------------------------------------------------------
0043 // Main script
0044 // ----------------------------------------------------------------------------------------------------------------
0045
0046 void L1TrackNtuplePlot(TString type,
0047                        TString type_dir = "",
0048                        TString treeName = "",
0049                        int TP_select_injet = 0,
0050                        int TP_select_pdgid = 0,
0051                        int TP_select_eventid = 0,
0052                        bool useTightCuts = false,
0053                        bool useDeadRegion = false,
0054                        float TP_minPt = 2.0,
0055                        float TP_maxPt = 100.0,
0056                        float TP_maxEta = 2.4,
0057                        float TP_maxDxy = 1.0,
0058                        float TP_maxD0 = 1.0,
0059                        bool doDetailedPlots = false) {
0060   // type:              this is the name of the input file you want to process (minus ".root" extension)
0061   // type_dir:          this is the directory containing the input file you want to process. Note that this must end with a "/", as in "EventSets/"
0062   // TP_select_pdgid:   if non-zero, only select TPs with a given PDG ID
0063   // TP_select_eventid: if zero, only look at TPs from primary interaction, else, include TPs from pileup
0064   // TP_minPt:          only look at TPs with pt > X GeV
0065   // TP_maxPt:          only look at TPs with pt < X GeV
0066   // TP_maxEta:         only look at TPs with |eta| < X
0067   // doDetailedPlots:   includes extra plots, such as  performance vs d0.
0068
0069   // TP_select_injet: only look at TPs that are within a jet with pt > 30 GeV (==1) or within a jet with pt > 100 GeV (==2), >200 GeV (==3) or all TPs (==0)
0070
0071   //--  N.B. For standard displaced tracking plots, set TP_minPt=3.0, TP_maxEta=2.0, TP_maxDxy=10.0,
0072   //--  TO_maxD0=10.0, doDetailedPlots=true. (Efficiency plots vs eta also usually made for d0 < 5).
0073
0074   gROOT->SetBatch();
0075   gErrorIgnoreLevel = kWarning;
0076
0077   SetPlotStyle();
0078
0079   // ----------------------------------------------------------------------------------------------------------------
0080   // define input options
0081
0082   // these are the LOOSE cuts, baseline scenario for efficiency and rate plots ==> configure as appropriate
0083   int L1Tk_minNstub = 4;
0084   float L1Tk_maxChi2 = 999999;
0085   float L1Tk_maxChi2dof = 999999.;
0086
0087   // TIGHT cuts (separate plots / rates) ==> configure as appropriate
0088   // this is currently set up as an either or for performance plots, to not duplicate a ton of code.
0089   int L1Tk_TIGHT_minNstub = 4;
0090   float L1Tk_TIGHT_maxChi2 = 999999;
0091   float L1Tk_TIGHT_maxChi2dof = 999999.;
0092   if (useTightCuts) {
0093     L1Tk_minNstub = L1Tk_TIGHT_minNstub;
0094     L1Tk_maxChi2 = L1Tk_TIGHT_maxChi2;
0095     L1Tk_maxChi2dof = L1Tk_TIGHT_maxChi2dof;
0096   }
0097
0098   bool doGausFit = false;     //do gaussian fit for resolution vs eta/pt plots
0099   bool doLooseMatch = false;  //looser MC truth matching
0100
0101   // tracklet variables
0102   int L1Tk_seed = 0;
0103
0104   //some counters for integrated efficiencies
0105   int n_all_eta2p5 = 0;
0106   int n_all_eta1p75 = 0;
0107   int n_all_eta1p0 = 0;
0108   int n_match_eta2p5 = 0;
0109   int n_match_eta1p75 = 0;
0110   int n_match_eta1p0 = 0;
0111   int n_all_ptg2 = 0;
0112   int n_all_ptg8 = 0;
0113   int n_all_pt2to8 = 0;
0114   int n_all_ptg40 = 0;
0115   int n_match_ptg2 = 0;
0116   int n_match_ptg8 = 0;
0117   int n_match_pt2to8 = 0;
0118   int n_match_ptg40 = 0;
0119
0120   // counters for total track rates
0121   int ntrk = 0;
0122   int ntrk_pt2 = 0;
0123   int ntrk_pt3 = 0;
0124   int ntrk_pt10 = 0;
0125   int ntp_pt2 = 0;
0126   int ntp_pt3 = 0;
0127   int ntp_pt10 = 0;
0128
0129   // ----------------------------------------------------------------------------------------------------------------
0130   // read ntuples
0131   TChain* tree = new TChain("L1TrackNtuple" + treeName + "/eventTree");
0132   tree->Add(type_dir + type + ".root");
0133
0134   if (tree->GetEntries() == 0) {
0135     cout << "File doesn't exist or is empty, returning..."
0136          << endl;  //cout's kept in this file as it is an example standalone plotting script, not running in central CMSSW
0137     return;
0138   }
0139
0140   // ----------------------------------------------------------------------------------------------------------------
0141   // define leafs & branches
0142
0143   // tracking particles
0144   vector<float>* tp_pt;
0145   vector<float>* tp_eta;
0146   vector<float>* tp_phi;
0147   vector<float>* tp_dxy;
0148   vector<float>* tp_z0;
0149   vector<float>* tp_d0;
0150   vector<int>* tp_pdgid;
0151   vector<int>* tp_nmatch;
0152   vector<int>* tp_nstub;
0153   vector<int>* tp_eventid;
0154   vector<int>* tp_injet;
0155   vector<int>* tp_injet_highpt;
0156   vector<int>* tp_injet_vhighpt;
0157
0158   // *L1 track* properties, for tracking particles matched to a L1 track
0159   vector<float>* matchtrk_pt;
0160   vector<float>* matchtrk_eta;
0161   vector<float>* matchtrk_phi;
0162   vector<float>* matchtrk_d0;
0163   vector<float>* matchtrk_z0;
0164   vector<float>* matchtrk_chi2;
0165   vector<float>* matchtrk_chi2rphi;
0166   vector<float>* matchtrk_chi2rz;
0167   vector<int>* matchtrk_nstub;
0168   vector<int>* matchtrk_lhits;
0169   vector<int>* matchtrk_dhits;
0170   vector<int>* matchtrk_seed;
0171   vector<int>* matchtrk_hitpattern;
0172   vector<int>* matchtrk_injet;
0173   vector<int>* matchtrk_injet_highpt;
0174   vector<int>* matchtrk_injet_vhighpt;
0175
0176   // all L1 tracks
0177   vector<float>* trk_pt;
0178   vector<float>* trk_eta;
0179   vector<float>* trk_phi;
0180   vector<float>* trk_chi2;
0181   vector<float>* trk_chi2rphi;
0182   vector<float>* trk_chi2rz;
0183   vector<int>* trk_nstub;
0184   vector<int>* trk_lhits;
0185   vector<int>* trk_dhits;
0186   vector<int>* trk_seed;
0187   vector<int>* trk_hitpattern;
0188   vector<unsigned int>* trk_phiSector;
0189   vector<int>* trk_injet;
0190   vector<int>* trk_injet_highpt;
0191   vector<int>* trk_injet_vhighpt;
0192   vector<int>* trk_fake;
0193   vector<int>* trk_genuine;
0194   vector<int>* trk_loose;
0195
0196   TBranch* b_tp_pt;
0197   TBranch* b_tp_eta;
0198   TBranch* b_tp_phi;
0199   TBranch* b_tp_dxy;
0200   TBranch* b_tp_z0;
0201   TBranch* b_tp_d0;
0202   TBranch* b_tp_pdgid;
0203   TBranch* b_tp_nmatch;
0204   TBranch* b_tp_nstub;
0205   TBranch* b_tp_eventid;
0206   TBranch* b_tp_injet;
0207   TBranch* b_tp_injet_highpt;
0208   TBranch* b_tp_injet_vhighpt;
0209
0210   TBranch* b_matchtrk_pt;
0211   TBranch* b_matchtrk_eta;
0212   TBranch* b_matchtrk_phi;
0213   TBranch* b_matchtrk_d0;
0214   TBranch* b_matchtrk_z0;
0215   TBranch* b_matchtrk_chi2;
0216   TBranch* b_matchtrk_chi2rphi;
0217   TBranch* b_matchtrk_chi2rz;
0218   TBranch* b_matchtrk_nstub;
0219   TBranch* b_matchtrk_lhits;
0220   TBranch* b_matchtrk_dhits;
0221   TBranch* b_matchtrk_seed;
0222   TBranch* b_matchtrk_hitpattern;
0223   TBranch* b_matchtrk_injet;
0224   TBranch* b_matchtrk_injet_highpt;
0225   TBranch* b_matchtrk_injet_vhighpt;
0226
0227   TBranch* b_trk_pt;
0228   TBranch* b_trk_eta;
0229   TBranch* b_trk_phi;
0230   TBranch* b_trk_chi2;
0231   TBranch* b_trk_chi2rphi;
0232   TBranch* b_trk_chi2rz;
0233   TBranch* b_trk_nstub;
0234   TBranch* b_trk_lhits;
0235   TBranch* b_trk_dhits;
0236   TBranch* b_trk_phiSector;
0237   TBranch* b_trk_seed;
0238   TBranch* b_trk_hitpattern;
0239   TBranch* b_trk_injet;
0240   TBranch* b_trk_injet_highpt;
0241   TBranch* b_trk_injet_vhighpt;
0242   TBranch* b_trk_fake;
0243   TBranch* b_trk_genuine;
0244   TBranch* b_trk_loose;
0245
0246   tp_pt = 0;
0247   tp_eta = 0;
0248   tp_phi = 0;
0249   tp_dxy = 0;
0250   tp_z0 = 0;
0251   tp_d0 = 0;
0252   tp_pdgid = 0;
0253   tp_nmatch = 0;
0254   tp_nstub = 0;
0255   tp_eventid = 0;
0256   tp_injet = 0;
0257   tp_injet_highpt = 0;
0258   tp_injet_vhighpt = 0;
0259
0260   matchtrk_pt = 0;
0261   matchtrk_eta = 0;
0262   matchtrk_phi = 0;
0263   matchtrk_d0 = 0;
0264   matchtrk_z0 = 0;
0265   matchtrk_chi2 = 0;
0266   matchtrk_chi2rphi = 0;
0267   matchtrk_chi2rz = 0;
0268   matchtrk_nstub = 0;
0269   matchtrk_lhits = 0;
0270   matchtrk_dhits = 0;
0271   matchtrk_seed = 0;
0272   matchtrk_hitpattern = 0;
0273   matchtrk_injet = 0;
0274   matchtrk_injet_highpt = 0;
0275   matchtrk_injet_vhighpt = 0;
0276
0277   trk_pt = 0;
0278   trk_eta = 0;
0279   trk_phi = 0;
0280   trk_chi2 = 0;
0281   trk_chi2rphi = 0;
0282   trk_chi2rz = 0;
0283   trk_nstub = 0;
0284   trk_lhits = 0;
0285   trk_dhits = 0;
0286   trk_phiSector = 0;
0287   trk_seed = 0;
0288   trk_hitpattern = 0;
0289   trk_injet = 0;
0290   trk_injet_highpt = 0;
0291   trk_injet_vhighpt = 0;
0292   trk_fake = 0;
0293   trk_genuine = 0;
0294   trk_loose = 0;
0295
0296   tree->SetBranchAddress("tp_pt", &tp_pt, &b_tp_pt);
0297   tree->SetBranchAddress("tp_eta", &tp_eta, &b_tp_eta);
0298   tree->SetBranchAddress("tp_phi", &tp_phi, &b_tp_phi);
0299   tree->SetBranchAddress("tp_dxy", &tp_dxy, &b_tp_dxy);
0300   tree->SetBranchAddress("tp_z0", &tp_z0, &b_tp_z0);
0301   tree->SetBranchAddress("tp_d0", &tp_d0, &b_tp_d0);
0302   tree->SetBranchAddress("tp_pdgid", &tp_pdgid, &b_tp_pdgid);
0303   if (doLooseMatch)
0304     tree->SetBranchAddress("tp_nloosematch", &tp_nmatch, &b_tp_nmatch);
0305   else
0306     tree->SetBranchAddress("tp_nmatch", &tp_nmatch, &b_tp_nmatch);
0307   tree->SetBranchAddress("tp_nstub", &tp_nstub, &b_tp_nstub);
0308   tree->SetBranchAddress("tp_eventid", &tp_eventid, &b_tp_eventid);
0309   if (TP_select_injet > 0) {
0310     tree->SetBranchAddress("tp_injet", &tp_injet, &b_tp_injet);
0311     tree->SetBranchAddress("tp_injet_highpt", &tp_injet_highpt, &b_tp_injet_highpt);
0312     tree->SetBranchAddress("tp_injet_vhighpt", &tp_injet_vhighpt, &b_tp_injet_vhighpt);
0313   }
0314
0315   if (doLooseMatch) {
0316     tree->SetBranchAddress("loosematchtrk_pt", &matchtrk_pt, &b_matchtrk_pt);
0317     tree->SetBranchAddress("loosematchtrk_eta", &matchtrk_eta, &b_matchtrk_eta);
0318     tree->SetBranchAddress("loosematchtrk_phi", &matchtrk_phi, &b_matchtrk_phi);
0319     tree->SetBranchAddress("loosematchtrk_d0", &matchtrk_d0, &b_matchtrk_d0);
0320     tree->SetBranchAddress("loosematchtrk_z0", &matchtrk_z0, &b_matchtrk_z0);
0321     tree->SetBranchAddress("loosematchtrk_chi2", &matchtrk_chi2, &b_matchtrk_chi2);
0322     tree->SetBranchAddress("loosematchtrk_chi2rphi", &matchtrk_chi2rphi, &b_matchtrk_chi2rphi);
0323     tree->SetBranchAddress("loosematchtrk_chi2rz", &matchtrk_chi2rz, &b_matchtrk_chi2rz);
0324     tree->SetBranchAddress("loosematchtrk_nstub", &matchtrk_nstub, &b_matchtrk_nstub);
0325     tree->SetBranchAddress("loosematchtrk_seed", &matchtrk_seed, &b_matchtrk_seed);
0326     tree->SetBranchAddress("loosematchtrk_hitpattern", &matchtrk_hitpattern, &b_matchtrk_hitpattern);
0327     if (TP_select_injet > 0) {
0328       tree->SetBranchAddress("loosematchtrk_injet", &matchtrk_injet, &b_matchtrk_injet);
0329       tree->SetBranchAddress("loosematchtrk_injet_highpt", &matchtrk_injet_highpt, &b_matchtrk_injet_highpt);
0330       tree->SetBranchAddress("loosematchtrk_injet_vhighpt", &matchtrk_injet_vhighpt, &b_matchtrk_injet_vhighpt);
0331     }
0332   } else {
0333     tree->SetBranchAddress("matchtrk_pt", &matchtrk_pt, &b_matchtrk_pt);
0334     tree->SetBranchAddress("matchtrk_eta", &matchtrk_eta, &b_matchtrk_eta);
0335     tree->SetBranchAddress("matchtrk_phi", &matchtrk_phi, &b_matchtrk_phi);
0336     tree->SetBranchAddress("matchtrk_d0", &matchtrk_d0, &b_matchtrk_d0);
0337     tree->SetBranchAddress("matchtrk_z0", &matchtrk_z0, &b_matchtrk_z0);
0338     tree->SetBranchAddress("matchtrk_chi2", &matchtrk_chi2, &b_matchtrk_chi2);
0339     tree->SetBranchAddress("matchtrk_chi2rphi", &matchtrk_chi2rphi, &b_matchtrk_chi2rphi);
0340     tree->SetBranchAddress("matchtrk_chi2rz", &matchtrk_chi2rz, &b_matchtrk_chi2rz);
0341     tree->SetBranchAddress("matchtrk_nstub", &matchtrk_nstub, &b_matchtrk_nstub);
0342     tree->SetBranchAddress("matchtrk_lhits", &matchtrk_lhits, &b_matchtrk_lhits);
0343     tree->SetBranchAddress("matchtrk_dhits", &matchtrk_dhits, &b_matchtrk_dhits);
0344     tree->SetBranchAddress("matchtrk_seed", &matchtrk_seed, &b_matchtrk_seed);
0345     tree->SetBranchAddress("matchtrk_hitpattern", &matchtrk_hitpattern, &b_matchtrk_hitpattern);
0346     if (TP_select_injet > 0) {
0347       tree->SetBranchAddress("matchtrk_injet", &matchtrk_injet, &b_matchtrk_injet);
0348       tree->SetBranchAddress("matchtrk_injet_highpt", &matchtrk_injet_highpt, &b_matchtrk_injet_highpt);
0349       tree->SetBranchAddress("matchtrk_injet_vhighpt", &matchtrk_injet_vhighpt, &b_matchtrk_injet_vhighpt);
0350     }
0351   }
0352
0353   tree->SetBranchAddress("trk_pt", &trk_pt, &b_trk_pt);
0354   tree->SetBranchAddress("trk_eta", &trk_eta, &b_trk_eta);
0355   tree->SetBranchAddress("trk_phi", &trk_phi, &b_trk_phi);
0356   tree->SetBranchAddress("trk_chi2", &trk_chi2, &b_trk_chi2);
0357   tree->SetBranchAddress("trk_chi2rphi", &trk_chi2rphi, &b_trk_chi2rphi);
0358   tree->SetBranchAddress("trk_chi2rz", &trk_chi2rz, &b_trk_chi2rz);
0359   tree->SetBranchAddress("trk_nstub", &trk_nstub, &b_trk_nstub);
0360   tree->SetBranchAddress("trk_lhits", &trk_lhits, &b_trk_lhits);
0361   tree->SetBranchAddress("trk_dhits", &trk_dhits, &b_trk_dhits);
0362   tree->SetBranchAddress("trk_phiSector", &trk_phiSector, &b_trk_phiSector);
0363   tree->SetBranchAddress("trk_seed", &trk_seed, &b_trk_seed);
0364   tree->SetBranchAddress("trk_hitpattern", &trk_hitpattern, &b_trk_hitpattern);
0365   tree->SetBranchAddress("trk_fake", &trk_fake, &b_trk_fake);
0366   tree->SetBranchAddress("trk_genuine", &trk_genuine, &b_trk_genuine);
0367   tree->SetBranchAddress("trk_loose", &trk_loose, &b_trk_loose);
0368   if (TP_select_injet > 0) {
0369     tree->SetBranchAddress("trk_injet", &trk_injet, &b_trk_injet);
0370     tree->SetBranchAddress("trk_injet_highpt", &trk_injet_highpt, &b_trk_injet_highpt);
0371     tree->SetBranchAddress("trk_injet_vhighpt", &trk_injet_vhighpt, &b_trk_injet_vhighpt);
0372   }
0373
0374   // ----------------------------------------------------------------------------------------------------------------
0375   // histograms
0376   // ----------------------------------------------------------------------------------------------------------------
0377
0378   /////////////////////////////////////////////////
0379   // NOTATION:                                   //
0380   // 'C' - Central eta range, |eta|<0.8          //
0381   // 'I' - Intermediate eta range, 0.8<|eta|<1.6 //
0382   // 'F' - Forward eta range, |eta|>1.6          //
0383   //                                             //
0384   // 'L' - Low pt range,  pt<8 GeV               //
0385   // 'H' - High pt range, pt>8 GeV               //
0386   /////////////////////////////////////////////////
0387
0388   // ----------------------------------------------------------------------------------------------------------------
0389   // for efficiencies
0390
0391   TH1F* h_tp_pt = new TH1F("tp_pt", ";Tracking particle p_{T} [GeV]; Tracking particles / 1.0 GeV", 100, 0, 100.0);
0392   TH1F* h_tp_pt_L = new TH1F("tp_pt_L", ";Tracking particle p_{T} [GeV]; Tracking particles / 0.1 GeV", 80, 0, 8.0);
0393   TH1F* h_tp_pt_LC = new TH1F("tp_pt_LC", ";Tracking particle p_{T} [GeV]; Tracking particles / 0.1 GeV", 80, 0, 8.0);
0394   TH1F* h_tp_pt_H = new TH1F("tp_pt_H", ";Tracking particle p_{T} [GeV]; Tracking particles / 1.0 GeV", 92, 8.0, 100.0);
0395   TH1F* h_tp_eta = new TH1F("tp_eta", ";Tracking particle #eta; Tracking particles / 0.1", 50, -2.5, 2.5);
0396   TH1F* h_tp_eta_L = new TH1F("tp_eta_L", ";Tracking particle #eta; Tracking particles / 0.1", 50, -2.5, 2.5);
0397   TH1F* h_tp_eta_H = new TH1F("tp_eta_H", ";Tracking particle #eta; Tracking particles / 0.1", 50, -2.5, 2.5);
0398   TH1F* h_tp_eta_23 = new TH1F("tp_eta_23", ";Tracking particle #eta; Tracking particles / 0.1", 50, -2.5, 2.5);
0399   TH1F* h_tp_eta_35 = new TH1F("tp_eta_35", ";Tracking particle #eta; Tracking particles / 0.1", 50, -2.5, 2.5);
0400   TH1F* h_tp_eta_5 = new TH1F("tp_eta_5", ";Tracking particle #eta; Tracking particles / 0.1", 50, -2.5, 2.5);
0401
0402   TH1F* h_match_tp_pt =
0403       new TH1F("match_tp_pt", ";Tracking particle p_{T} [GeV]; Tracking particles / 1.0 GeV", 100, 0, 100.0);
0404   TH1F* h_match_tp_pt_L =
0405       new TH1F("match_tp_pt_L", ";Tracking particle p_{T} [GeV]; Tracking particles / 0.1 GeV", 80, 0, 8.0);
0406   TH1F* h_match_tp_pt_LC =
0407       new TH1F("match_tp_pt_LC", ";Tracking particle p_{T} [GeV]; Tracking particles / 0.1 GeV", 80, 0, 8.0);
0408   TH1F* h_match_tp_pt_H =
0409       new TH1F("match_tp_pt_H", ";Tracking particle p_{T} [GeV]; Tracking particles / 0.1 GeV", 92, 8.0, 100.0);
0410   TH1F* h_match_tp_eta = new TH1F("match_tp_eta", ";Tracking particle #eta; Tracking particles / 0.1", 50, -2.5, 2.5);
0411   TH1F* h_match_tp_eta_L =
0412       new TH1F("match_tp_eta_L", ";Tracking particle #eta; Tracking particles / 0.1", 50, -2.5, 2.5);
0413   TH1F* h_match_tp_eta_H =
0414       new TH1F("match_tp_eta_H", ";Tracking particle #eta; Tracking particles / 0.1", 50, -2.5, 2.5);
0415   TH1F* h_match_tp_eta_23 =
0416       new TH1F("match_tp_eta_23", ";Tracking particle #eta; Tracking particles / 0.1", 50, -2.5, 2.5);
0417   TH1F* h_match_tp_eta_35 =
0418       new TH1F("match_tp_eta_35", ";Tracking particle #eta; Tracking particles / 0.1", 50, -2.5, 2.5);
0419   TH1F* h_match_tp_eta_5 =
0420       new TH1F("match_tp_eta_5", ";Tracking particle #eta; Tracking particles / 0.1", 50, -2.5, 2.5);
0421
0422   // ----------------------------------------------------------------------------------------------------------------
0423   // Tracklet propogation efficiencies vs. eta for seeding layers
0424
0425   int trackletEffEtaBins = 24;
0426   double trackletEffMaxEta = 2.4;
0427   int numLayers = 11;
0428   TH2F* h_trk_tracklet_hits = new TH2F("trk_tracklet_hits",
0429                                        ";Track |#eta|; Layer index (0-5 = L1-6, 6-10 = D1-5)",
0430                                        trackletEffEtaBins,
0431                                        0,
0432                                        trackletEffMaxEta,
0433                                        11,
0434                                        0,
0435                                        11);  //used to create below hist
0436   TH2F* h_trk_tracklet_eff = new TH2F("trk_tracklet_eff",
0437                                       ";Track |#eta|; Layer index (0-5 = L1-6, 6-10 = D1-5)",
0438                                       trackletEffEtaBins,
0439                                       0,
0440                                       trackletEffMaxEta,
0441                                       11,
0442                                       0,
0443                                       11);
0444
0445   // ----------------------------------------------------------------------------------------------------------------
0446   // resolution vs. pt histograms
0447
0448   // ----------------------------------------------
0449   // for new versions of resolution vs pt/eta plots
0450   unsigned int nBinsPtRes = 500;
0451   double maxPtRes = 30.;
0452
0453   unsigned int nBinsPtRelRes = 1000;
0454   double maxPtRelRes = 10.;
0455
0456   unsigned int nBinsEtaRes = 500;
0457   double maxEtaRes = 0.1;
0458
0459   unsigned int nBinsPhiRes = 500;
0460   double maxPhiRes = 0.2;
0461
0462   unsigned int nBinsZ0Res = 100;
0463   double maxZ0Res = 4.0;
0464   // ----------------------------------------------
0465
0466   const int nRANGE = 20;
0467   TString ptrange[nRANGE] = {"0-5",   "5-10",  "10-15", "15-20", "20-25", "25-30", "30-35", "35-40", "40-45", "45-50",
0468                              "50-55", "55-60", "60-65", "65-70", "70-75", "75-80", "80-85", "85-90", "90-95", "95-100"};
0469
0470   const float pt_resmin = 1.5;
0471   const int nRANGE_L = 13;
0472   TString ptrange_L[nRANGE_L] = {"1.5-2",
0473                                  "2-2.5",
0474                                  "2.5-3",
0475                                  "3-3.5",
0476                                  "3.5-4",
0477                                  "4-4.5",
0478                                  "4.5-5",
0479                                  "5-5.5",
0480                                  "5.5-6",
0481                                  "6-6.5",
0482                                  "6.5-7",
0483                                  "7-7.5",
0484                                  "7.5-8"};
0485
0486   TH1F* h_absResVsPt_pt[nRANGE];
0487   TH1F* h_absResVsPt_ptRel[nRANGE];
0488   TH1F* h_absResVsPt_z0[nRANGE];
0489   TH1F* h_absResVsPt_phi[nRANGE];
0490   TH1F* h_absResVsPt_eta[nRANGE];
0491   TH1F* h_absResVsPt_d0[nRANGE];
0492
0493   TH1F* h_absResVsPt_pt_L[nRANGE_L];
0494   TH1F* h_absResVsPt_ptRel_L[nRANGE_L];
0495   TH1F* h_absResVsPt_z0_L[nRANGE_L];
0496   TH1F* h_absResVsPt_phi_L[nRANGE_L];
0497   TH1F* h_absResVsPt_eta_L[nRANGE_L];
0498   TH1F* h_absResVsPt_d0_L[nRANGE_L];
0499
0500   for (int i = 0; i < nRANGE; i++) {
0501     h_absResVsPt_pt[i] = new TH1F(
0502         "absResVsPt_pt_" + ptrange[i], ";p_{T} residual (L1 - sim) [GeV]; L1 tracks / 0.1", nBinsPtRes, 0, maxPtRes);
0503     h_absResVsPt_ptRel[i] = new TH1F("absResVsPt_ptRel_" + ptrange[i],
0504                                      ";p_{T} residual (L1 - sim) / p_{T}; L1 tracks / 0.02",
0505                                      nBinsPtRelRes,
0506                                      0,
0507                                      maxPtRelRes);
0508     h_absResVsPt_z0[i] = new TH1F(
0509         "absResVsPt_z0_" + ptrange[i], ";z_{0} residual (L1 - sim) [GeV]; L1 tracks / 0.1", nBinsZ0Res, 0, maxZ0Res);
0510     h_absResVsPt_phi[i] = new TH1F(
0511         "absResVsPt_phi_" + ptrange[i], ";#phi residual (L1 - sim) [GeV]; L1 tracks / 0.1", nBinsPhiRes, 0, maxPhiRes);
0512     h_absResVsPt_eta[i] = new TH1F(
0513         "absResVsPt_eta_" + ptrange[i], ";#eta residual (L1 - sim) [GeV]; L1 tracks / 0.1", nBinsEtaRes, 0, maxEtaRes);
0514     h_absResVsPt_d0[i] =
0515         new TH1F("absResVsPt_d0_" + ptrange[i], ";d_{0}residual (L1 - sim) [GeV]; L1 tracks / 0.1", 100, 0, 0.02);
0516   }
0517
0518   for (int i = 0; i < nRANGE_L; i++) {
0519     h_absResVsPt_pt_L[i] = new TH1F(
0520         "absResVsPt_L_pt_" + ptrange_L[i], ";p_{T} residual (L1 - sim) [GeV]; L1 tracks / 0.1", nBinsPtRes, 0, maxPtRes);
0521     h_absResVsPt_ptRel_L[i] = new TH1F("absResVsPt_L_ptRel_" + ptrange_L[i],
0522                                        ";p_{T} residual (L1 - sim) / p_{T}; L1 tracks / 0.02",
0523                                        nBinsPtRelRes,
0524                                        0,
0525                                        maxPtRelRes);
0526     h_absResVsPt_z0_L[i] = new TH1F(
0527         "absResVsPt_L_z0_" + ptrange_L[i], ";z_{0} residual (L1 - sim) [GeV]; L1 tracks / 0.1", nBinsZ0Res, 0, maxZ0Res);
0528     h_absResVsPt_phi_L[i] = new TH1F("absResVsPt_L_phi_" + ptrange_L[i],
0529                                      ";#phi residual (L1 - sim) [GeV]; L1 tracks / 0.1",
0530                                      nBinsPhiRes,
0531                                      0,
0532                                      maxPhiRes);
0533     h_absResVsPt_eta_L[i] = new TH1F("absResVsPt_L_eta_" + ptrange_L[i],
0534                                      ";#eta residual (L1 - sim) [GeV]; L1 tracks / 0.1",
0535                                      nBinsEtaRes,
0536                                      0,
0537                                      maxEtaRes);
0538     h_absResVsPt_d0_L[i] =
0539         new TH1F("absResVsPt_L_d0_" + ptrange_L[i], ";d_{0}residual (L1 - sim) [GeV]; L1 tracks / 0.1", 100, 0, 0.02);
0540   }
0541
0542   // resolution vs. eta histograms
0543
0544   const float eta_resmax = 2.5;
0545   const int nETARANGE = 25;
0546   TString etarange[nETARANGE] = {"0.1", "0.2", "0.3", "0.4", "0.5", "0.6", "0.7", "0.8", "0.9",
0547                                  "1.0", "1.1", "1.2", "1.3", "1.4", "1.5", "1.6", "1.7", "1.8",
0548                                  "1.9", "2.0", "2.1", "2.2", "2.3", "2.4", "2.5"};
0549   /*
0550   const int nETARANGE = 12;
0551   TString etarange[nETARANGE] = {"0.2","0.4","0.6","0.8","1.0",
0552                  "1.2","1.4","1.6","1.8","2.0",
0553                  "2.2","2.4"};
0554   */
0555
0556   TH1F* h_absResVsEta_eta[nETARANGE];
0557   TH1F* h_absResVsEta_z0[nETARANGE];
0558   TH1F* h_absResVsEta_phi[nETARANGE];
0559   TH1F* h_absResVsEta_ptRel[nETARANGE];
0560   TH1F* h_absResVsEta_d0[nETARANGE];
0561
0562   TH1F* h_absResVsEta_eta_L[nETARANGE];
0563   TH1F* h_absResVsEta_z0_L[nETARANGE];
0564   TH1F* h_absResVsEta_phi_L[nETARANGE];
0565   TH1F* h_absResVsEta_ptRel_L[nETARANGE];
0566   TH1F* h_absResVsEta_d0_L[nETARANGE];
0567
0568   TH1F* h_absResVsEta_eta_H[nETARANGE];
0569   TH1F* h_absResVsEta_z0_H[nETARANGE];
0570   TH1F* h_absResVsEta_phi_H[nETARANGE];
0571   TH1F* h_absResVsEta_ptRel_H[nETARANGE];
0572   TH1F* h_absResVsEta_d0_H[nETARANGE];
0573
0574   for (int i = 0; i < nETARANGE; i++) {
0575     h_absResVsEta_eta[i] = new TH1F(
0576         "absResVsEta_eta_" + etarange[i], ";#eta residual (L1 - sim) [GeV]; L1 tracks / 0.1", nBinsEtaRes, 0, maxEtaRes);
0577     h_absResVsEta_z0[i] = new TH1F(
0578         "absResVsEta_z0_" + etarange[i], ";|z_{0} residual (L1 - sim)| [cm]; L1 tracks / 0.01", nBinsZ0Res, 0, maxZ0Res);
0579     h_absResVsEta_phi[i] = new TH1F(
0580         "absResVsEta_phi_" + etarange[i], ";#phi residual (L1 - sim) [GeV]; L1 tracks / 0.1", nBinsPhiRes, 0, maxPhiRes);
0581     h_absResVsEta_ptRel[i] = new TH1F("absResVsEta_ptRel_" + etarange[i],
0582                                       ";p_{T} residual (L1 - sim) / p_{T}; L1 tracks / 0.02",
0583                                       nBinsPtRelRes,
0584                                       0,
0585                                       maxPtRelRes);
0586     h_absResVsEta_d0[i] =
0587         new TH1F("absResVsEta_d0_" + etarange[i], ";d_{0}residual (L1 - sim) [GeV]; L1 tracks / 0.1", 100, 0, 0.02);
0588
0589     h_absResVsEta_eta_L[i] = new TH1F("absResVsEta_eta_L_" + etarange[i],
0590                                       ";#eta residual (L1 - sim) [GeV]; L1 tracks / 0.1",
0591                                       nBinsEtaRes,
0592                                       0,
0593                                       maxEtaRes);
0594     h_absResVsEta_z0_L[i] = new TH1F("absResVsEta_z0_L_" + etarange[i],
0595                                      ";|z_{0} residual (L1 - sim)| [cm]; L1 tracks / 0.01",
0596                                      nBinsZ0Res,
0597                                      0,
0598                                      maxZ0Res);
0599     h_absResVsEta_phi_L[i] = new TH1F("absResVsEta_phi_L_" + etarange[i],
0600                                       ";#phi residual (L1 - sim) [GeV]; L1 tracks / 0.1",
0601                                       nBinsPhiRes,
0602                                       0,
0603                                       maxPhiRes);
0604     h_absResVsEta_ptRel_L[i] = new TH1F("absResVsEta_ptRel_L_" + etarange[i],
0605                                         ";p_{T} residual (L1 - sim) / p_{T}; L1 tracks / 0.02",
0606                                         nBinsPtRelRes,
0607                                         0,
0608                                         maxPtRelRes);
0609     h_absResVsEta_d0_L[i] =
0610         new TH1F("absResVsEta_d0_L_" + etarange[i], ";d_{0}residual (L1 - sim) [GeV]; L1 tracks / 0.1", 100, 0, 0.02);
0611
0612     h_absResVsEta_eta_H[i] = new TH1F("absResVsEta_eta_H_" + etarange[i],
0613                                       ";#eta residual (L1 - sim) [GeV]; L1 tracks / 0.1",
0614                                       nBinsEtaRes,
0615                                       0,
0616                                       maxEtaRes);
0617     h_absResVsEta_z0_H[i] = new TH1F("absResVsEta_z0_H_" + etarange[i],
0618                                      ";|z_{0} residual (L1 - sim)| [cm]; L1 tracks / 0.01",
0619                                      nBinsZ0Res,
0620                                      0,
0621                                      maxZ0Res);
0622     h_absResVsEta_phi_H[i] = new TH1F("absResVsEta_phi_H_" + etarange[i],
0623                                       ";#phi residual (L1 - sim) [GeV]; L1 tracks / 0.1",
0624                                       nBinsPhiRes,
0625                                       0,
0626                                       maxPhiRes);
0627     h_absResVsEta_ptRel_H[i] = new TH1F("absResVsEta_ptRel_H_" + etarange[i],
0628                                         ";p_{T} residual (L1 - sim) / p_{T}; L1 tracks / 0.02",
0629                                         nBinsPtRelRes,
0630                                         0,
0631                                         maxPtRelRes);
0632     h_absResVsEta_d0_H[i] =
0633         new TH1F("absResVsEta_d0_H_" + etarange[i], ";d_{0}residual (L1 - sim) [GeV]; L1 tracks / 0.1", 100, 0, 0.02);
0634   }
0635
0636   // resolution vs phi
0637
0638   const int nPHIRANGE = 32;
0639   TString phirange[nPHIRANGE] = {"-3.0", "-2.8", "-2.6", "-2.4", "-2.2", "-2.0", "-1.8", "-1.6", "-1.4", "-1.2", "-1.0",
0640                                  "-0.8", "-0.6", "-0.4", "-0.2", "0.0",  "0.2",  "0.4",  "0.6",  "0.8",  "1.0",  "1.2",
0641                                  "1.4",  "1.6",  "1.8",  "2.0",  "2.2",  "2.4",  "2.6",  "2.8",  "3.0",  "3.2"};
0642
0643   TH1F* h_absResVsPhi_pt[nPHIRANGE];
0644   TH1F* h_absResVsPhi_ptRel[nPHIRANGE];
0645
0646   for (int i = 0; i < nPHIRANGE; i++) {
0647     h_absResVsPhi_pt[i] = new TH1F(
0648         "absResVsPt_pt_" + phirange[i], ";p_{T} residual (L1 - sim) [GeV]; L1 tracks / 0.1", nBinsPtRes, 0, maxPtRes);
0649     h_absResVsPhi_ptRel[i] = new TH1F("absResVsPt_ptRel_" + phirange[i],
0650                                       ";p_{T} residual (L1 - sim) / p_{T}; L1 tracks / 0.02",
0651                                       nBinsPtRelRes,
0652                                       0,
0653                                       maxPtRelRes);
0654   }
0655
0656   // ----------------------------------------------------------------------------------------------------------------
0657   // chi2 histograms (last bin is an overflow bin)
0658
0659   TH1F* h_trk_chi2 = new TH1F("trk_chi2", ";#chi^{2}; L1 tracks / 1.0", 100, 0, 100);
0660   TH1F* h_trk_chi2_dof = new TH1F("trk_chi2_dof", ";#chi^{2} / D.O.F.; L1 tracks / 0.2", 100, 0, 20);
0661   TH1F* h_match_trk_chi2 = new TH1F("match_trk_chi2", ";#chi^{2}; L1 tracks / 1.0", 100, 0, 100);
0662   TH1F* h_match_trk_chi2_dof = new TH1F("match_trk_chi2_dof", ";#chi^{2} / D.O.F.; L1 tracks / 0.2", 100, 0, 20);
0663
0664   TH1F* h_trk_chi2rphi = new TH1F("trk_chi2rphi", ";#chi^{2}_{r-#phi}; L1 tracks / 1.0", 100, 0, 100);
0665   TH1F* h_trk_chi2rphi_dof = new TH1F("trk_chi2rphi_dof", ";#chi^{2}_{r-#phi} / D.O.F.; L1 tracks / 0.2", 100, 0, 20);
0666   TH1F* h_match_trk_chi2rphi = new TH1F("match_trk_chi2rphi", ";#chi^{2}_{r-#phi}; L1 tracks / 1.0", 100, 0, 100);
0667   TH1F* h_match_trk_chi2rphi_dof =
0668       new TH1F("match_trk_chi2rphi_dof", ";#chi^{2}_{r-#phi} / D.O.F.; L1 tracks / 0.2", 100, 0, 20);
0669
0670   TH1F* h_trk_chi2rz = new TH1F("trk_chi2rz", ";#chi^{2}_{r-z}; L1 tracks / 1.0", 100, 0, 100);
0671   TH1F* h_trk_chi2rz_dof = new TH1F("trk_chi2rz_dof", ";#chi^{2}_{r-z} / D.O.F.; L1 tracks / 0.2", 100, 0, 20);
0672   TH1F* h_match_trk_chi2rz = new TH1F("match_trk_chi2rz", ";#chi^{2}_{r-z}; L1 tracks / 1.0", 100, 0, 100);
0673   TH1F* h_match_trk_chi2rz_dof =
0674       new TH1F("match_trk_chi2rz_dof", ";#chi^{2}_{r-z} / D.O.F.; L1 tracks / 0.2", 100, 0, 20);
0675
0676   // ----------------------------------------------------------------------------------------------------------------
0677   // total track rates
0678
0679   TH1F* h_trk_all_vspt = new TH1F("trk_all_vspt", ";Track p_{T} [GeV]; ", 50, 0, 25);
0680   TH1F* h_trk_loose_vspt = new TH1F("trk_loose_vspt", ";Track p_{T} [GeV]; ", 50, 0, 25);
0681   TH1F* h_trk_genuine_vspt = new TH1F("trk_genuine_vspt", ";Track p_{T} [GeV]; ", 50, 0, 25);
0682   TH1F* h_trk_notloose_vspt = new TH1F(
0683       "trk_notloose_vspt", ";Track p_{T} [GeV]; ", 50, 0, 25);  //(same as "fake" according to the trk_fake labeling)
0684   TH1F* h_trk_notgenuine_vspt = new TH1F("trk_notgenuine_vspt", ";Track p_{T} [GeV]; ", 50, 0, 25);
0685   TH1F* h_trk_duplicate_vspt = new TH1F("trk_duplicate_vspt",
0686                                         ";Track p_{T} [GeV]; ",
0687                                         50,
0688                                         0,
0689                                         25);  //where a TP is genuinely matched to more than one L1 track
0690   TH1F* h_tp_vspt = new TH1F("tp_vspt", ";TP p_{T} [GeV]; ", 50, 0, 25);
0691
0692   // ----------------------------------------------------------------------------------------------------------------
0693
0694   TH1F* h_tp_z0 = new TH1F("tp_z0", ";Tracking particle z_{0} [cm]; Tracking particles / 1.0 cm", 50, -25.0, 25.0);
0695   TH1F* h_tp_z0_L = new TH1F("tp_z0_L", ";Tracking particle z_{0} [cm]; Tracking particles / 1.0 cm", 50, -25.0, 25.0);
0696   TH1F* h_tp_z0_H = new TH1F("tp_z0_H", ";Tracking particle z_{0} [cm]; Tracking particles / 1.0 cm", 50, -25.0, 25.0);
0697
0698   TH1F* h_match_tp_z0 =
0699       new TH1F("match_tp_z0", ";Tracking particle z_{0} [cm]; Tracking particles / 1.0 cm", 50, -25.0, 25.0);
0700   TH1F* h_match_tp_z0_L =
0701       new TH1F("match_tp_z0_L", ";Tracking particle z_{0} [cm]; Tracking particles / 1.0 cm", 50, -25.0, 25.0);
0702   TH1F* h_match_tp_z0_H =
0703       new TH1F("match_tp_z0_H", ";Tracking particle z_{0} [cm]; Tracking particles / 1.0 cm", 50, -25.0, 25.0);
0704
0705   // ----------------------------------------------------------------------------------------------------------------
0706   //
0707   //       ******************   additional histograms drawn if using 'detailed' option   ******************
0708   //
0709   // ----------------------------------------------------------------------------------------------------------------
0710
0711   const float maxD0plot = TP_maxD0;
0712
0713   TH1F* h_tp_phi = new TH1F("tp_phi", ";Tracking particle #phi [rad]; Tracking particles / 0.1", 64, -3.2, 3.2);
0714   TH1F* h_tp_d0 =
0715       new TH1F("tp_d0", ";Tracking particle d_{0} [cm]; Tracking particles / 0.01 cm", 50, -maxD0plot, maxD0plot);
0716   TH1F* h_tp_absd0 =
0717       new TH1F("tp_absd0", ";Tracking particle |d_{0}| [cm]; Tracking particles / 0.04 cm", 50, 0, maxD0plot);
0718   TH1F* h_tp_absd0_eta2 =
0719       new TH1F("tp_absd0_eta2", ";Tracking particle |d_{0}| [cm]; Tracking particles / 0.04 cm", 50, 0, maxD0plot);
0720   TH1F* h_tp_absd0_eta2_pt3 =
0721       new TH1F("tp_absd0_eta2_pt3", ";Tracking particle |d_{0}| [cm]; Tracking particles / 0.04 cm", 50, 0, maxD0plot);
0722
0723   TH1F* h_match_tp_phi =
0724       new TH1F("match_tp_phi", ";Tracking particle #phi [rad]; Tracking particles / 0.1", 64, -3.2, 3.2);
0725   TH1F* h_match_tp_d0 =
0726       new TH1F("match_tp_d0", ";Tracking particle d_{0} [cm]; Tracking particles / 0.01 cm", 50, -maxD0plot, maxD0plot);
0727   TH1F* h_match_tp_absd0 =
0728       new TH1F("match_tp_absd0", ";Tracking particle d_{0} [cm]; Tracking particles / 0.04 cm", 50, 0, maxD0plot);
0729   TH1F* h_match_tp_absd0_eta2 =
0730       new TH1F("match_tp_absd0_eta2", ";Tracking particle d_{0} [cm]; Tracking particles / 0.04 cm", 50, 0, maxD0plot);
0731   TH1F* h_match_tp_absd0_eta2_pt3 = new TH1F(
0732       "match_tp_absd0_eta2_pt3", ";Tracking particle d_{0} [cm]; Tracking particles / 0.04 cm", 50, 0, maxD0plot);
0733
0734   TH1F* h_match_trk_nstub = new TH1F("match_trk_nstub", ";Number of stubs; L1 tracks / 1.0", 15, 0, 15);
0735   TH1F* h_match_trk_nstub_C = new TH1F("match_trk_nstub_C", ";Number of stubs; L1 tracks / 1.0", 15, 0, 15);
0736   TH1F* h_match_trk_nstub_I = new TH1F("match_trk_nstub_I", ";Number of stubs; L1 tracks / 1.0", 15, 0, 15);
0737   TH1F* h_match_trk_nstub_F = new TH1F("match_trk_nstub_F", ";Number of stubs; L1 tracks / 1.0", 15, 0, 15);
0738
0739   // note that we are only making the below chi2 histograms using the joint chi2, not the separate chi2rphi and chi2rz
0740
0741   // chi2 histograms
0742   // note: last bin is an overflow bin
0743   TH1F* h_match_trk_chi2_C_L = new TH1F("match_trk_chi2_C_L", ";#chi^{2}; L1 tracks / 1.0", 100, 0, 100);
0744   TH1F* h_match_trk_chi2_I_L = new TH1F("match_trk_chi2_I_L", ";#chi^{2}; L1 tracks / 1.0", 100, 0, 100);
0745   TH1F* h_match_trk_chi2_F_L = new TH1F("match_trk_chi2_F_L", ";#chi^{2}; L1 tracks / 1.0", 100, 0, 100);
0746   TH1F* h_match_trk_chi2_C_H = new TH1F("match_trk_chi2_C_H", ";#chi^{2}; L1 tracks / 1.0", 100, 0, 100);
0747   TH1F* h_match_trk_chi2_I_H = new TH1F("match_trk_chi2_I_H", ";#chi^{2}; L1 tracks / 1.0", 100, 0, 100);
0748   TH1F* h_match_trk_chi2_F_H = new TH1F("match_trk_chi2_F_H", ";#chi^{2}; L1 tracks / 1.0", 100, 0, 100);
0749
0750   // chi2/dof histograms
0751   // note: lastbin is an overflow bin
0752   TH1F* h_match_trk_chi2_dof_C_L =
0753       new TH1F("match_trk_chi2_dof_C_L", ";#chi^{2} / D.O.F.; L1 tracks / 0.2", 100, 0, 20);
0754   TH1F* h_match_trk_chi2_dof_I_L =
0755       new TH1F("match_trk_chi2_dof_I_L", ";#chi^{2} / D.O.F.; L1 tracks / 0.2", 100, 0, 20);
0756   TH1F* h_match_trk_chi2_dof_F_L =
0757       new TH1F("match_trk_chi2_dof_F_L", ";#chi^{2} / D.O.F.; L1 tracks / 0.2", 100, 0, 20);
0758   TH1F* h_match_trk_chi2_dof_C_H =
0759       new TH1F("match_trk_chi2_dof_C_H", ";#chi^{2} / D.O.F.; L1 tracks / 0.2", 100, 0, 20);
0760   TH1F* h_match_trk_chi2_dof_I_H =
0761       new TH1F("match_trk_chi2_dof_I_H", ";#chi^{2} / D.O.F.; L1 tracks / 0.2", 100, 0, 20);
0762   TH1F* h_match_trk_chi2_dof_F_H =
0763       new TH1F("match_trk_chi2_dof_F_H", ";#chi^{2} / D.O.F.; L1 tracks / 0.2", 100, 0, 20);
0764
0765   // ----------------------------------------------------------------------------------------------------------------
0766   // resolution histograms
0767   TH1F* h_res_pt = new TH1F("res_pt", ";p_{T} residual (L1 - sim) [GeV]; L1 tracks / 0.05", 200, -5.0, 5.0);
0768   TH1F* h_res_ptRel = new TH1F("res_ptRel", ";p_{T} residual (L1 - sim) / p_{T}; L1 tracks / 0.01", 200, -1.0, 1.0);
0769   TH1F* h_res_eta = new TH1F("res_eta", ";#eta residual (L1 - sim); L1 tracks / 0.0002", 100, -0.01, 0.01);
0770   TH1F* h_res_phi = new TH1F("res_phi", ";#phi residual (L1 - sim) [rad]; L1 tracks / 0.0001", 100, -0.005, 0.005);
0771
0772   TH1F* h_res_z0 = new TH1F("res_z0", ";z_{0} residual (L1 - sim) [cm]; L1 tracks / 0.02", 100, -1.0, 1.0);
0773   TH1F* h_res_z0_C = new TH1F("res_z0_C", ";z_{0} residual (L1 - sim) [cm]; L1 tracks / 0.02", 100, -1.0, 1.0);
0774   TH1F* h_res_z0_I = new TH1F("res_z0_I", ";z_{0} residual (L1 - sim) [cm]; L1 tracks / 0.02", 100, -1.0, 1.0);
0775   TH1F* h_res_z0_F = new TH1F("res_z0_F", ";z_{0} residual (L1 - sim) [cm]; L1 tracks / 0.02", 100, -1.0, 1.0);
0776   TH1F* h_res_z0_L = new TH1F("res_z0_L", ";z_{0} residual (L1 - sim) [cm]; L1 tracks / 0.02", 100, -1.0, 1.0);
0777   TH1F* h_res_z0_H = new TH1F("res_z0_H", ";z_{0} residual (L1 - sim) [cm]; L1 tracks / 0.02", 100, -1.0, 1.0);
0778
0779   TH1F* h_res_z0_C_L =
0780       new TH1F("res_z0_C_L", ";z_{0} residual (L1 - sim) [cm]; L1 tracks / 0.02", 100, (-1) * 1.0, 1.0);
0781   TH1F* h_res_z0_I_L =
0782       new TH1F("res_z0_I_L", ";z_{0} residual (L1 - sim) [cm]; L1 tracks / 0.02", 100, (-1) * 1.0, 1.0);
0783   TH1F* h_res_z0_F_L =
0784       new TH1F("res_z0_F_L", ";z_{0} residual (L1 - sim) [cm]; L1 tracks / 0.02", 100, (-1) * 1.0, 1.0);
0785   TH1F* h_res_z0_C_H =
0786       new TH1F("res_z0_C_H", ";z_{0} residual (L1 - sim) [cm]; L1 tracks / 0.02", 100, (-1) * 1.0, 1.0);
0787   TH1F* h_res_z0_I_H =
0788       new TH1F("res_z0_I_H", ";z_{0} residual (L1 - sim) [cm]; L1 tracks / 0.02", 100, (-1) * 1.0, 1.0);
0789   TH1F* h_res_z0_F_H =
0790       new TH1F("res_z0_F_H", ";z_{0} residual (L1 - sim) [cm]; L1 tracks / 0.02", 100, (-1) * 1.0, 1.0);
0791
0792   TH1F* h_res_d0 = new TH1F("res_d0", ";d_{0} residual (L1 - sim) [cm]; L1 tracks / 0.0002 cm", 200, -0.02, 0.02);
0793   TH1F* h_res_d0_C = new TH1F("res_d0_C", ";d_{0} residual (L1 - sim) [cm]; L1 tracks / 0.0001 cm", 200, -0.05, 0.05);
0794   TH1F* h_res_d0_I = new TH1F("res_d0_I", ";d_{0} residual (L1 - sim) [cm]; L1 tracks / 0.0001 cm", 200, -0.05, 0.05);
0795   TH1F* h_res_d0_F = new TH1F("res_d0_F", ";d_{0} residual (L1 - sim) [cm]; L1 tracks / 0.0001 cm", 200, -0.05, 0.05);
0796   TH1F* h_res_d0_L = new TH1F("res_d0_L", ";d_{0} residual (L1 - sim) [cm]; L1 tracks / 0.0001 cm", 200, -0.05, 0.05);
0797   TH1F* h_res_d0_H = new TH1F("res_d0_H", ";d_{0} residual (L1 - sim) [cm]; L1 tracks / 0.0001 cm", 200, -0.05, 0.05);
0798
0799   TH1F* h_res_d0_C_L =
0800       new TH1F("res_d0_C_L", ";d_{0} residual (L1 - sim) [cm]; L1 tracks / 0.0001 cm", 200, -0.05, 0.05);
0801   TH1F* h_res_d0_I_L =
0802       new TH1F("res_d0_I_L", ";d_{0} residual (L1 - sim) [cm]; L1 tracks / 0.0001 cm", 200, -0.05, 0.05);
0803   TH1F* h_res_d0_F_L =
0804       new TH1F("res_d0_F_L", ";d_{0} residual (L1 - sim) [cm]; L1 tracks / 0.0001 cm", 200, -0.05, 0.05);
0805   TH1F* h_res_d0_C_H =
0806       new TH1F("res_d0_C_H", ";d_{0} residual (L1 - sim) [cm]; L1 tracks / 0.0001 cm", 200, -0.05, 0.05);
0807   TH1F* h_res_d0_I_H =
0808       new TH1F("res_d0_I_H", ";d_{0} residual (L1 - sim) [cm]; L1 tracks / 0.0001 cm", 200, -0.05, 0.05);
0809   TH1F* h_res_d0_F_H =
0810       new TH1F("res_d0_F_H", ";d_{0} residual (L1 - sim) [cm]; L1 tracks / 0.0001 cm", 200, -0.05, 0.05);
0811
0812   // ----------------------------------------------------------------------------------------------------------------
0813   // more resolution vs pt
0814
0815   TH1F* h_resVsPt_pt[nRANGE];
0816   TH1F* h_resVsPt_pt_C[nRANGE];
0817   TH1F* h_resVsPt_pt_I[nRANGE];
0818   TH1F* h_resVsPt_pt_F[nRANGE];
0819
0820   TH1F* h_resVsPt_ptRel[nRANGE];
0821   TH1F* h_resVsPt_ptRel_C[nRANGE];
0822   TH1F* h_resVsPt_ptRel_I[nRANGE];
0823   TH1F* h_resVsPt_ptRel_F[nRANGE];
0824
0825   TH1F* h_resVsPt_z0[nRANGE];
0826   TH1F* h_resVsPt_z0_C[nRANGE];
0827   TH1F* h_resVsPt_z0_I[nRANGE];
0828   TH1F* h_resVsPt_z0_F[nRANGE];
0829
0830   TH1F* h_resVsPt_phi[nRANGE];
0831   TH1F* h_resVsPt_phi_C[nRANGE];
0832   TH1F* h_resVsPt_phi_I[nRANGE];
0833   TH1F* h_resVsPt_phi_F[nRANGE];
0834
0835   TH1F* h_resVsPt_eta[nRANGE];
0836   TH1F* h_resVsPt_d0[nRANGE];
0837
0838   for (int i = 0; i < nRANGE; i++) {
0839     h_resVsPt_pt[i] =
0840         new TH1F("resVsPt_pt_" + ptrange[i], ";p_{T} residual (L1 - sim) [GeV]; L1 tracks / 0.1", 100, -5.0, 5.0);
0841     h_resVsPt_pt_C[i] =
0842         new TH1F("resVsPt_pt_C_" + ptrange[i], ";p_{T} residual (L1 - sim) [GeV]; L1 tracks / 0.1", 100, -5.0, 5.0);
0843     h_resVsPt_pt_I[i] =
0844         new TH1F("resVsPt_pt_I_" + ptrange[i], ";p_{T} residual (L1 - sim) [GeV]; L1 tracks / 0.1", 100, -5.0, 5.0);
0845     h_resVsPt_pt_F[i] =
0846         new TH1F("resVsPt_pt_F_" + ptrange[i], ";p_{T} residual (L1 - sim) [GeV]; L1 tracks / 0.1", 100, -5.0, 5.0);
0847
0848     h_resVsPt_ptRel[i] = new TH1F(
0849         "resVsPt_ptRel_" + ptrange[i], ";p_{T} residual (L1 - sim) / p_{T}; L1 tracks / 0.02", 300, -0.15, 0.15);
0850     h_resVsPt_ptRel_C[i] = new TH1F(
0851         "resVsPt_ptRel_c_" + ptrange[i], ";p_{T} residual (L1 - sim) / p_{T}; L1 tracks / 0.02", 300, -0.15, 0.15);
0852     h_resVsPt_ptRel_I[i] = new TH1F(
0853         "resVsPt_ptRel_I_" + ptrange[i], ";p_{T} residual (L1 - sim) / p_{T}; L1 tracks / 0.02", 300, -0.15, 0.15);
0854     h_resVsPt_ptRel_F[i] = new TH1F(
0855         "resVsPt_ptRel_F_" + ptrange[i], ";p_{T} residual (L1 - sim) / p_{T}; L1 tracks / 0.02", 300, -0.15, 0.15);
0856
0857     h_resVsPt_z0[i] =
0858         new TH1F("resVsPt_z0_" + ptrange[i], ";z_{0} residual (L1 - sim) [cm]; L1 tracks / 0.02", 100, -1, 1);
0859     h_resVsPt_z0_C[i] =
0860         new TH1F("resVsPt_z0_C_" + ptrange[i], ";z_{0} residual (L1 - sim) [cm]; L1 tracks / 0.02", 100, -1, 1);
0861     h_resVsPt_z0_I[i] =
0862         new TH1F("resVsPt_z0_I_" + ptrange[i], ";z_{0} residual (L1 - sim) [cm]; L1 tracks / 0.02", 100, -1, 1);
0863     h_resVsPt_z0_F[i] =
0864         new TH1F("resVsPt_z0_F_" + ptrange[i], ";z_{0} residual (L1 - sim) [cm]; L1 tracks / 0.02", 100, -1, 1);
0865
0866     h_resVsPt_phi[i] = new TH1F(
0867         "resVsPt_phi_" + ptrange[i], ";#phi residual (L1 - sim) [rad]; L1 tracks / 0.0001", 100, -0.005, 0.005);
0868     h_resVsPt_phi_C[i] = new TH1F(
0869         "resVsPt_phi_C_" + ptrange[i], ";#phi residual (L1 - sim) [rad]; L1 tracks / 0.0001", 100, -0.005, 0.005);
0870     h_resVsPt_phi_I[i] = new TH1F(
0871         "resVsPt_phi_I_" + ptrange[i], ";#phi residual (L1 - sim) [rad]; L1 tracks / 0.0001", 100, -0.005, 0.005);
0872     h_resVsPt_phi_F[i] = new TH1F(
0873         "resVsPt_phi_F_" + ptrange[i], ";#phi residual (L1 - sim) [rad]; L1 tracks / 0.0001", 100, -0.005, 0.005);
0874
0875     h_resVsPt_eta[i] =
0876         new TH1F("resVsPt_eta_" + ptrange[i], ";#eta residual (L1 - sim); L1 tracks / 0.0002", 100, -0.01, 0.01);
0877
0878     h_resVsPt_d0[i] =
0879         new TH1F("resVsPt_d0_" + ptrange[i], ";d_{0} residual (L1 - sim) [cm]; L1 tracks / 0.0004", 100, -0.02, 0.02);
0880   }
0881
0882   // ----------------------------------------------------------------------------------------------------------------
0883   // more resolution vs eta
0884
0885   TH1F* h_resVsEta_eta[nETARANGE];
0886   TH1F* h_resVsEta_eta_L[nETARANGE];
0887   TH1F* h_resVsEta_eta_H[nETARANGE];
0888
0889   TH1F* h_resVsEta_z0[nETARANGE];
0890   TH1F* h_resVsEta_z0_L[nETARANGE];
0891   TH1F* h_resVsEta_z0_H[nETARANGE];
0892
0893   TH1F* h_resVsEta_phi[nETARANGE];
0894   TH1F* h_resVsEta_phi_L[nETARANGE];
0895   TH1F* h_resVsEta_phi_H[nETARANGE];
0896
0897   TH1F* h_resVsEta_ptRel[nETARANGE];
0898   TH1F* h_resVsEta_ptRel_L[nETARANGE];
0899   TH1F* h_resVsEta_ptRel_H[nETARANGE];
0900
0901   TH1F* h_resVsEta_d0[nETARANGE];
0902   TH1F* h_resVsEta_d0_L[nETARANGE];
0903   TH1F* h_resVsEta_d0_H[nETARANGE];
0904
0905   for (int i = 0; i < nETARANGE; i++) {
0906     h_resVsEta_eta[i] =
0907         new TH1F("resVsEta2_eta_" + etarange[i], ";#eta residual (L1 - sim); L1 tracks / 0.0002", 100, -0.01, 0.01);
0908     h_resVsEta_eta_L[i] =
0909         new TH1F("resVsEta2_eta_L_" + etarange[i], ";#eta residual (L1 - sim); L1 tracks / 0.0002", 100, -0.01, 0.01);
0910     h_resVsEta_eta_H[i] =
0911         new TH1F("resVsEta2_eta_H_" + etarange[i], ";#eta residual (L1 - sim); L1 tracks / 0.0002", 100, -0.01, 0.01);
0912
0913     h_resVsEta_z0[i] =
0914         new TH1F("resVsEta2_z0_" + etarange[i], ";z_{0} residual (L1 - sim) [cm]; L1 tracks / 0.01", 100, -1, 1);
0915     h_resVsEta_z0_L[i] =
0916         new TH1F("resVsEta2_z0_L_" + etarange[i], ";z_{0} residual (L1 - sim) [cm]; L1 tracks / 0.01", 100, -1, 1);
0917     h_resVsEta_z0_H[i] =
0918         new TH1F("resVsEta2_z0_H_" + etarange[i], ";z_{0} residual (L1 - sim) [cm]; L1 tracks / 0.01", 100, -1, 1);
0919
0920     h_resVsEta_phi[i] = new TH1F(
0921         "resVsEta2_phi_" + etarange[i], ";#phi residual (L1 - sim) [rad]; L1 tracks / 0.0001", 100, -0.005, 0.005);
0922     h_resVsEta_phi_L[i] = new TH1F(
0923         "resVsEta2_phi_L_" + etarange[i], ";#phi residual (L1 - sim) [rad]; L1 tracks / 0.0001", 100, -0.005, 0.005);
0924     h_resVsEta_phi_H[i] = new TH1F(
0925         "resVsEta2_phi_H_" + etarange[i], ";#phi residual (L1 - sim) [rad]; L1 tracks / 0.0001", 100, -0.005, 0.005);
0926
0927     h_resVsEta_ptRel[i] = new TH1F(
0928         "resVsEta2_ptRel_" + etarange[i], ";p_{T} residual (L1 - sim) / p_{T}; L1 tracks / 0.01", 100, -0.5, 0.5);
0929     h_resVsEta_ptRel_L[i] = new TH1F(
0930         "resVsEta2_ptRel_L_" + etarange[i], ";p_{T} residual (L1 - sim) / p_{T}; L1 tracks / 0.02", 100, -0.1, 0.1);
0931     h_resVsEta_ptRel_H[i] = new TH1F(
0932         "resVsEta2_ptRel_H_" + etarange[i], ";p_{T} residual (L1 - sim) / p_{T}; L1 tracks / 0.02", 100, -0.25, 0.25);
0933
0934     h_resVsEta_d0[i] =
0935         new TH1F("resVsEta2_d0_" + etarange[i], ";d_{0} residual (L1 - sim) [cm]; L1 tracks / 0.004", 100, -0.02, 0.02);
0936     h_resVsEta_d0_L[i] = new TH1F(
0937         "resVsEta2_d0_L_" + etarange[i], ";d_{0} residual (L1 - sim) [cm]; L1 tracks / 0.004", 100, -0.02, 0.02);
0938     h_resVsEta_d0_H[i] = new TH1F(
0939         "resVsEta2_d0_H_" + etarange[i], ";d_{0} residual (L1 - sim) [cm]; L1 tracks / 0.004", 100, -0.02, 0.02);
0940   }
0941   // ----------------------------------------------------------------------------------------------------------------
0942
0943   // ----------------------------------------------------------------------------------------------------------------
0944   // additional ones for sum pt in jets
0945   /*
0946   TH1F* h_jet_tp_sumpt_vspt  = new TH1F("jet_tp_sumpt_vspt",  ";sum(TP p_{T}) [GeV]; Gen jets / 5.0 GeV", 20, 0, 200.0);
0947   TH1F* h_jet_trk_sumpt_vspt = new TH1F("jet_trk_sumpt_vspt", ";sum(TP p_{T}) [GeV]; Gen jets / 5.0 GeV", 20, 0, 200.0);
0948   TH1F* h_jet_matchtrk_sumpt_vspt = new TH1F("jet_matchtrk_sumpt_vspt", ";sum(TP p_{T}) [GeV]; Gen jets / 5.0 GeV", 20, 0, 200.0);
0949
0950   TH1F* h_jet_tp_sumpt_vseta  = new TH1F("jet_tp_sumpt_vseta",  ";Gen jet #eta; Gen jets / 0.2", 24, -2.4, 2.4);
0951   TH1F* h_jet_trk_sumpt_vseta = new TH1F("jet_trk_sumpt_vseta", ";Gen jet #eta; Gen jets / 0.2", 24, -2.4, 2.4);
0952   TH1F* h_jet_matchtrk_sumpt_vseta = new TH1F("jet_matchtrk_sumpt_vseta", ";Gen jet #eta; Gen jets / 0.2", 24, -2.4, 2.4);
0953
0954   h_jet_tp_sumpt_vseta->Sumw2();
0955   h_jet_tp_sumpt_vspt->Sumw2();
0956   h_jet_trk_sumpt_vseta->Sumw2();
0957   h_jet_trk_sumpt_vspt->Sumw2();
0958   h_jet_matchtrk_sumpt_vseta->Sumw2();
0959   h_jet_matchtrk_sumpt_vspt->Sumw2();
0960   */
0961
0962   // ----------------------------------------------------------------------------------------------------------------
0963   // number of tracks per event
0964
0965   // all tracks
0966   TH1F* h_ntrk_pt2 = new TH1F("ntrk_pt2", ";# tracks (p_{T} > 2 GeV) / event; Events", 400, 0, 400.0);
0967   TH1F* h_ntrk_pt3 = new TH1F("ntrk_pt3", ";# tracks (p_{T} > 3 GeV) / event; Events", 300, 0, 300.0);
0968   TH1F* h_ntrk_pt10 = new TH1F("ntrk_pt10", ";# tracks (p_{T} > 10 GeV) / event; Events", 100, 0, 100.0);
0969
0970   // tracks flagged as genuine (this would include duplicates (?))
0971   TH1F* h_ntrk_genuine_pt2 =
0972       new TH1F("ntrk_genuine_pt2", ";# genuine tracks (p_{T} > 2 GeV) / event; Events", 400, 0, 400.0);
0973   TH1F* h_ntrk_genuine_pt3 =
0974       new TH1F("ntrk_genuine_pt3", ";# genuine tracks (p_{T} > 3 GeV) / event; Events", 300, 0, 300.0);
0975   TH1F* h_ntrk_genuine_pt10 =
0976       new TH1F("ntrk_genuine_pt10", ";# genuine tracks (p_{T} > 10 GeV) / event; Events", 100, 0, 100.0);
0977
0978   // Max N tracks from a sector per event
0979   TH1F* h_ntrkPerSector_all =
0980       new TH1F("ntrkPerSector_all", ";Max. # tracks from a sector / event; Events", 50, 0, 100.0);
0981   TH1F* h_ntrkPerSector_pt2 =
0982       new TH1F("ntrkPerSector_pt2", ";Max. # tracks from a sector (p_{T} > 2 GeV) / event; Events", 50, 0, 100.0);
0983   TH1F* h_ntrkPerSector_pt3 =
0984       new TH1F("ntrkPerSector_pt3", ";Max. # tracks from a sector (p_{T} > 3 GeV) / event; Events", 50, 0, 100.0);
0985   TH1F* h_ntrkPerSector_pt4 =
0986       new TH1F("ntrkPerSector_pt4", ";Max. # tracks from a sector (p_{T} > 10 GeV) / event; Events", 50, 0, 100.0);
0987
0988   // number of tracks vs. efficiency (eta, pT)
0989   TH1F* h_trk_pt = new TH1F("trk_pt", Form(";Track p_{T} (GeV);Tracks / 0.5 GeV"), 200, 0., 100.);
0990   TH1F* h_trk_eta = new TH1F("trk_eta", Form(";Track #eta;Tracks / 0.026"), 200, -2.6, 2.6);
0991
0992   // ----------------------------------------------------------------------------------------------------------------
0993   //        * * * * *     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     * * * * *
0994   // ----------------------------------------------------------------------------------------------------------------
0995
0996   int nevt = tree->GetEntries();
0997   cout << "number of events = " << nevt << endl;
0998
0999   // ----------------------------------------------------------------------------------------------------------------
1000   // event loop
1001   for (int i = 0; i < nevt; i++) {
1002     tree->GetEntry(i, 0);
1003
1004     /*
1005     // ----------------------------------------------------------------------------------------------------------------
1006     // sumpt in jets
1007     if (TP_select_injet > 0) {
1008       for (int ij=0; ij<(int)jet_tp_sumpt->size(); ij++) {
1009
1010     float fraction = 0;
1011     float fractionMatch = 0;
1012     if (jet_tp_sumpt->at(ij) > 0) {
1013       fraction = jet_trk_sumpt->at(ij)/jet_tp_sumpt->at(ij);
1014       fractionMatch = jet_matchtrk_sumpt->at(ij)/jet_tp_sumpt->at(ij);
1015     }
1016
1017     h_jet_tp_sumpt_vspt->Fill(jet_tp_sumpt->at(ij),1.0);
1018     h_jet_trk_sumpt_vspt->Fill(jet_tp_sumpt->at(ij),fraction);
1019     h_jet_matchtrk_sumpt_vspt->Fill(jet_tp_sumpt->at(ij),fractionMatch);
1020
1021     h_jet_tp_sumpt_vseta->Fill(jet_eta->at(ij),1.0);
1022     h_jet_trk_sumpt_vseta->Fill(jet_eta->at(ij),fraction);
1023     h_jet_matchtrk_sumpt_vseta->Fill(jet_eta->at(ij),fractionMatch);
1024       }
1025     }
1026     */
1027
1028     // ----------------------------------------------------------------------------------------------------------------
1029     // track loop for total rates & fake rates.
1030
1031     int ntrkevt_pt2 = 0;
1032     int ntrkevt_pt3 = 0;
1033     int ntrkevt_pt10 = 0;
1034
1035     int ntrkevt_genuine_pt2 = 0;
1036     int ntrkevt_genuine_pt3 = 0;
1037     int ntrkevt_genuine_pt10 = 0;
1038
1039     vector<unsigned int> nTrksPerSector_all(9, 0);
1040     vector<unsigned int> nTrksPerSector_pt2(9, 0);
1041     vector<unsigned int> nTrksPerSector_pt3(9, 0);
1042     vector<unsigned int> nTrksPerSector_pt4(9, 0);
1043
1044     for (int it = 0; it < (int)trk_pt->size(); it++) {
1045       // ----------------------------------------------------------------------------------------------------------------
1046       // track properties
1047
1048       // ----------------------------------------------------------------------------------------------------------------
1049       // Fill number of tracks vs track param
1050       h_trk_pt->Fill(trk_pt->at(it));
1051       h_trk_eta->Fill(trk_eta->at(it));
1052
1053       // fill all trk chi2 & chi2/dof histograms, including for chi2 r-phi and chi2 r-z
1054       int ndof = 2 * trk_nstub->at(it) - 4;
1055       float chi2 = trk_chi2->at(it);
1056       float chi2dof = (float)chi2 / ndof;
1057       float chi2rphi = trk_chi2rphi->at(it);
1058       float chi2rphidof = (float)chi2rphi / ndof;
1059       float chi2rz = trk_chi2rz->at(it);
1060       float chi2rzdof = (float)chi2rz / ndof;
1061
1062       // create overflow bins by restricting range of chi2
1063       int chi2Overflow = 100;
1064       int chi2DOFOverflow = 20;  //apprx chi2Overflow / avg. nstubs
1065       double buffer = 0.1;
1066
1067       if (chi2 > chi2Overflow)
1068         chi2 = chi2Overflow - buffer;
1069       if (chi2dof > chi2DOFOverflow)
1070         chi2dof = chi2DOFOverflow - buffer;
1071       if (chi2rphi > chi2Overflow)
1072         chi2rphi = chi2Overflow - buffer;
1073       if (chi2rphidof > chi2DOFOverflow)
1074         chi2rphidof = chi2DOFOverflow - buffer;
1075       if (chi2rz > chi2Overflow)
1076         chi2rz = chi2Overflow - buffer;
1077       if (chi2rzdof > chi2DOFOverflow)
1078         chi2rzdof = chi2DOFOverflow - buffer;
1079
1080       if (trk_pt->at(it) > TP_minPt) {  //TRK pt > TP_minPt
1081
1082         h_trk_chi2->Fill(chi2);
1083         h_trk_chi2_dof->Fill(chi2dof);
1084
1085         h_trk_chi2rphi->Fill(chi2rphi);
1086         h_trk_chi2rphi_dof->Fill(chi2rphidof);
1087
1088         h_trk_chi2rz->Fill(chi2rz);
1089         h_trk_chi2rz_dof->Fill(chi2rzdof);
1090
1091       }  //end TRK pt > TP_minPt
1092
1093       // ----------------------------------------------------------------------------------------------------------------
1094       // look at track rate and fake rate, etc.
1095
1096       // only look at tracks in (ttbar) jets ?
1097       if (TP_select_injet > 0) {
1098         if (TP_select_injet == 1 && trk_injet->at(it) == 0)
1099           continue;
1100         if (TP_select_injet == 2 && trk_injet_highpt->at(it) == 0)
1101           continue;
1102         if (TP_select_injet == 3 && trk_injet_vhighpt->at(it) == 0)
1103           continue;
1104       }
1105       ntrk++;
1106       if (trk_pt->at(it) >= 0.0)
1107         ++nTrksPerSector_all.at(trk_phiSector->at(it) % 9);
1108       if (std::abs(trk_eta->at(it)) > TP_maxEta)
1109         continue;
1110       if (trk_pt->at(it) < TP_minPt)
1111         continue;
1112
1113       // Uncomment these cuts to see effect on rate & fake rate.
1114       //int ndof = 2*trk_nstub->at(it)-4;
1115       //if (trk_chi2->at(it) > L1Tk_maxChi2) continue;
1116       //if (trk_chi2->at(it)/ndof > L1Tk_maxChi2dof) continue;
1117       //if (trk_nstub->at(it) < L1Tk_minNstub) continue;
1118
1119       // Tracklet & Hybrid have 9 sectors, but TMTT has 18 (with sectors 0 & 1 in nonant 0 etc).
1120       // As don't know here with algo used, "% 9" added to prevent crash, but not correct for TMTT.
1121       if (trk_pt->at(it) > 2.0)
1122         ++nTrksPerSector_pt2.at(trk_phiSector->at(it) % 9);
1123       if (trk_pt->at(it) > 3.0)
1124         ++nTrksPerSector_pt3.at(trk_phiSector->at(it) % 9);
1125       if (trk_pt->at(it) > 4.0)
1126         ++nTrksPerSector_pt4.at(trk_phiSector->at(it) % 9);
1127
1128       if (trk_pt->at(it) > 2.0) {
1129         ntrk_pt2++;
1130         ntrkevt_pt2++;
1131         h_trk_all_vspt->Fill(trk_pt->at(it));
1132         if (trk_genuine->at(it) == 1) {
1133           ntrkevt_genuine_pt2++;
1134           h_trk_genuine_vspt->Fill(trk_pt->at(it));
1135         } else
1136           h_trk_notgenuine_vspt->Fill(trk_pt->at(it));
1137         if (trk_loose->at(it) == 1)
1138           h_trk_loose_vspt->Fill(trk_pt->at(it));
1139         else
1140           h_trk_notloose_vspt->Fill(trk_pt->at(it));
1141       }
1142       if (trk_pt->at(it) > 3.0) {
1143         ntrk_pt3++;
1144         ntrkevt_pt3++;
1145         if (trk_genuine->at(it) == 1)
1146           ntrkevt_genuine_pt3++;
1147       }
1148       if (trk_pt->at(it) > 10.0) {
1149         ntrk_pt10++;
1150         ntrkevt_pt10++;
1151         if (trk_genuine->at(it) == 1)
1152           ntrkevt_genuine_pt10++;
1153       }
1154
1155       // ----------------------------------------------------------------------------------------------------------------
1156       // Fill tracklet propogation efficiency histo
1157
1158       // create an 11-bit long iterable from lhits and dhits
1159       int num_layers = 6;
1160       int num_discs = 5;
1161       int lhits = trk_lhits->at(it);
1162       int dhits = trk_dhits->at(it);
1163       std::vector<int> layers = {};
1164       for (int layer_index = 0; layer_index < num_layers + num_discs; layer_index++) {
1165         if (layer_index < num_layers) {
1166           layers.push_back(lhits % 10);
1167           lhits /= 10;
1168         } else {
1169           layers.push_back(dhits % 10);
1170           dhits /= 10;
1171         }
1172       }
1173
1174       for (unsigned int layer = 0; layer < layers.size(); layer++) {
1175         if (layers.at(layer)) {                                         // if there was a hit at this layer...
1176           h_trk_tracklet_hits->Fill(std::abs(trk_eta->at(it)), layer);  // ...fill this bin with the layer of the track.
1177         }
1178       }
1179     }
1180
1181     h_ntrk_pt2->Fill(ntrkevt_pt2);
1182     h_ntrk_pt3->Fill(ntrkevt_pt3);
1183     h_ntrk_pt10->Fill(ntrkevt_pt10);
1184
1185     h_ntrk_genuine_pt2->Fill(ntrkevt_genuine_pt2);
1186     h_ntrk_genuine_pt3->Fill(ntrkevt_genuine_pt3);
1187     h_ntrk_genuine_pt10->Fill(ntrkevt_genuine_pt10);
1188
1189     h_ntrkPerSector_all->Fill(*std::max_element(nTrksPerSector_all.begin(), nTrksPerSector_all.end()));
1190     h_ntrkPerSector_pt2->Fill(*std::max_element(nTrksPerSector_pt2.begin(), nTrksPerSector_pt2.end()));
1191     h_ntrkPerSector_pt3->Fill(*std::max_element(nTrksPerSector_pt3.begin(), nTrksPerSector_pt3.end()));
1192     h_ntrkPerSector_pt4->Fill(*std::max_element(nTrksPerSector_pt4.begin(), nTrksPerSector_pt4.end()));
1193
1194     // ----------------------------------------------------------------------------------------------------------------
1195     // tracking particle loop
1196     for (int it = 0; it < (int)tp_pt->size(); it++) {
1197       // only look at TPs in (ttbar) jets ?
1198       if (TP_select_injet > 0) {
1199         if (TP_select_injet == 1 && tp_injet->at(it) == 0)
1200           continue;
1201         if (TP_select_injet == 2 && tp_injet_highpt->at(it) == 0)
1202           continue;
1203         if (TP_select_injet == 3 && tp_injet_vhighpt->at(it) == 0)
1204           continue;
1205       }
1206
1207       // cut on PDG ID at plot stage?
1208       if (TP_select_pdgid != 0) {
1209         if (abs(tp_pdgid->at(it)) != abs(TP_select_pdgid))
1210           continue;
1211       }
1212
1213       // kinematic cuts
1214       if (std::abs(tp_dxy->at(it)) > TP_maxDxy)
1215         continue;
1216       if (std::abs(tp_d0->at(it)) > TP_maxD0)
1217         continue;
1218       if (tp_pt->at(it) < 0.2)
1219         continue;
1220       if (tp_pt->at(it) > TP_maxPt)
1221         continue;
1222       if (std::abs(tp_eta->at(it)) > TP_maxEta)
1223         continue;
1224
1225       // total track rates
1226       if (tp_pt->at(it) > TP_minPt) {
1227         if (tp_pt->at(it) > 2.0) {
1228           ntp_pt2++;
1229           h_tp_vspt->Fill(tp_pt->at(it));
1230           // duplicate rate
1231           if (tp_nmatch->at(it) > 1) {
1232             for (int inm = 1; inm < tp_nmatch->at(it); inm++)
1233               h_trk_duplicate_vspt->Fill(matchtrk_pt->at(it));
1234           }
1235         }
1236         if (tp_pt->at(it) > 3.0)
1237           ntp_pt3++;
1238         if (tp_pt->at(it) > 10.0)
1239           ntp_pt10++;
1240       }
1241
1242       // cut on event ID (eventid=0 means the TP is from the primary interaction, so *not* selecting only eventid=0 means including stuff from pileup)
1243       if (TP_select_eventid == 0 && tp_eventid->at(it) != 0)
1244         continue;
1245
1246       // look at failure scenarios?
1247       if (useDeadRegion) {
1248         if (tp_phi->at(it) < 0 || tp_phi->at(it) > 1)
1249           continue;
1250       }
1251
1252       h_tp_pt->Fill(tp_pt->at(it));
1253       if (tp_pt->at(it) < 8.0)
1254         h_tp_pt_L->Fill(tp_pt->at(it));
1255       else
1256         h_tp_pt_H->Fill(tp_pt->at(it));
1257       if (tp_pt->at(it) < 8.0 && std::abs(tp_eta->at(it)) < 1.0)
1258         h_tp_pt_LC->Fill(tp_pt->at(it));
1259
1260       if (tp_pt->at(it) > TP_minPt) {
1261         if (std::abs(tp_eta->at(it)) < 1.0)
1262           n_all_eta1p0++;
1263         else if (std::abs(tp_eta->at(it)) < 1.75)
1264           n_all_eta1p75++;
1265         else
1266           n_all_eta2p5++;
1267
1268         if (std::abs(tp_pt->at(it)) > 2.0)
1269           n_all_ptg2++;
1270         if (std::abs(tp_pt->at(it)) > 2.0 && std::abs(tp_pt->at(it)) < 8.0)
1271           n_all_pt2to8++;
1272         if (std::abs(tp_pt->at(it)) > 8.0)
1273           n_all_ptg8++;
1274         if (std::abs(tp_pt->at(it)) > 40.0)
1275           n_all_ptg40++;
1276
1277         h_tp_eta->Fill(tp_eta->at(it));
1278         h_tp_phi->Fill(tp_phi->at(it));
1279         h_tp_z0->Fill(tp_z0->at(it));
1280         h_tp_d0->Fill(tp_d0->at(it));
1281         h_tp_absd0->Fill(std::abs(tp_d0->at(it)));
1282         if (std::abs(tp_eta->at(it)) < 2.0)
1283           h_tp_absd0_eta2->Fill(std::abs(tp_d0->at(it)));
1284         if (std::abs(tp_eta->at(it)) < 2.0 && tp_pt->at(it) > 3.0)
1285           h_tp_absd0_eta2_pt3->Fill(std::abs(tp_d0->at(it)));
1286
1287         if (tp_pt->at(it) < 3.0)
1288           h_tp_eta_23->Fill(tp_eta->at(it));
1289         else if (tp_pt->at(it) < 5.0)
1290           h_tp_eta_35->Fill(tp_eta->at(it));
1291         else
1292           h_tp_eta_5->Fill(tp_eta->at(it));
1293
1294         if (tp_pt->at(it) < 8.0) {
1295           h_tp_eta_L->Fill(tp_eta->at(it));
1296           h_tp_z0_L->Fill(tp_z0->at(it));
1297         } else {
1298           h_tp_eta_H->Fill(tp_eta->at(it));
1299           h_tp_z0_H->Fill(tp_z0->at(it));
1300         }
1301       }
1302
1303       // ----------------------------------------------------------------------------------------------------------------
1304       // was the tracking particle matched to a L1 track?
1305       if (tp_nmatch->at(it) < 1)
1306         continue;
1307
1308       // ----------------------------------------------------------------------------------------------------------------
1309       // use only tracks with min X stubs
1310       if (matchtrk_nstub->at(it) < L1Tk_minNstub)
1311         continue;
1312
1313       int thisseed = matchtrk_seed->at(it);
1314       if (thisseed > 25)
1315         thisseed = thisseed - 20;
1316       if ((L1Tk_seed != 0) && (thisseed != L1Tk_seed))
1317         continue;
1318
1319       // ----------------------------------------------------------------------------------------------------------------
1320       // fill matchtrk chi2 & chi2/dof histograms before making chi2 cut
1321
1322       int ndof = 2 * matchtrk_nstub->at(it) - 4;
1323       float chi2 = matchtrk_chi2->at(it);
1324       float chi2dof = (float)chi2 / ndof;
1325       float chi2rphi = matchtrk_chi2rphi->at(it);
1326       float chi2rphidof = (float)chi2rphi / ndof;
1327       float chi2rz = matchtrk_chi2rz->at(it);
1328       float chi2rzdof = (float)chi2rz / ndof;
1329
1330       // create overflow bins by restricting range of chi2
1331       int chi2Overflow = 100;
1332       int chi2DOFOverflow = 20;  //apprx chi2Overflow / avg. nstubs
1333       double buffer = 0.1;
1334
1335       if (chi2 > chi2Overflow)
1336         chi2 = chi2Overflow - buffer;
1337       if (chi2dof > chi2DOFOverflow)
1338         chi2dof = chi2DOFOverflow - buffer;
1339       if (chi2rphi > chi2Overflow)
1340         chi2rphi = chi2Overflow - buffer;
1341       if (chi2rphidof > chi2DOFOverflow)
1342         chi2rphidof = chi2DOFOverflow - buffer;
1343       if (chi2rz > chi2Overflow)
1344         chi2rz = chi2Overflow - buffer;
1345       if (chi2rzdof > chi2DOFOverflow)
1346         chi2rzdof = chi2DOFOverflow - buffer;
1347
1348       if (tp_pt->at(it) > TP_minPt) {  //TP pt > TP_minPt
1349
1350         h_match_trk_chi2->Fill(chi2);
1351         h_match_trk_chi2_dof->Fill(chi2dof);
1352
1353         h_match_trk_chi2rphi->Fill(chi2rphi);
1354         h_match_trk_chi2rphi_dof->Fill(chi2rphidof);
1355
1356         h_match_trk_chi2rz->Fill(chi2rz);
1357         h_match_trk_chi2rz_dof->Fill(chi2rzdof);
1358
1359         // central eta
1360         if (std::abs(tp_eta->at(it)) < 0.8) {
1361           if (tp_pt->at(it) < 8.0) {
1362             h_match_trk_chi2_C_L->Fill(chi2);
1363             h_match_trk_chi2_dof_C_L->Fill(chi2dof);
1364           } else {
1365             h_match_trk_chi2_C_H->Fill(chi2);
1366             h_match_trk_chi2_dof_C_H->Fill(chi2dof);
1367           }
1368         }
1369         // intermediate eta
1370         else if (std::abs(tp_eta->at(it)) < 1.6 && std::abs(tp_eta->at(it)) >= 0.8) {
1371           if (tp_pt->at(it) < 8.0) {
1372             h_match_trk_chi2_I_L->Fill(chi2);
1373             h_match_trk_chi2_dof_I_L->Fill(chi2dof);
1374           } else {
1375             h_match_trk_chi2_I_H->Fill(chi2);
1376             h_match_trk_chi2_dof_I_H->Fill(chi2dof);
1377           }
1378         }
1379         // forward eta
1380         else if (std::abs(tp_eta->at(it)) >= 1.6) {
1381           if (tp_pt->at(it) < 8.0) {
1382             h_match_trk_chi2_F_L->Fill(chi2);
1383             h_match_trk_chi2_dof_F_L->Fill(chi2dof);
1384           } else {
1385             h_match_trk_chi2_F_H->Fill(chi2);
1386             h_match_trk_chi2_dof_F_H->Fill(chi2dof);
1387           }
1388         }
1389       }  //end TP pt > TP_minPt
1390
1391       // ----------------------------------------------------------------------------------------------------------------
1392       // cut on chi2?
1393       if (matchtrk_chi2->at(it) > L1Tk_maxChi2)
1394         continue;
1395       if (matchtrk_chi2->at(it) / ndof > L1Tk_maxChi2dof)
1396         continue;
1397
1398       // use tight quality cut selection?
1399       /*
1400       if (matchtrk_nstub->at(it)==4) {
1401     if (std::abs(matchtrk_eta->at(it))<2.2 && matchtrk_consistency->at(it)>10) continue;
1402     else if (std::abs(matchtrk_eta->at(it))>2.2 && chi2dof>5.0) continue;
1403       }
1404       if (matchtrk_pt->at(it)>10.0 && chi2dof>5.0) continue;
1405       */
1406
1407       // ----------------------------------------------------------------------------------------------------------------
1408       // more plots
1409
1410       // fill matched track histograms
1411       h_match_tp_pt->Fill(tp_pt->at(it));
1412       if (tp_pt->at(it) < 8.0)
1413         h_match_tp_pt_L->Fill(tp_pt->at(it));
1414       else
1415         h_match_tp_pt_H->Fill(tp_pt->at(it));
1416       if (tp_pt->at(it) < 8.0 && std::abs(tp_eta->at(it)) < 1.0)
1417         h_match_tp_pt_LC->Fill(tp_pt->at(it));
1418
1419       if (tp_pt->at(it) > TP_minPt) {
1420         h_match_tp_eta->Fill(tp_eta->at(it));
1421         h_match_tp_phi->Fill(tp_phi->at(it));
1422         h_match_tp_z0->Fill(tp_z0->at(it));
1423         h_match_tp_d0->Fill(tp_d0->at(it));
1424         h_match_tp_absd0->Fill(std::abs(tp_d0->at(it)));
1425         if (std::abs(tp_eta->at(it)) < 2.0)
1426           h_match_tp_absd0_eta2->Fill(std::abs(tp_d0->at(it)));
1427         if (std::abs(tp_eta->at(it)) < 2.0 && tp_pt->at(it) > 3.0)
1428           h_match_tp_absd0_eta2_pt3->Fill(std::abs(tp_d0->at(it)));
1429
1430         if (std::abs(tp_eta->at(it)) < 1.0)
1431           n_match_eta1p0++;
1432         else if (std::abs(tp_eta->at(it)) < 1.75)
1433           n_match_eta1p75++;
1434         else
1435           n_match_eta2p5++;
1436
1437         if (std::abs(tp_pt->at(it)) > 2.0)
1438           n_match_ptg2++;
1439         if (std::abs(tp_pt->at(it)) > 2.0 && std::abs(tp_pt->at(it)) < 8.0)
1440           n_match_pt2to8++;
1441         if (std::abs(tp_pt->at(it)) > 8.0)
1442           n_match_ptg8++;
1443         if (std::abs(tp_pt->at(it)) > 40.0)
1444           n_match_ptg40++;
1445
1446         if (tp_pt->at(it) < 3.0)
1447           h_match_tp_eta_23->Fill(tp_eta->at(it));
1448         else if (tp_pt->at(it) < 5.0)
1449           h_match_tp_eta_35->Fill(tp_eta->at(it));
1450         else
1451           h_match_tp_eta_5->Fill(tp_eta->at(it));
1452
1453         if (tp_pt->at(it) < 8.0) {
1454           h_match_tp_eta_L->Fill(tp_eta->at(it));
1455           h_match_tp_z0_L->Fill(tp_z0->at(it));
1456         } else {
1457           h_match_tp_z0_H->Fill(tp_z0->at(it));
1458           h_match_tp_eta_H->Fill(tp_eta->at(it));
1459         }
1460       }
1461
1462       // for the following, only consider TPs with pt > TP_minPt
1463       if (tp_pt->at(it) < TP_minPt)
1464         continue;
1465
1466       // fill nstub histograms
1467       h_match_trk_nstub->Fill(matchtrk_nstub->at(it));
1468       if (std::abs(tp_eta->at(it)) < 0.8)
1469         h_match_trk_nstub_C->Fill(matchtrk_nstub->at(it));
1470       else if (std::abs(tp_eta->at(it)) < 1.6 && std::abs(tp_eta->at(it)) >= 0.8)
1471         h_match_trk_nstub_I->Fill(matchtrk_nstub->at(it));
1472       else if (std::abs(tp_eta->at(it)) >= 1.6)
1473         h_match_trk_nstub_F->Fill(matchtrk_nstub->at(it));
1474
1475       // ----------------------------------------------------------------------------------------------------------------
1476       // fill resolution histograms
1477
1478       h_res_pt->Fill(matchtrk_pt->at(it) - tp_pt->at(it));
1479       h_res_ptRel->Fill((matchtrk_pt->at(it) - tp_pt->at(it)) / tp_pt->at(it));
1480       h_res_eta->Fill(matchtrk_eta->at(it) - tp_eta->at(it));
1481       h_res_phi->Fill(matchtrk_phi->at(it) - tp_phi->at(it));
1482       h_res_z0->Fill(matchtrk_z0->at(it) - tp_z0->at(it));
1483       if (matchtrk_d0->at(it) < 999.)
1484         h_res_d0->Fill(matchtrk_d0->at(it) - tp_d0->at(it));
1485
1486       if (std::abs(tp_eta->at(it)) < 0.8)
1487         h_res_z0_C->Fill(matchtrk_z0->at(it) - tp_z0->at(it));
1488       else if (std::abs(tp_eta->at(it)) < 1.6 && std::abs(tp_eta->at(it)) >= 0.8)
1489         h_res_z0_I->Fill(matchtrk_z0->at(it) - tp_z0->at(it));
1490       else if (std::abs(tp_eta->at(it)) >= 1.6)
1491         h_res_z0_F->Fill(matchtrk_z0->at(it) - tp_z0->at(it));
1492
1493       if (tp_pt->at(it) < 8.0) {
1494         h_res_z0_L->Fill(matchtrk_z0->at(it) - tp_z0->at(it));
1495         if (std::abs(tp_eta->at(it)) < 1.0)
1496           h_res_z0_C_L->Fill(matchtrk_z0->at(it) - tp_z0->at(it));
1497         else
1498           h_res_z0_F_L->Fill(matchtrk_z0->at(it) - tp_z0->at(it));
1499       } else {
1500         h_res_z0_H->Fill(matchtrk_z0->at(it) - tp_z0->at(it));
1501         if (std::abs(tp_eta->at(it)) < 1.0)
1502           h_res_z0_C_H->Fill(matchtrk_z0->at(it) - tp_z0->at(it));
1503         else
1504           h_res_z0_F_H->Fill(matchtrk_z0->at(it) - tp_z0->at(it));
1505       }
1506
1507       if (matchtrk_d0->at(it) < 999.) {
1508         if (std::abs(tp_eta->at(it)) < 0.8)
1509           h_res_d0_C->Fill(matchtrk_d0->at(it) - tp_d0->at(it));
1510         else if (std::abs(tp_eta->at(it)) < 1.6 && std::abs(tp_eta->at(it)) >= 0.8)
1511           h_res_d0_I->Fill(matchtrk_d0->at(it) - tp_d0->at(it));
1512         else if (std::abs(tp_eta->at(it)) >= 1.6)
1513           h_res_d0_F->Fill(matchtrk_d0->at(it) - tp_d0->at(it));
1514
1515         if (tp_pt->at(it) < 8.0) {
1516           h_res_d0_L->Fill(matchtrk_d0->at(it) - tp_d0->at(it));
1517           if (std::abs(tp_eta->at(it)) < 1.0)
1518             h_res_d0_C_L->Fill(matchtrk_d0->at(it) - tp_d0->at(it));
1519           else
1520             h_res_d0_F_L->Fill(matchtrk_d0->at(it) - tp_d0->at(it));
1521         } else {
1522           h_res_d0_H->Fill(matchtrk_d0->at(it) - tp_d0->at(it));
1523           if (std::abs(tp_eta->at(it)) < 1.0)
1524             h_res_d0_C_H->Fill(matchtrk_d0->at(it) - tp_d0->at(it));
1525           else
1526             h_res_d0_F_H->Fill(matchtrk_d0->at(it) - tp_d0->at(it));
1527         }
1528       }
1529
1530       // ----------------------------------------------------------------------------------------------------------------
1531       // fill resolution vs. pt histograms
1532       for (int im = 0; im < nRANGE; im++) {
1533         if ((tp_pt->at(it) > (float)im * 5.0) && (tp_pt->at(it) < (float)(im + 1) * 5.0)) {
1534           h_resVsPt_pt[im]->Fill(matchtrk_pt->at(it) - tp_pt->at(it));
1535           h_resVsPt_ptRel[im]->Fill((matchtrk_pt->at(it) - tp_pt->at(it)) / tp_pt->at(it));
1536           h_resVsPt_eta[im]->Fill(matchtrk_eta->at(it) - tp_eta->at(it));
1537           h_resVsPt_phi[im]->Fill(matchtrk_phi->at(it) - tp_phi->at(it));
1538           h_resVsPt_z0[im]->Fill(matchtrk_z0->at(it) - tp_z0->at(it));
1539
1540           h_absResVsPt_pt[im]->Fill(std::abs(matchtrk_pt->at(it) - tp_pt->at(it)));
1541           h_absResVsPt_ptRel[im]->Fill(std::abs((matchtrk_pt->at(it) - tp_pt->at(it))) / tp_pt->at(it));
1542           h_absResVsPt_z0[im]->Fill(std::abs(matchtrk_z0->at(it) - tp_z0->at(it)));
1543           h_absResVsPt_phi[im]->Fill(std::abs(matchtrk_phi->at(it) - tp_phi->at(it)));
1544           h_absResVsPt_eta[im]->Fill(std::abs(matchtrk_eta->at(it) - tp_eta->at(it)));
1545
1546           if (std::abs(tp_eta->at(it)) < 0.8) {
1547             h_resVsPt_pt_C[im]->Fill(matchtrk_pt->at(it) - tp_pt->at(it));
1548             h_resVsPt_ptRel_C[im]->Fill((matchtrk_pt->at(it) - tp_pt->at(it)) / tp_pt->at(it));
1549             h_resVsPt_z0_C[im]->Fill(matchtrk_z0->at(it) - tp_z0->at(it));
1550             h_resVsPt_phi_C[im]->Fill(matchtrk_phi->at(it) - tp_phi->at(it));
1551           } else if (std::abs(tp_eta->at(it)) < 1.6 && std::abs(tp_eta->at(it)) >= 0.8) {
1552             h_resVsPt_pt_I[im]->Fill(matchtrk_pt->at(it) - tp_pt->at(it));
1553             h_resVsPt_ptRel_I[im]->Fill((matchtrk_pt->at(it) - tp_pt->at(it)) / tp_pt->at(it));
1554             h_resVsPt_z0_I[im]->Fill(matchtrk_z0->at(it) - tp_z0->at(it));
1555             h_resVsPt_phi_I[im]->Fill(matchtrk_phi->at(it) - tp_phi->at(it));
1556           } else if (std::abs(tp_eta->at(it)) >= 1.6) {
1557             h_resVsPt_pt_F[im]->Fill(matchtrk_pt->at(it) - tp_pt->at(it));
1558             h_resVsPt_ptRel_F[im]->Fill((matchtrk_pt->at(it) - tp_pt->at(it)) / tp_pt->at(it));
1559             h_resVsPt_z0_F[im]->Fill(matchtrk_z0->at(it) - tp_z0->at(it));
1560             h_resVsPt_phi_F[im]->Fill(matchtrk_phi->at(it) - tp_phi->at(it));
1561           }
1562
1563           if (matchtrk_d0->at(it) < 999) {
1564             h_resVsPt_d0[im]->Fill(matchtrk_d0->at(it) - tp_d0->at(it));
1565             h_absResVsPt_d0[im]->Fill(std::abs(matchtrk_d0->at(it) - tp_d0->at(it)));
1566           }
1567         }
1568       }
1569
1570       for (int im = 3; im < nRANGE_L + 3; im++) {
1571         if ((tp_pt->at(it) > (float)im * 0.5) && (tp_pt->at(it) <= (float)(im + 1) * 0.5)) {
1572           h_absResVsPt_pt_L[im - 3]->Fill(std::abs(matchtrk_pt->at(it) - tp_pt->at(it)));
1573           h_absResVsPt_ptRel_L[im - 3]->Fill(std::abs((matchtrk_pt->at(it) - tp_pt->at(it))) / tp_pt->at(it));
1574           h_absResVsPt_z0_L[im - 3]->Fill(std::abs(matchtrk_z0->at(it) - tp_z0->at(it)));
1575           h_absResVsPt_phi_L[im - 3]->Fill(std::abs(matchtrk_phi->at(it) - tp_phi->at(it)));
1576           h_absResVsPt_eta_L[im - 3]->Fill(std::abs(matchtrk_eta->at(it) - tp_eta->at(it)));
1577           h_absResVsPt_d0_L[im - 3]->Fill(std::abs(matchtrk_d0->at(it) - tp_d0->at(it)));
1578         }
1579       }
1580
1581       // ----------------------------------------------------------------------------------------------------------------
1582       // fill resolution vs. eta histograms
1583       for (int im = 0; im < nETARANGE; im++) {
1584         if ((std::abs(tp_eta->at(it)) > (float)im * 0.1) && (std::abs(tp_eta->at(it)) < (float)(im + 1) * 0.1)) {
1585           h_resVsEta_ptRel[im]->Fill((matchtrk_pt->at(it) - tp_pt->at(it)) / tp_pt->at(it));
1586           h_resVsEta_eta[im]->Fill(matchtrk_eta->at(it) - tp_eta->at(it));
1587           h_resVsEta_phi[im]->Fill(matchtrk_phi->at(it) - tp_phi->at(it));
1588           h_resVsEta_z0[im]->Fill(matchtrk_z0->at(it) - tp_z0->at(it));
1589
1590           h_absResVsEta_ptRel[im]->Fill(std::abs((matchtrk_pt->at(it) - tp_pt->at(it))) / tp_pt->at(it));
1591           h_absResVsEta_eta[im]->Fill(std::abs(matchtrk_eta->at(it) - tp_eta->at(it)));
1592           h_absResVsEta_phi[im]->Fill(std::abs(matchtrk_phi->at(it) - tp_phi->at(it)));
1593           h_absResVsEta_z0[im]->Fill(std::abs(matchtrk_z0->at(it) - tp_z0->at(it)));
1594
1595           if (tp_pt->at(it) < 8.0) {
1596             h_resVsEta_ptRel_L[im]->Fill((matchtrk_pt->at(it) - tp_pt->at(it)) / tp_pt->at(it));
1597             h_resVsEta_eta_L[im]->Fill(matchtrk_eta->at(it) - tp_eta->at(it));
1598             h_resVsEta_z0_L[im]->Fill(matchtrk_z0->at(it) - tp_z0->at(it));
1599             h_resVsEta_phi_L[im]->Fill(matchtrk_phi->at(it) - tp_phi->at(it));
1600
1601             h_absResVsEta_ptRel_L[im]->Fill(std::abs((matchtrk_pt->at(it) - tp_pt->at(it))) / tp_pt->at(it));
1602             h_absResVsEta_eta_L[im]->Fill(std::abs(matchtrk_eta->at(it) - tp_eta->at(it)));
1603             h_absResVsEta_phi_L[im]->Fill(std::abs(matchtrk_phi->at(it) - tp_phi->at(it)));
1604             h_absResVsEta_z0_L[im]->Fill(std::abs(matchtrk_z0->at(it) - tp_z0->at(it)));
1605
1606             if (matchtrk_d0->at(it) < 999) {
1607               h_resVsEta_d0_L[im]->Fill(matchtrk_d0->at(it) - tp_d0->at(it));
1608               h_absResVsEta_d0_L[im]->Fill(std::abs(matchtrk_d0->at(it) - tp_d0->at(it)));
1609             }
1610           } else {
1611             h_resVsEta_ptRel_H[im]->Fill((matchtrk_pt->at(it) - tp_pt->at(it)) / tp_pt->at(it));
1612             h_resVsEta_eta_H[im]->Fill(matchtrk_eta->at(it) - tp_eta->at(it));
1613             h_resVsEta_z0_H[im]->Fill(matchtrk_z0->at(it) - tp_z0->at(it));
1614             h_resVsEta_phi_H[im]->Fill(matchtrk_phi->at(it) - tp_phi->at(it));
1615
1616             h_absResVsEta_ptRel_H[im]->Fill(std::abs((matchtrk_pt->at(it) - tp_pt->at(it))) / tp_pt->at(it));
1617             h_absResVsEta_eta_H[im]->Fill(std::abs(matchtrk_eta->at(it) - tp_eta->at(it)));
1618             h_absResVsEta_phi_H[im]->Fill(std::abs(matchtrk_phi->at(it) - tp_phi->at(it)));
1619             h_absResVsEta_z0_H[im]->Fill(std::abs(matchtrk_z0->at(it) - tp_z0->at(it)));
1620
1621             if (matchtrk_d0->at(it) < 999) {
1622               h_resVsEta_d0_H[im]->Fill(matchtrk_d0->at(it) - tp_d0->at(it));
1623               h_absResVsEta_d0_H[im]->Fill(std::abs(matchtrk_d0->at(it) - tp_d0->at(it)));
1624             }
1625           }
1626
1627           if (matchtrk_d0->at(it) < 999) {
1628             h_resVsEta_d0[im]->Fill(matchtrk_d0->at(it) - tp_d0->at(it));
1629             h_absResVsEta_d0[im]->Fill(std::abs(matchtrk_d0->at(it) - tp_d0->at(it)));
1630           }
1631         }
1632       }
1633
1634       // ----------------------------------------------------------------------------------------------------------------
1635       // fill resolution vs. phi histograms
1636       for (int im = 0; im < nPHIRANGE; im++) {
1637         if ((tp_phi->at(it) > (float)im * 0.2 - 3.2) && (tp_phi->at(it) < (float)(im + 1) * 0.2 - 3.2)) {
1638           h_absResVsPhi_pt[im]->Fill(std::abs(matchtrk_pt->at(it) - tp_pt->at(it)));
1639           h_absResVsPhi_ptRel[im]->Fill(std::abs((matchtrk_pt->at(it) - tp_pt->at(it))) / tp_pt->at(it));
1640         }
1641       }
1642
1643     }  // end of matched track loop
1644
1645   }  // end of event loop
1646
1647   // ----------------------------------------------------------------------------------------------------------------
1648
1649   // ----------------------------------------------------------------------------------------------------------------
1650   // Post-event-loop histogram normalization(s)
1651
1652   // Normalize tracklet efficiency by each eta slice
1653   double maxBinContents;
1654   double etaRes = trackletEffMaxEta / trackletEffEtaBins;
1655   for (double etaBin = 0; etaBin < trackletEffEtaBins;
1656        etaBin++) {  //loop through eta bin values (constants defined with relevant hist defs)
1657     maxBinContents = 0;
1658     std::vector<double> binContents = {};
1659     for (int layer = 0; layer < numLayers; layer++) {
1660       binContents.push_back(h_trk_tracklet_hits->GetBinContent(etaBin + 1, layer + 1));
1661       maxBinContents = std::max(maxBinContents, binContents.back());
1662     }
1663     float binWeight;
1664     for (int layer = 0; layer < numLayers; layer++) {
1665       binWeight = (maxBinContents == 0) ? 0 : binContents.at(layer) / maxBinContents;
1666       h_trk_tracklet_eff->Fill((etaBin + 0.5) * etaRes, layer, binWeight);
1667     }
1668   }
1669
1670   // Adjust tracklet hits and efficiency histograms for aesthetics
1671
1672   // ----------------------------------------------------------------------------------------------------------------
1673   // 2D plots
1674   // ----------------------------------------------------------------------------------------------------------------
1675
1676   TH1F* h2_resVsPt_pt = new TH1F("resVsPt2_pt", ";Tracking particle p_{T} [GeV]; p_{T} resolution [GeV]", 20, 0, 100);
1677   TH1F* h2_resVsPt_pt_C =
1678       new TH1F("resVsPt2_pt_C", ";Tracking particle p_{T} [GeV]; p_{T} resolution [GeV]", 20, 0, 100);
1679   TH1F* h2_resVsPt_pt_I =
1680       new TH1F("resVsPt2_pt_I", ";Tracking particle p_{T} [GeV]; p_{T} resolution [GeV]", 20, 0, 100);
1681   TH1F* h2_resVsPt_pt_F =
1682       new TH1F("resVsPt2_pt_F", ";Tracking particle p_{T} [GeV]; p_{T} resolution [GeV]", 20, 0, 100);
1683
1684   TH1F* h2_resVsPt_ptRel =
1685       new TH1F("resVsPt2_ptRel", ";Tracking particle p_{T} [GeV]; p_{T} resolution / p_{T}", 20, 0, 100);
1686   TH1F* h2_resVsPt_ptRel_C =
1687       new TH1F("resVsPt2_ptRel_C", ";Tracking particle p_{T} [GeV]; p_{T} resolution / p_{T}", 20, 0, 100);
1688   TH1F* h2_resVsPt_ptRel_I =
1689       new TH1F("resVsPt2_ptRel_I", ";Tracking particle p_{T} [GeV]; p_{T} resolution / p_{T}", 20, 0, 100);
1690   TH1F* h2_resVsPt_ptRel_F =
1691       new TH1F("resVsPt2_ptRel_F", ";Tracking particle p_{T} [GeV]; p_{T} resolution / p_{T}", 20, 0, 100);
1692
1693   TH1F* h2_mresVsPt_pt =
1694       new TH1F("mresVsPt2_pt", ";Tracking particle p_{T} [GeV]; Mean(p_{T} residual) [GeV]", 20, 0, 100);
1695   TH1F* h2_mresVsPt_pt_C =
1696       new TH1F("mresVsPt2_pt_C", ";Tracking particle p_{T} [GeV]; Mean(p_{T} residual) [GeV]", 20, 0, 100);
1697   TH1F* h2_mresVsPt_pt_I =
1698       new TH1F("mresVsPt2_pt_I", ";Tracking particle p_{T} [GeV]; Mean(p_{T} residual) [GeV]", 20, 0, 100);
1699   TH1F* h2_mresVsPt_pt_F =
1700       new TH1F("mresVsPt2_pt_F", ";Tracking particle p_{T} [GeV]; Mean(p_{T} residual) [GeV]", 20, 0, 100);
1701
1702   TH1F* h2_resVsPt_z0 = new TH1F("resVsPt2_z0", ";Tracking particle p_{T} [GeV]; z_{0} resolution [cm]", 20, 0, 100);
1703   TH1F* h2_resVsPt_z0_C =
1704       new TH1F("resVsPt2_z0_C", ";Tracking particle p_{T} [GeV]; z_{0} resolution [cm]", 20, 0, 100);
1705   TH1F* h2_resVsPt_z0_I =
1706       new TH1F("resVsPt2_z0_I", ";Tracking particle p_{T} [GeV]; z_{0} resolution [cm]", 20, 0, 100);
1707   TH1F* h2_resVsPt_z0_F =
1708       new TH1F("resVsPt2_z0_F", ";Tracking particle p_{T} [GeV]; z_{0} resolution [cm]", 20, 0, 100);
1709
1710   TH1F* h2_resVsPt_phi = new TH1F("resVsPt2_phi", ";Tracking particle p_{T} [GeV]; #phi resolution [rad]", 20, 0, 100);
1711   TH1F* h2_resVsPt_phi_C =
1712       new TH1F("resVsPt2_phi_C", ";Tracking particle p_{T} [GeV]; #phi resolution [rad]", 20, 0, 100);
1713   TH1F* h2_resVsPt_phi_I =
1714       new TH1F("resVsPt2_phi_I", ";Tracking particle p_{T} [GeV]; #phi resolution [rad]", 20, 0, 100);
1715   TH1F* h2_resVsPt_phi_F =
1716       new TH1F("resVsPt2_phi_F", ";Tracking particle p_{T} [GeV]; #phi resolution [rad]", 20, 0, 100);
1717
1718   TH1F* h2_resVsPt_eta = new TH1F("resVsPt2_eta", ";Tracking particle p_{T} [GeV]; #eta resolution", 20, 0, 100);
1719
1720   TH1F* h2_resVsPt_d0 = new TH1F("resVsPt2_d0", ";Tracking particle p_{T} [GeV]; d_{0} resolution [cm]", 20, 0, 100);
1721
1722   TH1F* h2_resVsPt_pt_68 =
1723       new TH1F("resVsPt2_pt_68", ";Tracking particle p_{T} [GeV]; p_{T} resolution [GeV]", 20, 0, 100);
1724   TH1F* h2_resVsPt_ptRel_68 =
1725       new TH1F("resVsPt2_ptRel_68", ";Tracking particle p_{T} [GeV]; p_{T} resolution / p_{T}", 20, 0, 100);
1726   TH1F* h2_resVsPt_z0_68 =
1727       new TH1F("resVsPt2_z0_68", ";Tracking particle p_{T} [GeV]; z_{0} resolution [cm]", 20, 0, 100);
1728   TH1F* h2_resVsPt_phi_68 =
1729       new TH1F("resVsPt2_phi_68", ";Tracking particle p_{T} [GeV]; #phi resolution [rad]", 20, 0, 100);
1730   TH1F* h2_resVsPt_eta_68 = new TH1F("resVsPt2_eta_68", ";Tracking particle p_{T} [GeV]; #eta resolution", 20, 0, 100);
1731   TH1F* h2_resVsPt_d0_68 =
1732       new TH1F("resVsPt2_d0_68", ";Tracking particle p_{T} [GeV]; d_{0} resolution [cm]", 20, 0, 100);
1733
1734   TH1F* h2_resVsPt_pt_90 =
1735       new TH1F("resVsPt2_pt_90", ";Tracking particle p_{T} [GeV]; p_{T} resolution [GeV]", 20, 0, 100);
1736   TH1F* h2_resVsPt_ptRel_90 =
1737       new TH1F("resVsPt2_ptRel_90", ";Tracking particle p_{T} [GeV]; p_{T} resolution / p_{T}", 20, 0, 100);
1738   TH1F* h2_resVsPt_z0_90 =
1739       new TH1F("resVsPt2_z0_90", ";Tracking particle p_{T} [GeV]; z_{0} resolution [cm]", 20, 0, 100);
1740   TH1F* h2_resVsPt_phi_90 =
1741       new TH1F("resVsPt2_phi_90", ";Tracking particle p_{T} [GeV]; #phi resolution [rad]", 20, 0, 100);
1742   TH1F* h2_resVsPt_eta_90 = new TH1F("resVsPt2_eta_90", ";Tracking particle p_{T} [GeV]; #eta resolution", 20, 0, 100);
1743   TH1F* h2_resVsPt_d0_90 =
1744       new TH1F("resVsPt2_d0_90", ";Tracking particle p_{T} [GeV]; d_{0} resolution [cm]", 20, 0, 100);
1745
1746   TH1F* h2_resVsPt_pt_99 =
1747       new TH1F("resVsPt2_pt_99", ";Tracking particle p_{T} [GeV]; p_{T} resolution [GeV]", 20, 0, 100);
1748   TH1F* h2_resVsPt_ptRel_99 =
1749       new TH1F("resVsPt2_ptRel_99", ";Tracking particle p_{T} [GeV]; p_{T} resolution / p_{T}", 20, 0, 100);
1750   TH1F* h2_resVsPt_z0_99 =
1751       new TH1F("resVsPt2_z0_99", ";Tracking particle p_{T} [GeV]; z_{0} resolution [cm]", 20, 0, 100);
1752   TH1F* h2_resVsPt_phi_99 =
1753       new TH1F("resVsPt2_phi_99", ";Tracking particle p_{T} [GeV]; #phi resolution [rad]", 20, 0, 100);
1754   TH1F* h2_resVsPt_eta_99 = new TH1F("resVsPt2_eta_99", ";Tracking particle p_{T} [GeV]; #eta resolution", 20, 0, 100);
1755   TH1F* h2_resVsPt_d0_99 =
1756       new TH1F("resVsPt2_d0_99", ";Tracking particle p_{T} [GeV]; d_{0} resolution [cm]", 20, 0, 100);
1757
1758   TH1F* h2_resVsPt_pt_L_68 =
1759       new TH1F("resVsPt2_pt_L_68", ";Tracking particle p_{T} [GeV]; p_{T} resolution [GeV]", nRANGE_L, pt_resmin, 8);
1760   TH1F* h2_resVsPt_ptRel_L_68 = new TH1F(
1761       "resVsPt2_ptRel_L_68", ";Tracking particle p_{T} [GeV]; p_{T} resolution / p_{T}", nRANGE_L, pt_resmin, 8);
1762   TH1F* h2_resVsPt_z0_L_68 =
1763       new TH1F("resVsPt2_z0_L_68", ";Tracking particle p_{T} [GeV]; z_{0} resolution [cm]", nRANGE_L, pt_resmin, 8);
1764   TH1F* h2_resVsPt_phi_L_68 =
1765       new TH1F("resVsPt2_phi_L_68", ";Tracking particle p_{T} [GeV]; #phi resolution [rad]", nRANGE_L, pt_resmin, 8);
1766   TH1F* h2_resVsPt_eta_L_68 =
1767       new TH1F("resVsPt2_eta_L_68", ";Tracking particle p_{T} [GeV]; #eta resolution", nRANGE_L, pt_resmin, 8);
1768   TH1F* h2_resVsPt_d0_L_68 =
1769       new TH1F("resVsPt2_d0_L_68", ";Tracking particle p_{T} [GeV]; d_{0} resolution [cm]", nRANGE_L, pt_resmin, 8);
1770
1771   TH1F* h2_resVsPt_pt_L_90 =
1772       new TH1F("resVsPt2_pt_L_90", ";Tracking particle p_{T} [GeV]; p_{T} resolution [GeV]", nRANGE_L, pt_resmin, 8);
1773   TH1F* h2_resVsPt_ptRel_L_90 = new TH1F(
1774       "resVsPt2_ptRel_L_90", ";Tracking particle p_{T} [GeV]; p_{T} resolution / p_{T}", nRANGE_L, pt_resmin, 8);
1775   TH1F* h2_resVsPt_z0_L_90 =
1776       new TH1F("resVsPt2_z0_L_90", ";Tracking particle p_{T} [GeV]; z_{0} resolution [cm]", nRANGE_L, pt_resmin, 8);
1777   TH1F* h2_resVsPt_phi_L_90 =
1778       new TH1F("resVsPt2_phi_L_90", ";Tracking particle p_{T} [GeV]; #phi resolution [rad]", nRANGE_L, pt_resmin, 8);
1779   TH1F* h2_resVsPt_eta_L_90 =
1780       new TH1F("resVsPt2_eta_L_90", ";Tracking particle p_{T} [GeV]; #eta resolution", nRANGE_L, pt_resmin, 8);
1781   TH1F* h2_resVsPt_d0_L_90 =
1782       new TH1F("resVsPt2_d0_L_90", ";Tracking particle p_{T} [GeV]; d_{0} resolution [cm]", nRANGE_L, pt_resmin, 8);
1783
1784   TH1F* h2_resVsPt_pt_L_99 =
1785       new TH1F("resVsPt2_pt_L_99", ";Tracking particle p_{T} [GeV]; p_{T} resolution [cm]", nRANGE_L, pt_resmin, 8);
1786   TH1F* h2_resVsPt_ptRel_L_99 = new TH1F(
1787       "resVsPt2_ptRel_L_99", ";Tracking particle p_{T} [GeV]; p_{T} resolution / p_{T}", nRANGE_L, pt_resmin, 8);
1788   TH1F* h2_resVsPt_z0_L_99 =
1789       new TH1F("resVsPt2_z0_L_99", ";Tracking particle p_{T} [GeV]; z_{0} resolution [cm]", nRANGE_L, pt_resmin, 8);
1790   TH1F* h2_resVsPt_phi_L_99 =
1791       new TH1F("resVsPt2_phi_L_99", ";Tracking particle p_{T} [GeV]; #phi resolution [rad]", nRANGE_L, pt_resmin, 8);
1792   TH1F* h2_resVsPt_eta_L_99 =
1793       new TH1F("resVsPt2_eta_L_99", ";Tracking particle p_{T} [GeV]; #eta resolution", nRANGE_L, pt_resmin, 8);
1794   TH1F* h2_resVsPt_d0_L_99 =
1795       new TH1F("resVsPt2_d0_L_99", ";Tracking particle p_{T} [GeV]; d_{0} resolution [cm]", nRANGE_L, pt_resmin, 8);
1796
1797   for (int i = 0; i < nRANGE; i++) {
1798     // set bin content and error
1799     h2_resVsPt_pt->SetBinContent(i + 1, h_resVsPt_pt[i]->GetRMS());
1800     h2_resVsPt_pt->SetBinError(i + 1, h_resVsPt_pt[i]->GetRMSError());
1801     h2_resVsPt_pt_C->SetBinContent(i + 1, h_resVsPt_pt_C[i]->GetRMS());
1802     h2_resVsPt_pt_C->SetBinError(i + 1, h_resVsPt_pt_C[i]->GetRMSError());
1803     h2_resVsPt_pt_I->SetBinContent(i + 1, h_resVsPt_pt_I[i]->GetRMS());
1804     h2_resVsPt_pt_I->SetBinError(i + 1, h_resVsPt_pt_I[i]->GetRMSError());
1805     h2_resVsPt_pt_F->SetBinContent(i + 1, h_resVsPt_pt_F[i]->GetRMS());
1806     h2_resVsPt_pt_F->SetBinError(i + 1, h_resVsPt_pt_F[i]->GetRMSError());
1807
1808     h2_resVsPt_ptRel->SetBinContent(i + 1, h_resVsPt_ptRel[i]->GetRMS());
1809     h2_resVsPt_ptRel->SetBinError(i + 1, h_resVsPt_ptRel[i]->GetRMSError());
1810     h2_resVsPt_ptRel_C->SetBinContent(i + 1, h_resVsPt_ptRel_C[i]->GetRMS());
1811     h2_resVsPt_ptRel_C->SetBinError(i + 1, h_resVsPt_ptRel_C[i]->GetRMSError());
1812     h2_resVsPt_ptRel_I->SetBinContent(i + 1, h_resVsPt_ptRel_I[i]->GetRMS());
1813     h2_resVsPt_ptRel_I->SetBinError(i + 1, h_resVsPt_ptRel_I[i]->GetRMSError());
1814     h2_resVsPt_ptRel_F->SetBinContent(i + 1, h_resVsPt_ptRel_F[i]->GetRMS());
1815     h2_resVsPt_ptRel_F->SetBinError(i + 1, h_resVsPt_ptRel_F[i]->GetRMSError());
1816
1817     h2_mresVsPt_pt->SetBinContent(i + 1, h_resVsPt_pt[i]->GetMean());
1818     h2_mresVsPt_pt->SetBinError(i + 1, h_resVsPt_pt[i]->GetMeanError());
1819     h2_mresVsPt_pt_C->SetBinContent(i + 1, h_resVsPt_pt_C[i]->GetMean());
1820     h2_mresVsPt_pt_C->SetBinError(i + 1, h_resVsPt_pt_C[i]->GetMeanError());
1821     h2_mresVsPt_pt_I->SetBinContent(i + 1, h_resVsPt_pt_I[i]->GetMean());
1822     h2_mresVsPt_pt_I->SetBinError(i + 1, h_resVsPt_pt_I[i]->GetMeanError());
1823     h2_mresVsPt_pt_F->SetBinContent(i + 1, h_resVsPt_pt_F[i]->GetMean());
1824     h2_mresVsPt_pt_F->SetBinError(i + 1, h_resVsPt_pt_F[i]->GetMeanError());
1825
1826     h2_resVsPt_z0->SetBinContent(i + 1, h_resVsPt_z0[i]->GetRMS());
1827     h2_resVsPt_z0->SetBinError(i + 1, h_resVsPt_z0[i]->GetRMSError());
1828     h2_resVsPt_z0_C->SetBinContent(i + 1, h_resVsPt_z0_C[i]->GetRMS());
1829     h2_resVsPt_z0_C->SetBinError(i + 1, h_resVsPt_z0_C[i]->GetRMSError());
1830     h2_resVsPt_z0_I->SetBinContent(i + 1, h_resVsPt_z0_I[i]->GetRMS());
1831     h2_resVsPt_z0_I->SetBinError(i + 1, h_resVsPt_z0_I[i]->GetRMSError());
1832     h2_resVsPt_z0_F->SetBinContent(i + 1, h_resVsPt_z0_F[i]->GetRMS());
1833     h2_resVsPt_z0_F->SetBinError(i + 1, h_resVsPt_z0_F[i]->GetRMSError());
1834
1835     h2_resVsPt_phi->SetBinContent(i + 1, h_resVsPt_phi[i]->GetRMS());
1836     h2_resVsPt_phi->SetBinError(i + 1, h_resVsPt_phi[i]->GetRMSError());
1837     h2_resVsPt_phi_C->SetBinContent(i + 1, h_resVsPt_phi_C[i]->GetRMS());
1838     h2_resVsPt_phi_C->SetBinError(i + 1, h_resVsPt_phi_C[i]->GetRMSError());
1839     h2_resVsPt_phi_I->SetBinContent(i + 1, h_resVsPt_phi_I[i]->GetRMS());
1840     h2_resVsPt_phi_I->SetBinError(i + 1, h_resVsPt_phi_I[i]->GetRMSError());
1841     h2_resVsPt_phi_F->SetBinContent(i + 1, h_resVsPt_phi_F[i]->GetRMS());
1842     h2_resVsPt_phi_F->SetBinError(i + 1, h_resVsPt_phi_F[i]->GetRMSError());
1843
1844     h2_resVsPt_eta->SetBinContent(i + 1, h_resVsPt_eta[i]->GetRMS());
1845     h2_resVsPt_eta->SetBinError(i + 1, h_resVsPt_eta[i]->GetRMSError());
1846
1847     h2_resVsPt_d0->SetBinContent(i + 1, h_resVsPt_d0[i]->GetRMS());
1848     h2_resVsPt_d0->SetBinError(i + 1, h_resVsPt_d0[i]->GetRMSError());
1849
1850     h2_resVsPt_pt_68->SetBinContent(i + 1, getIntervalContainingFractionOfEntries(h_absResVsPt_pt[i], 0.68));
1851     h2_resVsPt_pt_90->SetBinContent(i + 1, getIntervalContainingFractionOfEntries(h_absResVsPt_pt[i], 0.90));
1852     h2_resVsPt_pt_99->SetBinContent(i + 1, getIntervalContainingFractionOfEntries(h_absResVsPt_pt[i], 0.99));
1853
1854     h2_resVsPt_ptRel_68->SetBinContent(i + 1, getIntervalContainingFractionOfEntries(h_absResVsPt_ptRel[i], 0.68));
1855     h2_resVsPt_ptRel_90->SetBinContent(i + 1, getIntervalContainingFractionOfEntries(h_absResVsPt_ptRel[i], 0.90));
1856     h2_resVsPt_ptRel_99->SetBinContent(i + 1, getIntervalContainingFractionOfEntries(h_absResVsPt_ptRel[i], 0.99));
1857
1858     h2_resVsPt_eta_68->SetBinContent(i + 1, getIntervalContainingFractionOfEntries(h_absResVsPt_eta[i], 0.68));
1859     h2_resVsPt_eta_90->SetBinContent(i + 1, getIntervalContainingFractionOfEntries(h_absResVsPt_eta[i], 0.90));
1860     h2_resVsPt_eta_99->SetBinContent(i + 1, getIntervalContainingFractionOfEntries(h_absResVsPt_eta[i], 0.99));
1861
1862     h2_resVsPt_z0_68->SetBinContent(i + 1, getIntervalContainingFractionOfEntries(h_absResVsPt_z0[i], 0.68));
1863     h2_resVsPt_z0_90->SetBinContent(i + 1, getIntervalContainingFractionOfEntries(h_absResVsPt_z0[i], 0.90));
1864     h2_resVsPt_z0_99->SetBinContent(i + 1, getIntervalContainingFractionOfEntries(h_absResVsPt_z0[i], 0.99));
1865
1866     h2_resVsPt_phi_68->SetBinContent(i + 1, getIntervalContainingFractionOfEntries(h_absResVsPt_phi[i], 0.68));
1867     h2_resVsPt_phi_90->SetBinContent(i + 1, getIntervalContainingFractionOfEntries(h_absResVsPt_phi[i], 0.90));
1868     h2_resVsPt_phi_99->SetBinContent(i + 1, getIntervalContainingFractionOfEntries(h_absResVsPt_phi[i], 0.99));
1869
1870     h2_resVsPt_d0_68->SetBinContent(i + 1, getIntervalContainingFractionOfEntries(h_absResVsPt_d0[i], 0.68));
1871     h2_resVsPt_d0_90->SetBinContent(i + 1, getIntervalContainingFractionOfEntries(h_absResVsPt_d0[i], 0.90));
1872     h2_resVsPt_d0_99->SetBinContent(i + 1, getIntervalContainingFractionOfEntries(h_absResVsPt_d0[i], 0.99));
1873   }
1874
1875   for (int i = 0; i < nRANGE_L; i++) {
1876     h2_resVsPt_pt_L_68->SetBinContent(i + 1, getIntervalContainingFractionOfEntries(h_absResVsPt_pt_L[i], 0.68));
1877     h2_resVsPt_pt_L_90->SetBinContent(i + 1, getIntervalContainingFractionOfEntries(h_absResVsPt_pt_L[i], 0.90));
1878     h2_resVsPt_pt_L_99->SetBinContent(i + 1, getIntervalContainingFractionOfEntries(h_absResVsPt_pt_L[i], 0.99));
1879
1880     h2_resVsPt_ptRel_L_68->SetBinContent(i + 1, getIntervalContainingFractionOfEntries(h_absResVsPt_ptRel_L[i], 0.68));
1881     h2_resVsPt_ptRel_L_90->SetBinContent(i + 1, getIntervalContainingFractionOfEntries(h_absResVsPt_ptRel_L[i], 0.90));
1882     h2_resVsPt_ptRel_L_99->SetBinContent(i + 1, getIntervalContainingFractionOfEntries(h_absResVsPt_ptRel_L[i], 0.99));
1883
1884     h2_resVsPt_eta_L_68->SetBinContent(i + 1, getIntervalContainingFractionOfEntries(h_absResVsPt_eta_L[i], 0.68));
1885     h2_resVsPt_eta_L_90->SetBinContent(i + 1, getIntervalContainingFractionOfEntries(h_absResVsPt_eta_L[i], 0.90));
1886     h2_resVsPt_eta_L_99->SetBinContent(i + 1, getIntervalContainingFractionOfEntries(h_absResVsPt_eta_L[i], 0.99));
1887
1888     h2_resVsPt_z0_L_68->SetBinContent(i + 1, getIntervalContainingFractionOfEntries(h_absResVsPt_z0_L[i], 0.68));
1889     h2_resVsPt_z0_L_90->SetBinContent(i + 1, getIntervalContainingFractionOfEntries(h_absResVsPt_z0_L[i], 0.90));
1890     h2_resVsPt_z0_L_99->SetBinContent(i + 1, getIntervalContainingFractionOfEntries(h_absResVsPt_z0_L[i], 0.99));
1891
1892     h2_resVsPt_phi_L_68->SetBinContent(i + 1, getIntervalContainingFractionOfEntries(h_absResVsPt_phi_L[i], 0.68));
1893     h2_resVsPt_phi_L_90->SetBinContent(i + 1, getIntervalContainingFractionOfEntries(h_absResVsPt_phi_L[i], 0.90));
1894     h2_resVsPt_phi_L_99->SetBinContent(i + 1, getIntervalContainingFractionOfEntries(h_absResVsPt_phi_L[i], 0.99));
1895
1896     h2_resVsPt_d0_L_68->SetBinContent(i + 1, getIntervalContainingFractionOfEntries(h_absResVsPt_d0_L[i], 0.68));
1897     h2_resVsPt_d0_L_90->SetBinContent(i + 1, getIntervalContainingFractionOfEntries(h_absResVsPt_d0_L[i], 0.90));
1898     h2_resVsPt_d0_L_99->SetBinContent(i + 1, getIntervalContainingFractionOfEntries(h_absResVsPt_d0_L[i], 0.99));
1899   }
1900
1901   // resolution vs. eta histograms
1902   TH1F* h2_resVsEta_eta =
1903       new TH1F("resVsEta_eta", ";Tracking particle |#eta|; #eta resolution", nETARANGE, 0, eta_resmax);
1904   TH1F* h2_resVsEta_eta_L =
1905       new TH1F("resVsEta_eta_L", ";Tracking particle |#eta|; #eta resolution", nETARANGE, 0, eta_resmax);
1906   TH1F* h2_resVsEta_eta_H =
1907       new TH1F("resVsEta_eta_H", ";Tracking particle |#eta|; #eta resolution", nETARANGE, 0, eta_resmax);
1908
1909   TH1F* h2_resVsEta_z0 =
1910       new TH1F("resVsEta_z0", ";Tracking particle |#eta|; z_{0} resolution [cm]", nETARANGE, 0, eta_resmax);
1911   TH1F* h2_resVsEta_z0_L =
1912       new TH1F("resVsEta_z0_L", ";Tracking particle |#eta|; z_{0} resolution [cm]", nETARANGE, 0, eta_resmax);
1913   TH1F* h2_resVsEta_z0_H =
1914       new TH1F("resVsEta_z0_H", ";Tracking particle |#eta|; z_{0} resolution [cm]", nETARANGE, 0, eta_resmax);
1915
1916   TH1F* h2_resVsEta_phi =
1917       new TH1F("resVsEta_phi", ";Tracking particle |#eta|; #phi resolution [rad]", nETARANGE, 0, eta_resmax);
1918   TH1F* h2_resVsEta_phi_L =
1919       new TH1F("resVsEta_phi_L", ";Tracking particle |#eta|; #phi resolution [rad]", nETARANGE, 0, eta_resmax);
1920   TH1F* h2_resVsEta_phi_H =
1921       new TH1F("resVsEta_phi_H", ";Tracking particle |#eta|; #phi resolution [rad]", nETARANGE, 0, eta_resmax);
1922
1923   TH1F* h2_resVsEta_ptRel =
1924       new TH1F("resVsEta_ptRel", ";Tracking particle |#eta|; p_{T} resolution / p_{T}", nETARANGE, 0, eta_resmax);
1925   TH1F* h2_resVsEta_ptRel_L =
1926       new TH1F("resVsEta_ptRel_L", ";Tracking particle |#eta|; p_{T} resolution / p_{T}", nETARANGE, 0, eta_resmax);
1927   TH1F* h2_resVsEta_ptRel_H =
1928       new TH1F("resVsEta_ptRel_H", ";Tracking particle |#eta|; p_{T} resolution / p_{T}", nETARANGE, 0, eta_resmax);
1929
1930   TH1F* h2_resVsEta_d0 =
1931       new TH1F("resVsEta_d0", ";Tracking particle |#eta|; d_{0} resolution [cm]", nETARANGE, 0, eta_resmax);
1932   TH1F* h2_resVsEta_d0_L =
1933       new TH1F("resVsEta_d0_L", ";Tracking particle |#eta|; d_{0} resolution [cm]", nETARANGE, 0, eta_resmax);
1934   TH1F* h2_resVsEta_d0_H =
1935       new TH1F("resVsEta_d0_H", ";Tracking particle |#eta|; d_{0} resolution [cm]", nETARANGE, 0, eta_resmax);
1936
1937   // mean of residuals
1938   TH1F* h2_mresVsEta_eta =
1939       new TH1F("mresVsEta_eta", ";Tracking particle |#eta|; Mean(#eta residual)", nETARANGE, 0, eta_resmax);
1940   TH1F* h2_mresVsEta_z0 =
1941       new TH1F("mresVsEta_z0", ";Tracking particle |#eta|; Mean(z_{0} residual) [cm]", nETARANGE, 0, eta_resmax);
1942   TH1F* h2_mresVsEta_phi =
1943       new TH1F("mresVsEta_phi", ";Tracking particle |#eta|; Mean(phi residual) [rad]", nETARANGE, 0, eta_resmax);
1944   TH1F* h2_mresVsEta_ptRel =
1945       new TH1F("mresVsEta_ptRel", ";Tracking particle |#eta|; Mean(ptrel residual)", nETARANGE, 0, eta_resmax);
1946
1947   // 68 / 90 / 99% residuals
1948   TH1F* h2_resVsEta_eta_68 =
1949       new TH1F("resVsEta_eta_68", ";Tracking particle |#eta|; #eta resolution", nETARANGE, 0, eta_resmax);
1950   TH1F* h2_resVsEta_eta_90 =
1951       new TH1F("resVsEta_eta_90", ";Tracking particle |#eta|; #eta resolution", nETARANGE, 0, eta_resmax);
1952   TH1F* h2_resVsEta_eta_99 =
1953       new TH1F("resVsEta_eta_99", ";Tracking particle |#eta|; #eta resolution", nETARANGE, 0, eta_resmax);
1954
1955   TH1F* h2_resVsEta_z0_68 =
1956       new TH1F("resVsEta_z0_68", ";Tracking particle |#eta|; z_{0} resolution [cm]", nETARANGE, 0, eta_resmax);
1957   TH1F* h2_resVsEta_z0_90 =
1958       new TH1F("resVsEta_z0_90", ";Tracking particle |#eta|; z_{0} resolution [cm]", nETARANGE, 0, eta_resmax);
1959   TH1F* h2_resVsEta_z0_99 =
1960       new TH1F("resVsEta_z0_99", ";Tracking particle |#eta|; z_{0} resolution [cm]", nETARANGE, 0, eta_resmax);
1961
1962   TH1F* h2_resVsEta_phi_68 =
1963       new TH1F("resVsEta_phi_68", ";Tracking particle |#eta|; #phi resolution [rad]", nETARANGE, 0, eta_resmax);
1964   TH1F* h2_resVsEta_phi_90 =
1965       new TH1F("resVsEta_phi_90", ";Tracking particle |#eta|; #phi resolution [rad]", nETARANGE, 0, eta_resmax);
1966   TH1F* h2_resVsEta_phi_99 =
1967       new TH1F("resVsEta_phi_99", ";Tracking particle |#eta|; #phi resolution [rad]", nETARANGE, 0, eta_resmax);
1968
1969   TH1F* h2_resVsEta_ptRel_68 =
1970       new TH1F("resVsEta_ptRel_68", ";Tracking particle |#eta|; p_{T} resolution / p_{T}", nETARANGE, 0, eta_resmax);
1971   TH1F* h2_resVsEta_ptRel_90 =
1972       new TH1F("resVsEta_ptRel_90", ";Tracking particle |#eta|; p_{T} resolution / p_{T}", nETARANGE, 0, eta_resmax);
1973   TH1F* h2_resVsEta_ptRel_99 =
1974       new TH1F("resVsEta_ptRel_99", ";Tracking particle |#eta|; p_{T} resolution / p_{T}", nETARANGE, 0, eta_resmax);
1975
1976   TH1F* h2_resVsEta_d0_68 =
1977       new TH1F("resVsEta_d0_68", ";Tracking particle |#eta|; d_{0} resolution [cm]", nETARANGE, 0, eta_resmax);
1978   TH1F* h2_resVsEta_d0_90 =
1979       new TH1F("resVsEta_d0_90", ";Tracking particle |#eta|; d_{0} resolution [cm]", nETARANGE, 0, eta_resmax);
1980   TH1F* h2_resVsEta_d0_99 =
1981       new TH1F("resVsEta_d0_99", ";Tracking particle |#eta|; d_{0} resolution [cm]", nETARANGE, 0, eta_resmax);
1982
1983   TH1F* h2_resVsEta_eta_L_68 =
1984       new TH1F("resVsEta_eta_L_68", ";Tracking particle |#eta|; #eta resolution", nETARANGE, 0, eta_resmax);
1985   TH1F* h2_resVsEta_z0_L_68 =
1986       new TH1F("resVsEta_z0_L_68", ";Tracking particle |#eta|; z_{0} resolution [cm]", nETARANGE, 0, eta_resmax);
1987   TH1F* h2_resVsEta_phi_L_68 =
1988       new TH1F("resVsEta_phi_L_68", ";Tracking particle |#eta|; #phi resolution [rad]", nETARANGE, 0, eta_resmax);
1989   TH1F* h2_resVsEta_ptRel_L_68 =
1990       new TH1F("resVsEta_ptRel_L_68", ";Tracking particle |#eta|; p_{T} resolution / p_{T}", nETARANGE, 0, eta_resmax);
1991   TH1F* h2_resVsEta_d0_L_68 =
1992       new TH1F("resVsEta_d0_L_68", ";Tracking particle |#eta|; d_{0} resolution [cm]", nETARANGE, 0, eta_resmax);
1993
1994   TH1F* h2_resVsEta_eta_L_90 =
1995       new TH1F("resVsEta_eta_L_90", ";Tracking particle |#eta|; #eta resolution", nETARANGE, 0, eta_resmax);
1996   TH1F* h2_resVsEta_z0_L_90 =
1997       new TH1F("resVsEta_z0_L_90", ";Tracking particle |#eta|; z_{0} resolution [cm]", nETARANGE, 0, eta_resmax);
1998   TH1F* h2_resVsEta_phi_L_90 =
1999       new TH1F("resVsEta_phi_L_90", ";Tracking particle |#eta|; #phi resolution [rad]", nETARANGE, 0, eta_resmax);
2000   TH1F* h2_resVsEta_ptRel_L_90 =
2001       new TH1F("resVsEta_ptRel_L_90", ";Tracking particle |#eta|; p_{T} resolution / p_{T}", nETARANGE, 0, eta_resmax);
2002   TH1F* h2_resVsEta_d0_L_90 =
2003       new TH1F("resVsEta_d0_L_90", ";Tracking particle |#eta|; d_{0} resolution [cm]", nETARANGE, 0, eta_resmax);
2004
2005   TH1F* h2_resVsEta_eta_L_99 =
2006       new TH1F("resVsEta_eta_L_99", ";Tracking particle |#eta|; #eta resolution", nETARANGE, 0, eta_resmax);
2007   TH1F* h2_resVsEta_z0_L_99 =
2008       new TH1F("resVsEta_z0_L_99", ";Tracking particle |#eta|; z_{0} resolution [cm]", nETARANGE, 0, eta_resmax);
2009   TH1F* h2_resVsEta_phi_L_99 =
2010       new TH1F("resVsEta_phi_L_99", ";Tracking particle |#eta|; #phi resolution [rad]", nETARANGE, 0, eta_resmax);
2011   TH1F* h2_resVsEta_ptRel_L_99 =
2012       new TH1F("resVsEta_ptRel_L_99", ";Tracking particle |#eta|; p_{T} resolution / p_{T}", nETARANGE, 0, eta_resmax);
2013   TH1F* h2_resVsEta_d0_L_99 =
2014       new TH1F("resVsEta_d0_L_99", ";Tracking particle |#eta|; d_{0} resolution [cm]", nETARANGE, 0, eta_resmax);
2015
2016   TH1F* h2_resVsEta_eta_H_68 =
2017       new TH1F("resVsEta_eta_H_68", ";Tracking particle |#eta|; #eta resolution", nETARANGE, 0, eta_resmax);
2018   TH1F* h2_resVsEta_z0_H_68 =
2019       new TH1F("resVsEta_z0_H_68", ";Tracking particle |#eta|; z_{0} resolution [cm]", nETARANGE, 0, eta_resmax);
2020   TH1F* h2_resVsEta_phi_H_68 =
2021       new TH1F("resVsEta_phi_H_68", ";Tracking particle |#eta|; #phi resolution [rad]", nETARANGE, 0, eta_resmax);
2022   TH1F* h2_resVsEta_ptRel_H_68 =
2023       new TH1F("resVsEta_ptRel_H_68", ";Tracking particle |#eta|; p_{T} resolution / p_{T}", nETARANGE, 0, eta_resmax);
2024   TH1F* h2_resVsEta_d0_H_68 =
2025       new TH1F("resVsEta_d0_H_68", ";Tracking particle |#eta|; d_{0} resolution [cm]", nETARANGE, 0, eta_resmax);
2026
2027   TH1F* h2_resVsEta_eta_H_90 =
2028       new TH1F("resVsEta_eta_H_90", ";Tracking particle |#eta|; #eta resolution", nETARANGE, 0, eta_resmax);
2029   TH1F* h2_resVsEta_z0_H_90 =
2030       new TH1F("resVsEta_z0_H_90", ";Tracking particle |#eta|; z_{0} resolution [cm]", nETARANGE, 0, eta_resmax);
2031   TH1F* h2_resVsEta_phi_H_90 =
2032       new TH1F("resVsEta_phi_H_90", ";Tracking particle |#eta|; #phi resolution [rad]", nETARANGE, 0, eta_resmax);
2033   TH1F* h2_resVsEta_ptRel_H_90 =
2034       new TH1F("resVsEta_ptRel_H_90", ";Tracking particle |#eta|; p_{T} resolution / p_{T}", nETARANGE, 0, eta_resmax);
2035   TH1F* h2_resVsEta_d0_H_90 =
2036       new TH1F("resVsEta_d0_H_90", ";Tracking particle |#eta|; d_{0} resolution [cm]", nETARANGE, 0, eta_resmax);
2037
2038   TH1F* h2_resVsEta_eta_H_99 =
2039       new TH1F("resVsEta_eta_H_99", ";Tracking particle |#eta|; #eta resolution", nETARANGE, 0, eta_resmax);
2040   TH1F* h2_resVsEta_z0_H_99 =
2041       new TH1F("resVsEta_z0_H_99", ";Tracking particle |#eta|; z_{0} resolution [cm]", nETARANGE, 0, eta_resmax);
2042   TH1F* h2_resVsEta_phi_H_99 =
2043       new TH1F("resVsEta_phi_H_99", ";Tracking particle |#eta|; #phi resolution [rad]", nETARANGE, 0, eta_resmax);
2044   TH1F* h2_resVsEta_ptRel_H_99 =
2045       new TH1F("resVsEta_ptRel_H_99", ";Tracking particle |#eta|; p_{T} resolution / p_{T}", nETARANGE, 0, eta_resmax);
2046   TH1F* h2_resVsEta_d0_H_99 =
2047       new TH1F("resVsEta_d0_H_99", ";Tracking particle |#eta|; d_{0} resolution [cm]", nETARANGE, 0, eta_resmax);
2048
2049   // resolution vs. eta histograms (gaussian fit)
2050   TH1F* h3_resVsEta_eta_L = new TH1F("resVsEta_eta_L_gaus", ";|#eta|; #sigma(#eta)", nETARANGE, 0, eta_resmax);
2051   TH1F* h3_resVsEta_eta_H = new TH1F("resVsEta_eta_H_gaus", ";|#eta|; #sigma(#eta)", nETARANGE, 0, eta_resmax);
2052
2053   TH1F* h3_resVsEta_z0_L = new TH1F("resVsEta_z0_L_gaus", ";|#eta|; #sigma(z_{0}) [cm]", nETARANGE, 0, eta_resmax);
2054   TH1F* h3_resVsEta_z0_H = new TH1F("resVsEta_z0_H_gaus", ";|#eta|; #sigma(z_{0}) [cm]", nETARANGE, 0, eta_resmax);
2055
2056   TH1F* h3_resVsEta_phi_L = new TH1F("resVsEta_phi_L_gaus", ";|#eta|; #sigma(#phi) [rad]", nETARANGE, 0, eta_resmax);
2057   TH1F* h3_resVsEta_phi_H = new TH1F("resVsEta_phi_H_gaus", ";|#eta|; #sigma(#phi) [rad]", nETARANGE, 0, eta_resmax);
2058
2059   TH1F* h3_resVsEta_ptRel_L =
2060       new TH1F("resVsEta_ptRel_L_gaus", ";|#eta|; #sigma(p_{T}) / p_{T}", nETARANGE, 0, eta_resmax);
2061   TH1F* h3_resVsEta_ptRel_H =
2062       new TH1F("resVsEta_ptRel_H_gaus", ";|#eta|; #sigma(p_{T}) / p_{T}", nETARANGE, 0, eta_resmax);
2063
2064   gSystem->mkdir("FitResults");
2065   TString fitdir = "FitResults/";
2066
2067   for (int i = 0; i < nETARANGE; i++) {
2068     // set bin content and error
2069     h2_resVsEta_eta->SetBinContent(i + 1, h_resVsEta_eta[i]->GetRMS());
2070     h2_resVsEta_eta->SetBinError(i + 1, h_resVsEta_eta[i]->GetRMSError());
2071     h2_resVsEta_eta_L->SetBinContent(i + 1, h_resVsEta_eta_L[i]->GetRMS());
2072     h2_resVsEta_eta_L->SetBinError(i + 1, h_resVsEta_eta_L[i]->GetRMSError());
2073     h2_resVsEta_eta_H->SetBinContent(i + 1, h_resVsEta_eta_H[i]->GetRMS());
2074     h2_resVsEta_eta_H->SetBinError(i + 1, h_resVsEta_eta_H[i]->GetRMSError());
2075
2076     h2_resVsEta_z0->SetBinContent(i + 1, h_resVsEta_z0[i]->GetRMS());
2077     h2_resVsEta_z0->SetBinError(i + 1, h_resVsEta_z0[i]->GetRMSError());
2078     h2_resVsEta_z0_L->SetBinContent(i + 1, h_resVsEta_z0_L[i]->GetRMS());
2079     h2_resVsEta_z0_L->SetBinError(i + 1, h_resVsEta_z0_L[i]->GetRMSError());
2080     h2_resVsEta_z0_H->SetBinContent(i + 1, h_resVsEta_z0_H[i]->GetRMS());
2081     h2_resVsEta_z0_H->SetBinError(i + 1, h_resVsEta_z0_H[i]->GetRMSError());
2082
2083     h2_resVsEta_phi->SetBinContent(i + 1, h_resVsEta_phi[i]->GetRMS());
2084     h2_resVsEta_phi->SetBinError(i + 1, h_resVsEta_phi[i]->GetRMSError());
2085     h2_resVsEta_phi_L->SetBinContent(i + 1, h_resVsEta_phi_L[i]->GetRMS());
2086     h2_resVsEta_phi_L->SetBinError(i + 1, h_resVsEta_phi_L[i]->GetRMSError());
2087     h2_resVsEta_phi_H->SetBinContent(i + 1, h_resVsEta_phi_H[i]->GetRMS());
2088     h2_resVsEta_phi_H->SetBinError(i + 1, h_resVsEta_phi_H[i]->GetRMSError());
2089
2090     h2_resVsEta_ptRel->SetBinContent(i + 1, h_resVsEta_ptRel[i]->GetRMS());
2091     h2_resVsEta_ptRel->SetBinError(i + 1, h_resVsEta_ptRel[i]->GetRMSError());
2092     h2_resVsEta_ptRel_L->SetBinContent(i + 1, h_resVsEta_ptRel_L[i]->GetRMS());
2093     h2_resVsEta_ptRel_L->SetBinError(i + 1, h_resVsEta_ptRel_L[i]->GetRMSError());
2094     h2_resVsEta_ptRel_H->SetBinContent(i + 1, h_resVsEta_ptRel_H[i]->GetRMS());
2095     h2_resVsEta_ptRel_H->SetBinError(i + 1, h_resVsEta_ptRel_H[i]->GetRMSError());
2096
2097     h2_mresVsEta_eta->SetBinContent(i + 1, h_resVsEta_eta[i]->GetMean());
2098     h2_mresVsEta_eta->SetBinError(i + 1, h_resVsEta_eta[i]->GetMeanError());
2099     h2_mresVsEta_z0->SetBinContent(i + 1, h_resVsEta_z0[i]->GetMean());
2100     h2_mresVsEta_z0->SetBinError(i + 1, h_resVsEta_z0[i]->GetMeanError());
2101     h2_mresVsEta_phi->SetBinContent(i + 1, h_resVsEta_phi[i]->GetMean());
2102     h2_mresVsEta_phi->SetBinError(i + 1, h_resVsEta_phi[i]->GetMeanError());
2103     h2_mresVsEta_ptRel->SetBinContent(i + 1, h_resVsEta_ptRel[i]->GetMean());
2104     h2_mresVsEta_ptRel->SetBinError(i + 1, h_resVsEta_ptRel[i]->GetMeanError());
2105
2106     h2_resVsEta_d0->SetBinContent(i + 1, h_resVsEta_d0[i]->GetRMS());
2107     h2_resVsEta_d0->SetBinError(i + 1, h_resVsEta_d0[i]->GetRMSError());
2108     h2_resVsEta_d0_L->SetBinContent(i + 1, h_resVsEta_d0_L[i]->GetRMS());
2109     h2_resVsEta_d0_L->SetBinError(i + 1, h_resVsEta_d0_L[i]->GetRMSError());
2110     h2_resVsEta_d0_H->SetBinContent(i + 1, h_resVsEta_d0_H[i]->GetRMS());
2111     h2_resVsEta_d0_H->SetBinError(i + 1, h_resVsEta_d0_H[i]->GetRMSError());
2112
2113     h2_resVsEta_eta_68->SetBinContent(i + 1, getIntervalContainingFractionOfEntries(h_absResVsEta_eta[i], 0.68));
2114     h2_resVsEta_eta_90->SetBinContent(i + 1, getIntervalContainingFractionOfEntries(h_absResVsEta_eta[i], 0.90));
2115     h2_resVsEta_eta_99->SetBinContent(i + 1, getIntervalContainingFractionOfEntries(h_absResVsEta_eta[i], 0.99));
2116
2117     h2_resVsEta_z0_68->SetBinContent(i + 1, getIntervalContainingFractionOfEntries(h_absResVsEta_z0[i], 0.68));
2118     h2_resVsEta_z0_90->SetBinContent(i + 1, getIntervalContainingFractionOfEntries(h_absResVsEta_z0[i], 0.90));
2119     h2_resVsEta_z0_99->SetBinContent(i + 1, getIntervalContainingFractionOfEntries(h_absResVsEta_z0[i], 0.99));
2120
2121     h2_resVsEta_phi_68->SetBinContent(i + 1, getIntervalContainingFractionOfEntries(h_absResVsEta_phi[i], 0.68));
2122     h2_resVsEta_phi_90->SetBinContent(i + 1, getIntervalContainingFractionOfEntries(h_absResVsEta_phi[i], 0.90));
2123     h2_resVsEta_phi_99->SetBinContent(i + 1, getIntervalContainingFractionOfEntries(h_absResVsEta_phi[i], 0.99));
2124
2125     h2_resVsEta_ptRel_68->SetBinContent(i + 1, getIntervalContainingFractionOfEntries(h_absResVsEta_ptRel[i], 0.68));
2126     h2_resVsEta_ptRel_90->SetBinContent(i + 1, getIntervalContainingFractionOfEntries(h_absResVsEta_ptRel[i], 0.90));
2127     h2_resVsEta_ptRel_99->SetBinContent(i + 1, getIntervalContainingFractionOfEntries(h_absResVsEta_ptRel[i], 0.99));
2128
2129     h2_resVsEta_d0_68->SetBinContent(i + 1, getIntervalContainingFractionOfEntries(h_absResVsEta_d0[i], 0.68));
2130     h2_resVsEta_d0_90->SetBinContent(i + 1, getIntervalContainingFractionOfEntries(h_absResVsEta_d0[i], 0.90));
2131     h2_resVsEta_d0_99->SetBinContent(i + 1, getIntervalContainingFractionOfEntries(h_absResVsEta_d0[i], 0.99));
2132
2133     h2_resVsEta_eta_L_68->SetBinContent(i + 1, getIntervalContainingFractionOfEntries(h_absResVsEta_eta_L[i], 0.68));
2134     h2_resVsEta_z0_L_68->SetBinContent(i + 1, getIntervalContainingFractionOfEntries(h_absResVsEta_z0_L[i], 0.68));
2135     h2_resVsEta_phi_L_68->SetBinContent(i + 1, getIntervalContainingFractionOfEntries(h_absResVsEta_phi_L[i], 0.68));
2136     h2_resVsEta_ptRel_L_68->SetBinContent(i + 1,
2137                                           getIntervalContainingFractionOfEntries(h_absResVsEta_ptRel_L[i], 0.68));
2138     h2_resVsEta_d0_L_68->SetBinContent(i + 1, getIntervalContainingFractionOfEntries(h_absResVsEta_d0_L[i], 0.68));
2139
2140     h2_resVsEta_eta_L_90->SetBinContent(i + 1, getIntervalContainingFractionOfEntries(h_absResVsEta_eta_L[i], 0.90));
2141     h2_resVsEta_z0_L_90->SetBinContent(i + 1, getIntervalContainingFractionOfEntries(h_absResVsEta_z0_L[i], 0.90));
2142     h2_resVsEta_phi_L_90->SetBinContent(i + 1, getIntervalContainingFractionOfEntries(h_absResVsEta_phi_L[i], 0.90));
2143     h2_resVsEta_ptRel_L_90->SetBinContent(i + 1,
2144                                           getIntervalContainingFractionOfEntries(h_absResVsEta_ptRel_L[i], 0.90));
2145     h2_resVsEta_d0_L_90->SetBinContent(i + 1, getIntervalContainingFractionOfEntries(h_absResVsEta_d0_L[i], 0.90));
2146
2147     h2_resVsEta_eta_L_99->SetBinContent(i + 1, getIntervalContainingFractionOfEntries(h_absResVsEta_eta_L[i], 0.99));
2148     h2_resVsEta_z0_L_99->SetBinContent(i + 1, getIntervalContainingFractionOfEntries(h_absResVsEta_z0_L[i], 0.99));
2149     h2_resVsEta_phi_L_99->SetBinContent(i + 1, getIntervalContainingFractionOfEntries(h_absResVsEta_phi_L[i], 0.99));
2150     h2_resVsEta_ptRel_L_99->SetBinContent(i + 1,
2151                                           getIntervalContainingFractionOfEntries(h_absResVsEta_ptRel_L[i], 0.99));
2152     h2_resVsEta_d0_L_99->SetBinContent(i + 1, getIntervalContainingFractionOfEntries(h_absResVsEta_d0_L[i], 0.99));
2153
2154     h2_resVsEta_eta_H_68->SetBinContent(i + 1, getIntervalContainingFractionOfEntries(h_absResVsEta_eta_H[i], 0.68));
2155     h2_resVsEta_z0_H_68->SetBinContent(i + 1, getIntervalContainingFractionOfEntries(h_absResVsEta_z0_H[i], 0.68));
2156     h2_resVsEta_phi_H_68->SetBinContent(i + 1, getIntervalContainingFractionOfEntries(h_absResVsEta_phi_H[i], 0.68));
2157     h2_resVsEta_ptRel_H_68->SetBinContent(i + 1,
2158                                           getIntervalContainingFractionOfEntries(h_absResVsEta_ptRel_H[i], 0.68));
2159     h2_resVsEta_d0_H_68->SetBinContent(i + 1, getIntervalContainingFractionOfEntries(h_absResVsEta_d0_H[i], 0.68));
2160
2161     h2_resVsEta_eta_H_90->SetBinContent(i + 1, getIntervalContainingFractionOfEntries(h_absResVsEta_eta_H[i], 0.90));
2162     h2_resVsEta_z0_H_90->SetBinContent(i + 1, getIntervalContainingFractionOfEntries(h_absResVsEta_z0_H[i], 0.90));
2163     h2_resVsEta_phi_H_90->SetBinContent(i + 1, getIntervalContainingFractionOfEntries(h_absResVsEta_phi_H[i], 0.90));
2164     h2_resVsEta_ptRel_H_90->SetBinContent(i + 1,
2165                                           getIntervalContainingFractionOfEntries(h_absResVsEta_ptRel_H[i], 0.90));
2166     h2_resVsEta_d0_H_90->SetBinContent(i + 1, getIntervalContainingFractionOfEntries(h_absResVsEta_d0_H[i], 0.90));
2167
2168     h2_resVsEta_eta_H_99->SetBinContent(i + 1, getIntervalContainingFractionOfEntries(h_absResVsEta_eta_H[i], 0.99));
2169     h2_resVsEta_z0_H_99->SetBinContent(i + 1, getIntervalContainingFractionOfEntries(h_absResVsEta_z0_H[i], 0.99));
2170     h2_resVsEta_phi_H_99->SetBinContent(i + 1, getIntervalContainingFractionOfEntries(h_absResVsEta_phi_H[i], 0.99));
2171     h2_resVsEta_ptRel_H_99->SetBinContent(i + 1,
2172                                           getIntervalContainingFractionOfEntries(h_absResVsEta_ptRel_H[i], 0.99));
2173     h2_resVsEta_d0_H_99->SetBinContent(i + 1, getIntervalContainingFractionOfEntries(h_absResVsEta_d0_H[i], 0.99));
2174
2175     // ---------------------------------------------------------------------------------------------------
2176     // gaussian fit instead
2177     // ---------------------------------------------------------------------------------------------------
2178
2179     if (doGausFit) {
2180       TCanvas cfit;
2181       char text[500];
2182
2183       float sigma = 0;
2184       float esigma = 0;
2185       TF1* fit;
2186
2187       float rms = 0;
2188       float erms = 0;
2189
2190       fit = new TF1("fit", "gaus", -0.01, 0.01);
2191       h_resVsEta_eta_L[i]->Fit("fit", "R");
2192       sigma = fit->GetParameter(2);
2193       esigma = fit->GetParError(2);
2194       rms = h_resVsEta_eta_L[i]->GetRMS();
2195       erms = h_resVsEta_eta_L[i]->GetRMSError();
2196       h3_resVsEta_eta_L->SetBinContent(i + 1, sigma);
2197       h3_resVsEta_eta_L->SetBinError(i + 1, esigma);
2198       h_resVsEta_eta_L[i]->Draw();
2199       sprintf(text, "RMS: %.4f +/- %.4f", rms, erms);
2200       mySmallText(0.2, 0.86, 1, text);
2201       sprintf(text, "Fit: %.4f +/- %.4f", sigma, esigma);
2202       mySmallText(0.2, 0.8, 2, text);
2203       sprintf(text, "p_{T} < 5 GeV");
2204       mySmallText(0.2, 0.7, 1, text);
2205       cfit.SaveAs(fitdir + "resVsEta_eta_L_" + etarange[i] + ".pdf");
2206       delete fit;
2207
2208       fit = new TF1("fit", "gaus", -0.01, 0.01);
2209       h_resVsEta_eta_H[i]->Fit("fit", "R");
2210       sigma = fit->GetParameter(2);
2211       esigma = fit->GetParError(2);
2212       rms = h_resVsEta_eta_H[i]->GetRMS();
2213       erms = h_resVsEta_eta_H[i]->GetRMSError();
2214       h3_resVsEta_eta_H->SetBinContent(i + 1, sigma);
2215       h3_resVsEta_eta_H->SetBinError(i + 1, esigma);
2216       h_resVsEta_eta_H[i]->Draw();
2217       sprintf(text, "RMS: %.4f +/- %.4f", rms, erms);
2218       mySmallText(0.2, 0.86, 1, text);
2219       sprintf(text, "Fit: %.4f +/- %.4f", sigma, esigma);
2220       mySmallText(0.2, 0.8, 2, text);
2221       sprintf(text, "p_{T} > 15 GeV");
2222       mySmallText(0.2, 0.7, 1, text);
2223       cfit.SaveAs(fitdir + "resVsEta_eta_H_" + etarange[i] + ".pdf");
2224       delete fit;
2225
2226       fit = new TF1("fit", "gaus", -1, 1);
2227       h_resVsEta_z0_L[i]->Fit("fit", "R");
2228       sigma = fit->GetParameter(2);
2229       esigma = fit->GetParError(2);
2230       rms = h_resVsEta_z0_L[i]->GetRMS();
2231       erms = h_resVsEta_z0_L[i]->GetRMSError();
2232       h3_resVsEta_z0_L->SetBinContent(i + 1, sigma);
2233       h3_resVsEta_z0_L->SetBinError(i + 1, esigma);
2234       h_resVsEta_z0_L[i]->Draw();
2235       sprintf(text, "RMS: %.4f +/- %.4f", rms, erms);
2236       mySmallText(0.2, 0.86, 1, text);
2237       sprintf(text, "Fit: %.4f +/- %.4f", sigma, esigma);
2238       mySmallText(0.2, 0.8, 2, text);
2239       sprintf(text, "p_{T} < 5 GeV");
2240       mySmallText(0.2, 0.7, 1, text);
2241       cfit.SaveAs(fitdir + "resVsEta_z0_L_" + etarange[i] + ".pdf");
2242       delete fit;
2243
2244       fit = new TF1("fit", "gaus", -1, 1);
2245       h_resVsEta_z0_H[i]->Fit("fit", "R");
2246       sigma = fit->GetParameter(2);
2247       esigma = fit->GetParError(2);
2248       rms = h_resVsEta_z0_H[i]->GetRMS();
2249       erms = h_resVsEta_z0_H[i]->GetRMSError();
2250       h3_resVsEta_z0_H->SetBinContent(i + 1, sigma);
2251       h3_resVsEta_z0_H->SetBinError(i + 1, esigma);
2252       h_resVsEta_z0_H[i]->Draw();
2253       sprintf(text, "RMS: %.4f +/- %.4f", rms, erms);
2254       mySmallText(0.2, 0.86, 1, text);
2255       sprintf(text, "Fit: %.4f +/- %.4f", sigma, esigma);
2256       mySmallText(0.2, 0.8, 2, text);
2257       sprintf(text, "p_{T} > 15 GeV");
2258       mySmallText(0.2, 0.7, 1, text);
2259       cfit.SaveAs(fitdir + "resVsEta_z0_H_" + etarange[i] + ".pdf");
2260       delete fit;
2261
2262       fit = new TF1("fit", "gaus", -0.005, 0.005);
2263       h_resVsEta_phi_L[i]->Fit("fit", "R");
2264       sigma = fit->GetParameter(2);
2265       esigma = fit->GetParError(2);
2266       rms = h_resVsEta_phi_L[i]->GetRMS();
2267       erms = h_resVsEta_phi_L[i]->GetRMSError();
2268       h3_resVsEta_phi_L->SetBinContent(i + 1, sigma);
2269       h3_resVsEta_phi_L->SetBinError(i + 1, esigma);
2270       h_resVsEta_phi_L[i]->Draw();
2271       sprintf(text, "RMS: %.4f +/- %.4f", rms, erms);
2272       mySmallText(0.2, 0.86, 1, text);
2273       sprintf(text, "Fit: %.4f +/- %.4f", sigma, esigma);
2274       mySmallText(0.2, 0.8, 2, text);
2275       sprintf(text, "p_{T} < 5 GeV");
2276       mySmallText(0.2, 0.7, 1, text);
2277       cfit.SaveAs(fitdir + "resVsEta_phi_L_" + etarange[i] + ".pdf");
2278       delete fit;
2279
2280       fit = new TF1("fit", "gaus", -0.005, 0.005);
2281       h_resVsEta_phi_H[i]->Fit("fit", "R");
2282       sigma = fit->GetParameter(2);
2283       esigma = fit->GetParError(2);
2284       rms = h_resVsEta_phi_H[i]->GetRMS();
2285       erms = h_resVsEta_phi_H[i]->GetRMSError();
2286       h3_resVsEta_phi_H->SetBinContent(i + 1, sigma);
2287       h3_resVsEta_phi_H->SetBinError(i + 1, esigma);
2288       h_resVsEta_phi_H[i]->Draw();
2289       sprintf(text, "RMS: %.4f +/- %.4f", rms, erms);
2290       mySmallText(0.2, 0.86, 1, text);
2291       sprintf(text, "Fit: %.4f +/- %.4f", sigma, esigma);
2292       mySmallText(0.2, 0.8, 2, text);
2293       sprintf(text, "p_{T} > 15 GeV");
2294       mySmallText(0.2, 0.7, 1, text);
2295       cfit.SaveAs(fitdir + "resVsEta_phi_H_" + etarange[i] + ".pdf");
2296       delete fit;
2297
2298       fit = new TF1("fit", "gaus", -0.5, 0.5);
2299       h_resVsEta_ptRel_L[i]->Fit("fit", "R");
2300       sigma = fit->GetParameter(2);
2301       esigma = fit->GetParError(2);
2302       rms = h_resVsEta_ptRel_L[i]->GetRMS();
2303       erms = h_resVsEta_ptRel_L[i]->GetRMSError();
2304       h3_resVsEta_ptRel_L->SetBinContent(i + 1, sigma);
2305       h3_resVsEta_ptRel_L->SetBinError(i + 1, esigma);
2306       h_resVsEta_ptRel_L[i]->Draw();
2307       sprintf(text, "RMS: %.4f +/- %.4f", rms, erms);
2308       mySmallText(0.2, 0.86, 1, text);
2309       sprintf(text, "Fit: %.4f +/- %.4f", sigma, esigma);
2310       mySmallText(0.2, 0.8, 2, text);
2311       sprintf(text, "p_{T} < 5 GeV");
2312       mySmallText(0.2, 0.7, 1, text);
2313       cfit.SaveAs(fitdir + "resVsEta_ptRel_L_" + etarange[i] + ".pdf");
2314       delete fit;
2315
2316       fit = new TF1("fit", "gaus", -0.5, 0.5);
2317       h_resVsEta_ptRel_H[i]->Fit("fit", "R");
2318       sigma = fit->GetParameter(2);
2319       esigma = fit->GetParError(2);
2320       rms = h_resVsEta_ptRel_H[i]->GetRMS();
2321       erms = h_resVsEta_ptRel_H[i]->GetRMSError();
2322       h3_resVsEta_ptRel_H->SetBinContent(i + 1, sigma);
2323       h3_resVsEta_ptRel_H->SetBinError(i + 1, esigma);
2324       h_resVsEta_ptRel_H[i]->Draw();
2325       sprintf(text, "RMS: %.4f +/- %.4f", rms, erms);
2326       mySmallText(0.2, 0.86, 1, text);
2327       sprintf(text, "Fit: %.4f +/- %.4f", sigma, esigma);
2328       mySmallText(0.2, 0.8, 2, text);
2329       sprintf(text, "p_{T} > 15 GeV");
2330       mySmallText(0.2, 0.7, 1, text);
2331       cfit.SaveAs(fitdir + "resVsEta_ptRel_H_" + etarange[i] + ".pdf");
2332       delete fit;
2333
2334     }  //end doGausFit
2335   }
2336
2337   TH1F* h2_resVsPhi_pt_68 =
2338       new TH1F("resVsPhi2_pt_68", ";Tracking particle #phi; p_{T} resolution [GeV]", 32, -3.2, 3.2);
2339   TH1F* h2_resVsPhi_pt_90 =
2340       new TH1F("resVsPhi2_pt_90", ";Tracking particle #phi; p_{T} resolution [GeV]", 32, -3.2, 3.2);
2341   TH1F* h2_resVsPhi_pt_99 =
2342       new TH1F("resVsPhi2_pt_99", ";Tracking particle #phi; p_{T} resolution [GeV]", 32, -3.2, 3.2);
2343   TH1F* h2_resVsPhi_ptRel_68 =
2344       new TH1F("resVsPhi2_ptRel_68", ";Tracking particle #phi; p_{T} resolution / p_{T}", 32, -3.2, 3.2);
2345   TH1F* h2_resVsPhi_ptRel_90 =
2346       new TH1F("resVsPhi2_ptRel_90", ";Tracking particle #phi; p_{T} resolution / p_{T}", 32, -3.2, 3.2);
2347   TH1F* h2_resVsPhi_ptRel_99 =
2348       new TH1F("resVsPhi2_ptRel_99", ";Tracking particle #phi; p_{T} resolution / p_{T}", 32, -3.2, 3.2);
2349
2350   for (int i = 0; i < nPHIRANGE; i++) {
2351     h2_resVsPhi_pt_68->SetBinContent(i + 1, getIntervalContainingFractionOfEntries(h_absResVsPhi_pt[i], 0.68));
2352     h2_resVsPhi_pt_90->SetBinContent(i + 1, getIntervalContainingFractionOfEntries(h_absResVsPhi_pt[i], 0.90));
2353     h2_resVsPhi_pt_99->SetBinContent(i + 1, getIntervalContainingFractionOfEntries(h_absResVsPhi_pt[i], 0.99));
2354
2355     h2_resVsPhi_ptRel_68->SetBinContent(i + 1, getIntervalContainingFractionOfEntries(h_absResVsPhi_ptRel[i], 0.68));
2356     h2_resVsPhi_ptRel_90->SetBinContent(i + 1, getIntervalContainingFractionOfEntries(h_absResVsPhi_ptRel[i], 0.90));
2357     h2_resVsPhi_ptRel_99->SetBinContent(i + 1, getIntervalContainingFractionOfEntries(h_absResVsPhi_ptRel[i], 0.99));
2358   }
2359
2360   // set minimum to zero
2361   h2_resVsPt_pt->SetMinimum(0);
2362   h2_resVsPt_pt_C->SetMinimum(0);
2363   h2_resVsPt_pt_I->SetMinimum(0);
2364   h2_resVsPt_pt_F->SetMinimum(0);
2365
2366   h2_resVsPt_ptRel->SetMinimum(0);
2367   h2_resVsPt_ptRel_C->SetMinimum(0);
2368   h2_resVsPt_ptRel_I->SetMinimum(0);
2369   h2_resVsPt_ptRel_F->SetMinimum(0);
2370
2371   h2_resVsPt_z0->SetMinimum(0);
2372   h2_resVsPt_z0_C->SetMinimum(0);
2373   h2_resVsPt_z0_I->SetMinimum(0);
2374   h2_resVsPt_z0_F->SetMinimum(0);
2375
2376   h2_resVsPt_phi->SetMinimum(0);
2377   h2_resVsPt_phi_C->SetMinimum(0);
2378   h2_resVsPt_phi_I->SetMinimum(0);
2379   h2_resVsPt_phi_F->SetMinimum(0);
2380
2381   h2_resVsPt_eta->SetMinimum(0);
2382
2383   h2_resVsEta_eta->SetMinimum(0);
2384   h2_resVsEta_eta_L->SetMinimum(0);
2385   h2_resVsEta_eta_H->SetMinimum(0);
2386
2387   h2_resVsEta_z0->SetMinimum(0);
2388   h2_resVsEta_z0_L->SetMinimum(0);
2389   h2_resVsEta_z0_H->SetMinimum(0);
2390
2391   h2_resVsEta_phi->SetMinimum(0);
2392   h2_resVsEta_phi_L->SetMinimum(0);
2393   h2_resVsEta_phi_H->SetMinimum(0);
2394
2395   h2_resVsEta_ptRel->SetMinimum(0);
2396   h2_resVsEta_ptRel_L->SetMinimum(0);
2397   h2_resVsEta_ptRel_H->SetMinimum(0);
2398
2399   h2_resVsPt_d0->SetMinimum(0);
2400   h2_resVsEta_d0->SetMinimum(0);
2401   h2_resVsEta_d0_L->SetMinimum(0);
2402   h2_resVsEta_d0_H->SetMinimum(0);
2403
2404   // -------------------------------------------------------------------------------------------
2405   // output file for histograms
2406   // -------------------------------------------------------------------------------------------
2407
2408   if (TP_select_pdgid != 0) {
2409     char pdgidtxt[500];
2410     sprintf(pdgidtxt, "_pdgid%i", TP_select_pdgid);
2411     type = type + pdgidtxt;
2412   } else if (TP_select_injet == 1)
2413     type = type + "_injet";
2414   else if (TP_select_injet == 2)
2415     type = type + "_injet_highpt";
2416   else if (TP_select_injet == 3)
2417     type = type + "_injet_vhighpt";
2418
2419   if (TP_select_eventid != 0)
2420     type = type + "_wpu";
2421
2422   if (useTightCuts)
2423     type = type + "_tight";
2424   if (useDeadRegion)
2425     type = type + "_dead";
2426
2427   if (L1Tk_seed != 0) {
2428     char seedtxt[500];
2429     sprintf(seedtxt, "_seed%i", L1Tk_seed);
2430     type = type + seedtxt;
2431   }
2432
2433   if (TP_minPt > 2.0) {
2434     char pttxt[500];
2435     sprintf(pttxt, "_pt%.0f", TP_minPt);
2436     type = type + pttxt;
2437   }
2438
2439   TFile* fout;
2440   if (doLooseMatch)
2441     fout = new TFile("output_looseMatch_" + type + treeName + ".root", "recreate");
2442   else
2443     fout = new TFile(type_dir + "output_" + type + treeName + ".root", "recreate");
2444
2445   // -------------------------------------------------------------------------------------------
2446   // draw and save plots
2447   // -------------------------------------------------------------------------------------------
2448
2449   char ctxt[500];
2450   TCanvas c;
2451
2452   gSystem->mkdir("TrkPlots");
2453   TString DIR = "TrkPlots/";
2454
2455   // plots overlaying 68, 90, 99% confidence levels]
2456
2457   float max_eta_ptRel = 0.2;
2458   float max_pt_ptRel = 0.2;
2459   float max_pt_pt = 20;
2460   float max_z0 = 2.0;
2461   float max_phi = 0.01;
2462   float max_eta = 0.03;
2463
2464   if (type.Contains("El")) {
2465     max_pt_ptRel = 1.0;
2466     max_eta_ptRel = 1.0;
2467     max_phi = 0.1;
2468   }
2469
2470   // makeResidualIntervalPlot will save the individual plots to the root file
2471   makeResidualIntervalPlot(
2472       type, DIR, "resVsPt_ptRel", h2_resVsPt_ptRel_68, h2_resVsPt_ptRel_90, h2_resVsPt_ptRel_99, 0, max_pt_ptRel);
2473   makeResidualIntervalPlot(type, DIR, "resVsPt_pt", h2_resVsPt_pt_68, h2_resVsPt_pt_90, h2_resVsPt_pt_99, 0, max_pt_pt);
2474   makeResidualIntervalPlot(type, DIR, "resVsPt_z0", h2_resVsPt_z0_68, h2_resVsPt_z0_90, h2_resVsPt_z0_99, 0, max_z0);
2475   makeResidualIntervalPlot(
2476       type, DIR, "resVsPt_phi", h2_resVsPt_phi_68, h2_resVsPt_phi_90, h2_resVsPt_phi_99, 0, max_phi);
2477   makeResidualIntervalPlot(
2478       type, DIR, "resVsPt_eta", h2_resVsPt_eta_68, h2_resVsPt_eta_90, h2_resVsPt_eta_99, 0, max_eta);
2479   //makeResidualIntervalPlot( type, DIR, "resVsPt_d0", h2_resVsPt_d0_68, h2_resVsPt_d0_90, h2_resVsPt_d0_99, 0, 0.02 );
2480
2481   makeResidualIntervalPlot(type,
2482                            DIR,
2483                            "resVsPt_L_ptRel",
2484                            h2_resVsPt_ptRel_L_68,
2485                            h2_resVsPt_ptRel_L_90,
2486                            h2_resVsPt_ptRel_L_99,
2487                            0,
2488                            max_pt_ptRel);
2489   makeResidualIntervalPlot(type, DIR, "resVsPt_L_pt", h2_resVsPt_pt_L_68, h2_resVsPt_pt_L_90, h2_resVsPt_pt_L_99, 0, 4);
2490   makeResidualIntervalPlot(
2491       type, DIR, "resVsPt_L_z0", h2_resVsPt_z0_L_68, h2_resVsPt_z0_L_90, h2_resVsPt_z0_L_99, 0, max_z0);
2492   makeResidualIntervalPlot(
2493       type, DIR, "resVsPt_L_phi", h2_resVsPt_phi_L_68, h2_resVsPt_phi_L_90, h2_resVsPt_phi_L_99, 0, max_phi);
2494   makeResidualIntervalPlot(
2495       type, DIR, "resVsPt_L_eta", h2_resVsPt_eta_L_68, h2_resVsPt_eta_L_90, h2_resVsPt_eta_L_99, 0, max_eta);
2496   //makeResidualIntervalPlot( type, DIR, "resVsPt_L_d0", h2_resVsPt_d0_L_68, h2_resVsPt_d0_L_90, h2_resVsPt_d0_L_99, 0, 0.02 );
2497
2498   makeResidualIntervalPlot(
2499       type, DIR, "resVsEta_eta", h2_resVsEta_eta_68, h2_resVsEta_eta_90, h2_resVsEta_eta_99, 0, max_eta);
2500   makeResidualIntervalPlot(
2501       type, DIR, "resVsEta_z0", h2_resVsEta_z0_68, h2_resVsEta_z0_90, h2_resVsEta_z0_99, 0, max_z0);
2502   makeResidualIntervalPlot(
2503       type, DIR, "resVsEta_phi", h2_resVsEta_phi_68, h2_resVsEta_phi_90, h2_resVsEta_phi_99, 0, max_phi);
2504   makeResidualIntervalPlot(
2505       type, DIR, "resVsEta_ptRel", h2_resVsEta_ptRel_68, h2_resVsEta_ptRel_90, h2_resVsEta_ptRel_99, 0, max_eta_ptRel);
2506   //makeResidualIntervalPlot( type, DIR, "resVsEta_d0", h2_resVsEta_d0_68, h2_resVsEta_d0_90, h2_resVsEta_d0_99, 0, 0.02 );
2507
2508   makeResidualIntervalPlot(
2509       type, DIR, "resVsEta_L_eta", h2_resVsEta_eta_L_68, h2_resVsEta_eta_L_90, h2_resVsEta_eta_L_99, 0, max_eta);
2510   makeResidualIntervalPlot(
2511       type, DIR, "resVsEta_L_z0", h2_resVsEta_z0_L_68, h2_resVsEta_z0_L_90, h2_resVsEta_z0_L_99, 0, max_z0);
2512   makeResidualIntervalPlot(
2513       type, DIR, "resVsEta_L_phi", h2_resVsEta_phi_L_68, h2_resVsEta_phi_L_90, h2_resVsEta_phi_L_99, 0, max_phi);
2514   makeResidualIntervalPlot(type,
2515                            DIR,
2516                            "resVsEta_L_ptRel",
2517                            h2_resVsEta_ptRel_L_68,
2518                            h2_resVsEta_ptRel_L_90,
2519                            h2_resVsEta_ptRel_L_99,
2520                            0,
2521                            max_eta_ptRel);
2522   //makeResidualIntervalPlot( type, DIR, "resVsEta_L_d0", h2_resVsEta_d0_L_68, h2_resVsEta_d0_L_90, h2_resVsEta_d0_L_99, 0, 0.02 );
2523
2524   makeResidualIntervalPlot(
2525       type, DIR, "resVsEta_H_eta", h2_resVsEta_eta_H_68, h2_resVsEta_eta_H_90, h2_resVsEta_eta_H_99, 0, max_eta);
2526   makeResidualIntervalPlot(
2527       type, DIR, "resVsEta_H_z0", h2_resVsEta_z0_H_68, h2_resVsEta_z0_H_90, h2_resVsEta_z0_H_99, 0, max_z0);
2528   makeResidualIntervalPlot(
2529       type, DIR, "resVsEta_H_phi", h2_resVsEta_phi_H_68, h2_resVsEta_phi_H_90, h2_resVsEta_phi_H_99, 0, max_phi);
2530   makeResidualIntervalPlot(type,
2531                            DIR,
2532                            "resVsEta_H_ptRel",
2533                            h2_resVsEta_ptRel_H_68,
2534                            h2_resVsEta_ptRel_H_90,
2535                            h2_resVsEta_ptRel_H_99,
2536                            0,
2537                            max_eta_ptRel);
2538   //makeResidualIntervalPlot( type, DIR, "resVsEta_H_d0", h2_resVsEta_d0_H_68, h2_resVsEta_d0_H_90, h2_resVsEta_d0_H_99, 0, 0.02 );
2539
2540   if (doDetailedPlots) {
2541     makeResidualIntervalPlot(
2542         type, DIR, "resVsPhi_ptRel", h2_resVsPhi_ptRel_68, h2_resVsPhi_ptRel_90, h2_resVsPhi_ptRel_99, 0, 0.5);
2543     makeResidualIntervalPlot(type, DIR, "resVsPhi_pt", h2_resVsPhi_pt_68, h2_resVsPhi_pt_90, h2_resVsPhi_pt_99, 0, 20);
2544   }
2545
2546   // ----------------------------------------------------------------------------------------------------------
2547   // resoultion vs pt
2548   // ----------------------------------------------------------------------------------------------------------
2549
2550   h2_resVsPt_pt_90->SetMinimum(0);
2551   h2_resVsPt_pt_90->SetMarkerStyle(20);
2552   h2_resVsPt_pt_90->Draw("p");
2553   c.SaveAs(DIR + type + "_resVsPt_pt_90.pdf");
2554
2555   h2_resVsPt_ptRel_90->SetMinimum(0);
2556   h2_resVsPt_ptRel_90->SetMarkerStyle(20);
2557   h2_resVsPt_ptRel_90->Draw("p");
2558   c.SaveAs(DIR + type + "_resVsPt_ptRel_90.pdf");
2559
2560   h2_resVsPt_z0_90->SetMinimum(0);
2561   h2_resVsPt_z0_90->SetMarkerStyle(20);
2562   h2_resVsPt_z0_90->Draw("p");
2563   c.SaveAs(DIR + type + "_resVsPt_z0_90.pdf");
2564
2565   h2_resVsPt_phi_90->SetMinimum(0);
2566   h2_resVsPt_phi_90->SetMarkerStyle(20);
2567   h2_resVsPt_phi_90->Draw("p");
2568   c.SaveAs(DIR + type + "_resVsPt_phi_90.pdf");
2569
2570   h2_resVsPt_eta_90->SetMinimum(0);
2571   h2_resVsPt_eta_90->SetMarkerStyle(20);
2572   h2_resVsPt_eta_90->Draw("p");
2573   c.SaveAs(DIR + type + "_resVsPt_eta_90.pdf");
2574
2575   /*
2576   h2_resVsPt_phi_90->SetMinimum(0);
2577   h2_resVsPt_d0_90->SetMarkerStyle(20);
2578   h2_resVsPt_d0_90->Draw("p");
2579   c.SaveAs(DIR+type+"_resVsPt_d0_90.pdf");
2580   */
2581
2582   // Limit decimal places for doubles in normalized tracklet eff. graphs
2583   gStyle->SetPaintTextFormat("4.2f");
2584
2585   h_trk_tracklet_eff->SetMinimum(0);
2586   h_trk_tracklet_eff->Draw("col text");
2587   c.SaveAs(DIR + type + "_trk_tracklet_eff.pdf");
2588
2589   // Remove decimal places for ints in in tracklet hits graphs
2590   gStyle->SetPaintTextFormat("4.0f");
2591
2592   h_trk_tracklet_hits->SetMinimum(0);
2593   h_trk_tracklet_hits->Draw("col text");
2594   c.SaveAs(DIR + type + "_trk_tracklet_hits.pdf");
2595
2596   if (doDetailedPlots) {
2597     h2_resVsPt_eta->Write();
2598
2599     h2_resVsPt_pt->Write();
2600     h2_resVsPt_pt_C->Write();
2601     h2_resVsPt_pt_I->Write();
2602     h2_resVsPt_pt_F->Write();
2603
2604     h2_resVsPt_ptRel->Write();
2605     h2_resVsPt_ptRel_C->Write();
2606     h2_resVsPt_ptRel_I->Write();
2607     h2_resVsPt_ptRel_F->Write();
2608
2609     h2_mresVsPt_pt->Write();
2610     h2_mresVsPt_pt_C->Write();
2611     h2_mresVsPt_pt_I->Write();
2612     h2_mresVsPt_pt_F->Write();
2613
2614     h2_resVsPt_d0->Write();
2615
2616     h2_resVsPt_z0_C->Write();
2617     h2_resVsPt_z0_I->Write();
2618     h2_resVsPt_z0_F->Write();
2619
2620     h2_resVsPt_phi->Write();
2621     h2_resVsPt_phi_C->Write();
2622     h2_resVsPt_phi_I->Write();
2623     h2_resVsPt_phi_F->Write();
2624   }
2625
2626   // ----------------------------------------------------------------------------------------------------------
2627   // resolution vs eta
2628   // ----------------------------------------------------------------------------------------------------------
2629
2630   h2_resVsEta_eta_90->SetMinimum(0);
2631   h2_resVsEta_eta_90->SetMarkerStyle(20);
2632   h2_resVsEta_eta_90->Draw("p");
2633   c.SaveAs(DIR + type + "_resVsEta_eta_90.pdf");
2634
2635   h2_resVsEta_eta_68->SetMinimum(0);
2636   h2_resVsEta_eta_68->SetMarkerStyle(20);
2637   h2_resVsEta_eta_68->Draw("p");
2638   c.SaveAs(DIR + type + "_resVsEta_eta_68.pdf");
2639
2640   if (doDetailedPlots) {
2641     h2_resVsEta_eta_L_90->Draw("p");
2642     sprintf(ctxt, "p_{T} < 8 GeV");
2643     mySmallText(0.22, 0.82, 1, ctxt);
2644     c.SaveAs(DIR + type + "_resVsEta_eta_L_90.pdf");
2645
2646     h2_resVsEta_eta_H_90->Draw("p");
2647     sprintf(ctxt, "p_{T} > 8 GeV");
2648     mySmallText(0.22, 0.82, 1, ctxt);
2649     c.SaveAs(DIR + type + "_resVsEta_eta_H_90.pdf");
2650   }
2651
2652   h2_resVsEta_z0_90->SetMinimum(0);
2653   h2_resVsEta_z0_90->SetMarkerStyle(20);
2654   h2_resVsEta_z0_90->Draw("p");
2655   c.SaveAs(DIR + type + "_resVsEta_z0_90.pdf");
2656
2657   h2_resVsEta_z0_68->SetMinimum(0);
2658   h2_resVsEta_z0_68->SetMarkerStyle(20);
2659   h2_resVsEta_z0_68->Draw("p");
2660   c.SaveAs(DIR + type + "_resVsEta_z0_68.pdf");
2661
2662   if (doDetailedPlots) {
2663     h2_resVsEta_z0_L_90->Draw();
2664     sprintf(ctxt, "p_{T} < 8 GeV");
2665     mySmallText(0.22, 0.82, 1, ctxt);
2666     c.SaveAs(DIR + type + "_resVsEta_z0_L_90.pdf");
2667
2668     h2_resVsEta_z0_H_90->Draw();
2669     sprintf(ctxt, "p_{T} > 8 GeV");
2670     mySmallText(0.22, 0.82, 1, ctxt);
2671     c.SaveAs(DIR + type + "_resVsEta_z0_H_90.pdf");
2672   }
2673
2674   /*
2675   h2_resVsEta_d0_90->Draw();
2676   c.SaveAs(DIR+type+"_resVsEta_d0_90.pdf");
2677
2678   h2_resVsEta_d0_L_90->Draw();
2679   sprintf(ctxt,"p_{T} < 8 GeV");
2680   mySmallText(0.22,0.82,1,ctxt);
2681   c.SaveAs(DIR+type+"_resVsEta_d0_L_90.pdf");
2682
2683   h2_resVsEta_d0_H_90->Draw();
2684   sprintf(ctxt,"p_{T} > 8 GeV");
2685   mySmallText(0.22,0.82,1,ctxt);
2686   c.SaveAs(DIR+type+"_resVsEta_d0_H_90.pdf");
2687   */
2688
2689   h2_resVsEta_phi_90->SetMinimum(0);
2690   h2_resVsEta_phi_90->SetMarkerStyle(20);
2691   h2_resVsEta_phi_90->Draw("p");
2692   c.SaveAs(DIR + type + "_resVsEta_phi_90.pdf");
2693
2694   h2_resVsEta_phi_68->SetMinimum(0);
2695   h2_resVsEta_phi_68->SetMarkerStyle(20);
2696   h2_resVsEta_phi_68->Draw("p");
2697   c.SaveAs(DIR + type + "_resVsEta_phi_68.pdf");
2698
2699   if (doDetailedPlots) {
2700     h2_resVsEta_phi_L_90->Draw();
2701     sprintf(ctxt, "p_{T} < 8 GeV");
2702     mySmallText(0.22, 0.82, 1, ctxt);
2703     c.SaveAs(DIR + type + "_resVsEta_phi_L_90.pdf");
2704
2705     h2_resVsEta_phi_H_90->Draw();
2706     sprintf(ctxt, "p_{T} > 8 GeV");
2707     mySmallText(0.22, 0.82, 1, ctxt);
2708     c.SaveAs(DIR + type + "_resVsEta_phi_H_90.pdf");
2709   }
2710
2711   h2_resVsEta_ptRel_90->SetMinimum(0);
2712   h2_resVsEta_ptRel_90->SetMarkerStyle(20);
2713   h2_resVsEta_ptRel_90->Draw("p");
2714   c.SaveAs(DIR + type + "_resVsEta_ptRel_90.pdf");
2715
2716   h2_resVsEta_ptRel_68->SetMinimum(0);
2717   h2_resVsEta_ptRel_68->SetMarkerStyle(20);
2718   h2_resVsEta_ptRel_68->Draw("p");
2719   c.SaveAs(DIR + type + "_resVsEta_ptRel_68.pdf");
2720
2721   if (doDetailedPlots) {
2722     h2_resVsEta_ptRel_L_90->Draw();
2723     sprintf(ctxt, "p_{T} < 8 GeV");
2724     mySmallText(0.22, 0.82, 1, ctxt);
2725     c.SaveAs(DIR + type + "_resVsEta_ptRel_L_90.pdf");
2726
2727     h2_resVsEta_ptRel_H_90->Draw();
2728     sprintf(ctxt, "p_{T} > 8 GeV");
2729     mySmallText(0.22, 0.82, 1, ctxt);
2730     c.SaveAs(DIR + type + "_resVsEta_ptRel_H_90.pdf");
2731
2732     h2_resVsEta_eta->Write();
2733     h2_resVsEta_eta_L->Write();
2734     h2_resVsEta_eta_H->Write();
2735
2736     h2_resVsEta_z0->Draw();
2737     c.SaveAs(DIR + type + "_resVsEta_z0_rms.pdf");
2738     h2_resVsEta_eta->Draw();
2739     c.SaveAs(DIR + type + "_resVsEta_eta_rms.pdf");
2740     h2_resVsEta_ptRel->Draw();
2741     c.SaveAs(DIR + type + "_resVsEta_ptRel_rms.pdf");
2742     h2_resVsEta_phi->Draw();
2743     c.SaveAs(DIR + type + "_resVsEta_phi_rms.pdf");
2744
2745     // check residuals
2746     h2_mresVsEta_z0->Draw();
2747     h2_mresVsEta_z0->Write();
2748     c.SaveAs(DIR + type + "_mean-resVsEta_z0.pdf");
2749
2750     h2_mresVsEta_eta->Draw();
2751     h2_mresVsEta_eta->Write();
2752     c.SaveAs(DIR + type + "_mean-resVsEta_eta.pdf");
2753
2754     h2_mresVsEta_ptRel->Draw();
2755     h2_mresVsEta_ptRel->Write();
2756     c.SaveAs(DIR + type + "_mean-resVsEta_ptRel.pdf");
2757
2758     h2_mresVsEta_phi->Draw();
2759     h2_mresVsEta_phi->Write();
2760     c.SaveAs(DIR + type + "_mean-resVsEta_phi.pdf");
2761
2762     h2_resVsEta_z0->Write();
2763     h2_resVsEta_z0_L->Write();
2764     h2_resVsEta_z0_H->Write();
2765
2766     h2_resVsEta_d0->Write();
2767     h2_resVsEta_d0_L->Write();
2768     h2_resVsEta_d0_H->Write();
2769
2770     h2_resVsEta_phi->Write();
2771     h2_resVsEta_phi_L->Write();
2772     h2_resVsEta_phi_H->Write();
2773
2774     h2_resVsEta_ptRel->Write();
2775     h2_resVsEta_ptRel_L->Write();
2776     h2_resVsEta_ptRel_H->Write();
2777   }
2778
2779   if (doGausFit) {
2780     h3_resVsEta_eta_L->Write();
2781     h3_resVsEta_z0_L->Write();
2782     h3_resVsEta_phi_L->Write();
2783     h3_resVsEta_ptRel_L->Write();
2784
2785     h3_resVsEta_eta_H->Write();
2786     h3_resVsEta_z0_H->Write();
2787     h3_resVsEta_phi_H->Write();
2788     h3_resVsEta_ptRel_H->Write();
2789   }
2790
2791   // resolution vs phi
2792   if (doDetailedPlots) {
2793     h2_resVsPhi_pt_90->SetMinimum(0);
2794     h2_resVsPhi_pt_90->SetMarkerStyle(20);
2795     h2_resVsPhi_pt_90->Draw("p");
2796     c.SaveAs(DIR + type + "_resVsPhi_pt_90.pdf");
2797
2798     h2_resVsPhi_ptRel_90->SetMinimum(0);
2799     h2_resVsPhi_ptRel_90->SetMarkerStyle(20);
2800     h2_resVsPhi_ptRel_90->Draw("p");
2801     c.SaveAs(DIR + type + "_resVsPhi_ptRel_90.pdf");
2802   }
2803
2804   // ----------------------------------------------------------------------------------------------------------------
2805   // track quality plots
2806   // ----------------------------------------------------------------------------------------------------------------
2807
2808   if (doDetailedPlots) {
2809     h_match_trk_nstub->Write();
2810     h_match_trk_nstub_C->Write();
2811     h_match_trk_nstub_I->Write();
2812     h_match_trk_nstub_F->Write();
2813   }
2814
2815   h_trk_chi2->Draw();
2816   sprintf(ctxt, "|eta| < 2.4");
2817   mySmallText(0.52, 0.82, 1, ctxt);
2818   c.SaveAs(DIR + type + "_trk_chi2.pdf");
2819
2820   h_trk_chi2_dof->Draw();
2821   sprintf(ctxt, "|eta| < 2.4");
2822   mySmallText(0.52, 0.82, 1, ctxt);
2823   c.SaveAs(DIR + type + "_trk_chi2_dof.pdf");
2824
2825   h_trk_chi2rphi->Draw();
2826   sprintf(ctxt, "|eta| < 2.4");
2827   mySmallText(0.52, 0.82, 1, ctxt);
2828   c.SaveAs(DIR + type + "_trk_chi2rphi.pdf");
2829
2830   h_trk_chi2rphi_dof->Draw();
2831   sprintf(ctxt, "|eta| < 2.4");
2832   mySmallText(0.52, 0.82, 1, ctxt);
2833   c.SaveAs(DIR + type + "_trk_chi2rphi_dof.pdf");
2834
2835   h_trk_chi2rz->Draw();
2836   sprintf(ctxt, "|eta| < 2.4");
2837   mySmallText(0.52, 0.82, 1, ctxt);
2838   c.SaveAs(DIR + type + "_trk_chi2rz.pdf");
2839
2840   h_trk_chi2rz_dof->Draw();
2841   sprintf(ctxt, "|eta| < 2.4");
2842   mySmallText(0.52, 0.82, 1, ctxt);
2843   c.SaveAs(DIR + type + "_trk_chi2rz_dof.pdf");
2844
2845   h_trk_chi2->Write();
2846   h_trk_chi2rphi->Write();
2847   h_trk_chi2rz->Write();
2848   h_match_trk_chi2->Write();
2849   h_match_trk_chi2rphi->Write();
2850   h_match_trk_chi2rz->Write();
2851
2852   if (doDetailedPlots) {
2853     h_match_trk_chi2_C_L->Write();
2854     h_match_trk_chi2_I_L->Write();
2855     h_match_trk_chi2_F_L->Write();
2856     h_match_trk_chi2_C_H->Write();
2857     h_match_trk_chi2_I_H->Write();
2858     h_match_trk_chi2_F_H->Write();
2859
2860     h_match_trk_chi2_dof_C_L->Write();
2861     h_match_trk_chi2_dof_I_L->Write();
2862     h_match_trk_chi2_dof_F_L->Write();
2863     h_match_trk_chi2_dof_C_H->Write();
2864     h_match_trk_chi2_dof_I_H->Write();
2865     h_match_trk_chi2_dof_F_H->Write();
2866   }
2867
2868   // ----------------------------------------------------------------------------------------------------------------
2869   // efficiency plots
2870   // ----------------------------------------------------------------------------------------------------------------
2871
2872   // rebin pt/phi plots
2873   h_tp_pt->Rebin(2);
2874   h_match_tp_pt->Rebin(2);
2875   h_tp_phi->Rebin(2);
2876   h_match_tp_phi->Rebin(2);
2877
2878   h_tp_pt_L->Rebin(2);
2879   h_match_tp_pt_L->Rebin(2);
2880   h_tp_pt_LC->Rebin(2);
2881   h_match_tp_pt_LC->Rebin(2);
2882   h_tp_pt_H->Rebin(2);
2883   h_match_tp_pt_H->Rebin(2);
2884
2885   // h_tp_eta->Rebin(2);
2886   // h_match_tp_eta->Rebin(2);
2887   // h_tp_eta_L->Rebin(2);
2888   // h_match_tp_eta_L->Rebin(2);
2889   // h_tp_eta_H->Rebin(2);
2890   // h_match_tp_eta_H->Rebin(2);
2891
2892   // calculate the efficiency
2893   h_match_tp_pt->Sumw2();
2894   h_tp_pt->Sumw2();
2895   TH1F* h_eff_pt = (TH1F*)h_match_tp_pt->Clone();
2896   h_eff_pt->SetName("eff_pt");
2897   h_eff_pt->GetYaxis()->SetTitle("Efficiency");
2898   h_eff_pt->Divide(h_match_tp_pt, h_tp_pt, 1.0, 1.0, "B");
2899
2900   h_match_tp_pt_L->Sumw2();
2901   h_tp_pt_L->Sumw2();
2902   TH1F* h_eff_pt_L = (TH1F*)h_match_tp_pt_L->Clone();
2903   h_eff_pt_L->SetName("eff_pt_L");
2904   h_eff_pt_L->GetYaxis()->SetTitle("Efficiency");
2905   h_eff_pt_L->Divide(h_match_tp_pt_L, h_tp_pt_L, 1.0, 1.0, "B");
2906
2907   h_match_tp_pt_LC->Sumw2();
2908   h_tp_pt_LC->Sumw2();
2909   TH1F* h_eff_pt_LC = (TH1F*)h_match_tp_pt_LC->Clone();
2910   h_eff_pt_LC->SetName("eff_pt_LC");
2911   h_eff_pt_LC->GetYaxis()->SetTitle("Efficiency");
2912   h_eff_pt_LC->Divide(h_match_tp_pt_LC, h_tp_pt_LC, 1.0, 1.0, "B");
2913
2914   h_match_tp_pt_H->Sumw2();
2915   h_tp_pt_H->Sumw2();
2916   TH1F* h_eff_pt_H = (TH1F*)h_match_tp_pt_H->Clone();
2917   h_eff_pt_H->SetName("eff_pt_H");
2918   h_eff_pt_H->GetYaxis()->SetTitle("Efficiency");
2919   h_eff_pt_H->Divide(h_match_tp_pt_H, h_tp_pt_H, 1.0, 1.0, "B");
2920
2921   h_match_tp_eta->Sumw2();
2922   h_tp_eta->Sumw2();
2923   TH1F* h_eff_eta = (TH1F*)h_match_tp_eta->Clone();
2924   h_eff_eta->SetName("eff_eta");
2925   h_eff_eta->GetYaxis()->SetTitle("Efficiency");
2926   h_eff_eta->Divide(h_match_tp_eta, h_tp_eta, 1.0, 1.0, "B");
2927
2928   h_match_tp_eta_L->Sumw2();
2929   h_tp_eta_L->Sumw2();
2930   TH1F* h_eff_eta_L = (TH1F*)h_match_tp_eta_L->Clone();
2931   h_eff_eta_L->SetName("eff_eta_L");
2932   h_eff_eta_L->GetYaxis()->SetTitle("Efficiency");
2933   h_eff_eta_L->Divide(h_match_tp_eta_L, h_tp_eta_L, 1.0, 1.0, "B");
2934
2935   h_match_tp_eta_H->Sumw2();
2936   h_tp_eta_H->Sumw2();
2937   TH1F* h_eff_eta_H = (TH1F*)h_match_tp_eta_H->Clone();
2938   h_eff_eta_H->SetName("eff_eta_H");
2939   h_eff_eta_H->GetYaxis()->SetTitle("Efficiency");
2940   h_eff_eta_H->Divide(h_match_tp_eta_H, h_tp_eta_H, 1.0, 1.0, "B");
2941
2942   h_match_tp_eta_23->Sumw2();
2943   h_tp_eta_23->Sumw2();
2944   TH1F* h_eff_eta_23 = (TH1F*)h_match_tp_eta_23->Clone();
2945   h_eff_eta_23->SetName("eff_eta_23");
2946   h_eff_eta_23->GetYaxis()->SetTitle("Efficiency");
2947   h_eff_eta_23->Divide(h_match_tp_eta_23, h_tp_eta_23, 1.0, 1.0, "B");
2948
2949   h_match_tp_eta_35->Sumw2();
2950   h_tp_eta_35->Sumw2();
2951   TH1F* h_eff_eta_35 = (TH1F*)h_match_tp_eta_35->Clone();
2952   h_eff_eta_35->SetName("eff_eta_35");
2953   h_eff_eta_35->GetYaxis()->SetTitle("Efficiency");
2954   h_eff_eta_35->Divide(h_match_tp_eta_35, h_tp_eta_35, 1.0, 1.0, "B");
2955
2956   h_match_tp_eta_5->Sumw2();
2957   h_tp_eta_5->Sumw2();
2958   TH1F* h_eff_eta_5 = (TH1F*)h_match_tp_eta_5->Clone();
2959   h_eff_eta_5->SetName("eff_eta_5");
2960   h_eff_eta_5->GetYaxis()->SetTitle("Efficiency");
2961   h_eff_eta_5->Divide(h_match_tp_eta_5, h_tp_eta_5, 1.0, 1.0, "B");
2962
2963   h_match_tp_phi->Sumw2();
2964   h_tp_phi->Sumw2();
2965   TH1F* h_eff_phi = (TH1F*)h_match_tp_phi->Clone();
2966   h_eff_phi->SetName("eff_phi");
2967   h_eff_phi->GetYaxis()->SetTitle("Efficiency");
2968   h_eff_phi->Divide(h_match_tp_phi, h_tp_phi, 1.0, 1.0, "B");
2969
2970   h_match_tp_z0->Sumw2();
2971   h_tp_z0->Sumw2();
2972   TH1F* h_eff_z0 = (TH1F*)h_match_tp_z0->Clone();
2973   h_eff_z0->SetName("eff_z0");
2974   h_eff_z0->GetYaxis()->SetTitle("Efficiency");
2975   h_eff_z0->Divide(h_match_tp_z0, h_tp_z0, 1.0, 1.0, "B");
2976
2977   h_match_tp_z0_L->Sumw2();
2978   h_tp_z0_L->Sumw2();
2979   TH1F* h_eff_z0_L = (TH1F*)h_match_tp_z0_L->Clone();
2980   h_eff_z0_L->SetName("eff_z0_L");
2981   h_eff_z0_L->GetYaxis()->SetTitle("Efficiency");
2982   h_eff_z0_L->Divide(h_match_tp_z0_L, h_tp_z0_L, 1.0, 1.0, "B");
2983
2984   h_match_tp_z0_H->Sumw2();
2985   h_tp_z0_H->Sumw2();
2986   TH1F* h_eff_z0_H = (TH1F*)h_match_tp_z0_H->Clone();
2987   h_eff_z0_H->SetName("eff_z0_H");
2988   h_eff_z0_H->GetYaxis()->SetTitle("Efficiency");
2989   h_eff_z0_H->Divide(h_match_tp_z0_H, h_tp_z0_H, 1.0, 1.0, "B");
2990
2991   h_match_tp_d0->Sumw2();
2992   h_tp_d0->Sumw2();
2993   TH1F* h_eff_d0 = (TH1F*)h_match_tp_d0->Clone();
2994   h_eff_d0->SetName("eff_d0");
2995   h_eff_d0->GetYaxis()->SetTitle("Efficiency");
2996   h_eff_d0->Divide(h_match_tp_d0, h_tp_d0, 1.0, 1.0, "B");
2997
2998   h_match_tp_absd0->Sumw2();
2999   h_tp_absd0->Sumw2();
3000   TH1F* h_eff_absd0 = (TH1F*)h_match_tp_absd0->Clone();
3001   h_eff_absd0->SetName("eff_absd0");
3002   h_eff_absd0->GetYaxis()->SetTitle("Efficiency");
3003   h_eff_absd0->Divide(h_match_tp_absd0, h_tp_absd0, 1.0, 1.0, "B");
3004
3005   h_match_tp_absd0_eta2->Sumw2();
3006   h_tp_absd0_eta2->Sumw2();
3007   TH1F* h_eff_absd0_eta2 = (TH1F*)h_match_tp_absd0_eta2->Clone();
3008   h_eff_absd0_eta2->SetName("eff_absd0_eta2");
3009   h_eff_absd0_eta2->GetYaxis()->SetTitle("Efficiency");
3010   h_eff_absd0_eta2->Divide(h_match_tp_absd0_eta2, h_tp_absd0_eta2, 1.0, 1.0, "B");
3011
3012   h_match_tp_absd0_eta2_pt3->Sumw2();
3013   h_tp_absd0_eta2_pt3->Sumw2();
3014   TH1F* h_eff_absd0_eta2_pt3 = (TH1F*)h_match_tp_absd0_eta2_pt3->Clone();
3015   h_eff_absd0_eta2_pt3->SetName("eff_absd0_eta2_pt3");
3016   h_eff_absd0_eta2_pt3->GetYaxis()->SetTitle("Efficiency");
3017   h_eff_absd0_eta2_pt3->Divide(h_match_tp_absd0_eta2_pt3, h_tp_absd0_eta2_pt3, 1.0, 1.0, "B");
3018
3019   // set the axis range
3020   h_eff_pt->SetAxisRange(0, 1.1, "Y");
3021   h_eff_pt_L->SetAxisRange(0, 1.1, "Y");
3022   h_eff_pt_LC->SetAxisRange(0, 1.1, "Y");
3023   h_eff_pt_H->SetAxisRange(0, 1.1, "Y");
3024   h_eff_eta->SetAxisRange(0, 1.1, "Y");
3025   h_eff_eta_L->SetAxisRange(0, 1.1, "Y");
3026   h_eff_eta_H->SetAxisRange(0, 1.1, "Y");
3027   h_eff_eta_23->SetAxisRange(0, 1.1, "Y");
3028   h_eff_eta_35->SetAxisRange(0, 1.1, "Y");
3029   h_eff_eta_5->SetAxisRange(0, 1.1, "Y");
3030   h_eff_phi->SetAxisRange(0, 1.1, "Y");
3031   h_eff_z0->SetAxisRange(0, 1.1, "Y");
3032   h_eff_z0_L->SetAxisRange(0, 1.1, "Y");
3033   h_eff_z0_H->SetAxisRange(0, 1.1, "Y");
3034   h_eff_d0->SetAxisRange(0, 1.1, "Y");
3035   h_eff_absd0->SetAxisRange(0, 1.1, "Y");
3036   h_eff_absd0_eta2->SetAxisRange(0, 1.1, "Y");
3037   h_eff_absd0_eta2_pt3->SetAxisRange(0, 1.1, "Y");
3038
3039   gPad->SetGridx();
3040   gPad->SetGridy();
3041
3042   // draw and save plots
3043   h_eff_pt->Draw();
3044   h_eff_pt->Write();
3045   c.SaveAs(DIR + type + "_eff_pt.pdf");
3046
3047   if (type.Contains("Mu")) {
3048     h_eff_pt->GetYaxis()->SetRangeUser(0.8, 1.01);  // zoomed-in plot
3049     c.SaveAs(DIR + type + "_eff_pt_zoom.pdf");
3050   }
3051
3052   h_eff_pt_L->Draw();
3053   h_eff_pt_L->Write();
3054   sprintf(ctxt, "p_{T} < 8 GeV");
3055   mySmallText(0.45, 0.5, 1, ctxt);
3056   c.SaveAs(DIR + type + "_eff_pt_L.pdf");
3057
3058   if (doDetailedPlots) {
3059     h_eff_pt_LC->Draw();
3060     h_eff_pt_LC->Write();
3061     sprintf(ctxt, "p_{T} < 8 GeV, |#eta|<1.0");
3062     mySmallText(0.45, 0.5, 1, ctxt);
3063     c.SaveAs(DIR + type + "_eff_pt_LC.pdf");
3064   }
3065   h_eff_pt_H->Draw();
3066   h_eff_pt_H->Write();
3067   sprintf(ctxt, "p_{T} > 8 GeV");
3068   mySmallText(0.45, 0.5, 1, ctxt);
3069   c.SaveAs(DIR + type + "_eff_pt_H.pdf");
3070
3071   h_eff_eta->Draw();
3072   h_eff_eta->Write();
3073   c.SaveAs(DIR + type + "_eff_eta.pdf");
3074
3075   if (type.Contains("Mu")) {
3076     h_eff_eta->GetYaxis()->SetRangeUser(0.8, 1.01);  // zoomed-in plot
3077     c.SaveAs(DIR + type + "_eff_eta_zoom.pdf");
3078   }
3079
3080   h_eff_eta_L->Draw();
3081   h_eff_eta_L->Write();
3082   sprintf(ctxt, "p_{T} < 8 GeV");
3083   mySmallText(0.45, 0.5, 1, ctxt);
3084   c.SaveAs(DIR + type + "_eff_eta_L.pdf");
3085
3086   h_eff_eta_H->Draw();
3087   h_eff_eta_H->Write();
3088   sprintf(ctxt, "p_{T} > 8 GeV");
3089   mySmallText(0.45, 0.5, 1, ctxt);
3090   c.SaveAs(DIR + type + "_eff_eta_H.pdf");
3091
3092   h_eff_eta_23->Write();
3093   h_eff_eta_35->Write();
3094   h_eff_eta_5->Write();
3095
3096   if (doDetailedPlots) {
3097     h_eff_eta_23->Draw();
3098     sprintf(ctxt, "2 < p_{T} < 3 GeV");
3099     mySmallText(0.45, 0.5, 1, ctxt);
3100     c.SaveAs(DIR + type + "_eff_eta_23.pdf");
3101
3102     h_eff_eta_35->Draw();
3103     sprintf(ctxt, "3 < p_{T} < 5 GeV");
3104     mySmallText(0.45, 0.5, 1, ctxt);
3105     c.SaveAs(DIR + type + "_eff_eta_35.pdf");
3106
3107     h_eff_eta_5->Draw();
3108     sprintf(ctxt, "p_{T} > 5 GeV");
3109     mySmallText(0.45, 0.5, 1, ctxt);
3110     c.SaveAs(DIR + type + "_eff_eta_5.pdf");
3111
3112     h_eff_z0->Draw();
3113     h_eff_z0->Write();
3114     c.SaveAs(DIR + type + "_eff_z0.pdf");
3115
3116     h_eff_z0_L->Write();
3117     h_eff_z0_H->Write();
3118
3119     h_eff_phi->Draw();
3120     h_eff_phi->Write();
3121     c.SaveAs(DIR + type + "_eff_phi.pdf");
3122
3123     if (type.Contains("Mu")) {
3124       h_eff_phi->GetYaxis()->SetRangeUser(0.8, 1.01);  // zoomed-in plot
3125       c.SaveAs(DIR + type + "_eff_phi_zoom.pdf");
3126     }
3127   }
3128
3129   if (doDetailedPlots || TP_maxD0 > 1.0) {
3130     h_eff_d0->Write();
3131     h_eff_absd0->Write();
3132     h_eff_absd0_eta2->Write();
3133     h_eff_absd0_eta2_pt3->Write();
3134
3135     h_eff_d0->Draw();
3136     c.SaveAs(DIR + type + "_eff_d0.pdf");
3137
3138     h_eff_absd0->Draw();
3139     c.SaveAs(DIR + type + "_eff_absd0.pdf");
3140
3141     h_eff_absd0_eta2->Draw();
3142     c.SaveAs(DIR + type + "_eff_absd0_eta2.pdf");
3143
3144     h_eff_absd0_eta2_pt3->Draw();
3145     c.SaveAs(DIR + type + "_eff_absd0_eta2_pt3.pdf");
3146   }
3147
3148   gPad->SetGridx(0);
3149   gPad->SetGridy(0);
3150
3151   // ----------------------------------------------------------------------------------------------------------------
3152   // more resolution plots
3153   // ----------------------------------------------------------------------------------------------------------------
3154
3155   float rms = 0;
3156
3157   if (doDetailedPlots) {
3158     // draw and save plots
3159     h_res_pt->Draw();
3160     rms = h_res_pt->GetRMS();
3161     sprintf(ctxt, "RMS = %.4f", rms);
3162     mySmallText(0.22, 0.82, 1, ctxt);
3163     c.SaveAs(DIR + type + "_res_pt.pdf");
3164
3165     h_res_ptRel->Draw();
3166     rms = h_res_ptRel->GetRMS();
3167     sprintf(ctxt, "RMS = %.4f", rms);
3168     mySmallText(0.22, 0.82, 1, ctxt);
3169     c.SaveAs(DIR + type + "_res_ptRel.pdf");
3170
3171     h_res_eta->Draw();
3172     rms = h_res_eta->GetRMS();
3173     sprintf(ctxt, "RMS = %.3e", rms);
3174     mySmallText(0.22, 0.82, 1, ctxt);
3175     c.SaveAs(DIR + type + "_res_eta.pdf");
3176
3177     h_res_phi->Draw();
3178     rms = h_res_phi->GetRMS();
3179     sprintf(ctxt, "RMS = %.3e", rms);
3180     mySmallText(0.22, 0.82, 1, ctxt);
3181     c.SaveAs(DIR + type + "_res_phi.pdf");
3182
3183     h_res_z0->Draw();
3184     rms = h_res_z0->GetRMS();
3185     sprintf(ctxt, "RMS = %.4f", rms);
3186     mySmallText(0.22, 0.82, 1, ctxt);
3187     c.SaveAs(DIR + type + "_res_z0.pdf");
3188
3189     h_res_z0_C->Draw();
3190     rms = h_res_z0_C->GetRMS();
3191     sprintf(ctxt, "RMS = %.4f;", rms);
3192     mySmallText(0.22, 0.82, 1, ctxt);
3193     sprintf(ctxt, "|eta| < 0.8");
3194     mySmallText(0.22, 0.76, 1, ctxt);
3195     c.SaveAs(DIR + type + "_res_z0_C.pdf");
3196
3197     h_res_z0_I->Draw();
3198     rms = h_res_z0_I->GetRMS();
3199     sprintf(ctxt, "RMS = %.4f;", rms);
3200     mySmallText(0.22, 0.82, 1, ctxt);
3201     sprintf(ctxt, "0.8 < |eta| < 1.6");
3202     mySmallText(0.22, 0.76, 1, ctxt);
3203     c.SaveAs(DIR + type + "_res_z0_I.pdf");
3204
3205     h_res_z0_F->Draw();
3206     rms = h_res_z0_F->GetRMS();
3207     sprintf(ctxt, "RMS = %.4f;", rms);
3208     mySmallText(0.22, 0.82, 1, ctxt);
3209     sprintf(ctxt, "|eta| > 1.6");
3210     mySmallText(0.22, 0.76, 1, ctxt);
3211     c.SaveAs(DIR + type + "_res_z0_F.pdf");
3212
3213     h_res_z0_C_L->Draw();
3214     h_res_z0_C_L->Write();
3215     rms = h_res_z0_C_L->GetRMS();
3216     sprintf(ctxt, "RMS = %.4f;", rms);
3217     mySmallText(0.22, 0.82, 1, ctxt);
3218     sprintf(ctxt, "|eta| < 0.8");
3219     mySmallText(0.22, 0.76, 1, ctxt);
3220     c.SaveAs(DIR + type + "_res_z0_C_L.pdf");
3221
3222     h_res_z0_I_L->Draw();
3223     h_res_z0_I_L->Write();
3224     rms = h_res_z0_I_L->GetRMS();
3225     sprintf(ctxt, "RMS = %.4f;", rms);
3226     mySmallText(0.22, 0.82, 1, ctxt);
3227     sprintf(ctxt, "0.8 < |eta| < 1.6");
3228     mySmallText(0.22, 0.76, 1, ctxt);
3229     c.SaveAs(DIR + type + "_res_z0_I_L.pdf");
3230
3231     h_res_z0_F_L->Draw();
3232     h_res_z0_F_L->Write();
3233     rms = h_res_z0_F_L->GetRMS();
3234     sprintf(ctxt, "RMS = %.4f;", rms);
3235     mySmallText(0.22, 0.82, 1, ctxt);
3236     sprintf(ctxt, "|eta| > 1.6");
3237     mySmallText(0.22, 0.76, 1, ctxt);
3238     c.SaveAs(DIR + type + "_res_z0_F_L.pdf");
3239
3240     h_res_z0_C_H->Draw();
3241     h_res_z0_C_H->Write();
3242     rms = h_res_z0_C_H->GetRMS();
3243     sprintf(ctxt, "RMS = %.4f;", rms);
3244     mySmallText(0.22, 0.82, 1, ctxt);
3245     sprintf(ctxt, "|eta| < 0.8");
3246     mySmallText(0.22, 0.76, 1, ctxt);
3247     c.SaveAs(DIR + type + "_res_z0_C_H.pdf");
3248
3249     h_res_z0_I_H->Draw();
3250     h_res_z0_I_H->Write();
3251     rms = h_res_z0_I_H->GetRMS();
3252     sprintf(ctxt, "RMS = %.4f;", rms);
3253     mySmallText(0.22, 0.82, 1, ctxt);
3254     sprintf(ctxt, "0.8 < |eta| < 1.6");
3255     mySmallText(0.22, 0.76, 1, ctxt);
3256     c.SaveAs(DIR + type + "_res_z0_I_H.pdf");
3257
3258     h_res_z0_F_H->Draw();
3259     h_res_z0_F_H->Write();
3260     rms = h_res_z0_F_H->GetRMS();
3261     sprintf(ctxt, "RMS = %.4f;", rms);
3262     mySmallText(0.22, 0.82, 1, ctxt);
3263     sprintf(ctxt, "|eta| > 1.6");
3264     mySmallText(0.22, 0.76, 1, ctxt);
3265     c.SaveAs(DIR + type + "_res_z0_F_H.pdf");
3266
3267     h_res_z0_L->Draw();
3268     h_res_z0_L->Write();
3269     rms = h_res_z0_L->GetRMS();
3270     sprintf(ctxt, "RMS = %.4f;", rms);
3271     mySmallText(0.22, 0.82, 1, ctxt);
3272     sprintf(ctxt, "p_{T} < 5 GeV");
3273     mySmallText(0.22, 0.76, 1, ctxt);
3274     c.SaveAs(DIR + type + "_res_z0_L.pdf");
3275
3276     h_res_z0_H->Draw();
3277     h_res_z0_H->Write();
3278     rms = h_res_z0_H->GetRMS();
3279     sprintf(ctxt, "RMS = %.4f;", rms);
3280     mySmallText(0.22, 0.82, 1, ctxt);
3281     sprintf(ctxt, "p_{T} > 15 GeV");
3282     mySmallText(0.22, 0.76, 1, ctxt);
3283     c.SaveAs(DIR + type + "_res_z0_H.pdf");
3284
3285     if (h_res_d0->GetEntries() > 0) {
3286       h_res_d0->Draw();
3287       rms = h_res_d0->GetRMS();
3288       sprintf(ctxt, "RMS = %.4f", rms);
3289       mySmallText(0.22, 0.82, 1, ctxt);
3290       c.SaveAs(DIR + type + "_res_d0.pdf");
3291
3292       h_res_d0_C->Draw();
3293       h_res_d0_C->Write();
3294       rms = h_res_d0_C->GetRMS();
3295       sprintf(ctxt, "RMS = %.4f;", rms);
3296       mySmallText(0.22, 0.82, 1, ctxt);
3297       sprintf(ctxt, "|eta| < 0.8");
3298       mySmallText(0.22, 0.76, 1, ctxt);
3299       c.SaveAs(DIR + type + "_res_d0_C.pdf");
3300
3301       h_res_d0_I->Draw();
3302       h_res_d0_I->Write();
3303       rms = h_res_d0_I->GetRMS();
3304       sprintf(ctxt, "RMS = %.4f;", rms);
3305       mySmallText(0.22, 0.82, 1, ctxt);
3306       sprintf(ctxt, "0.8 < |eta| < 1.6");
3307       mySmallText(0.22, 0.76, 1, ctxt);
3308       c.SaveAs(DIR + type + "_res_d0_I.pdf");
3309
3310       h_res_d0_F->Draw();
3311       h_res_d0_F->Write();
3312       rms = h_res_d0_F->GetRMS();
3313       sprintf(ctxt, "RMS = %.4f;", rms);
3314       mySmallText(0.22, 0.82, 1, ctxt);
3315       sprintf(ctxt, "|eta| > 1.6");
3316       mySmallText(0.22, 0.76, 1, ctxt);
3317       c.SaveAs(DIR + type + "_res_d0_F.pdf");
3318
3319       h_res_d0_C_L->Draw();
3320       h_res_d0_C_L->Write();
3321       rms = h_res_d0_C_L->GetRMS();
3322       sprintf(ctxt, "RMS = %.4f;", rms);
3323       mySmallText(0.22, 0.82, 1, ctxt);
3324       sprintf(ctxt, "|eta| < 0.8");
3325       mySmallText(0.22, 0.76, 1, ctxt);
3326       c.SaveAs(DIR + type + "_res_d0_C_L.pdf");
3327
3328       h_res_d0_I_L->Draw();
3329       h_res_d0_I_L->Write();
3330       rms = h_res_d0_I_L->GetRMS();
3331       sprintf(ctxt, "RMS = %.4f;", rms);
3332       mySmallText(0.22, 0.82, 1, ctxt);
3333       sprintf(ctxt, "0.8 < |eta| < 1.6");
3334       mySmallText(0.22, 0.76, 1, ctxt);
3335       c.SaveAs(DIR + type + "_res_d0_I_L.pdf");
3336
3337       h_res_d0_F_L->Draw();
3338       h_res_d0_F_L->Write();
3339       rms = h_res_d0_F_L->GetRMS();
3340       sprintf(ctxt, "RMS = %.4f;", rms);
3341       mySmallText(0.22, 0.82, 1, ctxt);
3342       sprintf(ctxt, "|eta| > 1.6");
3343       mySmallText(0.22, 0.76, 1, ctxt);
3344       c.SaveAs(DIR + type + "_res_d0_F_L.pdf");
3345
3346       h_res_d0_C_H->Draw();
3347       h_res_d0_C_H->Write();
3348       rms = h_res_d0_C_H->GetRMS();
3349       sprintf(ctxt, "RMS = %.4f;", rms);
3350       mySmallText(0.22, 0.82, 1, ctxt);
3351       sprintf(ctxt, "|eta| < 0.8");
3352       mySmallText(0.22, 0.76, 1, ctxt);
3353       c.SaveAs(DIR + type + "_res_d0_C_H.pdf");
3354
3355       h_res_d0_I_H->Draw();
3356       h_res_d0_I_H->Write();
3357       rms = h_res_d0_I_H->GetRMS();
3358       sprintf(ctxt, "RMS = %.4f;", rms);
3359       mySmallText(0.22, 0.82, 1, ctxt);
3360       sprintf(ctxt, "0.8 < |eta| < 1.6");
3361       mySmallText(0.22, 0.76, 1, ctxt);
3362       c.SaveAs(DIR + type + "_res_d0_I_H.pdf");
3363
3364       h_res_d0_F_H->Draw();
3365       h_res_d0_F_H->Write();
3366       rms = h_res_d0_F_H->GetRMS();
3367       sprintf(ctxt, "RMS = %.4f;", rms);
3368       mySmallText(0.22, 0.82, 1, ctxt);
3369       sprintf(ctxt, "|eta| > 1.6");
3370       mySmallText(0.22, 0.76, 1, ctxt);
3371       c.SaveAs(DIR + type + "_res_d0_F_H.pdf");
3372
3373       h_res_d0_L->Draw();
3374       h_res_d0_L->Write();
3375       rms = h_res_d0_L->GetRMS();
3376       sprintf(ctxt, "RMS = %.4f;", rms);
3377       mySmallText(0.22, 0.82, 1, ctxt);
3378       sprintf(ctxt, "p_{T} < 5 GeV");
3379       mySmallText(0.22, 0.76, 1, ctxt);
3380       c.SaveAs(DIR + type + "_res_d0_L.pdf");
3381
3382       h_res_d0_H->Draw();
3383       h_res_d0_H->Write();
3384       rms = h_res_d0_H->GetRMS();
3385       sprintf(ctxt, "RMS = %.4f;", rms);
3386       mySmallText(0.22, 0.82, 1, ctxt);
3387       sprintf(ctxt, "p_{T} > 15 GeV");
3388       mySmallText(0.22, 0.76, 1, ctxt);
3389       c.SaveAs(DIR + type + "_res_d0_H.pdf");
3390     }
3391   }
3392
3393   // ---------------------------------------------------------------------------------------------------------
3394   // "fake rates"
3395
3396   h_trk_all_vspt->Sumw2();
3397   h_trk_loose_vspt->Sumw2();
3398   h_trk_genuine_vspt->Sumw2();
3399   h_trk_notloose_vspt->Sumw2();
3400   h_trk_notgenuine_vspt->Sumw2();
3401   h_trk_duplicate_vspt->Sumw2();
3402   h_tp_vspt->Sumw2();
3403
3404   // fraction of not genuine tracks
3405   TH1F* h_notgenuine_pt = (TH1F*)h_trk_notgenuine_vspt->Clone();
3406   h_notgenuine_pt->SetName("notgenuine_pt");
3407   h_notgenuine_pt->GetYaxis()->SetTitle("Not genuine fraction");
3408   h_notgenuine_pt->Divide(h_trk_notgenuine_vspt, h_trk_all_vspt, 1.0, 1.0, "B");
3409
3410   h_notgenuine_pt->Write();
3411   h_notgenuine_pt->Draw();
3412   c.SaveAs(DIR + type + "_notgenuine.pdf");
3413
3414   // fraction of not loosely genuine tracks
3415   TH1F* h_notloose_pt = (TH1F*)h_trk_notloose_vspt->Clone();
3416   h_notloose_pt->SetName("notloose_pt");
3417   h_notloose_pt->GetYaxis()->SetTitle("Not loose fraction");
3418   h_notloose_pt->Divide(h_trk_notloose_vspt, h_trk_all_vspt, 1.0, 1.0, "B");
3419
3420   h_notloose_pt->Write();
3421   h_notloose_pt->Draw();
3422   c.SaveAs(DIR + type + "_notloose.pdf");
3423
3424   // fraction of DUPLICATE tracks (genuine and not matched)
3425   TH1F* h_duplicatefrac_pt = (TH1F*)h_trk_duplicate_vspt->Clone();
3426   h_duplicatefrac_pt->SetName("duplicatefrac_pt");
3427   h_duplicatefrac_pt->GetYaxis()->SetTitle("Duplicate fraction");
3428   h_duplicatefrac_pt->Divide(h_trk_duplicate_vspt, h_trk_all_vspt, 1.0, 1.0, "B");
3429
3430   h_duplicatefrac_pt->Write();
3431   h_duplicatefrac_pt->Draw();
3432   c.SaveAs(DIR + type + "_duplicatefrac.pdf");
3433
3434   // ---------------------------------------------------------------------------------------------------------
3435   // total track rates vs pt
3436
3437   h_trk_all_vspt->Scale(1.0 / nevt);
3438   h_trk_loose_vspt->Scale(1.0 / nevt);
3439   h_trk_genuine_vspt->Scale(1.0 / nevt);
3440   h_trk_notloose_vspt->Scale(1.0 / nevt);
3441   h_trk_notgenuine_vspt->Scale(1.0 / nevt);
3442   h_trk_duplicate_vspt->Scale(1.0 / nevt);
3443   h_tp_vspt->Scale(1.0 / nevt);
3444
3445   h_tp_vspt->GetYaxis()->SetTitle("Tracks / event");
3446   h_tp_vspt->GetXaxis()->SetTitle("Track p_{T} [GeV]");
3447   h_tp_vspt->SetLineColor(4);
3448   h_tp_vspt->SetLineStyle(2);
3449
3450   h_trk_notgenuine_vspt->SetLineColor(2);
3451   h_trk_notgenuine_vspt->SetLineStyle(1);
3452
3453   h_trk_duplicate_vspt->SetLineColor(8);
3454   h_trk_duplicate_vspt->SetLineStyle(2);
3455
3456   float max = h_tp_vspt->GetMaximum();
3457   if (h_trk_all_vspt->GetMaximum() > max)
3458     max = h_trk_all_vspt->GetMaximum();
3459   h_tp_vspt->SetAxisRange(0.001, max * 1.5, "Y");
3460
3461   h_tp_vspt->Draw("hist");
3462   h_trk_all_vspt->Draw("same,hist");
3463   h_tp_vspt->Draw("same,hist");
3464   h_trk_notgenuine_vspt->Draw("same,hist");
3465   //h_trk_duplicate_vspt->Draw("same,hist");
3466
3467   h_trk_all_vspt->Write();
3468   h_trk_loose_vspt->Write();
3469   h_trk_genuine_vspt->Write();
3470   h_trk_notloose_vspt->Write();
3471   h_trk_notgenuine_vspt->Write();
3472   h_trk_duplicate_vspt->Write();
3473   h_tp_vspt->Write();
3474
3475   char txt[500];
3476   sprintf(txt, "# tracks/event = %.1f", h_trk_all_vspt->GetSum());
3477   mySmallText(0.5, 0.85, 1, txt);
3478   char txt3[500];
3479   sprintf(txt3, "# TPs(stubs in #geq 4 layers)/");
3480   char txt2[500];
3481   sprintf(txt2, "event = %.1f", h_tp_vspt->GetSum());
3482   mySmallText(0.5, 0.79, 4, txt3);
3483   mySmallText(0.5, 0.74, 4, txt2);
3484
3485   sprintf(txt, "# !genuine tracks/event = %.1f", h_trk_notgenuine_vspt->GetSum());
3486   mySmallText(0.5, 0.69, 2, txt);
3487   //sprintf(txt,"# duplicates/event = %.1f",h_trk_duplicate_vspt->GetSum());
3488   //mySmallText(0.5,0.64,8,txt);
3489
3490   c.SaveAs(DIR + type + "_trackrate_vspt.pdf");
3491
3492   gPad->SetLogy();
3493   c.SaveAs(DIR + type + "_trackrate_vspt_log.pdf");
3494   gPad->SetLogy(0);
3495
3496   // ---------------------------------------------------------------------------------------------------------
3497   // sum track/ TP pt in jets
3498   /*
3499   if (TP_select_injet > 0) {
3500
3501     TH1F* h_frac_sumpt_vspt = (TH1F*) h_jet_trk_sumpt_vspt->Clone();
3502     h_frac_sumpt_vspt->SetName("frac_sumpt_vspt");
3503     h_frac_sumpt_vspt->GetYaxis()->SetTitle("L1 sum(p_{T}) / TP sum(p_{T})");
3504     h_frac_sumpt_vspt->Divide(h_jet_trk_sumpt_vspt, h_jet_tp_sumpt_vspt, 1.0, 1.0, "B");
3505
3506     TH1F* h_frac_sumpt_vseta = (TH1F*) h_jet_trk_sumpt_vseta->Clone();
3507     h_frac_sumpt_vseta->SetName("frac_sumpt_vseta");
3508     h_frac_sumpt_vseta->GetYaxis()->SetTitle("L1 sum(p_{T}) / TP sum(p_{T})");
3509     h_frac_sumpt_vseta->Divide(h_jet_trk_sumpt_vseta, h_jet_tp_sumpt_vseta, 1.0, 1.0, "B");
3510
3511
3512     TH1F* h_matchfrac_sumpt_vspt = (TH1F*) h_jet_matchtrk_sumpt_vspt->Clone();
3513     h_matchfrac_sumpt_vspt->SetName("matchfrac_sumpt_vspt");
3514     h_matchfrac_sumpt_vspt->GetYaxis()->SetTitle("Matched L1 sum(p_{T}) / TP sum(p_{T})");
3515     h_matchfrac_sumpt_vspt->Divide(h_jet_matchtrk_sumpt_vspt, h_jet_tp_sumpt_vspt, 1.0, 1.0, "B");
3516
3517     TH1F* h_matchfrac_sumpt_vseta = (TH1F*) h_jet_matchtrk_sumpt_vseta->Clone();
3518     h_matchfrac_sumpt_vseta->SetName("matchfrac_sumpt_vseta");
3519     h_matchfrac_sumpt_vseta->GetYaxis()->SetTitle("Matched L1 sum(p_{T}) / TP sum(p_{T})");
3520     h_matchfrac_sumpt_vseta->Divide(h_jet_matchtrk_sumpt_vseta, h_jet_tp_sumpt_vseta, 1.0, 1.0, "B");
3521
3522
3523     h_frac_sumpt_vspt->Draw();
3524     c.SaveAs(DIR+type+"_sumpt_vspt.pdf");
3525
3526     h_frac_sumpt_vseta->Draw();
3527     c.SaveAs(DIR+type+"_sumpt_vseta.pdf");
3528
3529     h_matchfrac_sumpt_vspt->Draw();
3530     c.SaveAs(DIR+type+"_sumpt_match_vspt.pdf");
3531
3532     h_matchfrac_sumpt_vseta->Draw();
3533     c.SaveAs(DIR+type+"_sumpt_match_vseta.pdf");
3534   }
3535   */
3536
3537   // nbr tracks per event
3538
3539   h_ntrk_pt2->Write();
3540   h_ntrk_pt3->Write();
3541   h_ntrk_pt10->Write();
3542
3543   h_ntrkPerSector_all->Write();
3544   h_ntrkPerSector_pt2->Write();
3545   h_ntrkPerSector_pt3->Write();
3546   h_ntrkPerSector_pt4->Write();
3547
3548   h_ntrkPerSector_all->Scale(1.0 / nevt);
3549   h_ntrkPerSector_pt2->Scale(1.0 / nevt);
3550   h_ntrkPerSector_pt3->Scale(1.0 / nevt);
3551   h_ntrkPerSector_pt4->Scale(1.0 / nevt);
3552
3553   h_ntrkPerSector_all->GetYaxis()->SetTitle("Fraction of events");
3554   h_ntrkPerSector_all->GetXaxis()->SetTitle("Max number of transmitted tracks per #phi sector");
3555
3556   h_ntrkPerSector_all->SetLineColor(1);
3557   h_ntrkPerSector_pt2->SetLineColor(4);
3558   h_ntrkPerSector_pt3->SetLineColor(2);
3559   h_ntrkPerSector_pt4->SetLineColor(8);
3560
3561   max = h_ntrkPerSector_all->GetMaximum();
3562   h_ntrkPerSector_all->SetAxisRange(0.00001, max * 5, "Y");
3563   h_ntrkPerSector_all->SetAxisRange(0., 100, "X");
3564
3565   h_ntrkPerSector_all->Draw("hist");
3566   h_ntrkPerSector_pt2->Draw("same,hist");
3567   h_ntrkPerSector_pt3->Draw("same,hist");
3568   h_ntrkPerSector_pt4->Draw("same,hist");
3569   gPad->SetLogy();
3570
3571   TLegend* l = new TLegend(0.6, 0.55, 0.85, 0.85);
3572   l->SetFillStyle(0);
3573   l->SetBorderSize(0);
3574   l->SetTextSize(0.04);
3575   l->AddEntry(h_ntrkPerSector_all, "no p_{T}cut", "l");
3576   l->AddEntry(h_ntrkPerSector_pt2, "p_{T}^{track} > 2 GeV", "l");
3577   l->AddEntry(h_ntrkPerSector_pt3, "p_{T}^{track} > 3 GeV", "l");
3578   l->AddEntry(h_ntrkPerSector_pt4, "p_{T}^{track} > 4 GeV", "l");
3579   l->SetTextFont(42);
3580   l->Draw();
3581
3582   c.SaveAs(DIR + type + "_trackRatePerPhiSector_log.pdf");
3583   gPad->SetLogy(0);
3584
3585   h_ntrk_genuine_pt2->Write();
3586   h_ntrk_genuine_pt3->Write();
3587   h_ntrk_genuine_pt10->Write();
3588
3589   if (doDetailedPlots) {
3590     h_ntrk_pt2->Draw();
3591     c.SaveAs(DIR + type + "_trackrate_pt2_perevt.pdf");
3592
3593     h_ntrk_pt3->Draw();
3594     c.SaveAs(DIR + type + "_trackrate_pt3_perevt.pdf");
3595
3596     h_ntrk_pt10->Draw();
3597     c.SaveAs(DIR + type + "_trackrate_pt10_perevt.pdf");
3598   }
3599
3600   // number of tracks vs. eta, pT (trk_eta/trk_pt)
3601
3602   if (doDetailedPlots) {
3603     h_trk_eta->Write();
3604     h_trk_pt->Write();
3605
3606     h_trk_eta->Draw();
3607     c.SaveAs(DIR + type + "_trk_eta.pdf");
3608     h_trk_pt->Draw();
3609     c.SaveAs(DIR + type + "_trk_pt.pdf");
3610   }
3611
3612   fout->Close();
3613
3614   // ---------------------------------------------------------------------------------------------------------
3615   //some printouts
3616
3617   float k = (float)n_match_eta1p0;
3618   float N = (float)n_all_eta1p0;
3619   if (std::abs(N) > 0)
3620     cout << endl
3621          << "efficiency for |eta| < 1.0 = " << k / N * 100.0 << " +- " << 1.0 / N * sqrt(k * (1.0 - k / N)) * 100.0
3622          << endl;
3623   k = (float)n_match_eta1p75;
3624   N = (float)n_all_eta1p75;
3625   if (std::abs(N) > 0)
3626     cout << "efficiency for 1.0 < |eta| < 1.75 = " << k / N * 100.0 << " +- "
3627          << 1.0 / N * sqrt(k * (1.0 - k / N)) * 100.0 << endl;
3628   k = (float)n_match_eta2p5;
3629   N = (float)n_all_eta2p5;
3630   if (std::abs(N) > 0)
3631     cout << "efficiency for 1.75 < |eta| < " << std::min(TP_maxEta, 2.5f) << " = " << k / N * 100.0 << " +- "
3632          << 1.0 / N * sqrt(k * (1.0 - k / N)) * 100.0 << endl;
3633   N = (float)n_all_eta1p0 + n_all_eta1p75 + n_all_eta2p5;
3634   k = (float)n_match_eta1p0 + n_match_eta1p75 + n_match_eta2p5;
3635   if (std::abs(N) > 0)
3636     cout << "combined efficiency for |eta| < " << std::min(TP_maxEta, 2.5f) << " = " << k / N * 100.0 << " +- "
3637          << 1.0 / N * sqrt(k * (1.0 - k / N)) * 100.0 << " = " << k << "/" << N << endl
3638          << endl;
3639
3640   k = (float)n_match_ptg2;
3641   N = (float)n_all_ptg2;
3642   if (std::abs(N) > 0)
3643     cout << "efficiency for pt > " << std::max(TP_minPt, 2.0f) << " = " << k / N * 100.0 << " +- "
3644          << 1.0 / N * sqrt(k * (1.0 - k / N)) * 100.0 << endl;
3645   k = (float)n_match_pt2to8;
3646   N = (float)n_all_pt2to8;
3647   if (std::abs(N) > 0)
3648     cout << "efficiency for " << std::max(TP_minPt, 2.0f) << " < pt < 8.0 = " << k / N * 100.0 << " +- "
3649          << 1.0 / N * sqrt(k * (1.0 - k / N)) * 100.0 << endl;
3650   k = (float)n_match_ptg8;
3651   N = (float)n_all_ptg8;
3652   if (std::abs(N) > 0)
3653     cout << "efficiency for pt > 8.0 = " << k / N * 100.0 << " +- " << 1.0 / N * sqrt(k * (1.0 - k / N)) * 100.0
3654          << endl;
3655   k = (float)n_match_ptg40;
3656   N = (float)n_all_ptg40;
3657   if (std::abs(N) > 0)
3658     cout << "efficiency for pt > 40.0 = " << k / N * 100.0 << " +- " << 1.0 / N * sqrt(k * (1.0 - k / N)) * 100.0
3659          << endl
3660          << endl;
3661
3662   // track rates
3663   cout << "# TP/event (pt > " << std::max(TP_minPt, 2.0f) << ") = " << (float)ntp_pt2 / nevt << endl;
3664   cout << "# TP/event (pt > 3.0) = " << (float)ntp_pt3 / nevt << endl;
3665   cout << "# TP/event (pt > 10.0) = " << (float)ntp_pt10 / nevt << endl;
3666
3667   cout << "# tracks/event (no pt cut)= " << (float)ntrk / nevt << endl;
3668   cout << "# tracks/event (pt > " << std::max(TP_minPt, 2.0f) << ") = " << (float)ntrk_pt2 / nevt << endl;
3669   cout << "# tracks/event (pt > 3.0) = " << (float)ntrk_pt3 / nevt << endl;
3670   cout << "# tracks/event (pt > 10.0) = " << (float)ntrk_pt10 / nevt << endl;
3671 }
3672
3673 void SetPlotStyle() {
3674   // from ATLAS plot style macro
3675
3676   // use plain black on white colors
3677   gStyle->SetFrameBorderMode(0);
3678   gStyle->SetFrameFillColor(0);
3679   gStyle->SetCanvasBorderMode(0);
3680   gStyle->SetCanvasColor(0);
3681   gStyle->SetPadBorderMode(0);
3682   gStyle->SetPadColor(0);
3683   gStyle->SetStatColor(0);
3684   gStyle->SetHistLineColor(1);
3685
3686   gStyle->SetPalette(1);
3687
3688   // set the paper & margin sizes
3689   gStyle->SetPaperSize(20, 26);
3690   gStyle->SetPadTopMargin(0.05);
3691   gStyle->SetPadRightMargin(0.05);
3692   gStyle->SetPadBottomMargin(0.16);
3693   gStyle->SetPadLeftMargin(0.16);
3694
3695   // set title offsets (for axis label)
3696   gStyle->SetTitleXOffset(1.4);
3697   gStyle->SetTitleYOffset(1.4);
3698
3699   // use large fonts
3700   gStyle->SetTextFont(42);
3701   gStyle->SetTextSize(0.05);
3702   gStyle->SetLabelFont(42, "x");
3703   gStyle->SetTitleFont(42, "x");
3704   gStyle->SetLabelFont(42, "y");
3705   gStyle->SetTitleFont(42, "y");
3706   gStyle->SetLabelFont(42, "z");
3707   gStyle->SetTitleFont(42, "z");
3708   gStyle->SetLabelSize(0.05, "x");
3709   gStyle->SetTitleSize(0.05, "x");
3710   gStyle->SetLabelSize(0.05, "y");
3711   gStyle->SetTitleSize(0.05, "y");
3712   gStyle->SetLabelSize(0.05, "z");
3713   gStyle->SetTitleSize(0.05, "z");
3714
3715   // use bold lines and markers
3716   gStyle->SetMarkerStyle(20);
3717   gStyle->SetMarkerSize(1.2);
3718   gStyle->SetHistLineWidth(2.);
3719   gStyle->SetLineStyleString(2, "[12 12]");
3720
3721   // get rid of error bar caps
3722   gStyle->SetEndErrorSize(0.);
3723
3724   // do not display any of the standard histogram decorations
3725   gStyle->SetOptTitle(0);
3726   gStyle->SetOptStat(0);
3727   gStyle->SetOptFit(0);
3728
3729   // put tick marks on top and RHS of plots
3730   gStyle->SetPadTickX(1);
3731   gStyle->SetPadTickY(1);
3732 }
3733
3734 void mySmallText(Double_t x, Double_t y, Color_t color, char* text) {
3735   Double_t tsize = 0.044;
3736   TLatex l;
3737   l.SetTextSize(tsize);
3738   l.SetNDC();
3739   l.SetTextColor(color);
3740   l.DrawLatex(x, y, text);
3741 }
3742
3743 double getIntervalContainingFractionOfEntries(TH1* absResidualHistogram, double quantileToCalculate, int minEntries) {
3744   double totalIntegral = absResidualHistogram->Integral(0, absResidualHistogram->GetNbinsX() + 1);
3745   double numEntries = absResidualHistogram->GetEntries();
3746
3747   // Check that the interval is not somewhere in the overflow bin
3748   double maxAllowedEntriesInOverflow = totalIntegral * (1 - quantileToCalculate);
3749   double nEntriesInOverflow = absResidualHistogram->GetBinContent(absResidualHistogram->GetNbinsX() + 1);
3750   if (nEntriesInOverflow > maxAllowedEntriesInOverflow) {
3751     // cout << "WARNING : Cannot compute range corresponding to interval, as it is in the overflow bin" << endl;
3752     return absResidualHistogram->GetXaxis()->GetXmax() * 1.2;
3753   }
3754
3755   // Calculate quantile for given interval
3756   double interval[1];
3757   double quantile[1] = {quantileToCalculate};
3758   if (totalIntegral > 0 && numEntries >= minEntries) {
3759     absResidualHistogram->GetQuantiles(1, interval, quantile);
3760   } else {
3761     cout << "WARNING: histo " << absResidualHistogram->GetName()
3762          << " empty or with too few entries, so can't calc quantiles." << endl;
3763     interval[0] = 0.;
3764   }
3765
3766   return interval[0];
3767 }
3768
3769 void makeResidualIntervalPlot(
3770     TString type, TString dir, TString variable, TH1F* h_68, TH1F* h_90, TH1F* h_99, double minY, double maxY) {
3771   TCanvas c;
3772
3773   h_68->SetMinimum(minY);
3774   h_90->SetMinimum(minY);
3775   h_99->SetMinimum(minY);
3776
3777   h_68->SetMaximum(maxY);
3778   h_90->SetMaximum(maxY);
3779   h_99->SetMaximum(maxY);
3780
3781   h_68->SetMarkerStyle(20);
3782   h_90->SetMarkerStyle(26);
3783   h_99->SetMarkerStyle(24);
3784
3785   h_68->Draw("P");
3786   h_68->Write();
3787   h_90->Draw("P same");
3788   h_90->Write();
3789   h_99->Draw("P same");
3790   h_99->Write();
3791
3792   TLegend* l = new TLegend(0.65, 0.65, 0.85, 0.85);
3793   l->SetFillStyle(0);
3794   l->SetBorderSize(0);
3795   l->SetTextSize(0.04);
3796   l->AddEntry(h_99, "99%", "p");
3797   l->AddEntry(h_90, "90%", "p");
3798   l->AddEntry(h_68, "68%", "p");
3799   l->SetTextFont(42);
3800   l->Draw();
3801
3802   c.SaveAs(dir + type + "_" + variable + "_interval.pdf");
3803
3804   delete l;
3805 }