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