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