HcalTestBeam/macros/PlotTracks.C

0001 ////////////////////////////////////////////////////////////////////////////////
0002 //
0003 // void plotAll(std::string fname, std::string HLT, int var, int ien, int eta,
0004 //              bool varbin, int rebin, bool approve, bool logy, int pos,
0005 //              bool pv, int savePlot)
0006 // Gets a group of plots by calling plotEnergy or plotEnergyPV for a given
0007 // input file
0008 //   fname = name of the input file                             ["hlt.root"]
0009 //   HLT   = type of HLT used (to be given in the figure legend)["All HLTs"]
0010 //   var   = variable name (-1 means all variables)                 [-1]
0011 //   ien   = energy bin (-1 means all energy bins)                  [-1]
0012 //   eta   = eta bin (-1 means all eta bins)                        [-1]
0013 //   varbin= flag for using variable bin width                      [false]
0014 //   rebin = # of bins to be re-binned together                     [5]
0015 //   approve= If meant for approval talks                           [true]
0016 //   logy  = If y-axis scale shuld by linear/logarithmic            [true]
0017 //   pos   = position of the statistics boxes                       [0]
0018 //   pv    = flag deciding call to plotEnergyPV vs plotEnergy       [false]
0019 //   savePlot= Plot to be saved: no (-1), eps (0), gif (1), pdf (2) [-1]
0020 //
0021 // void plotEnergyAll(std::string fname, std::string hlt, int pv, int data,
0022 //                    bool varbin, int rebin, bool approve, bool logy, int pos,
0023 //                    int var, int ene, int eta, int savePlot)
0024 // Plots energy distributions for a number of variables, eta and energy bins
0025 //
0026 //   fname = name of the i/p root file                              [""]
0027 //   HLT   = type of HLT used (to be given in the figure legend)  ["All HLTs"]
0028 //   models= packed flag to select which files to be plotted        [15]
0029 //   pv    = selection of # of good PV's used                       [0]
0030 //   data  = flag to say if it is Data (1)/MC (2)/or comparison (4) [4]
0031 //   varbin= flag for using variable bin width                      [false]
0032 //   rebin = # of bins to be re-binned together                     [5]
0033 //   approve= If meant for approval talks                           [true]
0034 //   logy  = If y-axis scale shuld by linear/logarithmic            [true]
0035 //   pos   = position of the statistics boxes                       [0]
0036 //   var   = variable name (-1 means all variables)                 [-1]
0037 //   ene   = energy bin (-1 means all energy bins)                  [-1]
0038 //   eta   = eta bin (-1 means all eta bins)                        [-1]
0039 //   savePlot= Plot to be saved: no (-1), eps (0), gif (1), pdf (2) [-1]
0040 //
0041 // void plotEMeanAll(int data, int models, bool ratio, bool approve,
0042 //                   std::string postfix, int savePlot)
0043 //
0044 // Plots mean energy response as a function of track momentum or the ratio
0045 //
0046 //   data  = flag to say if it is Data (1)/MC (2)/or comparison (4) [4]
0047 //   models= packed flag to select which files to be plotted        [15]
0048 //   ratio = flag to say if raw mean will be shown or ratio         [false]
0049 //           wrt a reference
0050 //   approve= If meant for approval talks                           [true]
0051 //   postfix= String to be added in the name of saved file          [""]
0052 //   savePlot= Plot to be saved: no (-1), eps (0), gif (1), pdf (2) [-1]
0053 //
0054 // void plotEMean(std::string fname, std::string hlt, int models, int var,
0055 //                int eta, int pv, int data, bool ratio, bool approve,
0056 //                std::string postfix, int savePlot)
0057 //
0058 // Plots mean energy response as a function of track momentum or the ratio
0059 // MC/Data
0060 //
0061 //   fname = name of the i/p root file                              [""]
0062 //   HLT   = type of HLT used (to be given in the figure legend)   ["All HLTs"]
0063 //   models= packed flag to select which files to be plotted        [15]
0064 //   var   = type of energy reponse with values between 0:5         [0]
0065 //           E_{7x7}/p, H_{3x3}/p, (E_{7x7}+H_{3x3})/p,
0066 //           E_{11x11}/p, H_{5x5}/p, (E_{11x11}+H_{5x5})/p
0067 //   eta   = calorimeter cell where track will reach 0:3            [0]
0068 //           ieta (HCAL) values 1:7, 7-13, 13:17, 17:23
0069 //   pv    = selection of # of good PV's used                       [0]
0070 //   data  = flag to say if it is Data (1)/MC (2)/or comparison (4) [1]
0071 //   ratio = flag to say if raw mean will be shown or ratio wrt     [false]
0072 //           a reference
0073 //   approve= If meant for approval talks                           [true]
0074 //   postfix= String to be added in the name of saved file          [""]
0075 //   savePlot= Plot to be saved: no (-1), eps (0), gif (1), pdf (2) [-1]
0076 //
0077 // void plotEnergy(std::string fname, std::string HLT, int var, int ien,
0078 //                 int eta, bool varbin, int rebin, bool approve, bool logy,
0079 //                 int pos, int coloff)
0080 //
0081 // Plots energy response distribution measured energy/track momentum for tracks
0082 // of given momentum range in a eta window
0083 //
0084 //   fname = name of the i/p root file                            ["hlt.root"]
0085 //   HLT   = type of HLT used (to be given in the figure legend)  ["All HLTs"]
0086 //   var   = type of energy reponse with values between 0:5       [0]
0087 //           E_{7x7}/p, H_{3x3}/p, (E_{7x7}+H_{3x3})/p,
0088 //           E_{11x11}/p, H_{5x5}/p, (E_{11x11}+H_{5x5})/p
0089 //   ien   = Track momentum range with values between 0:9         [0]
0090 //           1:2,2:3,3:4,4:5,5:6,6:7,7:9,9:11,11:15,15:20
0091 //   eta   = calorimeter cell where track will reach 0:3          [0]
0092 //           ieta (HCAL) values 1:7, 7-13, 13:17, 17:23
0093 //   varbin= flag for using variable bin width                    [false]
0094 //   rebin = # of bins to be re-binned together                   [1]
0095 //   approve= If meant for approval talks                         [false]
0096 //   logy  = If y-axis scale shuld by linear/logarithmic          [true]
0097 //   pos   = position of the statistics boxes                     [0]
0098 //   coloff= color offset                                         [0]
0099 //
0100 // TCanvas* plotEnergyPV(std::string fname, std::string HLT, int var, int ien,
0101 //                       int eta, bool varbin, int rebin, bool approve,
0102 //                       bool logy, int pos)
0103 //
0104 // Plots energy response distribution measured energy/track momentum for tracks
0105 // of given momentum range in a eta window for 4 different selections of # of
0106 // primary vertex 1:2,2:3,3:5,5:100
0107 //
0108 //   fname = name of the i/p root file            ["StudyHLT_HLTZeroBias.root"]
0109 //   HLT   = type of HLT used (to be given in the figure legend) ["Zero Bias"]
0110 //   var   = type of energy reponse with values between 0:5      [0]
0111 //           E_{7x7}/p, H_{3x3}/p, (E_{7x7}+H_{3x3})/p,
0112 //           E_{11x11}/p, H_{5x5}/p, (E_{11x11}+H_{5x5})/p
0113 //   ien   = Track momentum range with values between 0:9        [0]
0114 //           1:2,2:3,3:4,4:5,5:6,6:7,7:9,9:11,11:15,15:20
0115 //   eta   = calorimeter cell where track will reach 0:3         [0]
0116 //           ieta (HCAL) values 1:7, 7-13, 13:17, 17:23
0117 //   varbin= flag for using variable bin width                   [false]
0118 //   rebin = # of bins to be re-binned together                  [1]
0119 //   approve= If meant for approval talks                        [false]
0120 //   logy  = If y-axis scale shuld by linear/logarithmic         [true]
0121 //   pos   = position of the statistics boxes                    [0]
0122 //
0123 // void plotTrack(std::string fname, std::string HLT, int var, bool varbin,
0124 //                int rebin, bool approve, bool logy, int pos)
0125 //
0126 // Plots kinematic propeties of the track
0127 //
0128 //   fname = name of the i/p root file                           ["hlt.root"]
0129 //   HLT   = type of HLT used (to be given in the figure legend  ["All HLTs"]
0130 //   var   = kinematic variable 0:3 --> p, pt, eta, phi          [0]
0131 //   varbin= flag for using variable bin width                   [false]
0132 //   rebin = # of bins to be re-binned together                  [1]
0133 //   approve= If meant for approval talks                        [false]
0134 //   logy  = If y-axis scale shuld by linear/logarithmic         [true]
0135 //   pos   = position of the statistics boxes                    [0]
0136 //
0137 // void plotIsolation(std::string fname, std::string HLT, int var, bool varbin,
0138 //                    int rebin, bool approve, bool logy, int pos)
0139 //
0140 // Plots variables used for deciding on track isolation
0141 //
0142 //   fname = name of the i/p root file                           ["hlt.root"]
0143 //   HLT   = type of HLT used (to be given in the figure legend  ["All HLTs"]
0144 //   var   = isolation variable 0:3 --> Charge isolation energy, [0]
0145 //     Neutral isolation energy, Energy in smaller cone,
0146 //     Energy in larger cone
0147 //   varbin= flag for using variable bin width                   [false]
0148 //   rebin = # of bins to be re-binned together                  [1]
0149 //   approve= If meant for approval talks                        [false]
0150 //   logy  = If y-axis scale shuld by linear/logarithmic         [true]
0151 //   pos   = position of the statistics boxes                    [0]
0152 //
0153 // void plotHLT(std::string fname, std::string HLT, int run, bool varbin,
0154 //                    int rebin, bool approve, bool logy, int pos)
0155 //
0156 // Plots HLT accept information for a given run or a summary
0157 //
0158 //   fname = name of the i/p root file                           ["hlt.root"]
0159 //   HLT   = type of HLT used (to be given in the figure legend  ["All HLTs"]
0160 //   run   = run number; if <=0 the overall summary              [-1]
0161 //   varbin= flag for using variable bin width                   [false]
0162 //   rebin = # of bins to be re-binned together                  [1]
0163 //   approve= If meant for approval talks                        [false]
0164 //   logy  = If y-axis scale shuld by linear/logarithmic         [true]
0165 //   pos   = position of the statistics boxes                    [0]
0166 ////////////////////////////////////////////////////////////////////////////////
0167
0168 #include "TCanvas.h"
0169 #include "TDirectory.h"
0170 #include "TF1.h"
0171 #include "TFile.h"
0172 #include "TFitResult.h"
0173 #include "TGraph.h"
0174 #include "TGraphAsymmErrors.h"
0175 #include "TH1D.h"
0176 #include "TH2D.h"
0177 #include "THStack.h"
0178 #include "TLegend.h"
0179 #include "TMath.h"
0180 #include "TProfile.h"
0181 #include "TPaveStats.h"
0182 #include "TPaveText.h"
0183 #include "TROOT.h"
0184 #include "TString.h"
0185 #include "TStyle.h"
0186
0187 #include <iostream>
0188 #include <iomanip>
0189 #include <vector>
0190 #include <string>
0191
0192 //const int nmodelm=15, nmodelx=16, nmodels=2;
0193 //const int nmodelm=3, nmodelx=4, nmodels=3;
0194 //const int nmodelm=4, nmodelx=5, nmodels=3;
0195 //const int nmodelm=4, nmodelx=5, nmodels=2;
0196 //const int nmodelm=2, nmodelx=3, nmodels=2;
0197 const int nmodelm=5, nmodelx=6, nmodels=3;
0198 int         styles[7]  = {20, 21, 24, 22, 23, 25, 33};
0199 int         colors[7]  = {1, 2, 4, 6, 7, 38, 3};
0200 std::string names[7]   = {"All", "Quality", "okEcal", "EcalCharIso",
0201               "HcalCharIso", "EcalNeutIso", "HcalNeutIso"};
0202 std::string namefull[7]= {"All tracks", "Good quality tracks",
0203               "Tracks reaching ECAL", "Charge isolation in ECAL",
0204               "Charge isolation in HCAL", "Isolated in ECAL",
0205               "Isolated in HCAL"};
0206 std::string nameEta[4] = {"i#eta 1:6","i#eta 7:12","i#eta 13:16","i#eta 17:22"};
0207 std::string nameEtas[4]= {"|#eta| < 0.52", "0.52 < |#eta| < 1.04",
0208               "1.04 < |#eta| < 1.39", "1.39 < |#eta| < 2.01"};
0209 std::string namePV[5]  = {"all PV","PV 1:1","PV 2:2","PV 3:5","PV > 5"};
0210 std::string varname[4] = {"p", "pt", "eta", "phi"};
0211 std::string vartitle[4]= {"p (GeV/c)", "p_{T} (GeV/c)", "#eta", "#phi"};
0212 std::string nameC[2]   = {"Ecal", "Hcal"};
0213 std::string nameCF[2]  = {"ECAL", "HCAL"};
0214 std::string varnameC[4]= {"maxNearP", "ediff", "ene1", "ene2"};
0215 std::string vartitlC[4]= {"Charge isolation energy",
0216               "Neutral isolation energy",
0217               "Energy in smaller cone",
0218               "Energy in larger cone"};
0219 const int NPT=10;
0220 double mom[NPT]={1.5,2.5,3.5,4.5,5.5,6.5,8.0,10.0,13.0,17.5};
0221 double dmom[NPT]={0.5,0.5,0.5,0.5,0.5,0.5,1.0,1.0,2.0,2.5};
0222 std::string varPs[NPT] = {"1:2","2:3","3:4","4:5","5:6","6:7","7:9",
0223               "9:11","11:15","15:20"};
0224 std::string varPs1[NPT] = {"1","2","3","4","5","6","7","9","11","15"};
0225 std::string varPPs[NPT] = {"1-2 GeV","2-3 GeV","3-4 GeV","4-5 GeV","5-6 GeV",
0226                "6-7 GeV","7-9 GeV","9-11 GeV","11-15 GeV",
0227                "15-20 GeV"};
0228 std::string varEta[4] = {"1:6", "7:12", "13:16", "17:23"};
0229 std::string varEta1[4]= {"1", "2", "3", "4"};
0230 std::string varEne[6] = {"E_{7x7}", "H_{3x3}", "(E_{7x7}+H_{3x3})",
0231              "E_{11x11}", "H_{5x5}", "(E_{11x11}+H_{5x5})"};
0232 std::string varEne1[6]= {"E7x7","H3x3","E7x7H3x3","E11x11","H5x5", "E11x11H5x5"};
0233 const int nbins=100;
0234 double xbins[nbins+1] = {0.00,0.01,0.02,0.03,0.04,0.05,0.06,0.07,0.08,0.09,
0235              0.10,0.11,0.12,0.13,0.14,0.15,0.16,0.17,0.18,0.19,
0236              0.20,0.21,0.22,0.23,0.24,0.25,0.26,0.27,0.28,0.29,
0237              0.30,0.31,0.32,0.33,0.34,0.35,0.36,0.37,0.38,0.39,
0238              0.40,0.41,0.42,0.43,0.44,0.45,0.46,0.47,0.48,0.49,
0239              0.50,0.52,0.54,0.56,0.58,0.60,0.62,0.64,0.66,0.68,
0240              0.70,0.72,0.74,0.76,0.78,0.80,0.82,0.84,0.86,0.88,
0241              0.90,0.92,0.94,0.96,0.98,1.00,1.05,1.10,1.15,1.20,
0242              1.25,1.30,1.35,1.40,1.45,1.50,1.55,1.60,1.65,1.70,
0243              1.75,1.80,1.85,1.90,1.95,2.00,2.10,2.20,2.30,2.40,
0244              2.50};
0245 int ibins[nbins+1] = { 11, 12, 13, 14, 15, 16, 17, 18, 19, 20,
0246                21, 22, 23, 24, 25, 26, 27, 28, 29, 30,
0247                31, 32, 33, 34, 35, 36, 37, 38, 39, 40,
0248                41, 42, 43, 44, 45, 46, 47, 48, 49, 50,
0249                51, 52, 53, 54, 55, 56, 57, 58, 59, 60,
0250                61, 63, 65, 67, 69, 71, 73, 75, 77, 79,
0251                81, 83, 85, 87, 89, 91, 93, 95, 97, 99,
0252               101,103,105,107,109,111,116,121,126,131,
0253               136,141,146,151,156,161,166,171,176,181,
0254               186,191,196,201,206,211,221,231,241,251,261};
0255 /*
0256 std::string files[nmodels]={"ZeroBiasStudyHLT.root","MinimumBiasStudyHLT.root"};
0257 std::string types[nmodels]={"Zero Bias Data","Minimum Bias Data"};
0258 */
0259 std::string files[nmodels]={"2017C_ZB.root","2017G_ZB.root","2017H_ZB.root"};
0260 std::string types[nmodels]={"Zero Bias (2017C)","Zero Bias (2017G)",
0261                 "Zero Bias (2017H)"};
0262
0263 /*
0264 std::string filem[nmodelm]={"pikp/FBE2p2StudyHLT.root",
0265                 "pikp/FBE3p3MixStudyHLT.root",
0266                 "pikp/FBE4bMixStudyHLT.root",
0267                 "pikp/FBE4r00MixStudyHLT.root",
0268                 "pikp/FBE4c01MixStudyHLT.root",
0269                 "pikp/FBE4c02MixStudyHLT.root",
0270                 "pikp/FBE3r11MixStudyHLT.root",
0271                 "pikp/FBE3r10MixStudyHLT.root",
0272                 "pikp/FBE3r8MixStudyHLT.root",
0273                 "pikp/FBE3r9MixStudyHLT.root",
0274                 "pikp/QFBE0p2StudyHLT.root",
0275                 "pikp/QFBE2p2StudyHLT.root",
0276                 "pikp/QFBE4bMixStudyHLT.root",
0277                 "pikp/FBAE2p2StudyHLT.root",
0278                 "pikp/FBAE4bMixStudyHLT.root"};
0279 std::string typem[nmodelm]={"10.2.p02 FTFP_BERT_EMM",
0280                 "10.3.p03 FTFP_BERT_EMM",
0281                 "10.4.beta FTFP_BERT_EMM",
0282                 "10.4 FTFP_BERT_EMM",
0283                 "10.4.cand01 FTFP_BERT_EMM",
0284                 "10.4.cand02 FTFP_BERT_EMM",
0285                 "10.3.ref11 FTFP_BERT_EMM",
0286                 "10.3.ref10 FTFP_BERT_EMM",
0287                 "10.3.ref08 FTFP_BERT_EMM",
0288                 "10.3.ref09 FTFP_BERT_EMM",
0289                 "10.0.p02 QGSP_FTFP_BERT_EML",
0290                 "10.2.p02 QGSP_FTFP_BERT_EMM",
0291                 "10.4.b01 QGSP_FTFP_BERT_EMM",
0292                 "10.2.p02 FTFP_BERT_ATL_EMM",
0293                 "10.4.b01 FTFP_BERT_ATL_EMM"};
0294 */
0295 /*
0296 std::string filem[nmodelm]={"pikp/FBE2p2MixStudyHLT10.root",
0297                 "pikp/FBE4MixStudyHLT10.root",
0298                 "pikp/FBE4vMixStudyHLT10.root",
0299                 "pikp/FBE4vcMixStudyHLT10.root"};
0300 std::string typem[nmodelm]={"10.2.p02 FTFP_BERT_EMM",
0301                 "10.4 FTFP_BERT_EMM",
0302                 "10.4 VecGeom FTFP_BERT_EMM",
0303                 "10.4 VecGeom+CLHEP FTFP_BERT_EMM"};
0304 */
0305 /*
0306 std::string filem[nmodelm]={"pikp/FBE4cMixStudyHLT10.root",
0307                 "pikp/FBE4vcMixStudyHLT10.root"};
0308 std::string typem[nmodelm]={"10.4 CLHEP FTFP_BERT_EMM",
0309                 "10.4 VecGeom+CLHEP FTFP_BERT_EMM"};
0310 */
0311 /*
0312 std::string filem[nmodelm]={"pikp/FBE3r9MixStudyHLT.root",
0313                 "pikp/FBE4r00MixStudyHLT.root",
0314                 "pikp/FBE4r01MixStudyHLT.root"};
0315 std::string typem[nmodelm]={"10.3.ref09 FTFP_BERT_EMM",
0316                 "10.4.ref00 FTFP_BERT_EMM",
0317                 "10.4.ref01 FTFP_BERT_EMM"};
0318 */
0319 /*
0320 std::string filem[nmodelm]={"pikp/QFB3r9MixStudyHLT.root",
0321                 "pikp/QFB4r00MixStudyHLT.root",
0322                 "pikp/QFB4r01MixStudyHLT.root"};
0323 std::string typem[nmodelm]={"10.3.ref09 QGSP_FTFP_BERT_EML",
0324                 "10.4.ref00 QGSP_FTFP_BERT_EML",
0325                 "10.4.ref01 QGSP_FTFP_BERT_EML"};
0326 */
0327 /*
0328 std::string filem[nmodelm]={"pikp/FBE2p2StudyHLT.root",
0329                 "pikp/FBE4r00MixStudyHLT.root",
0330                 "pikp/FBE4r00vMixStudyHLT.root"};
0331 std::string typem[nmodelm]={"10.2.p02 FTFP_BERT_EMM",
0332                 "10.4 FTFP_BERT_EMM (Native)",
0333                 "10.4 FTFP_BERT_EMM (VecGeom)"};
0334 */
0335 /*
0336 std::string filem[nmodelm]={"pikp/FBE3r6MixStudyHLT.root",
0337                 "pikp/FBE3r6vMixStudyHLT.root",
0338                 "pikp/FBE3r6vRMixStudyHLT.root"};
0339 std::string typem[nmodelm]={"10.3.ref06 FTFP_BERT_EMM (Native)",
0340                 "10.3.ref06 FTFP_BERT_EMM (VecGeom 4)",
0341                 "10.3.ref06 FTFP_BERT_EMM (VecGeom Corr)"};
0342 */
0343 /*
0344 std::string filem[nmodelm]={"pikp/FBE4r00vMixStudyHLT.root",
0345                 "pikp/FBE4r01MixStudyHLT.root",
0346                 "pikp/FBE4r01vMixStudyHLT.root"};
0347 std::string typem[nmodelm]={"10.4 VecGeom FTFP_BERT_EMM",
0348                 "10.4.ref01 FTFP_BERT_EMM",
0349                 "10.4.ref01 VecGeom FTFP_BERT_EMM"};
0350 */
0351 /*
0352 std::string filem[nmodelm]={"pikp/FBE4MixStudyHLT10d.root",
0353                 "pikp/FBE4MixStudyHLT10s.root",
0354                 "pikp/FBE4MixStudyHLT10t.root",
0355                 "pikp/FBE4MixStudyHLT10st.root"};
0356 std::string typem[nmodelm]={"10.4 FTFP_BERT_EMM (Default)",
0357                 "10.4 FTFP_BERT_EMM (New StepIntegrator)",
0358                 "10.4 FTFP_BERT_EMM (Smart Track)",
0359                 "10.4 FTFP_BERT_EMM (New Stepper + Smart Track)"};
0360 */
0361 /*
0362 std::string filem[nmodelm]={"pikp/FBE2p2StudyHLT.root",
0363                 "pikp/FBE4MixStudyHLT1025.root",
0364                 "pikp/FBE5c01vMixStudyHLT.root",
0365                 "pikp/FBE5c02vMixStudyHLT.root",
0366                 "pikp/QFB5c02vMixStudyHLT.root"};
0367 std::string typem[nmodelm]={"10.2.p02 FTFP_BERT_EMM",
0368                 "10.4 FTFP_BERT_EMM",
0369                 "10.5.cand01 FTFP_BERT_EMM",
0370                 "10.5.cand02 FTFP_BERT_EMM",
0371                 "10.5.cand02 QGSP_FTFP_BERT_EML"};
0372 */
0373 /*
0374 std::string filem[nmodelm]={"pikp/FBE4p0vMixStudyHLT.root",
0375                 "pikp/FBE4p3vMixStudyHLT.root",
0376                 "pikp/FBE5r00vMixStudyHLT.root",
0377                 "pikp/FBE5r05vMixStudyHLT.root",
0378                 "pikp/FBE5r06vMixStudyHLT.root"};
0379 std::string typem[nmodelm]={"10.4  FTFP_BERT_EMM",
0380                 "10.4.p03  FTFP_BERT_EMM",
0381                 "10.5  FTFP_BERT_EMM",
0382                 "10.5.ref05 FTFP_BERT_EMM",
0383                 "10.5.ref06 FTFP_BERT_EMM"};
0384 */
0385 std::string filem[nmodelm]={"pikp/QFB4p0vMixStudyHLT.root",
0386                 "pikp/QFB4p3vMixStudyHLT.root",
0387                 "pikp/QFB5r00vMixStudyHLT.root",
0388                 "pikp/QFB5r05vMixStudyHLT.root",
0389                 "pikp/QFB5r06vMixStudyHLT.root"};
0390 std::string typem[nmodelm]={"10.4  QGSP_FTFP_BERT_EML",
0391                 "10.4.p03  QGSP_FTFP_BERT_EML",
0392                 "10.5  QGSP_FTFP_BERT_EML",
0393                 "10.5.ref05 QGSP_FTFP_BERT_EML",
0394                 "10.5.ref06 QGSP_FTFP_BERT_EML"};
0395 /*
0396 std::string filex[nmodelx]={"AllDataStudyHLT.root",
0397                 "pikp/FBE2p2StudyHLT.root",
0398                 "pikp/FBE3p3MixStudyHLT.root",
0399                 "pikp/FBE4bMixStudyHLT.root",
0400                 "pikp/FBE4r00MixStudyHLT.root",
0401                 "pikp/FBE4c01MixStudyHLT.root",
0402                 "pikp/FBE4c02MixStudyHLT.root",
0403                 "pikp/FBE3r11MixStudyHLT.root",
0404                 "pikp/FBE3r10MixStudyHLT.root",
0405                 "pikp/FBE3r8MixStudyHLT.root",
0406                 "pikp/FBE3r9MixStudyHLT.root",
0407                 "pikp/QFBE0p2StudyHLT.root",
0408                 "pikp/QFBE2p2StudyHLT.root",
0409                 "pikp/QFBE4bMixStudyHLT.root",
0410                 "pikp/FBAE2p2StudyHLT.root",
0411                 "pikp/FBAE4bMixStudyHLT.root"};
0412 std::string typex[nmodelx]={"Data (2016B)",
0413                 "10.2.p02 FTFP_BERT_EMM",
0414                 "10.3.p03 FTFP_BERT_EMM",
0415                 "10.4.beta FTFP_BERT_EMM",
0416                 "10.4 FTFP_BERT_EMM",
0417                 "10.4.cand01 FTFP_BERT_EMM",
0418                 "10.4.cand02 FTFP_BERT_EMM",
0419                 "10.3.ref11 FTFP_BERT_EMM",
0420                 "10.3.ref10 FTFP_BERT_EMM",
0421                 "10.3.ref08 FTFP_BERT_EMM",
0422                 "10.3.ref09 FTFP_BERT_EMM",
0423                 "10.0.p02 QGSP_FTFP_BERT_EML",
0424                 "10.2.p02 QGSP_FTFP_BERT_EMM",
0425                 "10.4.b01 QGSP_FTFP_BERT_EMM",
0426                 "10.2.p02 FTFP_BERT_ATL_EMM",
0427                 "10.4.b01 FTFP_BERT_ATL_EMM"};
0428 */
0429 /*
0430 std::string filex[nmodelx]={"AllDataStudyHLT.root",
0431                 "pikp/FBE2p2MixStudyHLT10.root",
0432                 "pikp/FBE4MixStudyHLT10.root",
0433                 "pikp/FBE4vMixStudyHLT10.root",
0434                 "pikp/FBE4vcMixStudyHLT10.root"};
0435 std::string typex[nmodelx]={"Data (2016B)",
0436                 "10.2.p02 FTFP_BERT_EMM",
0437                 "10.4 FTFP_BERT_EMM",
0438                 "10.4 VecGeom FTFP_BERT_EMM",
0439                 "10.4 VecGeom+CLHEP FTFP_BERT_EMM"};
0440 */
0441 /*
0442 std::string filex[nmodelx]={"AllDataStudyHLT.root",
0443                             "pikp/FBE4cMixStudyHLT10.root",
0444                 "pikp/FBE4vcMixStudyHLT10.root"};
0445 std::string typex[nmodelx]={"Data (2016B)",
0446                             "10.4 CLHEP FTFP_BERT_EMM",
0447                 "10.4 VecGeom+CLHEP FTFP_BERT_EMM"};
0448 */
0449 /*
0450 std::string filex[nmodelx]={"AllDataStudyHLT.root",
0451                 "pikp/FBE3r9MixStudyHLT.root",
0452                 "pikp/FBE4r00MixStudyHLT.root",
0453                 "pikp/FBE4r01MixStudyHLT.root"};
0454 std::string typex[nmodelx]={"Data (2016B)",
0455                 "10.3.ref09 FTFP_BERT_EMM",
0456                 "10.4.ref00 FTFP_BERT_EMM",
0457                 "10.4.ref01 FTFP_BERT_EMM"};
0458 */
0459 /*
0460 std::string filex[nmodelx]={"AllDataStudyHLT.root",
0461                 "pikp/QFB3r9MixStudyHLT.root",
0462                 "pikp/QFB4r00MixStudyHLT.root",
0463                 "pikp/QFB4r01MixStudyHLT.root"};
0464 std::string typex[nmodelx]={"Data (2016B)",
0465                 "10.3.ref09 QGSP_FTFP_BERT_EML",
0466                 "10.4.ref00 QGSP_FTFP_BERT_EML",
0467                 "10.4.ref01 QGSP_FTFP_BERT_EML"};
0468 */
0469 /*
0470 std::string filex[nmodelx]={"AllDataStudyHLT.root",
0471                 "pikp/FBE2p2StudyHLT.root",
0472                 "pikp/FBE4r00MixStudyHLT.root",
0473                 "pikp/FBE4r00vMixStudyHLT.root"};
0474 std::string typex[nmodelx]={"Data (2016B)",
0475                 "10.2.p02 FTFP_BERT_EMM",
0476                 "10.4 FTFP_BERT_EMM (Native)",
0477                 "10.4 FTFP_BERT_EMM (VecGeom)"};
0478 */
0479 /*
0480 std::string filex[nmodelx]={"AllDataStudyHLT.root",
0481                 "pikp/FBE3r6MixStudyHLT.root",
0482                 "pikp/FBE3r6vMixStudyHLT.root",
0483                 "pikp/FBE3r6vRMixStudyHLT.root"};
0484 std::string typex[nmodelx]={"Data (2016B)",
0485                 "10.3.ref06 FTFP_BERT_EMM (Native)",
0486                 "10.3.ref06 FTFP_BERT_EMM (VecGeom 4)",
0487                 "10.3.ref06 FTFP_BERT_EMM (VecGeom Corr)"};
0488 */
0489 /*
0490 std::string filex[nmodelx]={"AllDataStudyHLT.root",
0491                 "pikp/FBE4r00vMixStudyHLT.root",
0492                 "pikp/FBE4r01MixStudyHLT.root",
0493                 "pikp/FBE4r01vMixStudyHLT.root"};
0494 std::string typex[nmodelx]={"Data (2016B)",
0495                 "10.4 VecGeom FTFP_BERT_EMM",
0496                 "10.4.ref01 FTFP_BERT_EMM",
0497                 "10.4.ref01 VecGeom FTFP_BERT_EMM"};
0498 */
0499 /*
0500 std::string filex[nmodelx]={"AllDataStudyHLT.root",
0501                 "pikp/FBE4MixStudyHLT10d.root",
0502                 "pikp/FBE4MixStudyHLT10s.root",
0503                 "pikp/FBE4MixStudyHLT10t.root",
0504                 "pikp/FBE4MixStudyHLT10st.root"};
0505 std::string typex[nmodelx]={"Data (2016B)",
0506                 "10.4 FTFP_BERT_EMM (Default)",
0507                 "10.4 FTFP_BERT_EMM (New StepIntegrator)",
0508                 "10.4 FTFP_BERT_EMM (Smart Track)",
0509                 "10.4 FTFP_BERT_EMM (New Stepper + Smart Track)"};
0510 */
0511 /*
0512 std::string filex[nmodelx]={"AllDataStudyHLT.root",
0513                 "pikp/FBE2p2StudyHLT.root",
0514                 "pikp/FBE4MixStudyHLT1025.root",
0515                 "pikp/FBE5c01vMixStudyHLT.root",
0516                 "pikp/FBE5c02vMixStudyHLT.root",
0517                 "pikp/QFB5c02vMixStudyHLT.root"};
0518 std::string typex[nmodelx]={"Data (2016B)",
0519                 "10.2.p02 FTFP_BERT_EMM",
0520                 "10.4 FTFP_BERT_EMM",
0521                 "10.5.cand01 FTFP_BERT_EMM",
0522                 "10.5.cand02 FTFP_BERT_EMM",
0523                 "10.5.cand02 QGSP_FTFP_BERT_EML"};
0524 */
0525 /*
0526 std::string filex[nmodelx]={"AllDataStudyHLT.root",
0527                 "pikp/FBE4p0vMixStudyHLT.root",
0528                 "pikp/FBE4p3vMixStudyHLT.root",
0529                 "pikp/FBE5r00vMixStudyHLT.root",
0530                 "pikp/FBE5r05vMixStudyHLT.root",
0531                 "pikp/FBE5r06vMixStudyHLT.root"};
0532 std::string typex[nmodelx]={"Data (2016B)",
0533                 "10.4  FTFP_BERT_EMM",
0534                 "10.4.p03  FTFP_BERT_EMM",
0535                 "10.5  FTFP_BERT_EMM",
0536                 "10.5.ref05 FTFP_BERT_EMM",
0537                 "10.5.ref06 FTFP_BERT_EMM"};
0538 */
0539 std::string filex[nmodelx]={"AllDataStudyHLT.root",
0540                             "pikp/QFB4p0vMixStudyHLT.root",
0541                 "pikp/QFB4p3vMixStudyHLT.root",
0542                 "pikp/QFB5r00vMixStudyHLT.root",
0543                 "pikp/QFB5r05vMixStudyHLT.root",
0544                 "pikp/QFB5r06vMixStudyHLT.root"};
0545 std::string typex[nmodelx]={"Data (2016B)",
0546                 "10.4  QGSP_FTFP_BERT_EML",
0547                 "10.4.p03  QGSP_FTFP_BERT_EML",
0548                 "10.5  QGSP_FTFP_BERT_EML",
0549                 "10.5.ref05 QGSP_FTFP_BERT_EML",
0550                 "10.5.ref06 QGSP_FTFP_BERT_EML"};
0551 /*
0552   std::string files[nmodelx]={"StudyHLT_ZeroBias_1PV.root","StudyHLT_PixelTrack_1PV.root","StudyHLT_1PV.root"};
0553   std::string types[nmodels]={"Zero Bias HLT","Pixel Track HLT","All HLTs"};
0554   std::string files[nmodels]={"StudyHLT_HLTZeroBias.root","StudyHLT_PixelTrack.root","StudyHLT_HLTJetE.root","StudyHLT_HLTPhysics.root","StudyHLT_All.root"};
0555   std::string types[nmodels]={"Zero Bias HLT","Pixel Track HLT","JetE HLT","Physics HLT","All HLTs"};
0556   std::string filem[nmodelm]={"StudyHLT_95p02_QGSP_FTFP_BERT.root", "StudyHLT_96p02_QGSP_FTFP_BERT.root", "StudyHLT_96p02_FTFP_BERT.root"};
0557   std::string typem[nmodelm]={"Pythia8 (9.5.p02 QGSP_FTFP_BERT)","Pythia8 (9.6.p02 QGSP_FTFP_BERT)","Pythia8 (9.6.p02 FTFP_BERT)"};
0558   std::string filex[nmodelx]={"AllDataStudyHLT.root","pikp/QFBE0p2StudyHLT.root", "pikp/FBE2p2StudyHLT.root"};
0559   std::string typex[nmodelx]={"Data (2016B)","10.0.p02 QGSP_FTFP_BERT_EML","10.2.p02 FTFP_BERT_EMM"};
0560 */
0561
0562 void plotAll(std::string fname="hlt.root", std::string HLT="All HLTs", int var=-1, int ien=-1, int eta=-1, bool varbin=false, int rebin=1, bool approve=true, bool logy=true, int pos=0, bool pv=false, int savePlot=-1);
0563 void plotEnergyAll(std::string fname="", std::string hlt="All HLTs", int models=15, int pv=0, int data=4, bool varbin=false, int rebin=5, bool approve=true, bool logy=true, int pos=0, int var=-1, int ene=-1, int eta=-1, int savePlot=-1);
0564 void plotEMeanAll(int data=4, int models=63, bool ratio=true, bool approve=true, std::string postfix="", int savePlot=-1);
0565 void plotEMean(std::string fname="", std::string hlt="All HLTs", int models=15, int var=0, int eta=0, int pv=0, int dataMC=1, bool raio=false, bool approve=true, std::string postfix="",int savePlot=-1);
0566 TCanvas* plotEMeanDraw(std::vector<std::string> fnames, std::vector<std::string> hlts, int var, int eta, int pv=0, bool approve=false, std::string dtype="Data", int coloff=0);
0567 TCanvas* plotEMeanRatioDraw(std::vector<std::string> fnames, std::vector<std::string> hlts, int var, int eta, int pv=0, bool approve=false, std::string dtype="Data", int coloff=0);
0568 TCanvas* plotEnergies(std::vector<std::string> fnames, std::vector<std::string> hlts, int var=0, int ien=0, int eta=0, int pv=0, bool varbin=false, int rebin=1, bool approve=false, std::string dtype="Data", bool logy=true, int pos=0, int coloff=0);
0569 TCanvas* plotEnergy(std::string fname="hlt.root", std::string HLT="All HLTs", int var=0, int ien=0, int eta=0, bool varbin=false, int rebin=1, bool approve=false, bool logy=true, int pos=0, int coloff=0);
0570 void plotEMeanPVAll(std::string fname="StudyHLT_ZeroBias_1PV.root", std::string HLT="Zero Bias", int var=-1, int eta=-1, bool approve=true);
0571 TCanvas* plotEMeanDrawPV(std::string fname="StudyHLT_ZeroBias_1PV.root", std::string HLT="Zero Bias", int var=0, int eta=0, bool approve=true);
0572 TCanvas* plotEnergyPV(std::string fnamee="StudyHLT_HLTZeroBias.root", std::string HLT="Zero Bias", int var=0, int ien=0, int eta=0, bool varbin=false, int rebin=1, bool approve=false, bool logy=true, int pos=0);
0573 TCanvas* plotTrack(std::string fname="hlt.root", std::string HLT="All HLTs", int var=0, bool varbin=false, int rebin=1, bool approve=false, bool logy=true, int pos=0);
0574 TCanvas* plotIsolation(std::string fname="hlt.root", std::string HLT="All HLTs", int var=0, bool varbin=false, int rebin=1, bool approve=false, bool logy=true, int pos=0);
0575 TCanvas* plotHLT(std::string fname="hlt.root", std::string HLT="All HLTs", int run=-1, bool varbin=false, int rebin=1, bool approve=false, bool logy=true, int pos=0);
0576 TCanvas* plotHisto(char* cname, std::string HLT, TObjArray& histArr, std::vector<std::string>& labels, std::vector<int>& color, char* name, double ymx0, bool logy, int pos, double yloff, double yhoff, double xmax=-1, bool varbin=false, int rebin=1, bool approve=false);
0577 void getHistStats(TH1D *h, int& entries, int& integral, double& mean, double& meanE, double& rms, double& rmsE, int& uflow, int& oflow);
0578 TFitResultPtr getHistFitStats(TH1F *h, const char* formula, double xlow, double xup, unsigned int& nPar, double* par, double* epar);
0579 void setHistAttr(TH1F *h, int icol, int lwid=1, int ltype=1);
0580 double getWeightedMean (int npt, int Start, std::vector<double>& mean, std::vector<double>& emean);
0581 TH1D* rebin(TH1D* histin, int);
0582
0583 void plotAll(std::string fname, std::string HLT, int var, int ien, int eta,
0584          bool varbin, int rebin, bool approve, bool logy, int pos, bool pv,
0585          int savePlot) {
0586   int varmin(0), varmax(5), enemin(0), enemax(9), etamin(0), etamax(3);
0587   if (var >= 0) varmin = varmax = var;
0588   if (ien >= 0) enemin = enemax = ien;
0589   if (eta >= 0) etamin = etamax = eta;
0590 //  std::cout << "Var " << varmin << ":" << varmax << " Ene " << enemin << ":"
0591 //      << enemax << " Eta " << etamin << ":" << etamax << std::endl;
0592   for (int var=varmin; var<=varmax; ++var) {
0593     for (int ene=enemin; ene<=enemax; ++ene) {
0594       for (int eta=etamin; eta<=etamax; ++eta) {
0595     TCanvas *c(0);
0596     if (pv) {
0597       c = plotEnergyPV(fname, HLT, var, ene, eta, varbin, rebin, approve, logy, pos);
0598     } else {
0599       c = plotEnergy(fname, HLT, var, ene, eta, varbin, rebin, approve, logy, pos);
0600     }
0601     if (c != 0 && savePlot >= 0 && savePlot < 3) {
0602       std::string ext[3] = {"eps", "gif", "pdf"};
0603       char name[200];
0604       sprintf (name, "%s.%s", c->GetName(), ext[savePlot].c_str());
0605       c->Print(name);
0606     }
0607       }
0608     }
0609   }
0610 }
0611
0612 void plotEnergyAll(std::string fname, std::string hlt, int models, int pv,
0613            int data, bool varbin, int rebin, bool approve, bool logy,
0614            int pos, int var, int ene, int eta, int savePlot) {
0615
0616   std::vector<std::string> fnames, hlts;
0617   std::string dtype = (data == 1) ? "Data" : "MC";
0618   int modeluse(models);
0619   int coloff = (data == 4) ? 0 : 1;
0620   if (fname == "") {
0621     if (data == 1) {
0622       for (int i=0; i<nmodels; ++i) {
0623     if (modeluse%2 == 1) {
0624       fnames.push_back(files[i]); hlts.push_back(types[i]);
0625     }
0626     modeluse /= 2;
0627       }
0628     } else if (data == 4) {
0629       for (int i=0; i<nmodelx; ++i) {
0630     if (modeluse%2 == 1) {
0631       fnames.push_back(filex[i]); hlts.push_back(typex[i]);
0632     }
0633     modeluse /= 2;
0634       }
0635     } else {
0636       for (int i=0; i<nmodelm; ++i) {
0637     if (modeluse%2 == 1) {
0638       fnames.push_back(filem[i]); hlts.push_back(typem[i]);
0639     }
0640     modeluse /= 2;
0641       }
0642     }
0643   } else {
0644     fnames.push_back(fname); hlts.push_back(hlt);
0645   }
0646   int varmin(0), varmax(5), enemin(0), enemax(9), etamin(0), etamax(3);
0647   if (var >= varmin && var <= varmax) varmin = varmax = var;
0648   if (ene >= enemin && ene <= enemax) enemin = enemax = ene;
0649   if (eta >= etamin && eta <= etamax) etamin = etamax = eta;
0650 //  std::cout << "Var " << varmin << ":" << varmax << " Ene " << enemin << ":"
0651 //      << enemax << " Eta " << etamin << ":" << etamax << std::endl;
0652   for (int var=varmin; var<=varmax; ++var) {
0653     for (int ene=enemin; ene<=enemax; ++ene) {
0654       for (int eta=etamin; eta<=etamax; ++eta) {
0655     TCanvas *c = plotEnergies(fnames, hlts, var, ene, eta, pv, varbin,
0656                   rebin, approve, dtype, logy, pos, coloff);
0657     if (c != 0 && savePlot >= 0 && savePlot < 3) {
0658       std::string ext[3] = {"eps", "gif", "pdf"};
0659       char name[200];
0660       sprintf (name, "%s.%s", c->GetName(), ext[savePlot].c_str());
0661       c->Print(name);
0662     }
0663       }
0664     }
0665   }
0666 }
0667
0668 void plotEMeanAll(int data, int models, bool ratio, bool approve,
0669           std::string postfix, int savePlot){
0670   int varmin(0), varmax(5), pvmin(0), pvmax(0), etamin(0), etamax(3);
0671 //  std::cout << "Var " << varmin << ":" << varmax << " PV " << pvmin << ":"
0672 //      << pvmax << " Eta " << etamin << ":" << etamax << std::endl;
0673   for (int var=varmin; var<=varmax; ++var) {
0674     for (int eta=etamin; eta<=etamax; ++eta) {
0675       for (int pv=pvmin; pv<=pvmax; ++pv) {
0676     plotEMean("", "", models, var, eta, pv, data, ratio, approve, postfix,
0677           savePlot);
0678       }
0679     }
0680   }
0681 }
0682
0683 void plotEMean(std::string fname, std::string hlt, int models, int var, int eta,
0684            int pv, int data, bool ratio, bool approve, std::string postfix,
0685            int savePlot) {
0686
0687   std::vector<std::string> fnames, hlts;
0688   std::string dtype = (data == 1) ? "Data" : "MC";
0689   int modeluse(models);
0690   int coloff = (data == 4 || data == 3) ? 0 : 1;
0691   if (fname == "") {
0692     if (data == 1) {
0693       for (int i=0; i<nmodels; ++i) {
0694     if (modeluse%2 == 1) {
0695       fnames.push_back(files[i]); hlts.push_back(types[i]);
0696     }
0697     modeluse /= 2;
0698       }
0699     } else if (data == 4) {
0700       for (int i=0; i<nmodelx; ++i) {
0701     if (modeluse%2 == 1) {
0702       fnames.push_back(filex[i]); hlts.push_back(typex[i]);
0703     }
0704     modeluse /= 2;
0705       }
0706     } else if (data == 3) {
0707       for (int i=0; i<nmodelx; ++i) {
0708     if (modeluse%2 == 1) {
0709       fnames.push_back(filex[i]); hlts.push_back(typex[i]);
0710     }
0711     modeluse /= 2;
0712       }
0713     } else {
0714       for (int i=0; i<nmodelm; ++i) {
0715     if (modeluse%2 == 1) {
0716       fnames.push_back(filem[i]); hlts.push_back(typem[i]);
0717     }
0718     modeluse /= 2;
0719       }
0720     }
0721   } else {
0722     fnames.push_back(fname); hlts.push_back(hlt);
0723   }
0724   int varmin(0), varmax(5), etamin(0), etamax(3), pvmin(0), pvmax(4);
0725   if (var >= 0) varmin = varmax = var;
0726   if (eta >= 0) etamin = etamax = eta;
0727   if (pv  >= 0) pvmin  = pvmax  = pv;
0728 //  std::cout << "Var " << varmin << ":" << varmax << " Eta " << etamin << ":"
0729 //      << etamax << std::endl;
0730   for (int var=varmin; var<=varmax; ++var) {
0731     for (int eta=etamin; eta<=etamax; ++eta) {
0732       for (int pv=pvmin; pv<=pvmax; ++pv) {
0733     TCanvas* c(0);
0734     if (ratio) {
0735       c = plotEMeanRatioDraw(fnames, hlts, var, eta, pv, approve, dtype, coloff);
0736     } else {
0737       c = plotEMeanDraw(fnames, hlts, var, eta, pv, approve, dtype, coloff);
0738     }
0739     if (c != 0 && savePlot >= 0 && savePlot < 3) {
0740       std::string ext[3] = {"eps", "gif", "pdf"};
0741       char name[200];
0742       sprintf (name, "%s%s.%s", c->GetName(), postfix.c_str(),
0743            ext[savePlot].c_str());
0744       c->Print(name);
0745     }
0746       }
0747     }
0748   }
0749 }
0750
0751 TCanvas* plotEMeanDraw(std::vector<std::string> fnames,
0752                std::vector<std::string> hlts, int var, int eta, int pv,
0753                bool approve, std::string dtype, int coloff) {
0754
0755   bool debug(false);
0756   std::vector<TGraphAsymmErrors*> graphs;
0757   TLegend*  legend = new TLegend(0.25, 0.80, 0.975, 0.95);
0758   legend->SetBorderSize(1); legend->SetFillColor(kWhite);
0759   legend->SetMargin(0.2);
0760   for (unsigned k=0; k<fnames.size(); ++k) {
0761     TFile *file = TFile::Open(fnames[k].c_str());
0762     double mean[NPT], dmean[NPT];
0763     for (int i=0; i<NPT; ++i) {
0764       char  name[100];
0765       sprintf (name, "h_energy_%d_%d_%d_%d", pv+3, i, eta, var);
0766       TH1D *histo = (TH1D*) file->FindObjectAny(name);
0767       if (histo) {
0768     mean[i] = histo->GetMean();
0769     dmean[i]= histo->GetMeanError();
0770       } else {
0771     mean[i] = -100.;
0772     dmean[i]= 0;
0773       }
0774     }
0775     if (debug) {
0776       std::cout << "Get mean for " << NPT << " points" << std::endl;
0777       for (int i=0; i<NPT; ++i)
0778     std::cout << "[" << i << "]" << " Momentum " << mom[i] << " +- "
0779           << dmom[i] << " Mean " << mean[i] << " +- " << dmean[i]
0780           << std::endl;
0781     }
0782     TGraphAsymmErrors *graph = new TGraphAsymmErrors(NPT, mom, mean, dmom,dmom, dmean,dmean);
0783     graph->SetMarkerStyle(styles[coloff+k]);
0784     graph->SetMarkerColor(colors[coloff+k]);
0785     graph->SetMarkerSize(1.6);
0786     graph->SetLineColor(colors[coloff+k]);
0787     graph->SetLineWidth(2);
0788     graphs.push_back(graph);
0789     legend->AddEntry(graph, hlts[k].c_str(), "lp");
0790     if (debug)  std::cout << "Complete " << hlts[k] << std::endl;
0791     file->Close();
0792   }
0793
0794   char cname[100], name[200];
0795   sprintf (cname, "c_%s_%d_%d_%s", varEne1[var].c_str(), eta, pv, dtype.c_str());
0796   gStyle->SetCanvasBorderMode(0); gStyle->SetCanvasColor(kWhite);
0797   gStyle->SetPadColor(kWhite);    gStyle->SetFillColor(kWhite);
0798   gStyle->SetOptTitle(kFALSE);    gStyle->SetPadBorderMode(0);
0799   gStyle->SetCanvasBorderMode(0);
0800   TCanvas *canvas = new TCanvas(cname, cname, 500, 400);
0801   gStyle->SetOptStat(0);     gPad->SetTopMargin(0.05);
0802   gPad->SetLeftMargin(0.15); gPad->SetRightMargin(0.025);
0803   gPad->SetBottomMargin(0.20);
0804   TH1F *vFrame = canvas->DrawFrame(0.0, 0.01, 50.0, 0.5);
0805   vFrame->GetYaxis()->SetRangeUser(0.0,1.6);
0806   vFrame->GetXaxis()->SetLabelSize(0.06);
0807   vFrame->GetYaxis()->SetLabelSize(0.05);
0808   vFrame->GetXaxis()->SetTitleSize(0.06);
0809   vFrame->GetYaxis()->SetTitleSize(0.06);
0810   vFrame->GetYaxis()->SetTitleOffset(0.9);
0811   vFrame->GetXaxis()->SetRangeUser(1.0,20.0);
0812   if (approve) {
0813     sprintf (name, "Mean of %s/p_{Track}", varEne[var].c_str());
0814   } else {
0815     sprintf (name, "<%s/p_{Track}>", varEne[var].c_str());
0816   }
0817   vFrame->GetYaxis()->SetTitle(name);
0818   sprintf (name, "p_{Track} (GeV/c)");
0819   vFrame->GetXaxis()->SetTitle(name);
0820   for (unsigned int ii=0; ii<graphs.size(); ++ii) graphs[ii]->Draw("P");
0821   legend->Draw();
0822   TLine *line = new TLine(1.0, 1.0, 20.0, 1.0);
0823   line->SetLineStyle(2); line->SetLineWidth(2); line->SetLineColor(kRed);
0824   line->Draw();
0825   TPaveText* text = new TPaveText(0.25, 0.74, 0.55, 0.79, "brNDC");
0826   if (approve) {
0827     sprintf (name, "(%s)", nameEtas[eta].c_str());
0828   } else {
0829     sprintf (name, "(%s, %s)", nameEta[eta].c_str(), namePV[pv].c_str());
0830   }
0831   if (debug) std::cout << "Name " << name << " |" << std::endl;
0832   text->AddText(name);
0833   text->Draw("same");
0834   TPaveText* text2 = new TPaveText(0.55, 0.74, 0.97, 0.79, "brNDC");
0835 //sprintf (name, "CMS Preliminary");
0836   sprintf (name, "CMS");
0837   text2->AddText(name);
0838   text2->Draw("same");
0839   return canvas;
0840 }
0841
0842 TCanvas* plotEMeanRatioDraw(std::vector<std::string> fnames,
0843                 std::vector<std::string> hlts, int var, int eta,
0844                 int pv, bool approve, std::string dtype,
0845                 int coloff) {
0846
0847   bool debug(false);
0848   std::vector<TGraphAsymmErrors*> graphs;
0849   TLegend*  legend = new TLegend(0.25, 0.80, 0.975, 0.95);
0850   legend->SetBorderSize(1); legend->SetFillColor(kWhite);
0851   legend->SetMargin(0.2);
0852   double mean0[NPT], dmean0[NPT];
0853   for (unsigned k=0; k<fnames.size(); ++k) {
0854     TFile *file = TFile::Open(fnames[k].c_str());
0855     double mean[NPT], dmean[NPT];
0856     for (int i=0; i<NPT; ++i) {
0857       char  name[100];
0858       sprintf (name, "h_energy_%d_%d_%d_%d", pv+3, i, eta, var);
0859       TH1D *histo = (TH1D*) file->FindObjectAny(name);
0860       if (histo) {
0861     mean[i] = histo->GetMean();
0862     dmean[i]= histo->GetMeanError();
0863       } else {
0864     mean[i] = -100.;
0865     dmean[i]= 0;
0866       }
0867     }
0868     if (debug) {
0869       std::cout << "Get mean for " << NPT << " points" << std::endl;
0870       for (int i=0; i<NPT; ++i)
0871     std::cout << "[" << i << "]" << " Momentum " << mom[i] << " +- "
0872           << dmom[i] << " Mean " << mean[i] << " +- " << dmean[i]
0873           << std::endl;
0874     }
0875     if (k == 0) {
0876       for (int i=0; i<NPT; ++i) {
0877     mean0[i] = mean[i];
0878     dmean0[i]= dmean[i];
0879       }
0880     } else {
0881       double sumNum(0), sumDen(0), sumNum1(0), sumDen1(0);
0882       for (int i=0; i<NPT; ++i) {
0883     if (dmean[i] > 0 && dmean0[i] > 0) {
0884       double er1 = dmean[i]/mean[i];
0885       double er2 = dmean0[i]/mean0[i];
0886       mean[i] = mean[i]/mean0[i];
0887       dmean[i]= mean[i]*sqrt(er1*er1+er2*er2);
0888       double temp1 = (mean[i]>1.0) ? 1.0/mean[i] : mean[i];
0889       double temp2 = (mean[i]>1.0) ? dmean[i]/(mean[i]*mean[i]) : dmean[i];
0890       if (i > 0) {
0891         sumNum += (fabs(1-temp1)/(temp2*temp2));
0892         sumDen += (1.0/(temp2*temp2));
0893       }
0894       sumNum1 += (fabs(1-temp1)/(temp2*temp2));
0895       sumDen1 += (1.0/(temp2*temp2));
0896     } else {
0897       mean[i] = -100.;
0898       dmean[i]= 0;
0899     }
0900       }
0901       sumNum  = (sumDen>0)  ? (sumNum/sumDen) : 0;
0902       sumDen  = (sumDen>0)  ? 1.0/sqrt(sumDen) : 0;
0903       sumNum1 = (sumDen1>0) ? (sumNum1/sumDen1) : 0;
0904       sumDen1 = (sumDen1>0) ? 1.0/sqrt(sumDen1) : 0;
0905       std::cout << "Get Ratio of mean for " << NPT << " points: Mean " << sumNum << " +- " << sumDen << " (" << sumNum1 << " +- " << sumDen1 << ") Input: " << fnames[k] << " var|eta|pv " << var << ":" << eta << ":" << pv << std::endl;
0906       if (debug) {
0907     std::cout << "Get Ratio of mean for " << NPT << " points: Mean "
0908           << sumNum << " +- " << sumDen << std::endl;
0909     for (int i=0; i<NPT; ++i)
0910       std::cout << "[" << i << "]" << " Momentum " << mom[i] << " +- "
0911             << dmom[i] << " Mean " << mean[i] << " +- " << dmean[i]
0912             << std::endl;
0913       }
0914       TGraphAsymmErrors *graph = new TGraphAsymmErrors(NPT, mom, mean, dmom,dmom, dmean,dmean);
0915       graph->SetMarkerStyle(styles[coloff+k]);
0916       graph->SetMarkerColor(colors[coloff+k]);
0917       graph->SetMarkerSize(1.6);
0918       graph->SetLineColor(colors[coloff+k]);
0919       graph->SetLineWidth(2);
0920       graphs.push_back(graph);
0921       char text[100];
0922       if (approve) {
0923     sprintf (text,"%s", hlts[k].c_str());
0924       } else {
0925     sprintf (text,"%5.3f #pm %5.3f  %s", sumNum, sumDen, hlts[k].c_str());
0926       }
0927       legend->AddEntry(graph, text, "lp");
0928       if (debug) std::cout << "Complete " << hlts[k] << std::endl;
0929     }
0930     file->Close();
0931   }
0932
0933   char cname[100], name[200];
0934   sprintf (cname, "cR_%s_%d_%d_%s", varEne1[var].c_str(), eta, pv, dtype.c_str());
0935   gStyle->SetCanvasBorderMode(0); gStyle->SetCanvasColor(kWhite);
0936   gStyle->SetPadColor(kWhite);    gStyle->SetFillColor(kWhite);
0937   gStyle->SetOptTitle(kFALSE);    gStyle->SetPadBorderMode(0);
0938   gStyle->SetCanvasBorderMode(0);
0939   TCanvas *canvas = new TCanvas(cname, cname, 500, 400);
0940   gStyle->SetOptStat(0);     gPad->SetTopMargin(0.05);
0941   gPad->SetLeftMargin(0.15); gPad->SetRightMargin(0.025);
0942   gPad->SetBottomMargin(0.20);
0943   TH1F *vFrame = canvas->DrawFrame(0.0, 0.01, 50.0, 0.5);
0944   vFrame->GetYaxis()->SetRangeUser(0.4,1.5);
0945   vFrame->GetXaxis()->SetLabelSize(0.06);
0946   vFrame->GetYaxis()->SetLabelSize(0.05);
0947   vFrame->GetXaxis()->SetTitleSize(0.06);
0948   vFrame->GetYaxis()->SetTitleSize(0.045);
0949   vFrame->GetYaxis()->SetTitleOffset(1.2);
0950   vFrame->GetXaxis()->SetRangeUser(1.0,20.0);
0951   if (approve) {
0952     sprintf (name, "%s/Data for mean of %s/p_{Track}", dtype.c_str(), varEne[var].c_str());
0953   } else {
0954     sprintf (name, "#frac{%s}{%s} for #frac{%s}{p_{Track}}", dtype.c_str(), hlts[0].c_str(), varEne[var].c_str());
0955   }
0956   vFrame->GetYaxis()->SetTitle(name);
0957   sprintf (name, "p_{Track} (GeV/c)");
0958   vFrame->GetXaxis()->SetTitle(name);
0959   for (unsigned int ii=0; ii<graphs.size(); ++ii) graphs[ii]->Draw("P");
0960   legend->Draw();
0961   TLine *line = new TLine(1.0, 1.0, 20.0, 1.0);
0962   line->SetLineStyle(2); line->SetLineWidth(2); line->SetLineColor(kRed);
0963   line->Draw();
0964   TPaveText* text = new TPaveText(0.55, 0.40, 0.95, 0.45, "brNDC");
0965   if (approve) {
0966     sprintf (name, "(%s)", nameEtas[eta].c_str());
0967   } else {
0968     sprintf (name, "(%s, %s)", nameEta[eta].c_str(), namePV[pv].c_str());
0969   }
0970   if (debug) std::cout << "Name " << name << " |" << std::endl;
0971   text->AddText(name);
0972   text->Draw("same");
0973   TPaveText* text2 = new TPaveText(0.55, 0.45, 0.95, 0.50, "brNDC");
0974 //sprintf (name, "CMS Preliminary");
0975   sprintf (name, "CMS");
0976   text2->AddText(name);
0977   text2->Draw("same");
0978   return canvas;
0979 }
0980
0981 TCanvas* plotEnergies(std::vector<std::string> fnames,
0982               std::vector<std::string> hlts, int var, int ien, int eta,
0983               int pv, bool varbin, int rebin, bool approve,
0984               std::string dtype, bool logy, int pos, int coloff) {
0985
0986   //  bool debug(false);
0987   TLegend*  legend = new TLegend(0.55, 0.70, 0.95, 0.85);
0988   legend->SetBorderSize(1); legend->SetFillColor(kWhite);
0989   legend->SetMargin(0.4);
0990   TObjArray                histArr;
0991   char                     name[100];
0992   std::vector<std::string> labels;
0993   std::vector<int>         color;
0994   double                   ymx0(0), ent0(0);
0995   sprintf (name, "h_energy_%d_%d_%d_%d", pv+3, ien, eta, var);
0996   for (unsigned k=0; k<fnames.size(); ++k) {
0997     TFile *file = TFile::Open(fnames[k].c_str());
0998     TH1D  *histo = (TH1D*) file->FindObjectAny(name);
0999     if (histo) {
1000       if (k == 0) {
1001     ent0 = histo->GetEntries();
1002       } else {
1003     double scale = (histo->GetEntries() > 0) ? ent0/histo->GetEntries() : 1;
1004     histo->Scale(scale);
1005       }
1006       histArr.AddLast(histo);
1007       labels.push_back(hlts[k]);
1008       color.push_back(colors[coloff+k]);
1009       int ibin = histo->GetMaximumBin();
1010       if (histo->GetBinContent(ibin) > ymx0) ymx0 = histo->GetBinContent(ibin);
1011     }
1012   }
1013   TCanvas* c(0);
1014   if (histArr.GetEntries()>0) {
1015     sprintf (name, "p=%s, %s", varPPs[ien].c_str(), nameEtas[eta].c_str());
1016     /*
1017     if (approve) {
1018       sprintf (name, "p=%s, %s", varPPs[ien].c_str(), nameEtas[eta].c_str());
1019     } else {
1020       sprintf (name, "p=%s, i#eta=%s, %s", varPs[ien].c_str(), varEta[eta].c_str(), namePV[pv].c_str());
1021     }
1022     */
1023     std::string clabel(name);
1024     char cname[50];
1025     sprintf (cname, "c_%s_%d_%d_%s", varEne1[var].c_str(), ien, eta, dtype.c_str());
1026     sprintf ( name, "%s/p", varEne[var].c_str());
1027     c = plotHisto(cname, clabel, histArr, labels, color, name, ymx0, logy, pos,
1028           0.10, 0.05, 2.5, varbin, rebin, approve);
1029   }
1030   return c;
1031 }
1032
1033 TCanvas* plotEnergy(std::string fname, std::string HLT, int var, int ien,
1034             int eta, bool varbin, int rebin, bool approve, bool logy,
1035             int pos, int coloff) {
1036
1037   TFile *file = TFile::Open(fname.c_str());
1038   char                     name[100];
1039   TObjArray                histArr;
1040   std::vector<std::string> labels;
1041   std::vector<int>         color;
1042   double                   ymx0(0);
1043   for (int i=0; i<4; ++i) {
1044     sprintf (name, "h_energy_%d_%d_%d_%d", i, ien, eta, var);
1045     TH1D *histo = (TH1D*) file->FindObjectAny(name);
1046     if (histo) {
1047       histArr.AddLast(histo);
1048       sprintf (name, "p=%s, #eta=%s %s", varPs[ien].c_str(), varEta[eta].c_str(), namefull[i+3].c_str());
1049       labels.push_back(name);
1050       color.push_back(colors[coloff+i]);
1051       int ibin = histo->GetMaximumBin();
1052       if (histo->GetBinContent(ibin) > ymx0) ymx0 = histo->GetBinContent(ibin);
1053     }
1054   }
1055   TCanvas* c(0);
1056   if (histArr.GetEntries()>0) {
1057     char cname[50];
1058     sprintf (cname, "c_%s_%d_%d", varEne1[var].c_str(), ien, eta);
1059     sprintf ( name, "%s/p", varEne[var].c_str());
1060     c = plotHisto(cname, HLT, histArr, labels, color, name, ymx0, logy, pos, 0.10, 0.05, 2.5, varbin, rebin, approve);
1061   }
1062   return c;
1063 }
1064
1065 void plotEMeanPVAll(std::string fname, std::string HLT, int var, int eta,
1066             bool approve) {
1067
1068   int varmin(0), varmax(5), etamin(0), etamax(3);
1069   if (var >= 0) varmin = varmax = var;
1070   if (eta >= 0) etamin = etamax = eta;
1071   for (int var=varmin; var<=varmax; ++var) {
1072     for (int eta=etamin; eta<=etamax; ++eta) {
1073       plotEMeanDrawPV(fname, HLT, var, eta, approve);
1074     }
1075   }
1076 }
1077
1078 TCanvas* plotEMeanDrawPV(std::string fname, std::string HLT, int var, int eta,
1079              bool approve) {
1080
1081   bool debug(false);
1082   std::vector<TGraphAsymmErrors*> graphs;
1083   TLegend*  legend = new TLegend(0.575, 0.80, 0.975, 0.95);
1084   legend->SetBorderSize(1); legend->SetFillColor(kWhite);
1085   legend->SetMargin(0.4);
1086   TFile *file = TFile::Open(fname.c_str());
1087   const int nPVBin=4;
1088   int    pvBins[nPVBin+1] = {1, 2, 3, 5, 100};
1089   for (int k=0; k<nPVBin; ++k) {
1090     char  name[100];
1091     double mean[NPT], dmean[NPT];
1092     for (int i=0; i<NPT; ++i) {
1093       sprintf (name, "h_energy_%d_%d_%d_%d", k, i, eta, var);
1094       TH1D *histo = (TH1D*) file->FindObjectAny(name);
1095       if (histo) {
1096     mean[i] = histo->GetMean();
1097     dmean[i]= histo->GetMeanError();
1098       } else {
1099     mean[i] = -100.;
1100     dmean[i]= 0;
1101       }
1102     }
1103     if (debug) {
1104       std::cout << "Get mean for " << NPT << " points" << std::endl;
1105       for (int i=0; i<NPT; ++i)
1106     std::cout << "[" << i << "]" << " Momentum " << mom[i] << " +- "
1107           << dmom[i] << " Mean " << mean[i] << " +- " << dmean[i]
1108           << std::endl;
1109     }
1110     TGraphAsymmErrors *graph = new TGraphAsymmErrors(NPT, mom, mean, dmom,dmom, dmean,dmean);
1111     graph->SetMarkerStyle(styles[k]);
1112     graph->SetMarkerColor(colors[k]);
1113     graph->SetMarkerSize(1.6);
1114     graph->SetLineColor(colors[k]);
1115     graph->SetLineWidth(2);
1116     graphs.push_back(graph);
1117     sprintf (name, "PV=%d:%d", pvBins[k], pvBins[k+1]-1);
1118     legend->AddEntry(graph, name, "lp");
1119     if (debug)  std::cout << "Complete " << name << std::endl;
1120   }
1121   file->Close();
1122
1123   char cname[100], name[200];
1124   sprintf (cname, "c_%s_PV_%d", varEne1[var].c_str(), eta);
1125   gStyle->SetCanvasBorderMode(0); gStyle->SetCanvasColor(kWhite);
1126   gStyle->SetPadColor(kWhite);    gStyle->SetFillColor(kWhite);
1127   gStyle->SetOptTitle(kFALSE);    gStyle->SetPadBorderMode(0);
1128   gStyle->SetCanvasBorderMode(0);
1129   TCanvas *canvas = new TCanvas(cname, cname, 500, 400);
1130   gStyle->SetOptStat(0);     gPad->SetTopMargin(0.05);
1131   gPad->SetLeftMargin(0.15); gPad->SetRightMargin(0.025);
1132   gPad->SetBottomMargin(0.20);
1133   TH1F *vFrame = canvas->DrawFrame(0.0, 0.01, 50.0, 0.5);
1134   vFrame->GetYaxis()->SetRangeUser(0.0,1.5);
1135   vFrame->GetXaxis()->SetLabelSize(0.06);
1136   vFrame->GetYaxis()->SetLabelSize(0.05);
1137   vFrame->GetXaxis()->SetTitleSize(0.06);
1138   vFrame->GetYaxis()->SetTitleSize(0.06);
1139   vFrame->GetYaxis()->SetTitleOffset(0.9);
1140   vFrame->GetXaxis()->SetRangeUser(1.0,20.0);
1141   if (approve) {
1142     sprintf (name, "Mean of %s/p_{Track}", varEne[var].c_str());
1143   } else {
1144     sprintf (name, "<%s/p_{Track}>", varEne[var].c_str());
1145   }
1146   vFrame->GetYaxis()->SetTitle(name);
1147   sprintf (name, "p_{Track} (GeV/c)");
1148   vFrame->GetXaxis()->SetTitle(name);
1149   for (unsigned int ii=0; ii<graphs.size(); ++ii) graphs[ii]->Draw("P");
1150   legend->Draw();
1151   TLine *line = new TLine(1.0, 1.0, 20.0, 1.0);
1152   line->SetLineStyle(2); line->SetLineWidth(2); line->SetLineColor(kRed);
1153   line->Draw();
1154   TPaveText* text = new TPaveText(0.575, 0.75, 0.97, 0.79, "brNDC");
1155   sprintf (name, "%s (%s)", HLT.c_str(), nameEtas[eta].c_str());
1156   if (debug) std::cout << "Name " << name << " |" << std::endl;
1157   text->AddText(name);
1158   text->Draw("same");
1159   TPaveText* text2 = new TPaveText(0.575, 0.71, 0.97, 0.75, "brNDC");
1160 //sprintf (name, "CMS Preliminary");
1161   sprintf (name, "CMS");
1162   text2->AddText(name);
1163   text2->Draw("same");
1164   return canvas;
1165 }
1166
1167 TCanvas* plotEnergyPV(std::string fname, std::string HLT, int var, int ien,
1168               int eta, bool varbin, int rebin, bool approve, bool logy,
1169               int pos) {
1170
1171   const int nPVBin=4;
1172   int    pvBins[nPVBin+1] = {1, 2, 3, 5, 100};
1173   TFile *file = TFile::Open(fname.c_str());
1174   char                     name[100];
1175   TObjArray                histArr;
1176   std::vector<std::string> labels;
1177   std::vector<int>         color;
1178   double                   ymx0(0);
1179   for (int i=4; i<nPVBin+4; ++i) {
1180     sprintf (name, "h_energy_%d_%d_%d_%d", i, ien, eta, var);
1181     TH1D *histo = (TH1D*) file->FindObjectAny(name);
1182     if (histo) {
1183       histArr.AddLast(histo);
1184       sprintf (name, "p=%s, #eta=%s, PV=%d:%d (%s)", varPs[ien].c_str(), varEta[eta].c_str(), pvBins[i-4], pvBins[i-3]-1, namefull[6].c_str());
1185       labels.push_back(name);
1186       color.push_back(colors[i-4]);
1187       int ibin = histo->GetMaximumBin();
1188       if (histo->GetBinContent(ibin) > ymx0) ymx0 = histo->GetBinContent(ibin);
1189     }
1190   }
1191   TCanvas *c(0);
1192   if (histArr.GetEntries()>0) {
1193     char cname[50];
1194     sprintf (cname, "c_%s_%d_%d", varEne1[var].c_str(), ien, eta);
1195     sprintf ( name, "%s/p", varEne[var].c_str());
1196     c = plotHisto(cname, HLT, histArr, labels, color, name, ymx0, logy, pos,
1197           0.10, 0.05, 2.5, varbin, rebin, approve);
1198   }
1199   return c;
1200 }
1201
1202 TCanvas* plotTrack(std::string fname, std::string HLT, int var, bool varbin,
1203            int rebin, bool approve, bool logy, int pos) {
1204
1205   TFile *file = TFile::Open(fname.c_str());
1206   char                     name[100];
1207   TObjArray                histArr;
1208   std::vector<std::string> labels;
1209   std::vector<int>         color;
1210   double                   ymx0(0);
1211   for (int i=0; i<7; ++i) {
1212     sprintf (name, "h_%s_%s", varname[var].c_str(), names[i].c_str());
1213     TH1D *histo = (TH1D*) file->FindObjectAny(name);
1214     if (histo) {
1215       histArr.AddLast(histo);
1216       labels.push_back(namefull[i]);
1217       color.push_back(colors[i]);
1218       int ibin = histo->GetMaximumBin();
1219       if (histo->GetBinContent(ibin) > ymx0) ymx0 = histo->GetBinContent(ibin);
1220     }
1221   }
1222   if (histArr.GetEntries()>0) {
1223     char cname[50];
1224     sprintf (cname, "c_%s", varname[var].c_str());
1225     sprintf ( name, "%s", vartitle[var].c_str());
1226     return plotHisto(cname, HLT, histArr, labels, color, name, ymx0, logy, pos,
1227              0.10, 0.05, -1, varbin, rebin, approve);
1228   } else {
1229     return 0;
1230   }
1231 }
1232
1233 TCanvas* plotIsolation(std::string fname, std::string HLT, int var, bool varbin,
1234                int rebin, bool approve, bool logy, int pos) {
1235
1236   TFile *file = TFile::Open(fname.c_str());
1237   char                     name[100];
1238   TObjArray                histArr;
1239   std::vector<std::string> labels;
1240   std::vector<int>         color;
1241   double                   ymx0(0);
1242   for (int i=0; i<2; ++i) {
1243     sprintf (name, "h_%s_%s", varnameC[var].c_str(), nameC[i].c_str());
1244     TH1D *histo = (TH1D*) file->FindObjectAny(name);
1245     if (histo) {
1246       histArr.AddLast(histo);
1247       labels.push_back(nameCF[i]);
1248       color.push_back(colors[i]);
1249       int ibin = histo->GetMaximumBin();
1250       if (histo->GetBinContent(ibin) > ymx0) ymx0 = histo->GetBinContent(ibin);
1251     }
1252   }
1253   if (histArr.GetEntries()>0) {
1254     char cname[50];
1255     sprintf (cname, "c_%s", varnameC[var].c_str());
1256     sprintf ( name, "%s (GeV)", vartitlC[var].c_str());
1257     return plotHisto(cname, HLT, histArr, labels, color, name, ymx0, logy, pos,
1258              0.10, 0.05, -1, varbin, rebin, approve);
1259   } else {
1260     return 0;
1261   }
1262 }
1263
1264 TCanvas* plotHLT(std::string fname, std::string HLT, int run, bool varbin,
1265          int rebin, bool approve, bool logy, int pos) {
1266
1267   TFile *file = TFile::Open(fname.c_str());
1268   char                     name[100];
1269   TObjArray                histArr;
1270   std::vector<std::string> labels;
1271   std::vector<int>         color;
1272   double                   ymx0(0);
1273   if (run > 0) sprintf (name, "h_HLTAccepts_%d", run);
1274   else         sprintf (name, "h_HLTAccept");
1275   TH1D *histo = (TH1D*) file->FindObjectAny(name);
1276   if (histo) {
1277     histArr.AddLast(histo);
1278     labels.push_back(HLT);
1279     color.push_back(colors[3]);
1280     int ibin = histo->GetMaximumBin();
1281     ymx0 = histo->GetBinContent(ibin);
1282   }
1283   if (histArr.GetEntries()>0) {
1284     char cname[50], hname[50];
1285     if (run > 0) {sprintf (cname,"c_HLT_%d",run); sprintf (name,"Run %d",run);
1286     }  else      {sprintf (cname,"c_HLTs");       sprintf (name, "All runs");}
1287     sprintf (hname, " ");
1288     return plotHisto(cname, "", histArr, labels, color, hname, ymx0, logy, pos,
1289              0.40, 0.01, -1, varbin, rebin, approve);
1290   } else {
1291     return 0;
1292   }
1293 }
1294
1295 TCanvas* plotHisto(char* cname, std::string HLT, TObjArray& histArr,
1296            std::vector<std::string>& labels, std::vector<int>& color,
1297            char* name, double ymx0, bool logy, int pos, double yloff,
1298            double yhoff, double xmax, bool varbin, int rebin0,
1299            bool approve) {
1300
1301   int nentry = histArr.GetEntries();
1302   double ymax = 10.;
1303   for (int i=0; i<10; ++i) {
1304     if (ymx0 < ymax) break;
1305     ymax *= 10.;
1306   }
1307   double ystep = ymax*0.1;
1308   for (int i=0; i<9; ++i) {
1309     if (ymax-ystep < ymx0) break;
1310     ymax -= ystep;
1311   }
1312   double ymin(0);
1313   if (logy) ymin = 1;
1314 //  std::cout << "ymin " << ymin << " ymax " << ymx0 << ":" << ymax  << std::endl;
1315
1316   gStyle->SetCanvasBorderMode(0); gStyle->SetCanvasColor(kWhite);
1317   gStyle->SetPadColor(kWhite);    gStyle->SetFillColor(kWhite);
1318   gStyle->SetOptTitle(kFALSE);    gStyle->SetPadBorderMode(0);
1319   gStyle->SetCanvasBorderMode(0);
1320   if (approve) gStyle->SetOptStat(0);
1321   else         gStyle->SetOptStat(1110);
1322   TCanvas *canvas = new TCanvas(cname, cname, 500, 500);
1323   gPad->SetTopMargin(yhoff);
1324   gPad->SetLeftMargin(0.15); gPad->SetRightMargin(0.025);
1325   gPad->SetBottomMargin(yloff);
1326   if (logy) canvas->SetLogy();
1327   double height = 0.08;
1328   double dx    = (nentry > 2) ? 0.50 : 0.30;
1329   double dy    = (nentry > 2) ? 0.12 : 0.04;
1330   double dy2   = 0.035;
1331   double xmin1 = (pos > 1) ? 0.375 : 0.75-dx;
1332   double xmin2 = (pos > 1) ? 0.12 : 0.65;
1333   double xmin3 = (pos > 1) ? 0.15 : 0.75;
1334   double ymin1 = (pos%2 == 0) ? (1.0-yhoff-dy) : (yloff+0.02);
1335   double ymin2 = (pos%2 == 0) ? (0.96-yhoff-nentry*height) : (yloff+0.025+nentry*height);
1336   double ymin3 = (pos%2 == 0) ? (ymin2-dy2) : (ymin2+dy2);
1337   double dx2   = (approve) ? dx : 0.32;
1338   double dx3   = (approve) ? dx : 0.22;
1339   if (approve) {
1340     xmin1 = xmin2 = xmin3 = 0.975-dx;
1341     ymin1 = (1.0-yhoff-dy);
1342     ymin2 = ymin1-dy2-0.01;
1343     ymin3 = ymin2-dy2;
1344   }
1345 //  std::cout << nentry << " " << pos << " " << dx << " " << dy << " " << xmin1 << " " << xmin2 << std::endl;
1346   TLegend  *legend = new TLegend(xmin1, ymin1, xmin1+dx, ymin1+dy);
1347   TPaveText *text, *text2;
1348   if (varbin || rebin0 == 1) {
1349     text  = new TPaveText(xmin2, ymin2, xmin2+dx2, ymin2+dy2, "brNDC");
1350     text2 = new TPaveText(xmin3, ymin3, xmin3+dx3, ymin3+dy2, "brNDC");
1351   } else {
1352     text  = new TPaveText(0.10,     0.95, dx2+0.10,     dy2+0.95, "brNDC");
1353     text2 = new TPaveText(dx2+0.10, 0.95, dx2+dx3+0.10, dy2+0.95, "brNDC");
1354   }
1355   legend->SetBorderSize(1); legend->SetFillColor(kWhite);
1356   char texts[200];
1357 //sprintf (texts, "CMS Preliminary");
1358   sprintf (texts, "CMS");
1359   text2->AddText(texts);
1360   THStack *Hs      = new THStack("hs2"," ");
1361   for (int i=0; i<nentry; i++) {
1362     TH1D *h =  (varbin) ? rebin((TH1D*)histArr[i], i) : (TH1D*)((TH1D*)histArr[i])->Rebin(rebin0);
1363     h->SetLineColor(color[i]);
1364     h->SetLineWidth(2);
1365     h->SetMarkerSize(0.8);
1366     double ymax0 = (varbin) ? ymax : rebin0*ymax;
1367     h->GetYaxis()->SetRangeUser(ymin,ymax0);
1368     if (xmax > 0 && (!varbin)) h->GetXaxis()->SetRangeUser(0,xmax);
1369     Hs->Add(h, "hist sames");
1370     legend->AddEntry(h,labels[i].c_str(),"l");
1371   }
1372   Hs->Draw("nostack");
1373   canvas->Update();
1374   Hs->GetHistogram()->GetXaxis()->SetTitle(name);
1375   Hs->GetHistogram()->GetXaxis()->SetLabelSize(0.035);
1376   Hs->GetHistogram()->GetYaxis()->SetTitleOffset(1.6);
1377   if (varbin) {
1378     Hs->GetHistogram()->GetYaxis()->SetTitle("Tracks/0.01");
1379   } else {
1380     Hs->GetHistogram()->GetYaxis()->SetTitle("Tracks");
1381     if (xmax > 0) Hs->GetHistogram()->GetXaxis()->SetRangeUser(0,xmax);
1382   }
1383   canvas->Modified();
1384
1385   canvas->Update();
1386   if (!approve) {
1387     for (int i=0; i<nentry; i++) {
1388       TH1D *h =  (TH1D*)histArr[i];
1389       if (h != NULL) {
1390     TPaveStats* st1 = (TPaveStats*)h->GetListOfFunctions()->FindObject("stats");
1391     if (st1 != NULL) {
1392       if (pos%2 == 0) {
1393         st1->SetY1NDC(1.0-yhoff-(i+1)*height); st1->SetY2NDC(1.0-yhoff-i*height);
1394       } else {
1395         st1->SetY1NDC(yloff+0.02+i*height); st1->SetY2NDC(yloff+0.02+(i+1)*height);
1396       }
1397       if (pos > 1) {
1398         st1->SetX1NDC(0.15); st1->SetX2NDC(.375);
1399       } else {
1400         st1->SetX1NDC(0.75); st1->SetX2NDC(.975);
1401       }
1402       st1->SetTextColor(colors[i]);
1403     }
1404       }
1405     }
1406   }
1407   legend->Draw("");
1408   if (HLT != "") {
1409     text->AddText(HLT.c_str());
1410     text->Draw("same");
1411   }
1412   text2->Draw("same");
1413
1414   return canvas;
1415 }
1416
1417 void getHistStats(TH1D *h, int& entries, int& integral, double& mean, double& meanE, double& rms, double& rmsE, int& uflow, int& oflow) {
1418   entries  = h->GetEntries();
1419   integral = h->Integral();
1420   mean     = h->GetMean();
1421   meanE    = h->GetMeanError();
1422   rms      = h->GetRMS();
1423   rmsE     = h->GetRMSError();
1424   uflow    = h->GetBinContent(0);
1425   oflow    = h->GetBinContent( h->GetNbinsX()+1 );
1426 }
1427
1428 TFitResultPtr getHistFitStats(TH1F *h, const char* formula, double xlow,
1429                   double xup, unsigned int& nPar, double* par,
1430                   double* epar) {
1431
1432   TFitResultPtr fit =  h->Fit(formula, "+qRB0", "", xlow, xup);
1433   nPar = fit->NPar();
1434   const std::vector<double> errors = fit->Errors();
1435   for (unsigned int i=0; i<nPar; i++) {
1436     par[i]  = fit->Value(i);
1437     epar[i] = errors[i];
1438   }
1439   return fit;
1440 }
1441
1442 void setHistAttr(TH1F *h, int icol, int lwid, int ltype) {
1443   h->SetLineColor(icol);
1444   h->SetLineStyle(ltype);
1445   h->SetLineWidth(lwid);
1446   TF1 *f = h->GetFunction("gaus");
1447   if (!f->IsZombie()) {
1448     f->SetLineColor(icol);
1449     f->SetLineStyle(2);
1450   }
1451 }
1452
1453 double getWeightedMean (int npt, int Start, std::vector<double>& mean, std::vector<double>& emean) {
1454   double sumDen=0, sumNum=0;
1455   for (int i=Start; i<npt; i++){
1456     if (mean[i]==0.0 || emean[i]==0.0) {
1457       sumNum += 0;
1458       sumDen += 0;
1459     } else {
1460       sumNum += mean[i]/(emean[i]*emean[i]);
1461       sumDen += 1.0/(emean[i]*emean[i]);
1462     }
1463   }
1464   double WeightedMean = sumNum/sumDen;
1465   return WeightedMean;
1466 }
1467
1468 TH1D* rebin(TH1D* histin, int indx) {
1469
1470   std::string nameIn(histin->GetName());
1471   char name[200];
1472   sprintf (name, "%sRebin%d", nameIn.c_str(), indx);
1473   TH1D* hist = new TH1D(name,histin->GetXaxis()->GetTitle(),nbins,xbins);
1474   std::vector<double> cont;
1475   for (int i=1; i<=histin->GetNbinsX(); ++i) {
1476     double value = histin->GetBinContent(i);
1477     cont.push_back(value);
1478 //    std::cout << "Input Bin [" << i << "] = " << value << std::endl;
1479   }
1480   for (int i=0; i<nbins; ++i) {
1481     double totl = 0;
1482     int    kount= 0;
1483     int    klow = ibins[i];
1484     int    khigh= ibins[i+1];
1485 //    std::cout << "i " << i << ":" << klow << ":" << khigh << std::endl;
1486     for (int k=klow; k<khigh; ++k) {
1487       totl += cont[k-1];
1488       kount++;
1489 //      std::cout << "K " << k << ":" << kount << ":" << cont[k-1] << ":" << totl << std::endl;
1490     }
1491     if (kount>0) totl /= kount;
1492 //    std::cout << "Output Bin [" << i+1 << "] Count " << kount << " Value " << totl << std::endl;
1493     hist->SetBinContent(i+1,totl);
1494   }
1495   return hist;
1496 }