Back to home page

Project CMSSW displayed by LXR
 
 

    

File indexing completed on 2021-12-15 05:16:54 

0001 //////////////////////////////////////////////////////////////////////////////
0002 // Usage:
0003 // .L CalibFitPlots.C+g
0004 //             For standard set of histograms from CalibMonitor
0005 //  FitHistStandard(infile, outfile, prefix, mode, type, append, saveAll,
0006 //                  debug);
0007 //      Defaults: mode=11111, type=0, append=true, saveAll=false, debug=false
0008 //
0009 //             For extended set of histograms from CalibMonitor
0010 //  FitHistExtended(infile, outfile, prefix, numb, type, append, fiteta, iname,
0011 //                  debug);
0012 //      Defaults: numb=50, type=3, append=true, fiteta=true, iname=3,
0013 //                debug=false
0014 //
0015 //             For RBX dependence in sets of histograms from CalibMonitor
0016 //  FitHistRBX(infile, outfile, prefix, append, iname);
0017 //      Defaults: append=true, iname=2
0018 //
0019 //             For plotting stored histograms from FitHist's
0020 //  PlotHist(infile, prefix, text, modePlot, kopt, lumi, ener, dataMC,
0021 //           drawStatBox, save);
0022 //      Defaults: modePlot=4, kopt=100, lumi=0, ener=13, dataMC=false,
0023 //                drawStatBox=true, save=0
0024 //
0025 //             For plotting histograms corresponding to individual ieta's
0026 //  PlotHistEta(infile, prefix, text, iene, numb, ieta, lumi, ener, dataMC,
0027 //           drawStatBox, save);
0028 //      Defaults iene=3, numb=50, ieta=0, lumi=0, ener=13.0, dataMC=false,
0029 //                drawStatBox=true, save=0
0030 //
0031 //             For plotting several histograms in the same plot
0032 //             (fits to different data sets for example)
0033 //  PlotHists(infile, prefix, text, drawStatBox, save)
0034 //      Defaults: drawStatBox=true; save=0;
0035 //      Note prefix is common part for all histograms
0036 //
0037 //             For plotting on the same canvas plots with different
0038 //             prefixes residing in the same file with approrprate text
0039 //  PlotTwoHists(infile, prefix1, text1, prefix2, text2, text0, type, iname,
0040 //               lumi, ener, drawStatBox, save);
0041 //      Defaults: type=0; iname=2; lumi=0; ener=13; drawStatBox=true;
0042 //                save=0;
0043 //      Note prefixN, textN have the same meaning as prefix and text for set N
0044 //           text0 is the text for general title added within ()
0045 //           type=0 plots response distributions and MPV of response vs ieta
0046 //               =1 plots MPV of response vs RBX #
0047 //
0048 //             For plotting stored histograms from CalibTree
0049 //  PlotFiveHists(infile, text0, prefix0, type, iname, drawStatBox, normalize,
0050 //                save, prefix1, text1, prefix2, text2, prefix3, text3,
0051 //                 prefix4, text4, prefix5, text5);
0052 //      Defaults: type=0; iname=0; drawStatBox=true; normalize=false;
0053 //                save=0; prefixN=""; textN=""; (for N > 0)
0054 //      Note prefixN, textN have the same meaning as prefix and text for set N
0055 //           text0 is the text for general title added within ()
0056 //           prefix0 is the tag attached to the canvas name
0057 //           type has the same meaning as in PlotTwoHists
0058 //
0059 //  PlotHistCorrResults(infile, text, prefixF, save);
0060 //      Defaults: save=0
0061 //
0062 //             For plotting correction factors
0063 //  PlotHistCorrFactor(infile, text, prefixF, scale, nmin, dataMC,
0064 //                    drawStatBox, iformat, save);
0065 //      Defaults: dataMC=true, drwaStatBox=false, nmin=100, iformat=0, save=0
0066 //
0067 //             For plotting (fractional) asymmetry in the correction factors
0068 //
0069 //  PlotHistCorrAsymmetry(infile, text, prefixF, iformat, save);
0070 //      Defaults: prefixF="", iformat=0, save=0
0071 //
0072 //             For plotting correction factors from upto 5 different runs
0073 //             on the same canvas
0074 //
0075 //  PlotHistCorrFactors(infile1, text1, infile2, text2, infile3, text3,
0076 //                      infile4, text4, infile5, text5, prefixF, ratio,
0077 //                      drawStatBox, nmin, dataMC, year, iformat, save)
0078 //      Defaults: ratio=false, drawStatBox=true, nmin=100, dataMC=false,
0079 //                year=2018, iformat=0, save=0
0080 //
0081 //             For plotting correction factors including systematics
0082 //  PlotHistCorrSys(infilec, conds, text, save)
0083 //      Defaults: save=0
0084 //
0085 //             For plotting uncertainties in correction factors with decreasing
0086 //             integrated luminpsoties starting from *lumi*
0087 //  PlotHistCorrLumis(infilec, conds, lumi, save)
0088 //      Defaults: save=0
0089 //
0090 //             For plotting correlation of correction factors
0091 //  PlotHistCorrRel(infile1, infile2, text1, text2, iformat1, iformat2, save)
0092 //      Defaults: iformat1=0, iformat2=0, save=0
0093 //
0094 //             For plotting difference of correction factors for a given depth
0095 //  PlotHistCorrDepth(infile1, infile2, text1, text2, depth, iformat1, iformat2,
0096 //                    save)
0097 //      Defaults: iformat1=0, iformat2=0, save=0
0098 //
0099 //             For plotting four histograms
0100 //  PlotFourHists(infile, prefix0, type, drawStatBox, normalize, save, prefix1,
0101 //                text1, prefix2, text2, prefix3, text3, prefix4, text4)
0102 //      Defaults: type=0, drawStatBox=0, normalize=false, save=0,
0103 //                prefixN="", textN=""
0104 //
0105 //            For plotting PU corrected histograms (o/p of CalibPlotCombine)
0106 //  PlotPUCorrHists(infile, prefix drawStatBox, approve, save)
0107 //      Defaults: infile = "corrfac.root", prefix = "", drawStatBox = 0,
0108 //                approve = true, save = 0
0109 //
0110 //             For plotting histograms obtained from fits to PU correction
0111 //             (o/p of CalibFitPU) for a given ieta using 2D/profile/Graphs
0112 //  PlotHistCorr(infile, prefix, text, eta, mode, drawStatBox, save)
0113 //      Defaults eta = 0 (all ieta values), mode = 1 (profile histograms),
0114 //               drawStatBox = true, save = 0
0115 //
0116 //             For plotting histograms created by CalibPlotProperties
0117 //  PlotPropertyHist(infile, prefix, text, etaMax, lumi, ener, dataMC,
0118 //           drawStatBox, save)
0119 //      Defaults etaMax = 25 (draws for eta = 1 .. etaMax), lumi = 0,
0120 //               ener = 13.0, dataMC = false,  drawStatBox = true, save = 0
0121 //
0122 //            For plotting mean response and resolution as a function of
0123 //            particle momentum
0124 //  PlotMeanError(infilest, region, resol, save, debug)
0125 //      Defaults region = 3 (overall), resol = false (response), save = 0,
0126 //               debug = false
0127 //      Format of the input file:
0128 //       # of energy points, # of types, # of regions
0129 //       Then for each type, energy point
0130 //         Type, lower and higher edge of momentum
0131 //         Mean response and its error for the 4 regions
0132 //         Width of response and uts error for the 4 regions
0133 //
0134 //  where:
0135 //  infile   (std::string)  = Name of the input ROOT file
0136 //  outfile  (std::string)  = Name of the output ROOT file
0137 //  prefix   (std::string)  = Prefix for the histogram names
0138 //  mode     (int)          = Flag to check which set of histograms to be
0139 //                            done. It has the format lthdo where each of
0140 //                            l, t,h,d,o can have a value 0 or 1 to select
0141 //                            or deselect. l,t,h,d,o for momentum range
0142 //                            60-100, 30-40, all, 20-30, 40-60 Gev (11111)
0143 //  type     (int)          = defines eta binning type (see CalibMonitor)
0144 //  append   (bool)         = Open the output in Update/Recreate mode (True)
0145 //  fiteta   (bool)         = fit the eta dependence with pol0
0146 //  iname    (int)          = choose the momentum bin (3: 40-60 GeV)
0147 //  saveAll  (bool)         = Flag to save intermediate plots (False)
0148 //  numb     (int)          = Number of eta bins (42 for -21:21)
0149 //  text     (std::string)  = Extra text to be put in the text title
0150 //  modePlot (int)          = Flag to plot E/p distribution (0);
0151 //                            <E/p> as a function of ieta (1);
0152 //                            <E/p> as a function of distance from L1 (2);
0153 //                            <E/p> as a function of number of vertex (3);
0154 //                            E/p for barrel, endcap and transition (4)
0155 //  kopt     (int)          = Option in format "hdo" where each of d, o can
0156 //                            have a value of 0 or 1 to select or deselect.
0157 //                            o>0 to carry out pol0 fit, o>1 to restrict
0158 //                            fit region between -20 & 20; d=1 to show grid;
0159 //                            h=0,1 to show plots with 2- or 1-Gaussian fit
0160 //  ieta     (int)          = specific ieta histogram to be plotted; if 0
0161 //                            histograms for all ieta's from -numb/2 to numb/2
0162 //                            will be plotted
0163 //  lumi     (double)       = Integrated luminosity of the dataset used which
0164 //                            needs to be drawn on the top of the canvas
0165 //                            along with CM energy (if lumi > 0)
0166 //  ener     (double)       = CM energy of the dataset used
0167 //  save     (int)          = if > 0 it saves the canvas as a pdf file; or
0168 //                            if < 0 it saves the canvas as a C file
0169 //  normalize(bool)         = if the histograms to be scaled to get
0170 //                            normalization to 1
0171 //  prefixF  (string)       = string to be included to the pad name
0172 //  infileN  (char*)        = name of the correction file and the corresponding
0173 //  textN    (string)         texts (if blank ignored)
0174 //  nmin     (int)          = minimum # of #ieta points needed to show the
0175 //                            fitted line
0176 //  scale    (double)       = constant scale factor applied to the factors
0177 //  ratio    (bool)         = set to show the ratio plot (false)
0178 //  drawStatBox (bool)      = set to show the statistical box (true)
0179 //  year     (int)          = Year of data taking (applicable to Data)
0180 //  infilc   (string)       = prefix of the file names of correction factors
0181 //                            (assumes file name would be the prefix followed
0182 //                            by _condX.txt where X=0 for the default version
0183 //                            and 1..conds for the variations)
0184 //  conds    (int)          = number of variations in estimating systematic
0185 //                            checks
0186 //  infilest (string)       = input file name containing the responses and
0187 //                            resolutions for barrel, transition, endcap,
0188 //                            overall regions at 5 energies using 3 methods
0189 //  region   (int)         = region to be selected: 0 = barrel, 1 = transition,
0190 //                           2 = endcap, 3 = overall (3)
0191 //  resol    (bool)        = parameter to be plotted: true = resolution,
0192 //                           false = response (false)
0193 //  iformat  (int)         = flag if it is created by standard (0) or by
0194 //                           Marina's (1) code
0195 //////////////////////////////////////////////////////////////////////////////
0196 
0197 #include <TCanvas.h>
0198 #include <TChain.h>
0199 #include <TProfile.h>
0200 #include <TF1.h>
0201 #include <TFile.h>
0202 #include <TFitResult.h>
0203 #include <TFitResultPtr.h>
0204 #include <TH1D.h>
0205 #include <TLegend.h>
0206 #include <TGraph.h>
0207 #include <TGraphErrors.h>
0208 #include <TGraphAsymmErrors.h>
0209 #include <TMath.h>
0210 #include <TPaveStats.h>
0211 #include <TPaveText.h>
0212 #include <TROOT.h>
0213 #include <TStyle.h>
0214 #include <cstdlib>
0215 #include <fstream>
0216 #include <iomanip>
0217 #include <iostream>
0218 #include <string>
0219 #include <vector>
0220 
0221 #include "CalibCorr.C"
0222 
0223 const double fitrangeFactor = 1.5;
0224 
0225 struct cfactors {
0226   int ieta, depth;
0227   double corrf, dcorr;
0228   cfactors(int ie = 0, int dp = 0, double cor = 1, double dc = 0) : ieta(ie), depth(dp), corrf(cor), dcorr(dc){};
0229 };
0230 
0231 struct results {
0232   double mean, errmean, width, errwidth;
0233   results(double v1 = 0, double er1 = 0, double v2 = 0, double er2 = 0)
0234       : mean(v1), errmean(er1), width(v2), errwidth(er2){};
0235 };
0236 
0237 std::pair<double, double> GetMean(TH1D* hist, double xmin, double xmax, double& rms) {
0238   double mean(0), err(0), wt(0);
0239   rms = 0;
0240   for (int i = 1; i <= hist->GetNbinsX(); ++i) {
0241     double xlow = hist->GetBinLowEdge(i);
0242     double xhigh = xlow + hist->GetBinWidth(i);
0243     if ((xlow >= xmin) && (xhigh <= xmax)) {
0244       double cont = hist->GetBinContent(i);
0245       double valu = 0.5 * (xlow + xhigh);
0246       wt += cont;
0247       mean += (valu * cont);
0248       rms += (valu * valu * cont);
0249     }
0250   }
0251   if (wt > 0) {
0252     mean /= wt;
0253     rms /= wt;
0254     err = std::sqrt((rms - mean * mean) / wt);
0255   }
0256   return std::pair<double, double>(mean, err);
0257 }
0258 
0259 std::pair<double, double> GetWidth(TH1D* hist, double xmin, double xmax) {
0260   double mean(0), mom2(0), rms(0), err(0), wt(0);
0261   for (int i = 1; i <= hist->GetNbinsX(); ++i) {
0262     double xlow = hist->GetBinLowEdge(i);
0263     double xhigh = xlow + hist->GetBinWidth(i);
0264     if ((xlow >= xmin) && (xhigh <= xmax)) {
0265       double cont = hist->GetBinContent(i);
0266       double valu = 0.5 * (xlow + xhigh);
0267       wt += cont;
0268       mean += (valu * cont);
0269       mom2 += (valu * valu * cont);
0270     }
0271   }
0272   if (wt > 0) {
0273     mean /= wt;
0274     mom2 /= wt;
0275     rms = std::sqrt(mom2 - mean * mean);
0276     err = rms / std::sqrt(2 * wt);
0277   }
0278   return std::pair<double, double>(rms, err);
0279 }
0280 
0281 Double_t langaufun(Double_t* x, Double_t* par) {
0282   //Fit parameters:
0283   //par[0]=Most Probable (MP, location) parameter of Landau density
0284   //par[1]=Total area (integral -inf to inf, normalization constant)
0285   //par[2]=Width (sigma) of convoluted Gaussian function
0286   //
0287   //In the Landau distribution (represented by the CERNLIB approximation),
0288   //the maximum is located at x=-0.22278298 with the location parameter=0.
0289   //This shift is corrected within this function, so that the actual
0290   //maximum is identical to the MP parameter.
0291 
0292   // Numeric constants
0293   Double_t invsq2pi = 0.3989422804014;  // (2 pi)^(-1/2)
0294   Double_t mpshift = -0.22278298;       // Landau maximum location
0295 
0296   // Control constants
0297   Double_t np = 100.0;  // number of convolution steps
0298   Double_t sc = 5.0;    // convolution extends to +-sc Gaussian sigmas
0299 
0300   // Variables
0301   Double_t xx;
0302   Double_t mpc;
0303   Double_t fland;
0304   Double_t sum = 0.0;
0305   Double_t xlow, xupp;
0306   Double_t step;
0307   Double_t i;
0308   Double_t par0 = 0.2;
0309 
0310   // MP shift correction
0311   mpc = par[0] - mpshift * par0 * par[0];
0312 
0313   // Range of convolution integral
0314   xlow = x[0] - sc * par[2];
0315   xupp = x[0] + sc * par[2];
0316 
0317   step = (xupp - xlow) / np;
0318 
0319   // Convolution integral of Landau and Gaussian by sum
0320   for (i = 1.0; i <= np / 2; i++) {
0321     xx = xlow + (i - .5) * step;
0322     fland = TMath::Landau(xx, mpc, par0 * par[0], kTRUE);  // / par[0];
0323     sum += fland * TMath::Gaus(x[0], xx, par[2]);
0324 
0325     xx = xupp - (i - .5) * step;
0326     fland = TMath::Landau(xx, mpc, par0 * par[0], kTRUE);  // / par[0];
0327     sum += fland * TMath::Gaus(x[0], xx, par[2]);
0328   }
0329 
0330   return (par[1] * step * sum * invsq2pi / par[2]);
0331 }
0332 
0333 Double_t doubleGauss(Double_t* x, Double_t* par) {
0334   double x1 = x[0] - par[1];
0335   double sig1 = par[2];
0336   double x2 = x[0] - par[4];
0337   double sig2 = par[5];
0338   double yval =
0339       (par[0] * std::exp(-0.5 * (x1 / sig1) * (x1 / sig1)) + par[3] * std::exp(-0.5 * (x2 / sig2) * (x2 / sig2)));
0340   return yval;
0341 }
0342 
0343 TFitResultPtr functionFit(TH1D* hist, double* fitrange, double* startvalues, double* parlimitslo, double* parlimitshi) {
0344   char FunName[100];
0345   sprintf(FunName, "Fitfcn_%s", hist->GetName());
0346   TF1* ffitold = (TF1*)gROOT->GetListOfFunctions()->FindObject(FunName);
0347   if (ffitold)
0348     delete ffitold;
0349 
0350   int npar(6);
0351   TF1* ffit = new TF1(FunName, doubleGauss, fitrange[0], fitrange[1], npar);
0352   ffit->SetParameters(startvalues);
0353   ffit->SetLineColor(kBlue);
0354   ffit->SetParNames("Area1", "Mean1", "Width1", "Area2", "Mean2", "Width2");
0355   for (int i = 0; i < npar; i++)
0356     ffit->SetParLimits(i, parlimitslo[i], parlimitshi[i]);
0357   TFitResultPtr Fit = hist->Fit(FunName, "QRWLS");
0358   return Fit;
0359 }
0360 
0361 std::pair<double, double> fitLanGau(TH1D* hist, bool debug) {
0362   double rms;
0363   std::pair<double, double> mrms = GetMean(hist, 0.005, 0.25, rms);
0364   double mean = mrms.first;
0365   double LowEdge = 0.005, HighEdge = mean + 3 * rms;
0366   TFitResultPtr Fit1 = hist->Fit("gaus", "+0wwqRS", "", LowEdge, HighEdge);
0367   if (debug)
0368     std::cout << hist->GetName() << " 0 " << Fit1->Value(0) << " 1 " << Fit1->Value(1) << " 2 " << Fit1->Value(2)
0369               << std::endl;
0370   double startvalues[3];
0371   startvalues[0] = Fit1->Value(1);
0372   startvalues[1] = Fit1->Value(0);
0373   startvalues[2] = Fit1->Value(2);
0374 
0375   char FunName[100];
0376   sprintf(FunName, "Fitfcn_%s", hist->GetName());
0377   TF1* ffitold = (TF1*)gROOT->GetListOfFunctions()->FindObject(FunName);
0378   if (ffitold)
0379     delete ffitold;
0380 
0381   TF1* ffit = new TF1(FunName, langaufun, LowEdge, HighEdge, 3);
0382   ffit->SetParameters(startvalues);
0383   ffit->SetParNames("MP", "Area", "GSigma");
0384   TFitResultPtr Fit2 = hist->Fit(FunName, "QRWLS");
0385   double value = Fit2->Value(0);
0386   double error = Fit2->FitResult::Error(0);
0387   if (debug)
0388     std::cout << hist->GetName() << " 0 " << Fit2->Value(0) << " 1 " << Fit2->Value(1) << " 2 " << Fit2->Value(2)
0389               << std::endl;
0390   return std::pair<double, double>(value, error);
0391 }
0392 
0393 results fitTwoGauss(TH1D* hist, bool debug) {
0394   double rms;
0395   std::pair<double, double> mrms = GetMean(hist, 0.2, 2.0, rms);
0396   double mean = mrms.first;
0397   double LowEdge = mean - fitrangeFactor * rms;
0398   double HighEdge = mean + fitrangeFactor * rms;
0399   if (LowEdge < 0.15)
0400     LowEdge = 0.15;
0401   std::string option = (hist->GetEntries() > 100) ? "QRS" : "QRWLS";
0402   TF1* g1 = new TF1("g1", "gaus", LowEdge, HighEdge);
0403   g1->SetLineColor(kGreen);
0404   TFitResultPtr Fit = hist->Fit(g1, option.c_str(), "");
0405 
0406   if (debug) {
0407     for (int k = 0; k < 3; ++k)
0408       std::cout << "Initial Parameter[" << k << "] = " << Fit->Value(k) << " +- " << Fit->FitResult::Error(k)
0409                 << std::endl;
0410   }
0411   double startvalues[6], fitrange[2], lowValue[6], highValue[6];
0412   startvalues[0] = Fit->Value(0);
0413   lowValue[0] = 0.5 * startvalues[0];
0414   highValue[0] = 2. * startvalues[0];
0415   startvalues[1] = Fit->Value(1);
0416   lowValue[1] = 0.5 * startvalues[1];
0417   highValue[1] = 2. * startvalues[1];
0418   startvalues[2] = Fit->Value(2);
0419   lowValue[2] = 0.5 * startvalues[2];
0420   highValue[2] = 2. * startvalues[2];
0421   startvalues[3] = 0.1 * Fit->Value(0);
0422   lowValue[3] = 0.0;
0423   highValue[3] = 10. * startvalues[3];
0424   startvalues[4] = Fit->Value(1);
0425   lowValue[4] = 0.5 * startvalues[4];
0426   highValue[4] = 2. * startvalues[4];
0427   startvalues[5] = 2.0 * Fit->Value(2);
0428   lowValue[5] = 0.5 * startvalues[5];
0429   highValue[5] = 100. * startvalues[5];
0430   //fitrange[0] = mean - 2.0*rms; fitrange[1] = mean + 2.0*rms;
0431   fitrange[0] = Fit->Value(1) - fitrangeFactor * Fit->Value(2);
0432   fitrange[1] = Fit->Value(1) + fitrangeFactor * Fit->Value(2);
0433   TFitResultPtr Fitfun = functionFit(hist, fitrange, startvalues, lowValue, highValue);
0434   double wt1 = (Fitfun->Value(0)) * (Fitfun->Value(2));
0435   double value1 = Fitfun->Value(1);
0436   double error1 = Fitfun->FitResult::Error(1);
0437   double wval1 = Fitfun->Value(2);
0438   double werr1 = Fitfun->FitResult::Error(2);
0439   double wt2 = (Fitfun->Value(3)) * (Fitfun->Value(5));
0440   double value2 = Fitfun->Value(4);
0441   double error2 = Fitfun->FitResult::Error(4);
0442   double wval2 = Fitfun->Value(5);
0443   double werr2 = Fitfun->FitResult::Error(5);
0444   double value = (wt1 * value1 + wt2 * value2) / (wt1 + wt2);
0445   double wval = (wt1 * wval1 + wt2 * wval2) / (wt1 + wt2);
0446   double error = (sqrt((wt1 * error1) * (wt1 * error1) + (wt2 * error2) * (wt2 * error2)) / (wt1 + wt2));
0447   double werror = (sqrt((wt1 * werr1) * (wt1 * werr1) + (wt2 * werr2) * (wt2 * werr2)) / (wt1 + wt2));
0448   std::cout << hist->GetName() << " Fit " << value << "+-" << error << " width " << wval << " +- " << werror
0449             << " First  " << value1 << "+-" << error1 << ":" << wval1 << "+-" << werr1 << ":" << wt1 << " Second "
0450             << value2 << "+-" << error2 << ":" << wval2 << "+-" << werr2 << ":" << wt2 << std::endl;
0451   if (debug) {
0452     for (int k = 0; k < 6; ++k)
0453       std::cout << hist->GetName() << ":Parameter[" << k << "] = " << Fitfun->Value(k) << " +- "
0454                 << Fitfun->FitResult::Error(k) << std::endl;
0455   }
0456   return results(value, error, wval, werror);
0457 }
0458 
0459 results fitOneGauss(TH1D* hist, bool fitTwice, bool debug) {
0460   double rms;
0461   std::pair<double, double> mrms = GetMean(hist, 0.2, 2.0, rms);
0462   double mean = mrms.first;
0463   double LowEdge = ((mean - fitrangeFactor * rms) < 0.5) ? 0.5 : (mean - fitrangeFactor * rms);
0464   double HighEdge = (mean + fitrangeFactor * rms);
0465   if (debug)
0466     std::cout << hist->GetName() << " Mean " << mean << " RMS " << rms << " Range " << LowEdge << ":" << HighEdge
0467               << "\n";
0468   std::string option = (hist->GetEntries() > 100) ? "QRS" : "QRWLS";
0469   TF1* g1 = new TF1("g1", "gaus", LowEdge, HighEdge);
0470   g1->SetLineColor(kGreen);
0471   TFitResultPtr Fit1 = hist->Fit(g1, option.c_str(), "");
0472   double value = Fit1->Value(1);
0473   double error = Fit1->FitResult::Error(1);
0474   double width = Fit1->Value(2);
0475   double werror = Fit1->FitResult::Error(2);
0476   if (fitTwice) {
0477     LowEdge = Fit1->Value(1) - fitrangeFactor * Fit1->Value(2);
0478     HighEdge = Fit1->Value(1) + fitrangeFactor * Fit1->Value(2);
0479     if (LowEdge < 0.5)
0480       LowEdge = 0.5;
0481     if (HighEdge > 5.0)
0482       HighEdge = 5.0;
0483     if (debug)
0484       std::cout << " Range for second Fit " << LowEdge << ":" << HighEdge << std::endl;
0485     TFitResultPtr Fit = hist->Fit("gaus", option.c_str(), "", LowEdge, HighEdge);
0486     value = Fit->Value(1);
0487     error = Fit->FitResult::Error(1);
0488     width = Fit->Value(2);
0489     werror = Fit->FitResult::Error(2);
0490   }
0491 
0492   std::pair<double, double> meaner = GetMean(hist, 0.2, 2.0, rms);
0493   std::pair<double, double> rmserr = GetWidth(hist, 0.2, 2.0);
0494   if (debug) {
0495     std::cout << "Fit " << value << ":" << error << ":" << hist->GetMeanError() << " Mean " << meaner.first << ":"
0496               << meaner.second << " Width " << rmserr.first << ":" << rmserr.second;
0497   }
0498   double minvalue(0.30);
0499   if (value < minvalue || value > 2.0 || error > 0.5) {
0500     value = meaner.first;
0501     error = meaner.second;
0502     width = rmserr.first;
0503     werror = rmserr.second;
0504   }
0505   if (debug) {
0506     std::cout << " Final " << value << "+-" << error << ":" << width << "+-" << werror << std::endl;
0507   }
0508   return results(value, error, width, werror);
0509 }
0510 
0511 void readCorrFactors(char* infile,
0512                      double scale,
0513                      std::map<int, cfactors>& cfacs,
0514                      int& etamin,
0515                      int& etamax,
0516                      int& maxdepth,
0517                      int iformat = 0,
0518                      bool debug = false) {
0519   cfacs.clear();
0520   std::ifstream fInput(infile);
0521   if (!fInput.good()) {
0522     std::cout << "Cannot open file " << infile << std::endl;
0523   } else {
0524     char buffer[1024];
0525     unsigned int all(0), good(0);
0526     while (fInput.getline(buffer, 1024)) {
0527       ++all;
0528       if (buffer[0] == '#')
0529         continue;  //ignore comment
0530       std::vector<std::string> items = splitString(std::string(buffer));
0531       if (items.size() != 5) {
0532         std::cout << "Ignore  line: " << buffer << std::endl;
0533       } else {
0534         ++good;
0535         int ieta = (iformat == 1) ? std::atoi(items[0].c_str()) : std::atoi(items[1].c_str());
0536         int depth = (iformat == 1) ? std::atoi(items[1].c_str()) : std::atoi(items[2].c_str());
0537         float corrf = std::atof(items[3].c_str());
0538         float dcorr = (iformat == 1) ? (0.02 * corrf) : std::atof(items[4].c_str());
0539         cfactors cfac(ieta, depth, scale * corrf, scale * dcorr);
0540         int detId = (iformat == 1) ? repackId(items[2], ieta, depth) : repackId(ieta, depth);
0541         cfacs[detId] = cfactors(ieta, depth, corrf, dcorr);
0542         if (ieta > etamax)
0543           etamax = ieta;
0544         if (ieta < etamin)
0545           etamin = ieta;
0546         if (depth > maxdepth)
0547           maxdepth = depth;
0548       }
0549     }
0550     fInput.close();
0551     std::cout << "Reads total of " << all << " and " << good << " good records"
0552               << " from " << infile << std::endl;
0553   }
0554   if (debug) {
0555     unsigned k(0);
0556     std::cout << "Eta Range " << etamin << ":" << etamax << " Max Depth " << maxdepth << std::endl;
0557     for (std::map<int, cfactors>::const_iterator itr = cfacs.begin(); itr != cfacs.end(); ++itr, ++k)
0558       std::cout << "[" << k << "] " << std::hex << itr->first << std::dec << ": " << (itr->second).ieta << " "
0559                 << (itr->second).depth << " " << (itr->second).corrf << " " << (itr->second).dcorr << std::endl;
0560   }
0561 }
0562 
0563 void FitHistStandard(std::string infile,
0564                      std::string outfile,
0565                      std::string prefix,
0566                      int mode = 11111,
0567                      int type = 0,
0568                      bool append = true,
0569                      bool saveAll = false,
0570                      bool debug = false) {
0571   int iname[6] = {0, 1, 2, 3, 4, 5};
0572   int checkmode[6] = {1000, 10, 10000, 1, 100000, 100};
0573   double xbin0[9] = {-21.0, -16.0, -12.0, -6.0, 0.0, 6.0, 12.0, 16.0, 21.0};
0574   double xbins[11] = {-25.0, -20.0, -15.0, -10.0, -5.0, 0.0, 5.0, 10.0, 15.0, 20.0, 25.0};
0575   double vbins[6] = {0.0, 7.0, 10.0, 13.0, 16.0, 50.0};
0576   double dlbins[9] = {0.0, 0.10, 0.20, 0.50, 1.0, 2.0, 2.5, 3.0, 10.0};
0577   std::string sname[4] = {"ratio", "etaR", "dl1R", "nvxR"};
0578   std::string lname[4] = {"Z", "E", "L", "V"};
0579   std::string wname[4] = {"W", "F", "M", "X"};
0580   std::string xname[4] = {"i#eta", "i#eta", "d_{L1}", "# Vertex"};
0581   int numb[4] = {10, 8, 8, 5};
0582 
0583   if (type == 0) {
0584     numb[0] = 8;
0585     for (int i = 0; i < 9; ++i)
0586       xbins[i] = xbin0[i];
0587   }
0588   TFile* file = new TFile(infile.c_str());
0589   std::vector<TH1D*> hists;
0590   char name[100], namw[100];
0591   if (file != nullptr) {
0592     for (int m1 = 0; m1 < 4; ++m1) {
0593       for (int m2 = 0; m2 < 6; ++m2) {
0594         sprintf(name, "%s%s%d0", prefix.c_str(), sname[m1].c_str(), iname[m2]);
0595         TH1D* hist0 = (TH1D*)file->FindObjectAny(name);
0596         bool ok = ((hist0 != nullptr) && (hist0->GetEntries() > 25));
0597         if ((mode / checkmode[m2]) % 10 > 0 && ok) {
0598           TH1D *histo(0), *histw(0);
0599           sprintf(name, "%s%s%d", prefix.c_str(), lname[m1].c_str(), iname[m2]);
0600           sprintf(namw, "%s%s%d", prefix.c_str(), wname[m1].c_str(), iname[m2]);
0601           if (m1 <= 1) {
0602             histo = new TH1D(name, hist0->GetTitle(), numb[m1], xbins);
0603             histw = new TH1D(namw, hist0->GetTitle(), numb[m1], xbins);
0604           } else if (m1 == 2) {
0605             histo = new TH1D(name, hist0->GetTitle(), numb[m1], dlbins);
0606             histw = new TH1D(namw, hist0->GetTitle(), numb[m1], dlbins);
0607           } else {
0608             histo = new TH1D(name, hist0->GetTitle(), numb[m1], vbins);
0609             histw = new TH1D(namw, hist0->GetTitle(), numb[m1], vbins);
0610           }
0611           int jmin(numb[m1]), jmax(0);
0612           for (int j = 0; j <= numb[m1]; ++j) {
0613             sprintf(name, "%s%s%d%d", prefix.c_str(), sname[m1].c_str(), iname[m2], j);
0614             TH1D* hist1 = (TH1D*)file->FindObjectAny(name);
0615             TH1D* hist = (TH1D*)hist1->Clone();
0616             double value(0), error(0), width(0), werror(0);
0617             if (hist->GetEntries() > 0) {
0618               value = hist->GetMean();
0619               error = hist->GetRMS();
0620               std::pair<double, double> rmserr = GetWidth(hist, 0.2, 2.0);
0621               width = rmserr.first;
0622               werror = rmserr.second;
0623             }
0624             if (hist->GetEntries() > 4) {
0625               bool flag = (j == 0) ? true : false;
0626               results meaner = fitOneGauss(hist, flag, debug);
0627               value = meaner.mean;
0628               error = meaner.errmean;
0629               width = meaner.width;
0630               werror = meaner.errwidth;
0631               if (j != 0) {
0632                 if (j < jmin)
0633                   jmin = j;
0634                 if (j > jmax)
0635                   jmax = j;
0636               }
0637             }
0638             if (j == 0) {
0639               hists.push_back(hist);
0640             } else {
0641               double wbyv = width / value;
0642               double wverr = wbyv * std::sqrt((werror * werror) / (width * width) + (error * error) / (value * value));
0643               if (saveAll)
0644                 hists.push_back(hist);
0645               histo->SetBinContent(j, value);
0646               histo->SetBinError(j, error);
0647               histw->SetBinContent(j, wbyv);
0648               histw->SetBinError(j, wverr);
0649             }
0650           }
0651           if (histo->GetEntries() > 2) {
0652             double LowEdge = histo->GetBinLowEdge(jmin);
0653             double HighEdge = histo->GetBinLowEdge(jmax) + histo->GetBinWidth(jmax);
0654             TFitResultPtr Fit = histo->Fit("pol0", "+QRWLS", "", LowEdge, HighEdge);
0655             if (debug) {
0656               std::cout << "Fit to Pol0: " << Fit->Value(0) << " +- " << Fit->FitResult::Error(0) << std::endl;
0657             }
0658             histo->GetXaxis()->SetTitle(xname[m1].c_str());
0659             histo->GetYaxis()->SetTitle("MPV(E_{HCAL}/(p-E_{ECAL}))");
0660             histo->GetYaxis()->SetRangeUser(0.4, 1.6);
0661             histw->GetXaxis()->SetTitle(xname[m1].c_str());
0662             histw->GetYaxis()->SetTitle("Width/MPV of (E_{HCAL}/(p-E_{ECAL}))");
0663             histw->GetYaxis()->SetRangeUser(0.0, 0.5);
0664           }
0665           hists.push_back(histo);
0666           hists.push_back(histw);
0667         }
0668       }
0669     }
0670     TFile* theFile(0);
0671     if (append) {
0672       theFile = new TFile(outfile.c_str(), "UPDATE");
0673     } else {
0674       theFile = new TFile(outfile.c_str(), "RECREATE");
0675     }
0676     theFile->cd();
0677     for (unsigned int i = 0; i < hists.size(); ++i) {
0678       TH1D* hnew = (TH1D*)hists[i]->Clone();
0679       hnew->Write();
0680     }
0681     theFile->Close();
0682     file->Close();
0683   }
0684 }
0685 
0686 void FitHistExtended(const char* infile,
0687                      const char* outfile,
0688                      std::string prefix,
0689                      int numb = 50,
0690                      int type = 3,
0691                      bool append = true,
0692                      bool fiteta = true,
0693                      int iname = 3,
0694                      bool debug = false) {
0695   std::string sname("ratio"), lname("Z"), wname("W"), ename("etaB");
0696   double xbins[99];
0697   double xbin[23] = {-23.0, -21.0, -19.0, -17.0, -15.0, -13.0, -11.0, -9.0, -7.0, -5.0, -3.0, 0.0,
0698                      3.0,   5.0,   7.0,   9.0,   11.0,  13.0,  15.0,  17.0, 19.0, 21.0, 23.0};
0699   if (type == 2) {
0700     numb = 22;
0701     for (int k = 0; k <= numb; ++k)
0702       xbins[k] = xbin[k];
0703   } else if (type == 1) {
0704     numb = 1;
0705     xbins[0] = -25;
0706     xbins[1] = 25;
0707   } else {
0708     int neta = numb / 2;
0709     for (int k = 0; k < neta; ++k) {
0710       xbins[k] = (k - neta) - 0.5;
0711       xbins[numb - k] = (neta - k) + 0.5;
0712     }
0713     xbins[neta] = 0;
0714   }
0715   if (debug) {
0716     for (int k = 0; k <= numb; ++k)
0717       std::cout << " " << xbins[k];
0718     std::cout << std::endl;
0719   }
0720   TFile* file = new TFile(infile);
0721   std::vector<TH1D*> hists;
0722   char name[200];
0723   if (debug) {
0724     std::cout << infile << " " << file << std::endl;
0725   }
0726   if (file != nullptr) {
0727     sprintf(name, "%s%s%d0", prefix.c_str(), sname.c_str(), iname);
0728     TH1D* hist0 = (TH1D*)file->FindObjectAny(name);
0729     bool ok = (hist0 != nullptr);
0730     if (debug) {
0731       std::cout << name << " Pointer " << hist0 << " " << ok << std::endl;
0732     }
0733     if (ok) {
0734       TH1D *histo(0), *histw(0);
0735       if (numb > 0) {
0736         sprintf(name, "%s%s%d", prefix.c_str(), lname.c_str(), iname);
0737         histo = new TH1D(name, hist0->GetTitle(), numb, xbins);
0738         sprintf(name, "%s%s%d", prefix.c_str(), wname.c_str(), iname);
0739         histw = new TH1D(name, hist0->GetTitle(), numb, xbins);
0740         if (debug)
0741           std::cout << name << " " << histo->GetNbinsX() << std::endl;
0742       }
0743       if (hist0->GetEntries() > 10) {
0744         double rms;
0745         results meaner0 = fitOneGauss(hist0, true, debug);
0746         std::pair<double, double> meaner1 = GetMean(hist0, 0.2, 2.0, rms);
0747         std::pair<double, double> meaner2 = GetWidth(hist0, 0.2, 2.0);
0748         if (debug) {
0749           std::cout << "Fit " << meaner0.mean << ":" << meaner0.errmean << " Mean1 " << hist0->GetMean() << ":"
0750                     << hist0->GetMeanError() << " Mean2 " << meaner1.first << ":" << meaner1.second << " Width "
0751                     << meaner2.first << ":" << meaner2.second << std::endl;
0752         }
0753       }
0754       int nv1(100), nv2(0);
0755       int jmin(numb), jmax(0);
0756       for (int j = 0; j <= numb; ++j) {
0757         sprintf(name, "%s%s%d%d", prefix.c_str(), sname.c_str(), iname, j);
0758         TH1D* hist1 = (TH1D*)file->FindObjectAny(name);
0759         if (debug) {
0760           std::cout << "Get Histogram for " << name << " at " << hist1 << std::endl;
0761         }
0762         double value(0), error(0), total(0), width(0), werror(0);
0763         if (hist1 == nullptr) {
0764           value = 1.0;
0765         } else {
0766           TH1D* hist = (TH1D*)hist1->Clone();
0767           if (hist->GetEntries() > 0) {
0768             value = hist->GetMean();
0769             error = hist->GetRMS();
0770             for (int i = 1; i <= hist->GetNbinsX(); ++i)
0771               total += hist->GetBinContent(i);
0772             std::pair<double, double> rmserr = GetWidth(hist, 0.2, 2.0);
0773             width = rmserr.first;
0774             werror = rmserr.second;
0775           }
0776           if (total > 4) {
0777             if (nv1 > j)
0778               nv1 = j;
0779             if (nv2 < j)
0780               nv2 = j;
0781             if (j == 0) {
0782               sprintf(name, "%sOne", hist1->GetName());
0783               TH1D* hist2 = (TH1D*)hist1->Clone(name);
0784               fitOneGauss(hist2, true, debug);
0785               hists.push_back(hist2);
0786               results meaner = fitOneGauss(hist, true, debug);
0787               value = meaner.mean;
0788               error = meaner.errmean;
0789               width = meaner.width;
0790               werror = meaner.errwidth;
0791             } else {
0792               results meaner = fitOneGauss(hist, true, debug);
0793               value = meaner.mean;
0794               error = meaner.errmean;
0795               width = meaner.width;
0796               werror = meaner.errwidth;
0797             }
0798             if (j != 0) {
0799               if (j < jmin)
0800                 jmin = j;
0801               if (j > jmax)
0802                 jmax = j;
0803             }
0804           }
0805           hists.push_back(hist);
0806         }
0807         if (debug) {
0808           std::cout << "Hist " << j << " Value " << value << " +- " << error << std::endl;
0809         }
0810         if (j != 0) {
0811           double wbyv = width / value;
0812           double wverr = wbyv * std::sqrt((werror * werror) / (width * width) + (error * error) / (value * value));
0813           histo->SetBinContent(j, value);
0814           histo->SetBinError(j, error);
0815           histw->SetBinContent(j, wbyv);
0816           histw->SetBinError(j, wverr);
0817         }
0818       }
0819       if (histo != nullptr) {
0820         if (histo->GetEntries() > 2 && fiteta) {
0821           if (debug) {
0822             std::cout << "Jmin/max " << jmin << ":" << jmax << ":" << histo->GetNbinsX() << std::endl;
0823           }
0824           double LowEdge = histo->GetBinLowEdge(jmin);
0825           double HighEdge = histo->GetBinLowEdge(jmax) + histo->GetBinWidth(jmax);
0826           TFitResultPtr Fit = histo->Fit("pol0", "+QRWLS", "", LowEdge, HighEdge);
0827           if (debug) {
0828             std::cout << "Fit to Pol0: " << Fit->Value(0) << " +- " << Fit->FitResult::Error(0) << " in range " << nv1
0829                       << ":" << xbins[nv1] << ":" << nv2 << ":" << xbins[nv2] << std::endl;
0830           }
0831           histo->GetXaxis()->SetTitle("i#eta");
0832           histo->GetYaxis()->SetTitle("MPV(E_{HCAL}/(p-E_{ECAL}))");
0833           histo->GetYaxis()->SetRangeUser(0.4, 1.6);
0834           histw->GetXaxis()->SetTitle("i#eta");
0835           histw->GetYaxis()->SetTitle("MPV/Width(E_{HCAL}/(p-E_{ECAL}))");
0836           histw->GetYaxis()->SetRangeUser(0.0, 0.5);
0837         }
0838         hists.push_back(histo);
0839         hists.push_back(histw);
0840       } else {
0841         hists.push_back(hist0);
0842       }
0843 
0844       //Barrel,Endcap
0845       for (int j = 1; j <= 4; ++j) {
0846         sprintf(name, "%s%s%d%d", prefix.c_str(), ename.c_str(), iname, j);
0847         TH1D* hist1 = (TH1D*)file->FindObjectAny(name);
0848         if (debug) {
0849           std::cout << "Get Histogram for " << name << " at " << hist1 << std::endl;
0850         }
0851         if (hist1 != nullptr) {
0852           TH1D* hist = (TH1D*)hist1->Clone();
0853           double value(0), error(0), total(0), width(0), werror(0);
0854           if (hist->GetEntries() > 0) {
0855             value = hist->GetMean();
0856             error = hist->GetRMS();
0857             for (int i = 1; i <= hist->GetNbinsX(); ++i)
0858               total += hist->GetBinContent(i);
0859           }
0860           if (total > 4) {
0861             sprintf(name, "%sOne", hist1->GetName());
0862             TH1D* hist2 = (TH1D*)hist1->Clone(name);
0863             results meanerr = fitOneGauss(hist2, true, debug);
0864             value = meanerr.mean;
0865             error = meanerr.errmean;
0866             width = meanerr.width;
0867             werror = meanerr.errwidth;
0868             double wbyv = width / value;
0869             double wverr = wbyv * std::sqrt((werror * werror) / (width * width) + (error * error) / (value * value));
0870             std::cout << hist2->GetName() << " MPV " << value << " +- " << error << " Width " << width << " +- "
0871                       << werror << " W/M " << wbyv << " +- " << wverr << std::endl;
0872             hists.push_back(hist2);
0873             if (hist1->GetBinLowEdge(1) < 0.1) {
0874               sprintf(name, "%sTwo", hist1->GetName());
0875               TH1D* hist3 = (TH1D*)hist1->Clone(name);
0876               fitLanGau(hist3, debug);
0877               hists.push_back(hist3);
0878             }
0879             //            results meaner0 = fitTwoGauss(hist, debug);
0880             results meaner0 = fitOneGauss(hist, true, debug);
0881             value = meaner0.mean;
0882             error = meaner0.errmean;
0883             double rms;
0884             std::pair<double, double> meaner = GetMean(hist, 0.2, 2.0, rms);
0885             if (debug) {
0886               std::cout << "Fit " << value << ":" << error << ":" << hist->GetMeanError() << " Mean " << meaner.first
0887                         << ":" << meaner.second << std::endl;
0888             }
0889           }
0890           hists.push_back(hist);
0891         }
0892       }
0893     }
0894     TFile* theFile(0);
0895     if (append) {
0896       if (debug) {
0897         std::cout << "Open file " << outfile << " in append mode" << std::endl;
0898       }
0899       theFile = new TFile(outfile, "UPDATE");
0900     } else {
0901       if (debug) {
0902         std::cout << "Open file " << outfile << " in recreate mode" << std::endl;
0903       }
0904       theFile = new TFile(outfile, "RECREATE");
0905     }
0906 
0907     theFile->cd();
0908     for (unsigned int i = 0; i < hists.size(); ++i) {
0909       TH1D* hnew = (TH1D*)hists[i]->Clone();
0910       if (debug) {
0911         std::cout << "Write Histogram " << hnew->GetTitle() << std::endl;
0912       }
0913       hnew->Write();
0914     }
0915     theFile->Close();
0916     file->Close();
0917   }
0918 }
0919 
0920 void FitHistRBX(const char* infile, const char* outfile, std::string prefix, bool append = true, int iname = 3) {
0921   std::string sname("RBX"), lname("R");
0922   int numb(18);
0923   bool debug(false);
0924   char name[200];
0925 
0926   TFile* file = new TFile(infile);
0927   std::vector<TH1D*> hists;
0928   sprintf(name, "%s%s%d", prefix.c_str(), lname.c_str(), iname);
0929   TH1D* histo = new TH1D(name, "", numb, 0, numb);
0930   if (debug) {
0931     std::cout << infile << " " << file << std::endl;
0932   }
0933   if (file != nullptr) {
0934     for (int j = 0; j < numb; ++j) {
0935       sprintf(name, "%s%s%d%d", prefix.c_str(), sname.c_str(), iname, j + 1);
0936       TH1D* hist1 = (TH1D*)file->FindObjectAny(name);
0937       if (debug) {
0938         std::cout << "Get Histogram for " << name << " at " << hist1 << std::endl;
0939       }
0940       TH1D* hist = (TH1D*)hist1->Clone();
0941       double value(0), error(0), total(0);
0942       if (hist->GetEntries() > 0) {
0943         value = hist->GetMean();
0944         error = hist->GetRMS();
0945         for (int i = 1; i <= hist->GetNbinsX(); ++i)
0946           total += hist->GetBinContent(i);
0947       }
0948       if (total > 4) {
0949         results meaner = fitOneGauss(hist, false, debug);
0950         value = meaner.mean;
0951         error = meaner.errmean;
0952       }
0953       hists.push_back(hist);
0954       histo->SetBinContent(j + 1, value);
0955       histo->SetBinError(j + 1, error);
0956     }
0957     histo->GetXaxis()->SetTitle("RBX #");
0958     histo->GetYaxis()->SetTitle("MPV(E_{HCAL}/(p-E_{ECAL}))");
0959     histo->GetYaxis()->SetRangeUser(0.75, 1.20);
0960     hists.push_back(histo);
0961 
0962     TFile* theFile(0);
0963     if (append) {
0964       if (debug) {
0965         std::cout << "Open file " << outfile << " in append mode" << std::endl;
0966       }
0967       theFile = new TFile(outfile, "UPDATE");
0968     } else {
0969       if (debug) {
0970         std::cout << "Open file " << outfile << " in recreate mode" << std::endl;
0971       }
0972       theFile = new TFile(outfile, "RECREATE");
0973     }
0974 
0975     theFile->cd();
0976     for (unsigned int i = 0; i < hists.size(); ++i) {
0977       TH1D* hnew = (TH1D*)hists[i]->Clone();
0978       if (debug) {
0979         std::cout << "Write Histogram " << hnew->GetTitle() << std::endl;
0980       }
0981       hnew->Write();
0982     }
0983     theFile->Close();
0984     file->Close();
0985   }
0986 }
0987 
0988 void PlotHist(const char* infile,
0989               std::string prefix,
0990               std::string text,
0991               int mode = 4,
0992               int kopt = 100,
0993               double lumi = 0,
0994               double ener = 13.0,
0995               bool dataMC = false,
0996               bool drawStatBox = true,
0997               int save = 0) {
0998   std::string name0[6] = {"ratio00", "ratio10", "ratio20", "ratio30", "ratio40", "ratio50"};
0999   std::string name1[5] = {"Z0", "Z1", "Z2", "Z3", "Z4"};
1000   std::string name2[5] = {"L0", "L1", "L2", "L3", "L4"};
1001   std::string name3[5] = {"V0", "V1", "V2", "V3", "V4"};
1002   std::string name4[20] = {"etaB41", "etaB42", "etaB43", "etaB44", "etaB31", "etaB32", "etaB33",
1003                            "etaB34", "etaB21", "etaB22", "etaB23", "etaB24", "etaB11", "etaB12",
1004                            "etaB13", "etaB14", "etaB01", "etaB02", "etaB03", "etaB04"};
1005   std::string name5[5] = {"W0", "W1", "W2", "W3", "W4"};
1006   std::string title[6] = {"Tracks with p = 10:20 GeV",
1007                           "Tracks with p = 20:30 GeV",
1008                           "Tracks with p = 30:40 GeV",
1009                           "Tracks with p = 40:60 GeV",
1010                           "Tracks with p = 60:100 GeV",
1011                           "Tracks with p = 20:100 GeV"};
1012   std::string title1[20] = {"Tracks with p = 60:100 GeV (Barrel)", "Tracks with p = 60:100 GeV (Transition)",
1013                             "Tracks with p = 60:100 GeV (Endcap)", "Tracks with p = 60:100 GeV",
1014                             "Tracks with p = 40:60 GeV (Barrel)",  "Tracks with p = 40:60 GeV (Transition)",
1015                             "Tracks with p = 40:60 GeV (Endcap)",  "Tracks with p = 40:60 GeV",
1016                             "Tracks with p = 30:40 GeV (Barrel)",  "Tracks with p = 30:40 GeV (Transition)",
1017                             "Tracks with p = 30:40 GeV (Endcap)",  "Tracks with p = 30:40 GeV",
1018                             "Tracks with p = 20:30 GeV (Barrel)",  "Tracks with p = 20:30 GeV (Transition)",
1019                             "Tracks with p = 20:30 GeV (Endcap)",  "Tracks with p = 20:30 GeV",
1020                             "Tracks with p = 10:20 GeV (Barrel)",  "Tracks with p = 10:20 GeV (Transition)",
1021                             "Tracks with p = 10:20 GeV (Endcap)",  "Tracks with p = 10:20 GeV"};
1022   std::string xtitl[5] = {"E_{HCAL}/(p-E_{ECAL})", "i#eta", "d_{L1}", "# Vertex", "E_{HCAL}/(p-E_{ECAL})"};
1023   std::string ytitl[5] = {
1024       "Tracks", "MPV(E_{HCAL}/(p-E_{ECAL}))", "MPV(E_{HCAL}/(p-E_{ECAL}))", "MPV(E_{HCAL}/(p-E_{ECAL}))", "Tracks"};
1025 
1026   gStyle->SetCanvasBorderMode(0);
1027   gStyle->SetCanvasColor(kWhite);
1028   gStyle->SetPadColor(kWhite);
1029   gStyle->SetFillColor(kWhite);
1030   gStyle->SetOptTitle(0);
1031   if (mode < 0 || mode > 5)
1032     mode = 0;
1033   if (drawStatBox) {
1034     int iopt(1110);
1035     if (mode != 0)
1036       iopt = 10;
1037     gStyle->SetOptStat(iopt);
1038     gStyle->SetOptFit(1);
1039   } else {
1040     gStyle->SetOptStat(0);
1041     gStyle->SetOptFit(0);
1042   }
1043   TFile* file = new TFile(infile);
1044   TLine* line(0);
1045   char name[100], namep[100];
1046   int kmax = (mode == 4) ? 20 : (((mode < 1) && (mode > 5)) ? 6 : 5);
1047   for (int k = 0; k < kmax; ++k) {
1048     if (mode == 1) {
1049       sprintf(name, "%s%s", prefix.c_str(), name1[k].c_str());
1050     } else if (mode == 2) {
1051       sprintf(name, "%s%s", prefix.c_str(), name2[k].c_str());
1052     } else if (mode == 3) {
1053       sprintf(name, "%s%s", prefix.c_str(), name3[k].c_str());
1054     } else if (mode == 4) {
1055       if ((kopt / 100) % 10 == 0) {
1056         sprintf(name, "%s%s", prefix.c_str(), name4[k].c_str());
1057       } else if ((kopt / 100) % 10 == 2) {
1058         sprintf(name, "%s%sTwo", prefix.c_str(), name4[k].c_str());
1059       } else {
1060         sprintf(name, "%s%sOne", prefix.c_str(), name4[k].c_str());
1061       }
1062     } else if (mode == 5) {
1063       sprintf(name, "%s%s", prefix.c_str(), name5[k].c_str());
1064     } else {
1065       if ((kopt / 100) % 10 == 0) {
1066         sprintf(name, "%s%s", prefix.c_str(), name0[k].c_str());
1067       } else if ((kopt / 100) % 10 == 2) {
1068         sprintf(name, "%s%sTwo", prefix.c_str(), name0[k].c_str());
1069       } else {
1070         sprintf(name, "%s%sOne", prefix.c_str(), name0[k].c_str());
1071       }
1072     }
1073     TH1D* hist1 = (TH1D*)file->FindObjectAny(name);
1074     if (hist1 != nullptr) {
1075       TH1D* hist = (TH1D*)(hist1->Clone());
1076       double p0(1);
1077       double ymin(0.90);
1078       sprintf(namep, "c_%s", name);
1079       TCanvas* pad = new TCanvas(namep, namep, 700, 500);
1080       pad->SetRightMargin(0.10);
1081       pad->SetTopMargin(0.10);
1082       if ((kopt / 10) % 10 > 0)
1083         gPad->SetGrid();
1084       hist->GetXaxis()->SetTitleSize(0.04);
1085       hist->GetXaxis()->SetTitle(xtitl[mode].c_str());
1086       hist->GetYaxis()->SetTitle(ytitl[mode].c_str());
1087       hist->GetYaxis()->SetLabelOffset(0.005);
1088       hist->GetYaxis()->SetTitleSize(0.04);
1089       hist->GetYaxis()->SetLabelSize(0.035);
1090       hist->GetYaxis()->SetTitleOffset(1.10);
1091       if (mode == 0 || mode == 4) {
1092         if ((kopt / 100) % 10 == 2)
1093           hist->GetXaxis()->SetRangeUser(0.0, 0.30);
1094         else
1095           hist->GetXaxis()->SetRangeUser(0.25, 2.25);
1096       } else {
1097         if (mode == 5)
1098           hist->GetYaxis()->SetRangeUser(0.1, 0.50);
1099         else
1100           hist->GetYaxis()->SetRangeUser(0.8, 1.20);
1101         if (kopt % 10 > 0) {
1102           int nbin = hist->GetNbinsX();
1103           double LowEdge = (kopt % 10 == 1) ? hist->GetBinLowEdge(1) : -20;
1104           double HighEdge = (kopt % 10 == 1) ? hist->GetBinLowEdge(nbin) + hist->GetBinWidth(nbin) : 20;
1105           TFitResultPtr Fit = hist->Fit("pol0", "+QRWLS", "", LowEdge, HighEdge);
1106           p0 = Fit->Value(0);
1107         }
1108       }
1109       hist->SetMarkerStyle(20);
1110       hist->SetMarkerColor(2);
1111       hist->SetLineColor(2);
1112       hist->Draw();
1113       pad->Update();
1114       TPaveStats* st1 = (TPaveStats*)hist->GetListOfFunctions()->FindObject("stats");
1115       if (st1 != nullptr) {
1116         ymin = (mode == 0 || mode == 4) ? 0.70 : 0.80;
1117         st1->SetY1NDC(ymin);
1118         st1->SetY2NDC(0.90);
1119         st1->SetX1NDC(0.65);
1120         st1->SetX2NDC(0.90);
1121       }
1122       if (mode != 0 && mode != 4 && kopt % 10 > 0) {
1123         double xmin = hist->GetBinLowEdge(1);
1124         int nbin = hist->GetNbinsX();
1125         double xmax = hist->GetBinLowEdge(nbin) + hist->GetBinWidth(nbin);
1126         double mean(0), rms(0), total(0);
1127         int kount(0);
1128         for (int k = 2; k < nbin; ++k) {
1129           double x = hist->GetBinContent(k);
1130           double w = hist->GetBinError(k);
1131           mean += (x * w);
1132           rms += (x * x * w);
1133           total += w;
1134           ++kount;
1135         }
1136         mean /= total;
1137         rms /= total;
1138         double error = sqrt(rms - mean * mean) / total;
1139         line = new TLine(xmin, p0, xmax, p0);  //etamin,1.0,etamax,1.0);
1140         std::cout << xmin << ":" << xmax << ":" << p0 << " Mean " << nbin << ":" << kount << ":" << total << ":" << mean
1141                   << ":" << rms << ":" << error << std::endl;
1142         line->SetLineWidth(2);
1143         line->SetLineStyle(2);
1144         line->Draw("same");
1145       }
1146       pad->Update();
1147       double ymx(0.96), xmi(0.25), xmx(0.90);
1148       char txt[100];
1149       if (lumi > 0.1) {
1150         ymx = ymin - 0.005;
1151         xmi = 0.45;
1152         TPaveText* txt0 = new TPaveText(0.65, 0.91, 0.90, 0.96, "blNDC");
1153         txt0->SetFillColor(0);
1154         sprintf(txt, "%4.1f TeV %5.1f fb^{-1}", ener, lumi);
1155         txt0->AddText(txt);
1156         txt0->Draw("same");
1157       }
1158       double ymi = ymx - 0.05;
1159       TPaveText* txt1 = new TPaveText(xmi, ymi, xmx, ymx, "blNDC");
1160       txt1->SetFillColor(0);
1161       if (text == "") {
1162         if (mode == 4)
1163           sprintf(txt, "%s", title1[k].c_str());
1164         else
1165           sprintf(txt, "%s", title[k].c_str());
1166       } else {
1167         if (mode == 4)
1168           sprintf(txt, "%s (%s)", title1[k].c_str(), text.c_str());
1169         else
1170           sprintf(txt, "%s (%s)", title[k].c_str(), text.c_str());
1171       }
1172       txt1->AddText(txt);
1173       txt1->Draw("same");
1174       double xmax = (dataMC) ? 0.33 : 0.44;
1175       ymi = (lumi > 0.1) ? 0.91 : 0.84;
1176       ymx = ymi + 0.05;
1177       TPaveText* txt2 = new TPaveText(0.11, ymi, xmax, ymx, "blNDC");
1178       txt2->SetFillColor(0);
1179       if (dataMC)
1180         sprintf(txt, "CMS Preliminary");
1181       else
1182         sprintf(txt, "CMS Simulation Preliminary");
1183       txt2->AddText(txt);
1184       txt2->Draw("same");
1185       pad->Modified();
1186       pad->Update();
1187       if (save > 0) {
1188         sprintf(name, "%s.pdf", pad->GetName());
1189         pad->Print(name);
1190       } else if (save < 0) {
1191         sprintf(name, "%s.C", pad->GetName());
1192         pad->Print(name);
1193       }
1194     }
1195   }
1196 }
1197 
1198 void PlotHistEta(const char* infile,
1199                  std::string prefix,
1200                  std::string text,
1201                  int iene = 3,
1202                  int numb = 50,
1203                  int ieta = 0,
1204                  double lumi = 0,
1205                  double ener = 13.0,
1206                  bool dataMC = false,
1207                  bool drawStatBox = true,
1208                  int save = 0) {
1209   std::string name0 = "ratio";
1210   std::string title[5] = {"10:20", "20:30", "30:40", "40:60", "60:100"};
1211   std::string xtitl = "E_{HCAL}/(p-E_{ECAL})";
1212   std::string ytitl = "Tracks";
1213 
1214   gStyle->SetCanvasBorderMode(0);
1215   gStyle->SetCanvasColor(kWhite);
1216   gStyle->SetPadColor(kWhite);
1217   gStyle->SetFillColor(kWhite);
1218   gStyle->SetOptTitle(0);
1219   if (drawStatBox) {
1220     int iopt(1110);
1221     gStyle->SetOptStat(iopt);
1222     gStyle->SetOptFit(1);
1223   } else {
1224     gStyle->SetOptStat(0);
1225     gStyle->SetOptFit(0);
1226   }
1227   if (iene < 0 || iene >= 5)
1228     iene = 3;
1229   int numb2 = numb / 2;
1230   if (ieta < -numb2 || ieta > numb2)
1231     ieta = 0;
1232   int ietaMin = ((ieta == 0) ? 1 : ((ieta > 0) ? (numb2 + ieta) : (numb2 + ieta + 1)));
1233   int ietaMax = (ieta == 0) ? numb : ietaMin;
1234   TFile* file = new TFile(infile);
1235   char name[100], namep[100];
1236   for (int k = ietaMin; k <= ietaMax; ++k) {
1237     int eta = (k > numb2) ? (k - numb2) : (k - numb2 - 1);
1238     sprintf(name, "%s%s%d%d", prefix.c_str(), name0.c_str(), iene, k);
1239     TH1D* hist1 = (TH1D*)file->FindObjectAny(name);
1240     std::cout << name << " at " << hist1 << std::endl;
1241     if (hist1 != nullptr) {
1242       TH1D* hist = (TH1D*)(hist1->Clone());
1243       double ymin(0.90);
1244       sprintf(namep, "c_%s", name);
1245       TCanvas* pad = new TCanvas(namep, namep, 700, 500);
1246       pad->SetRightMargin(0.10);
1247       pad->SetTopMargin(0.10);
1248       hist->GetXaxis()->SetTitleSize(0.04);
1249       hist->GetXaxis()->SetTitle(xtitl.c_str());
1250       hist->GetYaxis()->SetTitle(ytitl.c_str());
1251       hist->GetYaxis()->SetLabelOffset(0.005);
1252       hist->GetYaxis()->SetTitleSize(0.04);
1253       hist->GetYaxis()->SetLabelSize(0.035);
1254       hist->GetYaxis()->SetTitleOffset(1.10);
1255       hist->GetXaxis()->SetRangeUser(0.25, 2.25);
1256       hist->SetMarkerStyle(20);
1257       hist->SetMarkerColor(2);
1258       hist->SetLineColor(2);
1259       hist->Draw();
1260       pad->Update();
1261       TPaveStats* st1 = (TPaveStats*)hist->GetListOfFunctions()->FindObject("stats");
1262       if (st1 != nullptr) {
1263         ymin = 0.70;
1264         st1->SetY1NDC(ymin);
1265         st1->SetY2NDC(0.90);
1266         st1->SetX1NDC(0.65);
1267         st1->SetX2NDC(0.90);
1268       }
1269       double ymx(0.96), xmi(0.25), xmx(0.90);
1270       char txt[100];
1271       if (lumi > 0.1) {
1272         ymx = ymin - 0.005;
1273         xmi = 0.45;
1274         TPaveText* txt0 = new TPaveText(0.65, 0.91, 0.90, 0.96, "blNDC");
1275         txt0->SetFillColor(0);
1276         sprintf(txt, "%4.1f TeV %5.1f fb^{-1}", ener, lumi);
1277         txt0->AddText(txt);
1278         txt0->Draw("same");
1279       }
1280       double ymi = ymx - 0.05;
1281       TPaveText* txt1 = new TPaveText(xmi, ymi, xmx, ymx, "blNDC");
1282       txt1->SetFillColor(0);
1283       if (text == "") {
1284         sprintf(txt, "Tracks with p = %s GeV at i#eta = %d", title[iene].c_str(), eta);
1285       } else {
1286         sprintf(txt, "Tracks with p = %s GeV at i#eta = %d (%s)", title[iene].c_str(), eta, text.c_str());
1287       }
1288       txt1->AddText(txt);
1289       txt1->Draw("same");
1290       double xmax = (dataMC) ? 0.33 : 0.44;
1291       ymi = (lumi > 0.1) ? 0.91 : 0.84;
1292       ymx = ymi + 0.05;
1293       TPaveText* txt2 = new TPaveText(0.11, ymi, xmax, ymx, "blNDC");
1294       txt2->SetFillColor(0);
1295       if (dataMC)
1296         sprintf(txt, "CMS Preliminary");
1297       else
1298         sprintf(txt, "CMS Simulation Preliminary");
1299       txt2->AddText(txt);
1300       txt2->Draw("same");
1301       pad->Modified();
1302       pad->Update();
1303       if (save > 0) {
1304         sprintf(name, "%s.pdf", pad->GetName());
1305         pad->Print(name);
1306       } else if (save < 0) {
1307         sprintf(name, "%s.C", pad->GetName());
1308         pad->Print(name);
1309       }
1310     }
1311   }
1312 }
1313 
1314 void PlotHists(std::string infile, std::string prefix, std::string text, bool drawStatBox = true, int save = 0) {
1315   int colors[6] = {1, 6, 4, 7, 2, 9};
1316   std::string types[6] = {"B", "C", "D", "E", "F", "G"};
1317   std::string names[3] = {"ratio20", "Z2", "W2"};
1318   std::string xtitl[3] = {"E_{HCAL}/(p-E_{ECAL})", "i#eta", "i#eta"};
1319   std::string ytitl[3] = {"Tracks", "MPV(E_{HCAL}/(p-E_{ECAL}))", "MPV/Width(E_{HCAL}/(p-E_{ECAL}))"};
1320 
1321   gStyle->SetCanvasBorderMode(0);
1322   gStyle->SetCanvasColor(kWhite);
1323   gStyle->SetPadColor(kWhite);
1324   gStyle->SetFillColor(kWhite);
1325   gStyle->SetOptTitle(0);
1326   if (drawStatBox)
1327     gStyle->SetOptFit(10);
1328   else
1329     gStyle->SetOptFit(0);
1330 
1331   char name[100], namep[100];
1332   TFile* file = new TFile(infile.c_str());
1333   for (int i = 0; i < 3; ++i) {
1334     std::vector<TH1D*> hists;
1335     std::vector<int> kks;
1336     double ymax(0.77);
1337     if (drawStatBox) {
1338       if (i == 0)
1339         gStyle->SetOptStat(1100);
1340       else
1341         gStyle->SetOptStat(10);
1342     } else {
1343       gStyle->SetOptStat(0);
1344       ymax = 0.89;
1345     }
1346     for (int k = 0; k < 6; ++k) {
1347       sprintf(name, "%s%s%s", prefix.c_str(), types[k].c_str(), names[i].c_str());
1348       TH1D* hist1 = (TH1D*)file->FindObjectAny(name);
1349       if (hist1 != nullptr) {
1350         hists.push_back((TH1D*)(hist1->Clone()));
1351         kks.push_back(k);
1352       }
1353     }
1354     if (hists.size() > 0) {
1355       sprintf(namep, "c_%s%s", prefix.c_str(), names[i].c_str());
1356       TCanvas* pad = new TCanvas(namep, namep, 700, 500);
1357       TLegend* legend = new TLegend(0.44, 0.89 - 0.055 * hists.size(), 0.69, ymax);
1358       legend->SetFillColor(kWhite);
1359       pad->SetRightMargin(0.10);
1360       pad->SetTopMargin(0.10);
1361       double ymax(0.90);
1362       double dy = (i == 0) ? 0.13 : 0.08;
1363       for (unsigned int jk = 0; jk < hists.size(); ++jk) {
1364         int k = kks[jk];
1365         hists[jk]->GetXaxis()->SetTitle(xtitl[i].c_str());
1366         hists[jk]->GetYaxis()->SetTitle(ytitl[i].c_str());
1367         hists[jk]->GetYaxis()->SetLabelOffset(0.005);
1368         hists[jk]->GetYaxis()->SetLabelSize(0.035);
1369         hists[jk]->GetYaxis()->SetTitleOffset(1.15);
1370         if (i == 0) {
1371           hists[jk]->GetXaxis()->SetRangeUser(0.0, 2.5);
1372         } else if (i == 1) {
1373           hists[jk]->GetYaxis()->SetRangeUser(0.5, 2.0);
1374         } else {
1375           hists[jk]->GetYaxis()->SetRangeUser(0.0, 0.5);
1376         }
1377         hists[jk]->SetMarkerStyle(20);
1378         hists[jk]->SetMarkerColor(colors[k]);
1379         hists[jk]->SetLineColor(colors[k]);
1380         if (jk == 0)
1381           hists[jk]->Draw();
1382         else
1383           hists[jk]->Draw("sames");
1384         pad->Update();
1385         TPaveStats* st1 = (TPaveStats*)hists[jk]->GetListOfFunctions()->FindObject("stats");
1386         if (st1 != nullptr) {
1387           double ymin = ymax - dy;
1388           st1->SetLineColor(colors[k]);
1389           st1->SetTextColor(colors[k]);
1390           st1->SetY1NDC(ymin);
1391           st1->SetY2NDC(ymax);
1392           st1->SetX1NDC(0.70);
1393           st1->SetX2NDC(0.90);
1394           ymax = ymin;
1395         }
1396         sprintf(name, "%s%s", text.c_str(), types[k].c_str());
1397         legend->AddEntry(hists[jk], name, "lp");
1398       }
1399       legend->Draw("same");
1400       pad->Update();
1401       TPaveText* txt1 = new TPaveText(0.34, 0.825, 0.69, 0.895, "blNDC");
1402       txt1->SetFillColor(0);
1403       char txt[100];
1404       sprintf(txt, "Tracks with p = 40:60 GeV");
1405       txt1->AddText(txt);
1406       txt1->Draw("same");
1407       TPaveText* txt2 = new TPaveText(0.11, 0.825, 0.33, 0.895, "blNDC");
1408       txt2->SetFillColor(0);
1409       sprintf(txt, "CMS Preliminary");
1410       txt2->AddText(txt);
1411       txt2->Draw("same");
1412       if (!drawStatBox && i == 1) {
1413         double xmin = hists[0]->GetBinLowEdge(1);
1414         int nbin = hists[0]->GetNbinsX();
1415         double xmax = hists[0]->GetBinLowEdge(nbin) + hists[0]->GetBinWidth(nbin);
1416         TLine line = TLine(xmin, 1.0, xmax, 1.0);  //etamin,1.0,etamax,1.0);
1417         line.SetLineWidth(4);
1418         line.Draw("same");
1419         pad->Update();
1420       }
1421       pad->Modified();
1422       pad->Update();
1423       if (save > 0) {
1424         sprintf(name, "%s.pdf", pad->GetName());
1425         pad->Print(name);
1426       } else if (save < 0) {
1427         sprintf(name, "%s.C", pad->GetName());
1428         pad->Print(name);
1429       }
1430     }
1431   }
1432 }
1433 
1434 void PlotTwoHists(std::string infile,
1435                   std::string prefix1,
1436                   std::string text1,
1437                   std::string prefix2,
1438                   std::string text2,
1439                   std::string text0,
1440                   int type = 0,
1441                   int iname = 3,
1442                   double lumi = 0,
1443                   double ener = 13.0,
1444                   int drawStatBox = 0,
1445                   int save = 0) {
1446   int colors[2] = {2, 4};
1447   int numb[2] = {5, 1};
1448   std::string names0[5] = {"ratio00", "ratio00One", "etaB04One", "Z0", "W0"};
1449   std::string names1[5] = {"ratio10", "ratio10One", "etaB14One", "Z1", "W1"};
1450   std::string names2[5] = {"ratio30", "ratio30One", "etaB34One", "Z3", "W3"};
1451   std::string xtitl1[5] = {"E_{HCAL}/(p-E_{ECAL})", "E_{HCAL}/(p-E_{ECAL})", "E_{HCAL}/(p-E_{ECAL})", "i#eta", "i#eta"};
1452   std::string ytitl1[5] = {
1453       "Tracks", "Tracks", "Tracks", "MPV(E_{HCAL}/(p-E_{ECAL}))", "MPV/Width(E_{HCAL}/(p-E_{ECAL}))"};
1454   std::string names3[1] = {"R"};
1455   std::string xtitl2[1] = {"RBX #"};
1456   std::string ytitl2[1] = {"MPV(E_{HCAL}/(p-E_{ECAL}))"};
1457 
1458   gStyle->SetCanvasBorderMode(0);
1459   gStyle->SetCanvasColor(kWhite);
1460   gStyle->SetPadColor(kWhite);
1461   gStyle->SetFillColor(kWhite);
1462   gStyle->SetOptTitle(0);
1463   if ((drawStatBox / 10) % 10 > 0)
1464     gStyle->SetOptFit(10);
1465   else
1466     gStyle->SetOptFit(0);
1467 
1468   if (type != 1)
1469     type = 0;
1470   char name[100], namep[100];
1471   TFile* file = new TFile(infile.c_str());
1472   for (int i = 0; i < numb[type]; ++i) {
1473     std::vector<TH1D*> hists;
1474     std::vector<int> kks;
1475     double ymax(0.77);
1476     if (drawStatBox % 10 > 0) {
1477       if (i != 2)
1478         gStyle->SetOptStat(1100);
1479       else
1480         gStyle->SetOptStat(10);
1481     } else {
1482       gStyle->SetOptStat(0);
1483       ymax = 0.82;
1484     }
1485     for (int k = 0; k < 2; ++k) {
1486       std::string prefix = (k == 0) ? prefix1 : prefix2;
1487       if (type == 0) {
1488         if (iname == 0)
1489           sprintf(name, "%s%s", prefix.c_str(), names0[i].c_str());
1490         else if (iname == 1)
1491           sprintf(name, "%s%s", prefix.c_str(), names1[i].c_str());
1492         else
1493           sprintf(name, "%s%s", prefix.c_str(), names2[i].c_str());
1494       } else {
1495         sprintf(name, "%s%s%d", prefix.c_str(), names3[i].c_str(), iname);
1496       }
1497       TH1D* hist1 = (TH1D*)file->FindObjectAny(name);
1498       if (hist1 != nullptr) {
1499         hists.push_back((TH1D*)(hist1->Clone()));
1500         kks.push_back(k);
1501       }
1502     }
1503     if (hists.size() == 2) {
1504       if (type == 0) {
1505         if (iname == 0)
1506           sprintf(namep, "c_%s%s%s", prefix1.c_str(), prefix2.c_str(), names0[i].c_str());
1507         else if (iname == 1)
1508           sprintf(namep, "c_%s%s%s", prefix1.c_str(), prefix2.c_str(), names1[i].c_str());
1509         else
1510           sprintf(namep, "c_%s%s%s", prefix1.c_str(), prefix2.c_str(), names2[i].c_str());
1511       } else {
1512         sprintf(namep, "c_%s%s%s%d", prefix1.c_str(), prefix2.c_str(), names3[i].c_str(), iname);
1513       }
1514       double ymax(0.90);
1515       double dy = (i == 0) ? 0.13 : 0.08;
1516       double ymx0 = (drawStatBox == 0) ? (ymax - .01) : (ymax - dy * hists.size() - .01);
1517       TCanvas* pad = new TCanvas(namep, namep, 700, 500);
1518       TLegend* legend = new TLegend(0.64, ymx0 - 0.05 * hists.size(), 0.89, ymx0);
1519       legend->SetFillColor(kWhite);
1520       pad->SetRightMargin(0.10);
1521       pad->SetTopMargin(0.10);
1522       for (unsigned int jk = 0; jk < hists.size(); ++jk) {
1523         int k = kks[jk];
1524         hists[jk]->GetXaxis()->SetTitleSize(0.040);
1525         if (type == 0) {
1526           hists[jk]->GetXaxis()->SetTitle(xtitl1[i].c_str());
1527           hists[jk]->GetYaxis()->SetTitle(ytitl1[i].c_str());
1528         } else {
1529           hists[jk]->GetXaxis()->SetTitle(xtitl2[i].c_str());
1530           hists[jk]->GetYaxis()->SetTitle(ytitl2[i].c_str());
1531         }
1532         hists[jk]->GetYaxis()->SetLabelOffset(0.005);
1533         hists[jk]->GetYaxis()->SetLabelSize(0.035);
1534         hists[jk]->GetYaxis()->SetTitleSize(0.040);
1535         hists[jk]->GetYaxis()->SetTitleOffset(1.15);
1536         if ((type == 0) && (i != 3) && (i != 4))
1537           hists[jk]->GetXaxis()->SetRangeUser(0.0, 5.0);
1538         if (type == 0) {
1539           if (i == 3)
1540             hists[jk]->GetYaxis()->SetRangeUser(0.8, 1.2);
1541           else if (i == 4)
1542             hists[jk]->GetYaxis()->SetRangeUser(0.0, 0.5);
1543         }
1544         if (type != 0)
1545           hists[jk]->GetYaxis()->SetRangeUser(0.75, 1.2);
1546         hists[jk]->SetMarkerStyle(20);
1547         hists[jk]->SetMarkerColor(colors[k]);
1548         hists[jk]->SetLineColor(colors[k]);
1549         if (jk == 0)
1550           hists[jk]->Draw();
1551         else
1552           hists[jk]->Draw("sames");
1553         pad->Update();
1554         TPaveStats* st1 = (TPaveStats*)hists[jk]->GetListOfFunctions()->FindObject("stats");
1555         if (st1 != nullptr) {
1556           double ymin = ymax - dy;
1557           st1->SetLineColor(colors[k]);
1558           st1->SetTextColor(colors[k]);
1559           st1->SetY1NDC(ymin);
1560           st1->SetY2NDC(ymax);
1561           st1->SetX1NDC(0.70);
1562           st1->SetX2NDC(0.90);
1563           ymax = ymin;
1564         }
1565         if (k == 0)
1566           sprintf(name, "%s", text1.c_str());
1567         else
1568           sprintf(name, "%s", text2.c_str());
1569         legend->AddEntry(hists[jk], name, "lp");
1570       }
1571       legend->Draw("same");
1572       pad->Update();
1573       char txt[100];
1574       double xmi(0.10), xmx(0.895), ymx(0.95);
1575       if (lumi > 0.01) {
1576         xmx = 0.70;
1577         xmi = 0.30;
1578         TPaveText* txt0 = new TPaveText(0.705, 0.905, 0.90, 0.95, "blNDC");
1579         txt0->SetFillColor(0);
1580         sprintf(txt, "%4.1f TeV %5.1f fb^{-1}", ener, lumi);
1581         txt0->AddText(txt);
1582         txt0->Draw("same");
1583       }
1584       double ymi = ymx - 0.045;
1585       TPaveText* txt1 = new TPaveText(xmi, ymi, xmx, ymx, "blNDC");
1586       txt1->SetFillColor(0);
1587       if (iname == 0)
1588         sprintf(txt, "p = 20:30 GeV %s", text0.c_str());
1589       else
1590         sprintf(txt, "p = 40:60 GeV %s", text0.c_str());
1591       txt1->AddText(txt);
1592       txt1->Draw("same");
1593       ymi = (lumi > 0.1) ? 0.905 : 0.85;
1594       ymx = ymi + 0.045;
1595       TPaveText* txt2 = new TPaveText(0.105, ymi, 0.295, ymx, "blNDC");
1596       txt2->SetFillColor(0);
1597       sprintf(txt, "CMS Preliminary");
1598       txt2->AddText(txt);
1599       txt2->Draw("same");
1600       pad->Modified();
1601       if ((drawStatBox == 0) && (i == 3)) {
1602         double xmin = hists[0]->GetBinLowEdge(1);
1603         int nbin = hists[0]->GetNbinsX();
1604         double xmax = hists[0]->GetBinLowEdge(nbin) + hists[0]->GetBinWidth(nbin);
1605         TLine* line = new TLine(xmin, 1.0, xmax, 1.0);  //etamin,1.0,etamax,1.0);
1606         line->SetLineWidth(2);
1607         line->Draw("same");
1608         pad->Update();
1609       }
1610       pad->Update();
1611       if (save > 0) {
1612         sprintf(name, "%s.pdf", pad->GetName());
1613         pad->Print(name);
1614       } else if (save < 0) {
1615         sprintf(name, "%s.C", pad->GetName());
1616         pad->Print(name);
1617       }
1618     }
1619   }
1620 }
1621 
1622 void PlotFiveHists(std::string infile,
1623                    std::string text0,
1624                    std::string prefix0,
1625                    int type = 0,
1626                    int iname = 3,
1627                    int drawStatBox = 0,
1628                    bool normalize = false,
1629                    int save = 0,
1630                    std::string prefix1 = "",
1631                    std::string text1 = "",
1632                    std::string prefix2 = "",
1633                    std::string text2 = "",
1634                    std::string prefix3 = "",
1635                    std::string text3 = "",
1636                    std::string prefix4 = "",
1637                    std::string text4 = "",
1638                    std::string prefix5 = "",
1639                    std::string text5 = "") {
1640   int colors[5] = {2, 4, 6, 1, 7};
1641   int numb[3] = {5, 1, 4};
1642   std::string names0[5] = {"ratio00", "ratio00One", "etaB04", "Z0", "W0"};
1643   std::string names1[5] = {"ratio10", "ratio10One", "etaB14", "Z1", "W1"};
1644   std::string names2[5] = {"ratio30", "ratio30One", "etaB34", "Z3", "W3"};
1645   std::string xtitl1[5] = {"E_{HCAL}/(p-E_{ECAL})", "E_{HCAL}/(p-E_{ECAL})", "E_{HCAL}/(p-E_{ECAL})", "i#eta", "i#eta"};
1646   std::string ytitl1[5] = {
1647       "Tracks", "Tracks", "Tracks", "MPV(E_{HCAL}/(p-E_{ECAL}))", "MPV/Width(E_{HCAL}/(p-E_{ECAL}))"};
1648   std::string names3[1] = {"R"};
1649   std::string xtitl2[1] = {"RBX #"};
1650   std::string ytitl2[1] = {"MPV(E_{HCAL}/(p-E_{ECAL}))"};
1651   std::string names4[4] = {"pp21", "pp22", "pp23", "pp24"};
1652   std::string xtitl3[4] = {"p (GeV)", "p (GeV)", "p (GeV)", "p (GeV)"};
1653   std::string ytitl3[4] = {"Tracks", "Tracks", "Tracks", "Tracks"};
1654   std::string title3[4] = {"Barrel", "Transition", "Endcap", "Combined"};
1655 
1656   gStyle->SetCanvasBorderMode(0);
1657   gStyle->SetCanvasColor(kWhite);
1658   gStyle->SetPadColor(kWhite);
1659   gStyle->SetFillColor(kWhite);
1660   gStyle->SetOptTitle(0);
1661   if ((drawStatBox / 10) % 10 > 0)
1662     gStyle->SetOptFit(10);
1663   else
1664     gStyle->SetOptFit(0);
1665 
1666   if (type != 1 && type != 2)
1667     type = 0;
1668   char name[100], namep[100];
1669   TFile* file = new TFile(infile.c_str());
1670   for (int i = 0; i < numb[type]; ++i) {
1671     std::vector<TH1D*> hists;
1672     std::vector<int> kks;
1673     std::vector<std::string> texts;
1674     double ymax(0.77);
1675     if (drawStatBox % 10 > 0) {
1676       if (type == 2)
1677         gStyle->SetOptStat(1110);
1678       else if (i != 3)
1679         gStyle->SetOptStat(1100);
1680       else
1681         gStyle->SetOptStat(10);
1682     } else {
1683       gStyle->SetOptStat(0);
1684       ymax = 0.82;
1685     }
1686     for (int k = 0; k < 5; ++k) {
1687       std::string prefix, text;
1688       if (k == 0) {
1689         prefix = prefix1;
1690         text = text1;
1691       } else if (k == 1) {
1692         prefix = prefix2;
1693         text = text2;
1694       } else if (k == 2) {
1695         prefix = prefix3;
1696         text = text3;
1697       } else if (k == 3) {
1698         prefix = prefix4;
1699         text = text4;
1700       } else {
1701         prefix = prefix5;
1702         text = text5;
1703       }
1704       if (prefix != "") {
1705         if (type == 0) {
1706           if (iname == 0)
1707             sprintf(name, "%s%s", prefix.c_str(), names0[i].c_str());
1708           else if (iname == 1)
1709             sprintf(name, "%s%s", prefix.c_str(), names1[i].c_str());
1710           else
1711             sprintf(name, "%s%s", prefix.c_str(), names2[i].c_str());
1712         } else if (type == 1) {
1713           sprintf(name, "%s%s%d", prefix.c_str(), names3[i].c_str(), iname);
1714         } else {
1715           sprintf(name, "%s%s", prefix.c_str(), names4[i].c_str());
1716         }
1717         TH1D* hist1 = (TH1D*)file->FindObjectAny(name);
1718         if (hist1 != nullptr) {
1719           hists.push_back((TH1D*)(hist1->Clone()));
1720           kks.push_back(k);
1721           texts.push_back(text);
1722         }
1723       }
1724     }
1725     if (hists.size() > 0) {
1726       if (type == 0) {
1727         if (iname == 0)
1728           sprintf(namep, "c_%s%s", prefix0.c_str(), names0[i].c_str());
1729         else if (iname == 1)
1730           sprintf(namep, "c_%s%s", prefix0.c_str(), names1[i].c_str());
1731         else
1732           sprintf(namep, "c_%s%s", prefix0.c_str(), names2[i].c_str());
1733       } else if (type == 1) {
1734         sprintf(namep, "c_%s%s%d", prefix0.c_str(), names3[i].c_str(), iname);
1735       } else {
1736         sprintf(namep, "c_%s%s", prefix0.c_str(), names4[i].c_str());
1737       }
1738       double ymax(0.90);
1739       double dy = (i == 0 && type == 0) ? 0.13 : 0.08;
1740       double ymx0 = (drawStatBox == 0) ? (ymax - .01) : (ymax - dy * hists.size() - .01);
1741       TCanvas* pad = new TCanvas(namep, namep, 700, 500);
1742       TLegend* legend = new TLegend(0.64, ymx0 - 0.05 * hists.size(), 0.89, ymx0);
1743       legend->SetFillColor(kWhite);
1744       pad->SetRightMargin(0.10);
1745       pad->SetTopMargin(0.10);
1746       for (unsigned int jk = 0; jk < hists.size(); ++jk) {
1747         int k = kks[jk];
1748         hists[jk]->GetXaxis()->SetTitleSize(0.040);
1749         if (type == 0) {
1750           hists[jk]->GetXaxis()->SetTitle(xtitl1[i].c_str());
1751           hists[jk]->GetYaxis()->SetTitle(ytitl1[i].c_str());
1752         } else if (type == 1) {
1753           hists[jk]->GetXaxis()->SetTitle(xtitl2[i].c_str());
1754           hists[jk]->GetYaxis()->SetTitle(ytitl2[i].c_str());
1755         } else {
1756           hists[jk]->GetXaxis()->SetTitle(xtitl3[i].c_str());
1757           hists[jk]->GetYaxis()->SetTitle(ytitl3[i].c_str());
1758         }
1759         hists[jk]->GetYaxis()->SetLabelOffset(0.005);
1760         hists[jk]->GetYaxis()->SetLabelSize(0.035);
1761         hists[jk]->GetYaxis()->SetTitleSize(0.040);
1762         hists[jk]->GetYaxis()->SetTitleOffset(1.15);
1763         if ((type == 0) && (i != 3) && (i != 4))
1764           hists[jk]->GetXaxis()->SetRangeUser(0.0, 2.5);
1765         if (type == 0) {
1766           if (i == 3)
1767             hists[jk]->GetYaxis()->SetRangeUser(0.8, 1.2);
1768           else if (i == 4)
1769             hists[jk]->GetYaxis()->SetRangeUser(0.0, 0.5);
1770         }
1771         if (type == 1)
1772           hists[jk]->GetYaxis()->SetRangeUser(0.75, 1.2);
1773         hists[jk]->SetMarkerStyle(20);
1774         hists[jk]->SetMarkerColor(colors[k]);
1775         hists[jk]->SetLineColor(colors[k]);
1776         if (normalize && ((type == 2) || ((type == 0) && (i != 3)))) {
1777           if (jk == 0)
1778             hists[jk]->DrawNormalized();
1779           else
1780             hists[jk]->DrawNormalized("sames");
1781         } else {
1782           if (jk == 0)
1783             hists[jk]->Draw();
1784           else
1785             hists[jk]->Draw("sames");
1786         }
1787         pad->Update();
1788         TPaveStats* st1 = (TPaveStats*)hists[jk]->GetListOfFunctions()->FindObject("stats");
1789         if (st1 != nullptr) {
1790           double ymin = ymax - dy;
1791           st1->SetLineColor(colors[k]);
1792           st1->SetTextColor(colors[k]);
1793           st1->SetY1NDC(ymin);
1794           st1->SetY2NDC(ymax);
1795           st1->SetX1NDC(0.70);
1796           st1->SetX2NDC(0.90);
1797           ymax = ymin;
1798         }
1799         sprintf(name, "%s", texts[jk].c_str());
1800         legend->AddEntry(hists[jk], name, "lp");
1801       }
1802       legend->Draw("same");
1803       pad->Update();
1804       TPaveText* txt1 = new TPaveText(0.10, 0.905, 0.80, 0.95, "blNDC");
1805       txt1->SetFillColor(0);
1806       char txt[100];
1807       if (type == 2) {
1808         sprintf(txt, "p = 40:60 GeV (%s)", title3[i].c_str());
1809       } else if (((type == 0) && (iname == 0))) {
1810         sprintf(txt, "p = 20:30 GeV %s", text0.c_str());
1811       } else {
1812         sprintf(txt, "p = 40:60 GeV %s", text0.c_str());
1813       }
1814       txt1->AddText(txt);
1815       txt1->Draw("same");
1816       TPaveText* txt2 = new TPaveText(0.11, 0.825, 0.33, 0.895, "blNDC");
1817       txt2->SetFillColor(0);
1818       sprintf(txt, "CMS Preliminary");
1819       txt2->AddText(txt);
1820       txt2->Draw("same");
1821       pad->Modified();
1822       if ((drawStatBox == 0) && (i == 3)) {
1823         double xmin = hists[0]->GetBinLowEdge(1);
1824         int nbin = hists[0]->GetNbinsX();
1825         double xmax = hists[0]->GetBinLowEdge(nbin) + hists[0]->GetBinWidth(nbin);
1826         TLine* line = new TLine(xmin, 1.0, xmax, 1.0);  //etamin,1.0,etamax,1.0);
1827         line->SetLineWidth(2);
1828         line->Draw("same");
1829         pad->Update();
1830       }
1831       pad->Update();
1832       if (save > 0) {
1833         sprintf(name, "%s.pdf", pad->GetName());
1834         pad->Print(name);
1835       } else if (save < 0) {
1836         sprintf(name, "%s.C", pad->GetName());
1837         pad->Print(name);
1838       }
1839     }
1840   }
1841 }
1842 
1843 void PlotHistCorrResults(std::string infile, std::string text, std::string prefixF, int save = 0) {
1844   std::string name[5] = {"Eta1Bf", "Eta2Bf", "Eta1Af", "Eta2Af", "Cvg0"};
1845   std::string title[5] = {"Mean at the start of itertions",
1846                           "Median at the start of itertions",
1847                           "Mean at the end of itertions",
1848                           "Median at the end of itertions",
1849                           ""};
1850   int type[5] = {0, 0, 0, 0, 1};
1851 
1852   gStyle->SetCanvasBorderMode(0);
1853   gStyle->SetCanvasColor(kWhite);
1854   gStyle->SetPadColor(kWhite);
1855   gStyle->SetFillColor(kWhite);
1856   gStyle->SetOptTitle(0);
1857   gStyle->SetOptStat(10);
1858   gStyle->SetOptFit(10);
1859   TFile* file = new TFile(infile.c_str());
1860   char namep[100];
1861   for (int k = 0; k < 5; ++k) {
1862     TH1D* hist1 = (TH1D*)file->FindObjectAny(name[k].c_str());
1863     if (hist1 != nullptr) {
1864       TH1D* hist = (TH1D*)(hist1->Clone());
1865       sprintf(namep, "c_%s%s", prefixF.c_str(), name[k].c_str());
1866       TCanvas* pad = new TCanvas(namep, namep, 700, 500);
1867       pad->SetRightMargin(0.10);
1868       pad->SetTopMargin(0.10);
1869       hist->GetYaxis()->SetLabelOffset(0.005);
1870       hist->GetYaxis()->SetTitleOffset(1.20);
1871       double xmin = hist->GetBinLowEdge(1);
1872       int nbin = hist->GetNbinsX();
1873       double xmax = hist->GetBinLowEdge(nbin) + hist->GetBinWidth(nbin);
1874       std::cout << hist->GetTitle() << " Bins " << nbin << ":" << xmin << ":" << xmax << std::endl;
1875       double xlow(0.12), ylow(0.82);
1876       char txt[100], option[2];
1877       if (type[k] == 0) {
1878         sprintf(namep, "f_%s", name[k].c_str());
1879         TF1* func = new TF1(namep, "pol0", xmin, xmax);
1880         hist->Fit(func, "+QWL", "");
1881         if (text == "")
1882           sprintf(txt, "%s", title[k].c_str());
1883         else
1884           sprintf(txt, "%s (balancing the %s)", title[k].c_str(), text.c_str());
1885         sprintf(option, "%s", "");
1886       } else {
1887         hist->GetXaxis()->SetTitle("Iterations");
1888         hist->GetYaxis()->SetTitle("Convergence");
1889         hist->SetMarkerStyle(20);
1890         hist->SetMarkerColor(2);
1891         hist->SetMarkerSize(0.8);
1892         xlow = 0.50;
1893         ylow = 0.86;
1894         sprintf(txt, "(%s)", text.c_str());
1895         sprintf(option, "%s", "p");
1896       }
1897       hist->Draw(option);
1898       pad->Update();
1899       TPaveStats* st1 = (TPaveStats*)hist->GetListOfFunctions()->FindObject("stats");
1900       if (st1 != nullptr) {
1901         st1->SetY1NDC(ylow);
1902         st1->SetY2NDC(0.90);
1903         st1->SetX1NDC(0.70);
1904         st1->SetX2NDC(0.90);
1905       }
1906       TPaveText* txt1 = new TPaveText(xlow, 0.82, 0.68, 0.88, "blNDC");
1907       txt1->SetFillColor(0);
1908       txt1->AddText(txt);
1909       txt1->Draw("same");
1910       pad->Modified();
1911       pad->Update();
1912       if (save > 0) {
1913         sprintf(namep, "%s.pdf", pad->GetName());
1914         pad->Print(namep);
1915       } else if (save < 0) {
1916         sprintf(namep, "%s.C", pad->GetName());
1917         pad->Print(namep);
1918       }
1919     }
1920   }
1921 }
1922 
1923 void PlotHistCorrFactor(char* infile,
1924                         std::string text,
1925                         std::string prefixF = "",
1926                         double scale = 1.0,
1927                         int nmin = 100,
1928                         bool dataMC = false,
1929                         bool drawStatBox = true,
1930                         int iformat = 0,
1931                         int save = 0) {
1932   std::map<int, cfactors> cfacs;
1933   int etamin(100), etamax(-100), maxdepth(0);
1934   readCorrFactors(infile, scale, cfacs, etamin, etamax, maxdepth, iformat);
1935 
1936   gStyle->SetCanvasBorderMode(0);
1937   gStyle->SetCanvasColor(kWhite);
1938   gStyle->SetPadColor(kWhite);
1939   gStyle->SetFillColor(kWhite);
1940   gStyle->SetOptTitle(0);
1941   if (drawStatBox) {
1942     gStyle->SetOptStat(10);
1943     gStyle->SetOptFit(10);
1944   } else {
1945     gStyle->SetOptStat(0);
1946     gStyle->SetOptFit(0);
1947   }
1948   int colors[7] = {1, 6, 4, 7, 2, 9, 3};
1949   int mtype[7] = {20, 21, 22, 23, 24, 33, 25};
1950   int nbin = etamax - etamin + 1;
1951   std::vector<TH1D*> hists;
1952   std::vector<int> entries;
1953   char name[100];
1954   double dy(0);
1955   int fits(0);
1956   for (int j = 0; j < maxdepth; ++j) {
1957     sprintf(name, "hd%d", j + 1);
1958     TH1D* h = new TH1D(name, name, nbin, etamin, etamax);
1959     int nent(0);
1960     for (std::map<int, cfactors>::const_iterator itr = cfacs.begin(); itr != cfacs.end(); ++itr) {
1961       if ((itr->second).depth == j + 1) {
1962         int ieta = (itr->second).ieta;
1963         int bin = ieta - etamin + 1;
1964         float val = (itr->second).corrf;
1965         float dvl = (itr->second).dcorr;
1966         h->SetBinContent(bin, val);
1967         h->SetBinError(bin, dvl);
1968         nent++;
1969       }
1970     }
1971     if (nent > nmin) {
1972       fits++;
1973       dy += 0.025;
1974       sprintf(name, "hdf%d", j + 1);
1975       TF1* func = new TF1(name, "pol0", etamin, etamax);
1976       h->Fit(func, "+QWLR", "");
1977     }
1978     h->SetLineColor(colors[j]);
1979     h->SetMarkerColor(colors[j]);
1980     h->SetMarkerStyle(mtype[j]);
1981     h->GetXaxis()->SetTitle("i#eta");
1982     h->GetYaxis()->SetTitle("Correction Factor");
1983     h->GetYaxis()->SetLabelOffset(0.005);
1984     h->GetYaxis()->SetTitleOffset(1.20);
1985     h->GetYaxis()->SetRangeUser(0.0, 2.0);
1986     hists.push_back(h);
1987     entries.push_back(nent);
1988     dy += 0.025;
1989   }
1990   sprintf(name, "c_%sCorrFactor", prefixF.c_str());
1991   TCanvas* pad = new TCanvas(name, name, 700, 500);
1992   pad->SetRightMargin(0.10);
1993   pad->SetTopMargin(0.10);
1994   double yh = 0.90;
1995   // double yl = yh - 0.025 * hists.size() - dy - 0.01;
1996   double yl = 0.15;
1997   TLegend* legend = new TLegend(0.35, yl, 0.65, yl + 0.04 * hists.size());
1998   legend->SetFillColor(kWhite);
1999   for (unsigned int k = 0; k < hists.size(); ++k) {
2000     if (k == 0)
2001       hists[k]->Draw("");
2002     else
2003       hists[k]->Draw("sames");
2004     pad->Update();
2005     TPaveStats* st1 = (TPaveStats*)hists[k]->GetListOfFunctions()->FindObject("stats");
2006     if (st1 != nullptr) {
2007       dy = (entries[k] > nmin) ? 0.05 : 0.025;
2008       st1->SetLineColor(colors[k]);
2009       st1->SetTextColor(colors[k]);
2010       st1->SetY1NDC(yh - dy);
2011       st1->SetY2NDC(yh);
2012       st1->SetX1NDC(0.70);
2013       st1->SetX2NDC(0.90);
2014       yh -= dy;
2015     }
2016     sprintf(name, "Depth %d (%s)", k + 1, text.c_str());
2017     legend->AddEntry(hists[k], name, "lp");
2018   }
2019   legend->Draw("same");
2020   pad->Update();
2021   if (fits < 1) {
2022     double xmin = hists[0]->GetBinLowEdge(1);
2023     int nbin = hists[0]->GetNbinsX();
2024     double xmax = hists[0]->GetBinLowEdge(nbin) + hists[0]->GetBinWidth(nbin);
2025     TLine* line = new TLine(xmin, 1.0, xmax, 1.0);
2026     line->SetLineColor(9);
2027     line->SetLineWidth(2);
2028     line->SetLineStyle(2);
2029     line->Draw("same");
2030     pad->Modified();
2031     pad->Update();
2032   }
2033   char txt1[30];
2034   double xmax = (dataMC) ? 0.33 : 0.44;
2035   TPaveText* txt2 = new TPaveText(0.11, 0.85, xmax, 0.89, "blNDC");
2036   txt2->SetFillColor(0);
2037   if (dataMC)
2038     sprintf(txt1, "CMS Preliminary");
2039   else
2040     sprintf(txt1, "CMS Simulation Preliminary");
2041   txt2->AddText(txt1);
2042   txt2->Draw("same");
2043   pad->Modified();
2044   pad->Update();
2045   if (save > 0) {
2046     sprintf(name, "%s.pdf", pad->GetName());
2047     pad->Print(name);
2048   } else if (save < 0) {
2049     sprintf(name, "%s.C", pad->GetName());
2050     pad->Print(name);
2051   }
2052 }
2053 
2054 void PlotHistCorrAsymmetry(char* infile, std::string text, std::string prefixF = "", int iformat = 0, int save = 0) {
2055   std::map<int, cfactors> cfacs;
2056   int etamin(100), etamax(-100), maxdepth(0);
2057   double scale(1.0);
2058   readCorrFactors(infile, scale, cfacs, etamin, etamax, maxdepth, iformat);
2059 
2060   gStyle->SetCanvasBorderMode(0);
2061   gStyle->SetCanvasColor(kWhite);
2062   gStyle->SetPadColor(kWhite);
2063   gStyle->SetFillColor(kWhite);
2064   gStyle->SetOptTitle(0);
2065   gStyle->SetOptStat(0);
2066   gStyle->SetOptFit(10);
2067   int colors[6] = {1, 6, 4, 7, 2, 9};
2068   int mtype[6] = {20, 21, 22, 23, 24, 33};
2069   int nbin = etamax + 1;
2070   std::vector<TH1D*> hists;
2071   std::vector<int> entries;
2072   char name[100];
2073   double dy(0);
2074   for (int j = 0; j < maxdepth; ++j) {
2075     sprintf(name, "hd%d", j + 1);
2076     TH1D* h = new TH1D(name, name, nbin, 0, etamax);
2077     int nent(0);
2078     for (std::map<int, cfactors>::const_iterator itr = cfacs.begin(); itr != cfacs.end(); ++itr) {
2079       if ((itr->second).depth == j + 1) {
2080         int ieta = (itr->second).ieta;
2081         float vl1 = (itr->second).corrf;
2082         float dv1 = (itr->second).dcorr;
2083         if (ieta > 0) {
2084           for (std::map<int, cfactors>::const_iterator ktr = cfacs.begin(); ktr != cfacs.end(); ++ktr) {
2085             if (((ktr->second).depth == j + 1) && ((ktr->second).ieta == -ieta)) {
2086               float vl2 = (ktr->second).corrf;
2087               float dv2 = (ktr->second).dcorr;
2088               float val = 2.0 * (vl1 - vl2) / (vl1 + vl2);
2089               float dvl = (4.0 * sqrt(vl1 * vl1 * dv2 * dv2 + vl2 * vl2 * dv1 * dv1) / ((vl1 + vl2) * (vl1 + vl2)));
2090               int bin = ieta;
2091               h->SetBinContent(bin, val);
2092               h->SetBinError(bin, dvl);
2093               nent++;
2094             }
2095           }
2096         }
2097       }
2098     }
2099     h->SetLineColor(colors[j]);
2100     h->SetMarkerColor(colors[j]);
2101     h->SetMarkerStyle(mtype[j]);
2102     h->GetXaxis()->SetTitle("i#eta");
2103     h->GetYaxis()->SetTitle("Asymmetry in Correction Factor");
2104     h->GetYaxis()->SetLabelOffset(0.005);
2105     h->GetYaxis()->SetTitleOffset(1.20);
2106     h->GetYaxis()->SetRangeUser(-0.25, 0.25);
2107     hists.push_back(h);
2108     entries.push_back(nent);
2109     dy += 0.025;
2110   }
2111   sprintf(name, "c_%sCorrAsymmetry", prefixF.c_str());
2112   TCanvas* pad = new TCanvas(name, name, 700, 500);
2113   pad->SetRightMargin(0.10);
2114   pad->SetTopMargin(0.10);
2115   double yh = 0.90;
2116   double yl = yh - 0.035 * hists.size() - dy - 0.01;
2117   TLegend* legend = new TLegend(0.60, yl, 0.90, yl + 0.035 * hists.size());
2118   legend->SetFillColor(kWhite);
2119   for (unsigned int k = 0; k < hists.size(); ++k) {
2120     if (k == 0)
2121       hists[k]->Draw("");
2122     else
2123       hists[k]->Draw("sames");
2124     pad->Update();
2125     sprintf(name, "Depth %d (%s)", k + 1, text.c_str());
2126     legend->AddEntry(hists[k], name, "lp");
2127   }
2128   legend->Draw("same");
2129   pad->Update();
2130 
2131   TLine* line = new TLine(0.0, 0.0, etamax, 0.0);
2132   line->SetLineColor(9);
2133   line->SetLineWidth(2);
2134   line->SetLineStyle(2);
2135   line->Draw("same");
2136   pad->Update();
2137 
2138   if (save > 0) {
2139     sprintf(name, "%s.pdf", pad->GetName());
2140     pad->Print(name);
2141   } else if (save < 0) {
2142     sprintf(name, "%s.C", pad->GetName());
2143     pad->Print(name);
2144   }
2145 }
2146 
2147 void PlotHistCorrFactors(char* infile1,
2148                          std::string text1,
2149                          char* infile2,
2150                          std::string text2,
2151                          char* infile3,
2152                          std::string text3,
2153                          char* infile4,
2154                          std::string text4,
2155                          char* infile5,
2156                          std::string text5,
2157                          std::string prefixF,
2158                          bool ratio = false,
2159                          bool drawStatBox = true,
2160                          int nmin = 100,
2161                          bool dataMC = false,
2162                          int year = 2018,
2163                          int iformat = 0,
2164                          int save = 0) {
2165   std::map<int, cfactors> cfacs[5];
2166   std::vector<std::string> texts;
2167   int nfile(0), etamin(100), etamax(-100), maxdepth(0);
2168   const char* blank("");
2169   if (infile1 != blank) {
2170     readCorrFactors(infile1, 1.0, cfacs[nfile], etamin, etamax, maxdepth, iformat);
2171     if (cfacs[nfile].size() > 0) {
2172       texts.push_back(text1);
2173       ++nfile;
2174     }
2175   }
2176   if (infile2 != blank) {
2177     readCorrFactors(infile2, 1.0, cfacs[nfile], etamin, etamax, maxdepth, iformat);
2178     if (cfacs[nfile].size() > 0) {
2179       texts.push_back(text2);
2180       ++nfile;
2181     }
2182   }
2183   if (infile3 != blank) {
2184     readCorrFactors(infile3, 1.0, cfacs[nfile], etamin, etamax, maxdepth, iformat);
2185     if (cfacs[nfile].size() > 0) {
2186       texts.push_back(text3);
2187       ++nfile;
2188     }
2189   }
2190   if (infile4 != blank) {
2191     readCorrFactors(infile4, 1.0, cfacs[nfile], etamin, etamax, maxdepth, iformat);
2192     if (cfacs[nfile].size() > 0) {
2193       texts.push_back(text4);
2194       ++nfile;
2195     }
2196   }
2197   if (infile5 != blank) {
2198     readCorrFactors(infile5, 1.0, cfacs[nfile], etamin, etamax, maxdepth, iformat);
2199     if (cfacs[nfile].size() > 0) {
2200       texts.push_back(text5);
2201       ++nfile;
2202     }
2203   }
2204 
2205   if (nfile > 0) {
2206     gStyle->SetCanvasBorderMode(0);
2207     gStyle->SetCanvasColor(kWhite);
2208     gStyle->SetPadColor(kWhite);
2209     gStyle->SetFillColor(kWhite);
2210     gStyle->SetOptTitle(0);
2211     if ((!ratio) && drawStatBox) {
2212       gStyle->SetOptStat(10);
2213       gStyle->SetOptFit(10);
2214     } else {
2215       gStyle->SetOptStat(0);
2216       gStyle->SetOptFit(0);
2217     }
2218     int colors[6] = {1, 6, 4, 2, 7, 9};
2219     int mtype[6] = {20, 24, 22, 23, 21, 33};
2220     int nbin = etamax - etamin + 1;
2221     std::vector<TH1D*> hists;
2222     std::vector<int> entries, htype, depths;
2223     std::vector<double> fitr;
2224     char name[100];
2225     double dy(0);
2226     int fits(0);
2227     int nline(0);
2228     if (ratio) {
2229       for (int ih = 1; ih < nfile; ++ih) {
2230         for (int j = 0; j < maxdepth; ++j) {
2231           sprintf(name, "h%dd%d", ih, j + 1);
2232           TH1D* h = new TH1D(name, name, nbin, etamin, etamax);
2233           double sumNum(0), sumDen(0);
2234           std::map<int, cfactors>::const_iterator ktr = cfacs[ih].begin();
2235           for (std::map<int, cfactors>::const_iterator itr = cfacs[0].begin(); itr != cfacs[0].end(); ++itr, ++ktr) {
2236             int dep = (itr->second).depth;
2237             if (dep == j + 1) {
2238               int ieta = (itr->second).ieta;
2239               int bin = ieta - etamin + 1;
2240               float val = (itr->second).corrf / (ktr->second).corrf;
2241               float dvl =
2242                   val *
2243                   sqrt((((itr->second).dcorr * (itr->second).dcorr) / ((itr->second).corrf * (itr->second).corrf)) +
2244                        (((ktr->second).dcorr * (ktr->second).dcorr) / ((ktr->second).corrf * (ktr->second).corrf)));
2245               h->SetBinContent(bin, val);
2246               h->SetBinError(bin, dvl);
2247               sumNum += (val / (dvl * dvl));
2248               sumDen += (1.0 / (dvl * dvl));
2249             }
2250           }
2251           double fit = (sumDen > 0) ? (sumNum / sumDen) : 1.0;
2252           std::cout << "Fit to Pol0: " << fit << std::endl;
2253           h->SetLineColor(colors[ih]);
2254           h->SetMarkerColor(colors[ih]);
2255           h->SetMarkerStyle(mtype[j]);
2256           h->SetMarkerSize(0.9);
2257           h->GetXaxis()->SetTitle("i#eta");
2258           sprintf(name, "CF_{%s}/CF_{Set}", texts[0].c_str());
2259           h->GetYaxis()->SetTitle(name);
2260           h->GetYaxis()->SetLabelOffset(0.005);
2261           h->GetYaxis()->SetTitleSize(0.036);
2262           h->GetYaxis()->SetTitleOffset(1.20);
2263           h->GetYaxis()->SetRangeUser(0.80, 1.20);
2264           hists.push_back(h);
2265           fitr.push_back(fit);
2266           htype.push_back(ih);
2267           depths.push_back(j + 1);
2268         }
2269         if (ih == 1)
2270           nline += hists.size();
2271         else
2272           ++nline;
2273       }
2274     } else {
2275       for (int k1 = 0; k1 < nfile; ++k1) {
2276         for (int j = 0; j < maxdepth; ++j) {
2277           sprintf(name, "h%dd%d", k1, j + 1);
2278           TH1D* h = new TH1D(name, name, nbin, etamin, etamax);
2279           int nent(0);
2280           for (std::map<int, cfactors>::const_iterator itr = cfacs[k1].begin(); itr != cfacs[k1].end(); ++itr) {
2281             int dep = (itr->second).depth;
2282             if (dep == j + 1) {
2283               int ieta = (itr->second).ieta;
2284               int bin = ieta - etamin + 1;
2285               float val = (itr->second).corrf;
2286               float dvl = (itr->second).dcorr;
2287               h->SetBinContent(bin, val);
2288               h->SetBinError(bin, dvl);
2289               nent++;
2290             }
2291           }
2292           if (nent > nmin) {
2293             fits++;
2294             if (drawStatBox)
2295               dy += 0.025;
2296             sprintf(name, "h%ddf%d", k1, j + 1);
2297             TF1* func = new TF1(name, "pol0", etamin, etamax);
2298             h->Fit(func, "+QWLR", "");
2299           }
2300           h->SetLineColor(colors[k1]);
2301           h->SetMarkerColor(colors[k1]);
2302           h->SetMarkerStyle(mtype[j]);
2303           h->SetMarkerSize(0.9);
2304           h->GetXaxis()->SetTitle("i#eta");
2305           h->GetYaxis()->SetTitle("Correction Factor");
2306           h->GetYaxis()->SetLabelOffset(0.005);
2307           h->GetYaxis()->SetTitleOffset(1.20);
2308           h->GetYaxis()->SetRangeUser(0.5, 1.5);
2309           hists.push_back(h);
2310           entries.push_back(nent);
2311           if (drawStatBox)
2312             dy += 0.025;
2313           htype.push_back(k1);
2314           depths.push_back(j + 1);
2315         }
2316         if (k1 == 0)
2317           nline += hists.size();
2318         else
2319           ++nline;
2320       }
2321     }
2322     if (ratio)
2323       sprintf(name, "c_Corr%sRatio", prefixF.c_str());
2324     else
2325       sprintf(name, "c_Corr%s", prefixF.c_str());
2326     TCanvas* pad = new TCanvas(name, name, 700, 500);
2327     pad->SetRightMargin(0.10);
2328     pad->SetTopMargin(0.10);
2329     double yh = 0.90;
2330     double yl = yh - 0.035 * hists.size() - dy - 0.01;
2331     TLegend* legend = new TLegend(0.55, yl, 0.90, yl + 0.035 * nline);
2332     legend->SetFillColor(kWhite);
2333     for (unsigned int k = 0; k < hists.size(); ++k) {
2334       if (k == 0)
2335         hists[k]->Draw("");
2336       else
2337         hists[k]->Draw("sames");
2338       pad->Update();
2339       int k1 = htype[k];
2340       if (!ratio) {
2341         TPaveStats* st1 = (TPaveStats*)hists[k]->GetListOfFunctions()->FindObject("stats");
2342         if (st1 != nullptr) {
2343           dy = (entries[k] > nmin) ? 0.05 : 0.025;
2344           st1->SetLineColor(colors[k1]);
2345           st1->SetTextColor(colors[k1]);
2346           st1->SetY1NDC(yh - dy);
2347           st1->SetY2NDC(yh);
2348           st1->SetX1NDC(0.70);
2349           st1->SetX2NDC(0.90);
2350           yh -= dy;
2351         }
2352         sprintf(name, "Depth %d (%s)", depths[k], texts[k1].c_str());
2353       } else {
2354         sprintf(name, "Depth %d (%s Mean = %5.3f)", depths[k], texts[k1].c_str(), fitr[k]);
2355       }
2356       if ((depths[k] == 1) || (k1 == 0))
2357         legend->AddEntry(hists[k], name, "lp");
2358     }
2359     legend->Draw("same");
2360     TPaveText* txt0 = new TPaveText(0.12, 0.84, 0.49, 0.89, "blNDC");
2361     txt0->SetFillColor(0);
2362     char txt[40];
2363     if (dataMC)
2364       sprintf(txt, "CMS Preliminary (%d)", year);
2365     else
2366       sprintf(txt, "CMS Simulation Preliminary (%d)", year);
2367     txt0->AddText(txt);
2368     txt0->Draw("same");
2369     pad->Update();
2370     if (fits < 1) {
2371       double xmin = hists[0]->GetBinLowEdge(1);
2372       int nbin = hists[0]->GetNbinsX();
2373       double xmax = hists[0]->GetBinLowEdge(nbin) + hists[0]->GetBinWidth(nbin);
2374       TLine* line = new TLine(xmin, 1.0, xmax, 1.0);
2375       line->SetLineColor(9);
2376       line->SetLineWidth(2);
2377       line->SetLineStyle(2);
2378       line->Draw("same");
2379       pad->Update();
2380     }
2381     if (save > 0) {
2382       sprintf(name, "%s.pdf", pad->GetName());
2383       pad->Print(name);
2384     } else if (save < 0) {
2385       sprintf(name, "%s.C", pad->GetName());
2386       pad->Print(name);
2387     }
2388   }
2389 }
2390 
2391 void PlotHistCorrSys(std::string infilec, int conds, std::string text, int save = 0) {
2392   char fname[100];
2393   int iformat(0);
2394   sprintf(fname, "%s_cond0.txt", infilec.c_str());
2395   int etamin(100), etamax(-100), maxdepth(0);
2396   std::map<int, cfactors> cfacs;
2397   readCorrFactors(fname, 1.0, cfacs, etamin, etamax, maxdepth, iformat);
2398   // There are good records from the master file
2399   if (cfacs.size() > 0) {
2400     // Now read the other files
2401     std::map<int, cfactors> errfacs;
2402     for (int i = 0; i < conds; ++i) {
2403       sprintf(fname, "%s_cond%d.txt", infilec.c_str(), i + 1);
2404       std::map<int, cfactors> cfacx;
2405       int etamin1(100), etamax1(-100), maxdepth1(0);
2406       readCorrFactors(fname, 1.0, cfacx, etamin1, etamax1, maxdepth1, iformat);
2407       for (std::map<int, cfactors>::const_iterator itr1 = cfacx.begin(); itr1 != cfacx.end(); ++itr1) {
2408         std::map<int, cfactors>::iterator itr2 = errfacs.find(itr1->first);
2409         float corrf = (itr1->second).corrf;
2410         if (itr2 == errfacs.end()) {
2411           errfacs[itr1->first] = cfactors(1, 0, corrf, corrf * corrf);
2412         } else {
2413           int nent = (itr2->second).ieta + 1;
2414           float c1 = (itr2->second).corrf + corrf;
2415           float c2 = (itr2->second).dcorr + (corrf * corrf);
2416           errfacs[itr1->first] = cfactors(nent, 0, c1, c2);
2417         }
2418       }
2419     }
2420     // find the RMS from the distributions
2421     for (std::map<int, cfactors>::iterator itr = errfacs.begin(); itr != errfacs.end(); ++itr) {
2422       int nent = (itr->second).ieta;
2423       float mean = (itr->second).corrf / nent;
2424       float rms2 = (itr->second).dcorr / nent - (mean * mean);
2425       float rms = rms2 > 0 ? sqrt(rms2) : 0;
2426       errfacs[itr->first] = cfactors(nent, 0, mean, rms);
2427     }
2428     // Now combine the errors and plot
2429     int k(0);
2430     gStyle->SetCanvasBorderMode(0);
2431     gStyle->SetCanvasColor(kWhite);
2432     gStyle->SetPadColor(kWhite);
2433     gStyle->SetFillColor(kWhite);
2434     gStyle->SetOptTitle(0);
2435     gStyle->SetOptStat(10);
2436     gStyle->SetOptFit(10);
2437     int colors[6] = {1, 6, 4, 7, 2, 9};
2438     int mtype[6] = {20, 21, 22, 23, 24, 33};
2439     std::vector<TH1D*> hists;
2440     char name[100];
2441     int nbin = etamax - etamin + 1;
2442     for (int j = 0; j < maxdepth; ++j) {
2443       sprintf(name, "hd%d", j + 1);
2444       TH1D* h = new TH1D(name, name, nbin, etamin, etamax);
2445       h->SetLineColor(colors[j]);
2446       h->SetMarkerColor(colors[j]);
2447       h->SetMarkerStyle(mtype[j]);
2448       h->GetXaxis()->SetTitle("i#eta");
2449       h->GetYaxis()->SetTitle("Correction Factor");
2450       h->GetYaxis()->SetLabelOffset(0.005);
2451       h->GetYaxis()->SetTitleOffset(1.20);
2452       h->GetYaxis()->SetRangeUser(0.0, 2.0);
2453       hists.push_back(h);
2454     }
2455     for (std::map<int, cfactors>::iterator itr = cfacs.begin(); itr != cfacs.end(); ++itr, ++k) {
2456       std::map<int, cfactors>::iterator itr2 = errfacs.find(itr->first);
2457       float mean2 = (itr2 == errfacs.end()) ? 0 : (itr2->second).corrf;
2458       float ersys = (itr2 == errfacs.end()) ? 0 : (itr2->second).dcorr;
2459       float erstt = (itr->second).dcorr;
2460       float ertot = sqrt(erstt * erstt + ersys * ersys);
2461       float mean = (itr->second).corrf;
2462       int ieta = (itr->second).ieta;
2463       int depth = (itr->second).depth;
2464       std::cout << "[" << k << "] " << ieta << " " << depth << " " << mean << ":" << mean2 << " " << erstt << ":"
2465                 << ersys << ":" << ertot << std::endl;
2466       int bin = ieta - etamin + 1;
2467       hists[depth - 1]->SetBinContent(bin, mean);
2468       hists[depth - 1]->SetBinError(bin, ertot);
2469     }
2470     TCanvas* pad = new TCanvas("CorrFactorSys", "CorrFactorSys", 700, 500);
2471     pad->SetRightMargin(0.10);
2472     pad->SetTopMargin(0.10);
2473     double yh = 0.90;
2474     double yl = yh - 0.050 * hists.size() - 0.01;
2475     TLegend* legend = new TLegend(0.60, yl, 0.90, yl + 0.025 * hists.size());
2476     legend->SetFillColor(kWhite);
2477     for (unsigned int k = 0; k < hists.size(); ++k) {
2478       if (k == 0)
2479         hists[k]->Draw("");
2480       else
2481         hists[k]->Draw("sames");
2482       pad->Update();
2483       TPaveStats* st1 = (TPaveStats*)hists[k]->GetListOfFunctions()->FindObject("stats");
2484       if (st1 != nullptr) {
2485         st1->SetLineColor(colors[k]);
2486         st1->SetTextColor(colors[k]);
2487         st1->SetY1NDC(yh - 0.025);
2488         st1->SetY2NDC(yh);
2489         st1->SetX1NDC(0.70);
2490         st1->SetX2NDC(0.90);
2491         yh -= 0.025;
2492       }
2493       sprintf(name, "Depth %d (%s)", k + 1, text.c_str());
2494       legend->AddEntry(hists[k], name, "lp");
2495     }
2496     legend->Draw("same");
2497     pad->Update();
2498     if (save > 0) {
2499       sprintf(name, "%s.pdf", pad->GetName());
2500       pad->Print(name);
2501     } else if (save < 0) {
2502       sprintf(name, "%s.C", pad->GetName());
2503       pad->Print(name);
2504     }
2505   }
2506 }
2507 
2508 void PlotHistCorrLumis(std::string infilec, int conds, double lumi, int save = 0) {
2509   char fname[100];
2510   int iformat(0);
2511   sprintf(fname, "%s_0.txt", infilec.c_str());
2512   std::map<int, cfactors> cfacs;
2513   int etamin(100), etamax(-100), maxdepth(0);
2514   readCorrFactors(fname, 1.0, cfacs, etamin, etamax, maxdepth, iformat);
2515   int nbin = etamax - etamin + 1;
2516   std::cout << "Max Depth " << maxdepth << " and " << nbin << " eta bins for " << etamin << ":" << etamax << std::endl;
2517 
2518   // There are good records from the master file
2519   int colors[8] = {4, 2, 6, 7, 1, 9, 3, 5};
2520   int mtype[8] = {20, 21, 22, 23, 24, 25, 26, 27};
2521   if (cfacs.size() > 0) {
2522     // Now read the other files
2523     std::vector<TH1D*> hists;
2524     char name[100];
2525     for (int i = 0; i < conds; ++i) {
2526       int ih = (int)(hists.size());
2527       for (int j = 0; j < maxdepth; ++j) {
2528         sprintf(name, "hd%d%d", j + 1, i);
2529         TH1D* h = new TH1D(name, name, nbin, etamin, etamax);
2530         h->SetLineColor(colors[j]);
2531         h->SetMarkerColor(colors[j]);
2532         h->SetMarkerStyle(mtype[i]);
2533         h->SetMarkerSize(0.9);
2534         h->GetXaxis()->SetTitle("i#eta");
2535         h->GetYaxis()->SetTitle("Fractional Error");
2536         h->GetYaxis()->SetLabelOffset(0.005);
2537         h->GetYaxis()->SetTitleOffset(1.20);
2538         h->GetYaxis()->SetRangeUser(0.0, 0.10);
2539         hists.push_back(h);
2540       }
2541       sprintf(fname, "%s_%d.txt", infilec.c_str(), i);
2542       int etamin1(100), etamax1(-100), maxdepth1(0);
2543       readCorrFactors(fname, 1.0, cfacs, etamin1, etamax1, maxdepth1, iformat);
2544       for (std::map<int, cfactors>::const_iterator itr = cfacs.begin(); itr != cfacs.end(); ++itr) {
2545         double value = (itr->second).dcorr / (itr->second).corrf;
2546         int bin = (itr->second).ieta - etamin + 1;
2547         hists[ih + (itr->second).depth - 1]->SetBinContent(bin, value);
2548         hists[ih + (itr->second).depth - 1]->SetBinError(bin, 0.0001);
2549       }
2550     }
2551 
2552     // Now plot the histograms
2553     gStyle->SetCanvasBorderMode(0);
2554     gStyle->SetCanvasColor(kWhite);
2555     gStyle->SetPadColor(kWhite);
2556     gStyle->SetFillColor(kWhite);
2557     gStyle->SetOptTitle(0);
2558     gStyle->SetOptStat(0);
2559     gStyle->SetOptFit(0);
2560     TCanvas* pad = new TCanvas("CorrFactorErr", "CorrFactorErr", 700, 500);
2561     pad->SetRightMargin(0.10);
2562     pad->SetTopMargin(0.10);
2563     double yh(0.89);
2564     TLegend* legend = new TLegend(0.60, yh - 0.04 * conds, 0.89, yh);
2565     legend->SetFillColor(kWhite);
2566     double lumic(lumi);
2567     for (unsigned int k = 0; k < hists.size(); ++k) {
2568       if (k == 0)
2569         hists[k]->Draw("");
2570       else
2571         hists[k]->Draw("sames");
2572       pad->Update();
2573       if (k % maxdepth == 0) {
2574         sprintf(name, "L = %5.2f fb^{-1}", lumic);
2575         legend->AddEntry(hists[k], name, "lp");
2576         lumic *= 0.5;
2577       }
2578     }
2579     legend->Draw("same");
2580     pad->Update();
2581     if (save > 0) {
2582       sprintf(name, "%s.pdf", pad->GetName());
2583       pad->Print(name);
2584     } else if (save < 0) {
2585       sprintf(name, "%s.C", pad->GetName());
2586       pad->Print(name);
2587     }
2588   }
2589 }
2590 
2591 void PlotHistCorrRel(char* infile1,
2592                      char* infile2,
2593                      std::string text1,
2594                      std::string text2,
2595                      int iformat1 = 0,
2596                      int iformat2 = 0,
2597                      int save = 0) {
2598   std::map<int, cfactors> cfacs1, cfacs2;
2599   int etamin(100), etamax(-100), maxdepth(0);
2600   readCorrFactors(infile1, 1.0, cfacs1, etamin, etamax, maxdepth, iformat1);
2601   readCorrFactors(infile2, 1.0, cfacs2, etamin, etamax, maxdepth, iformat2);
2602   std::map<int, std::pair<cfactors, cfactors> > cfacs;
2603   for (std::map<int, cfactors>::iterator itr = cfacs1.begin(); itr != cfacs1.end(); ++itr) {
2604     std::map<int, cfactors>::iterator ktr = cfacs2.find(itr->first);
2605     if (ktr == cfacs2.end()) {
2606       cfactors fac2(((itr->second).ieta), ((itr->second).depth), 0, -1);
2607       cfacs[itr->first] = std::pair<cfactors, cfactors>((itr->second), fac2);
2608     } else {
2609       cfactors fac2(ktr->second);
2610       cfacs[itr->first] = std::pair<cfactors, cfactors>((itr->second), fac2);
2611     }
2612   }
2613   for (std::map<int, cfactors>::iterator itr = cfacs2.begin(); itr != cfacs2.end(); ++itr) {
2614     std::map<int, cfactors>::const_iterator ktr = cfacs1.find(itr->first);
2615     if (ktr == cfacs1.end()) {
2616       cfactors fac1(((itr->second).ieta), ((itr->second).depth), 0, -1);
2617       cfacs[itr->first] = std::pair<cfactors, cfactors>(fac1, (itr->second));
2618     }
2619   }
2620 
2621   // There are good records in bothe the files
2622   if ((cfacs1.size() > 0) && (cfacs2.size() > 0)) {
2623     int k(0);
2624     gStyle->SetCanvasBorderMode(0);
2625     gStyle->SetCanvasColor(kWhite);
2626     gStyle->SetPadColor(kWhite);
2627     gStyle->SetFillColor(kWhite);
2628     gStyle->SetOptTitle(0);
2629     gStyle->SetOptStat(10);
2630     gStyle->SetOptFit(10);
2631     int colors[6] = {1, 6, 4, 7, 2, 9};
2632     int mtype[6] = {20, 21, 22, 23, 24, 33};
2633     std::vector<TH1D*> hists;
2634     char name[100];
2635     int nbin = etamax - etamin + 1;
2636     for (int i = 0; i < 2; ++i) {
2637       for (int j = 0; j < maxdepth; ++j) {
2638         int j1 = (i == 0) ? j : maxdepth + j;
2639         sprintf(name, "hd%d%d", i, j + 1);
2640         TH1D* h = new TH1D(name, name, nbin, etamin, etamax);
2641         h->SetLineColor(colors[j1]);
2642         h->SetMarkerColor(colors[j1]);
2643         h->SetMarkerStyle(mtype[i]);
2644         h->GetXaxis()->SetTitle("i#eta");
2645         h->GetYaxis()->SetTitle("Correction Factor");
2646         h->GetYaxis()->SetLabelOffset(0.005);
2647         h->GetYaxis()->SetTitleOffset(1.20);
2648         h->GetYaxis()->SetRangeUser(0.0, 2.0);
2649         hists.push_back(h);
2650       }
2651     }
2652     for (std::map<int, std::pair<cfactors, cfactors> >::iterator it = cfacs.begin(); it != cfacs.end(); ++it, ++k) {
2653       float mean1 = (it->second).first.corrf;
2654       float error1 = (it->second).first.dcorr;
2655       float mean2 = (it->second).second.corrf;
2656       float error2 = (it->second).second.dcorr;
2657       int ieta = (it->second).first.ieta;
2658       int depth = (it->second).first.depth;
2659       /*
2660       std::cout << "[" << k << "] " << ieta << " " << depth << " " 
2661         << mean1 << ":" << mean2 << " " << error1 << ":" << error2
2662         << std::endl;
2663       */
2664       int bin = ieta - etamin + 1;
2665       if (error1 >= 0) {
2666         hists[depth - 1]->SetBinContent(bin, mean1);
2667         hists[depth - 1]->SetBinError(bin, error1);
2668       }
2669       if (error2 >= 0) {
2670         hists[maxdepth + depth - 1]->SetBinContent(bin, mean2);
2671         hists[maxdepth + depth - 1]->SetBinError(bin, error2);
2672       }
2673     }
2674     TCanvas* pad = new TCanvas("CorrFactors", "CorrFactors", 700, 500);
2675     pad->SetRightMargin(0.10);
2676     pad->SetTopMargin(0.10);
2677     double yh = 0.90;
2678     double yl = yh - 0.050 * hists.size() - 0.01;
2679     TLegend* legend = new TLegend(0.60, yl, 0.90, yl + 0.025 * hists.size());
2680     legend->SetFillColor(kWhite);
2681     for (unsigned int k = 0; k < hists.size(); ++k) {
2682       if (k == 0)
2683         hists[k]->Draw("");
2684       else
2685         hists[k]->Draw("sames");
2686       pad->Update();
2687       TPaveStats* st1 = (TPaveStats*)hists[k]->GetListOfFunctions()->FindObject("stats");
2688       if (st1 != nullptr) {
2689         st1->SetLineColor(colors[k]);
2690         st1->SetTextColor(colors[k]);
2691         st1->SetY1NDC(yh - 0.025);
2692         st1->SetY2NDC(yh);
2693         st1->SetX1NDC(0.70);
2694         st1->SetX2NDC(0.90);
2695         yh -= 0.025;
2696       }
2697       if (k < (unsigned int)(maxdepth)) {
2698         sprintf(name, "Depth %d (%s)", k + 1, text1.c_str());
2699       } else {
2700         sprintf(name, "Depth %d (%s)", k - maxdepth + 1, text2.c_str());
2701       }
2702       legend->AddEntry(hists[k], name, "lp");
2703     }
2704     legend->Draw("same");
2705     pad->Update();
2706     if (save > 0) {
2707       sprintf(name, "%s.pdf", pad->GetName());
2708       pad->Print(name);
2709     } else if (save < 0) {
2710       sprintf(name, "%s.C", pad->GetName());
2711       pad->Print(name);
2712     }
2713   }
2714 }
2715 
2716 void PlotHistCorrDepth(char* infile1,
2717                        char* infile2,
2718                        std::string text1,
2719                        std::string text2,
2720                        int depth,
2721                        int ietamax,
2722                        int iformat1 = 0,
2723                        int iformat2 = 0,
2724                        int save = 0,
2725                        int debug = 1) {
2726   std::map<int, cfactors> cfacs1, cfacs2;
2727   int etamin(100), etamax(-100), maxdepth(0);
2728   readCorrFactors(infile1, 1.0, cfacs1, etamin, etamax, maxdepth, iformat1, (debug > 1));
2729   readCorrFactors(infile2, 1.0, cfacs2, etamin, etamax, maxdepth, iformat2, (debug > 1));
2730 
2731   double sumNum(0), sumDen(0), ratMax(0);
2732   int npt0(0), npt1(0);
2733   for (std::map<int, cfactors>::iterator itr = cfacs1.begin(); itr != cfacs1.end(); ++itr) {
2734     std::map<int, cfactors>::iterator ktr = cfacs2.find(itr->first);
2735     if ((ktr != cfacs2.end()) && ((itr->second).depth == depth)) {
2736       double er1 = (itr->second).dcorr / (itr->second).corrf;
2737       double er2 = (ktr->second).dcorr / (ktr->second).corrf;
2738       double tmp = (ktr->second).corrf / (itr->second).corrf;
2739       double dtmp = tmp * sqrt(er1 * er1 + er2 * er2);
2740       double rat = (tmp > 1.0) ? 1.0 / tmp : tmp;
2741       double drt = (tmp > 1.0) ? dtmp / (tmp * tmp) : dtmp;
2742       rat = fabs(1.0 - rat);
2743       ratMax = std::max(ratMax, rat);
2744       ++npt0;
2745       if (debug > 0)
2746         std::cout << std::hex << (itr->first) << std::dec << " ieta:depth" << (itr->second).ieta << ":"
2747                   << (itr->second).depth << " Corr " << (itr->second).corrf << ":" << (ktr->second).corrf << " Ratio "
2748                   << rat << ":" << drt << std::endl;
2749       if (std::abs((itr->second).ieta) <= ietamax) {
2750         sumNum += (rat / (drt * drt));
2751         sumDen += (1.0 / (drt * drt));
2752         ++npt1;
2753       }
2754     }
2755   }
2756   sumNum = (sumDen > 0) ? (sumNum / sumDen) : 0;
2757   sumDen = (sumDen > 0) ? 1.0 / sqrt(sumDen) : 0;
2758   std::cout << "Get Ratio of mean for " << npt0 << ":" << npt1 << " points for depth " << depth << " Mean " << sumNum
2759             << " +- " << sumDen << " (Maximum " << ratMax << ")" << std::endl;
2760 
2761   gStyle->SetCanvasBorderMode(0);
2762   gStyle->SetCanvasColor(kWhite);
2763   gStyle->SetPadColor(kWhite);
2764   gStyle->SetFillColor(kWhite);
2765   gStyle->SetOptTitle(0);
2766   gStyle->SetOptStat(0);
2767   gStyle->SetOptFit(0);
2768   int colors[2] = {1, 2};
2769   int mtype[2] = {20, 24};
2770   int nbin = etamax - etamin + 1;
2771   std::vector<TH1D*> hists;
2772   char name[100];
2773   for (int j = 0; j < 2; ++j) {
2774     sprintf(name, "hd%d", (j + 1));
2775     TH1D* h = new TH1D(name, name, nbin, etamin, etamax);
2776     if (j == 0) {
2777       for (std::map<int, cfactors>::const_iterator itr = cfacs1.begin(); itr != cfacs1.end(); ++itr) {
2778         if ((itr->second).depth == depth) {
2779           int ieta = (itr->second).ieta;
2780           int bin = ieta - etamin + 1;
2781           float val = (itr->second).corrf;
2782           float dvl = (itr->second).dcorr;
2783           h->SetBinContent(bin, val);
2784           h->SetBinError(bin, dvl);
2785         }
2786       }
2787     } else {
2788       for (std::map<int, cfactors>::const_iterator itr = cfacs2.begin(); itr != cfacs2.end(); ++itr) {
2789         if ((itr->second).depth == depth) {
2790           int ieta = (itr->second).ieta;
2791           int bin = ieta - etamin + 1;
2792           float val = (itr->second).corrf;
2793           float dvl = (itr->second).dcorr;
2794           h->SetBinContent(bin, val);
2795           h->SetBinError(bin, dvl);
2796         }
2797       }
2798     }
2799     h->SetLineColor(colors[j]);
2800     h->SetMarkerColor(colors[j]);
2801     h->SetMarkerStyle(mtype[j]);
2802     h->GetXaxis()->SetTitle("i#eta");
2803     h->GetYaxis()->SetTitle("Correction Factor");
2804     h->GetYaxis()->SetLabelOffset(0.005);
2805     h->GetYaxis()->SetTitleOffset(1.20);
2806     h->GetYaxis()->SetRangeUser(0.0, 2.0);
2807     hists.push_back(h);
2808   }
2809   sprintf(name, "c_CorrFactorDepth%d", depth);
2810   TCanvas* pad = new TCanvas(name, name, 700, 500);
2811   pad->SetRightMargin(0.10);
2812   pad->SetTopMargin(0.10);
2813   double yl = 0.25;
2814   TLegend* legend = new TLegend(0.25, yl, 0.85, yl + 0.04 * hists.size());
2815   legend->SetFillColor(kWhite);
2816   for (unsigned int k = 0; k < hists.size(); ++k) {
2817     if (k == 0) {
2818       hists[k]->Draw("");
2819       sprintf(name, "Depth %d (%s)", depth, text1.c_str());
2820     } else {
2821       hists[k]->Draw("sames");
2822       sprintf(name, "Depth %d (%s)", depth, text2.c_str());
2823     }
2824     pad->Update();
2825     legend->AddEntry(hists[k], name, "lp");
2826   }
2827   legend->Draw("same");
2828   pad->Update();
2829   if (save > 0) {
2830     sprintf(name, "%s.pdf", pad->GetName());
2831     pad->Print(name);
2832   } else if (save < 0) {
2833     sprintf(name, "%s.C", pad->GetName());
2834     pad->Print(name);
2835   }
2836 }
2837 
2838 void PlotFourHists(std::string infile,
2839                    std::string prefix0,
2840                    int type = 0,
2841                    int drawStatBox = 0,
2842                    bool normalize = false,
2843                    int save = 0,
2844                    std::string prefix1 = "",
2845                    std::string text1 = "",
2846                    std::string prefix2 = "",
2847                    std::string text2 = "",
2848                    std::string prefix3 = "",
2849                    std::string text3 = "",
2850                    std::string prefix4 = "",
2851                    std::string text4 = "") {
2852   int colors[4] = {2, 4, 6, 1};
2853   std::string names[5] = {"eta03", "eta13", "eta23", "eta33", "eta43"};
2854   std::string xtitle[5] = {"i#eta", "i#eta", "i#eta", "i#eta", "i#eta"};
2855   std::string ytitle[5] = {"Tracks", "Tracks", "Tracks", "Tracks", "Tracks"};
2856   std::string title[5] = {"All Tracks (p = 40:60 GeV)",
2857                           "Good Quality Tracks (p = 40:60 GeV)",
2858                           "Selected Tracks (p = 40:60 GeV)",
2859                           "Isolated Tracks (p = 40:60 GeV)",
2860                           "Isolated MIP Tracks (p = 40:60 GeV)"};
2861 
2862   gStyle->SetCanvasBorderMode(0);
2863   gStyle->SetCanvasColor(kWhite);
2864   gStyle->SetPadColor(kWhite);
2865   gStyle->SetFillColor(kWhite);
2866   gStyle->SetOptTitle(0);
2867   gStyle->SetOptFit(0);
2868   if (drawStatBox == 0)
2869     gStyle->SetOptStat(0);
2870   else
2871     gStyle->SetOptStat(1110);
2872 
2873   if (type < 0 || type > 4)
2874     type = 0;
2875   char name[100], namep[100];
2876   TFile* file = new TFile(infile.c_str());
2877 
2878   std::vector<TH1D*> hists;
2879   std::vector<int> kks;
2880   std::vector<std::string> texts;
2881   for (int k = 0; k < 4; ++k) {
2882     std::string prefix, text;
2883     if (k == 0) {
2884       prefix = prefix1;
2885       text = text1;
2886     } else if (k == 1) {
2887       prefix = prefix2;
2888       text = text2;
2889     } else if (k == 2) {
2890       prefix = prefix3;
2891       text = text3;
2892     } else {
2893       prefix = prefix4;
2894       text = text4;
2895     }
2896     if (prefix != "") {
2897       sprintf(name, "%s%s", prefix.c_str(), names[type].c_str());
2898       TH1D* hist1 = (TH1D*)file->FindObjectAny(name);
2899       if (hist1 != nullptr) {
2900         hists.push_back((TH1D*)(hist1->Clone()));
2901         kks.push_back(k);
2902         texts.push_back(text);
2903       }
2904     }
2905   }
2906   if (hists.size() > 0) {
2907     sprintf(namep, "c_%s%s", prefix0.c_str(), names[type].c_str());
2908     double ymax(0.90), dy(0.13);
2909     double ymx0 = (drawStatBox == 0) ? (ymax - .01) : (ymax - dy * hists.size() - .01);
2910     TCanvas* pad = new TCanvas(namep, namep, 700, 500);
2911     TLegend* legend = new TLegend(0.64, ymx0 - 0.05 * hists.size(), 0.89, ymx0);
2912     legend->SetFillColor(kWhite);
2913     pad->SetRightMargin(0.10);
2914     pad->SetTopMargin(0.10);
2915     for (unsigned int jk = 0; jk < hists.size(); ++jk) {
2916       int k = kks[jk];
2917       hists[jk]->GetXaxis()->SetTitleSize(0.040);
2918       hists[jk]->GetXaxis()->SetTitle(xtitle[type].c_str());
2919       hists[jk]->GetYaxis()->SetTitle(ytitle[type].c_str());
2920       hists[jk]->GetYaxis()->SetLabelOffset(0.005);
2921       hists[jk]->GetYaxis()->SetLabelSize(0.035);
2922       hists[jk]->GetYaxis()->SetTitleSize(0.040);
2923       hists[jk]->GetYaxis()->SetTitleOffset(1.15);
2924       hists[jk]->SetMarkerStyle(20);
2925       hists[jk]->SetMarkerColor(colors[k]);
2926       hists[jk]->SetLineColor(colors[k]);
2927       if (normalize) {
2928         if (jk == 0)
2929           hists[jk]->DrawNormalized();
2930         else
2931           hists[jk]->DrawNormalized("sames");
2932       } else {
2933         if (jk == 0)
2934           hists[jk]->Draw();
2935         else
2936           hists[jk]->Draw("sames");
2937       }
2938       pad->Update();
2939       TPaveStats* st1 = (TPaveStats*)hists[jk]->GetListOfFunctions()->FindObject("stats");
2940       if (st1 != nullptr) {
2941         double ymin = ymax - dy;
2942         st1->SetLineColor(colors[k]);
2943         st1->SetTextColor(colors[k]);
2944         st1->SetY1NDC(ymin);
2945         st1->SetY2NDC(ymax);
2946         st1->SetX1NDC(0.70);
2947         st1->SetX2NDC(0.90);
2948         ymax = ymin;
2949       }
2950       sprintf(name, "%s", texts[jk].c_str());
2951       legend->AddEntry(hists[jk], name, "lp");
2952     }
2953     legend->Draw("same");
2954     pad->Update();
2955     TPaveText* txt1 = new TPaveText(0.10, 0.905, 0.80, 0.95, "blNDC");
2956     txt1->SetFillColor(0);
2957     char txt[100];
2958     sprintf(txt, "%s", title[type].c_str());
2959     txt1->AddText(txt);
2960     txt1->Draw("same");
2961     /*
2962     TPaveText *txt2 = new TPaveText(0.11,0.825,0.33,0.895,"blNDC");
2963     txt2->SetFillColor(0);
2964     sprintf (txt, "CMS Preliminary");
2965     txt2->AddText(txt);
2966     txt2->Draw("same");
2967     */
2968     pad->Modified();
2969     pad->Update();
2970     if (save > 0) {
2971       sprintf(name, "%s.pdf", pad->GetName());
2972       pad->Print(name);
2973     } else if (save < 0) {
2974       sprintf(name, "%s.C", pad->GetName());
2975       pad->Print(name);
2976     }
2977   }
2978 }
2979 
2980 void PlotPUCorrHists(std::string infile = "corrfac.root",
2981                      std::string prefix = "",
2982                      int drawStatBox = 0,
2983                      bool approve = true,
2984                      int save = 0) {
2985   std::string name1[4] = {"W0", "W1", "W2", "P"};
2986   std::string name2[4] = {"All", "Barrel", "Endcap", ""};
2987   std::string name3[2] = {"", "p = 40:60 GeV"};
2988   std::string name4[2] = {"Loose Isolation", "Tight Isolation"};
2989   std::string xtitle[4] = {"Correction Factor", "Correction Factor", "Correction Factor", "i#eta"};
2990   std::string ytitle[4] = {"Tracks", "Tracks", "Tracks", "Correction Factor"};
2991 
2992   gStyle->SetCanvasBorderMode(0);
2993   gStyle->SetCanvasColor(kWhite);
2994   gStyle->SetPadColor(kWhite);
2995   gStyle->SetFillColor(kWhite);
2996   gStyle->SetOptTitle(0);
2997   gStyle->SetOptFit(0);
2998   if (drawStatBox == 0)
2999     gStyle->SetOptStat(0);
3000   else
3001     gStyle->SetOptStat(1110);
3002 
3003   char name[100], namep[100], title[100];
3004   TFile* file = new TFile(infile.c_str());
3005 
3006   if (file != nullptr) {
3007     for (int i1 = 0; i1 < 4; ++i1) {
3008       for (int i2 = 0; i2 < 2; ++i2) {
3009         for (int i3 = 0; i3 < 2; ++i3) {
3010           sprintf(name, "%s%d%d", name1[i1].c_str(), i2, i3);
3011           if (i2 == 0)
3012             sprintf(title, "%s Tracks Selected with %s", name2[i1].c_str(), name4[i3].c_str());
3013           else
3014             sprintf(title, "%s Tracks Selected with %s (%s)", name2[i1].c_str(), name4[i3].c_str(), name3[i2].c_str());
3015           TH1D* hist1(nullptr);
3016           TProfile* hist2(nullptr);
3017           if (i1 != 3) {
3018             TH1D* hist = (TH1D*)file->FindObjectAny(name);
3019             if (hist != nullptr) {
3020               hist1 = (TH1D*)(hist->Clone());
3021               hist1->GetXaxis()->SetTitleSize(0.040);
3022               hist1->GetXaxis()->SetTitle(xtitle[i1].c_str());
3023               hist1->GetYaxis()->SetTitle(ytitle[i1].c_str());
3024               hist1->GetYaxis()->SetLabelOffset(0.005);
3025               hist1->GetYaxis()->SetLabelSize(0.035);
3026               hist1->GetYaxis()->SetTitleSize(0.040);
3027               hist1->GetYaxis()->SetTitleOffset(1.15);
3028             }
3029           } else {
3030             TProfile* hist = (TProfile*)file->FindObjectAny(name);
3031             if (hist != nullptr) {
3032               hist2 = (TProfile*)(hist->Clone());
3033               hist2->GetXaxis()->SetTitleSize(0.040);
3034               hist2->GetXaxis()->SetTitle(xtitle[i1].c_str());
3035               hist2->GetYaxis()->SetTitle(ytitle[i1].c_str());
3036               hist2->GetYaxis()->SetLabelOffset(0.005);
3037               hist2->GetYaxis()->SetLabelSize(0.035);
3038               hist2->GetYaxis()->SetTitleSize(0.040);
3039               hist2->GetYaxis()->SetTitleOffset(1.15);
3040               //          hist2->GetYaxis()->SetRangeUser(0.0, 1.5);
3041               hist2->SetMarkerStyle(20);
3042             }
3043           }
3044           if ((hist1 != nullptr) || (hist2 != nullptr)) {
3045             sprintf(namep, "c_%s%s", name, prefix.c_str());
3046             TCanvas* pad = new TCanvas(namep, namep, 700, 500);
3047             pad->SetRightMargin(0.10);
3048             pad->SetTopMargin(0.10);
3049             if (hist1 != nullptr) {
3050               pad->SetLogy();
3051               hist1->Draw();
3052               pad->Update();
3053               TPaveStats* st1 = (TPaveStats*)hist1->GetListOfFunctions()->FindObject("stats");
3054               if (st1 != nullptr) {
3055                 st1->SetY1NDC(0.77);
3056                 st1->SetY2NDC(0.90);
3057                 st1->SetX1NDC(0.70);
3058                 st1->SetX2NDC(0.90);
3059               }
3060             } else {
3061               hist2->Draw();
3062               pad->Update();
3063             }
3064             TPaveText* txt1 = new TPaveText(0.10, 0.905, 0.80, 0.95, "blNDC");
3065             txt1->SetFillColor(0);
3066             char txt[100];
3067             sprintf(txt, "%s", title);
3068             txt1->AddText(txt);
3069             txt1->Draw("same");
3070             if (approve) {
3071               double xoff = (i1 == 3) ? 0.11 : 0.22;
3072               TPaveText* txt2 = new TPaveText(xoff, 0.825, xoff + 0.22, 0.895, "blNDC");
3073               txt2->SetFillColor(0);
3074               sprintf(txt, "CMS Preliminary");
3075               txt2->AddText(txt);
3076               txt2->Draw("same");
3077             }
3078             pad->Modified();
3079             pad->Update();
3080             if (save > 0) {
3081               sprintf(name, "%s.pdf", pad->GetName());
3082               pad->Print(name);
3083             } else if (save < 0) {
3084               sprintf(name, "%s.C", pad->GetName());
3085               pad->Print(name);
3086             }
3087           }
3088         }
3089       }
3090     }
3091   }
3092 }
3093 
3094 void PlotHistCorr(const char* infile,
3095                   std::string prefix,
3096                   std::string text0,
3097                   int eta = 0,
3098                   int mode = 1,
3099                   bool drawStatBox = true,
3100                   int save = 0) {
3101   gStyle->SetCanvasBorderMode(0);
3102   gStyle->SetCanvasColor(kWhite);
3103   gStyle->SetPadColor(kWhite);
3104   gStyle->SetFillColor(kWhite);
3105   gStyle->SetOptTitle(0);
3106   if (drawStatBox)
3107     gStyle->SetOptStat(1100);
3108   else
3109     gStyle->SetOptStat(0);
3110 
3111   std::string tags[3] = {"UnNoPU", "UnPU", "Cor"};
3112   std::string text[3] = {"Uncorrected no PU", "Uncorrected PU", "Corrected PU"};
3113   int colors[3] = {1, 4, 2};
3114   int styles[3] = {1, 3, 2};
3115   TFile* file = new TFile(infile);
3116   if (mode < 0 || mode > 2)
3117     mode = 1;
3118   int etamin = (eta == 0) ? -27 : eta;
3119   int etamax = (eta == 0) ? 27 : eta;
3120   for (int ieta = etamin; ieta <= etamax; ++ieta) {
3121     char name[20];
3122     double yh(0.90), dy(0.09);
3123     double yh1 = drawStatBox ? (yh - 3 * dy - 0.01) : (yh - 0.01);
3124     TLegend* legend = new TLegend(0.55, yh1 - 0.15, 0.89, yh1);
3125     legend->SetFillColor(kWhite);
3126     sprintf(name, "c_%sEovp%d", prefix.c_str(), ieta);
3127     TCanvas* pad = new TCanvas(name, name, 700, 500);
3128     pad->SetRightMargin(0.10);
3129     pad->SetTopMargin(0.10);
3130     TH1D* hist[3];
3131     double ymax(0);
3132     for (int k = 0; k < 3; ++k) {
3133       if (k < 2)
3134         sprintf(name, "EovP_ieta%d%s", ieta, tags[k].c_str());
3135       else
3136         sprintf(name, "EovP_ieta%dCor%dPU", ieta, mode);
3137       TH1D* hist1 = (TH1D*)file->FindObjectAny(name);
3138       if (hist1 != nullptr) {
3139         hist[k] = (TH1D*)(hist1->Clone());
3140         ymax = std::max(ymax, (hist1->GetMaximum()));
3141       }
3142     }
3143     int imax = 10 * (2 + int(0.1 * ymax));
3144     for (int k = 0; k < 3; ++k) {
3145       hist[k]->GetYaxis()->SetLabelOffset(0.005);
3146       hist[k]->GetYaxis()->SetTitleOffset(1.20);
3147       hist[k]->GetXaxis()->SetTitle("E/p");
3148       hist[k]->GetYaxis()->SetTitle("Tracks");
3149       hist[k]->SetLineColor(colors[k]);
3150       hist[k]->SetLineStyle(styles[k]);
3151       hist[k]->GetYaxis()->SetRangeUser(0.0, imax);
3152       if (k == 0)
3153         hist[k]->Draw();
3154       else
3155         hist[k]->Draw("sames");
3156       legend->AddEntry(hist[k], text[k].c_str(), "lp");
3157       pad->Update();
3158       if (drawStatBox) {
3159         TPaveStats* st1 = (TPaveStats*)hist[k]->GetListOfFunctions()->FindObject("stats");
3160         if (st1 != nullptr) {
3161           st1->SetLineColor(colors[k]);
3162           st1->SetTextColor(colors[k]);
3163           st1->SetY1NDC(yh - dy);
3164           st1->SetY2NDC(yh);
3165           st1->SetX1NDC(0.70);
3166           st1->SetX2NDC(0.90);
3167           yh -= dy;
3168         }
3169       }
3170     }
3171     pad->Update();
3172     legend->Draw("same");
3173     pad->Update();
3174     TPaveText* txt1 = new TPaveText(0.10, 0.905, 0.80, 0.95, "blNDC");
3175     txt1->SetFillColor(0);
3176     char title[100];
3177     sprintf(title, "%s for i#eta = %d", text0.c_str(), ieta);
3178     txt1->AddText(title);
3179     txt1->Draw("same");
3180     pad->Modified();
3181     pad->Update();
3182     if (save > 0) {
3183       sprintf(name, "%s.pdf", pad->GetName());
3184       pad->Print(name);
3185     } else if (save < 0) {
3186       sprintf(name, "%s.C", pad->GetName());
3187       pad->Print(name);
3188     }
3189   }
3190 }
3191 
3192 void PlotPropertyHist(const char* infile,
3193                       std::string prefix,
3194                       std::string text,
3195                       int etaMax = 25,
3196                       double lumi = 0,
3197                       double ener = 13.0,
3198                       bool dataMC = false,
3199                       bool drawStatBox = true,
3200                       int save = 0) {
3201   std::string name0[3] = {"energyE2", "energyH2", "energyP2"};
3202   std::string title0[3] = {"Energy in ECAL", "Energy in HCAL", "Track Momentum"};
3203   std::string xtitl0[3] = {"Energy (GeV)", "Energy (GeV)", "p (GeV)"};
3204   std::string name1[5] = {"eta02", "eta12", "eta22", "eta32", "eta42"};
3205   std::string name10[5] = {"eta0", "eta1", "eta2", "eta3", "eta4"};
3206   std::string xtitl1 = "i#eta";
3207   std::string name2[5] = {"p0", "p1", "p2", "p3", "p4"};
3208   std::string xtitl2 = "p (GeV)";
3209   std::string title1[5] = {"Tracks with p=40:60 GeV",
3210                            "Good Quality Tracks with p=40:60 GeV",
3211                            "Selected Tracks with p=40:60 GeV",
3212                            "Isolated Tracks with p=40:60 GeV",
3213                            "Isolated Tracks with p=40:60 GeV and MIPS in ECAL"};
3214   std::string title2[5] = {
3215       "All Tracks", "Good Quality Tracks", "Selected Tracks", "Isolated Tracks", "Isolated Tracks with MIPS in ECAL"};
3216   std::string ytitle = "Tracks";
3217 
3218   gStyle->SetCanvasBorderMode(0);
3219   gStyle->SetCanvasColor(kWhite);
3220   gStyle->SetPadColor(kWhite);
3221   gStyle->SetFillColor(kWhite);
3222   gStyle->SetOptTitle(0);
3223   if (drawStatBox)
3224     gStyle->SetOptStat(1110);
3225   else
3226     gStyle->SetOptStat(0);
3227   gStyle->SetOptFit(0);
3228 
3229   TFile* file = new TFile(infile);
3230   char name[100], namep[100];
3231   for (int k = 1; k <= etaMax; ++k) {
3232     for (int j = 0; j < 3; ++j) {
3233       sprintf(name, "%s%s%d", prefix.c_str(), name0[j].c_str(), k);
3234       TH1D* hist1 = (TH1D*)file->FindObjectAny(name);
3235       if (hist1 != nullptr) {
3236         TH1D* hist = (TH1D*)(hist1->Clone());
3237         double ymin(0.90);
3238         sprintf(namep, "c_%s", name);
3239         TCanvas* pad = new TCanvas(namep, namep, 700, 500);
3240         pad->SetLogy();
3241         pad->SetRightMargin(0.10);
3242         pad->SetTopMargin(0.10);
3243         hist->GetXaxis()->SetTitleSize(0.04);
3244         hist->GetXaxis()->SetTitle(xtitl0[j].c_str());
3245         hist->GetYaxis()->SetTitle(ytitle.c_str());
3246         hist->GetYaxis()->SetLabelOffset(0.005);
3247         hist->GetYaxis()->SetTitleSize(0.04);
3248         hist->GetYaxis()->SetLabelSize(0.035);
3249         hist->GetYaxis()->SetTitleOffset(1.10);
3250         hist->SetMarkerStyle(20);
3251         hist->Draw();
3252         pad->Update();
3253         TPaveStats* st1 = (TPaveStats*)hist->GetListOfFunctions()->FindObject("stats");
3254         if (st1 != nullptr) {
3255           st1->SetY1NDC(0.78);
3256           st1->SetY2NDC(0.90);
3257           st1->SetX1NDC(0.65);
3258           st1->SetX2NDC(0.90);
3259         }
3260         pad->Update();
3261         double ymx(0.96), xmi(0.35), xmx(0.95);
3262         char txt[100];
3263         if (lumi > 0.1) {
3264           ymx = ymin - 0.005;
3265           xmi = 0.45;
3266           TPaveText* txt0 = new TPaveText(0.65, 0.91, 0.90, 0.96, "blNDC");
3267           txt0->SetFillColor(0);
3268           sprintf(txt, "%4.1f TeV %5.1f fb^{-1}", ener, lumi);
3269           txt0->AddText(txt);
3270           txt0->Draw("same");
3271         }
3272         double ymi = ymx - 0.05;
3273         TPaveText* txt1 = new TPaveText(xmi, ymi, xmx, ymx, "blNDC");
3274         txt1->SetFillColor(0);
3275         if (text == "") {
3276           sprintf(txt, "%s", title0[j].c_str());
3277         } else {
3278           sprintf(txt, "%s (%s)", title0[j].c_str(), text.c_str());
3279         }
3280         txt1->AddText(txt);
3281         txt1->Draw("same");
3282         double xmax = (dataMC) ? 0.24 : 0.35;
3283         ymi = 0.91;
3284         ymx = ymi + 0.05;
3285         TPaveText* txt2 = new TPaveText(0.02, ymi, xmax, ymx, "blNDC");
3286         txt2->SetFillColor(0);
3287         if (dataMC)
3288           sprintf(txt, "CMS Preliminary");
3289         else
3290           sprintf(txt, "CMS Simulation Preliminary");
3291         txt2->AddText(txt);
3292         txt2->Draw("same");
3293         pad->Modified();
3294         pad->Update();
3295         if (save > 0) {
3296           sprintf(name, "%s.pdf", pad->GetName());
3297           pad->Print(name);
3298         } else if (save < 0) {
3299           sprintf(name, "%s.C", pad->GetName());
3300           pad->Print(name);
3301         }
3302       }
3303     }
3304   }
3305 
3306   for (int k = 0; k < 3; ++k) {
3307     for (int j = 0; j < 5; ++j) {
3308       if (k == 0)
3309         sprintf(name, "%s%s", prefix.c_str(), name1[j].c_str());
3310       else if (k == 1)
3311         sprintf(name, "%s%s", prefix.c_str(), name10[j].c_str());
3312       else
3313         sprintf(name, "%s%s", prefix.c_str(), name2[j].c_str());
3314       TH1D* hist1 = (TH1D*)file->FindObjectAny(name);
3315       if (hist1 != nullptr) {
3316         TH1D* hist = (TH1D*)(hist1->Clone());
3317         double ymin(0.90);
3318         sprintf(namep, "c_%s", name);
3319         TCanvas* pad = new TCanvas(namep, namep, 700, 500);
3320         pad->SetLogy();
3321         pad->SetRightMargin(0.10);
3322         pad->SetTopMargin(0.10);
3323         hist->GetXaxis()->SetTitleSize(0.04);
3324         if (k <= 1)
3325           hist->GetXaxis()->SetTitle(xtitl1.c_str());
3326         else
3327           hist->GetXaxis()->SetTitle(xtitl2.c_str());
3328         hist->GetYaxis()->SetTitle(ytitle.c_str());
3329         hist->GetYaxis()->SetLabelOffset(0.005);
3330         hist->GetYaxis()->SetTitleSize(0.04);
3331         hist->GetYaxis()->SetLabelSize(0.035);
3332         hist->GetYaxis()->SetTitleOffset(1.10);
3333         hist->SetMarkerStyle(20);
3334         hist->Draw();
3335         pad->Update();
3336         TPaveStats* st1 = (TPaveStats*)hist->GetListOfFunctions()->FindObject("stats");
3337         if (st1 != nullptr) {
3338           st1->SetY1NDC(0.78);
3339           st1->SetY2NDC(0.90);
3340           st1->SetX1NDC(0.65);
3341           st1->SetX2NDC(0.90);
3342         }
3343         pad->Update();
3344         double ymx(0.96), xmi(0.35), xmx(0.95);
3345         char txt[100];
3346         if (lumi > 0.1) {
3347           ymx = ymin - 0.005;
3348           xmi = 0.45;
3349           TPaveText* txt0 = new TPaveText(0.65, 0.91, 0.90, 0.96, "blNDC");
3350           txt0->SetFillColor(0);
3351           sprintf(txt, "%4.1f TeV %5.1f fb^{-1}", ener, lumi);
3352           txt0->AddText(txt);
3353           txt0->Draw("same");
3354         }
3355         double ymi = ymx - 0.05;
3356         TPaveText* txt1 = new TPaveText(xmi, ymi, xmx, ymx, "blNDC");
3357         txt1->SetFillColor(0);
3358         if (text == "") {
3359           if (k == 0)
3360             sprintf(txt, "%s", title1[j].c_str());
3361           else
3362             sprintf(txt, "%s", title2[j].c_str());
3363         } else {
3364           if (k == 0)
3365             sprintf(txt, "%s (%s)", title1[j].c_str(), text.c_str());
3366           else
3367             sprintf(txt, "%s (%s)", title2[j].c_str(), text.c_str());
3368         }
3369         txt1->AddText(txt);
3370         txt1->Draw("same");
3371         double xmax = (dataMC) ? 0.24 : 0.35;
3372         ymi = 0.91;
3373         ymx = ymi + 0.05;
3374         TPaveText* txt2 = new TPaveText(0.02, ymi, xmax, ymx, "blNDC");
3375         txt2->SetFillColor(0);
3376         if (dataMC)
3377           sprintf(txt, "CMS Preliminary");
3378         else
3379           sprintf(txt, "CMS Simulation Preliminary");
3380         txt2->AddText(txt);
3381         txt2->Draw("same");
3382         pad->Modified();
3383         pad->Update();
3384         if (save > 0) {
3385           sprintf(name, "%s.pdf", pad->GetName());
3386           pad->Print(name);
3387         } else if (save < 0) {
3388           sprintf(name, "%s.C", pad->GetName());
3389           pad->Print(name);
3390         }
3391       }
3392     }
3393   }
3394 }
3395 
3396 void PlotMeanError(const std::string infilest, int reg = 3, bool resol = false, int save = 0, bool debug = false) {
3397   bool ok(false);
3398   const int ntypmx = 3;
3399   const int nregmx = 4;
3400   if (reg < 0 || reg >= nregmx)
3401     reg = nregmx - 1;
3402   int nEner(0), nType(0), nPts(0);
3403   std::vector<double> energy[ntypmx], denergy[ntypmx], value[ntypmx], dvalue[ntypmx];
3404   // First read the data
3405   std::ifstream fInput(infilest.c_str());
3406   if (!fInput.good()) {
3407     std::cout << "Cannot open file " << infilest << std::endl;
3408   } else {
3409     ok = true;
3410     fInput >> nEner >> nType >> nPts;
3411     int nmax = nEner * nType;
3412     int type, elow, ehigh;
3413     double v1[4], e1[4], v2[4], e2[4];
3414     for (int n = 0; n < nmax; ++n) {
3415       fInput >> type >> elow >> ehigh;
3416       fInput >> v1[0] >> e1[0] >> v1[1] >> e1[1] >> v1[2] >> e1[2] >> v1[3] >> e1[3];
3417       fInput >> v2[0] >> e2[0] >> v2[1] >> e2[1] >> v2[2] >> e2[2] >> v2[3] >> e2[3];
3418       double ener = 0.5 * (ehigh + elow);
3419       double dene = 0.5 * (ehigh - elow);
3420       energy[type].push_back(ener);
3421       denergy[type].push_back(dene);
3422       if (resol) {
3423         value[type].push_back(v2[reg]);
3424         dvalue[type].push_back(e2[reg]);
3425       } else {
3426         value[type].push_back(v1[reg]);
3427         dvalue[type].push_back(e1[reg]);
3428       }
3429     }
3430     fInput.close();
3431     std::cout << "Reads " << (nmax + 1) << " cards from " << infilest << " with measurements for " << nEner
3432               << " energies and " << nType << " types" << std::endl;
3433     if (debug) {
3434       for (int n = 0; n < nType; ++n) {
3435         std::cout << "Type " << n << " with " << energy[n].size() << " points\n";
3436         for (unsigned int k = 0; k < energy[n].size(); ++k)
3437           std::cout << " [" << k << "] " << energy[n][k] << " +- " << denergy[n][k] << " Value " << value[n][k]
3438                     << " +- " << dvalue[n][k] << std::endl;
3439       }
3440     }
3441   }
3442 
3443   // Now the plots
3444   if (ok) {
3445     int mvsres = (resol) ? 1 : 0;
3446     std::string names[2] = {"Mean", "Resol"};
3447     std::string namet[nregmx] = {"Barrel", "Transition", "Endcap", "Combined"};
3448     char cname[100];
3449     sprintf(cname, "c_%s%s", names[mvsres].c_str(), namet[reg].c_str());
3450     int color[ntypmx] = {2, 4, 1};
3451     int mtype[ntypmx] = {20, 21, 22};
3452     double ymin[2] = {0.65, 0.10};
3453     double ymax[2] = {1.30, 0.50};
3454     gStyle->SetCanvasBorderMode(0);
3455     gStyle->SetCanvasColor(kWhite);
3456     gStyle->SetPadColor(kWhite);
3457     gStyle->SetFillColor(kWhite);
3458     gStyle->SetOptTitle(kFALSE);
3459     gStyle->SetPadBorderMode(0);
3460     gStyle->SetCanvasBorderMode(0);
3461     gStyle->SetOptStat(0);
3462     TCanvas* canvas = new TCanvas(cname, cname, 500, 500);
3463     canvas->SetTopMargin(0.05);
3464     canvas->SetBottomMargin(0.14);
3465     canvas->SetLeftMargin(0.15);
3466     canvas->SetRightMargin(0.10);
3467     TH1F* vFrame = canvas->DrawFrame(0.0, ymin[mvsres], 120.0, ymax[mvsres]);
3468     vFrame->GetXaxis()->SetRangeUser(0.0, 120.0);
3469     vFrame->GetYaxis()->SetRangeUser(ymin[mvsres], ymax[mvsres]);
3470     vFrame->GetXaxis()->SetLabelSize(0.04);
3471     vFrame->GetYaxis()->SetLabelSize(0.04);
3472     vFrame->GetXaxis()->SetTitleSize(0.04);
3473     vFrame->GetYaxis()->SetTitleSize(0.04);
3474     vFrame->GetXaxis()->SetTitleOffset(1.2);
3475     vFrame->GetYaxis()->SetTitleOffset(1.6);
3476     vFrame->GetXaxis()->SetLabelOffset(0.0);
3477     vFrame->GetXaxis()->SetTitle("p_{Beam} (GeV/c)");
3478     if (resol) {
3479       vFrame->GetYaxis()->SetTitle("Width(E_{HCAL}/(p-E_{ECAL}))");
3480     } else {
3481       vFrame->GetYaxis()->SetTitle("MPV(E_{HCAL}/(p-E_{ECAL}))");
3482     }
3483     TLegend* legend = new TLegend(0.70, 0.80, 0.90, 0.94);
3484     legend->SetFillColor(kWhite);
3485     std::string nameg[ntypmx] = {"MAHI", "M0", "M2"};
3486     for (int n = 0; n < nType; ++n) {
3487       unsigned int nmax0 = energy[n].size();
3488       double mom[nmax0], dmom[nmax0], mean[nmax0], dmean[nmax0];
3489       for (unsigned int k = 0; k < nmax0; ++k) {
3490         mom[k] = energy[n][k];
3491         dmom[k] = denergy[n][k];
3492         mean[k] = value[n][k];
3493         dmean[k] = dvalue[n][k];
3494       }
3495       TGraphErrors* graph = new TGraphErrors(nmax0, mom, mean, dmom, dmean);
3496       graph->SetMarkerStyle(mtype[n]);
3497       graph->SetMarkerColor(color[n]);
3498       graph->SetMarkerSize(1.4);
3499       graph->SetLineColor(color[n]);
3500       graph->SetLineWidth(2);
3501       sprintf(cname, "%s", nameg[n].c_str());
3502       legend->AddEntry(graph, cname, "lp");
3503       graph->Draw("P");
3504     }
3505     legend->Draw("same");
3506     std::string regions[nregmx] = {"20118B Barrel", "2018B Transition", "2018B Endcap", "2018B"};
3507     sprintf(cname, "%s", regions[reg].c_str());
3508     TPaveText* txt0 = new TPaveText(0.16, 0.90, 0.40, 0.94, "blNDC");
3509     txt0->SetFillColor(0);
3510     txt0->AddText(cname);
3511     txt0->Draw("same");
3512     TPaveText* txt1 = new TPaveText(0.15, 0.95, 0.40, 0.99, "blNDC");
3513     txt1->SetFillColor(0);
3514     sprintf(cname, "CMS Preliminary");
3515     txt1->AddText(cname);
3516     txt1->Draw("same");
3517     canvas->Update();
3518     if (save > 0) {
3519       sprintf(cname, "%s.pdf", canvas->GetName());
3520       canvas->Print(cname);
3521     } else if (save < 0) {
3522       sprintf(cname, "%s.C", canvas->GetName());
3523       canvas->Print(cname);
3524     }
3525   }
3526 }
[ Source navigation ] [ Diff markup ] [ Identifier search ] [ general search ]

This page was automatically generated by the 2.2.1 LXR engine.

The LXR team
Valid CSS 2.1! Valid HTML 4.01!