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