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