Back to home page

Project CMSSW displayed by LXR

 
 

    


File indexing completed on 2024-07-16 02:42:25

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 //             Same as FitHistExtend with summary of each ieta/depth
0016 //  FitHistExtended2(infile, outfile, prefix, numb, append, iname, debug);
0017 //      Defaults: numb=50, append=true, iname=3, debug=false
0018 //
0019 //             For RBX dependence in sets of histograms from CalibMonitor
0020 //  FitHistRBX(infile, outfile, prefix, append, iname);
0021 //      Defaults: append=true, iname=2
0022 //
0023 //             For plotting stored histograms from FitHist's
0024 //  PlotHist(infile, prefix, text, modePlot, kopt, lumi, ener, isRealData,
0025 //           drawStatBox, save);
0026 //      Defaults: modePlot=4, kopt=100, lumi=0, ener=13, isRealData=false,
0027 //                drawStatBox=true, save=0
0028 //
0029 //             For plotting histograms corresponding to individual ieta's
0030 //  PlotHistEta(infile, prefix, text, iene, numb, ieta, lumi, ener, isRealData,
0031 //           drawStatBox, save);
0032 //      Defaults iene=3, numb=50, ieta=0, lumi=0, ener=13.0, isRealData=false,
0033 //                drawStatBox=true, save=0
0034 //
0035 //             For plotting several histograms in the same plot
0036 //             (fits to different data sets for example)
0037 //  PlotHists(infile, prefix, text, drawStatBox, save)
0038 //      Defaults: drawStatBox=true; save=0;
0039 //      Note prefix is common part for all histograms
0040 //
0041 //             For plotting on the same canvas plots with different
0042 //             prefixes residing in the same file with approrprate text
0043 //  PlotTwoHists(infile, prefix1, text1, prefix2, text2, text0, type, iname,
0044 //               lumi, ener, drawStatBox, save);
0045 //      Defaults: type=0; iname=2; lumi=0; ener=13; drawStatBox=true;
0046 //                save=0;
0047 //      Note prefixN, textN have the same meaning as prefix and text for set N
0048 //           text0 is the text for general title added within ()
0049 //           type=0 plots response distributions and MPV of response vs ieta
0050 //               =1 plots MPV of response vs RBX #
0051 //
0052 //             For plotting stored histograms from CalibTree
0053 //  PlotFiveHists(infile, text0, prefix0, type, iname, drawStatBox, normalize,
0054 //                save, prefix1, text1, prefix2, text2, prefix3, text3,
0055 //                 prefix4, text4, prefix5, text5);
0056 //      Defaults: type=0; iname=0; drawStatBox=true; normalize=false;
0057 //                save=0; prefixN=""; textN=""; (for N > 0)
0058 //      Note prefixN, textN have the same meaning as prefix and text for set N
0059 //           text0 is the text for general title added within ()
0060 //           prefix0 is the tag attached to the canvas name
0061 //           type has the same meaning as in PlotTwoHists
0062 //
0063 //  PlotHistCorrResults(infile, text, prefixF, save);
0064 //      Defaults: save=0
0065 //
0066 //             For plotting correction factors
0067 //  PlotHistCorrFactor(infile, text, prefixF, scale, nmin, isRealData,
0068 //                     drawStatBox, iformat, save);
0069 //      Defaults: isRealData=true, drwaStatBox=false, nmin=100, iformat=0,
0070 //                save=0
0071 //
0072 //             For plotting correction factors for a sigle depth
0073 //  PlotHistCorrFactor(infile, text, depth, prefixF, scale, nmin, isRealData,
0074 //                     drawStatBox, iformat, save);
0075 //      Defaults: isRealData=true, drwaStatBox=false, nmin=100, iformat=0,
0076 //                save=0
0077 //
0078 //             For plotting (fractional) asymmetry in the correction factors
0079 //
0080 //  PlotHistCorrAsymmetry(infile, text, prefixF, iformat, save);
0081 //      Defaults: prefixF="", iformat=0, save=0
0082 //
0083 //             For plotting correction factors from upto 5 different runs
0084 //             on the same canvas
0085 //
0086 //  PlotHistCorrFactors(infile1, text1, infile2, text2, infile3, text3,
0087 //                      infile4, text4, infile5, text5, prefixF, ratio,
0088 //                      drawStatBox, nmin, isRealData, year, iformat, save)
0089 //      Defaults: ratio=false, drawStatBox=true, nmin=100, isRealData=false,
0090 //                year=2018, iformat=0, save=0
0091 //
0092 //  PlotHistCorr2Factors(infile1, text1, infile2, text2, depth, prefixF, ratio,
0093 //                      drawStatBox, nmin, isRealData, year, iformat, save)
0094 //      Defaults: ratio=true, drawStatBox=false, nmin=100, isRealData=true,
0095 //                year=2023, iformat=0, save=0
0096 //
0097 //  PlotHistCorrDFactors(infile1, text1, infile2, text2, infile3, text3,
0098 //                      infile4, text4, infile5, text5, depth, prefixF, ratio,
0099 //                      drawStatBox, nmin, isRealData, year, iformat, save)
0100 //      Defaults: ratio=true, drawStatBox=false, nmin=100, isRealData=true,
0101 //                year=2024, iformat=0, save=0
0102 //
0103 //             For plotting correction factors including systematics
0104 //  PlotHistCorrSys(infilec, conds, text, save)
0105 //      Defaults: save=0
0106 //
0107 //             For plotting uncertainties in correction factors with decreasing
0108 //             integrated luminpsoties starting from *lumi*
0109 //  PlotHistCorrLumis(infilec, conds, lumi, save)
0110 //      Defaults: save=0
0111 //
0112 //             For plotting correlation of correction factors
0113 //  PlotHistCorrRel(infile1, infile2, text1, text2, iformat1, iformat2, save)
0114 //      Defaults: iformat1=0, iformat2=0, save=0
0115 //
0116 //             For plotting difference of correction factors for a given depth
0117 //  PlotHistCorrDepth(infile1, infile2, text1, text2, depth, iformat1, iformat2,
0118 //                    save)
0119 //      Defaults: iformat1=0, iformat2=0, save=0
0120 //
0121 //             For plotting four histograms
0122 //  PlotFourHists(infile, prefix0, type, drawStatBox, normalize, save, prefix1,
0123 //                text1, prefix2, text2, prefix3, text3, prefix4, text4)
0124 //      Defaults: type=0, drawStatBox=0, normalize=false, save=0,
0125 //                prefixN="", textN=""
0126 //
0127 //            For plotting PU corrected histograms (o/p of CalibPlotCombine)
0128 //  PlotPUCorrHists(infile, prefix drawStatBox, approve, save)
0129 //      Defaults: infile = "corrfac.root", prefix = "", drawStatBox = 0,
0130 //                approve = true, save = 0
0131 //
0132 //             For plotting histograms obtained from fits to PU correction
0133 //             (o/p of CalibFitPU) for a given ieta using 2D/profile/Graphs
0134 //  PlotHistCorr(infile, prefix, text, eta, mode, drawStatBox, save)
0135 //      Defaults eta = 0 (all ieta values), mode = 1 (profile histograms),
0136 //               drawStatBox = true, save = 0
0137 //
0138 //             For plotting histograms created by CalibPlotProperties
0139 //  PlotPropertyHist(infile, prefix, text, etaMax, lumi, ener, isRealData,
0140 //           drawStatBox, save)
0141 //      Defaults etaMax = 25 (draws for eta = 1 .. etaMax), lumi = 0,
0142 //               ener = 13.0, isRealData = false,  drawStatBox = true, save = 0
0143 //
0144 //            For plotting mean response and resolution as a function of
0145 //            particle momentum
0146 //  PlotMeanError(infilest, region, resol, save, debug)
0147 //      Defaults region = 3 (overall), resol = false (response), save = 0,
0148 //               debug = false
0149 //      Format of the input file:
0150 //       # of energy points, # of types, # of regions
0151 //       Then for each type, energy point
0152 //         Type, lower and higher edge of momentum
0153 //         Mean response and its error for the 4 regions
0154 //         Width of response and uts error for the 4 regions
0155 //
0156 //            For plotting depth dependent correction factors from muon study
0157 //  PlotDepthCorrFactor(infile, text, prefix, isRealData, drawStatBox, save)
0158 //      Defaults prefix = "", isRealData = true, drawStatBox = true, save = 0
0159 //      Format for the input file: ieta and correcrion factor with its
0160 //             uncertainty for each depth
0161 //
0162 //            For plotting ratio of correction factors as defined in a file
0163 //            give by infileX for 2 depths (depth1, depth2) as a function of
0164 //            ieta obaned from 2 sources of data (defined by text1 and text2)
0165 //  PlotHistCorrRatio(infile1, text1, infile2, text2, depth1, depth2, prefix,
0166 //                    text0, etaMin, etaMax, doFit, isRealData, year, iformat,
0167 //                    save)
0168 //      Defaults etaMin = -1, etaMax = -1, doFit = true, isRealData = true,
0169 //               year = 2022, iformat = 0, save = 0
0170 //      text0 is a general description common to both sets of corr factors
0171 //      etaMin < 0 and etaMax > 0 will take ieta range from -etaMax to +etaMax;
0172 //      etaMin > 0 will select ieta's where |ieta| is greater than etaMin
0173 //      with the plot either between -etaMax to etaMax if etaMax > 0 otherwise
0174 //      determined from data files;
0175 //      doFit determines if a Pol0 fit is to be done
0176 //
0177 //  where:
0178 //  infile   (std::string)  = Name of the input ROOT file
0179 //  outfile  (std::string)  = Name of the output ROOT file
0180 //  prefix   (std::string)  = Prefix for the histogram names
0181 //  mode     (int)          = Flag to check which set of histograms to be
0182 //                            done. It has the format lthdo where each of
0183 //                            l, t,h,d,o can have a value 0 or 1 to select
0184 //                            or deselect. l,t,h,d,o for momentum range
0185 //                            60-100, 30-40, all, 20-30, 40-60 Gev (11111)
0186 //  type     (int)          = defines eta binning type (see CalibMonitor)
0187 //  append   (bool)         = Open the output in Update/Recreate mode (True)
0188 //  fiteta   (bool)         = fit the eta dependence with pol0
0189 //  iname    (int)          = choose the momentum bin (3: 40-60 GeV)
0190 //  saveAll  (bool)         = Flag to save intermediate plots (False)
0191 //  numb     (int)          = Number of eta bins (42 for -21:21)
0192 //  text     (std::string)  = Extra text to be put in the text title
0193 //  modePlot (int)          = Flag to plot E/p distribution (0);
0194 //                            <E/p> as a function of ieta (1);
0195 //                            <E/p> as a function of distance from L1 (2);
0196 //                            <E/p> as a function of number of vertex (3);
0197 //                            E/p for barrel, endcap and transition (4)
0198 //  kopt     (int)          = Option in format "hdo" where each of d, o can
0199 //                            have a value of 0 or 1 to select or deselect.
0200 //                            o>0 to carry out pol0 fit, o>1 to restrict
0201 //                            fit region between -20 & 20; d=1 to show grid;
0202 //                            h=0,1 to show plots with 2- or 1-Gaussian fit
0203 //  ieta     (int)          = specific ieta histogram to be plotted; if 0
0204 //                            histograms for all ieta's from -numb/2 to numb/2
0205 //                            will be plotted
0206 //  lumi     (double)       = Integrated luminosity of the dataset used which
0207 //                            needs to be drawn on the top of the canvas
0208 //                            along with CM energy (if lumi > 0)
0209 //  ener     (double)       = CM energy of the dataset used
0210 //  save     (int)          = if > 0 it saves the canvas as a pdf file; or
0211 //                            if < 0 it saves the canvas as a C file
0212 //  normalize(bool)         = if the histograms to be scaled to get
0213 //                            normalization to 1
0214 //  prefixF  (string)       = string to be included to the pad name
0215 //  infileN  (char*)        = name of the correction file and the corresponding
0216 //  textN    (string)         texts (if blank ignored)
0217 //  nmin     (int)          = minimum # of #ieta points needed to show the
0218 //                            fitted line
0219 //  scale    (double)       = constant scale factor applied to the factors
0220 //  ratio    (bool)         = set to show the ratio plot (false)
0221 //  drawStatBox (bool)      = set to show the statistical box (true)
0222 //  year     (int)          = Year of data taking (applicable to Data)
0223 //  infilc   (string)       = prefix of the file names of correction factors
0224 //                            (assumes file name would be the prefix followed
0225 //                            by _condX.txt where X=0 for the default version
0226 //                            and 1..conds for the variations)
0227 //  conds    (int)          = number of variations in estimating systematic
0228 //                            checks
0229 //  infilest (string)       = input file name containing the responses and
0230 //                            resolutions for barrel, transition, endcap,
0231 //                            overall regions at 5 energies using 3 methods
0232 //  region   (int)         = region to be selected: 0 = barrel, 1 = transition,
0233 //                           2 = endcap, 3 = overall (3)
0234 //  resol    (bool)        = parameter to be plotted: true = resolution,
0235 //                           false = response (false)
0236 //  iformat  (int)         = flag if it is created by standard (0) or by
0237 //                           Marina's (1) code
0238 //////////////////////////////////////////////////////////////////////////////
0239 
0240 #include <TCanvas.h>
0241 #include <TChain.h>
0242 #include <TProfile.h>
0243 #include <TF1.h>
0244 #include <TFile.h>
0245 #include <TFitResult.h>
0246 #include <TFitResultPtr.h>
0247 #include <TH1D.h>
0248 #include <TLegend.h>
0249 #include <TLine.h>
0250 #include <TGraph.h>
0251 #include <TGraphErrors.h>
0252 #include <TGraphAsymmErrors.h>
0253 #include <TMath.h>
0254 #include <TPaveStats.h>
0255 #include <TPaveText.h>
0256 #include <TROOT.h>
0257 #include <TStyle.h>
0258 #include <cstdlib>
0259 #include <fstream>
0260 #include <iomanip>
0261 #include <iostream>
0262 #include <string>
0263 #include <vector>
0264 
0265 #include "CalibCorr.C"
0266 
0267 const double fitrangeFactor = 1.5;
0268 
0269 struct cfactors {
0270   int ieta, depth;
0271   double corrf, dcorr;
0272   cfactors(int ie = 0, int dp = 0, double cor = 1, double dc = 0) : ieta(ie), depth(dp), corrf(cor), dcorr(dc){};
0273 };
0274 
0275 struct results {
0276   double mean, errmean, width, errwidth;
0277   results(double v1 = 0, double er1 = 0, double v2 = 0, double er2 = 0)
0278       : mean(v1), errmean(er1), width(v2), errwidth(er2){};
0279 };
0280 
0281 std::pair<double, double> GetMean(TH1D* hist, double xmin, double xmax, double& rms) {
0282   double mean(0), err(0), wt(0);
0283   rms = 0;
0284   for (int i = 1; i <= hist->GetNbinsX(); ++i) {
0285     double xlow = hist->GetBinLowEdge(i);
0286     double xhigh = xlow + hist->GetBinWidth(i);
0287     if ((xlow >= xmin) && (xhigh <= xmax)) {
0288       double cont = hist->GetBinContent(i);
0289       double valu = 0.5 * (xlow + xhigh);
0290       wt += cont;
0291       mean += (valu * cont);
0292       rms += (valu * valu * cont);
0293     }
0294   }
0295   if (wt > 0) {
0296     mean /= wt;
0297     rms /= wt;
0298     err = std::sqrt((rms - mean * mean) / wt);
0299   }
0300   return std::pair<double, double>(mean, err);
0301 }
0302 
0303 std::pair<double, double> GetWidth(TH1D* hist, double xmin, double xmax) {
0304   double mean(0), mom2(0), rms(0), err(0), wt(0);
0305   for (int i = 1; i <= hist->GetNbinsX(); ++i) {
0306     double xlow = hist->GetBinLowEdge(i);
0307     double xhigh = xlow + hist->GetBinWidth(i);
0308     if ((xlow >= xmin) && (xhigh <= xmax)) {
0309       double cont = hist->GetBinContent(i);
0310       double valu = 0.5 * (xlow + xhigh);
0311       wt += cont;
0312       mean += (valu * cont);
0313       mom2 += (valu * valu * cont);
0314     }
0315   }
0316   if (wt > 0) {
0317     mean /= wt;
0318     mom2 /= wt;
0319     rms = std::sqrt(mom2 - mean * mean);
0320     err = rms / std::sqrt(2 * wt);
0321   }
0322   return std::pair<double, double>(rms, err);
0323 }
0324 
0325 Double_t langaufun(Double_t* x, Double_t* par) {
0326   //Fit parameters:
0327   //par[0]=Most Probable (MP, location) parameter of Landau density
0328   //par[1]=Total area (integral -inf to inf, normalization constant)
0329   //par[2]=Width (sigma) of convoluted Gaussian function
0330   //
0331   //In the Landau distribution (represented by the CERNLIB approximation),
0332   //the maximum is located at x=-0.22278298 with the location parameter=0.
0333   //This shift is corrected within this function, so that the actual
0334   //maximum is identical to the MP parameter.
0335 
0336   // Numeric constants
0337   Double_t invsq2pi = 0.3989422804014;  // (2 pi)^(-1/2)
0338   Double_t mpshift = -0.22278298;       // Landau maximum location
0339 
0340   // Control constants
0341   Double_t np = 100.0;  // number of convolution steps
0342   Double_t sc = 5.0;    // convolution extends to +-sc Gaussian sigmas
0343 
0344   // Variables
0345   Double_t xx;
0346   Double_t mpc;
0347   Double_t fland;
0348   Double_t sum = 0.0;
0349   Double_t xlow, xupp;
0350   Double_t step;
0351   Double_t i;
0352   Double_t par0 = 0.2;
0353 
0354   // MP shift correction
0355   mpc = par[0] - mpshift * par0 * par[0];
0356 
0357   // Range of convolution integral
0358   xlow = x[0] - sc * par[2];
0359   xupp = x[0] + sc * par[2];
0360 
0361   step = (xupp - xlow) / np;
0362 
0363   // Convolution integral of Landau and Gaussian by sum
0364   for (i = 1.0; i <= np / 2; i++) {
0365     xx = xlow + (i - .5) * step;
0366     fland = TMath::Landau(xx, mpc, par0 * par[0], kTRUE);  // / par[0];
0367     sum += fland * TMath::Gaus(x[0], xx, par[2]);
0368 
0369     xx = xupp - (i - .5) * step;
0370     fland = TMath::Landau(xx, mpc, par0 * par[0], kTRUE);  // / par[0];
0371     sum += fland * TMath::Gaus(x[0], xx, par[2]);
0372   }
0373 
0374   return (par[1] * step * sum * invsq2pi / par[2]);
0375 }
0376 
0377 Double_t doubleGauss(Double_t* x, Double_t* par) {
0378   double x1 = x[0] - par[1];
0379   double sig1 = par[2];
0380   double x2 = x[0] - par[4];
0381   double sig2 = par[5];
0382   double yval =
0383       (par[0] * std::exp(-0.5 * (x1 / sig1) * (x1 / sig1)) + par[3] * std::exp(-0.5 * (x2 / sig2) * (x2 / sig2)));
0384   return yval;
0385 }
0386 
0387 TFitResultPtr functionFit(TH1D* hist, double* fitrange, double* startvalues, double* parlimitslo, double* parlimitshi) {
0388   char FunName[100];
0389   sprintf(FunName, "Fitfcn_%s", hist->GetName());
0390   TF1* ffitold = (TF1*)gROOT->GetListOfFunctions()->FindObject(FunName);
0391   if (ffitold)
0392     delete ffitold;
0393 
0394   int npar(6);
0395   TObject* ob = gROOT->FindObject(FunName);
0396   if (ob)
0397     ob->Delete();
0398   TF1* ffit = new TF1(FunName, doubleGauss, fitrange[0], fitrange[1], npar);
0399   ffit->SetParameters(startvalues);
0400   ffit->SetLineColor(kBlue);
0401   ffit->SetParNames("Area1", "Mean1", "Width1", "Area2", "Mean2", "Width2");
0402   for (int i = 0; i < npar; i++)
0403     ffit->SetParLimits(i, parlimitslo[i], parlimitshi[i]);
0404   TFitResultPtr Fit = hist->Fit(FunName, "QRWLS");
0405   return Fit;
0406 }
0407 
0408 std::pair<double, double> fitLanGau(TH1D* hist, bool debug) {
0409   double rms;
0410   std::pair<double, double> mrms = GetMean(hist, 0.005, 0.25, rms);
0411   double mean = mrms.first;
0412   double LowEdge = 0.005, HighEdge = mean + 3 * rms;
0413   TFitResultPtr Fit1 = hist->Fit("gaus", "+0wwqRS", "", LowEdge, HighEdge);
0414   if (debug)
0415     std::cout << hist->GetName() << " 0 " << Fit1->Value(0) << " 1 " << Fit1->Value(1) << " 2 " << Fit1->Value(2)
0416               << std::endl;
0417   double startvalues[3];
0418   startvalues[0] = Fit1->Value(1);
0419   startvalues[1] = Fit1->Value(0);
0420   startvalues[2] = Fit1->Value(2);
0421 
0422   char FunName[100];
0423   sprintf(FunName, "Fitfcn_%s", hist->GetName());
0424   TF1* ffitold = (TF1*)gROOT->GetListOfFunctions()->FindObject(FunName);
0425   if (ffitold)
0426     delete ffitold;
0427 
0428   TObject* ob = gROOT->FindObject(FunName);
0429   if (ob)
0430     ob->Delete();
0431   TF1* ffit = new TF1(FunName, langaufun, LowEdge, HighEdge, 3);
0432   ffit->SetParameters(startvalues);
0433   ffit->SetParNames("MP", "Area", "GSigma");
0434   TFitResultPtr Fit2 = hist->Fit(FunName, "QRWLS");
0435   double value = Fit2->Value(0);
0436   double error = Fit2->FitResult::Error(0);
0437   if (debug)
0438     std::cout << hist->GetName() << " 0 " << Fit2->Value(0) << " 1 " << Fit2->Value(1) << " 2 " << Fit2->Value(2)
0439               << std::endl;
0440   return std::pair<double, double>(value, error);
0441 }
0442 
0443 results fitTwoGauss(TH1D* hist, bool debug) {
0444   double rms;
0445   std::pair<double, double> mrms = GetMean(hist, 0.2, 2.0, rms);
0446   double mean = mrms.first;
0447   double LowEdge = mean - fitrangeFactor * rms;
0448   double HighEdge = mean + fitrangeFactor * rms;
0449   if (LowEdge < 0.15)
0450     LowEdge = 0.15;
0451   std::string option = (hist->GetEntries() > 100) ? "QRS" : "QRWLS";
0452   TObject* ob = gROOT->FindObject("g1");
0453   if (ob)
0454     ob->Delete();
0455   TF1* g1 = new TF1("g1", "gaus", LowEdge, HighEdge);
0456   g1->SetLineColor(kGreen);
0457   TFitResultPtr Fit = hist->Fit(g1, option.c_str(), "");
0458 
0459   if (debug) {
0460     for (int k = 0; k < 3; ++k)
0461       std::cout << "Initial Parameter[" << k << "] = " << Fit->Value(k) << " +- " << Fit->FitResult::Error(k)
0462                 << std::endl;
0463   }
0464   double startvalues[6], fitrange[2], lowValue[6], highValue[6];
0465   startvalues[0] = Fit->Value(0);
0466   lowValue[0] = 0.5 * startvalues[0];
0467   highValue[0] = 2. * startvalues[0];
0468   startvalues[1] = Fit->Value(1);
0469   lowValue[1] = 0.5 * startvalues[1];
0470   highValue[1] = 2. * startvalues[1];
0471   startvalues[2] = Fit->Value(2);
0472   lowValue[2] = 0.5 * startvalues[2];
0473   highValue[2] = 2. * startvalues[2];
0474   startvalues[3] = 0.1 * Fit->Value(0);
0475   lowValue[3] = 0.0;
0476   highValue[3] = 10. * startvalues[3];
0477   startvalues[4] = Fit->Value(1);
0478   lowValue[4] = 0.5 * startvalues[4];
0479   highValue[4] = 2. * startvalues[4];
0480   startvalues[5] = 2.0 * Fit->Value(2);
0481   lowValue[5] = 0.5 * startvalues[5];
0482   highValue[5] = 100. * startvalues[5];
0483   //fitrange[0] = mean - 2.0*rms; fitrange[1] = mean + 2.0*rms;
0484   fitrange[0] = Fit->Value(1) - fitrangeFactor * Fit->Value(2);
0485   fitrange[1] = Fit->Value(1) + fitrangeFactor * Fit->Value(2);
0486   TFitResultPtr Fitfun = functionFit(hist, fitrange, startvalues, lowValue, highValue);
0487   double wt1 = (Fitfun->Value(0)) * (Fitfun->Value(2));
0488   double value1 = Fitfun->Value(1);
0489   double error1 = Fitfun->FitResult::Error(1);
0490   double wval1 = Fitfun->Value(2);
0491   double werr1 = Fitfun->FitResult::Error(2);
0492   double wt2 = (Fitfun->Value(3)) * (Fitfun->Value(5));
0493   double value2 = Fitfun->Value(4);
0494   double error2 = Fitfun->FitResult::Error(4);
0495   double wval2 = Fitfun->Value(5);
0496   double werr2 = Fitfun->FitResult::Error(5);
0497   double value = (wt1 * value1 + wt2 * value2) / (wt1 + wt2);
0498   double wval = (wt1 * wval1 + wt2 * wval2) / (wt1 + wt2);
0499   double error = (sqrt((wt1 * error1) * (wt1 * error1) + (wt2 * error2) * (wt2 * error2)) / (wt1 + wt2));
0500   double werror = (sqrt((wt1 * werr1) * (wt1 * werr1) + (wt2 * werr2) * (wt2 * werr2)) / (wt1 + wt2));
0501   std::cout << hist->GetName() << " Fit " << value << "+-" << error << " width " << wval << " +- " << werror
0502             << " First  " << value1 << "+-" << error1 << ":" << wval1 << "+-" << werr1 << ":" << wt1 << " Second "
0503             << value2 << "+-" << error2 << ":" << wval2 << "+-" << werr2 << ":" << wt2 << std::endl;
0504   if (debug) {
0505     for (int k = 0; k < 6; ++k)
0506       std::cout << hist->GetName() << ":Parameter[" << k << "] = " << Fitfun->Value(k) << " +- "
0507                 << Fitfun->FitResult::Error(k) << std::endl;
0508   }
0509   return results(value, error, wval, werror);
0510 }
0511 
0512 results fitOneGauss(TH1D* hist, bool fitTwice, bool debug) {
0513   double rms;
0514   std::pair<double, double> mrms = GetMean(hist, 0.2, 2.0, rms);
0515   double mean = mrms.first;
0516   double LowEdge = ((mean - fitrangeFactor * rms) < 0.5) ? 0.5 : (mean - fitrangeFactor * rms);
0517   double HighEdge = (mean + fitrangeFactor * rms);
0518   if (debug)
0519     std::cout << hist->GetName() << " Mean " << mean << " RMS " << rms << " Range " << LowEdge << ":" << HighEdge
0520               << "\n";
0521   std::string option = (hist->GetEntries() > 100) ? "QRS" : "QRWLS";
0522   TObject* ob = gROOT->FindObject("g1");
0523   if (ob)
0524     ob->Delete();
0525   TF1* g1 = new TF1("g1", "gaus", LowEdge, HighEdge);
0526   g1->SetLineColor(kGreen);
0527   TFitResultPtr Fit1 = hist->Fit(g1, option.c_str(), "");
0528   double value = Fit1->Value(1);
0529   double error = Fit1->FitResult::Error(1);
0530   double width = Fit1->Value(2);
0531   double werror = Fit1->FitResult::Error(2);
0532   if (fitTwice) {
0533     LowEdge = Fit1->Value(1) - fitrangeFactor * Fit1->Value(2);
0534     HighEdge = Fit1->Value(1) + fitrangeFactor * Fit1->Value(2);
0535     if (LowEdge < 0.5)
0536       LowEdge = 0.5;
0537     if (HighEdge > 5.0)
0538       HighEdge = 5.0;
0539     if (debug)
0540       std::cout << " Range for second Fit " << LowEdge << ":" << HighEdge << std::endl;
0541     TFitResultPtr Fit = hist->Fit("gaus", option.c_str(), "", LowEdge, HighEdge);
0542     value = Fit->Value(1);
0543     error = Fit->FitResult::Error(1);
0544     width = Fit->Value(2);
0545     werror = Fit->FitResult::Error(2);
0546   }
0547 
0548   std::pair<double, double> meaner = GetMean(hist, 0.2, 2.0, rms);
0549   std::pair<double, double> rmserr = GetWidth(hist, 0.2, 2.0);
0550   if (debug) {
0551     std::cout << "Fit " << value << ":" << error << ":" << hist->GetMeanError() << " Mean " << meaner.first << ":"
0552               << meaner.second << " Width " << rmserr.first << ":" << rmserr.second;
0553   }
0554   double minvalue(0.30);
0555   if (value < minvalue || value > 2.0 || error > 0.5) {
0556     value = meaner.first;
0557     error = meaner.second;
0558     width = rmserr.first;
0559     werror = rmserr.second;
0560   }
0561   if (debug) {
0562     std::cout << " Final " << value << "+-" << error << ":" << width << "+-" << werror << std::endl;
0563   }
0564   return results(value, error, width, werror);
0565 }
0566 
0567 void readCorrFactors(char* infile,
0568                      double scale,
0569                      std::map<int, cfactors>& cfacs,
0570                      int& etamin,
0571                      int& etamax,
0572                      int& maxdepth,
0573                      int iformat = 0,
0574                      bool debug = false) {
0575   cfacs.clear();
0576   std::ifstream fInput(infile);
0577   if (!fInput.good()) {
0578     std::cout << "Cannot open file " << infile << std::endl;
0579   } else {
0580     char buffer[1024];
0581     unsigned int all(0), good(0);
0582     while (fInput.getline(buffer, 1024)) {
0583       ++all;
0584       if (buffer[0] == '#')
0585         continue;  //ignore comment
0586       std::vector<std::string> items = splitString(std::string(buffer));
0587       if (items.size() != 5) {
0588         std::cout << "Ignore  line: " << buffer << std::endl;
0589       } else {
0590         ++good;
0591         int ieta = (iformat == 1) ? std::atoi(items[0].c_str()) : std::atoi(items[1].c_str());
0592         int depth = (iformat == 1) ? std::atoi(items[1].c_str()) : std::atoi(items[2].c_str());
0593         float corrf = std::atof(items[3].c_str());
0594         float dcorr = (iformat == 1) ? (0.02 * corrf) : std::atof(items[4].c_str());
0595         cfactors cfac(ieta, depth, scale * corrf, scale * dcorr);
0596         int detId = (iformat == 1) ? repackId(items[2], ieta, depth) : repackId(ieta, depth);
0597         cfacs[detId] = cfactors(ieta, depth, corrf, dcorr);
0598         if (ieta > etamax)
0599           etamax = ieta;
0600         if (ieta < etamin)
0601           etamin = ieta;
0602         if (depth > maxdepth)
0603           maxdepth = depth;
0604       }
0605     }
0606     fInput.close();
0607     std::cout << "Reads total of " << all << " and " << good << " good records"
0608               << " from " << infile << std::endl;
0609   }
0610   if (debug) {
0611     unsigned k(0);
0612     std::cout << "Eta Range " << etamin << ":" << etamax << " Max Depth " << maxdepth << std::endl;
0613     for (std::map<int, cfactors>::const_iterator itr = cfacs.begin(); itr != cfacs.end(); ++itr, ++k)
0614       std::cout << "[" << k << "] " << std::hex << itr->first << std::dec << ": " << (itr->second).ieta << " "
0615                 << (itr->second).depth << " " << (itr->second).corrf << " " << (itr->second).dcorr << std::endl;
0616   }
0617 }
0618 
0619 void FitHistStandard(std::string infile,
0620                      std::string outfile,
0621                      std::string prefix,
0622                      int mode = 11111,
0623                      int type = 0,
0624                      bool append = true,
0625                      bool saveAll = false,
0626                      bool debug = false) {
0627   int iname[6] = {0, 1, 2, 3, 4, 5};
0628   int checkmode[6] = {1000, 10, 10000, 1, 100000, 100};
0629   double xbin0[9] = {-21.0, -16.0, -12.0, -6.0, 0.0, 6.0, 12.0, 16.0, 21.0};
0630   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};
0631   double vbins[6] = {0.0, 7.0, 10.0, 13.0, 16.0, 50.0};
0632   double dlbins[9] = {0.0, 0.10, 0.20, 0.50, 1.0, 2.0, 2.5, 3.0, 10.0};
0633   std::string sname[4] = {"ratio", "etaR", "dl1R", "nvxR"};
0634   std::string lname[4] = {"Z", "E", "L", "V"};
0635   std::string wname[4] = {"W", "F", "M", "X"};
0636   std::string xname[4] = {"i#eta", "i#eta", "d_{L1}", "# Vertex"};
0637   int numb[4] = {10, 8, 8, 5};
0638 
0639   if (type == 0) {
0640     numb[0] = 8;
0641     for (int i = 0; i < 9; ++i)
0642       xbins[i] = xbin0[i];
0643   }
0644   TFile* file = new TFile(infile.c_str());
0645   std::vector<TH1D*> hists;
0646   char name[100], namw[100];
0647   if (file != nullptr) {
0648     for (int m1 = 0; m1 < 4; ++m1) {
0649       for (int m2 = 0; m2 < 6; ++m2) {
0650         sprintf(name, "%s%s%d0", prefix.c_str(), sname[m1].c_str(), iname[m2]);
0651         TH1D* hist0 = (TH1D*)file->FindObjectAny(name);
0652         bool ok = ((hist0 != nullptr) && (hist0->GetEntries() > 25));
0653         if ((mode / checkmode[m2]) % 10 > 0 && ok) {
0654           TH1D *histo(0), *histw(0);
0655           sprintf(name, "%s%s%d", prefix.c_str(), lname[m1].c_str(), iname[m2]);
0656           sprintf(namw, "%s%s%d", prefix.c_str(), wname[m1].c_str(), iname[m2]);
0657           if (m1 <= 1) {
0658             histo = new TH1D(name, hist0->GetTitle(), numb[m1], xbins);
0659             histw = new TH1D(namw, hist0->GetTitle(), numb[m1], xbins);
0660           } else if (m1 == 2) {
0661             histo = new TH1D(name, hist0->GetTitle(), numb[m1], dlbins);
0662             histw = new TH1D(namw, hist0->GetTitle(), numb[m1], dlbins);
0663           } else {
0664             histo = new TH1D(name, hist0->GetTitle(), numb[m1], vbins);
0665             histw = new TH1D(namw, hist0->GetTitle(), numb[m1], vbins);
0666           }
0667           int jmin(numb[m1]), jmax(0);
0668           for (int j = 0; j <= numb[m1]; ++j) {
0669             sprintf(name, "%s%s%d%d", prefix.c_str(), sname[m1].c_str(), iname[m2], j);
0670             TH1D* hist1 = (TH1D*)file->FindObjectAny(name);
0671             TH1D* hist = (TH1D*)hist1->Clone();
0672             double value(0), error(0), width(0), werror(0);
0673             if (hist->GetEntries() > 0) {
0674               value = hist->GetMean();
0675               error = hist->GetRMS();
0676               std::pair<double, double> rmserr = GetWidth(hist, 0.2, 2.0);
0677               width = rmserr.first;
0678               werror = rmserr.second;
0679             }
0680             if (hist->GetEntries() > 4) {
0681               bool flag = (j == 0) ? true : false;
0682               results meaner = fitOneGauss(hist, flag, debug);
0683               value = meaner.mean;
0684               error = meaner.errmean;
0685               width = meaner.width;
0686               werror = meaner.errwidth;
0687               if (j != 0) {
0688                 if (j < jmin)
0689                   jmin = j;
0690                 if (j > jmax)
0691                   jmax = j;
0692               }
0693             }
0694             if (j == 0) {
0695               hists.push_back(hist);
0696             } else {
0697               double wbyv = width / value;
0698               double wverr = wbyv * std::sqrt((werror * werror) / (width * width) + (error * error) / (value * value));
0699               if (saveAll)
0700                 hists.push_back(hist);
0701               histo->SetBinContent(j, value);
0702               histo->SetBinError(j, error);
0703               histw->SetBinContent(j, wbyv);
0704               histw->SetBinError(j, wverr);
0705             }
0706           }
0707           if (histo->GetEntries() > 2) {
0708             double LowEdge = histo->GetBinLowEdge(jmin);
0709             double HighEdge = histo->GetBinLowEdge(jmax) + histo->GetBinWidth(jmax);
0710             TFitResultPtr Fit = histo->Fit("pol0", "+QRWLS", "", LowEdge, HighEdge);
0711             if (debug) {
0712               std::cout << "Fit to Pol0: " << Fit->Value(0) << " +- " << Fit->FitResult::Error(0) << std::endl;
0713             }
0714             histo->GetXaxis()->SetTitle(xname[m1].c_str());
0715             histo->GetYaxis()->SetTitle("MPV(E_{HCAL}/(p-E_{ECAL}))");
0716             histo->GetYaxis()->SetRangeUser(0.4, 1.6);
0717             histw->GetXaxis()->SetTitle(xname[m1].c_str());
0718             histw->GetYaxis()->SetTitle("Width/MPV of (E_{HCAL}/(p-E_{ECAL}))");
0719             histw->GetYaxis()->SetRangeUser(0.0, 0.5);
0720           }
0721           hists.push_back(histo);
0722           hists.push_back(histw);
0723         }
0724       }
0725     }
0726     TFile* theFile(0);
0727     if (append) {
0728       theFile = new TFile(outfile.c_str(), "UPDATE");
0729     } else {
0730       theFile = new TFile(outfile.c_str(), "RECREATE");
0731     }
0732     theFile->cd();
0733     for (unsigned int i = 0; i < hists.size(); ++i) {
0734       TH1D* hnew = (TH1D*)hists[i]->Clone();
0735       hnew->Write();
0736     }
0737     theFile->Close();
0738     file->Close();
0739   }
0740 }
0741 
0742 void FitHistExtended(const char* infile,
0743                      const char* outfile,
0744                      std::string prefix,
0745                      int numb = 50,
0746                      int type = 3,
0747                      bool append = true,
0748                      bool fiteta = true,
0749                      int iname = 3,
0750                      bool debug = false) {
0751   std::string sname("ratio"), lname("Z"), wname("W"), ename("etaB");
0752   double xbins[99];
0753   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,
0754                      3.0,   5.0,   7.0,   9.0,   11.0,  13.0,  15.0,  17.0, 19.0, 21.0, 23.0};
0755   if (type == 2) {
0756     numb = 22;
0757     for (int k = 0; k <= numb; ++k)
0758       xbins[k] = xbin[k];
0759   } else if (type == 1) {
0760     numb = 1;
0761     xbins[0] = -25;
0762     xbins[1] = 25;
0763   } else {
0764     int neta = numb / 2;
0765     for (int k = 0; k < neta; ++k) {
0766       xbins[k] = (k - neta) - 0.5;
0767       xbins[numb - k] = (neta - k) + 0.5;
0768     }
0769     xbins[neta] = 0;
0770   }
0771   if (debug) {
0772     for (int k = 0; k <= numb; ++k)
0773       std::cout << " " << xbins[k];
0774     std::cout << std::endl;
0775   }
0776   TFile* file = new TFile(infile);
0777   std::vector<TH1D*> hists;
0778   char name[200];
0779   if (debug) {
0780     std::cout << infile << " " << file << std::endl;
0781   }
0782   if (file != nullptr) {
0783     sprintf(name, "%s%s%d0", prefix.c_str(), sname.c_str(), iname);
0784     TH1D* hist0 = (TH1D*)file->FindObjectAny(name);
0785     bool ok = (hist0 != nullptr);
0786     if (debug) {
0787       std::cout << name << " Pointer " << hist0 << " " << ok << std::endl;
0788     }
0789     if (ok) {
0790       TH1D *histo(0), *histw(0);
0791       if (numb > 0) {
0792         sprintf(name, "%s%s%d", prefix.c_str(), lname.c_str(), iname);
0793         histo = new TH1D(name, hist0->GetTitle(), numb, xbins);
0794         sprintf(name, "%s%s%d", prefix.c_str(), wname.c_str(), iname);
0795         histw = new TH1D(name, hist0->GetTitle(), numb, xbins);
0796         if (debug)
0797           std::cout << name << " " << histo->GetNbinsX() << std::endl;
0798       }
0799       if (hist0->GetEntries() > 10) {
0800         double rms;
0801         results meaner0 = fitOneGauss(hist0, true, debug);
0802         std::pair<double, double> meaner1 = GetMean(hist0, 0.2, 2.0, rms);
0803         std::pair<double, double> meaner2 = GetWidth(hist0, 0.2, 2.0);
0804         if (debug) {
0805           std::cout << "Fit " << meaner0.mean << ":" << meaner0.errmean << " Mean1 " << hist0->GetMean() << ":"
0806                     << hist0->GetMeanError() << " Mean2 " << meaner1.first << ":" << meaner1.second << " Width "
0807                     << meaner2.first << ":" << meaner2.second << std::endl;
0808         }
0809       }
0810       int nv1(100), nv2(0);
0811       int jmin(numb), jmax(0);
0812       for (int j = 0; j <= numb; ++j) {
0813         sprintf(name, "%s%s%d%d", prefix.c_str(), sname.c_str(), iname, j);
0814         TH1D* hist1 = (TH1D*)file->FindObjectAny(name);
0815         if (debug) {
0816           std::cout << "Get Histogram for " << name << " at " << hist1 << std::endl;
0817         }
0818         double value(0), error(0), total(0), width(0), werror(0);
0819         if (hist1 == nullptr) {
0820           value = 1.0;
0821         } else {
0822           TH1D* hist = (TH1D*)hist1->Clone();
0823           if (debug)
0824             std::cout << "Histogram " << name << ":" << (hist->GetName()) << " with " << (hist->GetEntries())
0825                       << " entries" << std::endl;
0826           if (hist->GetEntries() > 0) {
0827             value = hist->GetMean();
0828             error = hist->GetRMS();
0829             for (int i = 1; i <= hist->GetNbinsX(); ++i)
0830               total += hist->GetBinContent(i);
0831             std::pair<double, double> rmserr = GetWidth(hist, 0.2, 2.0);
0832             width = rmserr.first;
0833             werror = rmserr.second;
0834           }
0835           if (total > 4) {
0836             if (nv1 > j)
0837               nv1 = j;
0838             if (nv2 < j)
0839               nv2 = j;
0840             if (j == 0) {
0841               sprintf(name, "%sOne", hist1->GetName());
0842               TH1D* hist2 = (TH1D*)hist1->Clone(name);
0843               fitOneGauss(hist2, true, debug);
0844               hists.push_back(hist2);
0845               results meaner = fitOneGauss(hist, true, debug);
0846               value = meaner.mean;
0847               error = meaner.errmean;
0848               width = meaner.width;
0849               werror = meaner.errwidth;
0850             } else {
0851               results meaner = fitOneGauss(hist, true, debug);
0852               value = meaner.mean;
0853               error = meaner.errmean;
0854               width = meaner.width;
0855               werror = meaner.errwidth;
0856             }
0857             if (j != 0) {
0858               if (j < jmin)
0859                 jmin = j;
0860               if (j > jmax)
0861                 jmax = j;
0862             }
0863           }
0864           hists.push_back(hist);
0865         }
0866         if (debug) {
0867           std::cout << "Hist " << j << " Value " << value << " +- " << error << std::endl;
0868         }
0869         if (j != 0) {
0870           double wbyv = width / value;
0871           double wverr = wbyv * std::sqrt((werror * werror) / (width * width) + (error * error) / (value * value));
0872           histo->SetBinContent(j, value);
0873           histo->SetBinError(j, error);
0874           histw->SetBinContent(j, wbyv);
0875           histw->SetBinError(j, wverr);
0876         }
0877       }
0878       if (histo != nullptr) {
0879         if (histo->GetEntries() > 2 && fiteta) {
0880           if (debug) {
0881             std::cout << "Jmin/max " << jmin << ":" << jmax << ":" << histo->GetNbinsX() << std::endl;
0882           }
0883           double LowEdge = histo->GetBinLowEdge(jmin);
0884           double HighEdge = histo->GetBinLowEdge(jmax) + histo->GetBinWidth(jmax);
0885           TFitResultPtr Fit = histo->Fit("pol0", "+QRWLS", "", LowEdge, HighEdge);
0886           if (debug) {
0887             std::cout << "Fit to Pol0: " << Fit->Value(0) << " +- " << Fit->FitResult::Error(0) << " in range " << nv1
0888                       << ":" << xbins[nv1] << ":" << nv2 << ":" << xbins[nv2] << std::endl;
0889           }
0890           histo->GetXaxis()->SetTitle("i#eta");
0891           histo->GetYaxis()->SetTitle("MPV(E_{HCAL}/(p-E_{ECAL}))");
0892           histo->GetYaxis()->SetRangeUser(0.4, 1.6);
0893           histw->GetXaxis()->SetTitle("i#eta");
0894           histw->GetYaxis()->SetTitle("MPV/Width(E_{HCAL}/(p-E_{ECAL}))");
0895           histw->GetYaxis()->SetRangeUser(0.0, 0.5);
0896         }
0897         hists.push_back(histo);
0898         hists.push_back(histw);
0899       } else {
0900         hists.push_back(hist0);
0901       }
0902 
0903       //Barrel,Endcap
0904       for (int j = 1; j <= 4; ++j) {
0905         sprintf(name, "%s%s%d%d", prefix.c_str(), ename.c_str(), iname, j);
0906         TH1D* hist1 = (TH1D*)file->FindObjectAny(name);
0907         if (debug) {
0908           std::cout << "Get Histogram for " << name << " at " << hist1 << std::endl;
0909         }
0910         if (hist1 != nullptr) {
0911           TH1D* hist = (TH1D*)hist1->Clone();
0912           double value(0), error(0), total(0), width(0), werror(0);
0913           if (hist->GetEntries() > 0) {
0914             value = hist->GetMean();
0915             error = hist->GetRMS();
0916             for (int i = 1; i <= hist->GetNbinsX(); ++i)
0917               total += hist->GetBinContent(i);
0918           }
0919           if (total > 4) {
0920             sprintf(name, "%sOne", hist1->GetName());
0921             TH1D* hist2 = (TH1D*)hist1->Clone(name);
0922             results meanerr = fitOneGauss(hist2, true, debug);
0923             value = meanerr.mean;
0924             error = meanerr.errmean;
0925             width = meanerr.width;
0926             werror = meanerr.errwidth;
0927             double wbyv = width / value;
0928             double wverr = wbyv * std::sqrt((werror * werror) / (width * width) + (error * error) / (value * value));
0929             std::cout << hist2->GetName() << " MPV " << value << " +- " << error << " Width " << width << " +- "
0930                       << werror << " W/M " << wbyv << " +- " << wverr << std::endl;
0931             hists.push_back(hist2);
0932             if (hist1->GetBinLowEdge(1) < 0.1) {
0933               sprintf(name, "%sTwo", hist1->GetName());
0934               TH1D* hist3 = (TH1D*)hist1->Clone(name);
0935               fitLanGau(hist3, debug);
0936               hists.push_back(hist3);
0937             }
0938             //            results meaner0 = fitTwoGauss(hist, debug);
0939             results meaner0 = fitOneGauss(hist, true, debug);
0940             value = meaner0.mean;
0941             error = meaner0.errmean;
0942             double rms;
0943             std::pair<double, double> meaner = GetMean(hist, 0.2, 2.0, rms);
0944             if (debug) {
0945               std::cout << "Fit " << value << ":" << error << ":" << hist->GetMeanError() << " Mean " << meaner.first
0946                         << ":" << meaner.second << std::endl;
0947             }
0948           }
0949           hists.push_back(hist);
0950         }
0951       }
0952     }
0953     TFile* theFile(0);
0954     if (append) {
0955       if (debug) {
0956         std::cout << "Open file " << outfile << " in append mode" << std::endl;
0957       }
0958       theFile = new TFile(outfile, "UPDATE");
0959     } else {
0960       if (debug) {
0961         std::cout << "Open file " << outfile << " in recreate mode" << std::endl;
0962       }
0963       theFile = new TFile(outfile, "RECREATE");
0964     }
0965 
0966     theFile->cd();
0967     for (unsigned int i = 0; i < hists.size(); ++i) {
0968       TH1D* hnew = (TH1D*)hists[i]->Clone();
0969       if (debug) {
0970         std::cout << "Write Histogram " << hnew->GetTitle() << std::endl;
0971       }
0972       hnew->Write();
0973     }
0974     theFile->Close();
0975     file->Close();
0976   }
0977 }
0978 
0979 void FitHistExtended2(const char* infile,
0980                       const char* outfile,
0981                       std::string prefix,
0982                       int numb = 50,
0983                       bool append = true,
0984                       int iname = 3,
0985                       bool debug = false) {
0986   std::string sname("ratio"), lname("Z"), wname("W"), ename("etaB");
0987   double xbins[99];
0988   int ixbin[99];
0989   int neta = numb / 2;
0990   for (int k = 0; k < neta; ++k) {
0991     xbins[k] = (k - neta) - 0.5;
0992     xbins[numb - k] = (neta - k) + 0.5;
0993     ixbin[k] = (k - neta);
0994     ixbin[numb - k] = (neta - k);
0995   }
0996   xbins[neta] = 0;
0997   ixbin[neta] = 0;
0998   TFile* file = new TFile(infile);
0999   std::vector<TH1D*> hists;
1000   char name[200];
1001   if (debug) {
1002     std::cout << infile << " " << file << std::endl;
1003   }
1004   if (file != nullptr) {
1005     sprintf(name, "%s%s%d0", prefix.c_str(), sname.c_str(), iname);
1006     TH1D* hist0 = (TH1D*)file->FindObjectAny(name);
1007     bool ok = (hist0 != nullptr);
1008     if (debug) {
1009       std::cout << name << " Pointer " << hist0 << " " << ok << std::endl;
1010     }
1011     if (ok) {
1012       TH1D *histo(0), *histw(0);
1013       if (numb > 0) {
1014         sprintf(name, "%s%s%d", prefix.c_str(), lname.c_str(), iname);
1015         histo = new TH1D(name, hist0->GetTitle(), numb, xbins);
1016         sprintf(name, "%s%s%d", prefix.c_str(), wname.c_str(), iname);
1017         histw = new TH1D(name, hist0->GetTitle(), numb, xbins);
1018         if (debug)
1019           std::cout << name << " " << histo->GetNbinsX() << std::endl;
1020       }
1021       int nv1(100), nv2(0);
1022       int jmin(numb), jmax(0);
1023       for (int j = 0; j <= numb; ++j) {
1024         sprintf(name, "%s%s%d%d", prefix.c_str(), sname.c_str(), iname, j);
1025         TH1D* hist1 = (TH1D*)file->FindObjectAny(name);
1026         if (debug) {
1027           std::cout << "Get Histogram for " << name << " at " << hist1 << std::endl;
1028         }
1029         double value(0), error(0), total(0), width(0), werror(0);
1030         if (hist1 == nullptr) {
1031           value = 1.0;
1032         } else {
1033           TH1D* hist = (TH1D*)hist1->Clone();
1034           if (debug)
1035             std::cout << "Histogram " << name << ":" << (hist->GetName()) << " with " << (hist->GetEntries())
1036                       << " entries" << std::endl;
1037           if (hist->GetEntries() > 0) {
1038             value = hist->GetMean();
1039             error = hist->GetRMS();
1040             for (int i = 1; i <= hist->GetNbinsX(); ++i)
1041               total += hist->GetBinContent(i);
1042             std::pair<double, double> rmserr = GetWidth(hist, 0.2, 2.0);
1043             width = rmserr.first;
1044             werror = rmserr.second;
1045           }
1046           if (total > 4) {
1047             if (nv1 > j)
1048               nv1 = j;
1049             if (nv2 < j)
1050               nv2 = j;
1051             if (j == 0) {
1052               sprintf(name, "%sOne", hist1->GetName());
1053               TH1D* hist2 = (TH1D*)hist1->Clone(name);
1054               fitOneGauss(hist2, true, debug);
1055               hists.push_back(hist2);
1056               results meaner = fitOneGauss(hist, true, debug);
1057               value = meaner.mean;
1058               error = meaner.errmean;
1059               width = meaner.width;
1060               werror = meaner.errwidth;
1061             } else {
1062               results meaner = fitOneGauss(hist, true, debug);
1063               value = meaner.mean;
1064               error = meaner.errmean;
1065               width = meaner.width;
1066               werror = meaner.errwidth;
1067             }
1068             if (j != 0) {
1069               if (j < jmin)
1070                 jmin = j;
1071               if (j > jmax)
1072                 jmax = j;
1073             }
1074             double wbyv0 = width / value;
1075             double wverr0 = wbyv0 * std::sqrt((werror * werror) / (width * width) + (error * error) / (value * value));
1076             if (ixbin[j] != 0)
1077               std::cout << ixbin[j] << " MPV " << value << " +- " << error << " Width/MPV " << wbyv0 << " +- " << wverr0
1078                         << std::endl;
1079           }
1080           hists.push_back(hist);
1081         }
1082         if (debug) {
1083           std::cout << "Hist " << j << " Value " << value << " +- " << error << std::endl;
1084         }
1085         if (j != 0) {
1086           double wbyv = width / value;
1087           double wverr = wbyv * std::sqrt((werror * werror) / (width * width) + (error * error) / (value * value));
1088           histo->SetBinContent(j, value);
1089           histo->SetBinError(j, error);
1090           histw->SetBinContent(j, wbyv);
1091           histw->SetBinError(j, wverr);
1092         }
1093       }
1094       if (histo != nullptr) {
1095         if (histo->GetEntries() > 2) {
1096           if (debug) {
1097             std::cout << "Jmin/max " << jmin << ":" << jmax << ":" << histo->GetNbinsX() << std::endl;
1098           }
1099           double LowEdge = histo->GetBinLowEdge(jmin);
1100           double HighEdge = histo->GetBinLowEdge(jmax) + histo->GetBinWidth(jmax);
1101           TFitResultPtr Fit = histo->Fit("pol0", "+QRWLS", "", LowEdge, HighEdge);
1102           if (debug) {
1103             std::cout << "Fit to Pol0: " << Fit->Value(0) << " +- " << Fit->FitResult::Error(0) << " in range " << nv1
1104                       << ":" << xbins[nv1] << ":" << nv2 << ":" << xbins[nv2] << std::endl;
1105           }
1106           histo->GetXaxis()->SetTitle("i#eta");
1107           histo->GetYaxis()->SetTitle("MPV(E_{HCAL}/(p-E_{ECAL}))");
1108           histo->GetYaxis()->SetRangeUser(0.4, 1.6);
1109           histw->GetXaxis()->SetTitle("i#eta");
1110           histw->GetYaxis()->SetTitle("MPV/Width(E_{HCAL}/(p-E_{ECAL}))");
1111           histw->GetYaxis()->SetRangeUser(0.0, 0.5);
1112         }
1113         hists.push_back(histo);
1114         hists.push_back(histw);
1115       } else {
1116         hists.push_back(hist0);
1117       }
1118 
1119       //Barrel,Endcap
1120       for (int j = 1; j <= 4; ++j) {
1121         sprintf(name, "%s%s%d%d", prefix.c_str(), ename.c_str(), iname, j);
1122         TH1D* hist1 = (TH1D*)file->FindObjectAny(name);
1123         if (debug) {
1124           std::cout << "Get Histogram for " << name << " at " << hist1 << std::endl;
1125         }
1126         if (hist1 != nullptr) {
1127           TH1D* hist = (TH1D*)hist1->Clone();
1128           double value(0), error(0), total(0), width(0), werror(0);
1129           if (hist->GetEntries() > 0) {
1130             value = hist->GetMean();
1131             error = hist->GetRMS();
1132             for (int i = 1; i <= hist->GetNbinsX(); ++i)
1133               total += hist->GetBinContent(i);
1134           }
1135           if (total > 4) {
1136             sprintf(name, "%sOne", hist1->GetName());
1137             TH1D* hist2 = (TH1D*)hist1->Clone(name);
1138             results meanerr = fitOneGauss(hist2, true, debug);
1139             value = meanerr.mean;
1140             error = meanerr.errmean;
1141             width = meanerr.width;
1142             werror = meanerr.errwidth;
1143             double wbyv = width / value;
1144             double wverr = wbyv * std::sqrt((werror * werror) / (width * width) + (error * error) / (value * value));
1145             std::cout << hist2->GetName() << " MPV " << value << " +- " << error << " Width " << width << " +- "
1146                       << werror << " W/M " << wbyv << " +- " << wverr << std::endl;
1147             hists.push_back(hist2);
1148             results meaner0 = fitOneGauss(hist, true, debug);
1149             value = meaner0.mean;
1150             error = meaner0.errmean;
1151             double rms;
1152             std::pair<double, double> meaner = GetMean(hist, 0.2, 2.0, rms);
1153             if (debug) {
1154               std::cout << "Fit " << value << ":" << error << ":" << hist->GetMeanError() << " Mean " << meaner.first
1155                         << ":" << meaner.second << std::endl;
1156             }
1157           }
1158           hists.push_back(hist);
1159         }
1160       }
1161     }
1162     TFile* theFile(0);
1163     if (append) {
1164       if (debug) {
1165         std::cout << "Open file " << outfile << " in append mode" << std::endl;
1166       }
1167       theFile = new TFile(outfile, "UPDATE");
1168     } else {
1169       if (debug) {
1170         std::cout << "Open file " << outfile << " in recreate mode" << std::endl;
1171       }
1172       theFile = new TFile(outfile, "RECREATE");
1173     }
1174 
1175     theFile->cd();
1176     for (unsigned int i = 0; i < hists.size(); ++i) {
1177       TH1D* hnew = (TH1D*)hists[i]->Clone();
1178       if (debug) {
1179         std::cout << "Write Histogram " << hnew->GetTitle() << std::endl;
1180       }
1181       hnew->Write();
1182     }
1183     theFile->Close();
1184     file->Close();
1185   }
1186 }
1187 
1188 void FitHistRBX(const char* infile, const char* outfile, std::string prefix, bool append = true, int iname = 3) {
1189   std::string sname("RBX"), lname("R");
1190   int numb(18);
1191   bool debug(false);
1192   char name[200];
1193 
1194   TFile* file = new TFile(infile);
1195   std::vector<TH1D*> hists;
1196   sprintf(name, "%s%s%d", prefix.c_str(), lname.c_str(), iname);
1197   TH1D* histo = new TH1D(name, "", numb, 0, numb);
1198   if (debug) {
1199     std::cout << infile << " " << file << std::endl;
1200   }
1201   if (file != nullptr) {
1202     for (int j = 0; j < numb; ++j) {
1203       sprintf(name, "%s%s%d%d", prefix.c_str(), sname.c_str(), iname, j + 1);
1204       TH1D* hist1 = (TH1D*)file->FindObjectAny(name);
1205       if (debug) {
1206         std::cout << "Get Histogram for " << name << " at " << hist1 << std::endl;
1207       }
1208       TH1D* hist = (TH1D*)hist1->Clone();
1209       double value(0), error(0), total(0);
1210       if (hist->GetEntries() > 0) {
1211         value = hist->GetMean();
1212         error = hist->GetRMS();
1213         for (int i = 1; i <= hist->GetNbinsX(); ++i)
1214           total += hist->GetBinContent(i);
1215       }
1216       if (total > 4) {
1217         results meaner = fitOneGauss(hist, false, debug);
1218         value = meaner.mean;
1219         error = meaner.errmean;
1220       }
1221       hists.push_back(hist);
1222       histo->SetBinContent(j + 1, value);
1223       histo->SetBinError(j + 1, error);
1224     }
1225     histo->GetXaxis()->SetTitle("RBX #");
1226     histo->GetYaxis()->SetTitle("MPV(E_{HCAL}/(p-E_{ECAL}))");
1227     histo->GetYaxis()->SetRangeUser(0.75, 1.20);
1228     hists.push_back(histo);
1229 
1230     TFile* theFile(0);
1231     if (append) {
1232       if (debug) {
1233         std::cout << "Open file " << outfile << " in append mode" << std::endl;
1234       }
1235       theFile = new TFile(outfile, "UPDATE");
1236     } else {
1237       if (debug) {
1238         std::cout << "Open file " << outfile << " in recreate mode" << std::endl;
1239       }
1240       theFile = new TFile(outfile, "RECREATE");
1241     }
1242 
1243     theFile->cd();
1244     for (unsigned int i = 0; i < hists.size(); ++i) {
1245       TH1D* hnew = (TH1D*)hists[i]->Clone();
1246       if (debug) {
1247         std::cout << "Write Histogram " << hnew->GetTitle() << std::endl;
1248       }
1249       hnew->Write();
1250     }
1251     theFile->Close();
1252     file->Close();
1253   }
1254 }
1255 
1256 void PlotHist(const char* infile,
1257               std::string prefix,
1258               std::string text,
1259               int mode = 4,
1260               int kopt = 100,
1261               double lumi = 0,
1262               double ener = 13.0,
1263               bool isRealData = false,
1264               bool drawStatBox = true,
1265               int save = 0) {
1266   std::string name0[6] = {"ratio00", "ratio10", "ratio20", "ratio30", "ratio40", "ratio50"};
1267   std::string name1[5] = {"Z0", "Z1", "Z2", "Z3", "Z4"};
1268   std::string name2[5] = {"L0", "L1", "L2", "L3", "L4"};
1269   std::string name3[5] = {"V0", "V1", "V2", "V3", "V4"};
1270   std::string name4[20] = {"etaB41", "etaB42", "etaB43", "etaB44", "etaB31", "etaB32", "etaB33",
1271                            "etaB34", "etaB21", "etaB22", "etaB23", "etaB24", "etaB11", "etaB12",
1272                            "etaB13", "etaB14", "etaB01", "etaB02", "etaB03", "etaB04"};
1273   std::string name5[5] = {"W0", "W1", "W2", "W3", "W4"};
1274   std::string title[6] = {"Tracks with p = 10:20 GeV",
1275                           "Tracks with p = 20:30 GeV",
1276                           "Tracks with p = 30:40 GeV",
1277                           "Tracks with p = 40:60 GeV",
1278                           "Tracks with p = 60:100 GeV",
1279                           "Tracks with p = 20:100 GeV"};
1280   std::string title1[20] = {"Tracks with p = 60:100 GeV (Barrel)", "Tracks with p = 60:100 GeV (Transition)",
1281                             "Tracks with p = 60:100 GeV (Endcap)", "Tracks with p = 60:100 GeV",
1282                             "Tracks with p = 40:60 GeV (Barrel)",  "Tracks with p = 40:60 GeV (Transition)",
1283                             "Tracks with p = 40:60 GeV (Endcap)",  "Tracks with p = 40:60 GeV",
1284                             "Tracks with p = 30:40 GeV (Barrel)",  "Tracks with p = 30:40 GeV (Transition)",
1285                             "Tracks with p = 30:40 GeV (Endcap)",  "Tracks with p = 30:40 GeV",
1286                             "Tracks with p = 20:30 GeV (Barrel)",  "Tracks with p = 20:30 GeV (Transition)",
1287                             "Tracks with p = 20:30 GeV (Endcap)",  "Tracks with p = 20:30 GeV",
1288                             "Tracks with p = 10:20 GeV (Barrel)",  "Tracks with p = 10:20 GeV (Transition)",
1289                             "Tracks with p = 10:20 GeV (Endcap)",  "Tracks with p = 10:20 GeV"};
1290   std::string xtitl[5] = {"E_{HCAL}/(p-E_{ECAL})", "i#eta", "d_{L1}", "# Vertex", "E_{HCAL}/(p-E_{ECAL})"};
1291   std::string ytitl[5] = {
1292       "Tracks", "MPV(E_{HCAL}/(p-E_{ECAL}))", "MPV(E_{HCAL}/(p-E_{ECAL}))", "MPV(E_{HCAL}/(p-E_{ECAL}))", "Tracks"};
1293 
1294   gStyle->SetCanvasBorderMode(0);
1295   gStyle->SetCanvasColor(kWhite);
1296   gStyle->SetPadColor(kWhite);
1297   gStyle->SetFillColor(kWhite);
1298   gStyle->SetOptTitle(0);
1299   if (mode < 0 || mode > 5)
1300     mode = 0;
1301   if (drawStatBox) {
1302     int iopt(1110);
1303     if (mode != 0)
1304       iopt = 10;
1305     gStyle->SetOptStat(iopt);
1306     gStyle->SetOptFit(1);
1307   } else {
1308     gStyle->SetOptStat(0);
1309     gStyle->SetOptFit(0);
1310   }
1311   TFile* file = new TFile(infile);
1312   TLine* line(0);
1313   char name[100], namep[100];
1314   int kmax = (mode == 4) ? 20 : (((mode < 1) && (mode > 5)) ? 6 : 5);
1315   for (int k = 0; k < kmax; ++k) {
1316     if (mode == 1) {
1317       sprintf(name, "%s%s", prefix.c_str(), name1[k].c_str());
1318     } else if (mode == 2) {
1319       sprintf(name, "%s%s", prefix.c_str(), name2[k].c_str());
1320     } else if (mode == 3) {
1321       sprintf(name, "%s%s", prefix.c_str(), name3[k].c_str());
1322     } else if (mode == 4) {
1323       if ((kopt / 100) % 10 == 0) {
1324         sprintf(name, "%s%s", prefix.c_str(), name4[k].c_str());
1325       } else if ((kopt / 100) % 10 == 2) {
1326         sprintf(name, "%s%sTwo", prefix.c_str(), name4[k].c_str());
1327       } else {
1328         sprintf(name, "%s%sOne", prefix.c_str(), name4[k].c_str());
1329       }
1330     } else if (mode == 5) {
1331       sprintf(name, "%s%s", prefix.c_str(), name5[k].c_str());
1332     } else {
1333       if ((kopt / 100) % 10 == 0) {
1334         sprintf(name, "%s%s", prefix.c_str(), name0[k].c_str());
1335       } else if ((kopt / 100) % 10 == 2) {
1336         sprintf(name, "%s%sTwo", prefix.c_str(), name0[k].c_str());
1337       } else {
1338         sprintf(name, "%s%sOne", prefix.c_str(), name0[k].c_str());
1339       }
1340     }
1341     TH1D* hist1 = (TH1D*)file->FindObjectAny(name);
1342     if (hist1 != nullptr) {
1343       TH1D* hist = (TH1D*)(hist1->Clone());
1344       double p0(1);
1345       double ymin(0.90);
1346       sprintf(namep, "c_%s", name);
1347       TCanvas* pad = new TCanvas(namep, namep, 700, 500);
1348       pad->SetRightMargin(0.10);
1349       pad->SetTopMargin(0.10);
1350       if ((kopt / 10) % 10 > 0)
1351         gPad->SetGrid();
1352       hist->GetXaxis()->SetTitleSize(0.04);
1353       hist->GetXaxis()->SetTitle(xtitl[mode].c_str());
1354       hist->GetYaxis()->SetTitle(ytitl[mode].c_str());
1355       hist->GetYaxis()->SetLabelOffset(0.005);
1356       hist->GetYaxis()->SetTitleSize(0.04);
1357       hist->GetYaxis()->SetLabelSize(0.035);
1358       hist->GetYaxis()->SetTitleOffset(1.10);
1359       if (mode == 0 || mode == 4) {
1360         if ((kopt / 100) % 10 == 2)
1361           hist->GetXaxis()->SetRangeUser(0.0, 0.30);
1362         else
1363           hist->GetXaxis()->SetRangeUser(0.25, 2.25);
1364       } else {
1365         if (mode == 5)
1366           hist->GetYaxis()->SetRangeUser(0.1, 0.50);
1367         else if (isRealData)
1368           hist->GetYaxis()->SetRangeUser(0.5, 1.50);
1369         else
1370           hist->GetYaxis()->SetRangeUser(0.8, 1.20);
1371         if (kopt % 10 > 0) {
1372           int nbin = hist->GetNbinsX();
1373           double LowEdge = (kopt % 10 == 1) ? hist->GetBinLowEdge(1) : -20;
1374           double HighEdge = (kopt % 10 == 1) ? hist->GetBinLowEdge(nbin) + hist->GetBinWidth(nbin) : 20;
1375           TFitResultPtr Fit = hist->Fit("pol0", "+QRWLS", "", LowEdge, HighEdge);
1376           p0 = Fit->Value(0);
1377         }
1378       }
1379       hist->SetMarkerStyle(20);
1380       hist->SetMarkerColor(2);
1381       hist->SetLineColor(2);
1382       hist->Draw();
1383       pad->Update();
1384       TPaveStats* st1 = (TPaveStats*)hist->GetListOfFunctions()->FindObject("stats");
1385       if (st1 != nullptr) {
1386         ymin = (mode == 0 || mode == 4) ? 0.70 : 0.80;
1387         st1->SetY1NDC(ymin);
1388         st1->SetY2NDC(0.90);
1389         st1->SetX1NDC(0.65);
1390         st1->SetX2NDC(0.90);
1391       }
1392       if (mode != 0 && mode != 4 && kopt % 10 > 0) {
1393         double xmin = hist->GetBinLowEdge(1);
1394         int nbin = hist->GetNbinsX();
1395         double xmax = hist->GetBinLowEdge(nbin) + hist->GetBinWidth(nbin);
1396         double mean(0), rms(0), total(0);
1397         int kount(0);
1398         for (int k = 2; k < nbin; ++k) {
1399           double x = hist->GetBinContent(k);
1400           double w = hist->GetBinError(k);
1401           mean += (x * w);
1402           rms += (x * x * w);
1403           total += w;
1404           ++kount;
1405         }
1406         mean /= total;
1407         rms /= total;
1408         double error = sqrt(rms - mean * mean) / total;
1409         line = new TLine(xmin, p0, xmax, p0);  //etamin,1.0,etamax,1.0);
1410         std::cout << xmin << ":" << xmax << ":" << p0 << " Mean " << nbin << ":" << kount << ":" << total << ":" << mean
1411                   << ":" << rms << ":" << error << std::endl;
1412         line->SetLineWidth(2);
1413         line->SetLineStyle(2);
1414         line->Draw("same");
1415       }
1416       pad->Update();
1417       double ymx(0.96), xmi(0.25), xmx(0.90);
1418       char txt[100];
1419       if (lumi > 0.1) {
1420         ymx = ymin - 0.005;
1421         xmi = 0.45;
1422         TPaveText* txt0 = new TPaveText(0.65, 0.91, 0.90, 0.96, "blNDC");
1423         txt0->SetFillColor(0);
1424         sprintf(txt, "%4.1f TeV %5.1f fb^{-1}", ener, lumi);
1425         txt0->AddText(txt);
1426         txt0->Draw("same");
1427       }
1428       double ymi = ymx - 0.05;
1429       TPaveText* txt1 = new TPaveText(xmi, ymi, xmx, ymx, "blNDC");
1430       txt1->SetFillColor(0);
1431       if (text == "") {
1432         if (mode == 4)
1433           sprintf(txt, "%s", title1[k].c_str());
1434         else
1435           sprintf(txt, "%s", title[k].c_str());
1436       } else {
1437         if (mode == 4)
1438           sprintf(txt, "%s (%s)", title1[k].c_str(), text.c_str());
1439         else
1440           sprintf(txt, "%s (%s)", title[k].c_str(), text.c_str());
1441       }
1442       txt1->AddText(txt);
1443       txt1->Draw("same");
1444       double xmax = (isRealData) ? 0.33 : 0.44;
1445       ymi = (lumi > 0.1) ? 0.91 : 0.84;
1446       ymx = ymi + 0.05;
1447       TPaveText* txt2 = new TPaveText(0.11, ymi, xmax, ymx, "blNDC");
1448       txt2->SetFillColor(0);
1449       if (isRealData)
1450         sprintf(txt, "CMS Preliminary");
1451       else
1452         sprintf(txt, "CMS Simulation Preliminary");
1453       txt2->AddText(txt);
1454       txt2->Draw("same");
1455       pad->Modified();
1456       pad->Update();
1457       if (save > 0) {
1458         sprintf(name, "%s.pdf", pad->GetName());
1459         pad->Print(name);
1460       } else if (save < 0) {
1461         sprintf(name, "%s.C", pad->GetName());
1462         pad->Print(name);
1463       }
1464     }
1465   }
1466 }
1467 
1468 void PlotHistEta(const char* infile,
1469                  std::string prefix,
1470                  std::string text,
1471                  int iene = 3,
1472                  int numb = 50,
1473                  int ieta = 0,
1474                  double lumi = 0,
1475                  double ener = 13.0,
1476                  bool isRealData = false,
1477                  bool drawStatBox = true,
1478                  int save = 0) {
1479   std::string name0 = "ratio";
1480   std::string title[5] = {"10:20", "20:30", "30:40", "40:60", "60:100"};
1481   std::string xtitl = "E_{HCAL}/(p-E_{ECAL})";
1482   std::string ytitl = "Tracks";
1483 
1484   gStyle->SetCanvasBorderMode(0);
1485   gStyle->SetCanvasColor(kWhite);
1486   gStyle->SetPadColor(kWhite);
1487   gStyle->SetFillColor(kWhite);
1488   gStyle->SetOptTitle(0);
1489   if (drawStatBox) {
1490     int iopt(1110);
1491     gStyle->SetOptStat(iopt);
1492     gStyle->SetOptFit(1);
1493   } else {
1494     gStyle->SetOptStat(0);
1495     gStyle->SetOptFit(0);
1496   }
1497   if (iene < 0 || iene >= 5)
1498     iene = 3;
1499   int numb2 = numb / 2;
1500   if (ieta < -numb2 || ieta > numb2)
1501     ieta = 0;
1502   int ietaMin = ((ieta == 0) ? 1 : ((ieta > 0) ? (numb2 + ieta) : (numb2 + ieta + 1)));
1503   int ietaMax = (ieta == 0) ? numb : ietaMin;
1504   TFile* file = new TFile(infile);
1505   char name[100], namep[100];
1506   for (int k = ietaMin; k <= ietaMax; ++k) {
1507     int eta = (k > numb2) ? (k - numb2) : (k - numb2 - 1);
1508     sprintf(name, "%s%s%d%d", prefix.c_str(), name0.c_str(), iene, k);
1509     TH1D* hist1 = (TH1D*)file->FindObjectAny(name);
1510     std::cout << name << " at " << hist1 << std::endl;
1511     if (hist1 != nullptr) {
1512       TH1D* hist = (TH1D*)(hist1->Clone());
1513       double ymin(0.90);
1514       sprintf(namep, "c_%s", name);
1515       TCanvas* pad = new TCanvas(namep, namep, 700, 500);
1516       pad->SetRightMargin(0.10);
1517       pad->SetTopMargin(0.10);
1518       hist->GetXaxis()->SetTitleSize(0.04);
1519       hist->GetXaxis()->SetTitle(xtitl.c_str());
1520       hist->GetYaxis()->SetTitle(ytitl.c_str());
1521       hist->GetYaxis()->SetLabelOffset(0.005);
1522       hist->GetYaxis()->SetTitleSize(0.04);
1523       hist->GetYaxis()->SetLabelSize(0.035);
1524       hist->GetYaxis()->SetTitleOffset(1.10);
1525       hist->GetXaxis()->SetRangeUser(0.25, 2.25);
1526       hist->SetMarkerStyle(20);
1527       hist->SetMarkerColor(2);
1528       hist->SetLineColor(2);
1529       hist->Draw();
1530       pad->Update();
1531       TPaveStats* st1 = (TPaveStats*)hist->GetListOfFunctions()->FindObject("stats");
1532       if (st1 != nullptr) {
1533         ymin = 0.70;
1534         st1->SetY1NDC(ymin);
1535         st1->SetY2NDC(0.90);
1536         st1->SetX1NDC(0.65);
1537         st1->SetX2NDC(0.90);
1538       }
1539       double ymx(0.96), xmi(0.25), xmx(0.90);
1540       char txt[100];
1541       if (lumi > 0.1) {
1542         ymx = ymin - 0.005;
1543         xmi = 0.45;
1544         TPaveText* txt0 = new TPaveText(0.65, 0.91, 0.90, 0.96, "blNDC");
1545         txt0->SetFillColor(0);
1546         sprintf(txt, "%4.1f TeV %5.1f fb^{-1}", ener, lumi);
1547         txt0->AddText(txt);
1548         txt0->Draw("same");
1549       }
1550       double ymi = ymx - 0.05;
1551       TPaveText* txt1 = new TPaveText(xmi, ymi, xmx, ymx, "blNDC");
1552       txt1->SetFillColor(0);
1553       if (text == "") {
1554         sprintf(txt, "Tracks with p = %s GeV at i#eta = %d", title[iene].c_str(), eta);
1555       } else {
1556         sprintf(txt, "Tracks with p = %s GeV at i#eta = %d (%s)", title[iene].c_str(), eta, text.c_str());
1557       }
1558       txt1->AddText(txt);
1559       txt1->Draw("same");
1560       double xmax = (isRealData) ? 0.33 : 0.44;
1561       ymi = (lumi > 0.1) ? 0.91 : 0.84;
1562       ymx = ymi + 0.05;
1563       TPaveText* txt2 = new TPaveText(0.11, ymi, xmax, ymx, "blNDC");
1564       txt2->SetFillColor(0);
1565       if (isRealData)
1566         sprintf(txt, "CMS Preliminary");
1567       else
1568         sprintf(txt, "CMS Simulation Preliminary");
1569       txt2->AddText(txt);
1570       txt2->Draw("same");
1571       pad->Modified();
1572       pad->Update();
1573       if (save > 0) {
1574         sprintf(name, "%s.pdf", pad->GetName());
1575         pad->Print(name);
1576       } else if (save < 0) {
1577         sprintf(name, "%s.C", pad->GetName());
1578         pad->Print(name);
1579       }
1580     }
1581   }
1582 }
1583 
1584 void PlotHists(std::string infile, std::string prefix, std::string text, bool drawStatBox = true, int save = 0) {
1585   int colors[6] = {1, 6, 4, 7, 2, 9};
1586   std::string types[6] = {"B", "C", "D", "E", "F", "G"};
1587   std::string names[3] = {"ratio20", "Z2", "W2"};
1588   std::string xtitl[3] = {"E_{HCAL}/(p-E_{ECAL})", "i#eta", "i#eta"};
1589   std::string ytitl[3] = {"Tracks", "MPV(E_{HCAL}/(p-E_{ECAL}))", "MPV/Width(E_{HCAL}/(p-E_{ECAL}))"};
1590 
1591   gStyle->SetCanvasBorderMode(0);
1592   gStyle->SetCanvasColor(kWhite);
1593   gStyle->SetPadColor(kWhite);
1594   gStyle->SetFillColor(kWhite);
1595   gStyle->SetOptTitle(0);
1596   if (drawStatBox)
1597     gStyle->SetOptFit(10);
1598   else
1599     gStyle->SetOptFit(0);
1600 
1601   char name[100], namep[100];
1602   TFile* file = new TFile(infile.c_str());
1603   for (int i = 0; i < 3; ++i) {
1604     std::vector<TH1D*> hists;
1605     std::vector<int> kks;
1606     double ymax(0.77);
1607     if (drawStatBox) {
1608       if (i == 0)
1609         gStyle->SetOptStat(1100);
1610       else
1611         gStyle->SetOptStat(10);
1612     } else {
1613       gStyle->SetOptStat(0);
1614       ymax = 0.89;
1615     }
1616     for (int k = 0; k < 6; ++k) {
1617       sprintf(name, "%s%s%s", prefix.c_str(), types[k].c_str(), names[i].c_str());
1618       TH1D* hist1 = (TH1D*)file->FindObjectAny(name);
1619       if (hist1 != nullptr) {
1620         hists.push_back((TH1D*)(hist1->Clone()));
1621         kks.push_back(k);
1622       }
1623     }
1624     if (hists.size() > 0) {
1625       sprintf(namep, "c_%s%s", prefix.c_str(), names[i].c_str());
1626       TCanvas* pad = new TCanvas(namep, namep, 700, 500);
1627       TLegend* legend = new TLegend(0.44, 0.89 - 0.055 * hists.size(), 0.69, ymax);
1628       legend->SetFillColor(kWhite);
1629       pad->SetRightMargin(0.10);
1630       pad->SetTopMargin(0.10);
1631       double ymax(0.90);
1632       double dy = (i == 0) ? 0.13 : 0.08;
1633       for (unsigned int jk = 0; jk < hists.size(); ++jk) {
1634         int k = kks[jk];
1635         hists[jk]->GetXaxis()->SetTitle(xtitl[i].c_str());
1636         hists[jk]->GetYaxis()->SetTitle(ytitl[i].c_str());
1637         hists[jk]->GetYaxis()->SetLabelOffset(0.005);
1638         hists[jk]->GetYaxis()->SetLabelSize(0.035);
1639         hists[jk]->GetYaxis()->SetTitleOffset(1.15);
1640         if (i == 0) {
1641           hists[jk]->GetXaxis()->SetRangeUser(0.0, 2.5);
1642         } else if (i == 1) {
1643           hists[jk]->GetYaxis()->SetRangeUser(0.5, 2.0);
1644         } else {
1645           hists[jk]->GetYaxis()->SetRangeUser(0.0, 0.5);
1646         }
1647         hists[jk]->SetMarkerStyle(20);
1648         hists[jk]->SetMarkerColor(colors[k]);
1649         hists[jk]->SetLineColor(colors[k]);
1650         if (jk == 0)
1651           hists[jk]->Draw();
1652         else
1653           hists[jk]->Draw("sames");
1654         pad->Update();
1655         TPaveStats* st1 = (TPaveStats*)hists[jk]->GetListOfFunctions()->FindObject("stats");
1656         if (st1 != nullptr) {
1657           double ymin = ymax - dy;
1658           st1->SetLineColor(colors[k]);
1659           st1->SetTextColor(colors[k]);
1660           st1->SetY1NDC(ymin);
1661           st1->SetY2NDC(ymax);
1662           st1->SetX1NDC(0.70);
1663           st1->SetX2NDC(0.90);
1664           ymax = ymin;
1665         }
1666         sprintf(name, "%s%s", text.c_str(), types[k].c_str());
1667         legend->AddEntry(hists[jk], name, "lp");
1668       }
1669       legend->Draw("same");
1670       pad->Update();
1671       TPaveText* txt1 = new TPaveText(0.34, 0.825, 0.69, 0.895, "blNDC");
1672       txt1->SetFillColor(0);
1673       char txt[100];
1674       sprintf(txt, "Tracks with p = 40:60 GeV");
1675       txt1->AddText(txt);
1676       txt1->Draw("same");
1677       TPaveText* txt2 = new TPaveText(0.11, 0.825, 0.33, 0.895, "blNDC");
1678       txt2->SetFillColor(0);
1679       sprintf(txt, "CMS Preliminary");
1680       txt2->AddText(txt);
1681       txt2->Draw("same");
1682       if (!drawStatBox && i == 1) {
1683         double xmin = hists[0]->GetBinLowEdge(1);
1684         int nbin = hists[0]->GetNbinsX();
1685         double xmax = hists[0]->GetBinLowEdge(nbin) + hists[0]->GetBinWidth(nbin);
1686         TLine line = TLine(xmin, 1.0, xmax, 1.0);  //etamin,1.0,etamax,1.0);
1687         line.SetLineWidth(4);
1688         line.Draw("same");
1689         pad->Update();
1690       }
1691       pad->Modified();
1692       pad->Update();
1693       if (save > 0) {
1694         sprintf(name, "%s.pdf", pad->GetName());
1695         pad->Print(name);
1696       } else if (save < 0) {
1697         sprintf(name, "%s.C", pad->GetName());
1698         pad->Print(name);
1699       }
1700     }
1701   }
1702 }
1703 
1704 void PlotTwoHists(std::string infile,
1705                   std::string prefix1,
1706                   std::string text1,
1707                   std::string prefix2,
1708                   std::string text2,
1709                   std::string text0,
1710                   int type = 0,
1711                   int iname = 3,
1712                   double lumi = 0,
1713                   double ener = 13.0,
1714                   int drawStatBox = 0,
1715                   int save = 0) {
1716   int colors[2] = {2, 4};
1717   int numb[2] = {5, 1};
1718   std::string names0[5] = {"ratio00", "ratio00One", "etaB04One", "Z0", "W0"};
1719   std::string names1[5] = {"ratio10", "ratio10One", "etaB14One", "Z1", "W1"};
1720   std::string names2[5] = {"ratio30", "ratio30One", "etaB34One", "Z3", "W3"};
1721   std::string xtitl1[5] = {"E_{HCAL}/(p-E_{ECAL})", "E_{HCAL}/(p-E_{ECAL})", "E_{HCAL}/(p-E_{ECAL})", "i#eta", "i#eta"};
1722   std::string ytitl1[5] = {
1723       "Tracks", "Tracks", "Tracks", "MPV(E_{HCAL}/(p-E_{ECAL}))", "MPV/Width(E_{HCAL}/(p-E_{ECAL}))"};
1724   std::string names3[1] = {"R"};
1725   std::string xtitl2[1] = {"RBX #"};
1726   std::string ytitl2[1] = {"MPV(E_{HCAL}/(p-E_{ECAL}))"};
1727 
1728   gStyle->SetCanvasBorderMode(0);
1729   gStyle->SetCanvasColor(kWhite);
1730   gStyle->SetPadColor(kWhite);
1731   gStyle->SetFillColor(kWhite);
1732   gStyle->SetOptTitle(0);
1733   if ((drawStatBox / 10) % 10 > 0)
1734     gStyle->SetOptFit(10);
1735   else
1736     gStyle->SetOptFit(0);
1737 
1738   if (type != 1)
1739     type = 0;
1740   char name[100], namep[100];
1741   TFile* file = new TFile(infile.c_str());
1742   for (int i = 0; i < numb[type]; ++i) {
1743     std::vector<TH1D*> hists;
1744     std::vector<int> kks;
1745     double ymax(0.77);
1746     if (drawStatBox % 10 > 0) {
1747       if (i != 2)
1748         gStyle->SetOptStat(1100);
1749       else
1750         gStyle->SetOptStat(10);
1751     } else {
1752       gStyle->SetOptStat(0);
1753       ymax = 0.82;
1754     }
1755     for (int k = 0; k < 2; ++k) {
1756       std::string prefix = (k == 0) ? prefix1 : prefix2;
1757       if (type == 0) {
1758         if (iname == 0)
1759           sprintf(name, "%s%s", prefix.c_str(), names0[i].c_str());
1760         else if (iname == 1)
1761           sprintf(name, "%s%s", prefix.c_str(), names1[i].c_str());
1762         else
1763           sprintf(name, "%s%s", prefix.c_str(), names2[i].c_str());
1764       } else {
1765         sprintf(name, "%s%s%d", prefix.c_str(), names3[i].c_str(), iname);
1766       }
1767       TH1D* hist1 = (TH1D*)file->FindObjectAny(name);
1768       if (hist1 != nullptr) {
1769         hists.push_back((TH1D*)(hist1->Clone()));
1770         kks.push_back(k);
1771       }
1772     }
1773     if (hists.size() == 2) {
1774       if (type == 0) {
1775         if (iname == 0)
1776           sprintf(namep, "c_%s%s%s", prefix1.c_str(), prefix2.c_str(), names0[i].c_str());
1777         else if (iname == 1)
1778           sprintf(namep, "c_%s%s%s", prefix1.c_str(), prefix2.c_str(), names1[i].c_str());
1779         else
1780           sprintf(namep, "c_%s%s%s", prefix1.c_str(), prefix2.c_str(), names2[i].c_str());
1781       } else {
1782         sprintf(namep, "c_%s%s%s%d", prefix1.c_str(), prefix2.c_str(), names3[i].c_str(), iname);
1783       }
1784       double ymax(0.90);
1785       double dy = (i == 0) ? 0.13 : 0.08;
1786       double ymx0 = (drawStatBox == 0) ? (ymax - .01) : (ymax - dy * hists.size() - .01);
1787       TCanvas* pad = new TCanvas(namep, namep, 700, 500);
1788       TLegend* legend = new TLegend(0.64, ymx0 - 0.05 * hists.size(), 0.89, ymx0);
1789       legend->SetFillColor(kWhite);
1790       pad->SetRightMargin(0.10);
1791       pad->SetTopMargin(0.10);
1792       for (unsigned int jk = 0; jk < hists.size(); ++jk) {
1793         int k = kks[jk];
1794         hists[jk]->GetXaxis()->SetTitleSize(0.040);
1795         if (type == 0) {
1796           hists[jk]->GetXaxis()->SetTitle(xtitl1[i].c_str());
1797           hists[jk]->GetYaxis()->SetTitle(ytitl1[i].c_str());
1798         } else {
1799           hists[jk]->GetXaxis()->SetTitle(xtitl2[i].c_str());
1800           hists[jk]->GetYaxis()->SetTitle(ytitl2[i].c_str());
1801         }
1802         hists[jk]->GetYaxis()->SetLabelOffset(0.005);
1803         hists[jk]->GetYaxis()->SetLabelSize(0.035);
1804         hists[jk]->GetYaxis()->SetTitleSize(0.040);
1805         hists[jk]->GetYaxis()->SetTitleOffset(1.15);
1806         if ((type == 0) && (i != 3) && (i != 4))
1807           hists[jk]->GetXaxis()->SetRangeUser(0.0, 5.0);
1808         if (type == 0) {
1809           if (i == 3)
1810             hists[jk]->GetYaxis()->SetRangeUser(0.8, 1.2);
1811           else if (i == 4)
1812             hists[jk]->GetYaxis()->SetRangeUser(0.0, 0.5);
1813         }
1814         if (type != 0)
1815           hists[jk]->GetYaxis()->SetRangeUser(0.75, 1.2);
1816         hists[jk]->SetMarkerStyle(20);
1817         hists[jk]->SetMarkerColor(colors[k]);
1818         hists[jk]->SetLineColor(colors[k]);
1819         if (jk == 0)
1820           hists[jk]->Draw();
1821         else
1822           hists[jk]->Draw("sames");
1823         pad->Update();
1824         TPaveStats* st1 = (TPaveStats*)hists[jk]->GetListOfFunctions()->FindObject("stats");
1825         if (st1 != nullptr) {
1826           double ymin = ymax - dy;
1827           st1->SetLineColor(colors[k]);
1828           st1->SetTextColor(colors[k]);
1829           st1->SetY1NDC(ymin);
1830           st1->SetY2NDC(ymax);
1831           st1->SetX1NDC(0.70);
1832           st1->SetX2NDC(0.90);
1833           ymax = ymin;
1834         }
1835         if (k == 0)
1836           sprintf(name, "%s", text1.c_str());
1837         else
1838           sprintf(name, "%s", text2.c_str());
1839         legend->AddEntry(hists[jk], name, "lp");
1840       }
1841       legend->Draw("same");
1842       pad->Update();
1843       char txt[100];
1844       double xmi(0.10), xmx(0.895), ymx(0.95);
1845       if (lumi > 0.01) {
1846         xmx = 0.70;
1847         xmi = 0.30;
1848         TPaveText* txt0 = new TPaveText(0.705, 0.905, 0.90, 0.95, "blNDC");
1849         txt0->SetFillColor(0);
1850         sprintf(txt, "%4.1f TeV %5.1f fb^{-1}", ener, lumi);
1851         txt0->AddText(txt);
1852         txt0->Draw("same");
1853       }
1854       double ymi = ymx - 0.045;
1855       TPaveText* txt1 = new TPaveText(xmi, ymi, xmx, ymx, "blNDC");
1856       txt1->SetFillColor(0);
1857       if (iname == 0)
1858         sprintf(txt, "p = 20:30 GeV %s", text0.c_str());
1859       else
1860         sprintf(txt, "p = 40:60 GeV %s", text0.c_str());
1861       txt1->AddText(txt);
1862       txt1->Draw("same");
1863       ymi = (lumi > 0.1) ? 0.905 : 0.85;
1864       ymx = ymi + 0.045;
1865       TPaveText* txt2 = new TPaveText(0.105, ymi, 0.295, ymx, "blNDC");
1866       txt2->SetFillColor(0);
1867       sprintf(txt, "CMS Preliminary");
1868       txt2->AddText(txt);
1869       txt2->Draw("same");
1870       pad->Modified();
1871       if ((drawStatBox == 0) && (i == 3)) {
1872         double xmin = hists[0]->GetBinLowEdge(1);
1873         int nbin = hists[0]->GetNbinsX();
1874         double xmax = hists[0]->GetBinLowEdge(nbin) + hists[0]->GetBinWidth(nbin);
1875         TLine* line = new TLine(xmin, 1.0, xmax, 1.0);  //etamin,1.0,etamax,1.0);
1876         line->SetLineWidth(2);
1877         line->Draw("same");
1878         pad->Update();
1879       }
1880       pad->Update();
1881       if (save > 0) {
1882         sprintf(name, "%s.pdf", pad->GetName());
1883         pad->Print(name);
1884       } else if (save < 0) {
1885         sprintf(name, "%s.C", pad->GetName());
1886         pad->Print(name);
1887       }
1888     }
1889   }
1890 }
1891 
1892 void PlotFiveHists(std::string infile,
1893                    std::string text0,
1894                    std::string prefix0,
1895                    int type = 0,
1896                    int iname = 3,
1897                    int drawStatBox = 0,
1898                    bool normalize = false,
1899                    int save = 0,
1900                    std::string prefix1 = "",
1901                    std::string text1 = "",
1902                    std::string prefix2 = "",
1903                    std::string text2 = "",
1904                    std::string prefix3 = "",
1905                    std::string text3 = "",
1906                    std::string prefix4 = "",
1907                    std::string text4 = "",
1908                    std::string prefix5 = "",
1909                    std::string text5 = "") {
1910   int colors[5] = {2, 4, 6, 1, 7};
1911   int numb[3] = {5, 1, 4};
1912   std::string names0[5] = {"ratio00", "ratio00One", "etaB04", "Z0", "W0"};
1913   std::string names1[5] = {"ratio10", "ratio10One", "etaB14", "Z1", "W1"};
1914   std::string names2[5] = {"ratio30", "ratio30One", "etaB34", "Z3", "W3"};
1915   std::string xtitl1[5] = {"E_{HCAL}/(p-E_{ECAL})", "E_{HCAL}/(p-E_{ECAL})", "E_{HCAL}/(p-E_{ECAL})", "i#eta", "i#eta"};
1916   std::string ytitl1[5] = {
1917       "Tracks", "Tracks", "Tracks", "MPV(E_{HCAL}/(p-E_{ECAL}))", "MPV/Width(E_{HCAL}/(p-E_{ECAL}))"};
1918   std::string names3[1] = {"R"};
1919   std::string xtitl2[1] = {"RBX #"};
1920   std::string ytitl2[1] = {"MPV(E_{HCAL}/(p-E_{ECAL}))"};
1921   std::string names4[4] = {"pp21", "pp22", "pp23", "pp24"};
1922   std::string xtitl3[4] = {"p (GeV)", "p (GeV)", "p (GeV)", "p (GeV)"};
1923   std::string ytitl3[4] = {"Tracks", "Tracks", "Tracks", "Tracks"};
1924   std::string title3[4] = {"Barrel", "Transition", "Endcap", "Combined"};
1925 
1926   gStyle->SetCanvasBorderMode(0);
1927   gStyle->SetCanvasColor(kWhite);
1928   gStyle->SetPadColor(kWhite);
1929   gStyle->SetFillColor(kWhite);
1930   gStyle->SetOptTitle(0);
1931   if ((drawStatBox / 10) % 10 > 0)
1932     gStyle->SetOptFit(10);
1933   else
1934     gStyle->SetOptFit(0);
1935 
1936   if (type != 1 && type != 2)
1937     type = 0;
1938   char name[100], namep[100];
1939   TFile* file = new TFile(infile.c_str());
1940   for (int i = 0; i < numb[type]; ++i) {
1941     std::vector<TH1D*> hists;
1942     std::vector<int> kks;
1943     std::vector<std::string> texts;
1944     double ymax(0.77);
1945     if (drawStatBox % 10 > 0) {
1946       if (type == 2)
1947         gStyle->SetOptStat(1110);
1948       else if (i != 3)
1949         gStyle->SetOptStat(1100);
1950       else
1951         gStyle->SetOptStat(10);
1952     } else {
1953       gStyle->SetOptStat(0);
1954       ymax = 0.82;
1955     }
1956     for (int k = 0; k < 5; ++k) {
1957       std::string prefix, text;
1958       if (k == 0) {
1959         prefix = prefix1;
1960         text = text1;
1961       } else if (k == 1) {
1962         prefix = prefix2;
1963         text = text2;
1964       } else if (k == 2) {
1965         prefix = prefix3;
1966         text = text3;
1967       } else if (k == 3) {
1968         prefix = prefix4;
1969         text = text4;
1970       } else {
1971         prefix = prefix5;
1972         text = text5;
1973       }
1974       if (prefix != "") {
1975         if (type == 0) {
1976           if (iname == 0)
1977             sprintf(name, "%s%s", prefix.c_str(), names0[i].c_str());
1978           else if (iname == 1)
1979             sprintf(name, "%s%s", prefix.c_str(), names1[i].c_str());
1980           else
1981             sprintf(name, "%s%s", prefix.c_str(), names2[i].c_str());
1982         } else if (type == 1) {
1983           sprintf(name, "%s%s%d", prefix.c_str(), names3[i].c_str(), iname);
1984         } else {
1985           sprintf(name, "%s%s", prefix.c_str(), names4[i].c_str());
1986         }
1987         TH1D* hist1 = (TH1D*)file->FindObjectAny(name);
1988         if (hist1 != nullptr) {
1989           hists.push_back((TH1D*)(hist1->Clone()));
1990           kks.push_back(k);
1991           texts.push_back(text);
1992         }
1993       }
1994     }
1995     if (hists.size() > 0) {
1996       if (type == 0) {
1997         if (iname == 0)
1998           sprintf(namep, "c_%s%s", prefix0.c_str(), names0[i].c_str());
1999         else if (iname == 1)
2000           sprintf(namep, "c_%s%s", prefix0.c_str(), names1[i].c_str());
2001         else
2002           sprintf(namep, "c_%s%s", prefix0.c_str(), names2[i].c_str());
2003       } else if (type == 1) {
2004         sprintf(namep, "c_%s%s%d", prefix0.c_str(), names3[i].c_str(), iname);
2005       } else {
2006         sprintf(namep, "c_%s%s", prefix0.c_str(), names4[i].c_str());
2007       }
2008       double ymax(0.90);
2009       double dy = (i == 0 && type == 0) ? 0.13 : 0.08;
2010       double ymx0 = (drawStatBox == 0) ? (ymax - .01) : (ymax - dy * hists.size() - .01);
2011       TCanvas* pad = new TCanvas(namep, namep, 700, 500);
2012       TLegend* legend = new TLegend(0.64, ymx0 - 0.05 * hists.size(), 0.89, ymx0);
2013       legend->SetFillColor(kWhite);
2014       pad->SetRightMargin(0.10);
2015       pad->SetTopMargin(0.10);
2016       for (unsigned int jk = 0; jk < hists.size(); ++jk) {
2017         int k = kks[jk];
2018         hists[jk]->GetXaxis()->SetTitleSize(0.040);
2019         if (type == 0) {
2020           hists[jk]->GetXaxis()->SetTitle(xtitl1[i].c_str());
2021           hists[jk]->GetYaxis()->SetTitle(ytitl1[i].c_str());
2022         } else if (type == 1) {
2023           hists[jk]->GetXaxis()->SetTitle(xtitl2[i].c_str());
2024           hists[jk]->GetYaxis()->SetTitle(ytitl2[i].c_str());
2025         } else {
2026           hists[jk]->GetXaxis()->SetTitle(xtitl3[i].c_str());
2027           hists[jk]->GetYaxis()->SetTitle(ytitl3[i].c_str());
2028         }
2029         hists[jk]->GetYaxis()->SetLabelOffset(0.005);
2030         hists[jk]->GetYaxis()->SetLabelSize(0.035);
2031         hists[jk]->GetYaxis()->SetTitleSize(0.040);
2032         hists[jk]->GetYaxis()->SetTitleOffset(1.15);
2033         if ((type == 0) && (i != 3) && (i != 4))
2034           hists[jk]->GetXaxis()->SetRangeUser(0.0, 2.5);
2035         if (type == 0) {
2036           if (i == 3)
2037             hists[jk]->GetYaxis()->SetRangeUser(0.8, 1.2);
2038           else if (i == 4)
2039             hists[jk]->GetYaxis()->SetRangeUser(0.0, 0.5);
2040         }
2041         if (type == 1)
2042           hists[jk]->GetYaxis()->SetRangeUser(0.75, 1.2);
2043         hists[jk]->SetMarkerStyle(20);
2044         hists[jk]->SetMarkerColor(colors[k]);
2045         hists[jk]->SetLineColor(colors[k]);
2046         if (normalize && ((type == 2) || ((type == 0) && (i != 3)))) {
2047           if (jk == 0)
2048             hists[jk]->DrawNormalized();
2049           else
2050             hists[jk]->DrawNormalized("sames");
2051         } else {
2052           if (jk == 0)
2053             hists[jk]->Draw();
2054           else
2055             hists[jk]->Draw("sames");
2056         }
2057         pad->Update();
2058         TPaveStats* st1 = (TPaveStats*)hists[jk]->GetListOfFunctions()->FindObject("stats");
2059         if (st1 != nullptr) {
2060           double ymin = ymax - dy;
2061           st1->SetLineColor(colors[k]);
2062           st1->SetTextColor(colors[k]);
2063           st1->SetY1NDC(ymin);
2064           st1->SetY2NDC(ymax);
2065           st1->SetX1NDC(0.70);
2066           st1->SetX2NDC(0.90);
2067           ymax = ymin;
2068         }
2069         sprintf(name, "%s", texts[jk].c_str());
2070         legend->AddEntry(hists[jk], name, "lp");
2071       }
2072       legend->Draw("same");
2073       pad->Update();
2074       TPaveText* txt1 = new TPaveText(0.10, 0.905, 0.80, 0.95, "blNDC");
2075       txt1->SetFillColor(0);
2076       char txt[100];
2077       if (type == 2) {
2078         sprintf(txt, "p = 40:60 GeV (%s)", title3[i].c_str());
2079       } else if (((type == 0) && (iname == 0))) {
2080         sprintf(txt, "p = 20:30 GeV %s", text0.c_str());
2081       } else {
2082         sprintf(txt, "p = 40:60 GeV %s", text0.c_str());
2083       }
2084       txt1->AddText(txt);
2085       txt1->Draw("same");
2086       TPaveText* txt2 = new TPaveText(0.11, 0.825, 0.33, 0.895, "blNDC");
2087       txt2->SetFillColor(0);
2088       sprintf(txt, "CMS Preliminary");
2089       txt2->AddText(txt);
2090       txt2->Draw("same");
2091       pad->Modified();
2092       if ((drawStatBox == 0) && (i == 3)) {
2093         double xmin = hists[0]->GetBinLowEdge(1);
2094         int nbin = hists[0]->GetNbinsX();
2095         double xmax = hists[0]->GetBinLowEdge(nbin) + hists[0]->GetBinWidth(nbin);
2096         TLine* line = new TLine(xmin, 1.0, xmax, 1.0);  //etamin,1.0,etamax,1.0);
2097         line->SetLineWidth(2);
2098         line->Draw("same");
2099         pad->Update();
2100       }
2101       pad->Update();
2102       if (save > 0) {
2103         sprintf(name, "%s.pdf", pad->GetName());
2104         pad->Print(name);
2105       } else if (save < 0) {
2106         sprintf(name, "%s.C", pad->GetName());
2107         pad->Print(name);
2108       }
2109     }
2110   }
2111 }
2112 
2113 void PlotHistCorrResults(std::string infile, std::string text, std::string prefixF, int save = 0) {
2114   std::string name[5] = {"Eta1Bf", "Eta2Bf", "Eta1Af", "Eta2Af", "Cvg0"};
2115   std::string title[5] = {"Mean at the start of itertions",
2116                           "Median at the start of itertions",
2117                           "Mean at the end of itertions",
2118                           "Median at the end of itertions",
2119                           ""};
2120   int type[5] = {0, 0, 0, 0, 1};
2121 
2122   gStyle->SetCanvasBorderMode(0);
2123   gStyle->SetCanvasColor(kWhite);
2124   gStyle->SetPadColor(kWhite);
2125   gStyle->SetFillColor(kWhite);
2126   gStyle->SetOptTitle(0);
2127   gStyle->SetOptStat(10);
2128   gStyle->SetOptFit(10);
2129   TFile* file = new TFile(infile.c_str());
2130   char namep[100];
2131   for (int k = 0; k < 5; ++k) {
2132     TH1D* hist1 = (TH1D*)file->FindObjectAny(name[k].c_str());
2133     if (hist1 != nullptr) {
2134       TH1D* hist = (TH1D*)(hist1->Clone());
2135       sprintf(namep, "c_%s%s", prefixF.c_str(), name[k].c_str());
2136       TCanvas* pad = new TCanvas(namep, namep, 700, 500);
2137       pad->SetRightMargin(0.10);
2138       pad->SetTopMargin(0.10);
2139       hist->GetYaxis()->SetLabelOffset(0.005);
2140       hist->GetYaxis()->SetTitleOffset(1.20);
2141       double xmin = hist->GetBinLowEdge(1);
2142       int nbin = hist->GetNbinsX();
2143       double xmax = hist->GetBinLowEdge(nbin) + hist->GetBinWidth(nbin);
2144       std::cout << hist->GetTitle() << " Bins " << nbin << ":" << xmin << ":" << xmax << std::endl;
2145       double xlow(0.12), ylow(0.82);
2146       char txt[100], option[2];
2147       if (type[k] == 0) {
2148         sprintf(namep, "f_%s", name[k].c_str());
2149         TF1* func = new TF1(namep, "pol0", xmin, xmax);
2150         hist->Fit(func, "+QWL", "");
2151         if (text == "")
2152           sprintf(txt, "%s", title[k].c_str());
2153         else
2154           sprintf(txt, "%s (balancing the %s)", title[k].c_str(), text.c_str());
2155         sprintf(option, "%s", "");
2156       } else {
2157         hist->GetXaxis()->SetTitle("Iterations");
2158         hist->GetYaxis()->SetTitle("Convergence");
2159         hist->SetMarkerStyle(20);
2160         hist->SetMarkerColor(2);
2161         hist->SetMarkerSize(0.8);
2162         xlow = 0.50;
2163         ylow = 0.86;
2164         sprintf(txt, "(%s)", text.c_str());
2165         sprintf(option, "%s", "p");
2166       }
2167       hist->Draw(option);
2168       pad->Update();
2169       TPaveStats* st1 = (TPaveStats*)hist->GetListOfFunctions()->FindObject("stats");
2170       if (st1 != nullptr) {
2171         st1->SetY1NDC(ylow);
2172         st1->SetY2NDC(0.90);
2173         st1->SetX1NDC(0.70);
2174         st1->SetX2NDC(0.90);
2175       }
2176       TPaveText* txt1 = new TPaveText(xlow, 0.82, 0.68, 0.88, "blNDC");
2177       txt1->SetFillColor(0);
2178       txt1->AddText(txt);
2179       txt1->Draw("same");
2180       pad->Modified();
2181       pad->Update();
2182       if (save > 0) {
2183         sprintf(namep, "%s.pdf", pad->GetName());
2184         pad->Print(namep);
2185       } else if (save < 0) {
2186         sprintf(namep, "%s.C", pad->GetName());
2187         pad->Print(namep);
2188       }
2189     }
2190   }
2191 }
2192 
2193 void PlotHistCorrFactor(char* infile,
2194                         std::string text,
2195                         std::string prefixF,
2196                         double scale = 1.0,
2197                         int nmin = 100,
2198                         bool isRealData = true,
2199                         bool drawStatBox = false,
2200                         int iformat = 0,
2201                         int save = 0) {
2202   std::map<int, cfactors> cfacs;
2203   int etamin(100), etamax(-100), maxdepth(0);
2204   readCorrFactors(infile, scale, cfacs, etamin, etamax, maxdepth, iformat);
2205 
2206   gStyle->SetCanvasBorderMode(0);
2207   gStyle->SetCanvasColor(kWhite);
2208   gStyle->SetPadColor(kWhite);
2209   gStyle->SetFillColor(kWhite);
2210   gStyle->SetOptTitle(0);
2211   if (drawStatBox) {
2212     gStyle->SetOptStat(10);
2213     gStyle->SetOptFit(10);
2214   } else {
2215     gStyle->SetOptStat(0);
2216     gStyle->SetOptFit(0);
2217   }
2218   int colors[7] = {1, 6, 4, 7, 2, 9, 3};
2219   int mtype[7] = {20, 21, 22, 23, 24, 33, 25};
2220   int nbin = etamax - etamin + 1;
2221   std::vector<TH1D*> hists;
2222   std::vector<int> entries;
2223   char name[100];
2224   double dy(0);
2225   int fits(0);
2226   for (int j = 0; j < maxdepth; ++j) {
2227     sprintf(name, "hd%d", j + 1);
2228     TObject* ob = gROOT->FindObject(name);
2229     if (ob)
2230       ob->Delete();
2231     TH1D* h = new TH1D(name, name, nbin, etamin, etamax);
2232     int nent(0);
2233     double chi2(0);
2234     for (std::map<int, cfactors>::const_iterator itr = cfacs.begin(); itr != cfacs.end(); ++itr) {
2235       if ((itr->second).depth == j + 1) {
2236         int ieta = (itr->second).ieta;
2237         int bin = ieta - etamin + 1;
2238         float val = (itr->second).corrf;
2239         float dvl = (itr->second).dcorr;
2240         h->SetBinContent(bin, val);
2241         h->SetBinError(bin, dvl);
2242         nent++;
2243         chi2 += (((val - 1.0) / dvl) * ((val - 1.0) / dvl));
2244       }
2245     }
2246     std::cout << "Depth = " << (j + 1) << " chi2 = " << chi2 << "/" << nent << std::endl;
2247     if (nent > nmin) {
2248       fits++;
2249       dy += 0.025;
2250       sprintf(name, "hdf%d", j + 1);
2251       TObject* ob = gROOT->FindObject(name);
2252       if (ob)
2253         ob->Delete();
2254       TF1* func = new TF1(name, "pol0", etamin, etamax);
2255       h->Fit(func, "+QWLR", "");
2256     }
2257     h->SetLineColor(colors[j]);
2258     h->SetMarkerColor(colors[j]);
2259     h->SetMarkerStyle(mtype[j]);
2260     h->GetXaxis()->SetTitle("i#eta");
2261     h->GetYaxis()->SetTitle("Correction Factor");
2262     h->GetYaxis()->SetLabelOffset(0.005);
2263     h->GetYaxis()->SetTitleOffset(1.20);
2264     h->GetYaxis()->SetRangeUser(0.0, 2.0);
2265     hists.push_back(h);
2266     entries.push_back(nent);
2267     dy += 0.025;
2268   }
2269   sprintf(name, "c_%sCorrFactor", prefixF.c_str());
2270   TCanvas* pad = new TCanvas(name, name, 700, 500);
2271   pad->SetRightMargin(0.10);
2272   pad->SetTopMargin(0.10);
2273   double yh = 0.90;
2274   // double yl = yh - 0.025 * hists.size() - dy - 0.01;
2275   double yl = 0.15;
2276   TLegend* legend = new TLegend(0.35, yl, 0.65, yl + 0.04 * hists.size());
2277   legend->SetFillColor(kWhite);
2278   for (unsigned int k = 0; k < hists.size(); ++k) {
2279     if (k == 0)
2280       hists[k]->Draw("");
2281     else
2282       hists[k]->Draw("sames");
2283     pad->Update();
2284     TPaveStats* st1 = (TPaveStats*)hists[k]->GetListOfFunctions()->FindObject("stats");
2285     if (st1 != nullptr) {
2286       dy = (entries[k] > nmin) ? 0.05 : 0.025;
2287       st1->SetLineColor(colors[k]);
2288       st1->SetTextColor(colors[k]);
2289       st1->SetY1NDC(yh - dy);
2290       st1->SetY2NDC(yh);
2291       st1->SetX1NDC(0.70);
2292       st1->SetX2NDC(0.90);
2293       yh -= dy;
2294     }
2295     sprintf(name, "Depth %d (%s)", k + 1, text.c_str());
2296     legend->AddEntry(hists[k], name, "lp");
2297   }
2298   legend->Draw("same");
2299   pad->Update();
2300   if (fits < 1) {
2301     double xmin = hists[0]->GetBinLowEdge(1);
2302     int nbin = hists[0]->GetNbinsX();
2303     double xmax = hists[0]->GetBinLowEdge(nbin) + hists[0]->GetBinWidth(nbin);
2304     TLine* line = new TLine(xmin, 1.0, xmax, 1.0);
2305     line->SetLineColor(9);
2306     line->SetLineWidth(2);
2307     line->SetLineStyle(2);
2308     line->Draw("same");
2309     pad->Modified();
2310     pad->Update();
2311   }
2312   char txt1[30];
2313   double xmax = (isRealData) ? 0.33 : 0.44;
2314   TPaveText* txt2 = new TPaveText(0.11, 0.85, xmax, 0.89, "blNDC");
2315   txt2->SetFillColor(0);
2316   if (isRealData)
2317     sprintf(txt1, "CMS Preliminary");
2318   else
2319     sprintf(txt1, "CMS Simulation Preliminary");
2320   txt2->AddText(txt1);
2321   txt2->Draw("same");
2322   pad->Modified();
2323   pad->Update();
2324   if (save > 0) {
2325     sprintf(name, "%s.pdf", pad->GetName());
2326     pad->Print(name);
2327   } else if (save < 0) {
2328     sprintf(name, "%s.C", pad->GetName());
2329     pad->Print(name);
2330   }
2331 }
2332 
2333 void PlotHistCorrFactor(char* infile,
2334                         std::string text,
2335                         int depth,
2336                         std::string prefixF,
2337                         double scale = 1.0,
2338                         int nmin = 100,
2339                         bool isRealData = true,
2340                         bool drawStatBox = false,
2341                         int iformat = 0,
2342                         int save = 0) {
2343   std::map<int, cfactors> cfacs;
2344   int etamin(100), etamax(-100), maxdepth(0);
2345   readCorrFactors(infile, scale, cfacs, etamin, etamax, maxdepth, iformat);
2346 
2347   gStyle->SetCanvasBorderMode(0);
2348   gStyle->SetCanvasColor(kWhite);
2349   gStyle->SetPadColor(kWhite);
2350   gStyle->SetFillColor(kWhite);
2351   gStyle->SetOptTitle(0);
2352   if (drawStatBox) {
2353     gStyle->SetOptStat(10);
2354     gStyle->SetOptFit(10);
2355   } else {
2356     gStyle->SetOptStat(0);
2357     gStyle->SetOptFit(0);
2358   }
2359   int colors[7] = {1, 6, 4, 7, 2, 9, 3};
2360   int mtype[7] = {20, 21, 22, 23, 24, 33, 25};
2361   int nbin = etamax - etamin + 1;
2362   std::vector<TH1D*> hists;
2363   std::vector<int> entries;
2364   char name[100];
2365   double dy(0);
2366   int fits(0);
2367   sprintf(name, "hd%d", depth);
2368   TObject* ob = gROOT->FindObject(name);
2369   if (ob)
2370     ob->Delete();
2371   TH1D* h = new TH1D(name, name, nbin, etamin, etamax);
2372   int nent(0);
2373   double chi2(0);
2374   for (std::map<int, cfactors>::const_iterator itr = cfacs.begin(); itr != cfacs.end(); ++itr) {
2375     if ((itr->second).depth == depth) {
2376       int ieta = (itr->second).ieta;
2377       int bin = ieta - etamin + 1;
2378       float val = (itr->second).corrf;
2379       float dvl = (itr->second).dcorr;
2380       h->SetBinContent(bin, val);
2381       h->SetBinError(bin, dvl);
2382       nent++;
2383       chi2 += (((val - 1.0) / dvl) * ((val - 1.0) / dvl));
2384     }
2385   }
2386   std::cout << "Depth = " << depth << " chi2 = " << chi2 << "/" << nent << std::endl;
2387   if (nent > nmin) {
2388     fits++;
2389     dy += 0.025;
2390     sprintf(name, "hdf%d", depth);
2391     TObject* ob = gROOT->FindObject(name);
2392     if (ob)
2393       ob->Delete();
2394     TF1* func = new TF1(name, "pol0", etamin, etamax);
2395     h->Fit(func, "+QWLR", "");
2396   }
2397   h->SetLineColor(colors[depth - 1]);
2398   h->SetMarkerColor(colors[depth - 1]);
2399   h->SetMarkerStyle(mtype[depth - 1]);
2400   h->GetXaxis()->SetTitle("i#eta");
2401   h->GetYaxis()->SetTitle("Correction Factor");
2402   h->GetYaxis()->SetLabelOffset(0.005);
2403   h->GetYaxis()->SetTitleOffset(1.20);
2404   h->GetYaxis()->SetRangeUser(0.0, 2.0);
2405   hists.push_back(h);
2406   entries.push_back(nent);
2407   dy += 0.025;
2408 
2409   sprintf(name, "c_%sD%dCorrFactor", prefixF.c_str(), depth);
2410   TCanvas* pad = new TCanvas(name, name, 700, 500);
2411   pad->SetRightMargin(0.10);
2412   pad->SetTopMargin(0.10);
2413   double yh = 0.90;
2414   // double yl = yh - 0.025 * hists.size() - dy - 0.01;
2415   double yl = 0.15;
2416   TLegend* legend = new TLegend(0.35, yl, 0.85, yl + 0.04 * hists.size());
2417   legend->SetFillColor(kWhite);
2418   for (unsigned int k = 0; k < hists.size(); ++k) {
2419     if (k == 0)
2420       hists[k]->Draw("");
2421     else
2422       hists[k]->Draw("sames");
2423     pad->Update();
2424     if (drawStatBox) {
2425       TPaveStats* st1 = (TPaveStats*)hists[k]->GetListOfFunctions()->FindObject("stats");
2426       if (st1 != nullptr) {
2427         dy = (entries[k] > nmin) ? 0.05 : 0.025;
2428         st1->SetLineColor(colors[k]);
2429         st1->SetTextColor(colors[k]);
2430         st1->SetY1NDC(yh - dy);
2431         st1->SetY2NDC(yh);
2432         st1->SetX1NDC(0.70);
2433         st1->SetX2NDC(0.90);
2434         yh -= dy;
2435       }
2436     }
2437     sprintf(name, "Depth %d (%s)", depth, text.c_str());
2438     legend->AddEntry(hists[k], name, "lp");
2439   }
2440   legend->Draw("same");
2441   pad->Update();
2442   if (fits < 1) {
2443     double xmin = hists[0]->GetBinLowEdge(1);
2444     int nbin = hists[0]->GetNbinsX();
2445     double xmax = hists[0]->GetBinLowEdge(nbin) + hists[0]->GetBinWidth(nbin);
2446     TLine* line = new TLine(xmin, 1.0, xmax, 1.0);
2447     line->SetLineColor(9);
2448     line->SetLineWidth(2);
2449     line->SetLineStyle(2);
2450     line->Draw("same");
2451     pad->Modified();
2452     pad->Update();
2453   }
2454   char txt1[30];
2455   double xmax = (isRealData) ? 0.33 : 0.44;
2456   TPaveText* txt2 = new TPaveText(0.11, 0.85, xmax, 0.89, "blNDC");
2457   txt2->SetFillColor(0);
2458   if (isRealData)
2459     sprintf(txt1, "CMS Preliminary");
2460   else
2461     sprintf(txt1, "CMS Simulation Preliminary");
2462   txt2->AddText(txt1);
2463   txt2->Draw("same");
2464   pad->Modified();
2465   pad->Update();
2466   if (save > 0) {
2467     sprintf(name, "%s.pdf", pad->GetName());
2468     pad->Print(name);
2469   } else if (save < 0) {
2470     sprintf(name, "%s.C", pad->GetName());
2471     pad->Print(name);
2472   }
2473 }
2474 
2475 void PlotHistCorrAsymmetry(char* infile, std::string text, std::string prefixF = "", int iformat = 0, int save = 0) {
2476   std::map<int, cfactors> cfacs;
2477   int etamin(100), etamax(-100), maxdepth(0);
2478   double scale(1.0);
2479   readCorrFactors(infile, scale, cfacs, etamin, etamax, maxdepth, iformat);
2480 
2481   gStyle->SetCanvasBorderMode(0);
2482   gStyle->SetCanvasColor(kWhite);
2483   gStyle->SetPadColor(kWhite);
2484   gStyle->SetFillColor(kWhite);
2485   gStyle->SetOptTitle(0);
2486   gStyle->SetOptStat(0);
2487   gStyle->SetOptFit(10);
2488   int colors[6] = {1, 6, 4, 7, 2, 9};
2489   int mtype[6] = {20, 21, 22, 23, 24, 33};
2490   int nbin = etamax + 1;
2491   std::vector<TH1D*> hists;
2492   std::vector<int> entries;
2493   char name[100];
2494   double dy(0);
2495   for (int j = 0; j < maxdepth; ++j) {
2496     sprintf(name, "hd%d", j + 1);
2497     TObject* ob = gROOT->FindObject(name);
2498     if (ob)
2499       ob->Delete();
2500     TH1D* h = new TH1D(name, name, nbin, 0, etamax);
2501     int nent(0);
2502     for (std::map<int, cfactors>::const_iterator itr = cfacs.begin(); itr != cfacs.end(); ++itr) {
2503       if ((itr->second).depth == j + 1) {
2504         int ieta = (itr->second).ieta;
2505         float vl1 = (itr->second).corrf;
2506         float dv1 = (itr->second).dcorr;
2507         if (ieta > 0) {
2508           for (std::map<int, cfactors>::const_iterator ktr = cfacs.begin(); ktr != cfacs.end(); ++ktr) {
2509             if (((ktr->second).depth == j + 1) && ((ktr->second).ieta == -ieta)) {
2510               float vl2 = (ktr->second).corrf;
2511               float dv2 = (ktr->second).dcorr;
2512               float val = 2.0 * (vl1 - vl2) / (vl1 + vl2);
2513               float dvl = (4.0 * sqrt(vl1 * vl1 * dv2 * dv2 + vl2 * vl2 * dv1 * dv1) / ((vl1 + vl2) * (vl1 + vl2)));
2514               int bin = ieta;
2515               h->SetBinContent(bin, val);
2516               h->SetBinError(bin, dvl);
2517               nent++;
2518             }
2519           }
2520         }
2521       }
2522     }
2523     h->SetLineColor(colors[j]);
2524     h->SetMarkerColor(colors[j]);
2525     h->SetMarkerStyle(mtype[j]);
2526     h->GetXaxis()->SetTitle("i#eta");
2527     h->GetYaxis()->SetTitle("Asymmetry in Correction Factor");
2528     h->GetYaxis()->SetLabelOffset(0.005);
2529     h->GetYaxis()->SetTitleOffset(1.20);
2530     h->GetYaxis()->SetRangeUser(-0.25, 0.25);
2531     hists.push_back(h);
2532     entries.push_back(nent);
2533     dy += 0.025;
2534   }
2535   sprintf(name, "c_%sCorrAsymmetry", prefixF.c_str());
2536   TCanvas* pad = new TCanvas(name, name, 700, 500);
2537   pad->SetRightMargin(0.10);
2538   pad->SetTopMargin(0.10);
2539   double yh = 0.90;
2540   double yl = yh - 0.035 * hists.size() - dy - 0.01;
2541   TLegend* legend = new TLegend(0.60, yl, 0.90, yl + 0.035 * hists.size());
2542   legend->SetFillColor(kWhite);
2543   for (unsigned int k = 0; k < hists.size(); ++k) {
2544     if (k == 0)
2545       hists[k]->Draw("");
2546     else
2547       hists[k]->Draw("sames");
2548     pad->Update();
2549     sprintf(name, "Depth %d (%s)", k + 1, text.c_str());
2550     legend->AddEntry(hists[k], name, "lp");
2551   }
2552   legend->Draw("same");
2553   pad->Update();
2554 
2555   TLine* line = new TLine(0.0, 0.0, etamax, 0.0);
2556   line->SetLineColor(9);
2557   line->SetLineWidth(2);
2558   line->SetLineStyle(2);
2559   line->Draw("same");
2560   pad->Update();
2561 
2562   if (save > 0) {
2563     sprintf(name, "%s.pdf", pad->GetName());
2564     pad->Print(name);
2565   } else if (save < 0) {
2566     sprintf(name, "%s.C", pad->GetName());
2567     pad->Print(name);
2568   }
2569 }
2570 
2571 void PlotHistCorrFactors(char* infile1,
2572                          std::string text1,
2573                          char* infile2,
2574                          std::string text2,
2575                          char* infile3,
2576                          std::string text3,
2577                          char* infile4,
2578                          std::string text4,
2579                          char* infile5,
2580                          std::string text5,
2581                          std::string prefixF,
2582                          bool ratio = false,
2583                          bool drawStatBox = true,
2584                          int nmin = 100,
2585                          bool isRealData = false,
2586                          int year = 2018,
2587                          int iformat = 0,
2588                          int save = 0) {
2589   std::map<int, cfactors> cfacs[5];
2590   std::vector<std::string> texts;
2591   int nfile(0), etamin(100), etamax(-100), maxdepth(0);
2592   const char* blank("");
2593   if (infile1 != blank) {
2594     readCorrFactors(infile1, 1.0, cfacs[nfile], etamin, etamax, maxdepth, iformat);
2595     if (cfacs[nfile].size() > 0) {
2596       texts.push_back(text1);
2597       ++nfile;
2598     }
2599   }
2600   if (infile2 != blank) {
2601     readCorrFactors(infile2, 1.0, cfacs[nfile], etamin, etamax, maxdepth, iformat);
2602     if (cfacs[nfile].size() > 0) {
2603       texts.push_back(text2);
2604       ++nfile;
2605     }
2606   }
2607   if (infile3 != blank) {
2608     readCorrFactors(infile3, 1.0, cfacs[nfile], etamin, etamax, maxdepth, iformat);
2609     if (cfacs[nfile].size() > 0) {
2610       texts.push_back(text3);
2611       ++nfile;
2612     }
2613   }
2614   if (infile4 != blank) {
2615     readCorrFactors(infile4, 1.0, cfacs[nfile], etamin, etamax, maxdepth, iformat);
2616     if (cfacs[nfile].size() > 0) {
2617       texts.push_back(text4);
2618       ++nfile;
2619     }
2620   }
2621   if (infile5 != blank) {
2622     readCorrFactors(infile5, 1.0, cfacs[nfile], etamin, etamax, maxdepth, iformat);
2623     if (cfacs[nfile].size() > 0) {
2624       texts.push_back(text5);
2625       ++nfile;
2626     }
2627   }
2628 
2629   if (nfile > 0) {
2630     gStyle->SetCanvasBorderMode(0);
2631     gStyle->SetCanvasColor(kWhite);
2632     gStyle->SetPadColor(kWhite);
2633     gStyle->SetFillColor(kWhite);
2634     gStyle->SetOptTitle(0);
2635     if ((!ratio) && drawStatBox) {
2636       gStyle->SetOptStat(10);
2637       gStyle->SetOptFit(10);
2638     } else {
2639       gStyle->SetOptStat(0);
2640       gStyle->SetOptFit(0);
2641     }
2642     int colors[7] = {1, 6, 4, 2, 7, 9, 46};
2643     int mtype[7] = {20, 24, 22, 23, 21, 25, 33};
2644     int nbin = etamax - etamin + 1;
2645     std::vector<TH1D*> hists;
2646     std::vector<int> entries, htype, depths;
2647     std::vector<double> fitr;
2648     char name[100];
2649     double dy(0);
2650     int fits(0);
2651     int nline(0);
2652     if (ratio) {
2653       for (int ih = 1; ih < nfile; ++ih) {
2654         for (int j = 0; j < maxdepth; ++j) {
2655           sprintf(name, "h%dd%d", ih, j + 1);
2656           TObject* ob = gROOT->FindObject(name);
2657           if (ob)
2658             ob->Delete();
2659           TH1D* h = new TH1D(name, name, nbin, etamin, etamax);
2660           double sumNum(0), sumDen(0), chi2(0);
2661           int npt(0);
2662           std::map<int, cfactors>::const_iterator ktr = cfacs[ih].begin();
2663           for (std::map<int, cfactors>::const_iterator itr = cfacs[0].begin(); itr != cfacs[0].end(); ++itr, ++ktr) {
2664             int dep = (itr->second).depth;
2665             if (dep == j + 1) {
2666               int ieta = (itr->second).ieta;
2667               int bin = ieta - etamin + 1;
2668               float val = (itr->second).corrf / (ktr->second).corrf;
2669               float dvl =
2670                   val *
2671                   sqrt((((itr->second).dcorr * (itr->second).dcorr) / ((itr->second).corrf * (itr->second).corrf)) +
2672                        (((ktr->second).dcorr * (ktr->second).dcorr) / ((ktr->second).corrf * (ktr->second).corrf)));
2673               h->SetBinContent(bin, val);
2674               h->SetBinError(bin, dvl);
2675               sumNum += (val / (dvl * dvl));
2676               sumDen += (1.0 / (dvl * dvl));
2677               ++npt;
2678               chi2 += (((val - 1.0) / dvl) * ((val - 1.0) / dvl));
2679             }
2680           }
2681           double fit = (sumDen > 0) ? (sumNum / sumDen) : 1.0;
2682           std::cout << texts[ih] << " Depth = " << (j + 1) << " Fit to Pol0: " << fit << " chi2: " << chi2 << "/" << npt
2683                     << std::endl;
2684           h->SetLineColor(colors[ih]);
2685           h->SetMarkerColor(colors[ih]);
2686           h->SetMarkerStyle(mtype[j]);
2687           h->SetMarkerSize(0.9);
2688           h->GetXaxis()->SetTitle("i#eta");
2689           if (nfile > 2)
2690             sprintf(name, "CF_{%s}/CF_{Set}", texts[0].c_str());
2691           else
2692             sprintf(name, "CF_{%s}/CF_{%s}", texts[0].c_str(), texts[ih].c_str());
2693           h->GetYaxis()->SetTitle(name);
2694           h->GetYaxis()->SetLabelOffset(0.005);
2695           h->GetYaxis()->SetTitleSize(0.036);
2696           h->GetYaxis()->SetTitleOffset(1.20);
2697           h->GetYaxis()->SetRangeUser(0.50, 1.50);
2698           hists.push_back(h);
2699           fitr.push_back(fit);
2700           htype.push_back(ih);
2701           depths.push_back(j + 1);
2702         }
2703         if ((ih == 1) || (maxdepth <= 4))
2704           nline += hists.size();
2705         else
2706           ++nline;
2707       }
2708     } else {
2709       for (int k1 = 0; k1 < nfile; ++k1) {
2710         for (int j = 0; j < maxdepth; ++j) {
2711           sprintf(name, "h%dd%d", k1, j + 1);
2712           TObject* ob = gROOT->FindObject(name);
2713           if (ob)
2714             ob->Delete();
2715           TH1D* h = new TH1D(name, name, nbin, etamin, etamax);
2716           int nent(0);
2717           for (std::map<int, cfactors>::const_iterator itr = cfacs[k1].begin(); itr != cfacs[k1].end(); ++itr) {
2718             int dep = (itr->second).depth;
2719             if (dep == j + 1) {
2720               int ieta = (itr->second).ieta;
2721               int bin = ieta - etamin + 1;
2722               float val = (itr->second).corrf;
2723               float dvl = (itr->second).dcorr;
2724               h->SetBinContent(bin, val);
2725               h->SetBinError(bin, dvl);
2726               nent++;
2727             }
2728           }
2729           if (nent > nmin) {
2730             fits++;
2731             if (drawStatBox)
2732               dy += 0.025;
2733             sprintf(name, "h%ddf%d", k1, j + 1);
2734             TObject* ob = gROOT->FindObject(name);
2735             if (ob)
2736               ob->Delete();
2737             TF1* func = new TF1(name, "pol0", etamin, etamax);
2738             h->Fit(func, "+QWLR", "");
2739           }
2740           h->SetLineColor(colors[k1]);
2741           h->SetMarkerColor(colors[k1]);
2742           h->SetMarkerStyle(mtype[j]);
2743           h->SetMarkerSize(0.9);
2744           h->GetXaxis()->SetTitle("i#eta");
2745           h->GetYaxis()->SetTitle("Correction Factor");
2746           h->GetYaxis()->SetLabelOffset(0.005);
2747           h->GetYaxis()->SetTitleOffset(1.20);
2748           h->GetYaxis()->SetRangeUser(0.5, 1.5);
2749           hists.push_back(h);
2750           entries.push_back(nent);
2751           if (drawStatBox)
2752             dy += 0.025;
2753           htype.push_back(k1);
2754           depths.push_back(j + 1);
2755         }
2756         if ((k1 <= 1) || (maxdepth <= 4))
2757           nline += hists.size();
2758         else
2759           ++nline;
2760       }
2761     }
2762     if (ratio)
2763       sprintf(name, "c_Corr%sRatio", prefixF.c_str());
2764     else
2765       sprintf(name, "c_Corr%s", prefixF.c_str());
2766     TCanvas* pad = new TCanvas(name, name, 700, 500);
2767     pad->SetRightMargin(0.10);
2768     pad->SetTopMargin(0.10);
2769     double yh = 0.90;
2770     double yl = yh - 0.035 * hists.size() - dy - 0.01;
2771     TLegend* legend = new TLegend(0.45, yl, 0.90, yl + 0.035 * nline);
2772     legend->SetFillColor(kWhite);
2773     for (unsigned int k = 0; k < hists.size(); ++k) {
2774       if (k == 0)
2775         hists[k]->Draw("");
2776       else
2777         hists[k]->Draw("sames");
2778       pad->Update();
2779       int k1 = htype[k];
2780       if (!ratio) {
2781         TPaveStats* st1 = (TPaveStats*)hists[k]->GetListOfFunctions()->FindObject("stats");
2782         if (st1 != nullptr) {
2783           dy = (entries[k] > nmin) ? 0.05 : 0.025;
2784           st1->SetLineColor(colors[k1]);
2785           st1->SetTextColor(colors[k1]);
2786           st1->SetY1NDC(yh - dy);
2787           st1->SetY2NDC(yh);
2788           st1->SetX1NDC(0.70);
2789           st1->SetX2NDC(0.90);
2790           yh -= dy;
2791         }
2792         sprintf(name, "Depth %d (%s)", depths[k], texts[k1].c_str());
2793       } else {
2794         sprintf(name, "Depth %d (Mean[CF_{%s}/CF_{%s}] = %5.3f)", depths[k], text1.c_str(), texts[k1].c_str(), fitr[k]);
2795       }
2796       if ((depths[k] == 1) || (k1 <= 1) || (maxdepth <= 4))
2797         legend->AddEntry(hists[k], name, "lp");
2798     }
2799     legend->Draw("same");
2800     TPaveText* txt0 = new TPaveText(0.11, 0.84, 0.45, 0.89, "blNDC");
2801     txt0->SetFillColor(0);
2802     char txt[40];
2803     if (isRealData)
2804       sprintf(txt, "CMS Preliminary (%d)", year);
2805     else
2806       sprintf(txt, "CMS Simulation Preliminary (%d)", year);
2807     txt0->AddText(txt);
2808     txt0->Draw("same");
2809     pad->Update();
2810     if (fits < 1) {
2811       double xmin = hists[0]->GetBinLowEdge(1);
2812       int nbin = hists[0]->GetNbinsX();
2813       double xmax = hists[0]->GetBinLowEdge(nbin) + hists[0]->GetBinWidth(nbin);
2814       TLine* line = new TLine(xmin, 1.0, xmax, 1.0);
2815       line->SetLineColor(9);
2816       line->SetLineWidth(2);
2817       line->SetLineStyle(2);
2818       line->Draw("same");
2819       pad->Update();
2820     }
2821     if (save > 0) {
2822       sprintf(name, "%s.pdf", pad->GetName());
2823       pad->Print(name);
2824     } else if (save < 0) {
2825       sprintf(name, "%s.C", pad->GetName());
2826       pad->Print(name);
2827     }
2828   }
2829 }
2830 
2831 void PlotHistCorr2Factors(char* infile1,
2832                           std::string text1,
2833                           char* infile2,
2834                           std::string text2,
2835                           int depth,
2836                           std::string prefixF,
2837                           bool ratio = true,
2838                           bool drawStatBox = false,
2839                           int nmin = 100,
2840                           bool isRealData = true,
2841                           int year = 2023,
2842                           int iformat = 0,
2843                           int save = 0) {
2844   std::map<int, cfactors> cfacs[5];
2845   std::vector<std::string> texts;
2846   int nfile(0), etamin(100), etamax(-100), maxdepth(0);
2847   const char* blank("");
2848   if (infile1 != blank) {
2849     readCorrFactors(infile1, 1.0, cfacs[nfile], etamin, etamax, maxdepth, iformat);
2850     if (cfacs[nfile].size() > 0) {
2851       texts.push_back(text1);
2852       ++nfile;
2853     }
2854   }
2855   if (infile2 != blank) {
2856     readCorrFactors(infile2, 1.0, cfacs[nfile], etamin, etamax, maxdepth, iformat);
2857     if (cfacs[nfile].size() > 0) {
2858       texts.push_back(text2);
2859       ++nfile;
2860     }
2861   }
2862 
2863   if (nfile > 0) {
2864     gStyle->SetCanvasBorderMode(0);
2865     gStyle->SetCanvasColor(kWhite);
2866     gStyle->SetPadColor(kWhite);
2867     gStyle->SetFillColor(kWhite);
2868     gStyle->SetOptTitle(0);
2869     if ((!ratio) && drawStatBox) {
2870       gStyle->SetOptStat(10);
2871       gStyle->SetOptFit(10);
2872     } else {
2873       gStyle->SetOptStat(0);
2874       gStyle->SetOptFit(0);
2875     }
2876     int colors[7] = {1, 6, 4, 2, 7, 9, 46};
2877     int mtype[7] = {20, 24, 22, 23, 21, 25, 33};
2878     int nbin = etamax - etamin + 1;
2879     std::vector<TH1D*> hists;
2880     std::vector<int> entries, htype;
2881     std::vector<double> fitr;
2882     char name[100];
2883     double dy(0);
2884     int fits(0);
2885     int nline(0);
2886     if (ratio) {
2887       for (int ih = 1; ih < nfile; ++ih) {
2888         sprintf(name, "h%dd%d", ih, depth);
2889         TObject* ob = gROOT->FindObject(name);
2890         if (ob)
2891           ob->Delete();
2892         TH1D* h = new TH1D(name, name, nbin, etamin, etamax);
2893         double sumNum(0), sumDen(0), chi2(0);
2894         int npt(0);
2895         std::map<int, cfactors>::const_iterator ktr = cfacs[ih].begin();
2896         for (std::map<int, cfactors>::const_iterator itr = cfacs[0].begin(); itr != cfacs[0].end(); ++itr, ++ktr) {
2897           int dep = (itr->second).depth;
2898           if (dep == depth) {
2899             int ieta = (itr->second).ieta;
2900             int bin = ieta - etamin + 1;
2901             float val = (itr->second).corrf / (ktr->second).corrf;
2902             float dvl =
2903                 val * sqrt((((itr->second).dcorr * (itr->second).dcorr) / ((itr->second).corrf * (itr->second).corrf)) +
2904                            (((ktr->second).dcorr * (ktr->second).dcorr) / ((ktr->second).corrf * (ktr->second).corrf)));
2905             h->SetBinContent(bin, val);
2906             h->SetBinError(bin, dvl);
2907             sumNum += (val / (dvl * dvl));
2908             sumDen += (1.0 / (dvl * dvl));
2909             ++npt;
2910             chi2 += (((val - 1.0) / dvl) * ((val - 1.0) / dvl));
2911           }
2912         }
2913         double fit = (sumDen > 0) ? (sumNum / sumDen) : 1.0;
2914         std::cout << "Depth = " << depth << " Fit to Pol0: " << fit << " chi2 = " << chi2 << "/" << npt << std::endl;
2915         h->SetLineColor(colors[ih]);
2916         h->SetMarkerColor(colors[ih]);
2917         h->SetMarkerStyle(mtype[depth - 1]);
2918         h->SetMarkerSize(0.9);
2919         h->GetXaxis()->SetTitle("i#eta");
2920         sprintf(name, "CF_{%s}/CF_{%s}", texts[0].c_str(), texts[1].c_str());
2921         h->GetYaxis()->SetTitle(name);
2922         h->GetYaxis()->SetLabelOffset(0.005);
2923         h->GetYaxis()->SetTitleSize(0.036);
2924         h->GetYaxis()->SetTitleOffset(1.20);
2925         h->GetYaxis()->SetRangeUser(0.50, 1.50);
2926         hists.push_back(h);
2927         fitr.push_back(fit);
2928         htype.push_back(ih);
2929         ++nline;
2930       }
2931     } else {
2932       for (int k1 = 0; k1 < nfile; ++k1) {
2933         sprintf(name, "h%dd%d", k1, depth);
2934         TObject* ob = gROOT->FindObject(name);
2935         if (ob)
2936           ob->Delete();
2937         TH1D* h = new TH1D(name, name, nbin, etamin, etamax);
2938         int nent(0);
2939         for (std::map<int, cfactors>::const_iterator itr = cfacs[k1].begin(); itr != cfacs[k1].end(); ++itr) {
2940           int dep = (itr->second).depth;
2941           if (dep == depth) {
2942             int ieta = (itr->second).ieta;
2943             int bin = ieta - etamin + 1;
2944             float val = (itr->second).corrf;
2945             float dvl = (itr->second).dcorr;
2946             h->SetBinContent(bin, val);
2947             h->SetBinError(bin, dvl);
2948             nent++;
2949           }
2950         }
2951         if (nent > nmin) {
2952           fits++;
2953           if (drawStatBox)
2954             dy += 0.025;
2955           sprintf(name, "h%ddf%d", k1, depth);
2956           TObject* ob = gROOT->FindObject(name);
2957           if (ob)
2958             ob->Delete();
2959           TF1* func = new TF1(name, "pol0", etamin, etamax);
2960           h->Fit(func, "+QWLR", "");
2961         }
2962         h->SetLineColor(colors[k1]);
2963         h->SetMarkerColor(colors[k1]);
2964         h->SetMarkerStyle(mtype[depth - 1]);
2965         h->SetMarkerSize(0.9);
2966         h->GetXaxis()->SetTitle("i#eta");
2967         h->GetYaxis()->SetTitle("Correction Factor");
2968         h->GetYaxis()->SetLabelOffset(0.005);
2969         h->GetYaxis()->SetTitleOffset(1.20);
2970         h->GetYaxis()->SetRangeUser(0.5, 1.5);
2971         hists.push_back(h);
2972         entries.push_back(nent);
2973         if (drawStatBox)
2974           dy += 0.025;
2975         htype.push_back(k1);
2976       }
2977       ++nline;
2978     }
2979     if (ratio)
2980       sprintf(name, "c_Corr%sD%dRatio", prefixF.c_str(), depth);
2981     else
2982       sprintf(name, "c_Corr%sD%d", prefixF.c_str(), depth);
2983     TCanvas* pad = new TCanvas(name, name, 700, 500);
2984     pad->SetRightMargin(0.10);
2985     pad->SetTopMargin(0.10);
2986     double yh = 0.90;
2987     double yl = yh - 0.035 * hists.size() - dy - 0.01;
2988     TLegend* legend = new TLegend(0.45, yl, 0.90, yl + 0.035 * nline);
2989     legend->SetFillColor(kWhite);
2990     for (unsigned int k = 0; k < hists.size(); ++k) {
2991       if (k == 0)
2992         hists[k]->Draw("");
2993       else
2994         hists[k]->Draw("sames");
2995       pad->Update();
2996       int k1 = htype[k];
2997       if (!ratio) {
2998         TPaveStats* st1 = (TPaveStats*)hists[k]->GetListOfFunctions()->FindObject("stats");
2999         if (st1 != nullptr) {
3000           dy = (entries[k] > nmin) ? 0.05 : 0.025;
3001           st1->SetLineColor(colors[k1]);
3002           st1->SetTextColor(colors[k1]);
3003           st1->SetY1NDC(yh - dy);
3004           st1->SetY2NDC(yh);
3005           st1->SetX1NDC(0.70);
3006           st1->SetX2NDC(0.90);
3007           yh -= dy;
3008         }
3009         sprintf(name, "Depth %d (%s)", depth, texts[k1].c_str());
3010       } else {
3011         sprintf(name, "Depth %d (Mean[CF_{%s}/CF_{%s}] = %5.3f)", depth, text1.c_str(), texts[k1].c_str(), fitr[k]);
3012       }
3013       legend->AddEntry(hists[k], name, "lp");
3014     }
3015     legend->Draw("same");
3016     TPaveText* txt0 = new TPaveText(0.11, 0.84, 0.45, 0.89, "blNDC");
3017     txt0->SetFillColor(0);
3018     char txt[40];
3019     if (isRealData)
3020       sprintf(txt, "CMS Preliminary (%d)", year);
3021     else
3022       sprintf(txt, "CMS Simulation Preliminary (%d)", year);
3023     txt0->AddText(txt);
3024     txt0->Draw("same");
3025     pad->Update();
3026     if (fits < 1) {
3027       double xmin = hists[0]->GetBinLowEdge(1);
3028       int nbin = hists[0]->GetNbinsX();
3029       double xmax = hists[0]->GetBinLowEdge(nbin) + hists[0]->GetBinWidth(nbin);
3030       TLine* line = new TLine(xmin, 1.0, xmax, 1.0);
3031       line->SetLineColor(9);
3032       line->SetLineWidth(2);
3033       line->SetLineStyle(2);
3034       line->Draw("same");
3035       pad->Update();
3036     }
3037     if (save > 0) {
3038       sprintf(name, "%s.pdf", pad->GetName());
3039       pad->Print(name);
3040     } else if (save < 0) {
3041       sprintf(name, "%s.C", pad->GetName());
3042       pad->Print(name);
3043     }
3044   }
3045 }
3046 
3047 void PlotHistCorrDFactors(char* infile1,
3048                           std::string text1,
3049                           char* infile2,
3050                           std::string text2,
3051                           char* infile3,
3052                           std::string text3,
3053                           char* infile4,
3054                           std::string text4,
3055                           char* infile5,
3056                           std::string text5,
3057                           int depth,
3058                           std::string prefixF,
3059                           bool ratio = true,
3060                           bool drawStatBox = false,
3061                           int nmin = 100,
3062                           bool isRealData = true,
3063                           int year = 2024,
3064                           int iformat = 0,
3065                           int save = 0) {
3066   std::map<int, cfactors> cfacs[5];
3067   std::vector<std::string> texts;
3068   int nfile(0), etamin(100), etamax(-100), maxdepth(0);
3069   const char* blank("");
3070   if (infile1 != blank) {
3071     readCorrFactors(infile1, 1.0, cfacs[nfile], etamin, etamax, maxdepth, iformat);
3072     if (cfacs[nfile].size() > 0) {
3073       texts.push_back(text1);
3074       ++nfile;
3075     }
3076   }
3077   if (infile2 != blank) {
3078     readCorrFactors(infile2, 1.0, cfacs[nfile], etamin, etamax, maxdepth, iformat);
3079     if (cfacs[nfile].size() > 0) {
3080       texts.push_back(text2);
3081       ++nfile;
3082     }
3083   }
3084   if (infile3 != blank) {
3085     readCorrFactors(infile3, 1.0, cfacs[nfile], etamin, etamax, maxdepth, iformat);
3086     if (cfacs[nfile].size() > 0) {
3087       texts.push_back(text3);
3088       ++nfile;
3089     }
3090   }
3091   if (infile4 != blank) {
3092     readCorrFactors(infile4, 1.0, cfacs[nfile], etamin, etamax, maxdepth, iformat);
3093     if (cfacs[nfile].size() > 0) {
3094       texts.push_back(text4);
3095       ++nfile;
3096     }
3097   }
3098   if (infile5 != blank) {
3099     readCorrFactors(infile5, 1.0, cfacs[nfile], etamin, etamax, maxdepth, iformat);
3100     if (cfacs[nfile].size() > 0) {
3101       texts.push_back(text5);
3102       ++nfile;
3103     }
3104   }
3105 
3106   if (nfile > 0) {
3107     gStyle->SetCanvasBorderMode(0);
3108     gStyle->SetCanvasColor(kWhite);
3109     gStyle->SetPadColor(kWhite);
3110     gStyle->SetFillColor(kWhite);
3111     gStyle->SetOptTitle(0);
3112     if ((!ratio) && drawStatBox) {
3113       gStyle->SetOptStat(10);
3114       gStyle->SetOptFit(10);
3115     } else {
3116       gStyle->SetOptStat(0);
3117       gStyle->SetOptFit(0);
3118     }
3119     int colors[7] = {1, 6, 4, 2, 7, 9, 46};
3120     int mtype[7] = {20, 24, 22, 23, 21, 25, 33};
3121     int nbin = etamax - etamin + 1;
3122     std::vector<TH1D*> hists;
3123     std::vector<int> entries, htype;
3124     std::vector<double> fitr;
3125     char name[100];
3126     double dy(0);
3127     int fits(0);
3128     int nline(0);
3129     if (ratio) {
3130       for (int ih = 1; ih < nfile; ++ih) {
3131         sprintf(name, "h%dd%d", ih, depth);
3132         TObject* ob = gROOT->FindObject(name);
3133         if (ob)
3134           ob->Delete();
3135         TH1D* h = new TH1D(name, name, nbin, etamin, etamax);
3136         double sumNum(0), sumDen(0), chi2(0);
3137         int npt(0);
3138         std::map<int, cfactors>::const_iterator ktr = cfacs[ih].begin();
3139         for (std::map<int, cfactors>::const_iterator itr = cfacs[0].begin(); itr != cfacs[0].end(); ++itr, ++ktr) {
3140           int dep = (itr->second).depth;
3141           if (dep == depth) {
3142             int ieta = (itr->second).ieta;
3143             int bin = ieta - etamin + 1;
3144             float val = (itr->second).corrf / (ktr->second).corrf;
3145             float dvl =
3146                 val * sqrt((((itr->second).dcorr * (itr->second).dcorr) / ((itr->second).corrf * (itr->second).corrf)) +
3147                            (((ktr->second).dcorr * (ktr->second).dcorr) / ((ktr->second).corrf * (ktr->second).corrf)));
3148             h->SetBinContent(bin, val);
3149             h->SetBinError(bin, dvl);
3150             sumNum += (val / (dvl * dvl));
3151             sumDen += (1.0 / (dvl * dvl));
3152             ++npt;
3153             chi2 += (((val - 1.0) / dvl) * ((val - 1.0) / dvl));
3154           }
3155         }
3156         double fit = (sumDen > 0) ? (sumNum / sumDen) : 1.0;
3157         std::cout << texts[ih] << " Depth = " << depth << " Fit to Pol0: " << fit << " chi2 = " << chi2 << "/" << npt
3158                   << std::endl;
3159         h->SetLineColor(colors[ih]);
3160         h->SetMarkerColor(colors[ih]);
3161         h->SetMarkerStyle(mtype[depth - 1]);
3162         h->SetMarkerSize(0.9);
3163         h->GetXaxis()->SetTitle("i#eta");
3164         if (nfile > 2)
3165           sprintf(name, "CF_{%s}/CF_{Set}", texts[0].c_str());
3166         else
3167           sprintf(name, "CF_{%s}/CF_{%s}", texts[0].c_str(), texts[ih].c_str());
3168         h->GetYaxis()->SetTitle(name);
3169         h->GetYaxis()->SetLabelOffset(0.005);
3170         h->GetYaxis()->SetTitleSize(0.036);
3171         h->GetYaxis()->SetTitleOffset(1.20);
3172         h->GetYaxis()->SetRangeUser(0.50, 1.50);
3173         hists.push_back(h);
3174         fitr.push_back(fit);
3175         htype.push_back(ih);
3176         ++nline;
3177       }
3178     } else {
3179       for (int k1 = 0; k1 < nfile; ++k1) {
3180         sprintf(name, "h%dd%d", k1, depth);
3181         TObject* ob = gROOT->FindObject(name);
3182         if (ob)
3183           ob->Delete();
3184         TH1D* h = new TH1D(name, name, nbin, etamin, etamax);
3185         int nent(0);
3186         for (std::map<int, cfactors>::const_iterator itr = cfacs[k1].begin(); itr != cfacs[k1].end(); ++itr) {
3187           int dep = (itr->second).depth;
3188           if (dep == depth) {
3189             int ieta = (itr->second).ieta;
3190             int bin = ieta - etamin + 1;
3191             float val = (itr->second).corrf;
3192             float dvl = (itr->second).dcorr;
3193             h->SetBinContent(bin, val);
3194             h->SetBinError(bin, dvl);
3195             nent++;
3196           }
3197         }
3198         if (nent > nmin) {
3199           fits++;
3200           if (drawStatBox)
3201             dy += 0.025;
3202           sprintf(name, "h%ddf%d", k1, depth);
3203           TObject* ob = gROOT->FindObject(name);
3204           if (ob)
3205             ob->Delete();
3206           TF1* func = new TF1(name, "pol0", etamin, etamax);
3207           h->Fit(func, "+QWLR", "");
3208         }
3209         h->SetLineColor(colors[k1]);
3210         h->SetMarkerColor(colors[k1]);
3211         h->SetMarkerStyle(mtype[depth - 1]);
3212         h->SetMarkerSize(0.9);
3213         h->GetXaxis()->SetTitle("i#eta");
3214         h->GetYaxis()->SetTitle("Correction Factor");
3215         h->GetYaxis()->SetLabelOffset(0.005);
3216         h->GetYaxis()->SetTitleOffset(1.20);
3217         h->GetYaxis()->SetRangeUser(0.5, 1.5);
3218         hists.push_back(h);
3219         entries.push_back(nent);
3220         if (drawStatBox)
3221           dy += 0.025;
3222         htype.push_back(k1);
3223         ++nline;
3224       }
3225     }
3226     if (ratio)
3227       sprintf(name, "c_Corr%sRatioD%d", prefixF.c_str(), depth);
3228     else
3229       sprintf(name, "c_Corr%sD%d", prefixF.c_str(), depth);
3230     TCanvas* pad = new TCanvas(name, name, 700, 500);
3231     pad->SetRightMargin(0.10);
3232     pad->SetTopMargin(0.10);
3233     double yh = 0.90;
3234     double yl = yh - 0.035 * hists.size() - dy - 0.01;
3235     TLegend* legend = new TLegend(0.45, yl, 0.90, yl + 0.035 * nline);
3236     legend->SetFillColor(kWhite);
3237     for (unsigned int k = 0; k < hists.size(); ++k) {
3238       if (k == 0)
3239         hists[k]->Draw("");
3240       else
3241         hists[k]->Draw("sames");
3242       pad->Update();
3243       int k1 = htype[k];
3244       if (!ratio) {
3245         TPaveStats* st1 = (TPaveStats*)hists[k]->GetListOfFunctions()->FindObject("stats");
3246         if (st1 != nullptr) {
3247           dy = (entries[k] > nmin) ? 0.05 : 0.025;
3248           st1->SetLineColor(colors[k1]);
3249           st1->SetTextColor(colors[k1]);
3250           st1->SetY1NDC(yh - dy);
3251           st1->SetY2NDC(yh);
3252           st1->SetX1NDC(0.70);
3253           st1->SetX2NDC(0.90);
3254           yh -= dy;
3255         }
3256         sprintf(name, "Depth %d (%s)", depth, texts[k1].c_str());
3257       } else {
3258         sprintf(name, "Depth %d (Mean[CF_{%s}/CF_{%s}] = %5.3f)", depth, text1.c_str(), texts[k1].c_str(), fitr[k]);
3259       }
3260       legend->AddEntry(hists[k], name, "lp");
3261     }
3262     legend->Draw("same");
3263     TPaveText* txt0 = new TPaveText(0.11, 0.84, 0.45, 0.89, "blNDC");
3264     txt0->SetFillColor(0);
3265     char txt[40];
3266     if (isRealData)
3267       sprintf(txt, "CMS Preliminary (%d)", year);
3268     else
3269       sprintf(txt, "CMS Simulation Preliminary (%d)", year);
3270     txt0->AddText(txt);
3271     txt0->Draw("same");
3272     pad->Update();
3273     if (fits < 1) {
3274       double xmin = hists[0]->GetBinLowEdge(1);
3275       int nbin = hists[0]->GetNbinsX();
3276       double xmax = hists[0]->GetBinLowEdge(nbin) + hists[0]->GetBinWidth(nbin);
3277       TLine* line = new TLine(xmin, 1.0, xmax, 1.0);
3278       line->SetLineColor(9);
3279       line->SetLineWidth(2);
3280       line->SetLineStyle(2);
3281       line->Draw("same");
3282       pad->Update();
3283     }
3284     if (save > 0) {
3285       sprintf(name, "%s.pdf", pad->GetName());
3286       pad->Print(name);
3287     } else if (save < 0) {
3288       sprintf(name, "%s.C", pad->GetName());
3289       pad->Print(name);
3290     }
3291   }
3292 }
3293 
3294 void PlotHistCorrSys(std::string infilec, int conds, std::string text, int save = 0) {
3295   char fname[100];
3296   int iformat(0);
3297   sprintf(fname, "%s_cond0.txt", infilec.c_str());
3298   int etamin(100), etamax(-100), maxdepth(0);
3299   std::map<int, cfactors> cfacs;
3300   readCorrFactors(fname, 1.0, cfacs, etamin, etamax, maxdepth, iformat);
3301   // There are good records from the master file
3302   if (cfacs.size() > 0) {
3303     // Now read the other files
3304     std::map<int, cfactors> errfacs;
3305     for (int i = 0; i < conds; ++i) {
3306       sprintf(fname, "%s_cond%d.txt", infilec.c_str(), i + 1);
3307       std::map<int, cfactors> cfacx;
3308       int etamin1(100), etamax1(-100), maxdepth1(0);
3309       readCorrFactors(fname, 1.0, cfacx, etamin1, etamax1, maxdepth1, iformat);
3310       for (std::map<int, cfactors>::const_iterator itr1 = cfacx.begin(); itr1 != cfacx.end(); ++itr1) {
3311         std::map<int, cfactors>::iterator itr2 = errfacs.find(itr1->first);
3312         float corrf = (itr1->second).corrf;
3313         if (itr2 == errfacs.end()) {
3314           errfacs[itr1->first] = cfactors(1, 0, corrf, corrf * corrf);
3315         } else {
3316           int nent = (itr2->second).ieta + 1;
3317           float c1 = (itr2->second).corrf + corrf;
3318           float c2 = (itr2->second).dcorr + (corrf * corrf);
3319           errfacs[itr1->first] = cfactors(nent, 0, c1, c2);
3320         }
3321       }
3322     }
3323     // find the RMS from the distributions
3324     for (std::map<int, cfactors>::iterator itr = errfacs.begin(); itr != errfacs.end(); ++itr) {
3325       int nent = (itr->second).ieta;
3326       float mean = (itr->second).corrf / nent;
3327       float rms2 = (itr->second).dcorr / nent - (mean * mean);
3328       float rms = rms2 > 0 ? sqrt(rms2) : 0;
3329       errfacs[itr->first] = cfactors(nent, 0, mean, rms);
3330     }
3331     // Now combine the errors and plot
3332     int k(0);
3333     gStyle->SetCanvasBorderMode(0);
3334     gStyle->SetCanvasColor(kWhite);
3335     gStyle->SetPadColor(kWhite);
3336     gStyle->SetFillColor(kWhite);
3337     gStyle->SetOptTitle(0);
3338     gStyle->SetOptStat(10);
3339     gStyle->SetOptFit(10);
3340     int colors[6] = {1, 6, 4, 7, 2, 9};
3341     int mtype[6] = {20, 21, 22, 23, 24, 33};
3342     std::vector<TH1D*> hists;
3343     char name[100];
3344     int nbin = etamax - etamin + 1;
3345     for (int j = 0; j < maxdepth; ++j) {
3346       sprintf(name, "hd%d", j + 1);
3347       TObject* ob = gROOT->FindObject(name);
3348       if (ob)
3349         ob->Delete();
3350       TH1D* h = new TH1D(name, name, nbin, etamin, etamax);
3351       h->SetLineColor(colors[j]);
3352       h->SetMarkerColor(colors[j]);
3353       h->SetMarkerStyle(mtype[j]);
3354       h->GetXaxis()->SetTitle("i#eta");
3355       h->GetYaxis()->SetTitle("Correction Factor");
3356       h->GetYaxis()->SetLabelOffset(0.005);
3357       h->GetYaxis()->SetTitleOffset(1.20);
3358       h->GetYaxis()->SetRangeUser(0.0, 2.0);
3359       hists.push_back(h);
3360     }
3361     for (std::map<int, cfactors>::iterator itr = cfacs.begin(); itr != cfacs.end(); ++itr, ++k) {
3362       std::map<int, cfactors>::iterator itr2 = errfacs.find(itr->first);
3363       float mean2 = (itr2 == errfacs.end()) ? 0 : (itr2->second).corrf;
3364       float ersys = (itr2 == errfacs.end()) ? 0 : (itr2->second).dcorr;
3365       float erstt = (itr->second).dcorr;
3366       float ertot = sqrt(erstt * erstt + ersys * ersys);
3367       float mean = (itr->second).corrf;
3368       int ieta = (itr->second).ieta;
3369       int depth = (itr->second).depth;
3370       std::cout << "[" << k << "] " << ieta << " " << depth << " " << mean << ":" << mean2 << " " << erstt << ":"
3371                 << ersys << ":" << ertot << std::endl;
3372       int bin = ieta - etamin + 1;
3373       hists[depth - 1]->SetBinContent(bin, mean);
3374       hists[depth - 1]->SetBinError(bin, ertot);
3375     }
3376     TCanvas* pad = new TCanvas("CorrFactorSys", "CorrFactorSys", 700, 500);
3377     pad->SetRightMargin(0.10);
3378     pad->SetTopMargin(0.10);
3379     double yh = 0.90;
3380     double yl = yh - 0.050 * hists.size() - 0.01;
3381     TLegend* legend = new TLegend(0.60, yl, 0.90, yl + 0.025 * hists.size());
3382     legend->SetFillColor(kWhite);
3383     for (unsigned int k = 0; k < hists.size(); ++k) {
3384       if (k == 0)
3385         hists[k]->Draw("");
3386       else
3387         hists[k]->Draw("sames");
3388       pad->Update();
3389       TPaveStats* st1 = (TPaveStats*)hists[k]->GetListOfFunctions()->FindObject("stats");
3390       if (st1 != nullptr) {
3391         st1->SetLineColor(colors[k]);
3392         st1->SetTextColor(colors[k]);
3393         st1->SetY1NDC(yh - 0.025);
3394         st1->SetY2NDC(yh);
3395         st1->SetX1NDC(0.70);
3396         st1->SetX2NDC(0.90);
3397         yh -= 0.025;
3398       }
3399       sprintf(name, "Depth %d (%s)", k + 1, text.c_str());
3400       legend->AddEntry(hists[k], name, "lp");
3401     }
3402     legend->Draw("same");
3403     pad->Update();
3404     if (save > 0) {
3405       sprintf(name, "%s.pdf", pad->GetName());
3406       pad->Print(name);
3407     } else if (save < 0) {
3408       sprintf(name, "%s.C", pad->GetName());
3409       pad->Print(name);
3410     }
3411   }
3412 }
3413 
3414 void PlotHistCorrLumis(std::string infilec, int conds, double lumi, int save = 0) {
3415   char fname[100];
3416   int iformat(0);
3417   sprintf(fname, "%s_0.txt", infilec.c_str());
3418   std::map<int, cfactors> cfacs;
3419   int etamin(100), etamax(-100), maxdepth(0);
3420   readCorrFactors(fname, 1.0, cfacs, etamin, etamax, maxdepth, iformat);
3421   int nbin = etamax - etamin + 1;
3422   std::cout << "Max Depth " << maxdepth << " and " << nbin << " eta bins for " << etamin << ":" << etamax << std::endl;
3423 
3424   // There are good records from the master file
3425   int colors[8] = {4, 2, 6, 7, 1, 9, 3, 5};
3426   int mtype[8] = {20, 21, 22, 23, 24, 25, 26, 27};
3427   if (cfacs.size() > 0) {
3428     // Now read the other files
3429     std::vector<TH1D*> hists;
3430     char name[100];
3431     for (int i = 0; i < conds; ++i) {
3432       int ih = (int)(hists.size());
3433       for (int j = 0; j < maxdepth; ++j) {
3434         sprintf(name, "hd%d%d", j + 1, i);
3435         TObject* ob = gROOT->FindObject(name);
3436         if (ob)
3437           ob->Delete();
3438         TH1D* h = new TH1D(name, name, nbin, etamin, etamax);
3439         h->SetLineColor(colors[j]);
3440         h->SetMarkerColor(colors[j]);
3441         h->SetMarkerStyle(mtype[i]);
3442         h->SetMarkerSize(0.9);
3443         h->GetXaxis()->SetTitle("i#eta");
3444         h->GetYaxis()->SetTitle("Fractional Error");
3445         h->GetYaxis()->SetLabelOffset(0.005);
3446         h->GetYaxis()->SetTitleOffset(1.20);
3447         h->GetYaxis()->SetRangeUser(0.0, 0.10);
3448         hists.push_back(h);
3449       }
3450       sprintf(fname, "%s_%d.txt", infilec.c_str(), i);
3451       int etamin1(100), etamax1(-100), maxdepth1(0);
3452       readCorrFactors(fname, 1.0, cfacs, etamin1, etamax1, maxdepth1, iformat);
3453       for (std::map<int, cfactors>::const_iterator itr = cfacs.begin(); itr != cfacs.end(); ++itr) {
3454         double value = (itr->second).dcorr / (itr->second).corrf;
3455         int bin = (itr->second).ieta - etamin + 1;
3456         hists[ih + (itr->second).depth - 1]->SetBinContent(bin, value);
3457         hists[ih + (itr->second).depth - 1]->SetBinError(bin, 0.0001);
3458       }
3459     }
3460 
3461     // Now plot the histograms
3462     gStyle->SetCanvasBorderMode(0);
3463     gStyle->SetCanvasColor(kWhite);
3464     gStyle->SetPadColor(kWhite);
3465     gStyle->SetFillColor(kWhite);
3466     gStyle->SetOptTitle(0);
3467     gStyle->SetOptStat(0);
3468     gStyle->SetOptFit(0);
3469     TCanvas* pad = new TCanvas("CorrFactorErr", "CorrFactorErr", 700, 500);
3470     pad->SetRightMargin(0.10);
3471     pad->SetTopMargin(0.10);
3472     double yh(0.89);
3473     TLegend* legend = new TLegend(0.60, yh - 0.04 * conds, 0.89, yh);
3474     legend->SetFillColor(kWhite);
3475     double lumic(lumi);
3476     for (unsigned int k = 0; k < hists.size(); ++k) {
3477       if (k == 0)
3478         hists[k]->Draw("");
3479       else
3480         hists[k]->Draw("sames");
3481       pad->Update();
3482       if (k % maxdepth == 0) {
3483         sprintf(name, "L = %5.2f fb^{-1}", lumic);
3484         legend->AddEntry(hists[k], name, "lp");
3485         lumic *= 0.5;
3486       }
3487     }
3488     legend->Draw("same");
3489     pad->Update();
3490     if (save > 0) {
3491       sprintf(name, "%s.pdf", pad->GetName());
3492       pad->Print(name);
3493     } else if (save < 0) {
3494       sprintf(name, "%s.C", pad->GetName());
3495       pad->Print(name);
3496     }
3497   }
3498 }
3499 
3500 void PlotHistCorrRel(char* infile1,
3501                      char* infile2,
3502                      std::string text1,
3503                      std::string text2,
3504                      int iformat1 = 0,
3505                      int iformat2 = 0,
3506                      int save = 0) {
3507   std::map<int, cfactors> cfacs1, cfacs2;
3508   int etamin(100), etamax(-100), maxdepth(0);
3509   readCorrFactors(infile1, 1.0, cfacs1, etamin, etamax, maxdepth, iformat1);
3510   readCorrFactors(infile2, 1.0, cfacs2, etamin, etamax, maxdepth, iformat2);
3511   std::map<int, std::pair<cfactors, cfactors> > cfacs;
3512   for (std::map<int, cfactors>::iterator itr = cfacs1.begin(); itr != cfacs1.end(); ++itr) {
3513     std::map<int, cfactors>::iterator ktr = cfacs2.find(itr->first);
3514     if (ktr == cfacs2.end()) {
3515       cfactors fac2(((itr->second).ieta), ((itr->second).depth), 0, -1);
3516       cfacs[itr->first] = std::pair<cfactors, cfactors>((itr->second), fac2);
3517     } else {
3518       cfactors fac2(ktr->second);
3519       cfacs[itr->first] = std::pair<cfactors, cfactors>((itr->second), fac2);
3520     }
3521   }
3522   for (std::map<int, cfactors>::iterator itr = cfacs2.begin(); itr != cfacs2.end(); ++itr) {
3523     std::map<int, cfactors>::const_iterator ktr = cfacs1.find(itr->first);
3524     if (ktr == cfacs1.end()) {
3525       cfactors fac1(((itr->second).ieta), ((itr->second).depth), 0, -1);
3526       cfacs[itr->first] = std::pair<cfactors, cfactors>(fac1, (itr->second));
3527     }
3528   }
3529 
3530   // There are good records in bothe the files
3531   if ((cfacs1.size() > 0) && (cfacs2.size() > 0)) {
3532     int k(0);
3533     gStyle->SetCanvasBorderMode(0);
3534     gStyle->SetCanvasColor(kWhite);
3535     gStyle->SetPadColor(kWhite);
3536     gStyle->SetFillColor(kWhite);
3537     gStyle->SetOptTitle(0);
3538     gStyle->SetOptStat(10);
3539     gStyle->SetOptFit(10);
3540     int colors[6] = {1, 6, 4, 7, 2, 9};
3541     int mtype[6] = {20, 21, 22, 23, 24, 33};
3542     std::vector<TH1D*> hists;
3543     char name[100];
3544     int nbin = etamax - etamin + 1;
3545     for (int i = 0; i < 2; ++i) {
3546       for (int j = 0; j < maxdepth; ++j) {
3547         int j1 = (i == 0) ? j : maxdepth + j;
3548         sprintf(name, "hd%d%d", i, j + 1);
3549         TObject* ob = gROOT->FindObject(name);
3550         if (ob)
3551           ob->Delete();
3552         TH1D* h = new TH1D(name, name, nbin, etamin, etamax);
3553         h->SetLineColor(colors[j1]);
3554         h->SetMarkerColor(colors[j1]);
3555         h->SetMarkerStyle(mtype[i]);
3556         h->GetXaxis()->SetTitle("i#eta");
3557         h->GetYaxis()->SetTitle("Correction Factor");
3558         h->GetYaxis()->SetLabelOffset(0.005);
3559         h->GetYaxis()->SetTitleOffset(1.20);
3560         h->GetYaxis()->SetRangeUser(0.0, 2.0);
3561         hists.push_back(h);
3562       }
3563     }
3564     for (std::map<int, std::pair<cfactors, cfactors> >::iterator it = cfacs.begin(); it != cfacs.end(); ++it, ++k) {
3565       float mean1 = (it->second).first.corrf;
3566       float error1 = (it->second).first.dcorr;
3567       float mean2 = (it->second).second.corrf;
3568       float error2 = (it->second).second.dcorr;
3569       int ieta = (it->second).first.ieta;
3570       int depth = (it->second).first.depth;
3571       /*
3572       std::cout << "[" << k << "] " << ieta << " " << depth << " " 
3573         << mean1 << ":" << mean2 << " " << error1 << ":" << error2
3574         << std::endl;
3575       */
3576       int bin = ieta - etamin + 1;
3577       if (error1 >= 0) {
3578         hists[depth - 1]->SetBinContent(bin, mean1);
3579         hists[depth - 1]->SetBinError(bin, error1);
3580       }
3581       if (error2 >= 0) {
3582         hists[maxdepth + depth - 1]->SetBinContent(bin, mean2);
3583         hists[maxdepth + depth - 1]->SetBinError(bin, error2);
3584       }
3585     }
3586     TCanvas* pad = new TCanvas("CorrFactors", "CorrFactors", 700, 500);
3587     pad->SetRightMargin(0.10);
3588     pad->SetTopMargin(0.10);
3589     double yh = 0.90;
3590     double yl = yh - 0.050 * hists.size() - 0.01;
3591     TLegend* legend = new TLegend(0.60, yl, 0.90, yl + 0.025 * hists.size());
3592     legend->SetFillColor(kWhite);
3593     for (unsigned int k = 0; k < hists.size(); ++k) {
3594       if (k == 0)
3595         hists[k]->Draw("");
3596       else
3597         hists[k]->Draw("sames");
3598       pad->Update();
3599       TPaveStats* st1 = (TPaveStats*)hists[k]->GetListOfFunctions()->FindObject("stats");
3600       if (st1 != nullptr) {
3601         st1->SetLineColor(colors[k]);
3602         st1->SetTextColor(colors[k]);
3603         st1->SetY1NDC(yh - 0.025);
3604         st1->SetY2NDC(yh);
3605         st1->SetX1NDC(0.70);
3606         st1->SetX2NDC(0.90);
3607         yh -= 0.025;
3608       }
3609       if (k < (unsigned int)(maxdepth)) {
3610         sprintf(name, "Depth %d (%s)", k + 1, text1.c_str());
3611       } else {
3612         sprintf(name, "Depth %d (%s)", k - maxdepth + 1, text2.c_str());
3613       }
3614       legend->AddEntry(hists[k], name, "lp");
3615     }
3616     legend->Draw("same");
3617     pad->Update();
3618     if (save > 0) {
3619       sprintf(name, "%s.pdf", pad->GetName());
3620       pad->Print(name);
3621     } else if (save < 0) {
3622       sprintf(name, "%s.C", pad->GetName());
3623       pad->Print(name);
3624     }
3625   }
3626 }
3627 
3628 void PlotHistCorrDepth(char* infile1,
3629                        char* infile2,
3630                        std::string text1,
3631                        std::string text2,
3632                        int depth,
3633                        int ietamax,
3634                        int iformat1 = 0,
3635                        int iformat2 = 0,
3636                        int save = 0,
3637                        int debug = 1) {
3638   std::map<int, cfactors> cfacs1, cfacs2;
3639   int etamin(100), etamax(-100), maxdepth(0);
3640   readCorrFactors(infile1, 1.0, cfacs1, etamin, etamax, maxdepth, iformat1, (debug > 1));
3641   readCorrFactors(infile2, 1.0, cfacs2, etamin, etamax, maxdepth, iformat2, (debug > 1));
3642 
3643   double sumNum(0), sumDen(0), ratMax(0);
3644   int npt0(0), npt1(0);
3645   for (std::map<int, cfactors>::iterator itr = cfacs1.begin(); itr != cfacs1.end(); ++itr) {
3646     std::map<int, cfactors>::iterator ktr = cfacs2.find(itr->first);
3647     if ((ktr != cfacs2.end()) && ((itr->second).depth == depth)) {
3648       double er1 = (itr->second).dcorr / (itr->second).corrf;
3649       double er2 = (ktr->second).dcorr / (ktr->second).corrf;
3650       double tmp = (ktr->second).corrf / (itr->second).corrf;
3651       double dtmp = tmp * sqrt(er1 * er1 + er2 * er2);
3652       double rat = (tmp > 1.0) ? 1.0 / tmp : tmp;
3653       double drt = (tmp > 1.0) ? dtmp / (tmp * tmp) : dtmp;
3654       rat = fabs(1.0 - rat);
3655       ratMax = std::max(ratMax, rat);
3656       ++npt0;
3657       if (debug > 0)
3658         std::cout << std::hex << (itr->first) << std::dec << " ieta:depth" << (itr->second).ieta << ":"
3659                   << (itr->second).depth << " Corr " << (itr->second).corrf << ":" << (ktr->second).corrf << " Ratio "
3660                   << rat << ":" << drt << std::endl;
3661       if (std::abs((itr->second).ieta) <= ietamax) {
3662         sumNum += (rat / (drt * drt));
3663         sumDen += (1.0 / (drt * drt));
3664         ++npt1;
3665       }
3666     }
3667   }
3668   sumNum = (sumDen > 0) ? (sumNum / sumDen) : 0;
3669   sumDen = (sumDen > 0) ? 1.0 / sqrt(sumDen) : 0;
3670   std::cout << "Get Ratio of mean for " << npt0 << ":" << npt1 << " points for depth " << depth << " Mean " << sumNum
3671             << " +- " << sumDen << " (Maximum " << ratMax << ")" << std::endl;
3672 
3673   gStyle->SetCanvasBorderMode(0);
3674   gStyle->SetCanvasColor(kWhite);
3675   gStyle->SetPadColor(kWhite);
3676   gStyle->SetFillColor(kWhite);
3677   gStyle->SetOptTitle(0);
3678   gStyle->SetOptStat(0);
3679   gStyle->SetOptFit(0);
3680   int colors[2] = {1, 2};
3681   int mtype[2] = {20, 24};
3682   int nbin = etamax - etamin + 1;
3683   std::vector<TH1D*> hists;
3684   char name[100];
3685   for (int j = 0; j < 2; ++j) {
3686     sprintf(name, "hd%d", (j + 1));
3687     TObject* ob = gROOT->FindObject(name);
3688     if (ob)
3689       ob->Delete();
3690     TH1D* h = new TH1D(name, name, nbin, etamin, etamax);
3691     if (j == 0) {
3692       for (std::map<int, cfactors>::const_iterator itr = cfacs1.begin(); itr != cfacs1.end(); ++itr) {
3693         if ((itr->second).depth == depth) {
3694           int ieta = (itr->second).ieta;
3695           int bin = ieta - etamin + 1;
3696           float val = (itr->second).corrf;
3697           float dvl = (itr->second).dcorr;
3698           h->SetBinContent(bin, val);
3699           h->SetBinError(bin, dvl);
3700         }
3701       }
3702     } else {
3703       for (std::map<int, cfactors>::const_iterator itr = cfacs2.begin(); itr != cfacs2.end(); ++itr) {
3704         if ((itr->second).depth == depth) {
3705           int ieta = (itr->second).ieta;
3706           int bin = ieta - etamin + 1;
3707           float val = (itr->second).corrf;
3708           float dvl = (itr->second).dcorr;
3709           h->SetBinContent(bin, val);
3710           h->SetBinError(bin, dvl);
3711         }
3712       }
3713     }
3714     h->SetLineColor(colors[j]);
3715     h->SetMarkerColor(colors[j]);
3716     h->SetMarkerStyle(mtype[j]);
3717     h->GetXaxis()->SetTitle("i#eta");
3718     h->GetYaxis()->SetTitle("Correction Factor");
3719     h->GetYaxis()->SetLabelOffset(0.005);
3720     h->GetYaxis()->SetTitleOffset(1.20);
3721     h->GetYaxis()->SetRangeUser(0.0, 2.0);
3722     hists.push_back(h);
3723   }
3724   sprintf(name, "c_CorrFactorDepth%d", depth);
3725   TCanvas* pad = new TCanvas(name, name, 700, 500);
3726   pad->SetRightMargin(0.10);
3727   pad->SetTopMargin(0.10);
3728   double yl = 0.25;
3729   TLegend* legend = new TLegend(0.25, yl, 0.85, yl + 0.04 * hists.size());
3730   legend->SetFillColor(kWhite);
3731   for (unsigned int k = 0; k < hists.size(); ++k) {
3732     if (k == 0) {
3733       hists[k]->Draw("");
3734       sprintf(name, "Depth %d (%s)", depth, text1.c_str());
3735     } else {
3736       hists[k]->Draw("sames");
3737       sprintf(name, "Depth %d (%s)", depth, text2.c_str());
3738     }
3739     pad->Update();
3740     legend->AddEntry(hists[k], name, "lp");
3741   }
3742   legend->Draw("same");
3743   pad->Update();
3744   if (save > 0) {
3745     sprintf(name, "%s.pdf", pad->GetName());
3746     pad->Print(name);
3747   } else if (save < 0) {
3748     sprintf(name, "%s.C", pad->GetName());
3749     pad->Print(name);
3750   }
3751 }
3752 
3753 void PlotFourHists(std::string infile,
3754                    std::string prefix0,
3755                    int type = 0,
3756                    int drawStatBox = 0,
3757                    bool normalize = false,
3758                    int save = 0,
3759                    std::string prefix1 = "",
3760                    std::string text1 = "",
3761                    std::string prefix2 = "",
3762                    std::string text2 = "",
3763                    std::string prefix3 = "",
3764                    std::string text3 = "",
3765                    std::string prefix4 = "",
3766                    std::string text4 = "") {
3767   int colors[4] = {2, 4, 6, 1};
3768   std::string names[5] = {"eta03", "eta13", "eta23", "eta33", "eta43"};
3769   std::string xtitle[5] = {"i#eta", "i#eta", "i#eta", "i#eta", "i#eta"};
3770   std::string ytitle[5] = {"Tracks", "Tracks", "Tracks", "Tracks", "Tracks"};
3771   std::string title[5] = {"All Tracks (p = 40:60 GeV)",
3772                           "Good Quality Tracks (p = 40:60 GeV)",
3773                           "Selected Tracks (p = 40:60 GeV)",
3774                           "Isolated Tracks (p = 40:60 GeV)",
3775                           "Isolated MIP Tracks (p = 40:60 GeV)"};
3776 
3777   gStyle->SetCanvasBorderMode(0);
3778   gStyle->SetCanvasColor(kWhite);
3779   gStyle->SetPadColor(kWhite);
3780   gStyle->SetFillColor(kWhite);
3781   gStyle->SetOptTitle(0);
3782   gStyle->SetOptFit(0);
3783   if (drawStatBox == 0)
3784     gStyle->SetOptStat(0);
3785   else
3786     gStyle->SetOptStat(1110);
3787 
3788   if (type < 0 || type > 4)
3789     type = 0;
3790   char name[100], namep[100];
3791   TFile* file = new TFile(infile.c_str());
3792 
3793   std::vector<TH1D*> hists;
3794   std::vector<int> kks;
3795   std::vector<std::string> texts;
3796   for (int k = 0; k < 4; ++k) {
3797     std::string prefix, text;
3798     if (k == 0) {
3799       prefix = prefix1;
3800       text = text1;
3801     } else if (k == 1) {
3802       prefix = prefix2;
3803       text = text2;
3804     } else if (k == 2) {
3805       prefix = prefix3;
3806       text = text3;
3807     } else {
3808       prefix = prefix4;
3809       text = text4;
3810     }
3811     if (prefix != "") {
3812       sprintf(name, "%s%s", prefix.c_str(), names[type].c_str());
3813       TH1D* hist1 = (TH1D*)file->FindObjectAny(name);
3814       if (hist1 != nullptr) {
3815         hists.push_back((TH1D*)(hist1->Clone()));
3816         kks.push_back(k);
3817         texts.push_back(text);
3818       }
3819     }
3820   }
3821   if (hists.size() > 0) {
3822     sprintf(namep, "c_%s%s", prefix0.c_str(), names[type].c_str());
3823     double ymax(0.90), dy(0.13);
3824     double ymx0 = (drawStatBox == 0) ? (ymax - .01) : (ymax - dy * hists.size() - .01);
3825     TCanvas* pad = new TCanvas(namep, namep, 700, 500);
3826     TLegend* legend = new TLegend(0.64, ymx0 - 0.05 * hists.size(), 0.89, ymx0);
3827     legend->SetFillColor(kWhite);
3828     pad->SetRightMargin(0.10);
3829     pad->SetTopMargin(0.10);
3830     for (unsigned int jk = 0; jk < hists.size(); ++jk) {
3831       int k = kks[jk];
3832       hists[jk]->GetXaxis()->SetTitleSize(0.040);
3833       hists[jk]->GetXaxis()->SetTitle(xtitle[type].c_str());
3834       hists[jk]->GetYaxis()->SetTitle(ytitle[type].c_str());
3835       hists[jk]->GetYaxis()->SetLabelOffset(0.005);
3836       hists[jk]->GetYaxis()->SetLabelSize(0.035);
3837       hists[jk]->GetYaxis()->SetTitleSize(0.040);
3838       hists[jk]->GetYaxis()->SetTitleOffset(1.15);
3839       hists[jk]->SetMarkerStyle(20);
3840       hists[jk]->SetMarkerColor(colors[k]);
3841       hists[jk]->SetLineColor(colors[k]);
3842       if (normalize) {
3843         if (jk == 0)
3844           hists[jk]->DrawNormalized();
3845         else
3846           hists[jk]->DrawNormalized("sames");
3847       } else {
3848         if (jk == 0)
3849           hists[jk]->Draw();
3850         else
3851           hists[jk]->Draw("sames");
3852       }
3853       pad->Update();
3854       TPaveStats* st1 = (TPaveStats*)hists[jk]->GetListOfFunctions()->FindObject("stats");
3855       if (st1 != nullptr) {
3856         double ymin = ymax - dy;
3857         st1->SetLineColor(colors[k]);
3858         st1->SetTextColor(colors[k]);
3859         st1->SetY1NDC(ymin);
3860         st1->SetY2NDC(ymax);
3861         st1->SetX1NDC(0.70);
3862         st1->SetX2NDC(0.90);
3863         ymax = ymin;
3864       }
3865       sprintf(name, "%s", texts[jk].c_str());
3866       legend->AddEntry(hists[jk], name, "lp");
3867     }
3868     legend->Draw("same");
3869     pad->Update();
3870     TPaveText* txt1 = new TPaveText(0.10, 0.905, 0.80, 0.95, "blNDC");
3871     txt1->SetFillColor(0);
3872     char txt[100];
3873     sprintf(txt, "%s", title[type].c_str());
3874     txt1->AddText(txt);
3875     txt1->Draw("same");
3876     /*
3877     TPaveText *txt2 = new TPaveText(0.11,0.825,0.33,0.895,"blNDC");
3878     txt2->SetFillColor(0);
3879     sprintf (txt, "CMS Preliminary");
3880     txt2->AddText(txt);
3881     txt2->Draw("same");
3882     */
3883     pad->Modified();
3884     pad->Update();
3885     if (save > 0) {
3886       sprintf(name, "%s.pdf", pad->GetName());
3887       pad->Print(name);
3888     } else if (save < 0) {
3889       sprintf(name, "%s.C", pad->GetName());
3890       pad->Print(name);
3891     }
3892   }
3893 }
3894 
3895 void PlotPUCorrHists(std::string infile = "corrfac.root",
3896                      std::string prefix = "",
3897                      int drawStatBox = 0,
3898                      bool approve = true,
3899                      int save = 0) {
3900   std::string name1[4] = {"W0", "W1", "W2", "P"};
3901   std::string name2[4] = {"All", "Barrel", "Endcap", ""};
3902   std::string name3[2] = {"", "p = 40:60 GeV"};
3903   std::string name4[2] = {"Loose Isolation", "Tight Isolation"};
3904   std::string xtitle[4] = {"Correction Factor", "Correction Factor", "Correction Factor", "i#eta"};
3905   std::string ytitle[4] = {"Tracks", "Tracks", "Tracks", "Correction Factor"};
3906 
3907   gStyle->SetCanvasBorderMode(0);
3908   gStyle->SetCanvasColor(kWhite);
3909   gStyle->SetPadColor(kWhite);
3910   gStyle->SetFillColor(kWhite);
3911   gStyle->SetOptTitle(0);
3912   gStyle->SetOptFit(0);
3913   if (drawStatBox == 0)
3914     gStyle->SetOptStat(0);
3915   else
3916     gStyle->SetOptStat(1110);
3917 
3918   char name[100], namep[100], title[100];
3919   TFile* file = new TFile(infile.c_str());
3920 
3921   if (file != nullptr) {
3922     for (int i1 = 0; i1 < 4; ++i1) {
3923       for (int i2 = 0; i2 < 2; ++i2) {
3924         for (int i3 = 0; i3 < 2; ++i3) {
3925           sprintf(name, "%s%d%d", name1[i1].c_str(), i2, i3);
3926           if (i2 == 0)
3927             sprintf(title, "%s Tracks Selected with %s", name2[i1].c_str(), name4[i3].c_str());
3928           else
3929             sprintf(title, "%s Tracks Selected with %s (%s)", name2[i1].c_str(), name4[i3].c_str(), name3[i2].c_str());
3930           TH1D* hist1(nullptr);
3931           TProfile* hist2(nullptr);
3932           if (i1 != 3) {
3933             TH1D* hist = (TH1D*)file->FindObjectAny(name);
3934             if (hist != nullptr) {
3935               hist1 = (TH1D*)(hist->Clone());
3936               hist1->GetXaxis()->SetTitleSize(0.040);
3937               hist1->GetXaxis()->SetTitle(xtitle[i1].c_str());
3938               hist1->GetYaxis()->SetTitle(ytitle[i1].c_str());
3939               hist1->GetYaxis()->SetLabelOffset(0.005);
3940               hist1->GetYaxis()->SetLabelSize(0.035);
3941               hist1->GetYaxis()->SetTitleSize(0.040);
3942               hist1->GetYaxis()->SetTitleOffset(1.15);
3943             }
3944           } else {
3945             TProfile* hist = (TProfile*)file->FindObjectAny(name);
3946             if (hist != nullptr) {
3947               hist2 = (TProfile*)(hist->Clone());
3948               hist2->GetXaxis()->SetTitleSize(0.040);
3949               hist2->GetXaxis()->SetTitle(xtitle[i1].c_str());
3950               hist2->GetYaxis()->SetTitle(ytitle[i1].c_str());
3951               hist2->GetYaxis()->SetLabelOffset(0.005);
3952               hist2->GetYaxis()->SetLabelSize(0.035);
3953               hist2->GetYaxis()->SetTitleSize(0.040);
3954               hist2->GetYaxis()->SetTitleOffset(1.15);
3955               //          hist2->GetYaxis()->SetRangeUser(0.0, 1.5);
3956               hist2->SetMarkerStyle(20);
3957             }
3958           }
3959           if ((hist1 != nullptr) || (hist2 != nullptr)) {
3960             sprintf(namep, "c_%s%s", name, prefix.c_str());
3961             TCanvas* pad = new TCanvas(namep, namep, 700, 500);
3962             pad->SetRightMargin(0.10);
3963             pad->SetTopMargin(0.10);
3964             if (hist1 != nullptr) {
3965               pad->SetLogy();
3966               hist1->Draw();
3967               pad->Update();
3968               TPaveStats* st1 = (TPaveStats*)hist1->GetListOfFunctions()->FindObject("stats");
3969               if (st1 != nullptr) {
3970                 st1->SetY1NDC(0.77);
3971                 st1->SetY2NDC(0.90);
3972                 st1->SetX1NDC(0.70);
3973                 st1->SetX2NDC(0.90);
3974               }
3975             } else {
3976               hist2->Draw();
3977               pad->Update();
3978             }
3979             TPaveText* txt1 = new TPaveText(0.10, 0.905, 0.80, 0.95, "blNDC");
3980             txt1->SetFillColor(0);
3981             char txt[100];
3982             sprintf(txt, "%s", title);
3983             txt1->AddText(txt);
3984             txt1->Draw("same");
3985             if (approve) {
3986               double xoff = (i1 == 3) ? 0.11 : 0.22;
3987               TPaveText* txt2 = new TPaveText(xoff, 0.825, xoff + 0.22, 0.895, "blNDC");
3988               txt2->SetFillColor(0);
3989               sprintf(txt, "CMS Preliminary");
3990               txt2->AddText(txt);
3991               txt2->Draw("same");
3992             }
3993             pad->Modified();
3994             pad->Update();
3995             if (save > 0) {
3996               sprintf(name, "%s.pdf", pad->GetName());
3997               pad->Print(name);
3998             } else if (save < 0) {
3999               sprintf(name, "%s.C", pad->GetName());
4000               pad->Print(name);
4001             }
4002           }
4003         }
4004       }
4005     }
4006   }
4007 }
4008 
4009 void PlotHistCorr(const char* infile,
4010                   std::string prefix,
4011                   std::string text0,
4012                   int eta = 0,
4013                   int mode = 1,
4014                   bool drawStatBox = true,
4015                   int save = 0) {
4016   gStyle->SetCanvasBorderMode(0);
4017   gStyle->SetCanvasColor(kWhite);
4018   gStyle->SetPadColor(kWhite);
4019   gStyle->SetFillColor(kWhite);
4020   gStyle->SetOptTitle(0);
4021   if (drawStatBox)
4022     gStyle->SetOptStat(1100);
4023   else
4024     gStyle->SetOptStat(0);
4025 
4026   std::string tags[3] = {"UnNoPU", "UnPU", "Cor"};
4027   std::string text[3] = {"Uncorrected no PU", "Uncorrected PU", "Corrected PU"};
4028   int colors[3] = {1, 4, 2};
4029   int styles[3] = {1, 3, 2};
4030   TFile* file = new TFile(infile);
4031   if (mode < 0 || mode > 2)
4032     mode = 1;
4033   int etamin = (eta == 0) ? -27 : eta;
4034   int etamax = (eta == 0) ? 27 : eta;
4035   for (int ieta = etamin; ieta <= etamax; ++ieta) {
4036     char name[20];
4037     double yh(0.90), dy(0.09);
4038     double yh1 = drawStatBox ? (yh - 3 * dy - 0.01) : (yh - 0.01);
4039     TLegend* legend = new TLegend(0.55, yh1 - 0.15, 0.89, yh1);
4040     legend->SetFillColor(kWhite);
4041     sprintf(name, "c_%sEovp%d", prefix.c_str(), ieta);
4042     TCanvas* pad = new TCanvas(name, name, 700, 500);
4043     pad->SetRightMargin(0.10);
4044     pad->SetTopMargin(0.10);
4045     TH1D* hist[3];
4046     double ymax(0);
4047     for (int k = 0; k < 3; ++k) {
4048       if (k < 2)
4049         sprintf(name, "EovP_ieta%d%s", ieta, tags[k].c_str());
4050       else
4051         sprintf(name, "EovP_ieta%dCor%dPU", ieta, mode);
4052       TH1D* hist1 = (TH1D*)file->FindObjectAny(name);
4053       if (hist1 != nullptr) {
4054         hist[k] = (TH1D*)(hist1->Clone());
4055         ymax = std::max(ymax, (hist1->GetMaximum()));
4056       }
4057     }
4058     int imax = 10 * (2 + int(0.1 * ymax));
4059     for (int k = 0; k < 3; ++k) {
4060       hist[k]->GetYaxis()->SetLabelOffset(0.005);
4061       hist[k]->GetYaxis()->SetTitleOffset(1.20);
4062       hist[k]->GetXaxis()->SetTitle("E/p");
4063       hist[k]->GetYaxis()->SetTitle("Tracks");
4064       hist[k]->SetLineColor(colors[k]);
4065       hist[k]->SetLineStyle(styles[k]);
4066       hist[k]->GetYaxis()->SetRangeUser(0.0, imax);
4067       if (k == 0)
4068         hist[k]->Draw();
4069       else
4070         hist[k]->Draw("sames");
4071       legend->AddEntry(hist[k], text[k].c_str(), "lp");
4072       pad->Update();
4073       if (drawStatBox) {
4074         TPaveStats* st1 = (TPaveStats*)hist[k]->GetListOfFunctions()->FindObject("stats");
4075         if (st1 != nullptr) {
4076           st1->SetLineColor(colors[k]);
4077           st1->SetTextColor(colors[k]);
4078           st1->SetY1NDC(yh - dy);
4079           st1->SetY2NDC(yh);
4080           st1->SetX1NDC(0.70);
4081           st1->SetX2NDC(0.90);
4082           yh -= dy;
4083         }
4084       }
4085     }
4086     pad->Update();
4087     legend->Draw("same");
4088     pad->Update();
4089     TPaveText* txt1 = new TPaveText(0.10, 0.905, 0.80, 0.95, "blNDC");
4090     txt1->SetFillColor(0);
4091     char title[100];
4092     sprintf(title, "%s for i#eta = %d", text0.c_str(), ieta);
4093     txt1->AddText(title);
4094     txt1->Draw("same");
4095     pad->Modified();
4096     pad->Update();
4097     if (save > 0) {
4098       sprintf(name, "%s.pdf", pad->GetName());
4099       pad->Print(name);
4100     } else if (save < 0) {
4101       sprintf(name, "%s.C", pad->GetName());
4102       pad->Print(name);
4103     }
4104   }
4105 }
4106 
4107 void PlotPropertyHist(const char* infile,
4108                       std::string prefix,
4109                       std::string text,
4110                       int etaMax = 25,
4111                       double lumi = 0,
4112                       double ener = 13.0,
4113                       bool isRealData = false,
4114                       bool drawStatBox = true,
4115                       int save = 0) {
4116   std::string name0[3] = {"energyE2", "energyH2", "energyP2"};
4117   std::string title0[3] = {"Energy in ECAL", "Energy in HCAL", "Track Momentum"};
4118   std::string xtitl0[3] = {"Energy (GeV)", "Energy (GeV)", "p (GeV)"};
4119   std::string name1[5] = {"eta02", "eta12", "eta22", "eta32", "eta42"};
4120   std::string name10[5] = {"eta0", "eta1", "eta2", "eta3", "eta4"};
4121   std::string xtitl1 = "i#eta";
4122   std::string name2[5] = {"p0", "p1", "p2", "p3", "p4"};
4123   std::string xtitl2 = "p (GeV)";
4124   std::string title1[5] = {"Tracks with p=40:60 GeV",
4125                            "Good Quality Tracks with p=40:60 GeV",
4126                            "Selected Tracks with p=40:60 GeV",
4127                            "Isolated Tracks with p=40:60 GeV",
4128                            "Isolated Tracks with p=40:60 GeV and MIPS in ECAL"};
4129   std::string title2[5] = {
4130       "All Tracks", "Good Quality Tracks", "Selected Tracks", "Isolated Tracks", "Isolated Tracks with MIPS in ECAL"};
4131   std::string ytitle = "Tracks";
4132 
4133   gStyle->SetCanvasBorderMode(0);
4134   gStyle->SetCanvasColor(kWhite);
4135   gStyle->SetPadColor(kWhite);
4136   gStyle->SetFillColor(kWhite);
4137   gStyle->SetOptTitle(0);
4138   if (drawStatBox)
4139     gStyle->SetOptStat(1110);
4140   else
4141     gStyle->SetOptStat(0);
4142   gStyle->SetOptFit(0);
4143 
4144   TFile* file = new TFile(infile);
4145   char name[100], namep[100];
4146   for (int k = 1; k <= etaMax; ++k) {
4147     for (int j = 0; j < 3; ++j) {
4148       sprintf(name, "%s%s%d", prefix.c_str(), name0[j].c_str(), k);
4149       TH1D* hist1 = (TH1D*)file->FindObjectAny(name);
4150       if (hist1 != nullptr) {
4151         TH1D* hist = (TH1D*)(hist1->Clone());
4152         double ymin(0.90);
4153         sprintf(namep, "c_%s", name);
4154         TCanvas* pad = new TCanvas(namep, namep, 700, 500);
4155         pad->SetLogy();
4156         pad->SetRightMargin(0.10);
4157         pad->SetTopMargin(0.10);
4158         hist->GetXaxis()->SetTitleSize(0.04);
4159         hist->GetXaxis()->SetTitle(xtitl0[j].c_str());
4160         hist->GetYaxis()->SetTitle(ytitle.c_str());
4161         hist->GetYaxis()->SetLabelOffset(0.005);
4162         hist->GetYaxis()->SetTitleSize(0.04);
4163         hist->GetYaxis()->SetLabelSize(0.035);
4164         hist->GetYaxis()->SetTitleOffset(1.10);
4165         hist->SetMarkerStyle(20);
4166         hist->Draw();
4167         pad->Update();
4168         TPaveStats* st1 = (TPaveStats*)hist->GetListOfFunctions()->FindObject("stats");
4169         if (st1 != nullptr) {
4170           st1->SetY1NDC(0.78);
4171           st1->SetY2NDC(0.90);
4172           st1->SetX1NDC(0.65);
4173           st1->SetX2NDC(0.90);
4174         }
4175         pad->Update();
4176         double ymx(0.96), xmi(0.35), xmx(0.95);
4177         char txt[100];
4178         if (lumi > 0.1) {
4179           ymx = ymin - 0.005;
4180           xmi = 0.45;
4181           TPaveText* txt0 = new TPaveText(0.65, 0.91, 0.90, 0.96, "blNDC");
4182           txt0->SetFillColor(0);
4183           sprintf(txt, "%4.1f TeV %5.1f fb^{-1}", ener, lumi);
4184           txt0->AddText(txt);
4185           txt0->Draw("same");
4186         }
4187         double ymi = ymx - 0.05;
4188         TPaveText* txt1 = new TPaveText(xmi, ymi, xmx, ymx, "blNDC");
4189         txt1->SetFillColor(0);
4190         if (text == "") {
4191           sprintf(txt, "%s", title0[j].c_str());
4192         } else {
4193           sprintf(txt, "%s (%s)", title0[j].c_str(), text.c_str());
4194         }
4195         txt1->AddText(txt);
4196         txt1->Draw("same");
4197         double xmax = (isRealData) ? 0.24 : 0.35;
4198         ymi = 0.91;
4199         ymx = ymi + 0.05;
4200         TPaveText* txt2 = new TPaveText(0.02, ymi, xmax, ymx, "blNDC");
4201         txt2->SetFillColor(0);
4202         if (isRealData)
4203           sprintf(txt, "CMS Preliminary");
4204         else
4205           sprintf(txt, "CMS Simulation Preliminary");
4206         txt2->AddText(txt);
4207         txt2->Draw("same");
4208         pad->Modified();
4209         pad->Update();
4210         if (save > 0) {
4211           sprintf(name, "%s.pdf", pad->GetName());
4212           pad->Print(name);
4213         } else if (save < 0) {
4214           sprintf(name, "%s.C", pad->GetName());
4215           pad->Print(name);
4216         }
4217       }
4218     }
4219   }
4220 
4221   for (int k = 0; k < 3; ++k) {
4222     for (int j = 0; j < 5; ++j) {
4223       if (k == 0)
4224         sprintf(name, "%s%s", prefix.c_str(), name1[j].c_str());
4225       else if (k == 1)
4226         sprintf(name, "%s%s", prefix.c_str(), name10[j].c_str());
4227       else
4228         sprintf(name, "%s%s", prefix.c_str(), name2[j].c_str());
4229       TH1D* hist1 = (TH1D*)file->FindObjectAny(name);
4230       if (hist1 != nullptr) {
4231         TH1D* hist = (TH1D*)(hist1->Clone());
4232         double ymin(0.90);
4233         sprintf(namep, "c_%s", name);
4234         TCanvas* pad = new TCanvas(namep, namep, 700, 500);
4235         pad->SetLogy();
4236         pad->SetRightMargin(0.10);
4237         pad->SetTopMargin(0.10);
4238         hist->GetXaxis()->SetTitleSize(0.04);
4239         if (k <= 1)
4240           hist->GetXaxis()->SetTitle(xtitl1.c_str());
4241         else
4242           hist->GetXaxis()->SetTitle(xtitl2.c_str());
4243         hist->GetYaxis()->SetTitle(ytitle.c_str());
4244         hist->GetYaxis()->SetLabelOffset(0.005);
4245         hist->GetYaxis()->SetTitleSize(0.04);
4246         hist->GetYaxis()->SetLabelSize(0.035);
4247         hist->GetYaxis()->SetTitleOffset(1.10);
4248         hist->SetMarkerStyle(20);
4249         hist->Draw();
4250         pad->Update();
4251         TPaveStats* st1 = (TPaveStats*)hist->GetListOfFunctions()->FindObject("stats");
4252         if (st1 != nullptr) {
4253           st1->SetY1NDC(0.78);
4254           st1->SetY2NDC(0.90);
4255           st1->SetX1NDC(0.65);
4256           st1->SetX2NDC(0.90);
4257         }
4258         pad->Update();
4259         double ymx(0.96), xmi(0.35), xmx(0.95);
4260         char txt[100];
4261         if (lumi > 0.1) {
4262           ymx = ymin - 0.005;
4263           xmi = 0.45;
4264           TPaveText* txt0 = new TPaveText(0.65, 0.91, 0.90, 0.96, "blNDC");
4265           txt0->SetFillColor(0);
4266           sprintf(txt, "%4.1f TeV %5.1f fb^{-1}", ener, lumi);
4267           txt0->AddText(txt);
4268           txt0->Draw("same");
4269         }
4270         double ymi = ymx - 0.05;
4271         TPaveText* txt1 = new TPaveText(xmi, ymi, xmx, ymx, "blNDC");
4272         txt1->SetFillColor(0);
4273         if (text == "") {
4274           if (k == 0)
4275             sprintf(txt, "%s", title1[j].c_str());
4276           else
4277             sprintf(txt, "%s", title2[j].c_str());
4278         } else {
4279           if (k == 0)
4280             sprintf(txt, "%s (%s)", title1[j].c_str(), text.c_str());
4281           else
4282             sprintf(txt, "%s (%s)", title2[j].c_str(), text.c_str());
4283         }
4284         txt1->AddText(txt);
4285         txt1->Draw("same");
4286         double xmax = (isRealData) ? 0.24 : 0.35;
4287         ymi = 0.91;
4288         ymx = ymi + 0.05;
4289         TPaveText* txt2 = new TPaveText(0.02, ymi, xmax, ymx, "blNDC");
4290         txt2->SetFillColor(0);
4291         if (isRealData)
4292           sprintf(txt, "CMS Preliminary");
4293         else
4294           sprintf(txt, "CMS Simulation Preliminary");
4295         txt2->AddText(txt);
4296         txt2->Draw("same");
4297         pad->Modified();
4298         pad->Update();
4299         if (save > 0) {
4300           sprintf(name, "%s.pdf", pad->GetName());
4301           pad->Print(name);
4302         } else if (save < 0) {
4303           sprintf(name, "%s.C", pad->GetName());
4304           pad->Print(name);
4305         }
4306       }
4307     }
4308   }
4309 }
4310 
4311 void PlotMeanError(const std::string infilest, int reg = 3, bool resol = false, int save = 0, bool debug = false) {
4312   bool ok(false);
4313   const int ntypmx = 3;
4314   const int nregmx = 4;
4315   if (reg < 0 || reg >= nregmx)
4316     reg = nregmx - 1;
4317   int nEner(0), nType(0), nPts(0);
4318   std::vector<double> energy[ntypmx], denergy[ntypmx], value[ntypmx], dvalue[ntypmx];
4319   // First read the data
4320   std::ifstream fInput(infilest.c_str());
4321   if (!fInput.good()) {
4322     std::cout << "Cannot open file " << infilest << std::endl;
4323   } else {
4324     ok = true;
4325     fInput >> nEner >> nType >> nPts;
4326     int nmax = nEner * nType;
4327     int type, elow, ehigh;
4328     double v1[4], e1[4], v2[4], e2[4];
4329     for (int n = 0; n < nmax; ++n) {
4330       fInput >> type >> elow >> ehigh;
4331       fInput >> v1[0] >> e1[0] >> v1[1] >> e1[1] >> v1[2] >> e1[2] >> v1[3] >> e1[3];
4332       fInput >> v2[0] >> e2[0] >> v2[1] >> e2[1] >> v2[2] >> e2[2] >> v2[3] >> e2[3];
4333       double ener = 0.5 * (ehigh + elow);
4334       double dene = 0.5 * (ehigh - elow);
4335       energy[type].push_back(ener);
4336       denergy[type].push_back(dene);
4337       if (resol) {
4338         value[type].push_back(v2[reg]);
4339         dvalue[type].push_back(e2[reg]);
4340       } else {
4341         value[type].push_back(v1[reg]);
4342         dvalue[type].push_back(e1[reg]);
4343       }
4344     }
4345     fInput.close();
4346     std::cout << "Reads " << (nmax + 1) << " cards from " << infilest << " with measurements for " << nEner
4347               << " energies and " << nType << " types" << std::endl;
4348     if (debug) {
4349       for (int n = 0; n < nType; ++n) {
4350         std::cout << "Type " << n << " with " << energy[n].size() << " points\n";
4351         for (unsigned int k = 0; k < energy[n].size(); ++k)
4352           std::cout << " [" << k << "] " << energy[n][k] << " +- " << denergy[n][k] << " Value " << value[n][k]
4353                     << " +- " << dvalue[n][k] << std::endl;
4354       }
4355     }
4356   }
4357 
4358   // Now the plots
4359   if (ok) {
4360     int mvsres = (resol) ? 1 : 0;
4361     std::string names[2] = {"Mean", "Resol"};
4362     std::string namet[nregmx] = {"Barrel", "Transition", "Endcap", "Combined"};
4363     char cname[100];
4364     sprintf(cname, "c_%s%s", names[mvsres].c_str(), namet[reg].c_str());
4365     int color[ntypmx] = {2, 4, 1};
4366     int mtype[ntypmx] = {20, 21, 22};
4367     double ymin[2] = {0.65, 0.10};
4368     double ymax[2] = {1.30, 0.50};
4369     gStyle->SetCanvasBorderMode(0);
4370     gStyle->SetCanvasColor(kWhite);
4371     gStyle->SetPadColor(kWhite);
4372     gStyle->SetFillColor(kWhite);
4373     gStyle->SetOptTitle(kFALSE);
4374     gStyle->SetPadBorderMode(0);
4375     gStyle->SetCanvasBorderMode(0);
4376     gStyle->SetOptStat(0);
4377     TCanvas* canvas = new TCanvas(cname, cname, 500, 500);
4378     canvas->SetTopMargin(0.05);
4379     canvas->SetBottomMargin(0.14);
4380     canvas->SetLeftMargin(0.15);
4381     canvas->SetRightMargin(0.10);
4382     TH1F* vFrame = canvas->DrawFrame(0.0, ymin[mvsres], 120.0, ymax[mvsres]);
4383     vFrame->GetXaxis()->SetRangeUser(0.0, 120.0);
4384     vFrame->GetYaxis()->SetRangeUser(ymin[mvsres], ymax[mvsres]);
4385     vFrame->GetXaxis()->SetLabelSize(0.04);
4386     vFrame->GetYaxis()->SetLabelSize(0.04);
4387     vFrame->GetXaxis()->SetTitleSize(0.04);
4388     vFrame->GetYaxis()->SetTitleSize(0.04);
4389     vFrame->GetXaxis()->SetTitleOffset(1.2);
4390     vFrame->GetYaxis()->SetTitleOffset(1.6);
4391     vFrame->GetXaxis()->SetLabelOffset(0.0);
4392     vFrame->GetXaxis()->SetTitle("p_{Beam} (GeV/c)");
4393     if (resol) {
4394       vFrame->GetYaxis()->SetTitle("Width(E_{HCAL}/(p-E_{ECAL}))");
4395     } else {
4396       vFrame->GetYaxis()->SetTitle("MPV(E_{HCAL}/(p-E_{ECAL}))");
4397     }
4398     TLegend* legend = new TLegend(0.70, 0.80, 0.90, 0.94);
4399     legend->SetFillColor(kWhite);
4400     std::string nameg[ntypmx] = {"MAHI", "M0", "M2"};
4401     for (int n = 0; n < nType; ++n) {
4402       unsigned int nmax0 = energy[n].size();
4403       double mom[nmax0], dmom[nmax0], mean[nmax0], dmean[nmax0];
4404       for (unsigned int k = 0; k < nmax0; ++k) {
4405         mom[k] = energy[n][k];
4406         dmom[k] = denergy[n][k];
4407         mean[k] = value[n][k];
4408         dmean[k] = dvalue[n][k];
4409       }
4410       TGraphErrors* graph = new TGraphErrors(nmax0, mom, mean, dmom, dmean);
4411       graph->SetMarkerStyle(mtype[n]);
4412       graph->SetMarkerColor(color[n]);
4413       graph->SetMarkerSize(1.4);
4414       graph->SetLineColor(color[n]);
4415       graph->SetLineWidth(2);
4416       sprintf(cname, "%s", nameg[n].c_str());
4417       legend->AddEntry(graph, cname, "lp");
4418       graph->Draw("P");
4419     }
4420     legend->Draw("same");
4421     std::string regions[nregmx] = {"20118B Barrel", "2018B Transition", "2018B Endcap", "2018B"};
4422     sprintf(cname, "%s", regions[reg].c_str());
4423     TPaveText* txt0 = new TPaveText(0.16, 0.90, 0.40, 0.94, "blNDC");
4424     txt0->SetFillColor(0);
4425     txt0->AddText(cname);
4426     txt0->Draw("same");
4427     TPaveText* txt1 = new TPaveText(0.15, 0.95, 0.40, 0.99, "blNDC");
4428     txt1->SetFillColor(0);
4429     sprintf(cname, "CMS Preliminary");
4430     txt1->AddText(cname);
4431     txt1->Draw("same");
4432     canvas->Update();
4433     if (save > 0) {
4434       sprintf(cname, "%s.pdf", canvas->GetName());
4435       canvas->Print(cname);
4436     } else if (save < 0) {
4437       sprintf(cname, "%s.C", canvas->GetName());
4438       canvas->Print(cname);
4439     }
4440   }
4441 }
4442 
4443 void PlotDepthCorrFactor(char* infile,
4444                          std::string text,
4445                          std::string prefix = "",
4446                          bool isRealData = true,
4447                          bool drawStatBox = true,
4448                          int save = 0) {
4449   std::map<int, cfactors> cfacs;
4450   int etamin(100), etamax(-100), maxdepth(0);
4451   std::ifstream ifile(infile);
4452   if (!ifile.is_open()) {
4453     std::cout << "Cannot open duplicate file " << infile << std::endl;
4454   } else {
4455     unsigned int all(0), good(0);
4456     char buffer[1024];
4457     while (ifile.getline(buffer, 1024)) {
4458       ++all;
4459       std::string bufferString(buffer);
4460       if (bufferString.substr(0, 1) == "#") {
4461         continue;  //ignore other comments
4462       } else {
4463         std::vector<std::string> items = splitString(bufferString);
4464         if (items.size() < 3) {
4465           std::cout << "Ignore  line: " << buffer << " Size " << items.size();
4466           for (unsigned int k = 0; k < items.size(); ++k)
4467             std::cout << " [" << k << "] : " << items[k];
4468           std::cout << std::endl;
4469         } else {
4470           ++good;
4471           int ieta = std::atoi(items[0].c_str());
4472           if (ieta < etamin)
4473             etamin = ieta;
4474           if (ieta > etamax)
4475             etamax = ieta;
4476           unsigned int k(1);
4477           int depth(0);
4478           std::cout << "ieta " << ieta;
4479           while (k < items.size()) {
4480             ++depth;
4481             double corrf = std::atof(items[k].c_str());
4482             double dcorr = std::atof(items[k + 1].c_str());
4483             if (depth > maxdepth)
4484               maxdepth = depth;
4485             if ((depth == 1) && ((std::abs(ieta) == 17) || (std::abs(ieta) == 18))) {
4486             } else {
4487               int detId = repackId(ieta, depth);
4488               cfacs[detId] = cfactors(ieta, depth, corrf, dcorr);
4489               std::cout << " Depth " << depth << " " << corrf << " +- " << dcorr;
4490             }
4491             k += 2;
4492           }
4493           std::cout << std::endl;
4494         }
4495       }
4496     }
4497     ifile.close();
4498     std::cout << "Reads total of " << all << " and " << good << " good records of depth dependent factors from "
4499               << infile << " Eta range " << etamin << ":" << etamax << " maxdepth " << maxdepth << std::endl;
4500   }
4501 
4502   gStyle->SetCanvasBorderMode(0);
4503   gStyle->SetCanvasColor(kWhite);
4504   gStyle->SetPadColor(kWhite);
4505   gStyle->SetFillColor(kWhite);
4506   gStyle->SetOptTitle(0);
4507   if (drawStatBox) {
4508     gStyle->SetOptStat(10);
4509     gStyle->SetOptFit(10);
4510   } else {
4511     gStyle->SetOptStat(0);
4512     gStyle->SetOptFit(0);
4513   }
4514   int colors[7] = {1, 6, 4, 7, 2, 9, 3};
4515   int mtype[7] = {20, 21, 22, 23, 24, 33, 25};
4516   int nbin = etamax - etamin + 1;
4517   std::vector<TH1D*> hists;
4518   std::vector<int> entries;
4519   char name[100];
4520   double dy(0);
4521   int fits(0);
4522   for (int j = 0; j < maxdepth; ++j) {
4523     sprintf(name, "hd%d", j + 1);
4524     TObject* ob = gROOT->FindObject(name);
4525     if (ob)
4526       ob->Delete();
4527     TH1D* h = new TH1D(name, name, nbin, etamin, etamax);
4528     int nent(0);
4529     for (std::map<int, cfactors>::const_iterator itr = cfacs.begin(); itr != cfacs.end(); ++itr) {
4530       if ((itr->second).depth == j + 1) {
4531         int ieta = (itr->second).ieta;
4532         int bin = ieta - etamin + 1;
4533         float val = (itr->second).corrf;
4534         float dvl = (itr->second).dcorr;
4535         h->SetBinContent(bin, val);
4536         h->SetBinError(bin, dvl);
4537         nent++;
4538       }
4539     }
4540     h->SetLineColor(colors[j]);
4541     h->SetMarkerColor(colors[j]);
4542     h->SetMarkerStyle(mtype[j]);
4543     h->GetXaxis()->SetTitle("i#eta");
4544     h->GetYaxis()->SetTitle("Depth Dependent Correction Factor");
4545     h->GetYaxis()->SetLabelOffset(0.005);
4546     h->GetYaxis()->SetTitleOffset(1.20);
4547     h->GetYaxis()->SetRangeUser(0.0, 2.0);
4548     hists.push_back(h);
4549     entries.push_back(nent);
4550     dy += 0.025;
4551   }
4552   sprintf(name, "c_%sCorrFactor", prefix.c_str());
4553   TCanvas* pad = new TCanvas(name, name, 700, 500);
4554   pad->SetRightMargin(0.10);
4555   pad->SetTopMargin(0.10);
4556   double yh = 0.90;
4557   // double yl = yh - 0.025 * hists.size() - dy - 0.01;
4558   double yl = 0.15;
4559   TLegend* legend = new TLegend(0.35, yl, 0.65, yl + 0.04 * hists.size());
4560   legend->SetFillColor(kWhite);
4561   for (unsigned int k = 0; k < hists.size(); ++k) {
4562     if (k == 0)
4563       hists[k]->Draw("");
4564     else
4565       hists[k]->Draw("sames");
4566     pad->Update();
4567     TPaveStats* st1 = (TPaveStats*)hists[k]->GetListOfFunctions()->FindObject("stats");
4568     if (st1 != nullptr) {
4569       dy = 0.025;
4570       st1->SetLineColor(colors[k]);
4571       st1->SetTextColor(colors[k]);
4572       st1->SetY1NDC(yh - dy);
4573       st1->SetY2NDC(yh);
4574       st1->SetX1NDC(0.70);
4575       st1->SetX2NDC(0.90);
4576       yh -= dy;
4577     }
4578     sprintf(name, "Depth %d (%s)", k + 1, text.c_str());
4579     legend->AddEntry(hists[k], name, "lp");
4580   }
4581   legend->Draw("same");
4582   pad->Update();
4583   if (fits < 1) {
4584     double xmin = hists[0]->GetBinLowEdge(1);
4585     int nbin = hists[0]->GetNbinsX();
4586     double xmax = hists[0]->GetBinLowEdge(nbin) + hists[0]->GetBinWidth(nbin);
4587     TLine* line = new TLine(xmin, 1.0, xmax, 1.0);
4588     line->SetLineColor(9);
4589     line->SetLineWidth(2);
4590     line->SetLineStyle(2);
4591     line->Draw("same");
4592     pad->Modified();
4593     pad->Update();
4594   }
4595   char txt1[30];
4596   double xmax = (isRealData) ? 0.33 : 0.44;
4597   TPaveText* txt2 = new TPaveText(0.11, 0.85, xmax, 0.89, "blNDC");
4598   txt2->SetFillColor(0);
4599   if (isRealData)
4600     sprintf(txt1, "CMS Preliminary");
4601   else
4602     sprintf(txt1, "CMS Simulation Preliminary");
4603   txt2->AddText(txt1);
4604   txt2->Draw("same");
4605   pad->Modified();
4606   pad->Update();
4607   if (save > 0) {
4608     sprintf(name, "%s.pdf", pad->GetName());
4609     pad->Print(name);
4610   } else if (save < 0) {
4611     sprintf(name, "%s.C", pad->GetName());
4612     pad->Print(name);
4613   }
4614 }
4615 
4616 void DrawHistPhiSymmetry(TH1D* hist0, bool isRealData, bool drawStatBox, bool save) {
4617   char name[30], namep[30], txt1[30];
4618   TH1D* hist = (TH1D*)(hist0->Clone());
4619   sprintf(namep, "c_%s", hist->GetName());
4620   TCanvas* pad = new TCanvas(namep, namep, 700, 500);
4621   pad->SetRightMargin(0.10);
4622   pad->SetTopMargin(0.10);
4623   hist->GetXaxis()->SetTitleSize(0.04);
4624   hist->GetXaxis()->SetTitle(hist->GetTitle());
4625   hist->GetYaxis()->SetTitle("Channels");
4626   hist->GetYaxis()->SetLabelOffset(0.005);
4627   hist->GetYaxis()->SetTitleSize(0.04);
4628   hist->GetYaxis()->SetLabelSize(0.035);
4629   hist->GetYaxis()->SetTitleOffset(1.10);
4630   hist->Draw();
4631   pad->Update();
4632   if (drawStatBox) {
4633     TPaveStats* st1 = (TPaveStats*)hist->GetListOfFunctions()->FindObject("stats");
4634     if (st1 != nullptr) {
4635       st1->SetY1NDC(0.70);
4636       st1->SetY2NDC(0.90);
4637       st1->SetX1NDC(0.65);
4638       st1->SetX2NDC(0.90);
4639     }
4640     pad->Modified();
4641     pad->Update();
4642   }
4643   TPaveText* txt2 = new TPaveText(0.11, 0.85, 0.44, 0.89, "blNDC");
4644   txt2->SetFillColor(0);
4645   if (isRealData)
4646     sprintf(txt1, "CMS Preliminary");
4647   else
4648     sprintf(txt1, "CMS Simulation Preliminary");
4649   txt2->AddText(txt1);
4650   txt2->Draw("same");
4651   pad->Modified();
4652   pad->Update();
4653   if (save) {
4654     sprintf(name, "%s.pdf", pad->GetName());
4655     pad->Print(name);
4656   }
4657 }
4658 
4659 void PlotPhiSymmetryResults(
4660     char* infile, bool isRealData = true, bool drawStatBox = true, bool debug = false, bool save = false) {
4661   const int maxDepthHB(4), maxDepthHE(7);
4662   const double cfacMin(0.70), cfacMax(1.5);
4663   const int nbin = (100.0 * (cfacMax - cfacMin));
4664   std::vector<TH1D*> histHB, histHE;
4665   char name[20], title[80];
4666   for (int k = 0; k < maxDepthHB; ++k) {
4667     sprintf(name, "HB%d", k);
4668     sprintf(title, "Correction factor for depth %d of HB", k);
4669     TObject* ob = gROOT->FindObject(name);
4670     if (ob)
4671       ob->Delete();
4672     TH1D* h = new TH1D(name, title, nbin, cfacMin, cfacMax);
4673     histHB.push_back(h);
4674     if (debug)
4675       std::cout << "Book " << h->GetName() << " Title " << h->GetTitle() << " range " << nbin << ":" << cfacMin << ":"
4676                 << cfacMax << std::endl;
4677   }
4678   for (int k = 0; k < maxDepthHE; ++k) {
4679     sprintf(name, "HE%d", k);
4680     sprintf(title, "Correction factor for depth %d of HE", k);
4681     TObject* ob = gROOT->FindObject(name);
4682     if (ob)
4683       ob->Delete();
4684     TH1D* h = new TH1D(name, title, nbin, cfacMin, cfacMax);
4685     histHE.push_back(h);
4686     if (debug)
4687       std::cout << "Book " << h->GetName() << " Title " << h->GetTitle() << " range " << nbin << ":" << cfacMin << ":"
4688                 << cfacMax << std::endl;
4689   }
4690   std::ifstream fInput(infile);
4691   if (!fInput.good()) {
4692     std::cout << "Cannot open file " << infile << std::endl;
4693   } else {
4694     char buffer[1024];
4695     int all(0), good(0);
4696     std::map<int, int> kount;
4697     while (fInput.getline(buffer, 1024)) {
4698       ++all;
4699       std::string bufferString(buffer);
4700       if (bufferString.substr(0, 1) == "#") {
4701         continue;  //ignore other comments
4702       } else {
4703         std::vector<std::string> items = splitString(bufferString);
4704         if (items.size() < 5) {
4705           for (unsigned int k = 0; k < items.size(); ++k)
4706             std::cout << " [" << k << "] : " << items[k];
4707           std::cout << std::endl;
4708         } else {
4709           ++good;
4710           int subdet = std::atoi(items[0].c_str());
4711           int ieta = std::atoi(items[1].c_str());
4712           int depth = std::atoi(items[3].c_str());
4713           double corrf = std::atof(items[4].c_str());
4714           if ((debug) && (good % 100 == 0))
4715             std::cout << "Subdet " << subdet << " Depth " << depth << " Corr " << corrf << std::endl;
4716           if ((corrf < cfacMin) || (corrf > cfacMax))
4717             std::cout << "Outlier: " << buffer << std::endl;
4718           if (subdet == 1) {
4719             if (depth <= maxDepthHB)
4720               histHB[depth - 1]->Fill(corrf);
4721           } else if (subdet == 2) {
4722             if (depth <= maxDepthHE)
4723               histHE[depth - 1]->Fill(corrf);
4724           }
4725           if ((subdet == 1) || (subdet == 2)) {
4726             int indx = (ieta > 0) ? (ieta * 10 + depth) : (10000 + 10 * (-ieta) + depth);
4727             if (kount.find(indx) == kount.end())
4728               kount[indx] = 1;
4729             else
4730               ++kount[indx];
4731           }
4732         }
4733       }
4734     }
4735     fInput.close();
4736     std::cout << "Reads total of " << all << " and " << good << " good records of phi-symmetry factors from " << infile
4737               << std::endl;
4738     for (std::map<int, int>::iterator itr = kount.begin(); itr != kount.end(); ++itr) {
4739       int depth = (itr->first) % 10;
4740       int ieta = ((itr->first) / 10) % 100;
4741       int nphi = itr->second;
4742       bool miss = (((ieta <= 20) && (nphi != 72)) || ((ieta > 20) && (nphi != 36)));
4743       if (itr->first >= 10000)
4744         ieta = -ieta;
4745       if (debug || miss)
4746         std::cout << "Tower[" << ieta << ":" << depth << "] has " << nphi << " entries" << std::endl;
4747     }
4748   }
4749 
4750   // Now plot the histograms
4751   gStyle->SetCanvasBorderMode(0);
4752   gStyle->SetCanvasColor(kWhite);
4753   gStyle->SetPadColor(kWhite);
4754   gStyle->SetFillColor(kWhite);
4755   gStyle->SetOptTitle(0);
4756   gStyle->SetOptStat(111110);
4757   gStyle->SetOptFit(0);
4758 
4759   // HB first
4760   for (unsigned int k = 0; k < histHB.size(); ++k) {
4761     DrawHistPhiSymmetry(histHB[k], isRealData, drawStatBox, save);
4762   }
4763 
4764   // Then HE
4765   for (unsigned int k = 0; k < histHE.size(); ++k) {
4766     DrawHistPhiSymmetry(histHE[k], isRealData, drawStatBox, save);
4767   }
4768 }
4769 
4770 void PlotHistCorrRatio(char* infile1,
4771                        std::string text1,
4772                        char* infile2,
4773                        std::string text2,
4774                        int depth1,
4775                        int depth2,
4776                        std::string prefixF,
4777                        std::string text0,
4778                        int etaMin = -1,
4779                        int etaMax = -1,
4780                        bool doFit = true,
4781                        bool isRealData = true,
4782                        int year = 2022,
4783                        int iformat = 0,
4784                        int save = 0) {
4785   std::map<int, cfactors> cfacs[2];
4786   std::vector<std::string> texts;
4787   int nfile(0), etamin(100), etamax(-100), maxdepth(0);
4788   const char* blank("");
4789   if (infile1 != blank) {
4790     readCorrFactors(infile1, 1.0, cfacs[nfile], etamin, etamax, maxdepth, iformat);
4791     if (cfacs[nfile].size() > 0) {
4792       texts.push_back(text1);
4793       ++nfile;
4794     }
4795   }
4796   if (infile2 != blank) {
4797     readCorrFactors(infile2, 1.0, cfacs[nfile], etamin, etamax, maxdepth, iformat);
4798     if (cfacs[nfile].size() > 0) {
4799       texts.push_back(text2);
4800       ++nfile;
4801     }
4802   }
4803 
4804   if (etaMax > 0) {
4805     etamin = -etaMax;
4806     etamax = etaMax;
4807   }
4808   if (nfile == 2) {
4809     gStyle->SetCanvasBorderMode(0);
4810     gStyle->SetCanvasColor(kWhite);
4811     gStyle->SetPadColor(kWhite);
4812     gStyle->SetFillColor(kWhite);
4813     gStyle->SetOptTitle(0);
4814     if (doFit) {
4815       gStyle->SetOptStat(10);
4816       gStyle->SetOptFit(10);
4817     } else {
4818       gStyle->SetOptStat(0);
4819       gStyle->SetOptFit(0);
4820     }
4821     int colors[7] = {1, 6, 4, 2, 7, 9, 46};
4822     int mtype[7] = {20, 24, 22, 23, 21, 25, 33};
4823     int styles[7] = {2, 3, 1, 4, 1, 3, 2};
4824     int nbin = etamax - etamin + 1;
4825     std::vector<TH1D*> hists;
4826     std::vector<double> fitr, dfit;
4827     char name[100];
4828     for (int ih = 0; ih < nfile; ++ih) {
4829       sprintf(name, "h%d", ih);
4830       TObject* ob = gROOT->FindObject(name);
4831       if (ob)
4832         ob->Delete();
4833       TH1D* h = new TH1D(name, name, nbin, etamin, etamax);
4834       double sumNum(0), sumDen(0);
4835       int npt(0);
4836       for (std::map<int, cfactors>::const_iterator itr = cfacs[ih].begin(); itr != cfacs[ih].end(); ++itr) {
4837         int ieta = (itr->second).ieta;
4838         bool seleta = (etaMin > 0) ? (std::abs(ieta) > etaMin) : true;
4839         if ((ieta >= etamin) && (ieta <= etamax) && seleta && ((itr->second).depth == depth1)) {
4840           ++npt;
4841           int bin = ieta - etamin + 1;
4842           for (std::map<int, cfactors>::const_iterator ktr = cfacs[ih].begin(); ktr != cfacs[ih].end(); ++ktr) {
4843             if (((ktr->second).ieta == ieta) && ((ktr->second).depth == depth2)) {
4844               double er1 = (itr->second).dcorr / (itr->second).corrf;
4845               double er2 = (ktr->second).dcorr / (ktr->second).corrf;
4846               float val = (itr->second).corrf / (ktr->second).corrf;
4847               float dvl = val * sqrt(er1 * er1 + er2 * er2);
4848               double temp1 = ((itr->second).corrf > 1.0) ? 1.0 / (itr->second).corrf : (itr->second).corrf;
4849               double temp2 = ((itr->second).corrf > 1.0)
4850                                  ? (itr->second).dcorr / ((itr->second).corrf * (itr->second).corrf)
4851                                  : (itr->second).dcorr;
4852               h->SetBinContent(bin, val);
4853               h->SetBinError(bin, dvl);
4854               sumNum += (std::abs(1 - temp1) / (temp2 * temp2));
4855               sumDen += (1.0 / (temp2 * temp2));
4856               break;
4857             }
4858           }
4859         }
4860       }
4861       h->SetLineColor(colors[ih]);
4862       h->SetMarkerColor(colors[ih]);
4863       h->SetMarkerStyle(mtype[ih]);
4864       h->SetMarkerSize(0.9);
4865       h->GetXaxis()->SetTitle("i#eta");
4866       sprintf(name, "CF_{%d}/CF_{%d}", depth1, depth2);
4867       h->GetYaxis()->SetTitle(name);
4868       h->GetYaxis()->SetLabelOffset(0.005);
4869       h->GetYaxis()->SetTitleSize(0.036);
4870       h->GetYaxis()->SetTitleOffset(1.20);
4871       h->GetYaxis()->SetRangeUser(0.0, 3.0);
4872       if (doFit) {
4873         TObject* ob = gROOT->FindObject(name);
4874         if (ob)
4875           ob->Delete();
4876         TF1* func = new TF1(name, "pol0", etamin, etamax);
4877         func->SetLineColor(colors[ih]);
4878         func->SetLineStyle(styles[ih]);
4879         h->Fit(func, "+QWLR", "");
4880       }
4881       hists.push_back(h);
4882       sumNum = (sumDen > 0) ? (sumNum / sumDen) : 0;
4883       sumDen = (sumDen > 0) ? 1.0 / sqrt(sumDen) : 0;
4884       fitr.push_back(sumNum);
4885       dfit.push_back(sumDen);
4886       std::cout << "Get Ratio of mean for " << npt << " points: Mean " << sumNum << " +- " << sumDen << std::endl;
4887     }
4888     sprintf(name, "c_Corr%sRatio", prefixF.c_str());
4889     TCanvas* pad = new TCanvas(name, name, 700, 500);
4890     pad->SetRightMargin(0.10);
4891     pad->SetTopMargin(0.10);
4892     double yh = 0.90;
4893     double yl = yh - 0.035 * hists.size() - 0.01;
4894     TLegend* legend = new TLegend(0.11, yl, 0.50, yh - 0.01);
4895     legend->SetFillColor(kWhite);
4896     for (unsigned int k = 0; k < hists.size(); ++k) {
4897       if (k == 0)
4898         hists[k]->Draw("");
4899       else
4900         hists[k]->Draw("sames");
4901       pad->Update();
4902       if (doFit) {
4903         TPaveStats* st1 = (TPaveStats*)hists[k]->GetListOfFunctions()->FindObject("stats");
4904         if (st1 != nullptr) {
4905           st1->SetLineColor(colors[k]);
4906           st1->SetTextColor(colors[k]);
4907           yh = 0.90 - 0.070 * k;
4908           st1->SetY1NDC(yh - 0.07);
4909           st1->SetY2NDC(yh);
4910           st1->SetX1NDC(0.65);
4911           st1->SetX2NDC(0.90);
4912         }
4913       }
4914       pad->Update();
4915       sprintf(name, "%s (Mean dev. = %5.3f)", texts[k].c_str(), fitr[k]);
4916       legend->AddEntry(hists[k], name, "lp");
4917     }
4918     legend->Draw("same");
4919     TPaveText* txt0 = new TPaveText(0.12, 0.91, 0.49, 0.96, "blNDC");
4920     txt0->SetFillColor(0);
4921     char txt[40];
4922     if (isRealData)
4923       sprintf(txt, "CMS Preliminary (%d)", year);
4924     else
4925       sprintf(txt, "CMS Simulation Preliminary (%d)", year);
4926     txt0->AddText(txt);
4927     txt0->Draw("same");
4928     TPaveText* txt2 = new TPaveText(0.65, 0.91, 0.90, 0.96, "blNDC");
4929     txt2->SetFillColor(0);
4930     txt2->AddText(text0.c_str());
4931     txt2->Draw("same");
4932     pad->Update();
4933     if (save > 0) {
4934       sprintf(name, "%s.pdf", pad->GetName());
4935       pad->Print(name);
4936     } else if (save < 0) {
4937       sprintf(name, "%s.C", pad->GetName());
4938       pad->Print(name);
4939     }
4940   }
4941 }