Back to home page

Project CMSSW displayed by LXR

 
 

    


File indexing completed on 2022-11-02 03:49:48

0001 //////////////////////////////////////////////////////////////////////////////
0002 // Usage:
0003 // .L CalibMonitor.C+g
0004 //  CalibMonitor c1(fname, dirname, dupFileName, comFileName, outFileName,
0005 //                  prefix, corrFileName, rcorFileName, puCorr, flag, numb,
0006 //                  dataMC, truncateFlag, useGen, scale, useScale, etalo, etahi,
0007 //                  runlo, runhi, phimin, phimax, zside, nvxlo, nvxhi, rbx,
0008 //                  exclude, etamax);
0009 //  c1.Loop(nmax);
0010 //  c1.savePlot(histFileName,append,all);
0011 //
0012 //        This will prepare a set of histograms which can be used for a
0013 //        quick fit and display using the methods in CalibFitPlots.C
0014 //
0015 //  GetEntries g1(fname, dirname, dupFileName, bit1, bit2);
0016 //  g1.Loop(nmax);
0017 //
0018 //         This looks into the tree *EventInfo* and can provide a set
0019 //         of histograms with event statistics
0020 //
0021 //   where:
0022 //
0023 //   fname   (const char*)     = file name of the input ROOT tree
0024 //                               or name of the file containing a list of
0025 //                               file names of input ROOT trees
0026 //   dirname (const char*)     = name of the directory where Tree resides
0027 //                               (use "HcalIsoTrkAnalyzer")
0028 //   dupFileName (char*)       = name of the file containing list of entries
0029 //                               of duplicate events
0030 //   comFileName (char*)       = name of the file with list of run and event
0031 //                               number to be selected
0032 //   outFileName (char*)       = name of a text file to be created (under
0033 //                               control of value of flag) with information
0034 //                               about events
0035 //   prefix (std::string)      = String to be added to the name of histogram
0036 //                               (usually a 4 character string; default="")
0037 //   corrFileName (char*)      = name of the text file having the correction
0038 //                               factors to be used (default="", no corr.)
0039 //   rcorFileName (char*)      = name of the text file having the correction
0040 //                               factors as a function of run numbers or depth
0041 //                               or entry number to be used for raddam/depth/
0042 //                               pileup dependent correction  (default="",
0043 //                               no corr.)
0044 //   puCorr (int)              = PU correction to be applied or not: 0 no
0045 //                               correction; < 0 use eDelta; > 0 rho dependent
0046 //                               correction (-8)
0047 //   flag (int)                = 7 digit integer (xmlthdo) with control
0048 //                               information (x=3/2/1/0 for having 1000/500/50/
0049 //                               100 bins for response distribution in (0:5);
0050 //                               m=1/0 for (not) making plots for each RBX;
0051 //                               l=3/2/1/0 for type of rcorFileName (3 for
0052 //                               pileup correction using machine learning
0053 //                               method; 2 for overall response corrections;
0054 //                               1 for depth dependence corrections;
0055 //                               0 for raddam corrections);
0056 //                               t = bit information (lower bit set will
0057 //                               apply a cut on L1 closeness; and higher bit
0058 //                               set read correction file with Marina format);
0059 //                               h = 0/1/2 for not creating / creating in
0060 //                               recreate mode / creating in append mode
0061 //                               the output text file;
0062 //                               d = 0/1/2/3 produces 3 standard (0,1,2) or
0063 //                               extended (3) set of histograms;
0064 //                               o = 0/1/2 for tight / loose / flexible
0065 //                               selection). Default = 1031
0066 //   numb   (int)              = number of eta bins (50 for -25:25)
0067 //   dataMC (bool)             = true/false for data/MC (default true)
0068 //   truncateFlag    (int)     = Flag to treat different depths differently (0)
0069 //                               both depths of ieta 15, 16 of HB as depth 1 (1)
0070 //                               all depths as depth 1 (2), all depths in HE
0071 //                               with values > 1 as depth 2 (3), all depths in
0072 //                               HB with values > 1 as depth 2 (4), all depths
0073 //                               in HB and HE with values > 1 as depth 2 (5)
0074 //                               (default = 0)
0075 //   useGen (bool)             = true/false to use generator level momentum
0076 //                               or reconstruction level momentum
0077 //                               (default = false)
0078 //   scale (double)            = energy scale if correction factor to be used
0079 //                               (default = 1.0)
0080 //   useScale (int)            = application of scale factor (0: nowehere,
0081 //                               1: barrel; 2: endcap, 3: everywhere)
0082 //                               barrel => |ieta| < 16; endcap => |ieta| > 15
0083 //                               (default = 0)
0084 //   etalo/etahi (int,int)     = |eta| ranges (default = 0:30)
0085 //   runlo  (int)              = lower value of run number to be included (+ve)
0086 //                               or excluded (-ve) (default = 0)
0087 //   runhi  (int)              = higher value of run number to be included
0088 //                               (+ve) or excluded (-ve) (default = 9999999)
0089 //   phimin          (int)     = minimum iphi value (default = 1)
0090 //   phimax          (int)     = maximum iphi value (default = 72)
0091 //   zside           (int)     = the side of the detector if phimin and phimax
0092 //                               differ from 1-72 (default = 1)
0093 //   nvxlo           (int)     = minimum # of vertex in event to be used
0094 //                               (default = 0)
0095 //   nvxhi           (int)     = maximum # of vertex in event to be used
0096 //                               (default = 1000)
0097 //   rbx             (int)     = zside*(Subdet*100+RBX #) to be consdered
0098 //                               (default = 0). For HEP17 it will be 217
0099 //   exclude         (bool)    = RBX specified by *rbx* to be exluded or only
0100 //                               considered (default = false)
0101 //   etamax          (bool)    = if set and if the corr-factor not found in the
0102 //                               corrFactor table, the corr-factor for the
0103 //                               corresponding zside, depth=1 and maximum ieta
0104 //                               in the table is taken (default = false)
0105 //   nmax            (Long64_t)= maximum number of entries to be processed,
0106 //                               if -1, all entries to be processed (-1)
0107 //
0108 //   histFileName (std::string)= name of the file containing saved histograms
0109 //   append (bool)             = true/false if the histogram file to be opened
0110 //                               in append/output mode (default = true)
0111 //   all (bool)                = true/false if all histograms to be saved or
0112 //                               not (default = false)
0113 //
0114 //   bitX (unsigned int)       = bit number of the HLT used in the selection
0115 //        (X = 1, 2)             for example the bits of HLT_IsoTrackHB(HE)
0116 //////////////////////////////////////////////////////////////////////////////
0117 
0118 #include <TROOT.h>
0119 #include <TChain.h>
0120 #include <TFile.h>
0121 #include <TH1D.h>
0122 #include <TH2F.h>
0123 #include <TProfile.h>
0124 #include <TFitResult.h>
0125 #include <TFitResultPtr.h>
0126 #include <TStyle.h>
0127 #include <TCanvas.h>
0128 #include <TLegend.h>
0129 #include <TPaveStats.h>
0130 #include <TPaveText.h>
0131 
0132 #include <algorithm>
0133 #include <iomanip>
0134 #include <iostream>
0135 #include <fstream>
0136 #include <map>
0137 #include <vector>
0138 #include <string>
0139 
0140 void unpackDetId(unsigned int, int &, int &, int &, int &, int &);
0141 #include "CalibCorr.C"
0142 
0143 class CalibMonitor {
0144 public:
0145   TChain *fChain;  //!pointer to the analyzed TTree or TChain
0146   Int_t fCurrent;  //!current Tree number in a TChain
0147 
0148   // Declaration of leaf types
0149   Int_t t_Run;
0150   Int_t t_Event;
0151   Int_t t_DataType;
0152   Int_t t_ieta;
0153   Int_t t_iphi;
0154   Double_t t_EventWeight;
0155   Int_t t_nVtx;
0156   Int_t t_nTrk;
0157   Int_t t_goodPV;
0158   Double_t t_l1pt;
0159   Double_t t_l1eta;
0160   Double_t t_l1phi;
0161   Double_t t_l3pt;
0162   Double_t t_l3eta;
0163   Double_t t_l3phi;
0164   Double_t t_p;
0165   Double_t t_pt;
0166   Double_t t_phi;
0167   Double_t t_mindR1;
0168   Double_t t_mindR2;
0169   Double_t t_eMipDR;
0170   Double_t t_eHcal;
0171   Double_t t_eHcal10;
0172   Double_t t_eHcal30;
0173   Double_t t_hmaxNearP;
0174   Double_t t_rhoh;
0175   Bool_t t_selectTk;
0176   Bool_t t_qltyFlag;
0177   Bool_t t_qltyMissFlag;
0178   Bool_t t_qltyPVFlag;
0179   Double_t t_gentrackP;
0180   std::vector<unsigned int> *t_DetIds;
0181   std::vector<double> *t_HitEnergies;
0182   std::vector<bool> *t_trgbits;
0183   std::vector<unsigned int> *t_DetIds1;
0184   std::vector<unsigned int> *t_DetIds3;
0185   std::vector<double> *t_HitEnergies1;
0186   std::vector<double> *t_HitEnergies3;
0187 
0188   // List of branches
0189   TBranch *b_t_Run;           //!
0190   TBranch *b_t_Event;         //!
0191   TBranch *b_t_DataType;      //!
0192   TBranch *b_t_ieta;          //!
0193   TBranch *b_t_iphi;          //!
0194   TBranch *b_t_EventWeight;   //!
0195   TBranch *b_t_nVtx;          //!
0196   TBranch *b_t_nTrk;          //!
0197   TBranch *b_t_goodPV;        //!
0198   TBranch *b_t_l1pt;          //!
0199   TBranch *b_t_l1eta;         //!
0200   TBranch *b_t_l1phi;         //!
0201   TBranch *b_t_l3pt;          //!
0202   TBranch *b_t_l3eta;         //!
0203   TBranch *b_t_l3phi;         //!
0204   TBranch *b_t_p;             //!
0205   TBranch *b_t_pt;            //!
0206   TBranch *b_t_phi;           //!
0207   TBranch *b_t_mindR1;        //!
0208   TBranch *b_t_mindR2;        //!
0209   TBranch *b_t_eMipDR;        //!
0210   TBranch *b_t_eHcal;         //!
0211   TBranch *b_t_eHcal10;       //!
0212   TBranch *b_t_eHcal30;       //!
0213   TBranch *b_t_hmaxNearP;     //!
0214   TBranch *b_t_rhoh;          //!
0215   TBranch *b_t_selectTk;      //!
0216   TBranch *b_t_qltyFlag;      //!
0217   TBranch *b_t_qltyMissFlag;  //!
0218   TBranch *b_t_qltyPVFlag;    //!
0219   TBranch *b_t_gentrackP;     //!
0220   TBranch *b_t_DetIds;        //!
0221   TBranch *b_t_HitEnergies;   //!
0222   TBranch *b_t_trgbits;       //!
0223   TBranch *b_t_DetIds1;       //!
0224   TBranch *b_t_DetIds3;       //!
0225   TBranch *b_t_HitEnergies1;  //!
0226   TBranch *b_t_HitEnergies3;  //!
0227 
0228   struct counter {
0229     static const int npsize = 4;
0230     counter() {
0231       total = 0;
0232       for (int k = 0; k < npsize; ++k)
0233         count[k] = 0;
0234     };
0235     unsigned int total, count[npsize];
0236   };
0237 
0238   CalibMonitor(const char *fname,
0239                const char *dirname,
0240                const char *dupFileName,
0241                const char *comFileName,
0242                const char *outFileName,
0243                const std::string &prefix = "",
0244                const char *corrFileName = "",
0245                const char *rcorFileName = "",
0246                int puCorr = -8,
0247                int flag = 1031,
0248                int numb = 50,
0249                bool datMC = true,
0250                int truncateFlag = 0,
0251                bool useGen = false,
0252                double scale = 1.0,
0253                int useScale = 0,
0254                int etalo = 0,
0255                int etahi = 30,
0256                int runlo = 0,
0257                int runhi = 99999999,
0258                int phimin = 1,
0259                int phimax = 72,
0260                int zside = 1,
0261                int nvxlo = 0,
0262                int nvxhi = 1000,
0263                int rbx = 0,
0264                bool exclude = false,
0265                bool etamax = false);
0266   virtual ~CalibMonitor();
0267   virtual Int_t Cut(Long64_t entry);
0268   virtual Int_t GetEntry(Long64_t entry);
0269   virtual Long64_t LoadTree(Long64_t entry);
0270   virtual void Init(TChain *, const char *, const char *, const char *);
0271   virtual void Loop(Long64_t nmax = -1);
0272   virtual Bool_t Notify();
0273   virtual void Show(Long64_t entry = -1);
0274   bool goodTrack(double &eHcal, double &cut, const Long64_t &entry, bool debug);
0275   bool selectPhi(bool debug);
0276   void plotHist(int type, int num, bool save = false);
0277   template <class Hist>
0278   void drawHist(Hist *, TCanvas *);
0279   void savePlot(const std::string &theName, bool append = true, bool all = false);
0280   void correctEnergy(double &ener, const Long64_t &entry);
0281 
0282 private:
0283   static const unsigned int npbin = 6, kp50 = 3;
0284   CalibCorrFactor *corrFactor_;
0285   CalibCorr *cFactor_;
0286   CalibSelectRBX *cSelect_;
0287   const std::string fname_, dirnm_, prefix_, outFileName_;
0288   const int corrPU_, flag_, numb_;
0289   const bool dataMC_, useGen_;
0290   const int truncateFlag_;
0291   const int etalo_, etahi_;
0292   int runlo_, runhi_;
0293   const int phimin_, phimax_, zside_, nvxlo_, nvxhi_, rbx_;
0294   bool exclude_, corrE_, cutL1T_, selRBX_;
0295   bool includeRun_;
0296   int coarseBin_, etamp_, etamn_, plotType_;
0297   int flexibleSelect_, ifDepth_;
0298   double log2by18_;
0299   std::ofstream fileout_;
0300   std::vector<Long64_t> entries_;
0301   std::vector<std::pair<int, int> > events_;
0302   std::vector<double> etas_, ps_, dl1_;
0303   std::vector<int> nvx_, ietasL_, ietasH_;
0304   std::vector<TH1D *> h_rbx, h_etaF[npbin], h_etaB[npbin];
0305   std::vector<TProfile *> h_etaX[npbin];
0306   std::vector<TH1D *> h_etaR[npbin], h_nvxR[npbin], h_dL1R[npbin];
0307   std::vector<TH1D *> h_pp[npbin];
0308 };
0309 
0310 CalibMonitor::CalibMonitor(const char *fname,
0311                            const char *dirnm,
0312                            const char *dupFileName,
0313                            const char *comFileName,
0314                            const char *outFName,
0315                            const std::string &prefix,
0316                            const char *corrFileName,
0317                            const char *rcorFileName,
0318                            int puCorr,
0319                            int flag,
0320                            int numb,
0321                            bool dataMC,
0322                            int truncate,
0323                            bool useGen,
0324                            double scale,
0325                            int useScale,
0326                            int etalo,
0327                            int etahi,
0328                            int runlo,
0329                            int runhi,
0330                            int phimin,
0331                            int phimax,
0332                            int zside,
0333                            int nvxlo,
0334                            int nvxhi,
0335                            int rbx,
0336                            bool exc,
0337                            bool etam)
0338     : corrFactor_(nullptr),
0339       cFactor_(nullptr),
0340       cSelect_(nullptr),
0341       fname_(std::string(fname)),
0342       dirnm_(std::string(dirnm)),
0343       prefix_(prefix),
0344       outFileName_(std::string(outFName)),
0345       corrPU_(puCorr),
0346       flag_(flag),
0347       numb_(numb),
0348       dataMC_(dataMC),
0349       useGen_(useGen),
0350       truncateFlag_(truncate),
0351       etalo_(etalo),
0352       etahi_(etahi),
0353       runlo_(runlo),
0354       runhi_(runhi),
0355       phimin_(phimin),
0356       phimax_(phimax),
0357       zside_(zside),
0358       nvxlo_(nvxlo),
0359       nvxhi_(nvxhi),
0360       rbx_(rbx),
0361       exclude_(exc),
0362       includeRun_(true) {
0363   // if parameter tree is not specified (or zero), connect the file
0364   // used to generate this class and read the Tree
0365 
0366   plotType_ = ((flag_ / 10) % 10);
0367   if (plotType_ < 0 || plotType_ > 3)
0368     plotType_ = 3;
0369   flexibleSelect_ = (((flag_ / 1) % 10));
0370   int oneplace = ((flag_ / 1000) % 10);
0371   cutL1T_ = (oneplace % 2);
0372   bool marina = ((oneplace / 2) % 2);
0373   ifDepth_ = ((flag_ / 10000) % 10);
0374   selRBX_ = (((flag_ / 100000) % 10) > 0);
0375   coarseBin_ = ((flag_ / 1000000) % 10);
0376   log2by18_ = std::log(2.5) / 18.0;
0377   if (runlo_ < 0 || runhi_ < 0) {
0378     runlo_ = std::abs(runlo_);
0379     runhi_ = std::abs(runhi_);
0380     includeRun_ = false;
0381   }
0382   char treeName[400];
0383   sprintf(treeName, "%s/CalibTree", dirnm_.c_str());
0384   TChain *chain = new TChain(treeName);
0385   std::cout << "Create a chain for " << treeName << " from " << fname << " flags " << flexibleSelect_ << "|"
0386             << plotType_ << "|" << corrPU_ << "\n cons " << log2by18_ << " eta range " << etalo_ << ":" << etahi_
0387             << " run range " << runlo_ << ":" << runhi_ << " (inclusion flag " << includeRun_ << ")\n Selection of RBX "
0388             << selRBX_ << " Vertex Range " << nvxlo_ << ":" << nvxhi_ << "\n corrFileName: " << corrFileName
0389             << " useScale " << useScale << ":" << scale << ":" << etam << "\n rcorFileName: " << rcorFileName
0390             << " flag " << ifDepth_ << ":" << cutL1T_ << ":" << marina << std::endl;
0391   if (!fillChain(chain, fname)) {
0392     std::cout << "*****No valid tree chain can be obtained*****" << std::endl;
0393   } else {
0394     std::cout << "Proceed with a tree chain with " << chain->GetEntries() << " entries" << std::endl;
0395     corrFactor_ = new CalibCorrFactor(corrFileName, useScale, scale, etam, marina, false);
0396     Init(chain, dupFileName, comFileName, outFName);
0397     if (std::string(rcorFileName) != "") {
0398       cFactor_ = new CalibCorr(rcorFileName, ifDepth_, false);
0399     } else {
0400       ifDepth_ = 0;
0401     }
0402     if (rbx != 0)
0403       cSelect_ = new CalibSelectRBX(rbx, false);
0404   }
0405 }
0406 
0407 CalibMonitor::~CalibMonitor() {
0408   delete corrFactor_;
0409   delete cFactor_;
0410   delete cSelect_;
0411   if (!fChain)
0412     return;
0413   delete fChain->GetCurrentFile();
0414 }
0415 
0416 Int_t CalibMonitor::GetEntry(Long64_t entry) {
0417   // Read contents of entry.
0418   if (!fChain)
0419     return 0;
0420   return fChain->GetEntry(entry);
0421 }
0422 
0423 Long64_t CalibMonitor::LoadTree(Long64_t entry) {
0424   // Set the environment to read one entry
0425   if (!fChain)
0426     return -5;
0427   Long64_t centry = fChain->LoadTree(entry);
0428   if (centry < 0)
0429     return centry;
0430   if (!fChain->InheritsFrom(TChain::Class()))
0431     return centry;
0432   TChain *chain = (TChain *)fChain;
0433   if (chain->GetTreeNumber() != fCurrent) {
0434     fCurrent = chain->GetTreeNumber();
0435     Notify();
0436   }
0437   return centry;
0438 }
0439 
0440 void CalibMonitor::Init(TChain *tree, const char *dupFileName, const char *comFileName, const char *outFileName) {
0441   // The Init() function is called when the selector needs to initialize
0442   // a new tree or chain. Typically here the branch addresses and branch
0443   // pointers of the tree will be set.
0444   // It is normally not necessary to make changes to the generated
0445   // code, but the routine can be extended by the user if needed.
0446   // Init() will be called many times when running on PROOF
0447   // (once per file to be processed).
0448 
0449   // Set object pointer
0450   t_DetIds = 0;
0451   t_DetIds1 = 0;
0452   t_DetIds3 = 0;
0453   t_HitEnergies = 0;
0454   t_HitEnergies1 = 0;
0455   t_HitEnergies3 = 0;
0456   t_trgbits = 0;
0457   // Set branch addresses and branch pointers
0458   fChain = tree;
0459   fCurrent = -1;
0460   if (!tree)
0461     return;
0462   fChain->SetMakeClass(1);
0463 
0464   fChain->SetBranchAddress("t_Run", &t_Run, &b_t_Run);
0465   fChain->SetBranchAddress("t_Event", &t_Event, &b_t_Event);
0466   fChain->SetBranchAddress("t_DataType", &t_DataType, &b_t_DataType);
0467   fChain->SetBranchAddress("t_ieta", &t_ieta, &b_t_ieta);
0468   fChain->SetBranchAddress("t_iphi", &t_iphi, &b_t_iphi);
0469   fChain->SetBranchAddress("t_EventWeight", &t_EventWeight, &b_t_EventWeight);
0470   fChain->SetBranchAddress("t_nVtx", &t_nVtx, &b_t_nVtx);
0471   fChain->SetBranchAddress("t_nTrk", &t_nTrk, &b_t_nTrk);
0472   fChain->SetBranchAddress("t_goodPV", &t_goodPV, &b_t_goodPV);
0473   fChain->SetBranchAddress("t_l1pt", &t_l1pt, &b_t_l1pt);
0474   fChain->SetBranchAddress("t_l1eta", &t_l1eta, &b_t_l1eta);
0475   fChain->SetBranchAddress("t_l1phi", &t_l1phi, &b_t_l1phi);
0476   fChain->SetBranchAddress("t_l3pt", &t_l3pt, &b_t_l3pt);
0477   fChain->SetBranchAddress("t_l3eta", &t_l3eta, &b_t_l3eta);
0478   fChain->SetBranchAddress("t_l3phi", &t_l3phi, &b_t_l3phi);
0479   fChain->SetBranchAddress("t_p", &t_p, &b_t_p);
0480   fChain->SetBranchAddress("t_pt", &t_pt, &b_t_pt);
0481   fChain->SetBranchAddress("t_phi", &t_phi, &b_t_phi);
0482   fChain->SetBranchAddress("t_mindR1", &t_mindR1, &b_t_mindR1);
0483   fChain->SetBranchAddress("t_mindR2", &t_mindR2, &b_t_mindR2);
0484   fChain->SetBranchAddress("t_eMipDR", &t_eMipDR, &b_t_eMipDR);
0485   fChain->SetBranchAddress("t_eHcal", &t_eHcal, &b_t_eHcal);
0486   fChain->SetBranchAddress("t_eHcal10", &t_eHcal10, &b_t_eHcal10);
0487   fChain->SetBranchAddress("t_eHcal30", &t_eHcal30, &b_t_eHcal30);
0488   fChain->SetBranchAddress("t_hmaxNearP", &t_hmaxNearP, &b_t_hmaxNearP);
0489   fChain->SetBranchAddress("t_rhoh", &t_rhoh, &b_t_rhoh);
0490   fChain->SetBranchAddress("t_selectTk", &t_selectTk, &b_t_selectTk);
0491   fChain->SetBranchAddress("t_qltyFlag", &t_qltyFlag, &b_t_qltyFlag);
0492   fChain->SetBranchAddress("t_qltyMissFlag", &t_qltyMissFlag, &b_t_qltyMissFlag);
0493   fChain->SetBranchAddress("t_qltyPVFlag", &t_qltyPVFlag, &b_t_qltyPVFlag);
0494   fChain->SetBranchAddress("t_gentrackP", &t_gentrackP, &b_t_gentrackP);
0495   fChain->SetBranchAddress("t_DetIds", &t_DetIds, &b_t_DetIds);
0496   fChain->SetBranchAddress("t_HitEnergies", &t_HitEnergies, &b_t_HitEnergies);
0497   fChain->SetBranchAddress("t_trgbits", &t_trgbits, &b_t_trgbits);
0498   fChain->SetBranchAddress("t_DetIds1", &t_DetIds1, &b_t_DetIds1);
0499   fChain->SetBranchAddress("t_DetIds3", &t_DetIds3, &b_t_DetIds3);
0500   fChain->SetBranchAddress("t_HitEnergies1", &t_HitEnergies1, &b_t_HitEnergies1);
0501   fChain->SetBranchAddress("t_HitEnergies3", &t_HitEnergies3, &b_t_HitEnergies3);
0502   Notify();
0503 
0504   if (strcmp(dupFileName, "") != 0) {
0505     std::ifstream infil1(dupFileName);
0506     if (!infil1.is_open()) {
0507       std::cout << "Cannot open duplicate file " << dupFileName << std::endl;
0508     } else {
0509       while (1) {
0510         Long64_t jentry;
0511         infil1 >> jentry;
0512         if (!infil1.good())
0513           break;
0514         entries_.push_back(jentry);
0515       }
0516       infil1.close();
0517       std::cout << "Reads a list of " << entries_.size() << " events from " << dupFileName << std::endl;
0518     }
0519   } else {
0520     std::cout << "No duplicate events in the input file" << std::endl;
0521   }
0522 
0523   if (strcmp(comFileName, "") != 0) {
0524     std::ifstream infil2(comFileName);
0525     if (!infil2.is_open()) {
0526       std::cout << "Cannot open selection file " << comFileName << std::endl;
0527     } else {
0528       while (1) {
0529         int irun, ievt;
0530         infil2 >> irun >> ievt;
0531         if (!infil2.good())
0532           break;
0533         std::pair<int, int> key(irun, ievt);
0534         events_.push_back(key);
0535       }
0536       infil2.close();
0537       std::cout << "Reads a list of " << events_.size() << " events from " << comFileName << std::endl;
0538     }
0539   } else {
0540     std::cout << "No event list provided for selection" << std::endl;
0541   }
0542 
0543   if (((flag_ / 100) % 10) > 0) {
0544     if (((flag_ / 100) % 10) == 2) {
0545       fileout_.open(outFileName, std::ofstream::out);
0546       std::cout << "Opens " << outFileName << " in output mode" << std::endl;
0547     } else {
0548       fileout_.open(outFileName, std::ofstream::app);
0549       std::cout << "Opens " << outFileName << " in append mode" << std::endl;
0550     }
0551     fileout_ << "Input file: " << fname_ << " Directory: " << dirnm_ << " Prefix: " << prefix_ << std::endl;
0552   }
0553 
0554   double xbins[99];
0555   int nbins(-1);
0556   if (plotType_ == 0) {
0557     double xbin[9] = {-21.0, -16.0, -12.0, -6.0, 0.0, 6.0, 12.0, 16.0, 21.0};
0558     for (int i = 0; i < 9; ++i) {
0559       etas_.push_back(xbin[i]);
0560       xbins[i] = xbin[i];
0561     }
0562     nbins = 8;
0563   } else if (plotType_ == 1) {
0564     double xbin[11] = {-25.0, -20.0, -15.0, -10.0, -5.0, 0.0, 5.0, 10.0, 15.0, 20.0, 25.0};
0565     for (int i = 0; i < 11; ++i) {
0566       etas_.push_back(xbin[i]);
0567       xbins[i] = xbin[i];
0568     }
0569     nbins = 10;
0570   } else if (plotType_ == 2) {
0571     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,
0572                        3.0,   5.0,   7.0,   9.0,   11.0,  13.0,  15.0,  17.0, 19.0, 21.0, 23.0};
0573     for (int i = 0; i < 23; ++i) {
0574       etas_.push_back(xbin[i]);
0575       xbins[i] = xbin[i];
0576     }
0577     nbins = 22;
0578   } else {
0579     double xbina[99];
0580     int neta = numb_ / 2;
0581     for (int k = 0; k < neta; ++k) {
0582       xbina[k] = (k - neta) - 0.5;
0583       xbina[numb_ - k] = (neta - k) + 0.5;
0584     }
0585     xbina[neta] = 0;
0586     for (int i = 0; i < numb_ + 1; ++i) {
0587       etas_.push_back(xbina[i]);
0588       xbins[i] = xbina[i];
0589       ++nbins;
0590     }
0591   }
0592   int ipbin[npbin] = {10, 20, 30, 40, 60, 100};
0593   for (unsigned int i = 0; i < npbin; ++i)
0594     ps_.push_back((double)(ipbin[i]));
0595   int npvtx[6] = {0, 7, 10, 13, 16, 100};
0596   for (int i = 0; i < 6; ++i)
0597     nvx_.push_back(npvtx[i]);
0598   double dl1s[9] = {0, 0.10, 0.20, 0.50, 1.0, 2.0, 2.5, 3.0, 10.0};
0599   int ietasL[4] = {0, 13, 17, 0};
0600   int ietasH[4] = {14, 18, 25, 25};
0601   for (int i = 0; i < 4; ++i) {
0602     ietasL_.push_back(ietasL[i]);
0603     ietasH_.push_back(ietasH[i]);
0604   }
0605   int nxbin(100);
0606   double xlow(0.0), xhigh(5.0);
0607   if (coarseBin_ == 1) {
0608     xlow = 0.25;
0609     xhigh = 5.25;
0610     nxbin = 50;
0611   } else if (coarseBin_ > 1) {
0612     if (coarseBin_ == 2)
0613       nxbin = 500;
0614     else
0615       nxbin = 1000;
0616   }
0617 
0618   char name[20], title[200];
0619   std::string titl[5] = {
0620       "All tracks", "Good quality tracks", "Selected tracks", "Tracks with charge isolation", "Tracks MIP in ECAL"};
0621   for (int i = 0; i < 9; ++i)
0622     dl1_.push_back(dl1s[i]);
0623   unsigned int kp = (ps_.size() - 1);
0624   for (unsigned int k = 0; k < kp; ++k) {
0625     for (unsigned int j = 0; j <= ietasL_.size(); ++j) {
0626       sprintf(name, "%spp%d%d", prefix_.c_str(), k, j);
0627       if (j == 0)
0628         sprintf(title, "E/p for %s (p = %d:%d GeV All)", titl[4].c_str(), ipbin[k], ipbin[k + 1]);
0629       else
0630         sprintf(title,
0631                 "E/p for %s (p = %d:%d GeV |#eta| %d:%d)",
0632                 titl[4].c_str(),
0633                 ipbin[k],
0634                 ipbin[k + 1],
0635                 (ietasL_[j - 1] + 1),
0636                 ietasH_[j - 1]);
0637       h_pp[k].push_back(new TH1D(name, title, 100, 10.0, 110.0));
0638       int kk = h_pp[k].size() - 1;
0639       h_pp[k][kk]->Sumw2();
0640     }
0641   }
0642   if (plotType_ <= 1) {
0643     std::cout << "Book Histos for Standard" << std::endl;
0644     for (unsigned int k = 0; k < kp; ++k) {
0645       for (unsigned int i = 0; i < nvx_.size(); ++i) {
0646         sprintf(name, "%setaX%d%d", prefix_.c_str(), k, i);
0647         if (i == 0) {
0648           sprintf(title, "%s (p = %d:%d GeV all vertices)", titl[4].c_str(), ipbin[k], ipbin[k + 1]);
0649         } else {
0650           sprintf(
0651               title, "%s (p = %d:%d GeV # Vtx %d:%d)", titl[4].c_str(), ipbin[k], ipbin[k + 1], nvx_[i - 1], nvx_[i]);
0652         }
0653         h_etaX[k].push_back(new TProfile(name, title, nbins, xbins));
0654         unsigned int kk = h_etaX[k].size() - 1;
0655         h_etaX[k][kk]->Sumw2();
0656         sprintf(name, "%snvxR%d%d", prefix_.c_str(), k, i);
0657         if (i == 0) {
0658           sprintf(title, "E/p for %s (p = %d:%d GeV all vertices)", titl[4].c_str(), ipbin[k], ipbin[k + 1]);
0659         } else {
0660           sprintf(title,
0661                   "E/p for %s (p = %d:%d GeV # Vtx %d:%d)",
0662                   titl[4].c_str(),
0663                   ipbin[k],
0664                   ipbin[k + 1],
0665                   nvx_[i - 1],
0666                   nvx_[i]);
0667         }
0668         h_nvxR[k].push_back(new TH1D(name, title, nxbin, xlow, xhigh));
0669         kk = h_nvxR[k].size() - 1;
0670         h_nvxR[k][kk]->Sumw2();
0671       }
0672       for (unsigned int j = 0; j < etas_.size(); ++j) {
0673         sprintf(name, "%sratio%d%d", prefix_.c_str(), k, j);
0674         if (j == 0) {
0675           sprintf(title, "E/p for %s (p = %d:%d GeV)", titl[4].c_str(), ipbin[k], ipbin[k + 1]);
0676         } else {
0677           sprintf(title,
0678                   "E/p for %s (p = %d:%d GeV #eta %4.1f:%4.1f)",
0679                   titl[4].c_str(),
0680                   ipbin[k],
0681                   ipbin[k + 1],
0682                   etas_[j - 1],
0683                   etas_[j]);
0684         }
0685         h_etaF[k].push_back(new TH1D(name, title, nxbin, xlow, xhigh));
0686         unsigned int kk = h_etaF[k].size() - 1;
0687         h_etaF[k][kk]->Sumw2();
0688         sprintf(name, "%setaR%d%d", prefix_.c_str(), k, j);
0689         h_etaR[k].push_back(new TH1D(name, title, nxbin, xlow, xhigh));
0690         kk = h_etaR[k].size() - 1;
0691         h_etaR[k][kk]->Sumw2();
0692       }
0693       for (unsigned int j = 1; j <= ietasL_.size(); ++j) {
0694         sprintf(name, "%setaB%d%d", prefix_.c_str(), k, j);
0695         sprintf(title,
0696                 "E/p for %s (p = %d:%d GeV |#eta| %d:%d)",
0697                 titl[4].c_str(),
0698                 ipbin[k],
0699                 ipbin[k + 1],
0700                 (ietasL_[j - 1] + 1),
0701                 ietasH_[j - 1]);
0702         h_etaB[k].push_back(new TH1D(name, title, nxbin, xlow, xhigh));
0703         unsigned int kk = h_etaB[k].size() - 1;
0704         h_etaB[k][kk]->Sumw2();
0705       }
0706       for (unsigned int j = 0; j < dl1_.size(); ++j) {
0707         sprintf(name, "%sdl1R%d%d", prefix_.c_str(), k, j);
0708         if (j == 0) {
0709           sprintf(title, "E/p for %s (p = %d:%d GeV All d_{L1})", titl[4].c_str(), ipbin[k], ipbin[k + 1]);
0710         } else {
0711           sprintf(title,
0712                   "E/p for %s (p = %d:%d GeV d_{L1} %4.2f:%4.2f)",
0713                   titl[4].c_str(),
0714                   ipbin[k],
0715                   ipbin[k + 1],
0716                   dl1_[j - 1],
0717                   dl1_[j]);
0718         }
0719         h_dL1R[k].push_back(new TH1D(name, title, nxbin, xlow, xhigh));
0720         unsigned int kk = h_dL1R[k].size() - 1;
0721         h_dL1R[k][kk]->Sumw2();
0722       }
0723     }
0724     for (unsigned int i = 0; i < nvx_.size(); ++i) {
0725       sprintf(name, "%setaX%d%d", prefix_.c_str(), kp, i);
0726       if (i == 0) {
0727         sprintf(title, "%s (All Momentum all vertices)", titl[4].c_str());
0728       } else {
0729         sprintf(title, "%s (All Momentum # Vtx %d:%d)", titl[4].c_str(), nvx_[i - 1], nvx_[i]);
0730       }
0731       h_etaX[npbin - 1].push_back(new TProfile(name, title, nbins, xbins));
0732       unsigned int kk = h_etaX[npbin - 1].size() - 1;
0733       h_etaX[npbin - 1][kk]->Sumw2();
0734       sprintf(name, "%snvxR%d%d", prefix_.c_str(), kp, i);
0735       if (i == 0) {
0736         sprintf(title, "E/p for %s (All Momentum all vertices)", titl[4].c_str());
0737       } else {
0738         sprintf(title, "E/p for %s (All Momentum # Vtx %d:%d)", titl[4].c_str(), nvx_[i - 1], nvx_[i]);
0739       }
0740       h_nvxR[npbin - 1].push_back(new TH1D(name, title, 200, 0., 10.));
0741       kk = h_nvxR[npbin - 1].size() - 1;
0742       h_nvxR[npbin - 1][kk]->Sumw2();
0743     }
0744     for (unsigned int j = 0; j < etas_.size(); ++j) {
0745       sprintf(name, "%sratio%d%d", prefix_.c_str(), kp, j);
0746       if (j == 0) {
0747         sprintf(title, "E/p for %s (All momentum)", titl[4].c_str());
0748       } else {
0749         sprintf(title, "E/p for %s (All momentum #eta %4.1f:%4.1f)", titl[4].c_str(), etas_[j - 1], etas_[j]);
0750       }
0751       h_etaF[npbin - 1].push_back(new TH1D(name, title, 200, 0., 10.));
0752       unsigned int kk = h_etaF[npbin - 1].size() - 1;
0753       h_etaF[npbin - 1][kk]->Sumw2();
0754       sprintf(name, "%setaR%d%d", prefix_.c_str(), kp, j);
0755       h_etaR[npbin - 1].push_back(new TH1D(name, title, 200, 0., 10.));
0756       kk = h_etaR[npbin - 1].size() - 1;
0757       h_etaR[npbin - 1][kk]->Sumw2();
0758     }
0759     for (unsigned int j = 1; j <= ietasL_.size(); ++j) {
0760       sprintf(name, "%setaB%d%d", prefix_.c_str(), kp, j);
0761       sprintf(title, "E/p for %s (All momentum |#eta| %d:%d)", titl[4].c_str(), (ietasL_[j - 1] + 1), ietasH_[j - 1]);
0762       h_etaB[npbin - 1].push_back(new TH1D(name, title, nxbin, xlow, xhigh));
0763       unsigned int kk = h_etaB[npbin - 1].size() - 1;
0764       h_etaB[npbin - 1][kk]->Sumw2();
0765     }
0766     for (unsigned int j = 0; j < dl1_.size(); ++j) {
0767       sprintf(name, "%sdl1R%d%d", prefix_.c_str(), kp, j);
0768       if (j == 0) {
0769         sprintf(title, "E/p for %s (All momentum)", titl[4].c_str());
0770       } else {
0771         sprintf(title, "E/p for %s (All momentum d_{L1} %4.2f:%4.2f)", titl[4].c_str(), dl1_[j - 1], dl1_[j]);
0772       }
0773       h_dL1R[npbin - 1].push_back(new TH1D(name, title, 200, 0., 10.));
0774       unsigned int kk = h_dL1R[npbin - 1].size() - 1;
0775       h_dL1R[npbin - 1][kk]->Sumw2();
0776     }
0777   } else {
0778     std::cout << "Book Histos for Non-Standard " << etas_.size() << ":" << kp50 << std::endl;
0779     unsigned int kp = (ps_.size() - 1);
0780     for (unsigned int k = 0; k < kp; ++k) {
0781       for (unsigned int j = 0; j < etas_.size(); ++j) {
0782         sprintf(name, "%sratio%d%d", prefix_.c_str(), k, j);
0783         if (j == 0) {
0784           sprintf(title, "E/p for %s (p = %d:%d GeV)", titl[4].c_str(), ipbin[k], ipbin[k + 1]);
0785         } else {
0786           sprintf(title,
0787                   "E/p for %s (p = %d:%d GeV #eta %4.1f:%4.1f)",
0788                   titl[4].c_str(),
0789                   ipbin[k],
0790                   ipbin[k + 1],
0791                   etas_[j - 1],
0792                   etas_[j]);
0793         }
0794         h_etaF[k].push_back(new TH1D(name, title, nxbin, xlow, xhigh));
0795         unsigned int kk = h_etaF[k].size() - 1;
0796         h_etaF[k][kk]->Sumw2();
0797       }
0798       for (unsigned int j = 1; j <= ietasL_.size(); ++j) {
0799         sprintf(name, "%setaB%d%d", prefix_.c_str(), k, j);
0800         sprintf(title,
0801                 "E/p for %s (p = %d:%d GeV |#eta| %d:%d)",
0802                 titl[4].c_str(),
0803                 ipbin[k],
0804                 ipbin[k + 1],
0805                 (ietasL_[j - 1] + 1),
0806                 ietasH_[j - 1]);
0807         h_etaB[k].push_back(new TH1D(name, title, nxbin, xlow, xhigh));
0808         unsigned int kk = h_etaB[k].size() - 1;
0809         h_etaB[k][kk]->Sumw2();
0810       }
0811     }
0812   }
0813   if (selRBX_) {
0814     for (unsigned int j = 1; j <= 18; ++j) {
0815       sprintf(name, "%sRBX%d%d", prefix_.c_str(), kp50, j);
0816       sprintf(title, "E/p for RBX%d (p = %d:%d GeV |#eta| %d:%d)", j, ipbin[kp50], ipbin[kp50 + 1], etalo_, etahi_);
0817       h_rbx.push_back(new TH1D(name, title, nxbin, xlow, xhigh));
0818       h_rbx[j - 1]->Sumw2();
0819     }
0820   }
0821 }
0822 
0823 Bool_t CalibMonitor::Notify() {
0824   // The Notify() function is called when a new file is opened. This
0825   // can be either for a new TTree in a TChain or when when a new TTree
0826   // is started when using PROOF. It is normally not necessary to make changes
0827   // to the generated code, but the routine can be extended by the
0828   // user if needed. The return value is currently not used.
0829 
0830   return kTRUE;
0831 }
0832 
0833 void CalibMonitor::Show(Long64_t entry) {
0834   // Print contents of entry.
0835   // If entry is not specified, print current entry
0836   if (!fChain)
0837     return;
0838   fChain->Show(entry);
0839 }
0840 
0841 Int_t CalibMonitor::Cut(Long64_t) {
0842   // This function may be called from Loop.
0843   // returns  1 if entry is accepted.
0844   // returns -1 otherwise.
0845   return 1;
0846 }
0847 
0848 void CalibMonitor::Loop(Long64_t nmax) {
0849   //   In a ROOT session, you can do:
0850   //      Root > .L CalibMonitor.C
0851   //      Root > CalibMonitor t
0852   //      Root > t.GetEntry(12); // Fill t data members with entry number 12
0853   //      Root > t.Show();       // Show values of entry 12
0854   //      Root > t.Show(16);     // Read and show values of entry 16
0855   //      Root > t.Loop();       // Loop on all entries
0856   //
0857 
0858   //   This is the loop skeleton where:
0859   //      jentry is the global entry number in the chain
0860   //      ientry is the entry number in the current Tree
0861   //   Note that the argument to GetEntry must be:
0862   //      jentry for TChain::GetEntry
0863   //      ientry for TTree::GetEntry and TBranch::GetEntry
0864   //
0865   //       To read only selected branches, Insert statements like:
0866   // METHOD1:
0867   //    fChain->SetBranchStatus("*",0);  // disable all branches
0868   //    fChain->SetBranchStatus("branchname",1);  // activate branchname
0869   // METHOD2: replace line
0870   //    fChain->GetEntry(jentry);       //read all branches
0871   //by  b_branchname->GetEntry(ientry); //read only this branch
0872   const bool debug(false);
0873 
0874   if (fChain == 0)
0875     return;
0876   std::map<int, counter> runSum, runEn1, runEn2;
0877   if (debug) {
0878     std::cout << etas_.size() << " Eta Bins:";
0879     for (unsigned int j = 0; j < etas_.size(); ++j)
0880       std::cout << " " << etas_[j];
0881     std::cout << std::endl;
0882   }
0883   // Find list of duplicate events
0884   Long64_t nentries = fChain->GetEntriesFast();
0885   Long64_t entries = (nmax > 0) ? nmax : nentries;
0886   std::cout << "Total entries " << nentries << ":" << entries << std::endl;
0887   Long64_t nbytes(0), nb(0);
0888   unsigned int duplicate(0), good(0), kount(0);
0889   unsigned int kp1 = ps_.size() - 1;
0890   unsigned int kv1 = 0;
0891   std::vector<int> kounts(kp1, 0);
0892   std::vector<int> kount50(20, 0);
0893   std::vector<int> kount0(20, 0);
0894   std::vector<int> kount1(20, 0);
0895   std::vector<int> kount2(20, 0);
0896   std::vector<int> kount3(20, 0);
0897   std::vector<int> kount4(20, 0);
0898   std::vector<int> kount5(20, 0);
0899   for (Long64_t jentry = 0; jentry < entries; jentry++) {
0900     //for (Long64_t jentry=0; jentry<200;jentry++) {
0901     Long64_t ientry = LoadTree(jentry);
0902     if (ientry < 0)
0903       break;
0904     nb = fChain->GetEntry(jentry);
0905     nbytes += nb;
0906     if (jentry % 1000000 == 0)
0907       std::cout << "Entry " << jentry << " Run " << t_Run << " Event " << t_Event << std::endl;
0908     double pmom = (useGen_ && (t_gentrackP > 0)) ? t_gentrackP : t_p;
0909     int kp(-1);
0910     for (unsigned int k = 1; k < ps_.size(); ++k) {
0911       if (pmom >= ps_[k - 1] && pmom < ps_[k]) {
0912         kp = k - 1;
0913         break;
0914       }
0915     }
0916     bool p4060 = ((pmom >= 40.0) && (pmom <= 60.0));
0917     if (p4060)
0918       ++kount50[0];
0919     if (kp == 0) {
0920       ++kount0[0];
0921     } else if (kp == 1) {
0922       ++kount1[0];
0923     } else if (kp == 2) {
0924       ++kount2[0];
0925     } else if (kp == 3) {
0926       ++kount3[0];
0927     } else if (kp == 4) {
0928       ++kount4[0];
0929     } else if (kp == 5) {
0930       ++kount5[0];
0931     }
0932     bool select = (std::find(entries_.begin(), entries_.end(), jentry) == entries_.end());
0933     if (!select) {
0934       ++duplicate;
0935       if (debug)
0936         std::cout << "Duplicate event " << t_Run << " " << t_Event << " " << t_p << std::endl;
0937       continue;
0938     }
0939     if (p4060)
0940       ++kount50[1];
0941     if (kp == 0) {
0942       ++kount0[1];
0943     } else if (kp == 1) {
0944       ++kount1[1];
0945     } else if (kp == 2) {
0946       ++kount2[1];
0947     } else if (kp == 3) {
0948       ++kount3[1];
0949     } else if (kp == 4) {
0950       ++kount4[1];
0951     } else if (kp == 5) {
0952       ++kount5[1];
0953     }
0954     bool selRun = (includeRun_ ? ((t_Run >= runlo_) && (t_Run <= runhi_)) : ((t_Run < runlo_) || (t_Run > runhi_)));
0955     if (select) {
0956       if (p4060)
0957         ++kount50[2];
0958       if (kp == 0) {
0959         ++kount0[2];
0960       } else if (kp == 1) {
0961         ++kount1[2];
0962       } else if (kp == 2) {
0963         ++kount2[2];
0964       } else if (kp == 3) {
0965         ++kount3[2];
0966       } else if (kp == 4) {
0967         ++kount4[2];
0968       } else if (kp == 5) {
0969         ++kount5[2];
0970       }
0971     }
0972     select =
0973         (selRun && (fabs(t_ieta) >= etalo_) && (fabs(t_ieta) <= etahi_) && (t_nVtx >= nvxlo_) && (t_nVtx <= nvxhi_));
0974     if (select) {
0975       if (p4060)
0976         ++kount50[3];
0977       if (kp == 0) {
0978         ++kount0[3];
0979       } else if (kp == 1) {
0980         ++kount1[3];
0981       } else if (kp == 2) {
0982         ++kount2[3];
0983       } else if (kp == 3) {
0984         ++kount3[3];
0985       } else if (kp == 4) {
0986         ++kount4[3];
0987       } else if (kp == 5) {
0988         ++kount5[3];
0989       }
0990     }
0991     if (!select) {
0992       if (debug)
0993         std::cout << "Run # " << t_Run << " out of range of " << runlo_ << ":" << runhi_ << " or ieta " << t_ieta
0994                   << " (" << etalo_ << ":" << etahi_ << ") or nvtx " << t_nVtx << " (" << nvxlo_ << ":" << nvxhi_
0995                   << ") out of range" << std::endl;
0996       continue;
0997     }
0998     if (cSelect_ != nullptr) {
0999       if (exclude_) {
1000         if (cSelect_->isItRBX(t_DetIds))
1001           continue;
1002       } else {
1003         if (!(cSelect_->isItRBX(t_ieta, t_iphi)))
1004           continue;
1005       }
1006     }
1007     if (p4060)
1008       ++kount50[4];
1009     if (kp == 0) {
1010       ++kount0[4];
1011     } else if (kp == 1) {
1012       ++kount1[4];
1013     } else if (kp == 2) {
1014       ++kount2[4];
1015     } else if (kp == 3) {
1016       ++kount3[4];
1017     } else if (kp == 4) {
1018       ++kount4[4];
1019     } else if (kp == 5) {
1020       ++kount5[4];
1021     }
1022     select = (!cutL1T_ || (t_mindR1 >= 0.5));
1023     if (!select) {
1024       if (debug)
1025         std::cout << "Reject Run # " << t_Run << " Event # " << t_Event << " too close to L1 trigger " << t_mindR1
1026                   << std::endl;
1027       continue;
1028     }
1029     if (p4060)
1030       ++kount50[5];
1031     if (kp == 0) {
1032       ++kount0[5];
1033     } else if (kp == 1) {
1034       ++kount1[5];
1035     } else if (kp == 2) {
1036       ++kount2[5];
1037     } else if (kp == 3) {
1038       ++kount3[5];
1039     } else if (kp == 4) {
1040       ++kount4[5];
1041     } else if (kp == 5) {
1042       ++kount5[5];
1043     }
1044     select = ((events_.size() == 0) ||
1045               (std::find(events_.begin(), events_.end(), std::pair<int, int>(t_Run, t_Event)) != events_.end()));
1046     if (!select) {
1047       if (debug)
1048         std::cout << "Reject Run # " << t_Run << " Event # " << t_Event << " not in the selection list" << std::endl;
1049       continue;
1050     }
1051     if (p4060)
1052       ++kount50[6];
1053     if (kp == 0) {
1054       ++kount0[6];
1055     } else if (kp == 1) {
1056       ++kount1[6];
1057     } else if (kp == 2) {
1058       ++kount2[6];
1059     } else if (kp == 3) {
1060       ++kount3[6];
1061     } else if (kp == 4) {
1062       ++kount4[6];
1063     } else if (kp == 5) {
1064       ++kount5[6];
1065     }
1066     if ((ifDepth_ == 3) && (cFactor_ != nullptr) && (cFactor_->absent(ientry)))
1067       continue;
1068 
1069     // if (Cut(ientry) < 0) continue;
1070     int jp(-1), jp1(-1), jp2(-1);
1071     unsigned int kv = nvx_.size() - 1;
1072     for (unsigned int k = 1; k < nvx_.size(); ++k) {
1073       if (t_goodPV >= nvx_[k - 1] && t_goodPV < nvx_[k]) {
1074         kv = k;
1075         break;
1076       }
1077     }
1078     unsigned int kd1 = 0;
1079     unsigned int kd = dl1_.size() - 1;
1080     for (unsigned int k = 1; k < dl1_.size(); ++k) {
1081       if (t_mindR1 >= dl1_[k - 1] && t_mindR1 < dl1_[k]) {
1082         kd = k;
1083         break;
1084       }
1085     }
1086     double eta = (t_ieta > 0) ? ((double)(t_ieta)-0.001) : ((double)(t_ieta) + 0.001);
1087     for (unsigned int j = 1; j < etas_.size(); ++j) {
1088       if (eta > etas_[j - 1] && eta < etas_[j]) {
1089         jp = j;
1090         break;
1091       }
1092     }
1093     for (unsigned int j = 0; j < (ietasL_.size() - 1); ++j) {
1094       if (std::abs(t_ieta) > ietasL_[j] && std::abs(t_ieta) <= ietasH_[j]) {
1095         if (jp1 >= 0)
1096           jp2 = j;
1097         if (jp1 < 0)
1098           jp1 = j;
1099       }
1100     }
1101     int jp3 = (jp1 >= 0) ? (int)(ietasL_.size() - 1) : jp1;
1102     double cut = (pmom > 20) ? ((flexibleSelect_ == 0) ? 2.0 : 10.0) : 0.0;
1103     //  double rcut= (pmom > 20) ? 0.25: 0.1;
1104     double rcut(-1000.0);
1105     if (debug)
1106       std::cout << "Bin " << kp << ":" << kp1 << ":" << kv << ":" << kv1 << ":" << kd << ":" << kd1 << ":" << jp << ":"
1107                 << jp1 << ":" << jp2 << ":" << jp3 << " pmom:ieta:pv:mindR " << pmom << ":" << std::abs(t_ieta) << ":"
1108                 << t_ieta << ":" << t_goodPV << ":" << t_mindR1 << " Cuts " << cut << ":" << rcut << std::endl;
1109 
1110     // Selection of good track and energy measured in Hcal
1111     double rat(1.0), eHcal(t_eHcal);
1112     if (corrFactor_->doCorr() || (cFactor_ != nullptr)) {
1113       eHcal = 0;
1114       for (unsigned int k = 0; k < t_HitEnergies->size(); ++k) {
1115         // The masks are defined in DataFormats/HcalDetId/interface/HcalDetId.h
1116         double cfac(1.0);
1117         if (corrFactor_->doCorr()) {
1118           unsigned int id = truncateId((*t_DetIds)[k], truncateFlag_, false);
1119           cfac = corrFactor_->getCorr(id);
1120         }
1121         if ((cFactor_ != nullptr) && (ifDepth_ != 3))
1122           cfac *= cFactor_->getCorr(t_Run, (*t_DetIds)[k]);
1123         eHcal += (cfac * ((*t_HitEnergies)[k]));
1124         if (debug) {
1125           int subdet, zside, ieta, iphi, depth;
1126           unpackDetId((*t_DetIds)[k], subdet, zside, ieta, iphi, depth);
1127           std::cout << zside << ":" << ieta << ":" << depth << " Corr " << cfac << " " << (*t_HitEnergies)[k] << " Out "
1128                     << eHcal << std::endl;
1129         }
1130       }
1131     }
1132     bool goodTk = goodTrack(eHcal, cut, jentry, debug);
1133     bool selPhi = selectPhi(debug);
1134     if (t_qltyFlag) {
1135       if (p4060)
1136         ++kount50[7];
1137       if (kp == 0) {
1138         ++kount0[7];
1139       } else if (kp == 1) {
1140         ++kount1[7];
1141       } else if (kp == 2) {
1142         ++kount2[7];
1143       } else if (kp == 3) {
1144         ++kount3[7];
1145       } else if (kp == 4) {
1146         ++kount4[7];
1147       } else if (kp == 5) {
1148         ++kount5[7];
1149       }
1150       if (t_selectTk) {
1151         if (p4060)
1152           ++kount50[8];
1153         if (kp == 0) {
1154           ++kount0[8];
1155         } else if (kp == 1) {
1156           ++kount1[8];
1157         } else if (kp == 2) {
1158           ++kount2[8];
1159         } else if (kp == 3) {
1160           ++kount3[8];
1161         } else if (kp == 4) {
1162           ++kount4[8];
1163         } else if (kp == 5) {
1164           ++kount5[8];
1165         }
1166         if (t_hmaxNearP < cut) {
1167           if (p4060)
1168             ++kount50[9];
1169           if (kp == 0) {
1170             ++kount0[9];
1171           } else if (kp == 1) {
1172             ++kount1[9];
1173           } else if (kp == 2) {
1174             ++kount2[9];
1175           } else if (kp == 3) {
1176             ++kount3[9];
1177           } else if (kp == 4) {
1178             ++kount4[9];
1179           } else if (kp == 5) {
1180             ++kount5[9];
1181           }
1182           if (t_eMipDR < 1.0) {
1183             if (p4060)
1184               ++kount50[10];
1185             if (kp == 0) {
1186               ++kount0[10];
1187             } else if (kp == 1) {
1188               ++kount1[10];
1189             } else if (kp == 2) {
1190               ++kount2[10];
1191             } else if (kp == 3) {
1192               ++kount3[10];
1193             } else if (kp == 4) {
1194               ++kount4[10];
1195             } else if (kp == 5) {
1196               ++kount5[10];
1197             }
1198             if (eHcal > 0.001) {
1199               if (p4060)
1200                 ++kount50[11];
1201               if (kp == 0) {
1202                 ++kount0[11];
1203               } else if (kp == 1) {
1204                 ++kount1[11];
1205               } else if (kp == 2) {
1206                 ++kount2[11];
1207               } else if (kp == 3) {
1208                 ++kount3[11];
1209               } else if (kp == 4) {
1210                 ++kount4[11];
1211               } else if (kp == 5) {
1212                 ++kount5[11];
1213               }
1214               if (selPhi) {
1215                 if (p4060)
1216                   ++kount50[12];
1217                 if (kp == 0) {
1218                   ++kount0[12];
1219                 } else if (kp == 1) {
1220                   ++kount1[12];
1221                 } else if (kp == 2) {
1222                   ++kount2[12];
1223                 } else if (kp == 3) {
1224                   ++kount3[12];
1225                 } else if (kp == 4) {
1226                   ++kount4[12];
1227                 } else if (kp == 5) {
1228                   ++kount5[12];
1229                 }
1230               }
1231             }
1232           }
1233         }
1234       }
1235     }
1236     if (pmom > 0)
1237       rat = (eHcal / (pmom - t_eMipDR));
1238     if (debug) {
1239       std::cout << "Entry " << jentry << " p|eHcal|ratio " << pmom << "|" << t_eHcal << "|" << eHcal << "|" << rat
1240                 << "|" << kp << "|" << kv << "|" << jp << " Cuts " << t_qltyFlag << "|" << t_selectTk << "|"
1241                 << (t_hmaxNearP < cut) << "|" << (t_eMipDR < 1.0) << "|" << goodTk << "|" << (rat > rcut)
1242                 << " Select Phi " << selPhi << " hmaxNearP " << t_hmaxNearP << " eMipDR " << t_eMipDR << std::endl;
1243       std::cout << "D1 : " << kp << ":" << kp1 << ":" << kv << ":" << kv1 << ":" << kd << ":" << kd1 << ":" << jp
1244                 << std::endl;
1245     }
1246     if (goodTk && (kp >= 0) && selPhi) {
1247       if (p4060)
1248         ++kount50[13];
1249       if (kp == 0) {
1250         ++kount0[13];
1251       } else if (kp == 1) {
1252         ++kount1[13];
1253       } else if (kp == 2) {
1254         ++kount2[13];
1255       } else if (kp == 3) {
1256         ++kount3[13];
1257       } else if (kp == 4) {
1258         ++kount4[13];
1259       } else if (kp == 5) {
1260         ++kount5[13];
1261       }
1262       if (t_eHcal < 0.01) {
1263         std::map<int, counter>::const_iterator itr = runEn1.find(t_Run);
1264         if (itr == runEn1.end()) {
1265           counter knt;
1266           if ((kp >= 0) && (kp < counter::npsize))
1267             knt.count[kp] = 1;
1268           knt.total = 1;
1269           runEn1[t_Run] = knt;
1270         } else {
1271           counter knt = runEn1[t_Run];
1272           if ((kp >= 0) && (kp < counter::npsize))
1273             ++knt.count[kp];
1274           ++knt.total;
1275           runEn1[t_Run] = knt;
1276         }
1277       }
1278       if (t_eMipDR < 0.01 && t_eHcal < 0.01) {
1279         if (p4060)
1280           ++kount50[14];
1281         if (kp == 0) {
1282           ++kount0[14];
1283         } else if (kp == 1) {
1284           ++kount1[14];
1285         } else if (kp == 2) {
1286           ++kount2[14];
1287         } else if (kp == 3) {
1288           ++kount3[14];
1289         } else if (kp == 4) {
1290           ++kount4[14];
1291         } else if (kp == 5) {
1292           ++kount5[14];
1293         }
1294         std::map<int, counter>::const_iterator itr = runEn2.find(t_Run);
1295         if (itr == runEn2.end()) {
1296           counter knt;
1297           if ((kp >= 0) && (kp < counter::npsize))
1298             knt.count[kp] = 1;
1299           knt.total = 1;
1300           runEn2[t_Run] = knt;
1301         } else {
1302           counter knt = runEn2[t_Run];
1303           if ((kp >= 0) && (kp < counter::npsize))
1304             ++knt.count[kp];
1305           ++knt.total;
1306           runEn2[t_Run] = knt;
1307         }
1308       }
1309       if (rat > rcut) {
1310         if (p4060)
1311           ++kount50[15];
1312         if (kp == 0) {
1313           ++kount0[15];
1314         } else if (kp == 1) {
1315           ++kount1[15];
1316         } else if (kp == 2) {
1317           ++kount2[15];
1318         } else if (kp == 3) {
1319           ++kount3[15];
1320         } else if (kp == 4) {
1321           ++kount4[15];
1322         } else if (kp == 5) {
1323           ++kount5[15];
1324         }
1325         if (plotType_ <= 1) {
1326           h_etaX[kp][kv]->Fill(eta, rat, t_EventWeight);
1327           h_etaX[kp][kv1]->Fill(eta, rat, t_EventWeight);
1328           h_nvxR[kp][kv]->Fill(rat, t_EventWeight);
1329           h_nvxR[kp][kv1]->Fill(rat, t_EventWeight);
1330           h_dL1R[kp][kd]->Fill(rat, t_EventWeight);
1331           h_dL1R[kp][kd1]->Fill(rat, t_EventWeight);
1332           if (jp > 0)
1333             h_etaR[kp][jp]->Fill(rat, t_EventWeight);
1334           h_etaR[kp][0]->Fill(rat, t_EventWeight);
1335         }
1336         h_pp[kp][0]->Fill(pmom, t_EventWeight);
1337         if (jp1 >= 0) {
1338           h_pp[kp][jp1 + 1]->Fill(pmom, t_EventWeight);
1339           h_pp[kp][jp3 + 1]->Fill(pmom, t_EventWeight);
1340         }
1341         if (jp2 >= 0)
1342           h_pp[kp][jp2 + 1]->Fill(pmom, t_EventWeight);
1343         if (kp == (int)(kp50)) {
1344           std::map<int, counter>::const_iterator itr = runSum.find(t_Run);
1345           if (itr == runSum.end()) {
1346             counter knt;
1347             if ((kp >= 0) && (kp < counter::npsize))
1348               knt.count[kp] = 1;
1349             knt.total = 1;
1350             runSum[t_Run] = knt;
1351           } else {
1352             counter knt = runSum[t_Run];
1353             if ((kp >= 0) && (kp < counter::npsize))
1354               ++knt.count[kp];
1355             ++knt.total;
1356             runSum[t_Run] = knt;
1357           }
1358         }
1359         if ((!dataMC_) || (t_mindR1 > 0.5) || (t_DataType == 1)) {
1360           if (p4060)
1361             ++kount50[16];
1362           if (kp == 0) {
1363             ++kount0[16];
1364           } else if (kp == 1) {
1365             ++kount1[16];
1366           } else if (kp == 2) {
1367             ++kount2[16];
1368           } else if (kp == 3) {
1369             ++kount3[16];
1370           } else if (kp == 4) {
1371             ++kount4[16];
1372           } else if (kp == 5) {
1373             ++kount5[16];
1374           }
1375           ++kounts[kp];
1376           if (plotType_ <= 1) {
1377             if (jp > 0)
1378               h_etaF[kp][jp]->Fill(rat, t_EventWeight);
1379             h_etaF[kp][0]->Fill(rat, t_EventWeight);
1380           } else {
1381             if (debug) {
1382               std::cout << "kp " << kp << h_etaF[kp].size() << std::endl;
1383             }
1384             if (jp > 0)
1385               h_etaF[kp][jp]->Fill(rat, t_EventWeight);
1386             h_etaF[kp][0]->Fill(rat, t_EventWeight);
1387             if (jp1 >= 0) {
1388               h_etaB[kp][jp1]->Fill(rat, t_EventWeight);
1389               h_etaB[kp][jp3]->Fill(rat, t_EventWeight);
1390             }
1391             if (jp2 >= 0)
1392               h_etaB[kp][jp2]->Fill(rat, t_EventWeight);
1393           }
1394           if (selRBX_ && (kp == (int)(kp50)) && ((t_ieta * zside_) > 0)) {
1395             int rbx = (t_iphi > 70) ? 0 : (t_iphi + 1) / 4;
1396             h_rbx[rbx]->Fill(rat, t_EventWeight);
1397           }
1398         }
1399         if (pmom > 10.0) {
1400           if (plotType_ <= 1) {
1401             h_etaX[kp1][kv]->Fill(eta, rat, t_EventWeight);
1402             h_etaX[kp1][kv1]->Fill(eta, rat, t_EventWeight);
1403             h_nvxR[kp1][kv]->Fill(rat, t_EventWeight);
1404             h_nvxR[kp1][kv1]->Fill(rat, t_EventWeight);
1405             h_dL1R[kp1][kd]->Fill(rat, t_EventWeight);
1406             h_dL1R[kp1][kd1]->Fill(rat, t_EventWeight);
1407             if (jp > 0)
1408               h_etaR[kp1][jp]->Fill(rat, t_EventWeight);
1409             h_etaR[kp1][0]->Fill(rat, t_EventWeight);
1410             if (jp1 >= 0) {
1411               h_etaB[kp][jp1]->Fill(rat, t_EventWeight);
1412               h_etaB[kp][jp3]->Fill(rat, t_EventWeight);
1413             }
1414             if (jp2 >= 0)
1415               h_etaB[kp][jp2]->Fill(rat, t_EventWeight);
1416           }
1417           if (p4060)
1418             ++kount50[17];
1419           if (kp == 0) {
1420             ++kount0[17];
1421           } else if (kp == 1) {
1422             ++kount1[17];
1423           } else if (kp == 2) {
1424             ++kount2[17];
1425           } else if (kp == 3) {
1426             ++kount3[17];
1427           } else if (kp == 4) {
1428             ++kount4[17];
1429           } else if (kp == 5) {
1430             ++kount5[17];
1431           }
1432         }
1433       }
1434     }
1435     if (pmom > 10.0) {
1436       kount++;
1437       if (((flag_ / 100) % 10) != 0) {
1438         good++;
1439         fileout_ << good << " " << jentry << " " << t_Run << " " << t_Event << " " << t_ieta << " " << pmom
1440                  << std::endl;
1441       }
1442     }
1443   }
1444   unsigned int k(0);
1445   std::cout << "\nSummary of entries with " << runSum.size() << " runs\n";
1446   for (std::map<int, counter>::iterator itr = runSum.begin(); itr != runSum.end(); ++itr, ++k)
1447     std::cout << "[" << k << "] Run " << itr->first << " Total " << (itr->second).total << " in p-bins "
1448               << (itr->second).count[0] << ":" << (itr->second).count[1] << ":" << (itr->second).count[2] << ":"
1449               << (itr->second).count[3] << std::endl;
1450   k = 0;
1451   std::cout << runEn1.size() << " runs with 0 energy in HCAL\n";
1452   for (std::map<int, counter>::iterator itr = runEn1.begin(); itr != runEn1.end(); ++itr, ++k)
1453     std::cout << "[" << k << "] Run " << itr->first << " Total " << (itr->second).total << " in p-bins "
1454               << (itr->second).count[0] << ":" << (itr->second).count[1] << ":" << (itr->second).count[2] << ":"
1455               << (itr->second).count[3] << std::endl;
1456   k = 0;
1457   std::cout << runEn2.size() << " runs with 0 energy in ECAL and HCAL\n";
1458   for (std::map<int, counter>::iterator itr = runEn2.begin(); itr != runEn2.end(); ++itr, ++k)
1459     std::cout << "[" << k << "] Run " << itr->first << " Total " << (itr->second).total << " in p-bins "
1460               << (itr->second).count[0] << ":" << (itr->second).count[1] << ":" << (itr->second).count[2] << ":"
1461               << (itr->second).count[3] << std::endl;
1462   if (((flag_ / 100) % 10) > 0) {
1463     fileout_.close();
1464     std::cout << "Writes " << good << " events in the file " << outFileName_ << std::endl;
1465   }
1466   std::cout << "Finds " << duplicate << " Duplicate events out of " << kount << " events in this file with p>10 Gev"
1467             << std::endl;
1468   std::cout << "Number of selected events:" << std::endl;
1469   for (unsigned int k = 1; k < ps_.size(); ++k)
1470     std::cout << ps_[k - 1] << ":" << ps_[k] << "     " << kounts[k - 1] << std::endl;
1471   std::cout << "Number in each step for tracks of momentum 40-60 GeV: ";
1472   for (unsigned int k = 0; k < 18; ++k)
1473     std::cout << " [" << k << "] " << kount50[k];
1474   std::cout << std::endl;
1475   for (unsigned int k = 1; k < ps_.size(); ++k) {
1476     std::cout << "Number in each step for tracks of momentum " << ps_[k - 1] << "-" << ps_[k] << " Gev: ";
1477     for (unsigned int k1 = 0; k1 < 18; ++k1) {
1478       if (k == 1) {
1479         std::cout << " [" << k1 << "] " << kount0[k1];
1480       } else if (k == 2) {
1481         std::cout << " [" << k1 << "] " << kount1[k1];
1482       } else if (k == 3) {
1483         std::cout << " [" << k1 << "] " << kount2[k1];
1484       } else if (k == 4) {
1485         std::cout << " [" << k1 << "] " << kount3[k1];
1486       } else if (k == 5) {
1487         std::cout << " [" << k1 << "] " << kount4[k1];
1488       } else if (k == 6) {
1489         std::cout << " [" << k1 << "] " << kount5[k1];
1490       }
1491     }
1492     std::cout << std::endl;
1493   }
1494 }
1495 
1496 bool CalibMonitor::goodTrack(double &eHcal, double &cuti, const Long64_t &entry, bool debug) {
1497   bool select(true);
1498   double cut(cuti);
1499   if (debug) {
1500     std::cout << "goodTrack input " << eHcal << ":" << cut;
1501   }
1502   if (flexibleSelect_ > 1) {
1503     double eta = (t_ieta > 0) ? t_ieta : -t_ieta;
1504     cut = 8.0 * exp(eta * log2by18_);
1505   }
1506   correctEnergy(eHcal, entry);
1507   select = ((t_qltyFlag) && (t_selectTk) && (t_hmaxNearP < cut) && (t_eMipDR < 1.0) && (eHcal > 0.001));
1508   if (debug) {
1509     std::cout << " output " << select << " Based on " << t_qltyFlag << ":" << t_selectTk << ":" << t_hmaxNearP << ":"
1510               << cut << ":" << t_eMipDR << ":" << eHcal << std::endl;
1511   }
1512   return select;
1513 }
1514 
1515 bool CalibMonitor::selectPhi(bool debug) {
1516   bool select(true);
1517   if (phimin_ > 1 || phimax_ < 72) {
1518     double eTotal(0), eSelec(0);
1519     // The masks are defined in DataFormats/HcalDetId/interface/HcalDetId.h
1520     for (unsigned int k = 0; k < t_HitEnergies->size(); ++k) {
1521       int iphi = ((*t_DetIds)[k]) & (0x3FF);
1522       int zside = ((*t_DetIds)[k] & 0x80000) ? (1) : (-1);
1523       eTotal += ((*t_HitEnergies)[k]);
1524       if (iphi >= phimin_ && iphi <= phimax_ && zside == zside_)
1525         eSelec += ((*t_HitEnergies)[k]);
1526     }
1527     if (eSelec < 0.9 * eTotal)
1528       select = false;
1529     if (debug) {
1530       std::cout << "Etotal " << eTotal << " and ESelec " << eSelec << " (phi " << phimin_ << ":" << phimax_ << " z "
1531                 << zside_ << ") Selection " << select << std::endl;
1532     }
1533   }
1534   return select;
1535 }
1536 
1537 void CalibMonitor::plotHist(int itype, int inum, bool save) {
1538   gStyle->SetCanvasBorderMode(0);
1539   gStyle->SetCanvasColor(kWhite);
1540   gStyle->SetPadColor(kWhite);
1541   gStyle->SetFillColor(kWhite);
1542   gStyle->SetOptStat(111110);
1543   gStyle->SetOptFit(1);
1544   char name[100];
1545   int itmin = (itype >= 0 && itype < 4) ? itype : 0;
1546   int itmax = (itype >= 0 && itype < 4) ? itype : 3;
1547   std::string types[4] = {
1548       "E_{HCAL}/(p-E_{ECAL})", "E_{HCAL}/(p-E_{ECAL})", "E_{HCAL}/(p-E_{ECAL})", "E_{HCAL}/(p-E_{ECAL})"};
1549   int nmax[4] = {npbin, npbin, npbin, npbin};
1550   for (int type = itmin; type <= itmax; ++type) {
1551     int inmin = (inum >= 0 && inum < nmax[type]) ? inum : 0;
1552     int inmax = (inum >= 0 && inum < nmax[type]) ? inum : nmax[type] - 1;
1553     int kmax = 1;
1554     if (type == 0)
1555       kmax = (int)(etas_.size());
1556     else if (type == 1)
1557       kmax = (int)(etas_.size());
1558     else if (type == 2)
1559       kmax = (int)(nvx_.size());
1560     else
1561       kmax = (int)(dl1_.size());
1562     for (int num = inmin; num <= inmax; ++num) {
1563       for (int k = 0; k < kmax; ++k) {
1564         sprintf(name, "c_%d%d%d", type, num, k);
1565         TCanvas *pad = new TCanvas(name, name, 700, 500);
1566         pad->SetRightMargin(0.10);
1567         pad->SetTopMargin(0.10);
1568         sprintf(name, "%s", types[type].c_str());
1569         if (type != 7) {
1570           TH1D *hist(0);
1571           if (type == 0)
1572             hist = (TH1D *)(h_etaR[num][k]->Clone());
1573           else if (type == 1)
1574             hist = (TH1D *)(h_etaF[num][k]->Clone());
1575           else if (type == 2)
1576             hist = (TH1D *)(h_nvxR[num][k]->Clone());
1577           else
1578             hist = (TH1D *)(h_dL1R[num][k]->Clone());
1579           hist->GetXaxis()->SetTitle(name);
1580           hist->GetYaxis()->SetTitle("Tracks");
1581           drawHist(hist, pad);
1582           if (save) {
1583             sprintf(name, "c_%s%d%d%d.gif", prefix_.c_str(), type, num, k);
1584             pad->Print(name);
1585           }
1586         } else {
1587           TProfile *hist = (TProfile *)(h_etaX[num][k]->Clone());
1588           hist->GetXaxis()->SetTitle(name);
1589           hist->GetYaxis()->SetTitle("<E_{HCAL}/(p-E_{ECAL})>");
1590           hist->GetYaxis()->SetRangeUser(0.4, 1.6);
1591           hist->Fit("pol0", "q");
1592           drawHist(hist, pad);
1593           if (save) {
1594             sprintf(name, "c_%s%d%d%d.gif", prefix_.c_str(), type, num, k);
1595             pad->Print(name);
1596           }
1597         }
1598       }
1599     }
1600   }
1601 }
1602 
1603 template <class Hist>
1604 void CalibMonitor::drawHist(Hist *hist, TCanvas *pad) {
1605   hist->GetYaxis()->SetLabelOffset(0.005);
1606   hist->GetYaxis()->SetTitleOffset(1.20);
1607   hist->Draw();
1608   pad->Update();
1609   TPaveStats *st1 = (TPaveStats *)hist->GetListOfFunctions()->FindObject("stats");
1610   if (st1 != NULL) {
1611     st1->SetY1NDC(0.70);
1612     st1->SetY2NDC(0.90);
1613     st1->SetX1NDC(0.55);
1614     st1->SetX2NDC(0.90);
1615   }
1616   pad->Modified();
1617   pad->Update();
1618 }
1619 
1620 void CalibMonitor::savePlot(const std::string &theName, bool append, bool all) {
1621   TFile *theFile(0);
1622   if (append) {
1623     theFile = new TFile(theName.c_str(), "UPDATE");
1624   } else {
1625     theFile = new TFile(theName.c_str(), "RECREATE");
1626   }
1627 
1628   theFile->cd();
1629   for (unsigned int k = 0; k < ps_.size(); ++k) {
1630     for (unsigned int j = 0; j <= ietasL_.size(); ++j) {
1631       if ((all || k == kp50) && h_pp[k].size() > j && h_pp[k][j] != 0) {
1632         TH1D *hist = (TH1D *)h_pp[k][j]->Clone();
1633         hist->Write();
1634       }
1635     }
1636     if (plotType_ <= 1) {
1637       for (unsigned int i = 0; i < nvx_.size(); ++i) {
1638         if (h_etaX[k][i] != 0) {
1639           TProfile *hnew = (TProfile *)h_etaX[k][i]->Clone();
1640           hnew->Write();
1641         }
1642         if (h_nvxR[k].size() > i && h_nvxR[k][i] != 0) {
1643           TH1D *hist = (TH1D *)h_nvxR[k][i]->Clone();
1644           hist->Write();
1645         }
1646       }
1647     }
1648     for (unsigned int j = 0; j < etas_.size(); ++j) {
1649       if ((plotType_ <= 1) && (h_etaR[k][j] != 0)) {
1650         TH1D *hist = (TH1D *)h_etaR[k][j]->Clone();
1651         hist->Write();
1652       }
1653       if (h_etaF[k].size() > j && h_etaF[k][j] != 0 && (all || (k == kp50))) {
1654         TH1D *hist = (TH1D *)h_etaF[k][j]->Clone();
1655         hist->Write();
1656       }
1657     }
1658     if (plotType_ <= 1) {
1659       for (unsigned int j = 0; j < dl1_.size(); ++j) {
1660         if (h_dL1R[k][j] != 0) {
1661           TH1D *hist = (TH1D *)h_dL1R[k][j]->Clone();
1662           hist->Write();
1663         }
1664       }
1665     }
1666     for (unsigned int j = 0; j < ietasL_.size(); ++j) {
1667       if (h_etaB[k].size() > j && h_etaB[k][j] != 0 && (all || (k == kp50))) {
1668         TH1D *hist = (TH1D *)h_etaB[k][j]->Clone();
1669         hist->Write();
1670       }
1671     }
1672   }
1673   if (selRBX_) {
1674     for (unsigned int k = 0; k < 18; ++k) {
1675       if (h_rbx[k] != 0) {
1676         TH1D *h1 = (TH1D *)h_rbx[k]->Clone();
1677         h1->Write();
1678       }
1679     }
1680   }
1681   std::cout << "All done" << std::endl;
1682   theFile->Close();
1683 }
1684 
1685 void CalibMonitor::correctEnergy(double &eHcal, const Long64_t &entry) {
1686   bool debug(false);
1687   double pmom = (useGen_ && (t_gentrackP > 0)) ? t_gentrackP : t_p;
1688   if ((ifDepth_ == 3) && (cFactor_ != nullptr)) {
1689     double cfac = cFactor_->getCorr(entry);
1690     eHcal *= cfac;
1691     if (debug)
1692       std::cout << "PU Factor for " << ifDepth_ << ":" << entry << " = " << cfac << ":" << eHcal << std::endl;
1693   } else if ((corrPU_ < 0) && (pmom > 0)) {
1694     double ediff = (t_eHcal30 - t_eHcal10);
1695     if (t_DetIds1 != 0 && t_DetIds3 != 0) {
1696       double Etot1(0), Etot3(0);
1697       // The masks are defined in DataFormats/HcalDetId/interface/HcalDetId.h
1698       for (unsigned int idet = 0; idet < (*t_DetIds1).size(); idet++) {
1699         unsigned int id = truncateId((*t_DetIds1)[idet], truncateFlag_, false);
1700         double cfac = corrFactor_->getCorr(id);
1701         if (cFactor_ != 0)
1702           cfac *= cFactor_->getCorr(t_Run, (*t_DetIds1)[idet]);
1703         double hitEn = cfac * (*t_HitEnergies1)[idet];
1704         Etot1 += hitEn;
1705       }
1706       for (unsigned int idet = 0; idet < (*t_DetIds3).size(); idet++) {
1707         unsigned int id = truncateId((*t_DetIds3)[idet], truncateFlag_, false);
1708         double cfac = corrFactor_->getCorr(id);
1709         if (cFactor_ != 0)
1710           cfac *= cFactor_->getCorr(t_Run, (*t_DetIds3)[idet]);
1711         double hitEn = cfac * (*t_HitEnergies3)[idet];
1712         Etot3 += hitEn;
1713       }
1714       ediff = (Etot3 - Etot1);
1715     }
1716     double fac = puFactor(-corrPU_, t_ieta, pmom, eHcal, ediff, false);
1717     if (debug) {
1718       double fac1 = puFactor(-corrPU_, t_ieta, pmom, eHcal, ediff, true);
1719       double fac2 = puFactor(2, t_ieta, pmom, eHcal, ediff, true);
1720       std::cout << "PU Factor for " << -corrPU_ << " = " << fac1 << "; for 2 = " << fac2 << std::endl;
1721     }
1722     eHcal *= fac;
1723   } else if (corrPU_ > 0) {
1724     eHcal = puFactorRho(corrPU_, t_ieta, t_rhoh, eHcal);
1725   }
1726 }
1727 
1728 class GetEntries {
1729 public:
1730   TTree *fChain;   //!pointer to the analyzed TTree/TChain
1731   Int_t fCurrent;  //!current Tree number in a TChain
1732 
1733   // Declaration of leaf types
1734   UInt_t t_RunNo;
1735   UInt_t t_EventNo;
1736   Int_t t_Tracks;
1737   Int_t t_TracksProp;
1738   Int_t t_TracksSaved;
1739   Int_t t_TracksLoose;
1740   Int_t t_TracksTight;
1741   Int_t t_allvertex;
1742   Bool_t t_TrigPass;
1743   Bool_t t_TrigPassSel;
1744   Bool_t t_L1Bit;
1745   std::vector<Bool_t> *t_hltbits;
1746   std::vector<int> *t_ietaAll;
1747   std::vector<int> *t_ietaGood;
1748   std::vector<int> *t_trackType;
1749 
1750   // List of branches
1751   TBranch *b_t_RunNo;        //!
1752   TBranch *b_t_EventNo;      //!
1753   TBranch *b_t_Tracks;       //!
1754   TBranch *b_t_TracksProp;   //!
1755   TBranch *b_t_TracksSaved;  //!
1756   TBranch *b_t_TracksLoose;  //!
1757   TBranch *b_t_TracksTight;  //!
1758   TBranch *b_t_allvertex;    //!
1759   TBranch *b_t_TrigPass;     //!
1760   TBranch *b_t_TrigPassSel;  //!
1761   TBranch *b_t_L1Bit;        //!
1762   TBranch *b_t_hltbits;      //!
1763   TBranch *b_t_ietaAll;      //!
1764   TBranch *b_t_ietaGood;     //!
1765   TBranch *b_t_trackType;    //!
1766 
1767   GetEntries(const char *fname,
1768              const char *dirname,
1769              const char *dupFileName,
1770              const unsigned int bit1,
1771              const unsigned int bit2);
1772   virtual ~GetEntries();
1773   virtual Int_t Cut(Long64_t entry);
1774   virtual Int_t GetEntry(Long64_t entry);
1775   virtual Long64_t LoadTree(Long64_t entry);
1776   virtual void Init(TTree *tree, const char *dupFileName);
1777   virtual void Loop(Long64_t nmax = -1);
1778   virtual Bool_t Notify();
1779   virtual void Show(Long64_t entry = -1);
1780 
1781 private:
1782   unsigned int bit_[2];
1783   std::vector<Long64_t> entries_;
1784   TH1I *h_tk[3], *h_eta[4], *h_pvx[3];
1785   TH1D *h_eff[3];
1786 };
1787 
1788 GetEntries::GetEntries(
1789     const char *fname, const char *dirnm, const char *dupFileName, const unsigned int bit1, const unsigned int bit2) {
1790   TFile *file = new TFile(fname);
1791   TDirectory *dir = (TDirectory *)file->FindObjectAny(dirnm);
1792   std::cout << fname << " file " << file << " " << dirnm << " " << dir << std::endl;
1793   TTree *tree = (TTree *)dir->Get("EventInfo");
1794   std::cout << "CalibTree " << tree << std::endl;
1795   bit_[0] = bit1;
1796   bit_[1] = bit2;
1797   Init(tree, dupFileName);
1798 }
1799 
1800 GetEntries::~GetEntries() {
1801   if (!fChain)
1802     return;
1803   delete fChain->GetCurrentFile();
1804 }
1805 
1806 Int_t GetEntries::GetEntry(Long64_t entry) {
1807   // Read contents of entry.
1808   if (!fChain)
1809     return 0;
1810   return fChain->GetEntry(entry);
1811 }
1812 
1813 Long64_t GetEntries::LoadTree(Long64_t entry) {
1814   // Set the environment to read one entry
1815   if (!fChain)
1816     return -5;
1817   Long64_t centry = fChain->LoadTree(entry);
1818   if (centry < 0)
1819     return centry;
1820   if (!fChain->InheritsFrom(TChain::Class()))
1821     return centry;
1822   TChain *chain = (TChain *)fChain;
1823   if (chain->GetTreeNumber() != fCurrent) {
1824     fCurrent = chain->GetTreeNumber();
1825     Notify();
1826   }
1827   return centry;
1828 }
1829 
1830 void GetEntries::Init(TTree *tree, const char *dupFileName) {
1831   // The Init() function is called when the selector needs to initialize
1832   // a new tree or chain. Typically here the branch addresses and branch
1833   // pointers of the tree will be set.
1834   // It is normally not necessary to make changes to the generated
1835   // code, but the routine can be extended by the user if needed.
1836   // Init() will be called many times when running on PROOF
1837   // (once per file to be processed).
1838 
1839   // Set branch addresses and branch pointers
1840   // Set object pointer
1841   t_hltbits = 0;
1842   t_ietaAll = 0;
1843   t_ietaGood = 0;
1844   t_trackType = 0;
1845   t_L1Bit = false;
1846   fChain = tree;
1847   fCurrent = -1;
1848   if (!tree)
1849     return;
1850   fChain->SetMakeClass(1);
1851   fChain->SetBranchAddress("t_RunNo", &t_RunNo, &b_t_RunNo);
1852   fChain->SetBranchAddress("t_EventNo", &t_EventNo, &b_t_EventNo);
1853   fChain->SetBranchAddress("t_Tracks", &t_Tracks, &b_t_Tracks);
1854   fChain->SetBranchAddress("t_TracksProp", &t_TracksProp, &b_t_TracksProp);
1855   fChain->SetBranchAddress("t_TracksSaved", &t_TracksSaved, &b_t_TracksSaved);
1856   fChain->SetBranchAddress("t_TracksLoose", &t_TracksLoose, &b_t_TracksLoose);
1857   fChain->SetBranchAddress("t_TracksTight", &t_TracksTight, &b_t_TracksTight);
1858   fChain->SetBranchAddress("t_allvertex", &t_allvertex, &b_t_allvertex);
1859   fChain->SetBranchAddress("t_TrigPass", &t_TrigPass, &b_t_TrigPass);
1860   fChain->SetBranchAddress("t_TrigPassSel", &t_TrigPassSel, &b_t_TrigPassSel);
1861   fChain->SetBranchAddress("t_L1Bit", &t_L1Bit, &b_t_L1Bit);
1862   fChain->SetBranchAddress("t_hltbits", &t_hltbits, &b_t_hltbits);
1863   fChain->SetBranchAddress("t_ietaAll", &t_ietaAll, &b_t_ietaAll);
1864   fChain->SetBranchAddress("t_ietaGood", &t_ietaGood, &b_t_ietaGood);
1865   fChain->SetBranchAddress("t_trackType", &t_trackType, &b_t_trackType);
1866   Notify();
1867 
1868   if (std::string(dupFileName) != "") {
1869     std::ifstream infile(dupFileName);
1870     if (!infile.is_open()) {
1871       std::cout << "Cannot open " << dupFileName << std::endl;
1872     } else {
1873       while (1) {
1874         Long64_t jentry;
1875         infile >> jentry;
1876         if (!infile.good())
1877           break;
1878         entries_.push_back(jentry);
1879       }
1880       infile.close();
1881       std::cout << "Reads a list of " << entries_.size() << " events from " << dupFileName << std::endl;
1882     }
1883   }
1884 
1885   h_tk[0] = new TH1I("Track0", "# of tracks produced", 2000, 0, 2000);
1886   h_tk[1] = new TH1I("Track1", "# of tracks propagated", 2000, 0, 2000);
1887   h_tk[2] = new TH1I("Track2", "# of tracks saved in tree", 2000, 0, 2000);
1888   h_eta[0] = new TH1I("Eta0", "i#eta (All Tracks)", 60, -30, 30);
1889   h_eta[1] = new TH1I("Eta1", "i#eta (Good Tracks)", 60, -30, 30);
1890   h_eta[2] = new TH1I("Eta2", "i#eta (Loose Isolated Tracks)", 60, -30, 30);
1891   h_eta[3] = new TH1I("Eta3", "i#eta (Tight Isolated Tracks)", 60, -30, 30);
1892   h_eff[0] = new TH1D("Eff0", "i#eta (Selection Efficiency)", 60, -30, 30);
1893   h_eff[1] = new TH1D("Eff1", "i#eta (Loose Isolation Efficiency)", 60, -30, 30);
1894   h_eff[2] = new TH1D("Eff2", "i#eta (Tight Isolation Efficiency)", 60, -30, 30);
1895   h_pvx[0] = new TH1I("Pvx0", "Number of Good Vertex (All)", 100, 0, 100);
1896   h_pvx[1] = new TH1I("Pvx1", "Number of Good Vertex (Loose)", 100, 0, 100);
1897   h_pvx[2] = new TH1I("Pvx2", "Number of Good Vertex (Tight)", 100, 0, 100);
1898 }
1899 
1900 Bool_t GetEntries::Notify() {
1901   // The Notify() function is called when a new file is opened. This
1902   // can be either for a new TTree in a TChain or when when a new TTree
1903   // is started when using PROOF. It is normally not necessary to make changes
1904   // to the generated code, but the routine can be extended by the
1905   // user if needed. The return value is currently not used.
1906 
1907   return kTRUE;
1908 }
1909 
1910 void GetEntries::Show(Long64_t entry) {
1911   // Print contents of entry.
1912   // If entry is not specified, print current entry
1913   if (!fChain)
1914     return;
1915   fChain->Show(entry);
1916 }
1917 
1918 Int_t GetEntries::Cut(Long64_t) {
1919   // This function may be called from Loop.
1920   // returns  1 if entry is accepted.
1921   // returns -1 otherwise.
1922   return 1;
1923 }
1924 
1925 void GetEntries::Loop(Long64_t nmax) {
1926   //   In a ROOT session, you can do:
1927   //      Root > .L CalibMonitor.C+g
1928   //      Root > GetEntries t
1929   //      Root > t.GetEntry(12); // Fill t data members with entry number 12
1930   //      Root > t.Show();       // Show values of entry 12
1931   //      Root > t.Show(16);     // Read and show values of entry 16
1932   //      Root > t.Loop();       // Loop on all entries
1933   //
1934 
1935   //     This is the loop skeleton where:
1936   //    jentry is the global entry number in the chain
1937   //    ientry is the entry number in the current Tree
1938   //  Note that the argument to GetEntry must be:
1939   //    jentry for TChain::GetEntry
1940   //    ientry for TTree::GetEntry and TBranch::GetEntry
1941   //
1942   //       To read only selected branches, Insert statements like:
1943   // METHOD1:
1944   //    fChain->SetBranchStatus("*",0);  // disable all branches
1945   //    fChain->SetBranchStatus("branchname",1);  // activate branchname
1946   // METHOD2: replace line
1947   //    fChain->GetEntry(jentry);       //read all branches
1948   //by  b_branchname->GetEntry(ientry); //read only this branch
1949   if (fChain == 0)
1950     return;
1951 
1952   Long64_t nentries = fChain->GetEntriesFast();
1953   Long64_t nbytes = 0, nb = 0;
1954   unsigned int kount(0), duplicate(0), selected(0);
1955   int l1(0), hlt(0), loose(0), tight(0);
1956   int allHLT[3] = {0, 0, 0};
1957   int looseHLT[3] = {0, 0, 0};
1958   int tightHLT[3] = {0, 0, 0};
1959   Long64_t entries = (nmax > 0) ? nmax : nentries;
1960   for (Long64_t jentry = 0; jentry < entries; jentry++) {
1961     Long64_t ientry = LoadTree(jentry);
1962     if (ientry < 0)
1963       break;
1964     nb = fChain->GetEntry(jentry);
1965     nbytes += nb;
1966     bool select = (std::find(entries_.begin(), entries_.end(), jentry) == entries_.end());
1967     if (!select) {
1968       ++duplicate;
1969       continue;
1970     }
1971     h_tk[0]->Fill(t_Tracks);
1972     h_tk[1]->Fill(t_TracksProp);
1973     h_tk[2]->Fill(t_TracksSaved);
1974     h_pvx[0]->Fill(t_allvertex);
1975     if (t_TracksLoose > 0)
1976       h_pvx[1]->Fill(t_allvertex);
1977     if (t_TracksTight > 0)
1978       h_pvx[2]->Fill(t_allvertex);
1979     if (t_L1Bit) {
1980       ++l1;
1981       if (t_TracksLoose > 0)
1982         ++loose;
1983       if (t_TracksTight > 0)
1984         ++tight;
1985       if (t_TrigPass)
1986         ++hlt;
1987     }
1988     if (t_TrigPass) {
1989       ++kount;
1990       if (t_TrigPassSel)
1991         ++selected;
1992     }
1993     bool passt[2] = {false, false}, pass(false);
1994     for (unsigned k = 0; k < t_hltbits->size(); ++k) {
1995       if ((*t_hltbits)[k] > 0) {
1996         pass = true;
1997         for (int i = 0; i < 2; ++i)
1998           if (k == bit_[i])
1999             passt[i] = true;
2000       }
2001     }
2002     if (pass) {
2003       ++allHLT[0];
2004       for (int i = 0; i < 2; ++i)
2005         if (passt[i])
2006           ++allHLT[i + 1];
2007       if (t_TracksLoose > 0) {
2008         ++looseHLT[0];
2009         for (int i = 0; i < 2; ++i)
2010           if (passt[i])
2011             ++looseHLT[i + 1];
2012       }
2013       if (t_TracksTight > 0) {
2014         ++tightHLT[0];
2015         for (int i = 0; i < 2; ++i)
2016           if (passt[i])
2017             ++tightHLT[i + 1];
2018       }
2019     }
2020     for (unsigned int k = 0; k < t_ietaAll->size(); ++k)
2021       h_eta[0]->Fill((*t_ietaAll)[k]);
2022     for (unsigned int k = 0; k < t_ietaGood->size(); ++k) {
2023       h_eta[1]->Fill((*t_ietaGood)[k]);
2024       if (t_trackType->size() > 0) {
2025         if ((*t_trackType)[k] > 0)
2026           h_eta[2]->Fill((*t_ietaGood)[k]);
2027         if ((*t_trackType)[k] > 1)
2028           h_eta[3]->Fill((*t_ietaGood)[k]);
2029       }
2030     }
2031   }
2032   double ymaxk(0);
2033   for (int i = 1; i <= h_eff[0]->GetNbinsX(); ++i) {
2034     double rat(0), drat(0);
2035     if (h_eta[0]->GetBinContent(i) > ymaxk)
2036       ymaxk = h_eta[0]->GetBinContent(i);
2037     if ((h_eta[1]->GetBinContent(i) > 0) && (h_eta[0]->GetBinContent(i) > 0)) {
2038       rat = h_eta[1]->GetBinContent(i) / h_eta[0]->GetBinContent(i);
2039       drat = rat * std::sqrt(pow((h_eta[1]->GetBinError(i) / h_eta[1]->GetBinContent(i)), 2) +
2040                              pow((h_eta[0]->GetBinError(i) / h_eta[0]->GetBinContent(i)), 2));
2041     }
2042     h_eff[0]->SetBinContent(i, rat);
2043     h_eff[0]->SetBinError(i, drat);
2044     if ((h_eta[1]->GetBinContent(i) > 0) && (h_eta[2]->GetBinContent(i) > 0)) {
2045       rat = h_eta[2]->GetBinContent(i) / h_eta[1]->GetBinContent(i);
2046       drat = rat * std::sqrt(pow((h_eta[2]->GetBinError(i) / h_eta[2]->GetBinContent(i)), 2) +
2047                              pow((h_eta[1]->GetBinError(i) / h_eta[1]->GetBinContent(i)), 2));
2048     } else {
2049       rat = drat = 0;
2050     }
2051     h_eff[1]->SetBinContent(i, rat);
2052     h_eff[1]->SetBinError(i, drat);
2053     if ((h_eta[1]->GetBinContent(i) > 0) && (h_eta[3]->GetBinContent(i) > 0)) {
2054       rat = h_eta[3]->GetBinContent(i) / h_eta[1]->GetBinContent(i);
2055       drat = rat * std::sqrt(pow((h_eta[3]->GetBinError(i) / h_eta[3]->GetBinContent(i)), 2) +
2056                              pow((h_eta[1]->GetBinError(i) / h_eta[1]->GetBinContent(i)), 2));
2057     } else {
2058       rat = drat = 0;
2059     }
2060     h_eff[1]->SetBinContent(i, rat);
2061     h_eff[1]->SetBinError(i, drat);
2062   }
2063   std::cout << "===== Remove " << duplicate << " events from " << nentries << "\n===== " << kount
2064             << " events passed trigger of which " << selected << " events get selected =====\n"
2065             << std::endl;
2066   std::cout << "===== " << l1 << " events passed L1 " << hlt << " events passed HLT and " << loose << ":" << tight
2067             << " events have at least 1 track candidate with loose:tight"
2068             << " isolation cut =====\n"
2069             << std::endl;
2070   for (int i = 0; i < 3; ++i) {
2071     char tbit[20];
2072     if (i == 0)
2073       sprintf(tbit, "Any");
2074     else
2075       sprintf(tbit, "%3d", bit_[i - 1]);
2076     std::cout << "Satisfying HLT trigger bit " << tbit << " Kount " << allHLT[i] << " Loose " << looseHLT[i]
2077               << " Tight " << tightHLT[i] << std::endl;
2078   }
2079 
2080   gStyle->SetCanvasBorderMode(0);
2081   gStyle->SetCanvasColor(kWhite);
2082   gStyle->SetPadColor(kWhite);
2083   gStyle->SetFillColor(kWhite);
2084   gStyle->SetOptStat(1110);
2085   gStyle->SetOptTitle(0);
2086   int color[5] = {kBlack, kRed, kBlue, kMagenta, kCyan};
2087   int lines[5] = {1, 2, 3, 4, 5};
2088   TCanvas *pad1 = new TCanvas("c_track", "c_track", 500, 500);
2089   pad1->SetRightMargin(0.10);
2090   pad1->SetTopMargin(0.10);
2091   pad1->SetFillColor(kWhite);
2092   std::string titl1[3] = {"Reconstructed", "Propagated", "Saved"};
2093   TLegend *legend1 = new TLegend(0.11, 0.80, 0.50, 0.89);
2094   legend1->SetFillColor(kWhite);
2095   double ymax(0), xmax(0);
2096   for (int k = 0; k < 3; ++k) {
2097     int total(0), totaltk(0);
2098     for (int i = 1; i <= h_tk[k]->GetNbinsX(); ++i) {
2099       if (ymax < h_tk[k]->GetBinContent(i))
2100         ymax = h_tk[k]->GetBinContent(i);
2101       if (i > 1)
2102         total += (int)(h_tk[k]->GetBinContent(i));
2103       totaltk += (int)(h_tk[k]->GetBinContent(i)) * (i - 1);
2104       if (h_tk[k]->GetBinContent(i) > 0) {
2105         if (xmax < h_tk[k]->GetBinLowEdge(i) + h_tk[k]->GetBinWidth(i))
2106           xmax = h_tk[k]->GetBinLowEdge(i) + h_tk[k]->GetBinWidth(i);
2107       }
2108     }
2109     h_tk[k]->SetLineColor(color[k]);
2110     h_tk[k]->SetMarkerColor(color[k]);
2111     h_tk[k]->SetLineStyle(lines[k]);
2112     std::cout << h_tk[k]->GetTitle() << " Entries " << h_tk[k]->GetEntries() << " Events " << total << " Tracks "
2113               << totaltk << std::endl;
2114     legend1->AddEntry(h_tk[k], titl1[k].c_str(), "l");
2115   }
2116   int i1 = (int)(0.1 * xmax) + 1;
2117   xmax = 10.0 * i1;
2118   int i2 = (int)(0.01 * ymax) + 1;
2119 
2120   ymax = 100.0 * i2;
2121   for (int k = 0; k < 3; ++k) {
2122     h_tk[k]->GetXaxis()->SetRangeUser(0, xmax);
2123     h_tk[k]->GetYaxis()->SetRangeUser(0.1, ymax);
2124     h_tk[k]->GetXaxis()->SetTitle("# Tracks");
2125     h_tk[k]->GetYaxis()->SetTitle("Events");
2126     h_tk[k]->GetYaxis()->SetLabelOffset(0.005);
2127     h_tk[k]->GetYaxis()->SetTitleOffset(1.20);
2128     if (k == 0)
2129       h_tk[k]->Draw("hist");
2130     else
2131       h_tk[k]->Draw("hist sames");
2132   }
2133   pad1->Update();
2134   pad1->SetLogy();
2135   ymax = 0.90;
2136   for (int k = 0; k < 3; ++k) {
2137     TPaveStats *st1 = (TPaveStats *)h_tk[k]->GetListOfFunctions()->FindObject("stats");
2138     if (st1 != NULL) {
2139       st1->SetLineColor(color[k]);
2140       st1->SetTextColor(color[k]);
2141       st1->SetY1NDC(ymax - 0.09);
2142       st1->SetY2NDC(ymax);
2143       st1->SetX1NDC(0.55);
2144       st1->SetX2NDC(0.90);
2145       ymax -= 0.09;
2146     }
2147   }
2148   pad1->Modified();
2149   pad1->Update();
2150   legend1->Draw("same");
2151   pad1->Update();
2152 
2153   TCanvas *pad2 = new TCanvas("c_ieta", "c_ieta", 500, 500);
2154   pad2->SetRightMargin(0.10);
2155   pad2->SetTopMargin(0.10);
2156   pad2->SetFillColor(kWhite);
2157   pad2->SetLogy();
2158   std::string titl2[4] = {"All Tracks", "Selected Tracks", "Loose Isolation", "Tight Isolation"};
2159   TLegend *legend2 = new TLegend(0.11, 0.75, 0.50, 0.89);
2160   legend2->SetFillColor(kWhite);
2161   i2 = (int)(0.001 * ymaxk) + 1;
2162   ymax = 1000.0 * i2;
2163   for (int k = 0; k < 4; ++k) {
2164     h_eta[k]->GetYaxis()->SetRangeUser(1, ymax);
2165     h_eta[k]->SetLineColor(color[k]);
2166     h_eta[k]->SetMarkerColor(color[k]);
2167     h_eta[k]->SetLineStyle(lines[k]);
2168     h_eta[k]->GetXaxis()->SetTitle("i#eta");
2169     h_eta[k]->GetYaxis()->SetTitle("Tracks");
2170     h_eta[k]->GetYaxis()->SetLabelOffset(0.005);
2171     h_eta[k]->GetYaxis()->SetTitleOffset(1.20);
2172     legend2->AddEntry(h_eta[k], titl2[k].c_str(), "l");
2173     if (k == 0)
2174       h_eta[k]->Draw("hist");
2175     else
2176       h_eta[k]->Draw("hist sames");
2177   }
2178   pad2->Update();
2179   ymax = 0.90;
2180   //double ymin = 0.10;
2181   for (int k = 0; k < 4; ++k) {
2182     TPaveStats *st1 = (TPaveStats *)h_eta[k]->GetListOfFunctions()->FindObject("stats");
2183     if (st1 != NULL) {
2184       st1->SetLineColor(color[k]);
2185       st1->SetTextColor(color[k]);
2186       st1->SetY1NDC(ymax - 0.09);
2187       st1->SetY2NDC(ymax);
2188       st1->SetX1NDC(0.55);
2189       st1->SetX2NDC(0.90);
2190       ymax -= 0.09;
2191     }
2192   }
2193   pad2->Modified();
2194   pad2->Update();
2195   legend2->Draw("same");
2196   pad2->Update();
2197 
2198   std::string titl3[3] = {"Selection", "Loose Isolation", "Tight Isolation"};
2199   TCanvas *pad3 = new TCanvas("c_effi", "c_effi", 500, 500);
2200   TLegend *legend3 = new TLegend(0.11, 0.785, 0.50, 0.89);
2201   pad3->SetRightMargin(0.10);
2202   pad3->SetTopMargin(0.10);
2203   pad3->SetFillColor(kWhite);
2204   pad3->SetLogy();
2205   for (int k = 0; k < 3; ++k) {
2206     h_eff[k]->SetStats(0);
2207     h_eff[k]->SetMarkerStyle(20);
2208     h_eff[k]->SetLineColor(color[k]);
2209     h_eff[k]->SetMarkerColor(color[k]);
2210     h_eff[k]->GetXaxis()->SetTitle("i#eta");
2211     h_eff[k]->GetYaxis()->SetTitle("Efficiency");
2212     h_eff[k]->GetYaxis()->SetLabelOffset(0.005);
2213     h_eff[k]->GetYaxis()->SetTitleOffset(1.20);
2214     if (k == 0)
2215       h_eff[k]->Draw("");
2216     else
2217       h_eff[k]->Draw("same");
2218     legend3->AddEntry(h_eff[k], titl3[k].c_str(), "l");
2219   }
2220   pad3->Modified();
2221   pad3->Update();
2222   legend3->Draw("same");
2223   pad3->Update();
2224 
2225   std::string titl4[3] = {"All events", "Events with Loose Isolation", "Events with Tight Isolation"};
2226   TLegend *legend4 = new TLegend(0.11, 0.785, 0.50, 0.89);
2227   TCanvas *pad4 = new TCanvas("c_nvx", "c_nvx", 500, 500);
2228   pad4->SetRightMargin(0.10);
2229   pad4->SetTopMargin(0.10);
2230   pad4->SetFillColor(kWhite);
2231   pad4->SetLogy();
2232   for (int k = 0; k < 3; ++k) {
2233     h_pvx[k]->SetStats(1110);
2234     h_pvx[k]->SetMarkerStyle(20);
2235     h_pvx[k]->SetLineColor(color[k]);
2236     h_pvx[k]->SetMarkerColor(color[k]);
2237     h_pvx[k]->GetXaxis()->SetTitle("N_{PV}");
2238     h_pvx[k]->GetYaxis()->SetTitle("Events");
2239     h_pvx[k]->GetYaxis()->SetLabelOffset(0.005);
2240     h_pvx[k]->GetYaxis()->SetTitleOffset(1.20);
2241     if (k == 0)
2242       h_pvx[k]->Draw("");
2243     else
2244       h_pvx[k]->Draw("sames");
2245     legend4->AddEntry(h_pvx[k], titl4[k].c_str(), "l");
2246   }
2247   pad4->Update();
2248   ymax = 0.90;
2249   for (int k = 0; k < 3; ++k) {
2250     TPaveStats *st1 = (TPaveStats *)h_pvx[k]->GetListOfFunctions()->FindObject("stats");
2251     if (st1 != NULL) {
2252       st1->SetLineColor(color[k]);
2253       st1->SetTextColor(color[k]);
2254       st1->SetY1NDC(ymax - 0.09);
2255       st1->SetY2NDC(ymax);
2256       st1->SetX1NDC(0.55);
2257       st1->SetX2NDC(0.90);
2258       ymax -= 0.09;
2259     }
2260   }
2261   pad4->Modified();
2262   pad4->Update();
2263   legend4->Draw("same");
2264   pad4->Update();
2265 }