HcalCalibAlgos/macros/CalibMonitor.C

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