Back to home page

Project CMSSW displayed by LXR

 
 

    


File indexing completed on 2021-02-14 14:14:04

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