Back to home page

Project CMSSW displayed by LXR

 
 

    


File indexing completed on 2023-01-25 23:30:02

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 (debug)
0768             std::cout << "Histogram " << name << ":" << (hist->GetName()) << " with " << (hist->GetEntries())
0769                       << " entries" << std::endl;
0770           if (hist->GetEntries() > 0) {
0771             value = hist->GetMean();
0772             error = hist->GetRMS();
0773             for (int i = 1; i <= hist->GetNbinsX(); ++i)
0774               total += hist->GetBinContent(i);
0775             std::pair<double, double> rmserr = GetWidth(hist, 0.2, 2.0);
0776             width = rmserr.first;
0777             werror = rmserr.second;
0778           }
0779           if (total > 4) {
0780             if (nv1 > j)
0781               nv1 = j;
0782             if (nv2 < j)
0783               nv2 = j;
0784             if (j == 0) {
0785               sprintf(name, "%sOne", hist1->GetName());
0786               TH1D* hist2 = (TH1D*)hist1->Clone(name);
0787               fitOneGauss(hist2, true, debug);
0788               hists.push_back(hist2);
0789               results meaner = fitOneGauss(hist, true, debug);
0790               value = meaner.mean;
0791               error = meaner.errmean;
0792               width = meaner.width;
0793               werror = meaner.errwidth;
0794             } else {
0795               results meaner = fitOneGauss(hist, true, debug);
0796               value = meaner.mean;
0797               error = meaner.errmean;
0798               width = meaner.width;
0799               werror = meaner.errwidth;
0800             }
0801             if (j != 0) {
0802               if (j < jmin)
0803                 jmin = j;
0804               if (j > jmax)
0805                 jmax = j;
0806             }
0807           }
0808           hists.push_back(hist);
0809         }
0810         if (debug) {
0811           std::cout << "Hist " << j << " Value " << value << " +- " << error << std::endl;
0812         }
0813         if (j != 0) {
0814           double wbyv = width / value;
0815           double wverr = wbyv * std::sqrt((werror * werror) / (width * width) + (error * error) / (value * value));
0816           histo->SetBinContent(j, value);
0817           histo->SetBinError(j, error);
0818           histw->SetBinContent(j, wbyv);
0819           histw->SetBinError(j, wverr);
0820         }
0821       }
0822       if (histo != nullptr) {
0823         if (histo->GetEntries() > 2 && fiteta) {
0824           if (debug) {
0825             std::cout << "Jmin/max " << jmin << ":" << jmax << ":" << histo->GetNbinsX() << std::endl;
0826           }
0827           double LowEdge = histo->GetBinLowEdge(jmin);
0828           double HighEdge = histo->GetBinLowEdge(jmax) + histo->GetBinWidth(jmax);
0829           TFitResultPtr Fit = histo->Fit("pol0", "+QRWLS", "", LowEdge, HighEdge);
0830           if (debug) {
0831             std::cout << "Fit to Pol0: " << Fit->Value(0) << " +- " << Fit->FitResult::Error(0) << " in range " << nv1
0832                       << ":" << xbins[nv1] << ":" << nv2 << ":" << xbins[nv2] << std::endl;
0833           }
0834           histo->GetXaxis()->SetTitle("i#eta");
0835           histo->GetYaxis()->SetTitle("MPV(E_{HCAL}/(p-E_{ECAL}))");
0836           histo->GetYaxis()->SetRangeUser(0.4, 1.6);
0837           histw->GetXaxis()->SetTitle("i#eta");
0838           histw->GetYaxis()->SetTitle("MPV/Width(E_{HCAL}/(p-E_{ECAL}))");
0839           histw->GetYaxis()->SetRangeUser(0.0, 0.5);
0840         }
0841         hists.push_back(histo);
0842         hists.push_back(histw);
0843       } else {
0844         hists.push_back(hist0);
0845       }
0846 
0847       //Barrel,Endcap
0848       for (int j = 1; j <= 4; ++j) {
0849         sprintf(name, "%s%s%d%d", prefix.c_str(), ename.c_str(), iname, j);
0850         TH1D* hist1 = (TH1D*)file->FindObjectAny(name);
0851         if (debug) {
0852           std::cout << "Get Histogram for " << name << " at " << hist1 << std::endl;
0853         }
0854         if (hist1 != nullptr) {
0855           TH1D* hist = (TH1D*)hist1->Clone();
0856           double value(0), error(0), total(0), width(0), werror(0);
0857           if (hist->GetEntries() > 0) {
0858             value = hist->GetMean();
0859             error = hist->GetRMS();
0860             for (int i = 1; i <= hist->GetNbinsX(); ++i)
0861               total += hist->GetBinContent(i);
0862           }
0863           if (total > 4) {
0864             sprintf(name, "%sOne", hist1->GetName());
0865             TH1D* hist2 = (TH1D*)hist1->Clone(name);
0866             results meanerr = fitOneGauss(hist2, true, debug);
0867             value = meanerr.mean;
0868             error = meanerr.errmean;
0869             width = meanerr.width;
0870             werror = meanerr.errwidth;
0871             double wbyv = width / value;
0872             double wverr = wbyv * std::sqrt((werror * werror) / (width * width) + (error * error) / (value * value));
0873             std::cout << hist2->GetName() << " MPV " << value << " +- " << error << " Width " << width << " +- "
0874                       << werror << " W/M " << wbyv << " +- " << wverr << std::endl;
0875             hists.push_back(hist2);
0876             if (hist1->GetBinLowEdge(1) < 0.1) {
0877               sprintf(name, "%sTwo", hist1->GetName());
0878               TH1D* hist3 = (TH1D*)hist1->Clone(name);
0879               fitLanGau(hist3, debug);
0880               hists.push_back(hist3);
0881             }
0882             //            results meaner0 = fitTwoGauss(hist, debug);
0883             results meaner0 = fitOneGauss(hist, true, debug);
0884             value = meaner0.mean;
0885             error = meaner0.errmean;
0886             double rms;
0887             std::pair<double, double> meaner = GetMean(hist, 0.2, 2.0, rms);
0888             if (debug) {
0889               std::cout << "Fit " << value << ":" << error << ":" << hist->GetMeanError() << " Mean " << meaner.first
0890                         << ":" << meaner.second << std::endl;
0891             }
0892           }
0893           hists.push_back(hist);
0894         }
0895       }
0896     }
0897     TFile* theFile(0);
0898     if (append) {
0899       if (debug) {
0900         std::cout << "Open file " << outfile << " in append mode" << std::endl;
0901       }
0902       theFile = new TFile(outfile, "UPDATE");
0903     } else {
0904       if (debug) {
0905         std::cout << "Open file " << outfile << " in recreate mode" << std::endl;
0906       }
0907       theFile = new TFile(outfile, "RECREATE");
0908     }
0909 
0910     theFile->cd();
0911     for (unsigned int i = 0; i < hists.size(); ++i) {
0912       TH1D* hnew = (TH1D*)hists[i]->Clone();
0913       if (debug) {
0914         std::cout << "Write Histogram " << hnew->GetTitle() << std::endl;
0915       }
0916       hnew->Write();
0917     }
0918     theFile->Close();
0919     file->Close();
0920   }
0921 }
0922 
0923 void FitHistRBX(const char* infile, const char* outfile, std::string prefix, bool append = true, int iname = 3) {
0924   std::string sname("RBX"), lname("R");
0925   int numb(18);
0926   bool debug(false);
0927   char name[200];
0928 
0929   TFile* file = new TFile(infile);
0930   std::vector<TH1D*> hists;
0931   sprintf(name, "%s%s%d", prefix.c_str(), lname.c_str(), iname);
0932   TH1D* histo = new TH1D(name, "", numb, 0, numb);
0933   if (debug) {
0934     std::cout << infile << " " << file << std::endl;
0935   }
0936   if (file != nullptr) {
0937     for (int j = 0; j < numb; ++j) {
0938       sprintf(name, "%s%s%d%d", prefix.c_str(), sname.c_str(), iname, j + 1);
0939       TH1D* hist1 = (TH1D*)file->FindObjectAny(name);
0940       if (debug) {
0941         std::cout << "Get Histogram for " << name << " at " << hist1 << std::endl;
0942       }
0943       TH1D* hist = (TH1D*)hist1->Clone();
0944       double value(0), error(0), total(0);
0945       if (hist->GetEntries() > 0) {
0946         value = hist->GetMean();
0947         error = hist->GetRMS();
0948         for (int i = 1; i <= hist->GetNbinsX(); ++i)
0949           total += hist->GetBinContent(i);
0950       }
0951       if (total > 4) {
0952         results meaner = fitOneGauss(hist, false, debug);
0953         value = meaner.mean;
0954         error = meaner.errmean;
0955       }
0956       hists.push_back(hist);
0957       histo->SetBinContent(j + 1, value);
0958       histo->SetBinError(j + 1, error);
0959     }
0960     histo->GetXaxis()->SetTitle("RBX #");
0961     histo->GetYaxis()->SetTitle("MPV(E_{HCAL}/(p-E_{ECAL}))");
0962     histo->GetYaxis()->SetRangeUser(0.75, 1.20);
0963     hists.push_back(histo);
0964 
0965     TFile* theFile(0);
0966     if (append) {
0967       if (debug) {
0968         std::cout << "Open file " << outfile << " in append mode" << std::endl;
0969       }
0970       theFile = new TFile(outfile, "UPDATE");
0971     } else {
0972       if (debug) {
0973         std::cout << "Open file " << outfile << " in recreate mode" << std::endl;
0974       }
0975       theFile = new TFile(outfile, "RECREATE");
0976     }
0977 
0978     theFile->cd();
0979     for (unsigned int i = 0; i < hists.size(); ++i) {
0980       TH1D* hnew = (TH1D*)hists[i]->Clone();
0981       if (debug) {
0982         std::cout << "Write Histogram " << hnew->GetTitle() << std::endl;
0983       }
0984       hnew->Write();
0985     }
0986     theFile->Close();
0987     file->Close();
0988   }
0989 }
0990 
0991 void PlotHist(const char* infile,
0992               std::string prefix,
0993               std::string text,
0994               int mode = 4,
0995               int kopt = 100,
0996               double lumi = 0,
0997               double ener = 13.0,
0998               bool dataMC = false,
0999               bool drawStatBox = true,
1000               int save = 0) {
1001   std::string name0[6] = {"ratio00", "ratio10", "ratio20", "ratio30", "ratio40", "ratio50"};
1002   std::string name1[5] = {"Z0", "Z1", "Z2", "Z3", "Z4"};
1003   std::string name2[5] = {"L0", "L1", "L2", "L3", "L4"};
1004   std::string name3[5] = {"V0", "V1", "V2", "V3", "V4"};
1005   std::string name4[20] = {"etaB41", "etaB42", "etaB43", "etaB44", "etaB31", "etaB32", "etaB33",
1006                            "etaB34", "etaB21", "etaB22", "etaB23", "etaB24", "etaB11", "etaB12",
1007                            "etaB13", "etaB14", "etaB01", "etaB02", "etaB03", "etaB04"};
1008   std::string name5[5] = {"W0", "W1", "W2", "W3", "W4"};
1009   std::string title[6] = {"Tracks with p = 10:20 GeV",
1010                           "Tracks with p = 20:30 GeV",
1011                           "Tracks with p = 30:40 GeV",
1012                           "Tracks with p = 40:60 GeV",
1013                           "Tracks with p = 60:100 GeV",
1014                           "Tracks with p = 20:100 GeV"};
1015   std::string title1[20] = {"Tracks with p = 60:100 GeV (Barrel)", "Tracks with p = 60:100 GeV (Transition)",
1016                             "Tracks with p = 60:100 GeV (Endcap)", "Tracks with p = 60:100 GeV",
1017                             "Tracks with p = 40:60 GeV (Barrel)",  "Tracks with p = 40:60 GeV (Transition)",
1018                             "Tracks with p = 40:60 GeV (Endcap)",  "Tracks with p = 40:60 GeV",
1019                             "Tracks with p = 30:40 GeV (Barrel)",  "Tracks with p = 30:40 GeV (Transition)",
1020                             "Tracks with p = 30:40 GeV (Endcap)",  "Tracks with p = 30:40 GeV",
1021                             "Tracks with p = 20:30 GeV (Barrel)",  "Tracks with p = 20:30 GeV (Transition)",
1022                             "Tracks with p = 20:30 GeV (Endcap)",  "Tracks with p = 20:30 GeV",
1023                             "Tracks with p = 10:20 GeV (Barrel)",  "Tracks with p = 10:20 GeV (Transition)",
1024                             "Tracks with p = 10:20 GeV (Endcap)",  "Tracks with p = 10:20 GeV"};
1025   std::string xtitl[5] = {"E_{HCAL}/(p-E_{ECAL})", "i#eta", "d_{L1}", "# Vertex", "E_{HCAL}/(p-E_{ECAL})"};
1026   std::string ytitl[5] = {
1027       "Tracks", "MPV(E_{HCAL}/(p-E_{ECAL}))", "MPV(E_{HCAL}/(p-E_{ECAL}))", "MPV(E_{HCAL}/(p-E_{ECAL}))", "Tracks"};
1028 
1029   gStyle->SetCanvasBorderMode(0);
1030   gStyle->SetCanvasColor(kWhite);
1031   gStyle->SetPadColor(kWhite);
1032   gStyle->SetFillColor(kWhite);
1033   gStyle->SetOptTitle(0);
1034   if (mode < 0 || mode > 5)
1035     mode = 0;
1036   if (drawStatBox) {
1037     int iopt(1110);
1038     if (mode != 0)
1039       iopt = 10;
1040     gStyle->SetOptStat(iopt);
1041     gStyle->SetOptFit(1);
1042   } else {
1043     gStyle->SetOptStat(0);
1044     gStyle->SetOptFit(0);
1045   }
1046   TFile* file = new TFile(infile);
1047   TLine* line(0);
1048   char name[100], namep[100];
1049   int kmax = (mode == 4) ? 20 : (((mode < 1) && (mode > 5)) ? 6 : 5);
1050   for (int k = 0; k < kmax; ++k) {
1051     if (mode == 1) {
1052       sprintf(name, "%s%s", prefix.c_str(), name1[k].c_str());
1053     } else if (mode == 2) {
1054       sprintf(name, "%s%s", prefix.c_str(), name2[k].c_str());
1055     } else if (mode == 3) {
1056       sprintf(name, "%s%s", prefix.c_str(), name3[k].c_str());
1057     } else if (mode == 4) {
1058       if ((kopt / 100) % 10 == 0) {
1059         sprintf(name, "%s%s", prefix.c_str(), name4[k].c_str());
1060       } else if ((kopt / 100) % 10 == 2) {
1061         sprintf(name, "%s%sTwo", prefix.c_str(), name4[k].c_str());
1062       } else {
1063         sprintf(name, "%s%sOne", prefix.c_str(), name4[k].c_str());
1064       }
1065     } else if (mode == 5) {
1066       sprintf(name, "%s%s", prefix.c_str(), name5[k].c_str());
1067     } else {
1068       if ((kopt / 100) % 10 == 0) {
1069         sprintf(name, "%s%s", prefix.c_str(), name0[k].c_str());
1070       } else if ((kopt / 100) % 10 == 2) {
1071         sprintf(name, "%s%sTwo", prefix.c_str(), name0[k].c_str());
1072       } else {
1073         sprintf(name, "%s%sOne", prefix.c_str(), name0[k].c_str());
1074       }
1075     }
1076     TH1D* hist1 = (TH1D*)file->FindObjectAny(name);
1077     if (hist1 != nullptr) {
1078       TH1D* hist = (TH1D*)(hist1->Clone());
1079       double p0(1);
1080       double ymin(0.90);
1081       sprintf(namep, "c_%s", name);
1082       TCanvas* pad = new TCanvas(namep, namep, 700, 500);
1083       pad->SetRightMargin(0.10);
1084       pad->SetTopMargin(0.10);
1085       if ((kopt / 10) % 10 > 0)
1086         gPad->SetGrid();
1087       hist->GetXaxis()->SetTitleSize(0.04);
1088       hist->GetXaxis()->SetTitle(xtitl[mode].c_str());
1089       hist->GetYaxis()->SetTitle(ytitl[mode].c_str());
1090       hist->GetYaxis()->SetLabelOffset(0.005);
1091       hist->GetYaxis()->SetTitleSize(0.04);
1092       hist->GetYaxis()->SetLabelSize(0.035);
1093       hist->GetYaxis()->SetTitleOffset(1.10);
1094       if (mode == 0 || mode == 4) {
1095         if ((kopt / 100) % 10 == 2)
1096           hist->GetXaxis()->SetRangeUser(0.0, 0.30);
1097         else
1098           hist->GetXaxis()->SetRangeUser(0.25, 2.25);
1099       } else {
1100         if (mode == 5)
1101           hist->GetYaxis()->SetRangeUser(0.1, 0.50);
1102         else if (dataMC)
1103           hist->GetYaxis()->SetRangeUser(0.5, 1.50);
1104         else
1105           hist->GetYaxis()->SetRangeUser(0.8, 1.20);
1106         if (kopt % 10 > 0) {
1107           int nbin = hist->GetNbinsX();
1108           double LowEdge = (kopt % 10 == 1) ? hist->GetBinLowEdge(1) : -20;
1109           double HighEdge = (kopt % 10 == 1) ? hist->GetBinLowEdge(nbin) + hist->GetBinWidth(nbin) : 20;
1110           TFitResultPtr Fit = hist->Fit("pol0", "+QRWLS", "", LowEdge, HighEdge);
1111           p0 = Fit->Value(0);
1112         }
1113       }
1114       hist->SetMarkerStyle(20);
1115       hist->SetMarkerColor(2);
1116       hist->SetLineColor(2);
1117       hist->Draw();
1118       pad->Update();
1119       TPaveStats* st1 = (TPaveStats*)hist->GetListOfFunctions()->FindObject("stats");
1120       if (st1 != nullptr) {
1121         ymin = (mode == 0 || mode == 4) ? 0.70 : 0.80;
1122         st1->SetY1NDC(ymin);
1123         st1->SetY2NDC(0.90);
1124         st1->SetX1NDC(0.65);
1125         st1->SetX2NDC(0.90);
1126       }
1127       if (mode != 0 && mode != 4 && kopt % 10 > 0) {
1128         double xmin = hist->GetBinLowEdge(1);
1129         int nbin = hist->GetNbinsX();
1130         double xmax = hist->GetBinLowEdge(nbin) + hist->GetBinWidth(nbin);
1131         double mean(0), rms(0), total(0);
1132         int kount(0);
1133         for (int k = 2; k < nbin; ++k) {
1134           double x = hist->GetBinContent(k);
1135           double w = hist->GetBinError(k);
1136           mean += (x * w);
1137           rms += (x * x * w);
1138           total += w;
1139           ++kount;
1140         }
1141         mean /= total;
1142         rms /= total;
1143         double error = sqrt(rms - mean * mean) / total;
1144         line = new TLine(xmin, p0, xmax, p0);  //etamin,1.0,etamax,1.0);
1145         std::cout << xmin << ":" << xmax << ":" << p0 << " Mean " << nbin << ":" << kount << ":" << total << ":" << mean
1146                   << ":" << rms << ":" << error << std::endl;
1147         line->SetLineWidth(2);
1148         line->SetLineStyle(2);
1149         line->Draw("same");
1150       }
1151       pad->Update();
1152       double ymx(0.96), xmi(0.25), xmx(0.90);
1153       char txt[100];
1154       if (lumi > 0.1) {
1155         ymx = ymin - 0.005;
1156         xmi = 0.45;
1157         TPaveText* txt0 = new TPaveText(0.65, 0.91, 0.90, 0.96, "blNDC");
1158         txt0->SetFillColor(0);
1159         sprintf(txt, "%4.1f TeV %5.1f fb^{-1}", ener, lumi);
1160         txt0->AddText(txt);
1161         txt0->Draw("same");
1162       }
1163       double ymi = ymx - 0.05;
1164       TPaveText* txt1 = new TPaveText(xmi, ymi, xmx, ymx, "blNDC");
1165       txt1->SetFillColor(0);
1166       if (text == "") {
1167         if (mode == 4)
1168           sprintf(txt, "%s", title1[k].c_str());
1169         else
1170           sprintf(txt, "%s", title[k].c_str());
1171       } else {
1172         if (mode == 4)
1173           sprintf(txt, "%s (%s)", title1[k].c_str(), text.c_str());
1174         else
1175           sprintf(txt, "%s (%s)", title[k].c_str(), text.c_str());
1176       }
1177       txt1->AddText(txt);
1178       txt1->Draw("same");
1179       double xmax = (dataMC) ? 0.33 : 0.44;
1180       ymi = (lumi > 0.1) ? 0.91 : 0.84;
1181       ymx = ymi + 0.05;
1182       TPaveText* txt2 = new TPaveText(0.11, ymi, xmax, ymx, "blNDC");
1183       txt2->SetFillColor(0);
1184       if (dataMC)
1185         sprintf(txt, "CMS Preliminary");
1186       else
1187         sprintf(txt, "CMS Simulation Preliminary");
1188       txt2->AddText(txt);
1189       txt2->Draw("same");
1190       pad->Modified();
1191       pad->Update();
1192       if (save > 0) {
1193         sprintf(name, "%s.pdf", pad->GetName());
1194         pad->Print(name);
1195       } else if (save < 0) {
1196         sprintf(name, "%s.C", pad->GetName());
1197         pad->Print(name);
1198       }
1199     }
1200   }
1201 }
1202 
1203 void PlotHistEta(const char* infile,
1204                  std::string prefix,
1205                  std::string text,
1206                  int iene = 3,
1207                  int numb = 50,
1208                  int ieta = 0,
1209                  double lumi = 0,
1210                  double ener = 13.0,
1211                  bool dataMC = false,
1212                  bool drawStatBox = true,
1213                  int save = 0) {
1214   std::string name0 = "ratio";
1215   std::string title[5] = {"10:20", "20:30", "30:40", "40:60", "60:100"};
1216   std::string xtitl = "E_{HCAL}/(p-E_{ECAL})";
1217   std::string ytitl = "Tracks";
1218 
1219   gStyle->SetCanvasBorderMode(0);
1220   gStyle->SetCanvasColor(kWhite);
1221   gStyle->SetPadColor(kWhite);
1222   gStyle->SetFillColor(kWhite);
1223   gStyle->SetOptTitle(0);
1224   if (drawStatBox) {
1225     int iopt(1110);
1226     gStyle->SetOptStat(iopt);
1227     gStyle->SetOptFit(1);
1228   } else {
1229     gStyle->SetOptStat(0);
1230     gStyle->SetOptFit(0);
1231   }
1232   if (iene < 0 || iene >= 5)
1233     iene = 3;
1234   int numb2 = numb / 2;
1235   if (ieta < -numb2 || ieta > numb2)
1236     ieta = 0;
1237   int ietaMin = ((ieta == 0) ? 1 : ((ieta > 0) ? (numb2 + ieta) : (numb2 + ieta + 1)));
1238   int ietaMax = (ieta == 0) ? numb : ietaMin;
1239   TFile* file = new TFile(infile);
1240   char name[100], namep[100];
1241   for (int k = ietaMin; k <= ietaMax; ++k) {
1242     int eta = (k > numb2) ? (k - numb2) : (k - numb2 - 1);
1243     sprintf(name, "%s%s%d%d", prefix.c_str(), name0.c_str(), iene, k);
1244     TH1D* hist1 = (TH1D*)file->FindObjectAny(name);
1245     std::cout << name << " at " << hist1 << std::endl;
1246     if (hist1 != nullptr) {
1247       TH1D* hist = (TH1D*)(hist1->Clone());
1248       double ymin(0.90);
1249       sprintf(namep, "c_%s", name);
1250       TCanvas* pad = new TCanvas(namep, namep, 700, 500);
1251       pad->SetRightMargin(0.10);
1252       pad->SetTopMargin(0.10);
1253       hist->GetXaxis()->SetTitleSize(0.04);
1254       hist->GetXaxis()->SetTitle(xtitl.c_str());
1255       hist->GetYaxis()->SetTitle(ytitl.c_str());
1256       hist->GetYaxis()->SetLabelOffset(0.005);
1257       hist->GetYaxis()->SetTitleSize(0.04);
1258       hist->GetYaxis()->SetLabelSize(0.035);
1259       hist->GetYaxis()->SetTitleOffset(1.10);
1260       hist->GetXaxis()->SetRangeUser(0.25, 2.25);
1261       hist->SetMarkerStyle(20);
1262       hist->SetMarkerColor(2);
1263       hist->SetLineColor(2);
1264       hist->Draw();
1265       pad->Update();
1266       TPaveStats* st1 = (TPaveStats*)hist->GetListOfFunctions()->FindObject("stats");
1267       if (st1 != nullptr) {
1268         ymin = 0.70;
1269         st1->SetY1NDC(ymin);
1270         st1->SetY2NDC(0.90);
1271         st1->SetX1NDC(0.65);
1272         st1->SetX2NDC(0.90);
1273       }
1274       double ymx(0.96), xmi(0.25), xmx(0.90);
1275       char txt[100];
1276       if (lumi > 0.1) {
1277         ymx = ymin - 0.005;
1278         xmi = 0.45;
1279         TPaveText* txt0 = new TPaveText(0.65, 0.91, 0.90, 0.96, "blNDC");
1280         txt0->SetFillColor(0);
1281         sprintf(txt, "%4.1f TeV %5.1f fb^{-1}", ener, lumi);
1282         txt0->AddText(txt);
1283         txt0->Draw("same");
1284       }
1285       double ymi = ymx - 0.05;
1286       TPaveText* txt1 = new TPaveText(xmi, ymi, xmx, ymx, "blNDC");
1287       txt1->SetFillColor(0);
1288       if (text == "") {
1289         sprintf(txt, "Tracks with p = %s GeV at i#eta = %d", title[iene].c_str(), eta);
1290       } else {
1291         sprintf(txt, "Tracks with p = %s GeV at i#eta = %d (%s)", title[iene].c_str(), eta, text.c_str());
1292       }
1293       txt1->AddText(txt);
1294       txt1->Draw("same");
1295       double xmax = (dataMC) ? 0.33 : 0.44;
1296       ymi = (lumi > 0.1) ? 0.91 : 0.84;
1297       ymx = ymi + 0.05;
1298       TPaveText* txt2 = new TPaveText(0.11, ymi, xmax, ymx, "blNDC");
1299       txt2->SetFillColor(0);
1300       if (dataMC)
1301         sprintf(txt, "CMS Preliminary");
1302       else
1303         sprintf(txt, "CMS Simulation Preliminary");
1304       txt2->AddText(txt);
1305       txt2->Draw("same");
1306       pad->Modified();
1307       pad->Update();
1308       if (save > 0) {
1309         sprintf(name, "%s.pdf", pad->GetName());
1310         pad->Print(name);
1311       } else if (save < 0) {
1312         sprintf(name, "%s.C", pad->GetName());
1313         pad->Print(name);
1314       }
1315     }
1316   }
1317 }
1318 
1319 void PlotHists(std::string infile, std::string prefix, std::string text, bool drawStatBox = true, int save = 0) {
1320   int colors[6] = {1, 6, 4, 7, 2, 9};
1321   std::string types[6] = {"B", "C", "D", "E", "F", "G"};
1322   std::string names[3] = {"ratio20", "Z2", "W2"};
1323   std::string xtitl[3] = {"E_{HCAL}/(p-E_{ECAL})", "i#eta", "i#eta"};
1324   std::string ytitl[3] = {"Tracks", "MPV(E_{HCAL}/(p-E_{ECAL}))", "MPV/Width(E_{HCAL}/(p-E_{ECAL}))"};
1325 
1326   gStyle->SetCanvasBorderMode(0);
1327   gStyle->SetCanvasColor(kWhite);
1328   gStyle->SetPadColor(kWhite);
1329   gStyle->SetFillColor(kWhite);
1330   gStyle->SetOptTitle(0);
1331   if (drawStatBox)
1332     gStyle->SetOptFit(10);
1333   else
1334     gStyle->SetOptFit(0);
1335 
1336   char name[100], namep[100];
1337   TFile* file = new TFile(infile.c_str());
1338   for (int i = 0; i < 3; ++i) {
1339     std::vector<TH1D*> hists;
1340     std::vector<int> kks;
1341     double ymax(0.77);
1342     if (drawStatBox) {
1343       if (i == 0)
1344         gStyle->SetOptStat(1100);
1345       else
1346         gStyle->SetOptStat(10);
1347     } else {
1348       gStyle->SetOptStat(0);
1349       ymax = 0.89;
1350     }
1351     for (int k = 0; k < 6; ++k) {
1352       sprintf(name, "%s%s%s", prefix.c_str(), types[k].c_str(), names[i].c_str());
1353       TH1D* hist1 = (TH1D*)file->FindObjectAny(name);
1354       if (hist1 != nullptr) {
1355         hists.push_back((TH1D*)(hist1->Clone()));
1356         kks.push_back(k);
1357       }
1358     }
1359     if (hists.size() > 0) {
1360       sprintf(namep, "c_%s%s", prefix.c_str(), names[i].c_str());
1361       TCanvas* pad = new TCanvas(namep, namep, 700, 500);
1362       TLegend* legend = new TLegend(0.44, 0.89 - 0.055 * hists.size(), 0.69, ymax);
1363       legend->SetFillColor(kWhite);
1364       pad->SetRightMargin(0.10);
1365       pad->SetTopMargin(0.10);
1366       double ymax(0.90);
1367       double dy = (i == 0) ? 0.13 : 0.08;
1368       for (unsigned int jk = 0; jk < hists.size(); ++jk) {
1369         int k = kks[jk];
1370         hists[jk]->GetXaxis()->SetTitle(xtitl[i].c_str());
1371         hists[jk]->GetYaxis()->SetTitle(ytitl[i].c_str());
1372         hists[jk]->GetYaxis()->SetLabelOffset(0.005);
1373         hists[jk]->GetYaxis()->SetLabelSize(0.035);
1374         hists[jk]->GetYaxis()->SetTitleOffset(1.15);
1375         if (i == 0) {
1376           hists[jk]->GetXaxis()->SetRangeUser(0.0, 2.5);
1377         } else if (i == 1) {
1378           hists[jk]->GetYaxis()->SetRangeUser(0.5, 2.0);
1379         } else {
1380           hists[jk]->GetYaxis()->SetRangeUser(0.0, 0.5);
1381         }
1382         hists[jk]->SetMarkerStyle(20);
1383         hists[jk]->SetMarkerColor(colors[k]);
1384         hists[jk]->SetLineColor(colors[k]);
1385         if (jk == 0)
1386           hists[jk]->Draw();
1387         else
1388           hists[jk]->Draw("sames");
1389         pad->Update();
1390         TPaveStats* st1 = (TPaveStats*)hists[jk]->GetListOfFunctions()->FindObject("stats");
1391         if (st1 != nullptr) {
1392           double ymin = ymax - dy;
1393           st1->SetLineColor(colors[k]);
1394           st1->SetTextColor(colors[k]);
1395           st1->SetY1NDC(ymin);
1396           st1->SetY2NDC(ymax);
1397           st1->SetX1NDC(0.70);
1398           st1->SetX2NDC(0.90);
1399           ymax = ymin;
1400         }
1401         sprintf(name, "%s%s", text.c_str(), types[k].c_str());
1402         legend->AddEntry(hists[jk], name, "lp");
1403       }
1404       legend->Draw("same");
1405       pad->Update();
1406       TPaveText* txt1 = new TPaveText(0.34, 0.825, 0.69, 0.895, "blNDC");
1407       txt1->SetFillColor(0);
1408       char txt[100];
1409       sprintf(txt, "Tracks with p = 40:60 GeV");
1410       txt1->AddText(txt);
1411       txt1->Draw("same");
1412       TPaveText* txt2 = new TPaveText(0.11, 0.825, 0.33, 0.895, "blNDC");
1413       txt2->SetFillColor(0);
1414       sprintf(txt, "CMS Preliminary");
1415       txt2->AddText(txt);
1416       txt2->Draw("same");
1417       if (!drawStatBox && i == 1) {
1418         double xmin = hists[0]->GetBinLowEdge(1);
1419         int nbin = hists[0]->GetNbinsX();
1420         double xmax = hists[0]->GetBinLowEdge(nbin) + hists[0]->GetBinWidth(nbin);
1421         TLine line = TLine(xmin, 1.0, xmax, 1.0);  //etamin,1.0,etamax,1.0);
1422         line.SetLineWidth(4);
1423         line.Draw("same");
1424         pad->Update();
1425       }
1426       pad->Modified();
1427       pad->Update();
1428       if (save > 0) {
1429         sprintf(name, "%s.pdf", pad->GetName());
1430         pad->Print(name);
1431       } else if (save < 0) {
1432         sprintf(name, "%s.C", pad->GetName());
1433         pad->Print(name);
1434       }
1435     }
1436   }
1437 }
1438 
1439 void PlotTwoHists(std::string infile,
1440                   std::string prefix1,
1441                   std::string text1,
1442                   std::string prefix2,
1443                   std::string text2,
1444                   std::string text0,
1445                   int type = 0,
1446                   int iname = 3,
1447                   double lumi = 0,
1448                   double ener = 13.0,
1449                   int drawStatBox = 0,
1450                   int save = 0) {
1451   int colors[2] = {2, 4};
1452   int numb[2] = {5, 1};
1453   std::string names0[5] = {"ratio00", "ratio00One", "etaB04One", "Z0", "W0"};
1454   std::string names1[5] = {"ratio10", "ratio10One", "etaB14One", "Z1", "W1"};
1455   std::string names2[5] = {"ratio30", "ratio30One", "etaB34One", "Z3", "W3"};
1456   std::string xtitl1[5] = {"E_{HCAL}/(p-E_{ECAL})", "E_{HCAL}/(p-E_{ECAL})", "E_{HCAL}/(p-E_{ECAL})", "i#eta", "i#eta"};
1457   std::string ytitl1[5] = {
1458       "Tracks", "Tracks", "Tracks", "MPV(E_{HCAL}/(p-E_{ECAL}))", "MPV/Width(E_{HCAL}/(p-E_{ECAL}))"};
1459   std::string names3[1] = {"R"};
1460   std::string xtitl2[1] = {"RBX #"};
1461   std::string ytitl2[1] = {"MPV(E_{HCAL}/(p-E_{ECAL}))"};
1462 
1463   gStyle->SetCanvasBorderMode(0);
1464   gStyle->SetCanvasColor(kWhite);
1465   gStyle->SetPadColor(kWhite);
1466   gStyle->SetFillColor(kWhite);
1467   gStyle->SetOptTitle(0);
1468   if ((drawStatBox / 10) % 10 > 0)
1469     gStyle->SetOptFit(10);
1470   else
1471     gStyle->SetOptFit(0);
1472 
1473   if (type != 1)
1474     type = 0;
1475   char name[100], namep[100];
1476   TFile* file = new TFile(infile.c_str());
1477   for (int i = 0; i < numb[type]; ++i) {
1478     std::vector<TH1D*> hists;
1479     std::vector<int> kks;
1480     double ymax(0.77);
1481     if (drawStatBox % 10 > 0) {
1482       if (i != 2)
1483         gStyle->SetOptStat(1100);
1484       else
1485         gStyle->SetOptStat(10);
1486     } else {
1487       gStyle->SetOptStat(0);
1488       ymax = 0.82;
1489     }
1490     for (int k = 0; k < 2; ++k) {
1491       std::string prefix = (k == 0) ? prefix1 : prefix2;
1492       if (type == 0) {
1493         if (iname == 0)
1494           sprintf(name, "%s%s", prefix.c_str(), names0[i].c_str());
1495         else if (iname == 1)
1496           sprintf(name, "%s%s", prefix.c_str(), names1[i].c_str());
1497         else
1498           sprintf(name, "%s%s", prefix.c_str(), names2[i].c_str());
1499       } else {
1500         sprintf(name, "%s%s%d", prefix.c_str(), names3[i].c_str(), iname);
1501       }
1502       TH1D* hist1 = (TH1D*)file->FindObjectAny(name);
1503       if (hist1 != nullptr) {
1504         hists.push_back((TH1D*)(hist1->Clone()));
1505         kks.push_back(k);
1506       }
1507     }
1508     if (hists.size() == 2) {
1509       if (type == 0) {
1510         if (iname == 0)
1511           sprintf(namep, "c_%s%s%s", prefix1.c_str(), prefix2.c_str(), names0[i].c_str());
1512         else if (iname == 1)
1513           sprintf(namep, "c_%s%s%s", prefix1.c_str(), prefix2.c_str(), names1[i].c_str());
1514         else
1515           sprintf(namep, "c_%s%s%s", prefix1.c_str(), prefix2.c_str(), names2[i].c_str());
1516       } else {
1517         sprintf(namep, "c_%s%s%s%d", prefix1.c_str(), prefix2.c_str(), names3[i].c_str(), iname);
1518       }
1519       double ymax(0.90);
1520       double dy = (i == 0) ? 0.13 : 0.08;
1521       double ymx0 = (drawStatBox == 0) ? (ymax - .01) : (ymax - dy * hists.size() - .01);
1522       TCanvas* pad = new TCanvas(namep, namep, 700, 500);
1523       TLegend* legend = new TLegend(0.64, ymx0 - 0.05 * hists.size(), 0.89, ymx0);
1524       legend->SetFillColor(kWhite);
1525       pad->SetRightMargin(0.10);
1526       pad->SetTopMargin(0.10);
1527       for (unsigned int jk = 0; jk < hists.size(); ++jk) {
1528         int k = kks[jk];
1529         hists[jk]->GetXaxis()->SetTitleSize(0.040);
1530         if (type == 0) {
1531           hists[jk]->GetXaxis()->SetTitle(xtitl1[i].c_str());
1532           hists[jk]->GetYaxis()->SetTitle(ytitl1[i].c_str());
1533         } else {
1534           hists[jk]->GetXaxis()->SetTitle(xtitl2[i].c_str());
1535           hists[jk]->GetYaxis()->SetTitle(ytitl2[i].c_str());
1536         }
1537         hists[jk]->GetYaxis()->SetLabelOffset(0.005);
1538         hists[jk]->GetYaxis()->SetLabelSize(0.035);
1539         hists[jk]->GetYaxis()->SetTitleSize(0.040);
1540         hists[jk]->GetYaxis()->SetTitleOffset(1.15);
1541         if ((type == 0) && (i != 3) && (i != 4))
1542           hists[jk]->GetXaxis()->SetRangeUser(0.0, 5.0);
1543         if (type == 0) {
1544           if (i == 3)
1545             hists[jk]->GetYaxis()->SetRangeUser(0.8, 1.2);
1546           else if (i == 4)
1547             hists[jk]->GetYaxis()->SetRangeUser(0.0, 0.5);
1548         }
1549         if (type != 0)
1550           hists[jk]->GetYaxis()->SetRangeUser(0.75, 1.2);
1551         hists[jk]->SetMarkerStyle(20);
1552         hists[jk]->SetMarkerColor(colors[k]);
1553         hists[jk]->SetLineColor(colors[k]);
1554         if (jk == 0)
1555           hists[jk]->Draw();
1556         else
1557           hists[jk]->Draw("sames");
1558         pad->Update();
1559         TPaveStats* st1 = (TPaveStats*)hists[jk]->GetListOfFunctions()->FindObject("stats");
1560         if (st1 != nullptr) {
1561           double ymin = ymax - dy;
1562           st1->SetLineColor(colors[k]);
1563           st1->SetTextColor(colors[k]);
1564           st1->SetY1NDC(ymin);
1565           st1->SetY2NDC(ymax);
1566           st1->SetX1NDC(0.70);
1567           st1->SetX2NDC(0.90);
1568           ymax = ymin;
1569         }
1570         if (k == 0)
1571           sprintf(name, "%s", text1.c_str());
1572         else
1573           sprintf(name, "%s", text2.c_str());
1574         legend->AddEntry(hists[jk], name, "lp");
1575       }
1576       legend->Draw("same");
1577       pad->Update();
1578       char txt[100];
1579       double xmi(0.10), xmx(0.895), ymx(0.95);
1580       if (lumi > 0.01) {
1581         xmx = 0.70;
1582         xmi = 0.30;
1583         TPaveText* txt0 = new TPaveText(0.705, 0.905, 0.90, 0.95, "blNDC");
1584         txt0->SetFillColor(0);
1585         sprintf(txt, "%4.1f TeV %5.1f fb^{-1}", ener, lumi);
1586         txt0->AddText(txt);
1587         txt0->Draw("same");
1588       }
1589       double ymi = ymx - 0.045;
1590       TPaveText* txt1 = new TPaveText(xmi, ymi, xmx, ymx, "blNDC");
1591       txt1->SetFillColor(0);
1592       if (iname == 0)
1593         sprintf(txt, "p = 20:30 GeV %s", text0.c_str());
1594       else
1595         sprintf(txt, "p = 40:60 GeV %s", text0.c_str());
1596       txt1->AddText(txt);
1597       txt1->Draw("same");
1598       ymi = (lumi > 0.1) ? 0.905 : 0.85;
1599       ymx = ymi + 0.045;
1600       TPaveText* txt2 = new TPaveText(0.105, ymi, 0.295, ymx, "blNDC");
1601       txt2->SetFillColor(0);
1602       sprintf(txt, "CMS Preliminary");
1603       txt2->AddText(txt);
1604       txt2->Draw("same");
1605       pad->Modified();
1606       if ((drawStatBox == 0) && (i == 3)) {
1607         double xmin = hists[0]->GetBinLowEdge(1);
1608         int nbin = hists[0]->GetNbinsX();
1609         double xmax = hists[0]->GetBinLowEdge(nbin) + hists[0]->GetBinWidth(nbin);
1610         TLine* line = new TLine(xmin, 1.0, xmax, 1.0);  //etamin,1.0,etamax,1.0);
1611         line->SetLineWidth(2);
1612         line->Draw("same");
1613         pad->Update();
1614       }
1615       pad->Update();
1616       if (save > 0) {
1617         sprintf(name, "%s.pdf", pad->GetName());
1618         pad->Print(name);
1619       } else if (save < 0) {
1620         sprintf(name, "%s.C", pad->GetName());
1621         pad->Print(name);
1622       }
1623     }
1624   }
1625 }
1626 
1627 void PlotFiveHists(std::string infile,
1628                    std::string text0,
1629                    std::string prefix0,
1630                    int type = 0,
1631                    int iname = 3,
1632                    int drawStatBox = 0,
1633                    bool normalize = false,
1634                    int save = 0,
1635                    std::string prefix1 = "",
1636                    std::string text1 = "",
1637                    std::string prefix2 = "",
1638                    std::string text2 = "",
1639                    std::string prefix3 = "",
1640                    std::string text3 = "",
1641                    std::string prefix4 = "",
1642                    std::string text4 = "",
1643                    std::string prefix5 = "",
1644                    std::string text5 = "") {
1645   int colors[5] = {2, 4, 6, 1, 7};
1646   int numb[3] = {5, 1, 4};
1647   std::string names0[5] = {"ratio00", "ratio00One", "etaB04", "Z0", "W0"};
1648   std::string names1[5] = {"ratio10", "ratio10One", "etaB14", "Z1", "W1"};
1649   std::string names2[5] = {"ratio30", "ratio30One", "etaB34", "Z3", "W3"};
1650   std::string xtitl1[5] = {"E_{HCAL}/(p-E_{ECAL})", "E_{HCAL}/(p-E_{ECAL})", "E_{HCAL}/(p-E_{ECAL})", "i#eta", "i#eta"};
1651   std::string ytitl1[5] = {
1652       "Tracks", "Tracks", "Tracks", "MPV(E_{HCAL}/(p-E_{ECAL}))", "MPV/Width(E_{HCAL}/(p-E_{ECAL}))"};
1653   std::string names3[1] = {"R"};
1654   std::string xtitl2[1] = {"RBX #"};
1655   std::string ytitl2[1] = {"MPV(E_{HCAL}/(p-E_{ECAL}))"};
1656   std::string names4[4] = {"pp21", "pp22", "pp23", "pp24"};
1657   std::string xtitl3[4] = {"p (GeV)", "p (GeV)", "p (GeV)", "p (GeV)"};
1658   std::string ytitl3[4] = {"Tracks", "Tracks", "Tracks", "Tracks"};
1659   std::string title3[4] = {"Barrel", "Transition", "Endcap", "Combined"};
1660 
1661   gStyle->SetCanvasBorderMode(0);
1662   gStyle->SetCanvasColor(kWhite);
1663   gStyle->SetPadColor(kWhite);
1664   gStyle->SetFillColor(kWhite);
1665   gStyle->SetOptTitle(0);
1666   if ((drawStatBox / 10) % 10 > 0)
1667     gStyle->SetOptFit(10);
1668   else
1669     gStyle->SetOptFit(0);
1670 
1671   if (type != 1 && type != 2)
1672     type = 0;
1673   char name[100], namep[100];
1674   TFile* file = new TFile(infile.c_str());
1675   for (int i = 0; i < numb[type]; ++i) {
1676     std::vector<TH1D*> hists;
1677     std::vector<int> kks;
1678     std::vector<std::string> texts;
1679     double ymax(0.77);
1680     if (drawStatBox % 10 > 0) {
1681       if (type == 2)
1682         gStyle->SetOptStat(1110);
1683       else if (i != 3)
1684         gStyle->SetOptStat(1100);
1685       else
1686         gStyle->SetOptStat(10);
1687     } else {
1688       gStyle->SetOptStat(0);
1689       ymax = 0.82;
1690     }
1691     for (int k = 0; k < 5; ++k) {
1692       std::string prefix, text;
1693       if (k == 0) {
1694         prefix = prefix1;
1695         text = text1;
1696       } else if (k == 1) {
1697         prefix = prefix2;
1698         text = text2;
1699       } else if (k == 2) {
1700         prefix = prefix3;
1701         text = text3;
1702       } else if (k == 3) {
1703         prefix = prefix4;
1704         text = text4;
1705       } else {
1706         prefix = prefix5;
1707         text = text5;
1708       }
1709       if (prefix != "") {
1710         if (type == 0) {
1711           if (iname == 0)
1712             sprintf(name, "%s%s", prefix.c_str(), names0[i].c_str());
1713           else if (iname == 1)
1714             sprintf(name, "%s%s", prefix.c_str(), names1[i].c_str());
1715           else
1716             sprintf(name, "%s%s", prefix.c_str(), names2[i].c_str());
1717         } else if (type == 1) {
1718           sprintf(name, "%s%s%d", prefix.c_str(), names3[i].c_str(), iname);
1719         } else {
1720           sprintf(name, "%s%s", prefix.c_str(), names4[i].c_str());
1721         }
1722         TH1D* hist1 = (TH1D*)file->FindObjectAny(name);
1723         if (hist1 != nullptr) {
1724           hists.push_back((TH1D*)(hist1->Clone()));
1725           kks.push_back(k);
1726           texts.push_back(text);
1727         }
1728       }
1729     }
1730     if (hists.size() > 0) {
1731       if (type == 0) {
1732         if (iname == 0)
1733           sprintf(namep, "c_%s%s", prefix0.c_str(), names0[i].c_str());
1734         else if (iname == 1)
1735           sprintf(namep, "c_%s%s", prefix0.c_str(), names1[i].c_str());
1736         else
1737           sprintf(namep, "c_%s%s", prefix0.c_str(), names2[i].c_str());
1738       } else if (type == 1) {
1739         sprintf(namep, "c_%s%s%d", prefix0.c_str(), names3[i].c_str(), iname);
1740       } else {
1741         sprintf(namep, "c_%s%s", prefix0.c_str(), names4[i].c_str());
1742       }
1743       double ymax(0.90);
1744       double dy = (i == 0 && type == 0) ? 0.13 : 0.08;
1745       double ymx0 = (drawStatBox == 0) ? (ymax - .01) : (ymax - dy * hists.size() - .01);
1746       TCanvas* pad = new TCanvas(namep, namep, 700, 500);
1747       TLegend* legend = new TLegend(0.64, ymx0 - 0.05 * hists.size(), 0.89, ymx0);
1748       legend->SetFillColor(kWhite);
1749       pad->SetRightMargin(0.10);
1750       pad->SetTopMargin(0.10);
1751       for (unsigned int jk = 0; jk < hists.size(); ++jk) {
1752         int k = kks[jk];
1753         hists[jk]->GetXaxis()->SetTitleSize(0.040);
1754         if (type == 0) {
1755           hists[jk]->GetXaxis()->SetTitle(xtitl1[i].c_str());
1756           hists[jk]->GetYaxis()->SetTitle(ytitl1[i].c_str());
1757         } else if (type == 1) {
1758           hists[jk]->GetXaxis()->SetTitle(xtitl2[i].c_str());
1759           hists[jk]->GetYaxis()->SetTitle(ytitl2[i].c_str());
1760         } else {
1761           hists[jk]->GetXaxis()->SetTitle(xtitl3[i].c_str());
1762           hists[jk]->GetYaxis()->SetTitle(ytitl3[i].c_str());
1763         }
1764         hists[jk]->GetYaxis()->SetLabelOffset(0.005);
1765         hists[jk]->GetYaxis()->SetLabelSize(0.035);
1766         hists[jk]->GetYaxis()->SetTitleSize(0.040);
1767         hists[jk]->GetYaxis()->SetTitleOffset(1.15);
1768         if ((type == 0) && (i != 3) && (i != 4))
1769           hists[jk]->GetXaxis()->SetRangeUser(0.0, 2.5);
1770         if (type == 0) {
1771           if (i == 3)
1772             hists[jk]->GetYaxis()->SetRangeUser(0.8, 1.2);
1773           else if (i == 4)
1774             hists[jk]->GetYaxis()->SetRangeUser(0.0, 0.5);
1775         }
1776         if (type == 1)
1777           hists[jk]->GetYaxis()->SetRangeUser(0.75, 1.2);
1778         hists[jk]->SetMarkerStyle(20);
1779         hists[jk]->SetMarkerColor(colors[k]);
1780         hists[jk]->SetLineColor(colors[k]);
1781         if (normalize && ((type == 2) || ((type == 0) && (i != 3)))) {
1782           if (jk == 0)
1783             hists[jk]->DrawNormalized();
1784           else
1785             hists[jk]->DrawNormalized("sames");
1786         } else {
1787           if (jk == 0)
1788             hists[jk]->Draw();
1789           else
1790             hists[jk]->Draw("sames");
1791         }
1792         pad->Update();
1793         TPaveStats* st1 = (TPaveStats*)hists[jk]->GetListOfFunctions()->FindObject("stats");
1794         if (st1 != nullptr) {
1795           double ymin = ymax - dy;
1796           st1->SetLineColor(colors[k]);
1797           st1->SetTextColor(colors[k]);
1798           st1->SetY1NDC(ymin);
1799           st1->SetY2NDC(ymax);
1800           st1->SetX1NDC(0.70);
1801           st1->SetX2NDC(0.90);
1802           ymax = ymin;
1803         }
1804         sprintf(name, "%s", texts[jk].c_str());
1805         legend->AddEntry(hists[jk], name, "lp");
1806       }
1807       legend->Draw("same");
1808       pad->Update();
1809       TPaveText* txt1 = new TPaveText(0.10, 0.905, 0.80, 0.95, "blNDC");
1810       txt1->SetFillColor(0);
1811       char txt[100];
1812       if (type == 2) {
1813         sprintf(txt, "p = 40:60 GeV (%s)", title3[i].c_str());
1814       } else if (((type == 0) && (iname == 0))) {
1815         sprintf(txt, "p = 20:30 GeV %s", text0.c_str());
1816       } else {
1817         sprintf(txt, "p = 40:60 GeV %s", text0.c_str());
1818       }
1819       txt1->AddText(txt);
1820       txt1->Draw("same");
1821       TPaveText* txt2 = new TPaveText(0.11, 0.825, 0.33, 0.895, "blNDC");
1822       txt2->SetFillColor(0);
1823       sprintf(txt, "CMS Preliminary");
1824       txt2->AddText(txt);
1825       txt2->Draw("same");
1826       pad->Modified();
1827       if ((drawStatBox == 0) && (i == 3)) {
1828         double xmin = hists[0]->GetBinLowEdge(1);
1829         int nbin = hists[0]->GetNbinsX();
1830         double xmax = hists[0]->GetBinLowEdge(nbin) + hists[0]->GetBinWidth(nbin);
1831         TLine* line = new TLine(xmin, 1.0, xmax, 1.0);  //etamin,1.0,etamax,1.0);
1832         line->SetLineWidth(2);
1833         line->Draw("same");
1834         pad->Update();
1835       }
1836       pad->Update();
1837       if (save > 0) {
1838         sprintf(name, "%s.pdf", pad->GetName());
1839         pad->Print(name);
1840       } else if (save < 0) {
1841         sprintf(name, "%s.C", pad->GetName());
1842         pad->Print(name);
1843       }
1844     }
1845   }
1846 }
1847 
1848 void PlotHistCorrResults(std::string infile, std::string text, std::string prefixF, int save = 0) {
1849   std::string name[5] = {"Eta1Bf", "Eta2Bf", "Eta1Af", "Eta2Af", "Cvg0"};
1850   std::string title[5] = {"Mean at the start of itertions",
1851                           "Median at the start of itertions",
1852                           "Mean at the end of itertions",
1853                           "Median at the end of itertions",
1854                           ""};
1855   int type[5] = {0, 0, 0, 0, 1};
1856 
1857   gStyle->SetCanvasBorderMode(0);
1858   gStyle->SetCanvasColor(kWhite);
1859   gStyle->SetPadColor(kWhite);
1860   gStyle->SetFillColor(kWhite);
1861   gStyle->SetOptTitle(0);
1862   gStyle->SetOptStat(10);
1863   gStyle->SetOptFit(10);
1864   TFile* file = new TFile(infile.c_str());
1865   char namep[100];
1866   for (int k = 0; k < 5; ++k) {
1867     TH1D* hist1 = (TH1D*)file->FindObjectAny(name[k].c_str());
1868     if (hist1 != nullptr) {
1869       TH1D* hist = (TH1D*)(hist1->Clone());
1870       sprintf(namep, "c_%s%s", prefixF.c_str(), name[k].c_str());
1871       TCanvas* pad = new TCanvas(namep, namep, 700, 500);
1872       pad->SetRightMargin(0.10);
1873       pad->SetTopMargin(0.10);
1874       hist->GetYaxis()->SetLabelOffset(0.005);
1875       hist->GetYaxis()->SetTitleOffset(1.20);
1876       double xmin = hist->GetBinLowEdge(1);
1877       int nbin = hist->GetNbinsX();
1878       double xmax = hist->GetBinLowEdge(nbin) + hist->GetBinWidth(nbin);
1879       std::cout << hist->GetTitle() << " Bins " << nbin << ":" << xmin << ":" << xmax << std::endl;
1880       double xlow(0.12), ylow(0.82);
1881       char txt[100], option[2];
1882       if (type[k] == 0) {
1883         sprintf(namep, "f_%s", name[k].c_str());
1884         TF1* func = new TF1(namep, "pol0", xmin, xmax);
1885         hist->Fit(func, "+QWL", "");
1886         if (text == "")
1887           sprintf(txt, "%s", title[k].c_str());
1888         else
1889           sprintf(txt, "%s (balancing the %s)", title[k].c_str(), text.c_str());
1890         sprintf(option, "%s", "");
1891       } else {
1892         hist->GetXaxis()->SetTitle("Iterations");
1893         hist->GetYaxis()->SetTitle("Convergence");
1894         hist->SetMarkerStyle(20);
1895         hist->SetMarkerColor(2);
1896         hist->SetMarkerSize(0.8);
1897         xlow = 0.50;
1898         ylow = 0.86;
1899         sprintf(txt, "(%s)", text.c_str());
1900         sprintf(option, "%s", "p");
1901       }
1902       hist->Draw(option);
1903       pad->Update();
1904       TPaveStats* st1 = (TPaveStats*)hist->GetListOfFunctions()->FindObject("stats");
1905       if (st1 != nullptr) {
1906         st1->SetY1NDC(ylow);
1907         st1->SetY2NDC(0.90);
1908         st1->SetX1NDC(0.70);
1909         st1->SetX2NDC(0.90);
1910       }
1911       TPaveText* txt1 = new TPaveText(xlow, 0.82, 0.68, 0.88, "blNDC");
1912       txt1->SetFillColor(0);
1913       txt1->AddText(txt);
1914       txt1->Draw("same");
1915       pad->Modified();
1916       pad->Update();
1917       if (save > 0) {
1918         sprintf(namep, "%s.pdf", pad->GetName());
1919         pad->Print(namep);
1920       } else if (save < 0) {
1921         sprintf(namep, "%s.C", pad->GetName());
1922         pad->Print(namep);
1923       }
1924     }
1925   }
1926 }
1927 
1928 void PlotHistCorrFactor(char* infile,
1929                         std::string text,
1930                         std::string prefixF = "",
1931                         double scale = 1.0,
1932                         int nmin = 100,
1933                         bool dataMC = false,
1934                         bool drawStatBox = true,
1935                         int iformat = 0,
1936                         int save = 0) {
1937   std::map<int, cfactors> cfacs;
1938   int etamin(100), etamax(-100), maxdepth(0);
1939   readCorrFactors(infile, scale, cfacs, etamin, etamax, maxdepth, iformat);
1940 
1941   gStyle->SetCanvasBorderMode(0);
1942   gStyle->SetCanvasColor(kWhite);
1943   gStyle->SetPadColor(kWhite);
1944   gStyle->SetFillColor(kWhite);
1945   gStyle->SetOptTitle(0);
1946   if (drawStatBox) {
1947     gStyle->SetOptStat(10);
1948     gStyle->SetOptFit(10);
1949   } else {
1950     gStyle->SetOptStat(0);
1951     gStyle->SetOptFit(0);
1952   }
1953   int colors[7] = {1, 6, 4, 7, 2, 9, 3};
1954   int mtype[7] = {20, 21, 22, 23, 24, 33, 25};
1955   int nbin = etamax - etamin + 1;
1956   std::vector<TH1D*> hists;
1957   std::vector<int> entries;
1958   char name[100];
1959   double dy(0);
1960   int fits(0);
1961   for (int j = 0; j < maxdepth; ++j) {
1962     sprintf(name, "hd%d", j + 1);
1963     TH1D* h = new TH1D(name, name, nbin, etamin, etamax);
1964     int nent(0);
1965     for (std::map<int, cfactors>::const_iterator itr = cfacs.begin(); itr != cfacs.end(); ++itr) {
1966       if ((itr->second).depth == j + 1) {
1967         int ieta = (itr->second).ieta;
1968         int bin = ieta - etamin + 1;
1969         float val = (itr->second).corrf;
1970         float dvl = (itr->second).dcorr;
1971         h->SetBinContent(bin, val);
1972         h->SetBinError(bin, dvl);
1973         nent++;
1974       }
1975     }
1976     if (nent > nmin) {
1977       fits++;
1978       dy += 0.025;
1979       sprintf(name, "hdf%d", j + 1);
1980       TF1* func = new TF1(name, "pol0", etamin, etamax);
1981       h->Fit(func, "+QWLR", "");
1982     }
1983     h->SetLineColor(colors[j]);
1984     h->SetMarkerColor(colors[j]);
1985     h->SetMarkerStyle(mtype[j]);
1986     h->GetXaxis()->SetTitle("i#eta");
1987     h->GetYaxis()->SetTitle("Correction Factor");
1988     h->GetYaxis()->SetLabelOffset(0.005);
1989     h->GetYaxis()->SetTitleOffset(1.20);
1990     h->GetYaxis()->SetRangeUser(0.0, 2.0);
1991     hists.push_back(h);
1992     entries.push_back(nent);
1993     dy += 0.025;
1994   }
1995   sprintf(name, "c_%sCorrFactor", prefixF.c_str());
1996   TCanvas* pad = new TCanvas(name, name, 700, 500);
1997   pad->SetRightMargin(0.10);
1998   pad->SetTopMargin(0.10);
1999   double yh = 0.90;
2000   // double yl = yh - 0.025 * hists.size() - dy - 0.01;
2001   double yl = 0.15;
2002   TLegend* legend = new TLegend(0.35, yl, 0.65, yl + 0.04 * hists.size());
2003   legend->SetFillColor(kWhite);
2004   for (unsigned int k = 0; k < hists.size(); ++k) {
2005     if (k == 0)
2006       hists[k]->Draw("");
2007     else
2008       hists[k]->Draw("sames");
2009     pad->Update();
2010     TPaveStats* st1 = (TPaveStats*)hists[k]->GetListOfFunctions()->FindObject("stats");
2011     if (st1 != nullptr) {
2012       dy = (entries[k] > nmin) ? 0.05 : 0.025;
2013       st1->SetLineColor(colors[k]);
2014       st1->SetTextColor(colors[k]);
2015       st1->SetY1NDC(yh - dy);
2016       st1->SetY2NDC(yh);
2017       st1->SetX1NDC(0.70);
2018       st1->SetX2NDC(0.90);
2019       yh -= dy;
2020     }
2021     sprintf(name, "Depth %d (%s)", k + 1, text.c_str());
2022     legend->AddEntry(hists[k], name, "lp");
2023   }
2024   legend->Draw("same");
2025   pad->Update();
2026   if (fits < 1) {
2027     double xmin = hists[0]->GetBinLowEdge(1);
2028     int nbin = hists[0]->GetNbinsX();
2029     double xmax = hists[0]->GetBinLowEdge(nbin) + hists[0]->GetBinWidth(nbin);
2030     TLine* line = new TLine(xmin, 1.0, xmax, 1.0);
2031     line->SetLineColor(9);
2032     line->SetLineWidth(2);
2033     line->SetLineStyle(2);
2034     line->Draw("same");
2035     pad->Modified();
2036     pad->Update();
2037   }
2038   char txt1[30];
2039   double xmax = (dataMC) ? 0.33 : 0.44;
2040   TPaveText* txt2 = new TPaveText(0.11, 0.85, xmax, 0.89, "blNDC");
2041   txt2->SetFillColor(0);
2042   if (dataMC)
2043     sprintf(txt1, "CMS Preliminary");
2044   else
2045     sprintf(txt1, "CMS Simulation Preliminary");
2046   txt2->AddText(txt1);
2047   txt2->Draw("same");
2048   pad->Modified();
2049   pad->Update();
2050   if (save > 0) {
2051     sprintf(name, "%s.pdf", pad->GetName());
2052     pad->Print(name);
2053   } else if (save < 0) {
2054     sprintf(name, "%s.C", pad->GetName());
2055     pad->Print(name);
2056   }
2057 }
2058 
2059 void PlotHistCorrAsymmetry(char* infile, std::string text, std::string prefixF = "", int iformat = 0, int save = 0) {
2060   std::map<int, cfactors> cfacs;
2061   int etamin(100), etamax(-100), maxdepth(0);
2062   double scale(1.0);
2063   readCorrFactors(infile, scale, cfacs, etamin, etamax, maxdepth, iformat);
2064 
2065   gStyle->SetCanvasBorderMode(0);
2066   gStyle->SetCanvasColor(kWhite);
2067   gStyle->SetPadColor(kWhite);
2068   gStyle->SetFillColor(kWhite);
2069   gStyle->SetOptTitle(0);
2070   gStyle->SetOptStat(0);
2071   gStyle->SetOptFit(10);
2072   int colors[6] = {1, 6, 4, 7, 2, 9};
2073   int mtype[6] = {20, 21, 22, 23, 24, 33};
2074   int nbin = etamax + 1;
2075   std::vector<TH1D*> hists;
2076   std::vector<int> entries;
2077   char name[100];
2078   double dy(0);
2079   for (int j = 0; j < maxdepth; ++j) {
2080     sprintf(name, "hd%d", j + 1);
2081     TH1D* h = new TH1D(name, name, nbin, 0, etamax);
2082     int nent(0);
2083     for (std::map<int, cfactors>::const_iterator itr = cfacs.begin(); itr != cfacs.end(); ++itr) {
2084       if ((itr->second).depth == j + 1) {
2085         int ieta = (itr->second).ieta;
2086         float vl1 = (itr->second).corrf;
2087         float dv1 = (itr->second).dcorr;
2088         if (ieta > 0) {
2089           for (std::map<int, cfactors>::const_iterator ktr = cfacs.begin(); ktr != cfacs.end(); ++ktr) {
2090             if (((ktr->second).depth == j + 1) && ((ktr->second).ieta == -ieta)) {
2091               float vl2 = (ktr->second).corrf;
2092               float dv2 = (ktr->second).dcorr;
2093               float val = 2.0 * (vl1 - vl2) / (vl1 + vl2);
2094               float dvl = (4.0 * sqrt(vl1 * vl1 * dv2 * dv2 + vl2 * vl2 * dv1 * dv1) / ((vl1 + vl2) * (vl1 + vl2)));
2095               int bin = ieta;
2096               h->SetBinContent(bin, val);
2097               h->SetBinError(bin, dvl);
2098               nent++;
2099             }
2100           }
2101         }
2102       }
2103     }
2104     h->SetLineColor(colors[j]);
2105     h->SetMarkerColor(colors[j]);
2106     h->SetMarkerStyle(mtype[j]);
2107     h->GetXaxis()->SetTitle("i#eta");
2108     h->GetYaxis()->SetTitle("Asymmetry in Correction Factor");
2109     h->GetYaxis()->SetLabelOffset(0.005);
2110     h->GetYaxis()->SetTitleOffset(1.20);
2111     h->GetYaxis()->SetRangeUser(-0.25, 0.25);
2112     hists.push_back(h);
2113     entries.push_back(nent);
2114     dy += 0.025;
2115   }
2116   sprintf(name, "c_%sCorrAsymmetry", prefixF.c_str());
2117   TCanvas* pad = new TCanvas(name, name, 700, 500);
2118   pad->SetRightMargin(0.10);
2119   pad->SetTopMargin(0.10);
2120   double yh = 0.90;
2121   double yl = yh - 0.035 * hists.size() - dy - 0.01;
2122   TLegend* legend = new TLegend(0.60, yl, 0.90, yl + 0.035 * hists.size());
2123   legend->SetFillColor(kWhite);
2124   for (unsigned int k = 0; k < hists.size(); ++k) {
2125     if (k == 0)
2126       hists[k]->Draw("");
2127     else
2128       hists[k]->Draw("sames");
2129     pad->Update();
2130     sprintf(name, "Depth %d (%s)", k + 1, text.c_str());
2131     legend->AddEntry(hists[k], name, "lp");
2132   }
2133   legend->Draw("same");
2134   pad->Update();
2135 
2136   TLine* line = new TLine(0.0, 0.0, etamax, 0.0);
2137   line->SetLineColor(9);
2138   line->SetLineWidth(2);
2139   line->SetLineStyle(2);
2140   line->Draw("same");
2141   pad->Update();
2142 
2143   if (save > 0) {
2144     sprintf(name, "%s.pdf", pad->GetName());
2145     pad->Print(name);
2146   } else if (save < 0) {
2147     sprintf(name, "%s.C", pad->GetName());
2148     pad->Print(name);
2149   }
2150 }
2151 
2152 void PlotHistCorrFactors(char* infile1,
2153                          std::string text1,
2154                          char* infile2,
2155                          std::string text2,
2156                          char* infile3,
2157                          std::string text3,
2158                          char* infile4,
2159                          std::string text4,
2160                          char* infile5,
2161                          std::string text5,
2162                          std::string prefixF,
2163                          bool ratio = false,
2164                          bool drawStatBox = true,
2165                          int nmin = 100,
2166                          bool dataMC = false,
2167                          int year = 2018,
2168                          int iformat = 0,
2169                          int save = 0) {
2170   std::map<int, cfactors> cfacs[5];
2171   std::vector<std::string> texts;
2172   int nfile(0), etamin(100), etamax(-100), maxdepth(0);
2173   const char* blank("");
2174   if (infile1 != blank) {
2175     readCorrFactors(infile1, 1.0, cfacs[nfile], etamin, etamax, maxdepth, iformat);
2176     if (cfacs[nfile].size() > 0) {
2177       texts.push_back(text1);
2178       ++nfile;
2179     }
2180   }
2181   if (infile2 != blank) {
2182     readCorrFactors(infile2, 1.0, cfacs[nfile], etamin, etamax, maxdepth, iformat);
2183     if (cfacs[nfile].size() > 0) {
2184       texts.push_back(text2);
2185       ++nfile;
2186     }
2187   }
2188   if (infile3 != blank) {
2189     readCorrFactors(infile3, 1.0, cfacs[nfile], etamin, etamax, maxdepth, iformat);
2190     if (cfacs[nfile].size() > 0) {
2191       texts.push_back(text3);
2192       ++nfile;
2193     }
2194   }
2195   if (infile4 != blank) {
2196     readCorrFactors(infile4, 1.0, cfacs[nfile], etamin, etamax, maxdepth, iformat);
2197     if (cfacs[nfile].size() > 0) {
2198       texts.push_back(text4);
2199       ++nfile;
2200     }
2201   }
2202   if (infile5 != blank) {
2203     readCorrFactors(infile5, 1.0, cfacs[nfile], etamin, etamax, maxdepth, iformat);
2204     if (cfacs[nfile].size() > 0) {
2205       texts.push_back(text5);
2206       ++nfile;
2207     }
2208   }
2209 
2210   if (nfile > 0) {
2211     gStyle->SetCanvasBorderMode(0);
2212     gStyle->SetCanvasColor(kWhite);
2213     gStyle->SetPadColor(kWhite);
2214     gStyle->SetFillColor(kWhite);
2215     gStyle->SetOptTitle(0);
2216     if ((!ratio) && drawStatBox) {
2217       gStyle->SetOptStat(10);
2218       gStyle->SetOptFit(10);
2219     } else {
2220       gStyle->SetOptStat(0);
2221       gStyle->SetOptFit(0);
2222     }
2223     int colors[6] = {1, 6, 4, 2, 7, 9};
2224     int mtype[6] = {20, 24, 22, 23, 21, 33};
2225     int nbin = etamax - etamin + 1;
2226     std::vector<TH1D*> hists;
2227     std::vector<int> entries, htype, depths;
2228     std::vector<double> fitr;
2229     char name[100];
2230     double dy(0);
2231     int fits(0);
2232     int nline(0);
2233     if (ratio) {
2234       for (int ih = 1; ih < nfile; ++ih) {
2235         for (int j = 0; j < maxdepth; ++j) {
2236           sprintf(name, "h%dd%d", ih, j + 1);
2237           TH1D* h = new TH1D(name, name, nbin, etamin, etamax);
2238           double sumNum(0), sumDen(0);
2239           std::map<int, cfactors>::const_iterator ktr = cfacs[ih].begin();
2240           for (std::map<int, cfactors>::const_iterator itr = cfacs[0].begin(); itr != cfacs[0].end(); ++itr, ++ktr) {
2241             int dep = (itr->second).depth;
2242             if (dep == j + 1) {
2243               int ieta = (itr->second).ieta;
2244               int bin = ieta - etamin + 1;
2245               float val = (itr->second).corrf / (ktr->second).corrf;
2246               float dvl =
2247                   val *
2248                   sqrt((((itr->second).dcorr * (itr->second).dcorr) / ((itr->second).corrf * (itr->second).corrf)) +
2249                        (((ktr->second).dcorr * (ktr->second).dcorr) / ((ktr->second).corrf * (ktr->second).corrf)));
2250               h->SetBinContent(bin, val);
2251               h->SetBinError(bin, dvl);
2252               sumNum += (val / (dvl * dvl));
2253               sumDen += (1.0 / (dvl * dvl));
2254             }
2255           }
2256           double fit = (sumDen > 0) ? (sumNum / sumDen) : 1.0;
2257           std::cout << "Fit to Pol0: " << fit << std::endl;
2258           h->SetLineColor(colors[ih]);
2259           h->SetMarkerColor(colors[ih]);
2260           h->SetMarkerStyle(mtype[j]);
2261           h->SetMarkerSize(0.9);
2262           h->GetXaxis()->SetTitle("i#eta");
2263           sprintf(name, "CF_{%s}/CF_{Set}", texts[0].c_str());
2264           h->GetYaxis()->SetTitle(name);
2265           h->GetYaxis()->SetLabelOffset(0.005);
2266           h->GetYaxis()->SetTitleSize(0.036);
2267           h->GetYaxis()->SetTitleOffset(1.20);
2268           h->GetYaxis()->SetRangeUser(0.80, 1.20);
2269           hists.push_back(h);
2270           fitr.push_back(fit);
2271           htype.push_back(ih);
2272           depths.push_back(j + 1);
2273         }
2274         if (ih == 1)
2275           nline += hists.size();
2276         else
2277           ++nline;
2278       }
2279     } else {
2280       for (int k1 = 0; k1 < nfile; ++k1) {
2281         for (int j = 0; j < maxdepth; ++j) {
2282           sprintf(name, "h%dd%d", k1, j + 1);
2283           TH1D* h = new TH1D(name, name, nbin, etamin, etamax);
2284           int nent(0);
2285           for (std::map<int, cfactors>::const_iterator itr = cfacs[k1].begin(); itr != cfacs[k1].end(); ++itr) {
2286             int dep = (itr->second).depth;
2287             if (dep == j + 1) {
2288               int ieta = (itr->second).ieta;
2289               int bin = ieta - etamin + 1;
2290               float val = (itr->second).corrf;
2291               float dvl = (itr->second).dcorr;
2292               h->SetBinContent(bin, val);
2293               h->SetBinError(bin, dvl);
2294               nent++;
2295             }
2296           }
2297           if (nent > nmin) {
2298             fits++;
2299             if (drawStatBox)
2300               dy += 0.025;
2301             sprintf(name, "h%ddf%d", k1, j + 1);
2302             TF1* func = new TF1(name, "pol0", etamin, etamax);
2303             h->Fit(func, "+QWLR", "");
2304           }
2305           h->SetLineColor(colors[k1]);
2306           h->SetMarkerColor(colors[k1]);
2307           h->SetMarkerStyle(mtype[j]);
2308           h->SetMarkerSize(0.9);
2309           h->GetXaxis()->SetTitle("i#eta");
2310           h->GetYaxis()->SetTitle("Correction Factor");
2311           h->GetYaxis()->SetLabelOffset(0.005);
2312           h->GetYaxis()->SetTitleOffset(1.20);
2313           h->GetYaxis()->SetRangeUser(0.5, 1.5);
2314           hists.push_back(h);
2315           entries.push_back(nent);
2316           if (drawStatBox)
2317             dy += 0.025;
2318           htype.push_back(k1);
2319           depths.push_back(j + 1);
2320         }
2321         if (k1 == 0)
2322           nline += hists.size();
2323         else
2324           ++nline;
2325       }
2326     }
2327     if (ratio)
2328       sprintf(name, "c_Corr%sRatio", prefixF.c_str());
2329     else
2330       sprintf(name, "c_Corr%s", prefixF.c_str());
2331     TCanvas* pad = new TCanvas(name, name, 700, 500);
2332     pad->SetRightMargin(0.10);
2333     pad->SetTopMargin(0.10);
2334     double yh = 0.90;
2335     double yl = yh - 0.035 * hists.size() - dy - 0.01;
2336     TLegend* legend = new TLegend(0.55, yl, 0.90, yl + 0.035 * nline);
2337     legend->SetFillColor(kWhite);
2338     for (unsigned int k = 0; k < hists.size(); ++k) {
2339       if (k == 0)
2340         hists[k]->Draw("");
2341       else
2342         hists[k]->Draw("sames");
2343       pad->Update();
2344       int k1 = htype[k];
2345       if (!ratio) {
2346         TPaveStats* st1 = (TPaveStats*)hists[k]->GetListOfFunctions()->FindObject("stats");
2347         if (st1 != nullptr) {
2348           dy = (entries[k] > nmin) ? 0.05 : 0.025;
2349           st1->SetLineColor(colors[k1]);
2350           st1->SetTextColor(colors[k1]);
2351           st1->SetY1NDC(yh - dy);
2352           st1->SetY2NDC(yh);
2353           st1->SetX1NDC(0.70);
2354           st1->SetX2NDC(0.90);
2355           yh -= dy;
2356         }
2357         sprintf(name, "Depth %d (%s)", depths[k], texts[k1].c_str());
2358       } else {
2359         sprintf(name, "Depth %d (%s Mean = %5.3f)", depths[k], texts[k1].c_str(), fitr[k]);
2360       }
2361       if ((depths[k] == 1) || (k1 == 0))
2362         legend->AddEntry(hists[k], name, "lp");
2363     }
2364     legend->Draw("same");
2365     TPaveText* txt0 = new TPaveText(0.12, 0.84, 0.49, 0.89, "blNDC");
2366     txt0->SetFillColor(0);
2367     char txt[40];
2368     if (dataMC)
2369       sprintf(txt, "CMS Preliminary (%d)", year);
2370     else
2371       sprintf(txt, "CMS Simulation Preliminary (%d)", year);
2372     txt0->AddText(txt);
2373     txt0->Draw("same");
2374     pad->Update();
2375     if (fits < 1) {
2376       double xmin = hists[0]->GetBinLowEdge(1);
2377       int nbin = hists[0]->GetNbinsX();
2378       double xmax = hists[0]->GetBinLowEdge(nbin) + hists[0]->GetBinWidth(nbin);
2379       TLine* line = new TLine(xmin, 1.0, xmax, 1.0);
2380       line->SetLineColor(9);
2381       line->SetLineWidth(2);
2382       line->SetLineStyle(2);
2383       line->Draw("same");
2384       pad->Update();
2385     }
2386     if (save > 0) {
2387       sprintf(name, "%s.pdf", pad->GetName());
2388       pad->Print(name);
2389     } else if (save < 0) {
2390       sprintf(name, "%s.C", pad->GetName());
2391       pad->Print(name);
2392     }
2393   }
2394 }
2395 
2396 void PlotHistCorrSys(std::string infilec, int conds, std::string text, int save = 0) {
2397   char fname[100];
2398   int iformat(0);
2399   sprintf(fname, "%s_cond0.txt", infilec.c_str());
2400   int etamin(100), etamax(-100), maxdepth(0);
2401   std::map<int, cfactors> cfacs;
2402   readCorrFactors(fname, 1.0, cfacs, etamin, etamax, maxdepth, iformat);
2403   // There are good records from the master file
2404   if (cfacs.size() > 0) {
2405     // Now read the other files
2406     std::map<int, cfactors> errfacs;
2407     for (int i = 0; i < conds; ++i) {
2408       sprintf(fname, "%s_cond%d.txt", infilec.c_str(), i + 1);
2409       std::map<int, cfactors> cfacx;
2410       int etamin1(100), etamax1(-100), maxdepth1(0);
2411       readCorrFactors(fname, 1.0, cfacx, etamin1, etamax1, maxdepth1, iformat);
2412       for (std::map<int, cfactors>::const_iterator itr1 = cfacx.begin(); itr1 != cfacx.end(); ++itr1) {
2413         std::map<int, cfactors>::iterator itr2 = errfacs.find(itr1->first);
2414         float corrf = (itr1->second).corrf;
2415         if (itr2 == errfacs.end()) {
2416           errfacs[itr1->first] = cfactors(1, 0, corrf, corrf * corrf);
2417         } else {
2418           int nent = (itr2->second).ieta + 1;
2419           float c1 = (itr2->second).corrf + corrf;
2420           float c2 = (itr2->second).dcorr + (corrf * corrf);
2421           errfacs[itr1->first] = cfactors(nent, 0, c1, c2);
2422         }
2423       }
2424     }
2425     // find the RMS from the distributions
2426     for (std::map<int, cfactors>::iterator itr = errfacs.begin(); itr != errfacs.end(); ++itr) {
2427       int nent = (itr->second).ieta;
2428       float mean = (itr->second).corrf / nent;
2429       float rms2 = (itr->second).dcorr / nent - (mean * mean);
2430       float rms = rms2 > 0 ? sqrt(rms2) : 0;
2431       errfacs[itr->first] = cfactors(nent, 0, mean, rms);
2432     }
2433     // Now combine the errors and plot
2434     int k(0);
2435     gStyle->SetCanvasBorderMode(0);
2436     gStyle->SetCanvasColor(kWhite);
2437     gStyle->SetPadColor(kWhite);
2438     gStyle->SetFillColor(kWhite);
2439     gStyle->SetOptTitle(0);
2440     gStyle->SetOptStat(10);
2441     gStyle->SetOptFit(10);
2442     int colors[6] = {1, 6, 4, 7, 2, 9};
2443     int mtype[6] = {20, 21, 22, 23, 24, 33};
2444     std::vector<TH1D*> hists;
2445     char name[100];
2446     int nbin = etamax - etamin + 1;
2447     for (int j = 0; j < maxdepth; ++j) {
2448       sprintf(name, "hd%d", j + 1);
2449       TH1D* h = new TH1D(name, name, nbin, etamin, etamax);
2450       h->SetLineColor(colors[j]);
2451       h->SetMarkerColor(colors[j]);
2452       h->SetMarkerStyle(mtype[j]);
2453       h->GetXaxis()->SetTitle("i#eta");
2454       h->GetYaxis()->SetTitle("Correction Factor");
2455       h->GetYaxis()->SetLabelOffset(0.005);
2456       h->GetYaxis()->SetTitleOffset(1.20);
2457       h->GetYaxis()->SetRangeUser(0.0, 2.0);
2458       hists.push_back(h);
2459     }
2460     for (std::map<int, cfactors>::iterator itr = cfacs.begin(); itr != cfacs.end(); ++itr, ++k) {
2461       std::map<int, cfactors>::iterator itr2 = errfacs.find(itr->first);
2462       float mean2 = (itr2 == errfacs.end()) ? 0 : (itr2->second).corrf;
2463       float ersys = (itr2 == errfacs.end()) ? 0 : (itr2->second).dcorr;
2464       float erstt = (itr->second).dcorr;
2465       float ertot = sqrt(erstt * erstt + ersys * ersys);
2466       float mean = (itr->second).corrf;
2467       int ieta = (itr->second).ieta;
2468       int depth = (itr->second).depth;
2469       std::cout << "[" << k << "] " << ieta << " " << depth << " " << mean << ":" << mean2 << " " << erstt << ":"
2470                 << ersys << ":" << ertot << std::endl;
2471       int bin = ieta - etamin + 1;
2472       hists[depth - 1]->SetBinContent(bin, mean);
2473       hists[depth - 1]->SetBinError(bin, ertot);
2474     }
2475     TCanvas* pad = new TCanvas("CorrFactorSys", "CorrFactorSys", 700, 500);
2476     pad->SetRightMargin(0.10);
2477     pad->SetTopMargin(0.10);
2478     double yh = 0.90;
2479     double yl = yh - 0.050 * hists.size() - 0.01;
2480     TLegend* legend = new TLegend(0.60, yl, 0.90, yl + 0.025 * hists.size());
2481     legend->SetFillColor(kWhite);
2482     for (unsigned int k = 0; k < hists.size(); ++k) {
2483       if (k == 0)
2484         hists[k]->Draw("");
2485       else
2486         hists[k]->Draw("sames");
2487       pad->Update();
2488       TPaveStats* st1 = (TPaveStats*)hists[k]->GetListOfFunctions()->FindObject("stats");
2489       if (st1 != nullptr) {
2490         st1->SetLineColor(colors[k]);
2491         st1->SetTextColor(colors[k]);
2492         st1->SetY1NDC(yh - 0.025);
2493         st1->SetY2NDC(yh);
2494         st1->SetX1NDC(0.70);
2495         st1->SetX2NDC(0.90);
2496         yh -= 0.025;
2497       }
2498       sprintf(name, "Depth %d (%s)", k + 1, text.c_str());
2499       legend->AddEntry(hists[k], name, "lp");
2500     }
2501     legend->Draw("same");
2502     pad->Update();
2503     if (save > 0) {
2504       sprintf(name, "%s.pdf", pad->GetName());
2505       pad->Print(name);
2506     } else if (save < 0) {
2507       sprintf(name, "%s.C", pad->GetName());
2508       pad->Print(name);
2509     }
2510   }
2511 }
2512 
2513 void PlotHistCorrLumis(std::string infilec, int conds, double lumi, int save = 0) {
2514   char fname[100];
2515   int iformat(0);
2516   sprintf(fname, "%s_0.txt", infilec.c_str());
2517   std::map<int, cfactors> cfacs;
2518   int etamin(100), etamax(-100), maxdepth(0);
2519   readCorrFactors(fname, 1.0, cfacs, etamin, etamax, maxdepth, iformat);
2520   int nbin = etamax - etamin + 1;
2521   std::cout << "Max Depth " << maxdepth << " and " << nbin << " eta bins for " << etamin << ":" << etamax << std::endl;
2522 
2523   // There are good records from the master file
2524   int colors[8] = {4, 2, 6, 7, 1, 9, 3, 5};
2525   int mtype[8] = {20, 21, 22, 23, 24, 25, 26, 27};
2526   if (cfacs.size() > 0) {
2527     // Now read the other files
2528     std::vector<TH1D*> hists;
2529     char name[100];
2530     for (int i = 0; i < conds; ++i) {
2531       int ih = (int)(hists.size());
2532       for (int j = 0; j < maxdepth; ++j) {
2533         sprintf(name, "hd%d%d", j + 1, i);
2534         TH1D* h = new TH1D(name, name, nbin, etamin, etamax);
2535         h->SetLineColor(colors[j]);
2536         h->SetMarkerColor(colors[j]);
2537         h->SetMarkerStyle(mtype[i]);
2538         h->SetMarkerSize(0.9);
2539         h->GetXaxis()->SetTitle("i#eta");
2540         h->GetYaxis()->SetTitle("Fractional Error");
2541         h->GetYaxis()->SetLabelOffset(0.005);
2542         h->GetYaxis()->SetTitleOffset(1.20);
2543         h->GetYaxis()->SetRangeUser(0.0, 0.10);
2544         hists.push_back(h);
2545       }
2546       sprintf(fname, "%s_%d.txt", infilec.c_str(), i);
2547       int etamin1(100), etamax1(-100), maxdepth1(0);
2548       readCorrFactors(fname, 1.0, cfacs, etamin1, etamax1, maxdepth1, iformat);
2549       for (std::map<int, cfactors>::const_iterator itr = cfacs.begin(); itr != cfacs.end(); ++itr) {
2550         double value = (itr->second).dcorr / (itr->second).corrf;
2551         int bin = (itr->second).ieta - etamin + 1;
2552         hists[ih + (itr->second).depth - 1]->SetBinContent(bin, value);
2553         hists[ih + (itr->second).depth - 1]->SetBinError(bin, 0.0001);
2554       }
2555     }
2556 
2557     // Now plot the histograms
2558     gStyle->SetCanvasBorderMode(0);
2559     gStyle->SetCanvasColor(kWhite);
2560     gStyle->SetPadColor(kWhite);
2561     gStyle->SetFillColor(kWhite);
2562     gStyle->SetOptTitle(0);
2563     gStyle->SetOptStat(0);
2564     gStyle->SetOptFit(0);
2565     TCanvas* pad = new TCanvas("CorrFactorErr", "CorrFactorErr", 700, 500);
2566     pad->SetRightMargin(0.10);
2567     pad->SetTopMargin(0.10);
2568     double yh(0.89);
2569     TLegend* legend = new TLegend(0.60, yh - 0.04 * conds, 0.89, yh);
2570     legend->SetFillColor(kWhite);
2571     double lumic(lumi);
2572     for (unsigned int k = 0; k < hists.size(); ++k) {
2573       if (k == 0)
2574         hists[k]->Draw("");
2575       else
2576         hists[k]->Draw("sames");
2577       pad->Update();
2578       if (k % maxdepth == 0) {
2579         sprintf(name, "L = %5.2f fb^{-1}", lumic);
2580         legend->AddEntry(hists[k], name, "lp");
2581         lumic *= 0.5;
2582       }
2583     }
2584     legend->Draw("same");
2585     pad->Update();
2586     if (save > 0) {
2587       sprintf(name, "%s.pdf", pad->GetName());
2588       pad->Print(name);
2589     } else if (save < 0) {
2590       sprintf(name, "%s.C", pad->GetName());
2591       pad->Print(name);
2592     }
2593   }
2594 }
2595 
2596 void PlotHistCorrRel(char* infile1,
2597                      char* infile2,
2598                      std::string text1,
2599                      std::string text2,
2600                      int iformat1 = 0,
2601                      int iformat2 = 0,
2602                      int save = 0) {
2603   std::map<int, cfactors> cfacs1, cfacs2;
2604   int etamin(100), etamax(-100), maxdepth(0);
2605   readCorrFactors(infile1, 1.0, cfacs1, etamin, etamax, maxdepth, iformat1);
2606   readCorrFactors(infile2, 1.0, cfacs2, etamin, etamax, maxdepth, iformat2);
2607   std::map<int, std::pair<cfactors, cfactors> > cfacs;
2608   for (std::map<int, cfactors>::iterator itr = cfacs1.begin(); itr != cfacs1.end(); ++itr) {
2609     std::map<int, cfactors>::iterator ktr = cfacs2.find(itr->first);
2610     if (ktr == cfacs2.end()) {
2611       cfactors fac2(((itr->second).ieta), ((itr->second).depth), 0, -1);
2612       cfacs[itr->first] = std::pair<cfactors, cfactors>((itr->second), fac2);
2613     } else {
2614       cfactors fac2(ktr->second);
2615       cfacs[itr->first] = std::pair<cfactors, cfactors>((itr->second), fac2);
2616     }
2617   }
2618   for (std::map<int, cfactors>::iterator itr = cfacs2.begin(); itr != cfacs2.end(); ++itr) {
2619     std::map<int, cfactors>::const_iterator ktr = cfacs1.find(itr->first);
2620     if (ktr == cfacs1.end()) {
2621       cfactors fac1(((itr->second).ieta), ((itr->second).depth), 0, -1);
2622       cfacs[itr->first] = std::pair<cfactors, cfactors>(fac1, (itr->second));
2623     }
2624   }
2625 
2626   // There are good records in bothe the files
2627   if ((cfacs1.size() > 0) && (cfacs2.size() > 0)) {
2628     int k(0);
2629     gStyle->SetCanvasBorderMode(0);
2630     gStyle->SetCanvasColor(kWhite);
2631     gStyle->SetPadColor(kWhite);
2632     gStyle->SetFillColor(kWhite);
2633     gStyle->SetOptTitle(0);
2634     gStyle->SetOptStat(10);
2635     gStyle->SetOptFit(10);
2636     int colors[6] = {1, 6, 4, 7, 2, 9};
2637     int mtype[6] = {20, 21, 22, 23, 24, 33};
2638     std::vector<TH1D*> hists;
2639     char name[100];
2640     int nbin = etamax - etamin + 1;
2641     for (int i = 0; i < 2; ++i) {
2642       for (int j = 0; j < maxdepth; ++j) {
2643         int j1 = (i == 0) ? j : maxdepth + j;
2644         sprintf(name, "hd%d%d", i, j + 1);
2645         TH1D* h = new TH1D(name, name, nbin, etamin, etamax);
2646         h->SetLineColor(colors[j1]);
2647         h->SetMarkerColor(colors[j1]);
2648         h->SetMarkerStyle(mtype[i]);
2649         h->GetXaxis()->SetTitle("i#eta");
2650         h->GetYaxis()->SetTitle("Correction Factor");
2651         h->GetYaxis()->SetLabelOffset(0.005);
2652         h->GetYaxis()->SetTitleOffset(1.20);
2653         h->GetYaxis()->SetRangeUser(0.0, 2.0);
2654         hists.push_back(h);
2655       }
2656     }
2657     for (std::map<int, std::pair<cfactors, cfactors> >::iterator it = cfacs.begin(); it != cfacs.end(); ++it, ++k) {
2658       float mean1 = (it->second).first.corrf;
2659       float error1 = (it->second).first.dcorr;
2660       float mean2 = (it->second).second.corrf;
2661       float error2 = (it->second).second.dcorr;
2662       int ieta = (it->second).first.ieta;
2663       int depth = (it->second).first.depth;
2664       /*
2665       std::cout << "[" << k << "] " << ieta << " " << depth << " " 
2666         << mean1 << ":" << mean2 << " " << error1 << ":" << error2
2667         << std::endl;
2668       */
2669       int bin = ieta - etamin + 1;
2670       if (error1 >= 0) {
2671         hists[depth - 1]->SetBinContent(bin, mean1);
2672         hists[depth - 1]->SetBinError(bin, error1);
2673       }
2674       if (error2 >= 0) {
2675         hists[maxdepth + depth - 1]->SetBinContent(bin, mean2);
2676         hists[maxdepth + depth - 1]->SetBinError(bin, error2);
2677       }
2678     }
2679     TCanvas* pad = new TCanvas("CorrFactors", "CorrFactors", 700, 500);
2680     pad->SetRightMargin(0.10);
2681     pad->SetTopMargin(0.10);
2682     double yh = 0.90;
2683     double yl = yh - 0.050 * hists.size() - 0.01;
2684     TLegend* legend = new TLegend(0.60, yl, 0.90, yl + 0.025 * hists.size());
2685     legend->SetFillColor(kWhite);
2686     for (unsigned int k = 0; k < hists.size(); ++k) {
2687       if (k == 0)
2688         hists[k]->Draw("");
2689       else
2690         hists[k]->Draw("sames");
2691       pad->Update();
2692       TPaveStats* st1 = (TPaveStats*)hists[k]->GetListOfFunctions()->FindObject("stats");
2693       if (st1 != nullptr) {
2694         st1->SetLineColor(colors[k]);
2695         st1->SetTextColor(colors[k]);
2696         st1->SetY1NDC(yh - 0.025);
2697         st1->SetY2NDC(yh);
2698         st1->SetX1NDC(0.70);
2699         st1->SetX2NDC(0.90);
2700         yh -= 0.025;
2701       }
2702       if (k < (unsigned int)(maxdepth)) {
2703         sprintf(name, "Depth %d (%s)", k + 1, text1.c_str());
2704       } else {
2705         sprintf(name, "Depth %d (%s)", k - maxdepth + 1, text2.c_str());
2706       }
2707       legend->AddEntry(hists[k], name, "lp");
2708     }
2709     legend->Draw("same");
2710     pad->Update();
2711     if (save > 0) {
2712       sprintf(name, "%s.pdf", pad->GetName());
2713       pad->Print(name);
2714     } else if (save < 0) {
2715       sprintf(name, "%s.C", pad->GetName());
2716       pad->Print(name);
2717     }
2718   }
2719 }
2720 
2721 void PlotHistCorrDepth(char* infile1,
2722                        char* infile2,
2723                        std::string text1,
2724                        std::string text2,
2725                        int depth,
2726                        int ietamax,
2727                        int iformat1 = 0,
2728                        int iformat2 = 0,
2729                        int save = 0,
2730                        int debug = 1) {
2731   std::map<int, cfactors> cfacs1, cfacs2;
2732   int etamin(100), etamax(-100), maxdepth(0);
2733   readCorrFactors(infile1, 1.0, cfacs1, etamin, etamax, maxdepth, iformat1, (debug > 1));
2734   readCorrFactors(infile2, 1.0, cfacs2, etamin, etamax, maxdepth, iformat2, (debug > 1));
2735 
2736   double sumNum(0), sumDen(0), ratMax(0);
2737   int npt0(0), npt1(0);
2738   for (std::map<int, cfactors>::iterator itr = cfacs1.begin(); itr != cfacs1.end(); ++itr) {
2739     std::map<int, cfactors>::iterator ktr = cfacs2.find(itr->first);
2740     if ((ktr != cfacs2.end()) && ((itr->second).depth == depth)) {
2741       double er1 = (itr->second).dcorr / (itr->second).corrf;
2742       double er2 = (ktr->second).dcorr / (ktr->second).corrf;
2743       double tmp = (ktr->second).corrf / (itr->second).corrf;
2744       double dtmp = tmp * sqrt(er1 * er1 + er2 * er2);
2745       double rat = (tmp > 1.0) ? 1.0 / tmp : tmp;
2746       double drt = (tmp > 1.0) ? dtmp / (tmp * tmp) : dtmp;
2747       rat = fabs(1.0 - rat);
2748       ratMax = std::max(ratMax, rat);
2749       ++npt0;
2750       if (debug > 0)
2751         std::cout << std::hex << (itr->first) << std::dec << " ieta:depth" << (itr->second).ieta << ":"
2752                   << (itr->second).depth << " Corr " << (itr->second).corrf << ":" << (ktr->second).corrf << " Ratio "
2753                   << rat << ":" << drt << std::endl;
2754       if (std::abs((itr->second).ieta) <= ietamax) {
2755         sumNum += (rat / (drt * drt));
2756         sumDen += (1.0 / (drt * drt));
2757         ++npt1;
2758       }
2759     }
2760   }
2761   sumNum = (sumDen > 0) ? (sumNum / sumDen) : 0;
2762   sumDen = (sumDen > 0) ? 1.0 / sqrt(sumDen) : 0;
2763   std::cout << "Get Ratio of mean for " << npt0 << ":" << npt1 << " points for depth " << depth << " Mean " << sumNum
2764             << " +- " << sumDen << " (Maximum " << ratMax << ")" << std::endl;
2765 
2766   gStyle->SetCanvasBorderMode(0);
2767   gStyle->SetCanvasColor(kWhite);
2768   gStyle->SetPadColor(kWhite);
2769   gStyle->SetFillColor(kWhite);
2770   gStyle->SetOptTitle(0);
2771   gStyle->SetOptStat(0);
2772   gStyle->SetOptFit(0);
2773   int colors[2] = {1, 2};
2774   int mtype[2] = {20, 24};
2775   int nbin = etamax - etamin + 1;
2776   std::vector<TH1D*> hists;
2777   char name[100];
2778   for (int j = 0; j < 2; ++j) {
2779     sprintf(name, "hd%d", (j + 1));
2780     TH1D* h = new TH1D(name, name, nbin, etamin, etamax);
2781     if (j == 0) {
2782       for (std::map<int, cfactors>::const_iterator itr = cfacs1.begin(); itr != cfacs1.end(); ++itr) {
2783         if ((itr->second).depth == depth) {
2784           int ieta = (itr->second).ieta;
2785           int bin = ieta - etamin + 1;
2786           float val = (itr->second).corrf;
2787           float dvl = (itr->second).dcorr;
2788           h->SetBinContent(bin, val);
2789           h->SetBinError(bin, dvl);
2790         }
2791       }
2792     } else {
2793       for (std::map<int, cfactors>::const_iterator itr = cfacs2.begin(); itr != cfacs2.end(); ++itr) {
2794         if ((itr->second).depth == depth) {
2795           int ieta = (itr->second).ieta;
2796           int bin = ieta - etamin + 1;
2797           float val = (itr->second).corrf;
2798           float dvl = (itr->second).dcorr;
2799           h->SetBinContent(bin, val);
2800           h->SetBinError(bin, dvl);
2801         }
2802       }
2803     }
2804     h->SetLineColor(colors[j]);
2805     h->SetMarkerColor(colors[j]);
2806     h->SetMarkerStyle(mtype[j]);
2807     h->GetXaxis()->SetTitle("i#eta");
2808     h->GetYaxis()->SetTitle("Correction Factor");
2809     h->GetYaxis()->SetLabelOffset(0.005);
2810     h->GetYaxis()->SetTitleOffset(1.20);
2811     h->GetYaxis()->SetRangeUser(0.0, 2.0);
2812     hists.push_back(h);
2813   }
2814   sprintf(name, "c_CorrFactorDepth%d", depth);
2815   TCanvas* pad = new TCanvas(name, name, 700, 500);
2816   pad->SetRightMargin(0.10);
2817   pad->SetTopMargin(0.10);
2818   double yl = 0.25;
2819   TLegend* legend = new TLegend(0.25, yl, 0.85, yl + 0.04 * hists.size());
2820   legend->SetFillColor(kWhite);
2821   for (unsigned int k = 0; k < hists.size(); ++k) {
2822     if (k == 0) {
2823       hists[k]->Draw("");
2824       sprintf(name, "Depth %d (%s)", depth, text1.c_str());
2825     } else {
2826       hists[k]->Draw("sames");
2827       sprintf(name, "Depth %d (%s)", depth, text2.c_str());
2828     }
2829     pad->Update();
2830     legend->AddEntry(hists[k], name, "lp");
2831   }
2832   legend->Draw("same");
2833   pad->Update();
2834   if (save > 0) {
2835     sprintf(name, "%s.pdf", pad->GetName());
2836     pad->Print(name);
2837   } else if (save < 0) {
2838     sprintf(name, "%s.C", pad->GetName());
2839     pad->Print(name);
2840   }
2841 }
2842 
2843 void PlotFourHists(std::string infile,
2844                    std::string prefix0,
2845                    int type = 0,
2846                    int drawStatBox = 0,
2847                    bool normalize = false,
2848                    int save = 0,
2849                    std::string prefix1 = "",
2850                    std::string text1 = "",
2851                    std::string prefix2 = "",
2852                    std::string text2 = "",
2853                    std::string prefix3 = "",
2854                    std::string text3 = "",
2855                    std::string prefix4 = "",
2856                    std::string text4 = "") {
2857   int colors[4] = {2, 4, 6, 1};
2858   std::string names[5] = {"eta03", "eta13", "eta23", "eta33", "eta43"};
2859   std::string xtitle[5] = {"i#eta", "i#eta", "i#eta", "i#eta", "i#eta"};
2860   std::string ytitle[5] = {"Tracks", "Tracks", "Tracks", "Tracks", "Tracks"};
2861   std::string title[5] = {"All Tracks (p = 40:60 GeV)",
2862                           "Good Quality Tracks (p = 40:60 GeV)",
2863                           "Selected Tracks (p = 40:60 GeV)",
2864                           "Isolated Tracks (p = 40:60 GeV)",
2865                           "Isolated MIP Tracks (p = 40:60 GeV)"};
2866 
2867   gStyle->SetCanvasBorderMode(0);
2868   gStyle->SetCanvasColor(kWhite);
2869   gStyle->SetPadColor(kWhite);
2870   gStyle->SetFillColor(kWhite);
2871   gStyle->SetOptTitle(0);
2872   gStyle->SetOptFit(0);
2873   if (drawStatBox == 0)
2874     gStyle->SetOptStat(0);
2875   else
2876     gStyle->SetOptStat(1110);
2877 
2878   if (type < 0 || type > 4)
2879     type = 0;
2880   char name[100], namep[100];
2881   TFile* file = new TFile(infile.c_str());
2882 
2883   std::vector<TH1D*> hists;
2884   std::vector<int> kks;
2885   std::vector<std::string> texts;
2886   for (int k = 0; k < 4; ++k) {
2887     std::string prefix, text;
2888     if (k == 0) {
2889       prefix = prefix1;
2890       text = text1;
2891     } else if (k == 1) {
2892       prefix = prefix2;
2893       text = text2;
2894     } else if (k == 2) {
2895       prefix = prefix3;
2896       text = text3;
2897     } else {
2898       prefix = prefix4;
2899       text = text4;
2900     }
2901     if (prefix != "") {
2902       sprintf(name, "%s%s", prefix.c_str(), names[type].c_str());
2903       TH1D* hist1 = (TH1D*)file->FindObjectAny(name);
2904       if (hist1 != nullptr) {
2905         hists.push_back((TH1D*)(hist1->Clone()));
2906         kks.push_back(k);
2907         texts.push_back(text);
2908       }
2909     }
2910   }
2911   if (hists.size() > 0) {
2912     sprintf(namep, "c_%s%s", prefix0.c_str(), names[type].c_str());
2913     double ymax(0.90), dy(0.13);
2914     double ymx0 = (drawStatBox == 0) ? (ymax - .01) : (ymax - dy * hists.size() - .01);
2915     TCanvas* pad = new TCanvas(namep, namep, 700, 500);
2916     TLegend* legend = new TLegend(0.64, ymx0 - 0.05 * hists.size(), 0.89, ymx0);
2917     legend->SetFillColor(kWhite);
2918     pad->SetRightMargin(0.10);
2919     pad->SetTopMargin(0.10);
2920     for (unsigned int jk = 0; jk < hists.size(); ++jk) {
2921       int k = kks[jk];
2922       hists[jk]->GetXaxis()->SetTitleSize(0.040);
2923       hists[jk]->GetXaxis()->SetTitle(xtitle[type].c_str());
2924       hists[jk]->GetYaxis()->SetTitle(ytitle[type].c_str());
2925       hists[jk]->GetYaxis()->SetLabelOffset(0.005);
2926       hists[jk]->GetYaxis()->SetLabelSize(0.035);
2927       hists[jk]->GetYaxis()->SetTitleSize(0.040);
2928       hists[jk]->GetYaxis()->SetTitleOffset(1.15);
2929       hists[jk]->SetMarkerStyle(20);
2930       hists[jk]->SetMarkerColor(colors[k]);
2931       hists[jk]->SetLineColor(colors[k]);
2932       if (normalize) {
2933         if (jk == 0)
2934           hists[jk]->DrawNormalized();
2935         else
2936           hists[jk]->DrawNormalized("sames");
2937       } else {
2938         if (jk == 0)
2939           hists[jk]->Draw();
2940         else
2941           hists[jk]->Draw("sames");
2942       }
2943       pad->Update();
2944       TPaveStats* st1 = (TPaveStats*)hists[jk]->GetListOfFunctions()->FindObject("stats");
2945       if (st1 != nullptr) {
2946         double ymin = ymax - dy;
2947         st1->SetLineColor(colors[k]);
2948         st1->SetTextColor(colors[k]);
2949         st1->SetY1NDC(ymin);
2950         st1->SetY2NDC(ymax);
2951         st1->SetX1NDC(0.70);
2952         st1->SetX2NDC(0.90);
2953         ymax = ymin;
2954       }
2955       sprintf(name, "%s", texts[jk].c_str());
2956       legend->AddEntry(hists[jk], name, "lp");
2957     }
2958     legend->Draw("same");
2959     pad->Update();
2960     TPaveText* txt1 = new TPaveText(0.10, 0.905, 0.80, 0.95, "blNDC");
2961     txt1->SetFillColor(0);
2962     char txt[100];
2963     sprintf(txt, "%s", title[type].c_str());
2964     txt1->AddText(txt);
2965     txt1->Draw("same");
2966     /*
2967     TPaveText *txt2 = new TPaveText(0.11,0.825,0.33,0.895,"blNDC");
2968     txt2->SetFillColor(0);
2969     sprintf (txt, "CMS Preliminary");
2970     txt2->AddText(txt);
2971     txt2->Draw("same");
2972     */
2973     pad->Modified();
2974     pad->Update();
2975     if (save > 0) {
2976       sprintf(name, "%s.pdf", pad->GetName());
2977       pad->Print(name);
2978     } else if (save < 0) {
2979       sprintf(name, "%s.C", pad->GetName());
2980       pad->Print(name);
2981     }
2982   }
2983 }
2984 
2985 void PlotPUCorrHists(std::string infile = "corrfac.root",
2986                      std::string prefix = "",
2987                      int drawStatBox = 0,
2988                      bool approve = true,
2989                      int save = 0) {
2990   std::string name1[4] = {"W0", "W1", "W2", "P"};
2991   std::string name2[4] = {"All", "Barrel", "Endcap", ""};
2992   std::string name3[2] = {"", "p = 40:60 GeV"};
2993   std::string name4[2] = {"Loose Isolation", "Tight Isolation"};
2994   std::string xtitle[4] = {"Correction Factor", "Correction Factor", "Correction Factor", "i#eta"};
2995   std::string ytitle[4] = {"Tracks", "Tracks", "Tracks", "Correction Factor"};
2996 
2997   gStyle->SetCanvasBorderMode(0);
2998   gStyle->SetCanvasColor(kWhite);
2999   gStyle->SetPadColor(kWhite);
3000   gStyle->SetFillColor(kWhite);
3001   gStyle->SetOptTitle(0);
3002   gStyle->SetOptFit(0);
3003   if (drawStatBox == 0)
3004     gStyle->SetOptStat(0);
3005   else
3006     gStyle->SetOptStat(1110);
3007 
3008   char name[100], namep[100], title[100];
3009   TFile* file = new TFile(infile.c_str());
3010 
3011   if (file != nullptr) {
3012     for (int i1 = 0; i1 < 4; ++i1) {
3013       for (int i2 = 0; i2 < 2; ++i2) {
3014         for (int i3 = 0; i3 < 2; ++i3) {
3015           sprintf(name, "%s%d%d", name1[i1].c_str(), i2, i3);
3016           if (i2 == 0)
3017             sprintf(title, "%s Tracks Selected with %s", name2[i1].c_str(), name4[i3].c_str());
3018           else
3019             sprintf(title, "%s Tracks Selected with %s (%s)", name2[i1].c_str(), name4[i3].c_str(), name3[i2].c_str());
3020           TH1D* hist1(nullptr);
3021           TProfile* hist2(nullptr);
3022           if (i1 != 3) {
3023             TH1D* hist = (TH1D*)file->FindObjectAny(name);
3024             if (hist != nullptr) {
3025               hist1 = (TH1D*)(hist->Clone());
3026               hist1->GetXaxis()->SetTitleSize(0.040);
3027               hist1->GetXaxis()->SetTitle(xtitle[i1].c_str());
3028               hist1->GetYaxis()->SetTitle(ytitle[i1].c_str());
3029               hist1->GetYaxis()->SetLabelOffset(0.005);
3030               hist1->GetYaxis()->SetLabelSize(0.035);
3031               hist1->GetYaxis()->SetTitleSize(0.040);
3032               hist1->GetYaxis()->SetTitleOffset(1.15);
3033             }
3034           } else {
3035             TProfile* hist = (TProfile*)file->FindObjectAny(name);
3036             if (hist != nullptr) {
3037               hist2 = (TProfile*)(hist->Clone());
3038               hist2->GetXaxis()->SetTitleSize(0.040);
3039               hist2->GetXaxis()->SetTitle(xtitle[i1].c_str());
3040               hist2->GetYaxis()->SetTitle(ytitle[i1].c_str());
3041               hist2->GetYaxis()->SetLabelOffset(0.005);
3042               hist2->GetYaxis()->SetLabelSize(0.035);
3043               hist2->GetYaxis()->SetTitleSize(0.040);
3044               hist2->GetYaxis()->SetTitleOffset(1.15);
3045               //          hist2->GetYaxis()->SetRangeUser(0.0, 1.5);
3046               hist2->SetMarkerStyle(20);
3047             }
3048           }
3049           if ((hist1 != nullptr) || (hist2 != nullptr)) {
3050             sprintf(namep, "c_%s%s", name, prefix.c_str());
3051             TCanvas* pad = new TCanvas(namep, namep, 700, 500);
3052             pad->SetRightMargin(0.10);
3053             pad->SetTopMargin(0.10);
3054             if (hist1 != nullptr) {
3055               pad->SetLogy();
3056               hist1->Draw();
3057               pad->Update();
3058               TPaveStats* st1 = (TPaveStats*)hist1->GetListOfFunctions()->FindObject("stats");
3059               if (st1 != nullptr) {
3060                 st1->SetY1NDC(0.77);
3061                 st1->SetY2NDC(0.90);
3062                 st1->SetX1NDC(0.70);
3063                 st1->SetX2NDC(0.90);
3064               }
3065             } else {
3066               hist2->Draw();
3067               pad->Update();
3068             }
3069             TPaveText* txt1 = new TPaveText(0.10, 0.905, 0.80, 0.95, "blNDC");
3070             txt1->SetFillColor(0);
3071             char txt[100];
3072             sprintf(txt, "%s", title);
3073             txt1->AddText(txt);
3074             txt1->Draw("same");
3075             if (approve) {
3076               double xoff = (i1 == 3) ? 0.11 : 0.22;
3077               TPaveText* txt2 = new TPaveText(xoff, 0.825, xoff + 0.22, 0.895, "blNDC");
3078               txt2->SetFillColor(0);
3079               sprintf(txt, "CMS Preliminary");
3080               txt2->AddText(txt);
3081               txt2->Draw("same");
3082             }
3083             pad->Modified();
3084             pad->Update();
3085             if (save > 0) {
3086               sprintf(name, "%s.pdf", pad->GetName());
3087               pad->Print(name);
3088             } else if (save < 0) {
3089               sprintf(name, "%s.C", pad->GetName());
3090               pad->Print(name);
3091             }
3092           }
3093         }
3094       }
3095     }
3096   }
3097 }
3098 
3099 void PlotHistCorr(const char* infile,
3100                   std::string prefix,
3101                   std::string text0,
3102                   int eta = 0,
3103                   int mode = 1,
3104                   bool drawStatBox = true,
3105                   int save = 0) {
3106   gStyle->SetCanvasBorderMode(0);
3107   gStyle->SetCanvasColor(kWhite);
3108   gStyle->SetPadColor(kWhite);
3109   gStyle->SetFillColor(kWhite);
3110   gStyle->SetOptTitle(0);
3111   if (drawStatBox)
3112     gStyle->SetOptStat(1100);
3113   else
3114     gStyle->SetOptStat(0);
3115 
3116   std::string tags[3] = {"UnNoPU", "UnPU", "Cor"};
3117   std::string text[3] = {"Uncorrected no PU", "Uncorrected PU", "Corrected PU"};
3118   int colors[3] = {1, 4, 2};
3119   int styles[3] = {1, 3, 2};
3120   TFile* file = new TFile(infile);
3121   if (mode < 0 || mode > 2)
3122     mode = 1;
3123   int etamin = (eta == 0) ? -27 : eta;
3124   int etamax = (eta == 0) ? 27 : eta;
3125   for (int ieta = etamin; ieta <= etamax; ++ieta) {
3126     char name[20];
3127     double yh(0.90), dy(0.09);
3128     double yh1 = drawStatBox ? (yh - 3 * dy - 0.01) : (yh - 0.01);
3129     TLegend* legend = new TLegend(0.55, yh1 - 0.15, 0.89, yh1);
3130     legend->SetFillColor(kWhite);
3131     sprintf(name, "c_%sEovp%d", prefix.c_str(), ieta);
3132     TCanvas* pad = new TCanvas(name, name, 700, 500);
3133     pad->SetRightMargin(0.10);
3134     pad->SetTopMargin(0.10);
3135     TH1D* hist[3];
3136     double ymax(0);
3137     for (int k = 0; k < 3; ++k) {
3138       if (k < 2)
3139         sprintf(name, "EovP_ieta%d%s", ieta, tags[k].c_str());
3140       else
3141         sprintf(name, "EovP_ieta%dCor%dPU", ieta, mode);
3142       TH1D* hist1 = (TH1D*)file->FindObjectAny(name);
3143       if (hist1 != nullptr) {
3144         hist[k] = (TH1D*)(hist1->Clone());
3145         ymax = std::max(ymax, (hist1->GetMaximum()));
3146       }
3147     }
3148     int imax = 10 * (2 + int(0.1 * ymax));
3149     for (int k = 0; k < 3; ++k) {
3150       hist[k]->GetYaxis()->SetLabelOffset(0.005);
3151       hist[k]->GetYaxis()->SetTitleOffset(1.20);
3152       hist[k]->GetXaxis()->SetTitle("E/p");
3153       hist[k]->GetYaxis()->SetTitle("Tracks");
3154       hist[k]->SetLineColor(colors[k]);
3155       hist[k]->SetLineStyle(styles[k]);
3156       hist[k]->GetYaxis()->SetRangeUser(0.0, imax);
3157       if (k == 0)
3158         hist[k]->Draw();
3159       else
3160         hist[k]->Draw("sames");
3161       legend->AddEntry(hist[k], text[k].c_str(), "lp");
3162       pad->Update();
3163       if (drawStatBox) {
3164         TPaveStats* st1 = (TPaveStats*)hist[k]->GetListOfFunctions()->FindObject("stats");
3165         if (st1 != nullptr) {
3166           st1->SetLineColor(colors[k]);
3167           st1->SetTextColor(colors[k]);
3168           st1->SetY1NDC(yh - dy);
3169           st1->SetY2NDC(yh);
3170           st1->SetX1NDC(0.70);
3171           st1->SetX2NDC(0.90);
3172           yh -= dy;
3173         }
3174       }
3175     }
3176     pad->Update();
3177     legend->Draw("same");
3178     pad->Update();
3179     TPaveText* txt1 = new TPaveText(0.10, 0.905, 0.80, 0.95, "blNDC");
3180     txt1->SetFillColor(0);
3181     char title[100];
3182     sprintf(title, "%s for i#eta = %d", text0.c_str(), ieta);
3183     txt1->AddText(title);
3184     txt1->Draw("same");
3185     pad->Modified();
3186     pad->Update();
3187     if (save > 0) {
3188       sprintf(name, "%s.pdf", pad->GetName());
3189       pad->Print(name);
3190     } else if (save < 0) {
3191       sprintf(name, "%s.C", pad->GetName());
3192       pad->Print(name);
3193     }
3194   }
3195 }
3196 
3197 void PlotPropertyHist(const char* infile,
3198                       std::string prefix,
3199                       std::string text,
3200                       int etaMax = 25,
3201                       double lumi = 0,
3202                       double ener = 13.0,
3203                       bool dataMC = false,
3204                       bool drawStatBox = true,
3205                       int save = 0) {
3206   std::string name0[3] = {"energyE2", "energyH2", "energyP2"};
3207   std::string title0[3] = {"Energy in ECAL", "Energy in HCAL", "Track Momentum"};
3208   std::string xtitl0[3] = {"Energy (GeV)", "Energy (GeV)", "p (GeV)"};
3209   std::string name1[5] = {"eta02", "eta12", "eta22", "eta32", "eta42"};
3210   std::string name10[5] = {"eta0", "eta1", "eta2", "eta3", "eta4"};
3211   std::string xtitl1 = "i#eta";
3212   std::string name2[5] = {"p0", "p1", "p2", "p3", "p4"};
3213   std::string xtitl2 = "p (GeV)";
3214   std::string title1[5] = {"Tracks with p=40:60 GeV",
3215                            "Good Quality Tracks with p=40:60 GeV",
3216                            "Selected Tracks with p=40:60 GeV",
3217                            "Isolated Tracks with p=40:60 GeV",
3218                            "Isolated Tracks with p=40:60 GeV and MIPS in ECAL"};
3219   std::string title2[5] = {
3220       "All Tracks", "Good Quality Tracks", "Selected Tracks", "Isolated Tracks", "Isolated Tracks with MIPS in ECAL"};
3221   std::string ytitle = "Tracks";
3222 
3223   gStyle->SetCanvasBorderMode(0);
3224   gStyle->SetCanvasColor(kWhite);
3225   gStyle->SetPadColor(kWhite);
3226   gStyle->SetFillColor(kWhite);
3227   gStyle->SetOptTitle(0);
3228   if (drawStatBox)
3229     gStyle->SetOptStat(1110);
3230   else
3231     gStyle->SetOptStat(0);
3232   gStyle->SetOptFit(0);
3233 
3234   TFile* file = new TFile(infile);
3235   char name[100], namep[100];
3236   for (int k = 1; k <= etaMax; ++k) {
3237     for (int j = 0; j < 3; ++j) {
3238       sprintf(name, "%s%s%d", prefix.c_str(), name0[j].c_str(), k);
3239       TH1D* hist1 = (TH1D*)file->FindObjectAny(name);
3240       if (hist1 != nullptr) {
3241         TH1D* hist = (TH1D*)(hist1->Clone());
3242         double ymin(0.90);
3243         sprintf(namep, "c_%s", name);
3244         TCanvas* pad = new TCanvas(namep, namep, 700, 500);
3245         pad->SetLogy();
3246         pad->SetRightMargin(0.10);
3247         pad->SetTopMargin(0.10);
3248         hist->GetXaxis()->SetTitleSize(0.04);
3249         hist->GetXaxis()->SetTitle(xtitl0[j].c_str());
3250         hist->GetYaxis()->SetTitle(ytitle.c_str());
3251         hist->GetYaxis()->SetLabelOffset(0.005);
3252         hist->GetYaxis()->SetTitleSize(0.04);
3253         hist->GetYaxis()->SetLabelSize(0.035);
3254         hist->GetYaxis()->SetTitleOffset(1.10);
3255         hist->SetMarkerStyle(20);
3256         hist->Draw();
3257         pad->Update();
3258         TPaveStats* st1 = (TPaveStats*)hist->GetListOfFunctions()->FindObject("stats");
3259         if (st1 != nullptr) {
3260           st1->SetY1NDC(0.78);
3261           st1->SetY2NDC(0.90);
3262           st1->SetX1NDC(0.65);
3263           st1->SetX2NDC(0.90);
3264         }
3265         pad->Update();
3266         double ymx(0.96), xmi(0.35), xmx(0.95);
3267         char txt[100];
3268         if (lumi > 0.1) {
3269           ymx = ymin - 0.005;
3270           xmi = 0.45;
3271           TPaveText* txt0 = new TPaveText(0.65, 0.91, 0.90, 0.96, "blNDC");
3272           txt0->SetFillColor(0);
3273           sprintf(txt, "%4.1f TeV %5.1f fb^{-1}", ener, lumi);
3274           txt0->AddText(txt);
3275           txt0->Draw("same");
3276         }
3277         double ymi = ymx - 0.05;
3278         TPaveText* txt1 = new TPaveText(xmi, ymi, xmx, ymx, "blNDC");
3279         txt1->SetFillColor(0);
3280         if (text == "") {
3281           sprintf(txt, "%s", title0[j].c_str());
3282         } else {
3283           sprintf(txt, "%s (%s)", title0[j].c_str(), text.c_str());
3284         }
3285         txt1->AddText(txt);
3286         txt1->Draw("same");
3287         double xmax = (dataMC) ? 0.24 : 0.35;
3288         ymi = 0.91;
3289         ymx = ymi + 0.05;
3290         TPaveText* txt2 = new TPaveText(0.02, ymi, xmax, ymx, "blNDC");
3291         txt2->SetFillColor(0);
3292         if (dataMC)
3293           sprintf(txt, "CMS Preliminary");
3294         else
3295           sprintf(txt, "CMS Simulation Preliminary");
3296         txt2->AddText(txt);
3297         txt2->Draw("same");
3298         pad->Modified();
3299         pad->Update();
3300         if (save > 0) {
3301           sprintf(name, "%s.pdf", pad->GetName());
3302           pad->Print(name);
3303         } else if (save < 0) {
3304           sprintf(name, "%s.C", pad->GetName());
3305           pad->Print(name);
3306         }
3307       }
3308     }
3309   }
3310 
3311   for (int k = 0; k < 3; ++k) {
3312     for (int j = 0; j < 5; ++j) {
3313       if (k == 0)
3314         sprintf(name, "%s%s", prefix.c_str(), name1[j].c_str());
3315       else if (k == 1)
3316         sprintf(name, "%s%s", prefix.c_str(), name10[j].c_str());
3317       else
3318         sprintf(name, "%s%s", prefix.c_str(), name2[j].c_str());
3319       TH1D* hist1 = (TH1D*)file->FindObjectAny(name);
3320       if (hist1 != nullptr) {
3321         TH1D* hist = (TH1D*)(hist1->Clone());
3322         double ymin(0.90);
3323         sprintf(namep, "c_%s", name);
3324         TCanvas* pad = new TCanvas(namep, namep, 700, 500);
3325         pad->SetLogy();
3326         pad->SetRightMargin(0.10);
3327         pad->SetTopMargin(0.10);
3328         hist->GetXaxis()->SetTitleSize(0.04);
3329         if (k <= 1)
3330           hist->GetXaxis()->SetTitle(xtitl1.c_str());
3331         else
3332           hist->GetXaxis()->SetTitle(xtitl2.c_str());
3333         hist->GetYaxis()->SetTitle(ytitle.c_str());
3334         hist->GetYaxis()->SetLabelOffset(0.005);
3335         hist->GetYaxis()->SetTitleSize(0.04);
3336         hist->GetYaxis()->SetLabelSize(0.035);
3337         hist->GetYaxis()->SetTitleOffset(1.10);
3338         hist->SetMarkerStyle(20);
3339         hist->Draw();
3340         pad->Update();
3341         TPaveStats* st1 = (TPaveStats*)hist->GetListOfFunctions()->FindObject("stats");
3342         if (st1 != nullptr) {
3343           st1->SetY1NDC(0.78);
3344           st1->SetY2NDC(0.90);
3345           st1->SetX1NDC(0.65);
3346           st1->SetX2NDC(0.90);
3347         }
3348         pad->Update();
3349         double ymx(0.96), xmi(0.35), xmx(0.95);
3350         char txt[100];
3351         if (lumi > 0.1) {
3352           ymx = ymin - 0.005;
3353           xmi = 0.45;
3354           TPaveText* txt0 = new TPaveText(0.65, 0.91, 0.90, 0.96, "blNDC");
3355           txt0->SetFillColor(0);
3356           sprintf(txt, "%4.1f TeV %5.1f fb^{-1}", ener, lumi);
3357           txt0->AddText(txt);
3358           txt0->Draw("same");
3359         }
3360         double ymi = ymx - 0.05;
3361         TPaveText* txt1 = new TPaveText(xmi, ymi, xmx, ymx, "blNDC");
3362         txt1->SetFillColor(0);
3363         if (text == "") {
3364           if (k == 0)
3365             sprintf(txt, "%s", title1[j].c_str());
3366           else
3367             sprintf(txt, "%s", title2[j].c_str());
3368         } else {
3369           if (k == 0)
3370             sprintf(txt, "%s (%s)", title1[j].c_str(), text.c_str());
3371           else
3372             sprintf(txt, "%s (%s)", title2[j].c_str(), text.c_str());
3373         }
3374         txt1->AddText(txt);
3375         txt1->Draw("same");
3376         double xmax = (dataMC) ? 0.24 : 0.35;
3377         ymi = 0.91;
3378         ymx = ymi + 0.05;
3379         TPaveText* txt2 = new TPaveText(0.02, ymi, xmax, ymx, "blNDC");
3380         txt2->SetFillColor(0);
3381         if (dataMC)
3382           sprintf(txt, "CMS Preliminary");
3383         else
3384           sprintf(txt, "CMS Simulation Preliminary");
3385         txt2->AddText(txt);
3386         txt2->Draw("same");
3387         pad->Modified();
3388         pad->Update();
3389         if (save > 0) {
3390           sprintf(name, "%s.pdf", pad->GetName());
3391           pad->Print(name);
3392         } else if (save < 0) {
3393           sprintf(name, "%s.C", pad->GetName());
3394           pad->Print(name);
3395         }
3396       }
3397     }
3398   }
3399 }
3400 
3401 void PlotMeanError(const std::string infilest, int reg = 3, bool resol = false, int save = 0, bool debug = false) {
3402   bool ok(false);
3403   const int ntypmx = 3;
3404   const int nregmx = 4;
3405   if (reg < 0 || reg >= nregmx)
3406     reg = nregmx - 1;
3407   int nEner(0), nType(0), nPts(0);
3408   std::vector<double> energy[ntypmx], denergy[ntypmx], value[ntypmx], dvalue[ntypmx];
3409   // First read the data
3410   std::ifstream fInput(infilest.c_str());
3411   if (!fInput.good()) {
3412     std::cout << "Cannot open file " << infilest << std::endl;
3413   } else {
3414     ok = true;
3415     fInput >> nEner >> nType >> nPts;
3416     int nmax = nEner * nType;
3417     int type, elow, ehigh;
3418     double v1[4], e1[4], v2[4], e2[4];
3419     for (int n = 0; n < nmax; ++n) {
3420       fInput >> type >> elow >> ehigh;
3421       fInput >> v1[0] >> e1[0] >> v1[1] >> e1[1] >> v1[2] >> e1[2] >> v1[3] >> e1[3];
3422       fInput >> v2[0] >> e2[0] >> v2[1] >> e2[1] >> v2[2] >> e2[2] >> v2[3] >> e2[3];
3423       double ener = 0.5 * (ehigh + elow);
3424       double dene = 0.5 * (ehigh - elow);
3425       energy[type].push_back(ener);
3426       denergy[type].push_back(dene);
3427       if (resol) {
3428         value[type].push_back(v2[reg]);
3429         dvalue[type].push_back(e2[reg]);
3430       } else {
3431         value[type].push_back(v1[reg]);
3432         dvalue[type].push_back(e1[reg]);
3433       }
3434     }
3435     fInput.close();
3436     std::cout << "Reads " << (nmax + 1) << " cards from " << infilest << " with measurements for " << nEner
3437               << " energies and " << nType << " types" << std::endl;
3438     if (debug) {
3439       for (int n = 0; n < nType; ++n) {
3440         std::cout << "Type " << n << " with " << energy[n].size() << " points\n";
3441         for (unsigned int k = 0; k < energy[n].size(); ++k)
3442           std::cout << " [" << k << "] " << energy[n][k] << " +- " << denergy[n][k] << " Value " << value[n][k]
3443                     << " +- " << dvalue[n][k] << std::endl;
3444       }
3445     }
3446   }
3447 
3448   // Now the plots
3449   if (ok) {
3450     int mvsres = (resol) ? 1 : 0;
3451     std::string names[2] = {"Mean", "Resol"};
3452     std::string namet[nregmx] = {"Barrel", "Transition", "Endcap", "Combined"};
3453     char cname[100];
3454     sprintf(cname, "c_%s%s", names[mvsres].c_str(), namet[reg].c_str());
3455     int color[ntypmx] = {2, 4, 1};
3456     int mtype[ntypmx] = {20, 21, 22};
3457     double ymin[2] = {0.65, 0.10};
3458     double ymax[2] = {1.30, 0.50};
3459     gStyle->SetCanvasBorderMode(0);
3460     gStyle->SetCanvasColor(kWhite);
3461     gStyle->SetPadColor(kWhite);
3462     gStyle->SetFillColor(kWhite);
3463     gStyle->SetOptTitle(kFALSE);
3464     gStyle->SetPadBorderMode(0);
3465     gStyle->SetCanvasBorderMode(0);
3466     gStyle->SetOptStat(0);
3467     TCanvas* canvas = new TCanvas(cname, cname, 500, 500);
3468     canvas->SetTopMargin(0.05);
3469     canvas->SetBottomMargin(0.14);
3470     canvas->SetLeftMargin(0.15);
3471     canvas->SetRightMargin(0.10);
3472     TH1F* vFrame = canvas->DrawFrame(0.0, ymin[mvsres], 120.0, ymax[mvsres]);
3473     vFrame->GetXaxis()->SetRangeUser(0.0, 120.0);
3474     vFrame->GetYaxis()->SetRangeUser(ymin[mvsres], ymax[mvsres]);
3475     vFrame->GetXaxis()->SetLabelSize(0.04);
3476     vFrame->GetYaxis()->SetLabelSize(0.04);
3477     vFrame->GetXaxis()->SetTitleSize(0.04);
3478     vFrame->GetYaxis()->SetTitleSize(0.04);
3479     vFrame->GetXaxis()->SetTitleOffset(1.2);
3480     vFrame->GetYaxis()->SetTitleOffset(1.6);
3481     vFrame->GetXaxis()->SetLabelOffset(0.0);
3482     vFrame->GetXaxis()->SetTitle("p_{Beam} (GeV/c)");
3483     if (resol) {
3484       vFrame->GetYaxis()->SetTitle("Width(E_{HCAL}/(p-E_{ECAL}))");
3485     } else {
3486       vFrame->GetYaxis()->SetTitle("MPV(E_{HCAL}/(p-E_{ECAL}))");
3487     }
3488     TLegend* legend = new TLegend(0.70, 0.80, 0.90, 0.94);
3489     legend->SetFillColor(kWhite);
3490     std::string nameg[ntypmx] = {"MAHI", "M0", "M2"};
3491     for (int n = 0; n < nType; ++n) {
3492       unsigned int nmax0 = energy[n].size();
3493       double mom[nmax0], dmom[nmax0], mean[nmax0], dmean[nmax0];
3494       for (unsigned int k = 0; k < nmax0; ++k) {
3495         mom[k] = energy[n][k];
3496         dmom[k] = denergy[n][k];
3497         mean[k] = value[n][k];
3498         dmean[k] = dvalue[n][k];
3499       }
3500       TGraphErrors* graph = new TGraphErrors(nmax0, mom, mean, dmom, dmean);
3501       graph->SetMarkerStyle(mtype[n]);
3502       graph->SetMarkerColor(color[n]);
3503       graph->SetMarkerSize(1.4);
3504       graph->SetLineColor(color[n]);
3505       graph->SetLineWidth(2);
3506       sprintf(cname, "%s", nameg[n].c_str());
3507       legend->AddEntry(graph, cname, "lp");
3508       graph->Draw("P");
3509     }
3510     legend->Draw("same");
3511     std::string regions[nregmx] = {"20118B Barrel", "2018B Transition", "2018B Endcap", "2018B"};
3512     sprintf(cname, "%s", regions[reg].c_str());
3513     TPaveText* txt0 = new TPaveText(0.16, 0.90, 0.40, 0.94, "blNDC");
3514     txt0->SetFillColor(0);
3515     txt0->AddText(cname);
3516     txt0->Draw("same");
3517     TPaveText* txt1 = new TPaveText(0.15, 0.95, 0.40, 0.99, "blNDC");
3518     txt1->SetFillColor(0);
3519     sprintf(cname, "CMS Preliminary");
3520     txt1->AddText(cname);
3521     txt1->Draw("same");
3522     canvas->Update();
3523     if (save > 0) {
3524       sprintf(cname, "%s.pdf", canvas->GetName());
3525       canvas->Print(cname);
3526     } else if (save < 0) {
3527       sprintf(cname, "%s.C", canvas->GetName());
3528       canvas->Print(cname);
3529     }
3530   }
3531 }