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 //                  dataMC, truncateFlag, useGen, scale, useScale, etalo, etahi,
0007 //                  runlo, runhi, phimin, phimax, zside, nvxlo, nvxhi, rbx,
0008 //                  exclude, etamax);
0009 //  c1.Loop();
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();
0017 //
0018 //         This looks into the tree *EventInfo* and can provide a set
0019 //         of histograms with event statistics
0020 //
0021 //  CalibPlotProperties cc(fname, dirname, dupFileName, prefix, corrFileName,
0022 //                 rcorFileName, puCorr, flagC, dataMC, truncateFlag,
0023 //                         useGen, scale, useScale, etalo, etahi, runlo, runhi,
0024 //                 phimin, phimax, zside, nvxlo, nvxhi, rbx, exclude,
0025 //                         etamax);
0026 //  cc.Loop();
0027 //  cc.savePlot(histFileName,append,all);
0028 //
0029 //         This makes different style of plots to understand the quality of
0030 //         reconstruction rather than the response plots which are done by
0031 //         CalibMonitor
0032 //
0033 //   where:
0034 //
0035 //   fname   (const char*)     = file name of the input ROOT tree
0036 //                               or name of the file containing a list of
0037 //                               file names of input ROOT trees
0038 //   dirname (std::string)     = name of the directory where Tree resides
0039 //                               (use "HcalIsoTrkAnalyzer")
0040 //   dupFileName (char*)       = name of the file containing list of entries
0041 //                               of duplicate events
0042 //   comFileName (char*)       = name of the file with list of run and event
0043 //                               number to be selected
0044 //   outFileName (char*)       = name of a text file to be created (under
0045 //                               control of value of flag) with information
0046 //                               about events
0047 //   prefix (std::string)      = String to be added to the name of histogram
0048 //                               (usually a 4 character string; default="")
0049 //   corrFileName (char*)      = name of the text file having the correction
0050 //                               factors to be used (default="", no corr.)
0051 //   rcorFileName (char*)      = name of the text file having the correction
0052 //                               factors as a function of run numbers or depth
0053 //                               or entry number to be used for raddam/depth/
0054 //                               pileup dependent correction  (default="",
0055 //                               no corr.)
0056 //   puCorr (int)              = PU correction to be applied or not: 0 no
0057 //                               correction; < 0 use eDelta; > 0 rho dependent
0058 //                               correction (-2)
0059 //   flag (int)                = 7 digit integer (xmlthdo) with control
0060 //                               information (x=3/2/1/0 for having 1000/500/50/
0061 //                               100 bins for response distribution in (0:5);
0062 //                               m=1/0 for (not) making plots for each RBX;
0063 //                               l=3/2/1/0 for type of rcorFileName (3 for
0064 //                               pileup correction using machine learning
0065 //                               method; 2 for overall response corrections;
0066 //                               1 for depth dependence corrections;
0067 //                               0 for raddam corrections);
0068 //                               t = bit information (lower bit set will
0069 //                               apply a cut on L1 closeness; and higher bit
0070 //                               set read correction file with Marina format);
0071 //                               h = 0/1/2 for not creating / creating in
0072 //                               recreate mode / creating in append mode
0073 //                               the output text file;
0074 //                               d = 0/1/2/3 produces 3 standard (0,1,2) or
0075 //                               extended (3) set of histograms;
0076 //                               o = 0/1/2 for tight / loose / flexible
0077 //                               selection). Default = 1031
0078 //   flagC (int)               = 6 digit integer (mlthdo) with control
0079 //                               information (m=1/0 for making or not rechit
0080 //                               energy distributions;  l=2/1/0 for type of
0081 //                               correction (2 for overall response corrections;
0082 //                               1 for depth dependence corrections; 0 for
0083 //                               raddam corrections); t=bit information (lower
0084 //                               bit set will apply a cut on L1 closeness; and
0085 //                               higher bit set read correction file with
0086 //                               Marina format); h=1/0 for making or not
0087 //                               plots of momentum and total energies in the
0088 //                               two calorimeters ECAL/HCAL; d=1/0 for making
0089 //                               plots of momentum and eta distributions;
0090 //                               o = 0/1/2 for tight / loose / flexible
0091 //                               selection). Default = 1031
0092 //   numb   (int)              = number of eta bins (50 for -25:25)
0093 //   dataMC (bool)             = true/false for data/MC (default true)
0094 //   truncateFlag    (int)     = Flag to treat different depths differently (0)
0095 //                               both depths of ieta 15, 16 of HB as depth 1 (1)
0096 //                               all depths as depth 1 (2), all depths in HE
0097 //                               with values > 1 as depth 2 (3), all depths in
0098 //                               HB with values > 1 as depth 2 (4), all depths
0099 //                               in HB and HE with values > 1 as depth 2 (5)
0100 //                               (default = 0)
0101 //   useGen (bool)             = true/false to use generator level momentum
0102 //                               or reconstruction level momentum
0103 //                               (default = false)
0104 //   scale (double)            = energy scale if correction factor to be used
0105 //                               (default = 1.0)
0106 //   useScale (int)            = application of scale factor (0: nowehere,
0107 //                               1: barrel; 2: endcap, 3: everywhere)
0108 //                               barrel => |ieta| < 16; endcap => |ieta| > 15
0109 //                               (default = 0)
0110 //   etalo/etahi (int,int)     = |eta| ranges (default = 0:30)
0111 //   runlo  (int)              = lower value of run number to be included (+ve)
0112 //                               or excluded (-ve) (default = 0)
0113 //   runhi  (int)              = higher value of run number to be included
0114 //                               (+ve) or excluded (-ve) (default = 9999999)
0115 //   phimin          (int)     = minimum iphi value (default = 1)
0116 //   phimax          (int)     = maximum iphi value (default = 72)
0117 //   zside           (int)     = the side of the detector if phimin and phimax
0118 //                               differ from 1-72 (default = 1)
0119 //   nvxlo           (int)     = minimum # of vertex in event to be used
0120 //                               (default = 0)
0121 //   nvxhi           (int)     = maximum # of vertex in event to be used
0122 //                               (default = 1000)
0123 //   rbx             (int)     = zside*(Subdet*100+RBX #) to be consdered
0124 //                               (default = 0). For HEP17 it will be 217
0125 //   exclude         (bool)    = RBX specified by *rbx* to be exluded or only
0126 //                               considered (default = false)
0127 //   etamax          (bool)    = if set and if the corr-factor not found in the
0128 //                               corrFactor table, the corr-factor for the
0129 //                               corresponding zside, depth=1 and maximum ieta
0130 //                               in the table is taken (default = false)
0131 //
0132 //   histFileName (std::string)= name of the file containing saved histograms
0133 //   append (bool)             = true/false if the histogram file to be opened
0134 //                               in append/output mode (default = true)
0135 //   all (bool)                = true/false if all histograms to be saved or
0136 //                               not (default = false)
0137 //
0138 //   bitX (unsigned int)       = bit number of the HLT used in the selection
0139 //        (X = 1, 2)             for example the bits of HLT_IsoTrackHB(HE)
0140 //////////////////////////////////////////////////////////////////////////////
0141
0142 #include <TROOT.h>
0143 #include <TChain.h>
0144 #include <TFile.h>
0145 #include <TF1.h>
0146 #include <TH1D.h>
0147 #include <TH2F.h>
0148 #include <TProfile.h>
0149 #include <TFitResult.h>
0150 #include <TFitResultPtr.h>
0151 #include <TStyle.h>
0152 #include <TCanvas.h>
0153 #include <TLegend.h>
0154 #include <TPaveStats.h>
0155 #include <TPaveText.h>
0156
0157 #include <algorithm>
0158 #include <iomanip>
0159 #include <iostream>
0160 #include <fstream>
0161 #include <map>
0162 #include <vector>
0163 #include <string>
0164
0165 #include "CalibCorr.C"
0166
0167 class CalibMonitor {
0168 public:
0169   TChain *fChain;  //!pointer to the analyzed TTree or TChain
0170   Int_t fCurrent;  //!current Tree number in a TChain
0171
0172   // Declaration of leaf types
0173   Int_t t_Run;
0174   Int_t t_Event;
0175   Int_t t_DataType;
0176   Int_t t_ieta;
0177   Int_t t_iphi;
0178   Double_t t_EventWeight;
0179   Int_t t_nVtx;
0180   Int_t t_nTrk;
0181   Int_t t_goodPV;
0182   Double_t t_l1pt;
0183   Double_t t_l1eta;
0184   Double_t t_l1phi;
0185   Double_t t_l3pt;
0186   Double_t t_l3eta;
0187   Double_t t_l3phi;
0188   Double_t t_p;
0189   Double_t t_pt;
0190   Double_t t_phi;
0191   Double_t t_mindR1;
0192   Double_t t_mindR2;
0193   Double_t t_eMipDR;
0194   Double_t t_eHcal;
0195   Double_t t_eHcal10;
0196   Double_t t_eHcal30;
0197   Double_t t_hmaxNearP;
0198   Double_t t_rhoh;
0199   Bool_t t_selectTk;
0200   Bool_t t_qltyFlag;
0201   Bool_t t_qltyMissFlag;
0202   Bool_t t_qltyPVFlag;
0203   Double_t t_gentrackP;
0204   std::vector<unsigned int> *t_DetIds;
0205   std::vector<double> *t_HitEnergies;
0206   std::vector<bool> *t_trgbits;
0207   std::vector<unsigned int> *t_DetIds1;
0208   std::vector<unsigned int> *t_DetIds3;
0209   std::vector<double> *t_HitEnergies1;
0210   std::vector<double> *t_HitEnergies3;
0211
0212   // List of branches
0213   TBranch *b_t_Run;           //!
0214   TBranch *b_t_Event;         //!
0215   TBranch *b_t_DataType;      //!
0216   TBranch *b_t_ieta;          //!
0217   TBranch *b_t_iphi;          //!
0218   TBranch *b_t_EventWeight;   //!
0219   TBranch *b_t_nVtx;          //!
0220   TBranch *b_t_nTrk;          //!
0221   TBranch *b_t_goodPV;        //!
0222   TBranch *b_t_l1pt;          //!
0223   TBranch *b_t_l1eta;         //!
0224   TBranch *b_t_l1phi;         //!
0225   TBranch *b_t_l3pt;          //!
0226   TBranch *b_t_l3eta;         //!
0227   TBranch *b_t_l3phi;         //!
0228   TBranch *b_t_p;             //!
0229   TBranch *b_t_pt;            //!
0230   TBranch *b_t_phi;           //!
0231   TBranch *b_t_mindR1;        //!
0232   TBranch *b_t_mindR2;        //!
0233   TBranch *b_t_eMipDR;        //!
0234   TBranch *b_t_eHcal;         //!
0235   TBranch *b_t_eHcal10;       //!
0236   TBranch *b_t_eHcal30;       //!
0237   TBranch *b_t_hmaxNearP;     //!
0238   TBranch *b_t_rhoh;          //!
0239   TBranch *b_t_selectTk;      //!
0240   TBranch *b_t_qltyFlag;      //!
0241   TBranch *b_t_qltyMissFlag;  //!
0242   TBranch *b_t_qltyPVFlag;    //!
0243   TBranch *b_t_gentrackP;     //!
0244   TBranch *b_t_DetIds;        //!
0245   TBranch *b_t_HitEnergies;   //!
0246   TBranch *b_t_trgbits;       //!
0247   TBranch *b_t_DetIds1;       //!
0248   TBranch *b_t_DetIds3;       //!
0249   TBranch *b_t_HitEnergies1;  //!
0250   TBranch *b_t_HitEnergies3;  //!
0251
0252   struct counter {
0253     static const int npsize = 4;
0254     counter() {
0255       total = 0;
0256       for (int k = 0; k < npsize; ++k)
0257         count[k] = 0;
0258     };
0259     unsigned int total, count[npsize];
0260   };
0261
0262   CalibMonitor(const char *fname,
0263                const std::string &dirname,
0264                const char *dupFileName,
0265                const char *comFileName,
0266                const char *outFileName,
0267                const std::string &prefix = "",
0268                const char *corrFileName = "",
0269                const char *rcorFileName = "",
0270                int puCorr = -2,
0271                int flag = 1031,
0272                int numb = 50,
0273                bool datMC = true,
0274                int truncateFlag = 0,
0275                bool useGen = false,
0276                double scale = 1.0,
0277                int useScale = 0,
0278                int etalo = 0,
0279                int etahi = 30,
0280                int runlo = 0,
0281                int runhi = 99999999,
0282                int phimin = 1,
0283                int phimax = 72,
0284                int zside = 1,
0285                int nvxlo = 0,
0286                int nvxhi = 1000,
0287                int rbx = 0,
0288                bool exclude = false,
0289                bool etamax = false);
0290   virtual ~CalibMonitor();
0291   virtual Int_t Cut(Long64_t entry);
0292   virtual Int_t GetEntry(Long64_t entry);
0293   virtual Long64_t LoadTree(Long64_t entry);
0294   virtual void Init(TChain *, const char *, const char *, const char *);
0295   virtual void Loop();
0296   virtual Bool_t Notify();
0297   virtual void Show(Long64_t entry = -1);
0298   bool goodTrack(double &eHcal, double &cut, const Long64_t &entry, bool debug);
0299   bool selectPhi(bool debug);
0300   void plotHist(int type, int num, bool save = false);
0301   template <class Hist>
0302   void drawHist(Hist *, TCanvas *);
0303   void savePlot(const std::string &theName, bool append = true, bool all = false);
0304   void correctEnergy(double &ener, const Long64_t &entry);
0305
0306 private:
0307   static const unsigned int npbin = 6, kp50 = 3;
0308   CalibCorrFactor *corrFactor_;
0309   CalibCorr *cFactor_;
0310   CalibSelectRBX *cSelect_;
0311   const std::string fname_, dirnm_, prefix_, outFileName_;
0312   const int corrPU_, flag_, numb_;
0313   const bool dataMC_, useGen_;
0314   const int truncateFlag_;
0315   const int etalo_, etahi_;
0316   int runlo_, runhi_;
0317   const int phimin_, phimax_, zside_, nvxlo_, nvxhi_, rbx_;
0318   bool exclude_, corrE_, cutL1T_, selRBX_;
0319   bool includeRun_;
0320   int coarseBin_, etamp_, etamn_, plotType_;
0321   int flexibleSelect_, ifDepth_;
0322   double log2by18_;
0323   std::ofstream fileout_;
0324   std::vector<Long64_t> entries_;
0325   std::vector<std::pair<int, int> > events_;
0326   std::vector<double> etas_, ps_, dl1_;
0327   std::vector<int> nvx_, ietasL_, ietasH_;
0328   std::vector<TH1D *> h_rbx, h_etaF[npbin], h_etaB[npbin];
0329   std::vector<TProfile *> h_etaX[npbin];
0330   std::vector<TH1D *> h_etaR[npbin], h_nvxR[npbin], h_dL1R[npbin];
0331   std::vector<TH1D *> h_pp[npbin];
0332 };
0333
0334 CalibMonitor::CalibMonitor(const char *fname,
0335                            const std::string &dirnm,
0336                            const char *dupFileName,
0337                            const char *comFileName,
0338                            const char *outFName,
0339                            const std::string &prefix,
0340                            const char *corrFileName,
0341                            const char *rcorFileName,
0342                            int puCorr,
0343                            int flag,
0344                            int numb,
0345                            bool dataMC,
0346                            int truncate,
0347                            bool useGen,
0348                            double scale,
0349                            int useScale,
0350                            int etalo,
0351                            int etahi,
0352                            int runlo,
0353                            int runhi,
0354                            int phimin,
0355                            int phimax,
0356                            int zside,
0357                            int nvxlo,
0358                            int nvxhi,
0359                            int rbx,
0360                            bool exc,
0361                            bool etam)
0362     : corrFactor_(nullptr),
0363       cFactor_(nullptr),
0364       cSelect_(nullptr),
0365       fname_(fname),
0366       dirnm_(dirnm),
0367       prefix_(prefix),
0368       outFileName_(std::string(outFName)),
0369       corrPU_(puCorr),
0370       flag_(flag),
0371       numb_(numb),
0372       dataMC_(dataMC),
0373       useGen_(useGen),
0374       truncateFlag_(truncate),
0375       etalo_(etalo),
0376       etahi_(etahi),
0377       runlo_(runlo),
0378       runhi_(runhi),
0379       phimin_(phimin),
0380       phimax_(phimax),
0381       zside_(zside),
0382       nvxlo_(nvxlo),
0383       nvxhi_(nvxhi),
0384       rbx_(rbx),
0385       exclude_(exc),
0386       includeRun_(true) {
0387   // if parameter tree is not specified (or zero), connect the file
0388   // used to generate this class and read the Tree
0389
0390   plotType_ = ((flag_ / 10) % 10);
0391   if (plotType_ < 0 || plotType_ > 3)
0392     plotType_ = 3;
0393   flexibleSelect_ = (((flag_ / 1) % 10));
0394   int oneplace = ((flag_ / 1000) % 10);
0395   cutL1T_ = (oneplace % 2);
0396   bool marina = ((oneplace / 2) % 2);
0397   ifDepth_ = ((flag_ / 10000) % 10);
0398   selRBX_ = (((flag_ / 100000) % 10) > 0);
0399   coarseBin_ = ((flag_ / 1000000) % 10);
0400   log2by18_ = std::log(2.5) / 18.0;
0401   if (runlo_ < 0 || runhi_ < 0) {
0402     runlo_ = std::abs(runlo_);
0403     runhi_ = std::abs(runhi_);
0404     includeRun_ = false;
0405   }
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 << 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, useScale, scale, etam, marina, false);
0420     Init(chain, dupFileName, comFileName, outFName);
0421     if (std::string(rcorFileName) != "") {
0422       cFactor_ = new CalibCorr(rcorFileName, ifDepth_, false);
0423     } else {
0424       ifDepth_ = 0;
0425     }
0426     if (rbx != 0)
0427       cSelect_ = new CalibSelectRBX(rbx, false);
0428   }
0429 }
0430
0431 CalibMonitor::~CalibMonitor() {
0432   delete corrFactor_;
0433   delete cFactor_;
0434   delete cSelect_;
0435   if (!fChain)
0436     return;
0437   delete fChain->GetCurrentFile();
0438 }
0439
0440 Int_t CalibMonitor::GetEntry(Long64_t entry) {
0441   // Read contents of entry.
0442   if (!fChain)
0443     return 0;
0444   return fChain->GetEntry(entry);
0445 }
0446
0447 Long64_t CalibMonitor::LoadTree(Long64_t entry) {
0448   // Set the environment to read one entry
0449   if (!fChain)
0450     return -5;
0451   Long64_t centry = fChain->LoadTree(entry);
0452   if (centry < 0)
0453     return centry;
0454   if (!fChain->InheritsFrom(TChain::Class()))
0455     return centry;
0456   TChain *chain = (TChain *)fChain;
0457   if (chain->GetTreeNumber() != fCurrent) {
0458     fCurrent = chain->GetTreeNumber();
0459     Notify();
0460   }
0461   return centry;
0462 }
0463
0464 void CalibMonitor::Init(TChain *tree, const char *dupFileName, const char *comFileName, const char *outFileName) {
0465   // The Init() function is called when the selector needs to initialize
0466   // a new tree or chain. Typically here the branch addresses and branch
0467   // pointers of the tree will be set.
0468   // It is normally not necessary to make changes to the generated
0469   // code, but the routine can be extended by the user if needed.
0470   // Init() will be called many times when running on PROOF
0471   // (once per file to be processed).
0472
0473   // Set object pointer
0474   t_DetIds = 0;
0475   t_DetIds1 = 0;
0476   t_DetIds3 = 0;
0477   t_HitEnergies = 0;
0478   t_HitEnergies1 = 0;
0479   t_HitEnergies3 = 0;
0480   t_trgbits = 0;
0481   // Set branch addresses and branch pointers
0482   fChain = tree;
0483   fCurrent = -1;
0484   if (!tree)
0485     return;
0486   fChain->SetMakeClass(1);
0487
0488   fChain->SetBranchAddress("t_Run", &t_Run, &b_t_Run);
0489   fChain->SetBranchAddress("t_Event", &t_Event, &b_t_Event);
0490   fChain->SetBranchAddress("t_DataType", &t_DataType, &b_t_DataType);
0491   fChain->SetBranchAddress("t_ieta", &t_ieta, &b_t_ieta);
0492   fChain->SetBranchAddress("t_iphi", &t_iphi, &b_t_iphi);
0493   fChain->SetBranchAddress("t_EventWeight", &t_EventWeight, &b_t_EventWeight);
0494   fChain->SetBranchAddress("t_nVtx", &t_nVtx, &b_t_nVtx);
0495   fChain->SetBranchAddress("t_nTrk", &t_nTrk, &b_t_nTrk);
0496   fChain->SetBranchAddress("t_goodPV", &t_goodPV, &b_t_goodPV);
0497   fChain->SetBranchAddress("t_l1pt", &t_l1pt, &b_t_l1pt);
0498   fChain->SetBranchAddress("t_l1eta", &t_l1eta, &b_t_l1eta);
0499   fChain->SetBranchAddress("t_l1phi", &t_l1phi, &b_t_l1phi);
0500   fChain->SetBranchAddress("t_l3pt", &t_l3pt, &b_t_l3pt);
0501   fChain->SetBranchAddress("t_l3eta", &t_l3eta, &b_t_l3eta);
0502   fChain->SetBranchAddress("t_l3phi", &t_l3phi, &b_t_l3phi);
0503   fChain->SetBranchAddress("t_p", &t_p, &b_t_p);
0504   fChain->SetBranchAddress("t_pt", &t_pt, &b_t_pt);
0505   fChain->SetBranchAddress("t_phi", &t_phi, &b_t_phi);
0506   fChain->SetBranchAddress("t_mindR1", &t_mindR1, &b_t_mindR1);
0507   fChain->SetBranchAddress("t_mindR2", &t_mindR2, &b_t_mindR2);
0508   fChain->SetBranchAddress("t_eMipDR", &t_eMipDR, &b_t_eMipDR);
0509   fChain->SetBranchAddress("t_eHcal", &t_eHcal, &b_t_eHcal);
0510   fChain->SetBranchAddress("t_eHcal10", &t_eHcal10, &b_t_eHcal10);
0511   fChain->SetBranchAddress("t_eHcal30", &t_eHcal30, &b_t_eHcal30);
0512   fChain->SetBranchAddress("t_hmaxNearP", &t_hmaxNearP, &b_t_hmaxNearP);
0513   fChain->SetBranchAddress("t_rhoh", &t_rhoh, &b_t_rhoh);
0514   fChain->SetBranchAddress("t_selectTk", &t_selectTk, &b_t_selectTk);
0515   fChain->SetBranchAddress("t_qltyFlag", &t_qltyFlag, &b_t_qltyFlag);
0516   fChain->SetBranchAddress("t_qltyMissFlag", &t_qltyMissFlag, &b_t_qltyMissFlag);
0517   fChain->SetBranchAddress("t_qltyPVFlag", &t_qltyPVFlag, &b_t_qltyPVFlag);
0518   fChain->SetBranchAddress("t_gentrackP", &t_gentrackP, &b_t_gentrackP);
0519   fChain->SetBranchAddress("t_DetIds", &t_DetIds, &b_t_DetIds);
0520   fChain->SetBranchAddress("t_HitEnergies", &t_HitEnergies, &b_t_HitEnergies);
0521   fChain->SetBranchAddress("t_trgbits", &t_trgbits, &b_t_trgbits);
0522   fChain->SetBranchAddress("t_DetIds1", &t_DetIds1, &b_t_DetIds1);
0523   fChain->SetBranchAddress("t_DetIds3", &t_DetIds3, &b_t_DetIds3);
0524   fChain->SetBranchAddress("t_HitEnergies1", &t_HitEnergies1, &b_t_HitEnergies1);
0525   fChain->SetBranchAddress("t_HitEnergies3", &t_HitEnergies3, &b_t_HitEnergies3);
0526   Notify();
0527
0528   if (strcmp(dupFileName, "") != 0) {
0529     ifstream infil1(dupFileName);
0530     if (!infil1.is_open()) {
0531       std::cout << "Cannot open duplicate file " << dupFileName << std::endl;
0532     } else {
0533       while (1) {
0534         Long64_t jentry;
0535         infil1 >> jentry;
0536         if (!infil1.good())
0537           break;
0538         entries_.push_back(jentry);
0539       }
0540       infil1.close();
0541       std::cout << "Reads a list of " << entries_.size() << " events from " << dupFileName << std::endl;
0542     }
0543   } else {
0544     std::cout << "No duplicate events in the input file" << std::endl;
0545   }
0546
0547   if (strcmp(comFileName, "") != 0) {
0548     ifstream infil2(comFileName);
0549     if (!infil2.is_open()) {
0550       std::cout << "Cannot open selection file " << comFileName << std::endl;
0551     } else {
0552       while (1) {
0553         int irun, ievt;
0554         infil2 >> irun >> ievt;
0555         if (!infil2.good())
0556           break;
0557         std::pair<int, int> key(irun, ievt);
0558         events_.push_back(key);
0559       }
0560       infil2.close();
0561       std::cout << "Reads a list of " << events_.size() << " events from " << comFileName << std::endl;
0562     }
0563   } else {
0564     std::cout << "No event list provided for selection" << std::endl;
0565   }
0566
0567   if (((flag_ / 100) % 10) > 0) {
0568     if (((flag_ / 100) % 10) == 2) {
0569       fileout_.open(outFileName, std::ofstream::out);
0570       std::cout << "Opens " << outFileName << " in output mode" << std::endl;
0571     } else {
0572       fileout_.open(outFileName, std::ofstream::app);
0573       std::cout << "Opens " << outFileName << " in append mode" << std::endl;
0574     }
0575     fileout_ << "Input file: " << fname_ << " Directory: " << dirnm_ << " Prefix: " << prefix_ << std::endl;
0576   }
0577
0578   double xbins[99];
0579   int nbins(-1);
0580   if (plotType_ == 0) {
0581     double xbin[9] = {-21.0, -16.0, -12.0, -6.0, 0.0, 6.0, 12.0, 16.0, 21.0};
0582     for (int i = 0; i < 9; ++i) {
0583       etas_.push_back(xbin[i]);
0584       xbins[i] = xbin[i];
0585     }
0586     nbins = 8;
0587   } else if (plotType_ == 1) {
0588     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};
0589     for (int i = 0; i < 11; ++i) {
0590       etas_.push_back(xbin[i]);
0591       xbins[i] = xbin[i];
0592     }
0593     nbins = 10;
0594   } else if (plotType_ == 2) {
0595     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,
0596                        3.0,   5.0,   7.0,   9.0,   11.0,  13.0,  15.0,  17.0, 19.0, 21.0, 23.0};
0597     for (int i = 0; i < 23; ++i) {
0598       etas_.push_back(xbin[i]);
0599       xbins[i] = xbin[i];
0600     }
0601     nbins = 22;
0602   } else {
0603     double xbina[99];
0604     int neta = numb_ / 2;
0605     for (int k = 0; k < neta; ++k) {
0606       xbina[k] = (k - neta) - 0.5;
0607       xbina[numb_ - k] = (neta - k) + 0.5;
0608     }
0609     xbina[neta] = 0;
0610     for (int i = 0; i < numb_ + 1; ++i) {
0611       etas_.push_back(xbina[i]);
0612       xbins[i] = xbina[i];
0613       ++nbins;
0614     }
0615   }
0616   int ipbin[npbin] = {10, 20, 30, 40, 60, 100};
0617   for (unsigned int i = 0; i < npbin; ++i)
0618     ps_.push_back((double)(ipbin[i]));
0619   int npvtx[6] = {0, 7, 10, 13, 16, 100};
0620   for (int i = 0; i < 6; ++i)
0621     nvx_.push_back(npvtx[i]);
0622   double dl1s[9] = {0, 0.10, 0.20, 0.50, 1.0, 2.0, 2.5, 3.0, 10.0};
0623   int ietasL[4] = {0, 13, 17, 0};
0624   int ietasH[4] = {14, 18, 25, 25};
0625   for (int i = 0; i < 4; ++i) {
0626     ietasL_.push_back(ietasL[i]);
0627     ietasH_.push_back(ietasH[i]);
0628   }
0629   int nxbin(100);
0630   double xlow(0.0), xhigh(5.0);
0631   if (coarseBin_ == 1) {
0632     xlow = 0.25;
0633     xhigh = 5.25;
0634     nxbin = 50;
0635   } else if (coarseBin_ > 1) {
0636     if (coarseBin_ == 2)
0637       nxbin = 500;
0638     else
0639       nxbin = 1000;
0640   }
0641
0642   char name[20], title[200];
0643   std::string titl[5] = {
0644       "All tracks", "Good quality tracks", "Selected tracks", "Tracks with charge isolation", "Tracks MIP in ECAL"};
0645   for (int i = 0; i < 9; ++i)
0646     dl1_.push_back(dl1s[i]);
0647   unsigned int kp = (ps_.size() - 1);
0648   for (unsigned int k = 0; k < kp; ++k) {
0649     for (unsigned int j = 0; j <= ietasL_.size(); ++j) {
0650       sprintf(name, "%spp%d%d", prefix_.c_str(), k, j);
0651       if (j == 0)
0652         sprintf(title, "E/p for %s (p = %d:%d GeV All)", titl[4].c_str(), ipbin[k], ipbin[k + 1]);
0653       else
0654         sprintf(title,
0655                 "E/p for %s (p = %d:%d GeV |#eta| %d:%d)",
0656                 titl[4].c_str(),
0657                 ipbin[k],
0658                 ipbin[k + 1],
0659                 (ietasL_[j - 1] + 1),
0660                 ietasH_[j - 1]);
0661       h_pp[k].push_back(new TH1D(name, title, 100, 10.0, 110.0));
0662       int kk = h_pp[k].size() - 1;
0663       h_pp[k][kk]->Sumw2();
0664     }
0665   }
0666   if (plotType_ <= 1) {
0667     std::cout << "Book Histos for Standard" << std::endl;
0668     for (unsigned int k = 0; k < kp; ++k) {
0669       for (unsigned int i = 0; i < nvx_.size(); ++i) {
0670         sprintf(name, "%setaX%d%d", prefix_.c_str(), k, i);
0671         if (i == 0) {
0672           sprintf(title, "%s (p = %d:%d GeV all vertices)", titl[4].c_str(), ipbin[k], ipbin[k + 1]);
0673         } else {
0674           sprintf(
0675               title, "%s (p = %d:%d GeV # Vtx %d:%d)", titl[4].c_str(), ipbin[k], ipbin[k + 1], nvx_[i - 1], nvx_[i]);
0676         }
0677         h_etaX[k].push_back(new TProfile(name, title, nbins, xbins));
0678         unsigned int kk = h_etaX[k].size() - 1;
0679         h_etaX[k][kk]->Sumw2();
0680         sprintf(name, "%snvxR%d%d", prefix_.c_str(), k, i);
0681         if (i == 0) {
0682           sprintf(title, "E/p for %s (p = %d:%d GeV all vertices)", titl[4].c_str(), ipbin[k], ipbin[k + 1]);
0683         } else {
0684           sprintf(title,
0685                   "E/p for %s (p = %d:%d GeV # Vtx %d:%d)",
0686                   titl[4].c_str(),
0687                   ipbin[k],
0688                   ipbin[k + 1],
0689                   nvx_[i - 1],
0690                   nvx_[i]);
0691         }
0692         h_nvxR[k].push_back(new TH1D(name, title, nxbin, xlow, xhigh));
0693         kk = h_nvxR[k].size() - 1;
0694         h_nvxR[k][kk]->Sumw2();
0695       }
0696       for (unsigned int j = 0; j < etas_.size(); ++j) {
0697         sprintf(name, "%sratio%d%d", prefix_.c_str(), k, j);
0698         if (j == 0) {
0699           sprintf(title, "E/p for %s (p = %d:%d GeV)", titl[4].c_str(), ipbin[k], ipbin[k + 1]);
0700         } else {
0701           sprintf(title,
0702                   "E/p for %s (p = %d:%d GeV #eta %4.1f:%4.1f)",
0703                   titl[4].c_str(),
0704                   ipbin[k],
0705                   ipbin[k + 1],
0706                   etas_[j - 1],
0707                   etas_[j]);
0708         }
0709         h_etaF[k].push_back(new TH1D(name, title, nxbin, xlow, xhigh));
0710         unsigned int kk = h_etaF[k].size() - 1;
0711         h_etaF[k][kk]->Sumw2();
0712         sprintf(name, "%setaR%d%d", prefix_.c_str(), k, j);
0713         h_etaR[k].push_back(new TH1D(name, title, nxbin, xlow, xhigh));
0714         kk = h_etaR[k].size() - 1;
0715         h_etaR[k][kk]->Sumw2();
0716       }
0717       for (unsigned int j = 1; j <= ietasL_.size(); ++j) {
0718         sprintf(name, "%setaB%d%d", prefix_.c_str(), k, j);
0719         sprintf(title,
0720                 "E/p for %s (p = %d:%d GeV |#eta| %d:%d)",
0721                 titl[4].c_str(),
0722                 ipbin[k],
0723                 ipbin[k + 1],
0724                 (ietasL_[j - 1] + 1),
0725                 ietasH_[j - 1]);
0726         h_etaB[k].push_back(new TH1D(name, title, nxbin, xlow, xhigh));
0727         unsigned int kk = h_etaB[k].size() - 1;
0728         h_etaB[k][kk]->Sumw2();
0729       }
0730       for (unsigned int j = 0; j < dl1_.size(); ++j) {
0731         sprintf(name, "%sdl1R%d%d", prefix_.c_str(), k, j);
0732         if (j == 0) {
0733           sprintf(title, "E/p for %s (p = %d:%d GeV All d_{L1})", titl[4].c_str(), ipbin[k], ipbin[k + 1]);
0734         } else {
0735           sprintf(title,
0736                   "E/p for %s (p = %d:%d GeV d_{L1} %4.2f:%4.2f)",
0737                   titl[4].c_str(),
0738                   ipbin[k],
0739                   ipbin[k + 1],
0740                   dl1_[j - 1],
0741                   dl1_[j]);
0742         }
0743         h_dL1R[k].push_back(new TH1D(name, title, nxbin, xlow, xhigh));
0744         unsigned int kk = h_dL1R[k].size() - 1;
0745         h_dL1R[k][kk]->Sumw2();
0746       }
0747     }
0748     for (unsigned int i = 0; i < nvx_.size(); ++i) {
0749       sprintf(name, "%setaX%d%d", prefix_.c_str(), kp, i);
0750       if (i == 0) {
0751         sprintf(title, "%s (All Momentum all vertices)", titl[4].c_str());
0752       } else {
0753         sprintf(title, "%s (All Momentum # Vtx %d:%d)", titl[4].c_str(), nvx_[i - 1], nvx_[i]);
0754       }
0755       h_etaX[npbin - 1].push_back(new TProfile(name, title, nbins, xbins));
0756       unsigned int kk = h_etaX[npbin - 1].size() - 1;
0757       h_etaX[npbin - 1][kk]->Sumw2();
0758       sprintf(name, "%snvxR%d%d", prefix_.c_str(), kp, i);
0759       if (i == 0) {
0760         sprintf(title, "E/p for %s (All Momentum all vertices)", titl[4].c_str());
0761       } else {
0762         sprintf(title, "E/p for %s (All Momentum # Vtx %d:%d)", titl[4].c_str(), nvx_[i - 1], nvx_[i]);
0763       }
0764       h_nvxR[npbin - 1].push_back(new TH1D(name, title, 200, 0., 10.));
0765       kk = h_nvxR[npbin - 1].size() - 1;
0766       h_nvxR[npbin - 1][kk]->Sumw2();
0767     }
0768     for (unsigned int j = 0; j < etas_.size(); ++j) {
0769       sprintf(name, "%sratio%d%d", prefix_.c_str(), kp, j);
0770       if (j == 0) {
0771         sprintf(title, "E/p for %s (All momentum)", titl[4].c_str());
0772       } else {
0773         sprintf(title, "E/p for %s (All momentum #eta %4.1f:%4.1f)", titl[4].c_str(), etas_[j - 1], etas_[j]);
0774       }
0775       h_etaF[npbin - 1].push_back(new TH1D(name, title, 200, 0., 10.));
0776       unsigned int kk = h_etaF[npbin - 1].size() - 1;
0777       h_etaF[npbin - 1][kk]->Sumw2();
0778       sprintf(name, "%setaR%d%d", prefix_.c_str(), kp, j);
0779       h_etaR[npbin - 1].push_back(new TH1D(name, title, 200, 0., 10.));
0780       kk = h_etaR[npbin - 1].size() - 1;
0781       h_etaR[npbin - 1][kk]->Sumw2();
0782     }
0783     for (unsigned int j = 1; j <= ietasL_.size(); ++j) {
0784       sprintf(name, "%setaB%d%d", prefix_.c_str(), kp, j);
0785       sprintf(title, "E/p for %s (All momentum |#eta| %d:%d)", titl[4].c_str(), (ietasL_[j - 1] + 1), ietasH_[j - 1]);
0786       h_etaB[npbin - 1].push_back(new TH1D(name, title, nxbin, xlow, xhigh));
0787       unsigned int kk = h_etaB[npbin - 1].size() - 1;
0788       h_etaB[npbin - 1][kk]->Sumw2();
0789     }
0790     for (unsigned int j = 0; j < dl1_.size(); ++j) {
0791       sprintf(name, "%sdl1R%d%d", prefix_.c_str(), kp, j);
0792       if (j == 0) {
0793         sprintf(title, "E/p for %s (All momentum)", titl[4].c_str());
0794       } else {
0795         sprintf(title, "E/p for %s (All momentum d_{L1} %4.2f:%4.2f)", titl[4].c_str(), dl1_[j - 1], dl1_[j]);
0796       }
0797       h_dL1R[npbin - 1].push_back(new TH1D(name, title, 200, 0., 10.));
0798       unsigned int kk = h_dL1R[npbin - 1].size() - 1;
0799       h_dL1R[npbin - 1][kk]->Sumw2();
0800     }
0801   } else {
0802     std::cout << "Book Histos for Non-Standard " << etas_.size() << ":" << kp50 << std::endl;
0803     unsigned int kp = (ps_.size() - 1);
0804     for (unsigned int k = 0; k < kp; ++k) {
0805       for (unsigned int j = 0; j < etas_.size(); ++j) {
0806         sprintf(name, "%sratio%d%d", prefix_.c_str(), k, j);
0807         if (j == 0) {
0808           sprintf(title, "E/p for %s (p = %d:%d GeV)", titl[4].c_str(), ipbin[k], ipbin[k + 1]);
0809         } else {
0810           sprintf(title,
0811                   "E/p for %s (p = %d:%d GeV #eta %4.1f:%4.1f)",
0812                   titl[4].c_str(),
0813                   ipbin[k],
0814                   ipbin[k + 1],
0815                   etas_[j - 1],
0816                   etas_[j]);
0817         }
0818         h_etaF[k].push_back(new TH1D(name, title, nxbin, xlow, xhigh));
0819         unsigned int kk = h_etaF[k].size() - 1;
0820         h_etaF[k][kk]->Sumw2();
0821       }
0822       for (unsigned int j = 1; j <= ietasL_.size(); ++j) {
0823         sprintf(name, "%setaB%d%d", prefix_.c_str(), k, j);
0824         sprintf(title,
0825                 "E/p for %s (p = %d:%d GeV |#eta| %d:%d)",
0826                 titl[4].c_str(),
0827                 ipbin[k],
0828                 ipbin[k + 1],
0829                 (ietasL_[j - 1] + 1),
0830                 ietasH_[j - 1]);
0831         h_etaB[k].push_back(new TH1D(name, title, nxbin, xlow, xhigh));
0832         unsigned int kk = h_etaB[k].size() - 1;
0833         h_etaB[k][kk]->Sumw2();
0834       }
0835     }
0836   }
0837   if (selRBX_) {
0838     for (unsigned int j = 1; j <= 18; ++j) {
0839       sprintf(name, "%sRBX%d%d", prefix_.c_str(), kp50, j);
0840       sprintf(title, "E/p for RBX%d (p = %d:%d GeV |#eta| %d:%d)", j, ipbin[kp50], ipbin[kp50 + 1], etalo_, etahi_);
0841       h_rbx.push_back(new TH1D(name, title, nxbin, xlow, xhigh));
0842       h_rbx[j - 1]->Sumw2();
0843     }
0844   }
0845 }
0846
0847 Bool_t CalibMonitor::Notify() {
0848   // The Notify() function is called when a new file is opened. This
0849   // can be either for a new TTree in a TChain or when when a new TTree
0850   // is started when using PROOF. It is normally not necessary to make changes
0851   // to the generated code, but the routine can be extended by the
0852   // user if needed. The return value is currently not used.
0853
0854   return kTRUE;
0855 }
0856
0857 void CalibMonitor::Show(Long64_t entry) {
0858   // Print contents of entry.
0859   // If entry is not specified, print current entry
0860   if (!fChain)
0861     return;
0862   fChain->Show(entry);
0863 }
0864
0865 Int_t CalibMonitor::Cut(Long64_t) {
0866   // This function may be called from Loop.
0867   // returns  1 if entry is accepted.
0868   // returns -1 otherwise.
0869   return 1;
0870 }
0871
0872 void CalibMonitor::Loop() {
0873   //   In a ROOT session, you can do:
0874   //      Root > .L CalibMonitor.C
0875   //      Root > CalibMonitor t
0876   //      Root > t.GetEntry(12); // Fill t data members with entry number 12
0877   //      Root > t.Show();       // Show values of entry 12
0878   //      Root > t.Show(16);     // Read and show values of entry 16
0879   //      Root > t.Loop();       // Loop on all entries
0880   //
0881
0882   //   This is the loop skeleton where:
0883   //      jentry is the global entry number in the chain
0884   //      ientry is the entry number in the current Tree
0885   //   Note that the argument to GetEntry must be:
0886   //      jentry for TChain::GetEntry
0887   //      ientry for TTree::GetEntry and TBranch::GetEntry
0888   //
0889   //       To read only selected branches, Insert statements like:
0890   // METHOD1:
0891   //    fChain->SetBranchStatus("*",0);  // disable all branches
0892   //    fChain->SetBranchStatus("branchname",1);  // activate branchname
0893   // METHOD2: replace line
0894   //    fChain->GetEntry(jentry);       //read all branches
0895   //by  b_branchname->GetEntry(ientry); //read only this branch
0896   const bool debug(false);
0897
0898   if (fChain == 0)
0899     return;
0900   std::map<int, counter> runSum, runEn1, runEn2;
0901
0902   // Find list of duplicate events
0903   Long64_t nentries = fChain->GetEntriesFast();
0904   std::cout << "Total entries " << nentries << std::endl;
0905   Long64_t nbytes(0), nb(0);
0906   unsigned int duplicate(0), good(0), kount(0);
0907   unsigned int kp1 = ps_.size() - 1;
0908   unsigned int kv1 = 0;
0909   std::vector<int> kounts(kp1, 0);
0910   std::vector<int> kount50(20, 0);
0911   std::vector<int> kount0(20, 0);
0912   std::vector<int> kount1(20, 0);
0913   std::vector<int> kount2(20, 0);
0914   std::vector<int> kount3(20, 0);
0915   std::vector<int> kount4(20, 0);
0916   std::vector<int> kount5(20, 0);
0917   for (Long64_t jentry = 0; jentry < nentries; jentry++) {
0918     //for (Long64_t jentry=0; jentry<200;jentry++) {
0919     Long64_t ientry = LoadTree(jentry);
0920     if (ientry < 0)
0921       break;
0922     nb = fChain->GetEntry(jentry);
0923     nbytes += nb;
0924     if (jentry % 1000000 == 0)
0925       std::cout << "Entry " << jentry << " Run " << t_Run << " Event " << t_Event << std::endl;
0926     double pmom = (useGen_ && (t_gentrackP > 0)) ? t_gentrackP : t_p;
0927     int kp(-1);
0928     for (unsigned int k = 1; k < ps_.size(); ++k) {
0929       if (pmom >= ps_[k - 1] && pmom < ps_[k]) {
0930         kp = k - 1;
0931         break;
0932       }
0933     }
0934     bool p4060 = ((pmom >= 40.0) && (pmom <= 60.0));
0935     if (p4060)
0936       ++kount50[0];
0937     if (kp == 0) {
0938       ++kount0[0];
0939     } else if (kp == 1) {
0940       ++kount1[0];
0941     } else if (kp == 2) {
0942       ++kount2[0];
0943     } else if (kp == 3) {
0944       ++kount3[0];
0945     } else if (kp == 4) {
0946       ++kount4[0];
0947     } else if (kp == 5) {
0948       ++kount5[0];
0949     }
0950     bool select = (std::find(entries_.begin(), entries_.end(), jentry) == entries_.end());
0951     if (!select) {
0952       ++duplicate;
0953       if (debug)
0954         std::cout << "Duplicate event " << t_Run << " " << t_Event << " " << t_p << std::endl;
0955       continue;
0956     }
0957     if (p4060)
0958       ++kount50[1];
0959     if (kp == 0) {
0960       ++kount0[1];
0961     } else if (kp == 1) {
0962       ++kount1[1];
0963     } else if (kp == 2) {
0964       ++kount2[1];
0965     } else if (kp == 3) {
0966       ++kount3[1];
0967     } else if (kp == 4) {
0968       ++kount4[1];
0969     } else if (kp == 5) {
0970       ++kount5[1];
0971     }
0972     bool selRun = (includeRun_ ? ((t_Run >= runlo_) && (t_Run <= runhi_)) : ((t_Run < runlo_) || (t_Run > runhi_)));
0973     if (select) {
0974       if (p4060)
0975         ++kount50[2];
0976       if (kp == 0) {
0977         ++kount0[2];
0978       } else if (kp == 1) {
0979         ++kount1[2];
0980       } else if (kp == 2) {
0981         ++kount2[2];
0982       } else if (kp == 3) {
0983         ++kount3[2];
0984       } else if (kp == 4) {
0985         ++kount4[2];
0986       } else if (kp == 5) {
0987         ++kount5[2];
0988       }
0989     }
0990     select =
0991         (selRun && (fabs(t_ieta) >= etalo_) && (fabs(t_ieta) <= etahi_) && (t_nVtx >= nvxlo_) && (t_nVtx <= nvxhi_));
0992     if (select) {
0993       if (p4060)
0994         ++kount50[3];
0995       if (kp == 0) {
0996         ++kount0[3];
0997       } else if (kp == 1) {
0998         ++kount1[3];
0999       } else if (kp == 2) {
1000         ++kount2[3];
1001       } else if (kp == 3) {
1002         ++kount3[3];
1003       } else if (kp == 4) {
1004         ++kount4[3];
1005       } else if (kp == 5) {
1006         ++kount5[3];
1007       }
1008     }
1009     if (!select) {
1010       if (debug)
1011         std::cout << "Run # " << t_Run << " out of range of " << runlo_ << ":" << runhi_ << " or ieta " << t_ieta
1012                   << " (" << etalo_ << ":" << etahi_ << ") or nvtx " << t_nVtx << " (" << nvxlo_ << ":" << nvxhi_
1013                   << ") out of range" << std::endl;
1014       continue;
1015     }
1016     if (cSelect_ != nullptr) {
1017       if (exclude_) {
1018         if (cSelect_->isItRBX(t_DetIds))
1019           continue;
1020       } else {
1021         if (!(cSelect_->isItRBX(t_ieta, t_iphi)))
1022           continue;
1023       }
1024     }
1025     if (p4060)
1026       ++kount50[4];
1027     if (kp == 0) {
1028       ++kount0[4];
1029     } else if (kp == 1) {
1030       ++kount1[4];
1031     } else if (kp == 2) {
1032       ++kount2[4];
1033     } else if (kp == 3) {
1034       ++kount3[4];
1035     } else if (kp == 4) {
1036       ++kount4[4];
1037     } else if (kp == 5) {
1038       ++kount5[4];
1039     }
1040     select = (!cutL1T_ || (t_mindR1 >= 0.5));
1041     if (!select) {
1042       if (debug)
1043         std::cout << "Reject Run # " << t_Run << " Event # " << t_Event << " too close to L1 trigger " << t_mindR1
1044                   << std::endl;
1045       continue;
1046     }
1047     if (p4060)
1048       ++kount50[5];
1049     if (kp == 0) {
1050       ++kount0[5];
1051     } else if (kp == 1) {
1052       ++kount1[5];
1053     } else if (kp == 2) {
1054       ++kount2[5];
1055     } else if (kp == 3) {
1056       ++kount3[5];
1057     } else if (kp == 4) {
1058       ++kount4[5];
1059     } else if (kp == 5) {
1060       ++kount5[5];
1061     }
1062     select = ((events_.size() == 0) ||
1063               (std::find(events_.begin(), events_.end(), std::pair<int, int>(t_Run, t_Event)) != events_.end()));
1064     if (!select) {
1065       if (debug)
1066         std::cout << "Reject Run # " << t_Run << " Event # " << t_Event << " not in the selection list" << std::endl;
1067       continue;
1068     }
1069     if (p4060)
1070       ++kount50[6];
1071     if (kp == 0) {
1072       ++kount0[6];
1073     } else if (kp == 1) {
1074       ++kount1[6];
1075     } else if (kp == 2) {
1076       ++kount2[6];
1077     } else if (kp == 3) {
1078       ++kount3[6];
1079     } else if (kp == 4) {
1080       ++kount4[6];
1081     } else if (kp == 5) {
1082       ++kount5[6];
1083     }
1084     if ((ifDepth_ == 3) && (cFactor_ != nullptr) && (cFactor_->absent(ientry)))
1085       continue;
1086
1087     // if (Cut(ientry) < 0) continue;
1088     int jp(-1), jp1(-1), jp2(-1);
1089     unsigned int kv = nvx_.size() - 1;
1090     for (unsigned int k = 1; k < nvx_.size(); ++k) {
1091       if (t_goodPV >= nvx_[k - 1] && t_goodPV < nvx_[k]) {
1092         kv = k;
1093         break;
1094       }
1095     }
1096     unsigned int kd1 = 0;
1097     unsigned int kd = dl1_.size() - 1;
1098     for (unsigned int k = 1; k < dl1_.size(); ++k) {
1099       if (t_mindR1 >= dl1_[k - 1] && t_mindR1 < dl1_[k]) {
1100         kd = k;
1101         break;
1102       }
1103     }
1104     double eta = (t_ieta > 0) ? ((double)(t_ieta)-0.001) : ((double)(t_ieta) + 0.001);
1105     for (unsigned int j = 1; j < etas_.size(); ++j) {
1106       if (eta > etas_[j - 1] && eta < etas_[j]) {
1107         jp = j;
1108         break;
1109       }
1110     }
1111     for (unsigned int j = 0; j < (ietasL_.size() - 1); ++j) {
1112       if (std::abs(t_ieta) > ietasL_[j] && std::abs(t_ieta) <= ietasH_[j]) {
1113         if (jp1 >= 0)
1114           jp2 = j;
1115         if (jp1 < 0)
1116           jp1 = j;
1117       }
1118     }
1119     int jp3 = (jp1 >= 0) ? (int)(ietasL_.size() - 1) : jp1;
1120     if (debug)
1121       std::cout << "Bin " << kp << ":" << kp1 << ":" << kv << ":" << kv1 << ":" << kd << ":" << kd1 << ":" << jp << ":"
1122                 << jp1 << ":" << jp2 << ":" << jp3 << " pmom:ieta:pv:mindR " << pmom << ":" << std::abs(t_ieta) << ":"
1123                 << t_goodPV << ":" << t_mindR1 << std::endl;
1124     double cut = (pmom > 20) ? ((flexibleSelect_ == 0) ? 2.0 : 10.0) : 0.0;
1125     //  double rcut= (pmom > 20) ? 0.25: 0.1;
1126     double rcut(-1000.0);
1127
1128     // Selection of good track and energy measured in Hcal
1129     double rat(1.0), eHcal(t_eHcal);
1130     if (corrFactor_->doCorr() || (cFactor_ != nullptr)) {
1131       eHcal = 0;
1132       for (unsigned int k = 0; k < t_HitEnergies->size(); ++k) {
1133         // The masks are defined in DataFormats/HcalDetId/interface/HcalDetId.h
1134         double cfac(1.0);
1135         if (corrFactor_->doCorr()) {
1136           unsigned int id = truncateId((*t_DetIds)[k], truncateFlag_, false);
1137           cfac = corrFactor_->getCorr(id);
1138         }
1139         if ((cFactor_ != nullptr) && (ifDepth_ != 3))
1140           cfac *= cFactor_->getCorr(t_Run, (*t_DetIds)[k]);
1141         eHcal += (cfac * ((*t_HitEnergies)[k]));
1142         if (debug) {
1143           int subdet, zside, ieta, iphi, depth;
1144           unpackDetId((*t_DetIds)[k], subdet, zside, ieta, iphi, depth);
1145           std::cout << zside << ":" << ieta << ":" << depth << " Corr " << cfac << " " << (*t_HitEnergies)[k] << " Out "
1146                     << eHcal << std::endl;
1147         }
1148       }
1149     }
1150     bool goodTk = goodTrack(eHcal, cut, jentry, debug);
1151     bool selPhi = selectPhi(debug);
1152     if (t_qltyFlag) {
1153       if (p4060)
1154         ++kount50[7];
1155       if (kp == 0) {
1156         ++kount0[7];
1157       } else if (kp == 1) {
1158         ++kount1[7];
1159       } else if (kp == 2) {
1160         ++kount2[7];
1161       } else if (kp == 3) {
1162         ++kount3[7];
1163       } else if (kp == 4) {
1164         ++kount4[7];
1165       } else if (kp == 5) {
1166         ++kount5[7];
1167       }
1168       if (t_selectTk) {
1169         if (p4060)
1170           ++kount50[8];
1171         if (kp == 0) {
1172           ++kount0[8];
1173         } else if (kp == 1) {
1174           ++kount1[8];
1175         } else if (kp == 2) {
1176           ++kount2[8];
1177         } else if (kp == 3) {
1178           ++kount3[8];
1179         } else if (kp == 4) {
1180           ++kount4[8];
1181         } else if (kp == 5) {
1182           ++kount5[8];
1183         }
1184         if (t_hmaxNearP < cut) {
1185           if (p4060)
1186             ++kount50[9];
1187           if (kp == 0) {
1188             ++kount0[9];
1189           } else if (kp == 1) {
1190             ++kount1[9];
1191           } else if (kp == 2) {
1192             ++kount2[9];
1193           } else if (kp == 3) {
1194             ++kount3[9];
1195           } else if (kp == 4) {
1196             ++kount4[9];
1197           } else if (kp == 5) {
1198             ++kount5[9];
1199           }
1200           if (t_eMipDR < 1.0) {
1201             if (p4060)
1202               ++kount50[10];
1203             if (kp == 0) {
1204               ++kount0[10];
1205             } else if (kp == 1) {
1206               ++kount1[10];
1207             } else if (kp == 2) {
1208               ++kount2[10];
1209             } else if (kp == 3) {
1210               ++kount3[10];
1211             } else if (kp == 4) {
1212               ++kount4[10];
1213             } else if (kp == 5) {
1214               ++kount5[10];
1215             }
1216             if (eHcal > 0.001) {
1217               if (p4060)
1218                 ++kount50[11];
1219               if (kp == 0) {
1220                 ++kount0[11];
1221               } else if (kp == 1) {
1222                 ++kount1[11];
1223               } else if (kp == 2) {
1224                 ++kount2[11];
1225               } else if (kp == 3) {
1226                 ++kount3[11];
1227               } else if (kp == 4) {
1228                 ++kount4[11];
1229               } else if (kp == 5) {
1230                 ++kount5[11];
1231               }
1232               if (selPhi) {
1233                 if (p4060)
1234                   ++kount50[12];
1235                 if (kp == 0) {
1236                   ++kount0[12];
1237                 } else if (kp == 1) {
1238                   ++kount1[12];
1239                 } else if (kp == 2) {
1240                   ++kount2[12];
1241                 } else if (kp == 3) {
1242                   ++kount3[12];
1243                 } else if (kp == 4) {
1244                   ++kount4[12];
1245                 } else if (kp == 5) {
1246                   ++kount5[12];
1247                 }
1248               }
1249             }
1250           }
1251         }
1252       }
1253     }
1254     if (pmom > 0)
1255       rat = (eHcal / (pmom - t_eMipDR));
1256     if (debug) {
1257       std::cout << "Entry " << jentry << " p|eHcal|ratio " << pmom << "|" << t_eHcal << "|" << eHcal << "|" << rat
1258                 << "|" << kp << "|" << kv << "|" << jp << " Cuts " << t_qltyFlag << "|" << t_selectTk << "|"
1259                 << (t_hmaxNearP < cut) << "|" << (t_eMipDR < 1.0) << "|" << goodTk << "|" << (rat > rcut)
1260                 << " Select Phi " << selPhi << std::endl;
1261       std::cout << "D1 : " << kp << ":" << kp1 << ":" << kv << ":" << kv1 << ":" << kd << ":" << kd1 << ":" << jp
1262                 << std::endl;
1263     }
1264     if (goodTk && (kp >= 0) && selPhi) {
1265       if (p4060)
1266         ++kount50[13];
1267       if (kp == 0) {
1268         ++kount0[13];
1269       } else if (kp == 1) {
1270         ++kount1[13];
1271       } else if (kp == 2) {
1272         ++kount2[13];
1273       } else if (kp == 3) {
1274         ++kount3[13];
1275       } else if (kp == 4) {
1276         ++kount4[13];
1277       } else if (kp == 5) {
1278         ++kount5[13];
1279       }
1280       if (t_eHcal < 0.01) {
1281         std::map<int, counter>::const_iterator itr = runEn1.find(t_Run);
1282         if (itr == runEn1.end()) {
1283           counter knt;
1284           if ((kp >= 0) && (kp < counter::npsize))
1285             knt.count[kp] = 1;
1286           knt.total = 1;
1287           runEn1[t_Run] = knt;
1288         } else {
1289           counter knt = runEn1[t_Run];
1290           if ((kp >= 0) && (kp < counter::npsize))
1291             ++knt.count[kp];
1292           ++knt.total;
1293           runEn1[t_Run] = knt;
1294         }
1295       }
1296       if (t_eMipDR < 0.01 && t_eHcal < 0.01) {
1297         if (p4060)
1298           ++kount50[14];
1299         if (kp == 0) {
1300           ++kount0[14];
1301         } else if (kp == 1) {
1302           ++kount1[14];
1303         } else if (kp == 2) {
1304           ++kount2[14];
1305         } else if (kp == 3) {
1306           ++kount3[14];
1307         } else if (kp == 4) {
1308           ++kount4[14];
1309         } else if (kp == 5) {
1310           ++kount5[14];
1311         }
1312         std::map<int, counter>::const_iterator itr = runEn2.find(t_Run);
1313         if (itr == runEn2.end()) {
1314           counter knt;
1315           if ((kp >= 0) && (kp < counter::npsize))
1316             knt.count[kp] = 1;
1317           knt.total = 1;
1318           runEn2[t_Run] = knt;
1319         } else {
1320           counter knt = runEn2[t_Run];
1321           if ((kp >= 0) && (kp < counter::npsize))
1322             ++knt.count[kp];
1323           ++knt.total;
1324           runEn2[t_Run] = knt;
1325         }
1326       }
1327       if (rat > rcut) {
1328         if (p4060)
1329           ++kount50[15];
1330         if (kp == 0) {
1331           ++kount0[15];
1332         } else if (kp == 1) {
1333           ++kount1[15];
1334         } else if (kp == 2) {
1335           ++kount2[15];
1336         } else if (kp == 3) {
1337           ++kount3[15];
1338         } else if (kp == 4) {
1339           ++kount4[15];
1340         } else if (kp == 5) {
1341           ++kount5[15];
1342         }
1343         if (plotType_ <= 1) {
1344           h_etaX[kp][kv]->Fill(eta, rat, t_EventWeight);
1345           h_etaX[kp][kv1]->Fill(eta, rat, t_EventWeight);
1346           h_nvxR[kp][kv]->Fill(rat, t_EventWeight);
1347           h_nvxR[kp][kv1]->Fill(rat, t_EventWeight);
1348           h_dL1R[kp][kd]->Fill(rat, t_EventWeight);
1349           h_dL1R[kp][kd1]->Fill(rat, t_EventWeight);
1350           if (jp > 0)
1351             h_etaR[kp][jp]->Fill(rat, t_EventWeight);
1352           h_etaR[kp][0]->Fill(rat, t_EventWeight);
1353         }
1354         h_pp[kp][0]->Fill(pmom, t_EventWeight);
1355         if (jp1 >= 0) {
1356           h_pp[kp][jp1 + 1]->Fill(pmom, t_EventWeight);
1357           h_pp[kp][jp3 + 1]->Fill(pmom, t_EventWeight);
1358         }
1359         if (jp2 >= 0)
1360           h_pp[kp][jp2 + 1]->Fill(pmom, t_EventWeight);
1361         if (kp == (int)(kp50)) {
1362           std::map<int, counter>::const_iterator itr = runSum.find(t_Run);
1363           if (itr == runSum.end()) {
1364             counter knt;
1365             if ((kp >= 0) && (kp < counter::npsize))
1366               knt.count[kp] = 1;
1367             knt.total = 1;
1368             runSum[t_Run] = knt;
1369           } else {
1370             counter knt = runSum[t_Run];
1371             if ((kp >= 0) && (kp < counter::npsize))
1372               ++knt.count[kp];
1373             ++knt.total;
1374             runSum[t_Run] = knt;
1375           }
1376         }
1377         if ((!dataMC_) || (t_mindR1 > 0.5) || (t_DataType == 1)) {
1378           if (p4060)
1379             ++kount50[16];
1380           if (kp == 0) {
1381             ++kount0[16];
1382           } else if (kp == 1) {
1383             ++kount1[16];
1384           } else if (kp == 2) {
1385             ++kount2[16];
1386           } else if (kp == 3) {
1387             ++kount3[16];
1388           } else if (kp == 4) {
1389             ++kount4[16];
1390           } else if (kp == 5) {
1391             ++kount5[16];
1392           }
1393           ++kounts[kp];
1394           if (plotType_ <= 1) {
1395             if (jp > 0)
1396               h_etaF[kp][jp]->Fill(rat, t_EventWeight);
1397             h_etaF[kp][0]->Fill(rat, t_EventWeight);
1398           } else {
1399             if (debug) {
1400               std::cout << "kp " << kp << h_etaF[kp].size() << std::endl;
1401             }
1402             if (jp > 0)
1403               h_etaF[kp][jp]->Fill(rat, t_EventWeight);
1404             h_etaF[kp][0]->Fill(rat, t_EventWeight);
1405             if (jp1 >= 0) {
1406               h_etaB[kp][jp1]->Fill(rat, t_EventWeight);
1407               h_etaB[kp][jp3]->Fill(rat, t_EventWeight);
1408             }
1409             if (jp2 >= 0)
1410               h_etaB[kp][jp2]->Fill(rat, t_EventWeight);
1411           }
1412           if (selRBX_ && (kp == (int)(kp50)) && ((t_ieta * zside_) > 0)) {
1413             int rbx = (t_iphi > 70) ? 0 : (t_iphi + 1) / 4;
1414             h_rbx[rbx]->Fill(rat, t_EventWeight);
1415           }
1416         }
1417         if (pmom > 10.0) {
1418           if (plotType_ <= 1) {
1419             h_etaX[kp1][kv]->Fill(eta, rat, t_EventWeight);
1420             h_etaX[kp1][kv1]->Fill(eta, rat, t_EventWeight);
1421             h_nvxR[kp1][kv]->Fill(rat, t_EventWeight);
1422             h_nvxR[kp1][kv1]->Fill(rat, t_EventWeight);
1423             h_dL1R[kp1][kd]->Fill(rat, t_EventWeight);
1424             h_dL1R[kp1][kd1]->Fill(rat, t_EventWeight);
1425             if (jp > 0)
1426               h_etaR[kp1][jp]->Fill(rat, t_EventWeight);
1427             h_etaR[kp1][0]->Fill(rat, t_EventWeight);
1428             if (jp1 >= 0) {
1429               h_etaB[kp][jp1]->Fill(rat, t_EventWeight);
1430               h_etaB[kp][jp3]->Fill(rat, t_EventWeight);
1431             }
1432             if (jp2 >= 0)
1433               h_etaB[kp][jp2]->Fill(rat, t_EventWeight);
1434           }
1435           if (p4060)
1436             ++kount50[17];
1437           if (kp == 0) {
1438             ++kount0[17];
1439           } else if (kp == 1) {
1440             ++kount1[17];
1441           } else if (kp == 2) {
1442             ++kount2[17];
1443           } else if (kp == 3) {
1444             ++kount3[17];
1445           } else if (kp == 4) {
1446             ++kount4[17];
1447           } else if (kp == 5) {
1448             ++kount5[17];
1449           }
1450         }
1451       }
1452     }
1453     if (pmom > 10.0) {
1454       kount++;
1455       if (((flag_ / 100) % 10) != 0) {
1456         good++;
1457         fileout_ << good << " " << jentry << " " << t_Run << " " << t_Event << " " << t_ieta << " " << pmom
1458                  << std::endl;
1459       }
1460     }
1461   }
1462   unsigned int k(0);
1463   std::cout << "\nSummary of entries with " << runSum.size() << " runs\n";
1464   for (std::map<int, counter>::iterator itr = runSum.begin(); itr != runSum.end(); ++itr, ++k)
1465     std::cout << "[" << k << "] Run " << itr->first << " Total " << (itr->second).total << " in p-bins "
1466               << (itr->second).count[0] << ":" << (itr->second).count[1] << ":" << (itr->second).count[2] << ":"
1467               << (itr->second).count[3] << std::endl;
1468   k = 0;
1469   std::cout << runEn1.size() << " runs with 0 energy in HCAL\n";
1470   for (std::map<int, counter>::iterator itr = runEn1.begin(); itr != runEn1.end(); ++itr, ++k)
1471     std::cout << "[" << k << "] Run " << itr->first << " Total " << (itr->second).total << " in p-bins "
1472               << (itr->second).count[0] << ":" << (itr->second).count[1] << ":" << (itr->second).count[2] << ":"
1473               << (itr->second).count[3] << std::endl;
1474   k = 0;
1475   std::cout << runEn2.size() << " runs with 0 energy in ECAL and HCAL\n";
1476   for (std::map<int, counter>::iterator itr = runEn2.begin(); itr != runEn2.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   if (((flag_ / 100) % 10) > 0) {
1481     fileout_.close();
1482     std::cout << "Writes " << good << " events in the file " << outFileName_ << std::endl;
1483   }
1484   std::cout << "Finds " << duplicate << " Duplicate events out of " << kount << " events in this file with p>10 Gev"
1485             << std::endl;
1486   std::cout << "Number of selected events:" << std::endl;
1487   for (unsigned int k = 1; k < ps_.size(); ++k)
1488     std::cout << ps_[k - 1] << ":" << ps_[k] << "     " << kounts[k - 1] << std::endl;
1489   std::cout << "Number in each step for tracks of momentum 40-60 GeV: ";
1490   for (unsigned int k = 0; k < 18; ++k)
1491     std::cout << " [" << k << "] " << kount50[k];
1492   std::cout << std::endl;
1493   for (unsigned int k = 1; k < ps_.size(); ++k) {
1494     std::cout << "Number in each step for tracks of momentum " << ps_[k - 1] << "-" << ps_[k] << " Gev: ";
1495     for (unsigned int k1 = 0; k1 < 18; ++k1) {
1496       if (k == 1) {
1497         std::cout << " [" << k1 << "] " << kount0[k1];
1498       } else if (k == 2) {
1499         std::cout << " [" << k1 << "] " << kount1[k1];
1500       } else if (k == 3) {
1501         std::cout << " [" << k1 << "] " << kount2[k1];
1502       } else if (k == 4) {
1503         std::cout << " [" << k1 << "] " << kount3[k1];
1504       } else if (k == 5) {
1505         std::cout << " [" << k1 << "] " << kount4[k1];
1506       } else if (k == 6) {
1507         std::cout << " [" << k1 << "] " << kount5[k1];
1508       }
1509     }
1510     std::cout << std::endl;
1511   }
1512 }
1513
1514 bool CalibMonitor::goodTrack(double &eHcal, double &cuti, const Long64_t &entry, bool debug) {
1515   bool select(true);
1516   double cut(cuti);
1517   if (debug) {
1518     std::cout << "goodTrack input " << eHcal << ":" << cut;
1519   }
1520   if (flexibleSelect_ > 1) {
1521     double eta = (t_ieta > 0) ? t_ieta : -t_ieta;
1522     cut = 8.0 * exp(eta * log2by18_);
1523   }
1524   correctEnergy(eHcal, entry);
1525   select = ((t_qltyFlag) && (t_selectTk) && (t_hmaxNearP < cut) && (t_eMipDR < 1.0) && (eHcal > 0.001));
1526   if (debug) {
1527     std::cout << " output " << select << " Based on " << t_qltyFlag << ":" << t_selectTk << ":" << t_hmaxNearP << ":"
1528               << cut << ":" << t_eMipDR << ":" << eHcal << std::endl;
1529   }
1530   return select;
1531 }
1532
1533 bool CalibMonitor::selectPhi(bool debug) {
1534   bool select(true);
1535   if (phimin_ > 1 || phimax_ < 72) {
1536     double eTotal(0), eSelec(0);
1537     // The masks are defined in DataFormats/HcalDetId/interface/HcalDetId.h
1538     for (unsigned int k = 0; k < t_HitEnergies->size(); ++k) {
1539       int iphi = ((*t_DetIds)[k]) & (0x3FF);
1540       int zside = ((*t_DetIds)[k] & 0x80000) ? (1) : (-1);
1541       eTotal += ((*t_HitEnergies)[k]);
1542       if (iphi >= phimin_ && iphi <= phimax_ && zside == zside_)
1543         eSelec += ((*t_HitEnergies)[k]);
1544     }
1545     if (eSelec < 0.9 * eTotal)
1546       select = false;
1547     if (debug) {
1548       std::cout << "Etotal " << eTotal << " and ESelec " << eSelec << " (phi " << phimin_ << ":" << phimax_ << " z "
1549                 << zside_ << ") Selection " << select << std::endl;
1550     }
1551   }
1552   return select;
1553 }
1554
1555 void CalibMonitor::plotHist(int itype, int inum, bool save) {
1556   gStyle->SetCanvasBorderMode(0);
1557   gStyle->SetCanvasColor(kWhite);
1558   gStyle->SetPadColor(kWhite);
1559   gStyle->SetFillColor(kWhite);
1560   gStyle->SetOptStat(111110);
1561   gStyle->SetOptFit(1);
1562   char name[100];
1563   int itmin = (itype >= 0 && itype < 4) ? itype : 0;
1564   int itmax = (itype >= 0 && itype < 4) ? itype : 3;
1565   std::string types[4] = {
1566       "E_{HCAL}/(p-E_{ECAL})", "E_{HCAL}/(p-E_{ECAL})", "E_{HCAL}/(p-E_{ECAL})", "E_{HCAL}/(p-E_{ECAL})"};
1567   int nmax[4] = {npbin, npbin, npbin, npbin};
1568   for (int type = itmin; type <= itmax; ++type) {
1569     int inmin = (inum >= 0 && inum < nmax[type]) ? inum : 0;
1570     int inmax = (inum >= 0 && inum < nmax[type]) ? inum : nmax[type] - 1;
1571     int kmax = 1;
1572     if (type == 0)
1573       kmax = (int)(etas_.size());
1574     else if (type == 1)
1575       kmax = (int)(etas_.size());
1576     else if (type == 2)
1577       kmax = (int)(nvx_.size());
1578     else
1579       kmax = (int)(dl1_.size());
1580     for (int num = inmin; num <= inmax; ++num) {
1581       for (int k = 0; k < kmax; ++k) {
1582         sprintf(name, "c_%d%d%d", type, num, k);
1583         TCanvas *pad = new TCanvas(name, name, 700, 500);
1584         pad->SetRightMargin(0.10);
1585         pad->SetTopMargin(0.10);
1586         sprintf(name, "%s", types[type].c_str());
1587         if (type != 7) {
1588           TH1D *hist(0);
1589           if (type == 0)
1590             hist = (TH1D *)(h_etaR[num][k]->Clone());
1591           else if (type == 1)
1592             hist = (TH1D *)(h_etaF[num][k]->Clone());
1593           else if (type == 2)
1594             hist = (TH1D *)(h_nvxR[num][k]->Clone());
1595           else
1596             hist = (TH1D *)(h_dL1R[num][k]->Clone());
1597           hist->GetXaxis()->SetTitle(name);
1598           hist->GetYaxis()->SetTitle("Tracks");
1599           drawHist(hist, pad);
1600           if (save) {
1601             sprintf(name, "c_%s%d%d%d.gif", prefix_.c_str(), type, num, k);
1602             pad->Print(name);
1603           }
1604         } else {
1605           TProfile *hist = (TProfile *)(h_etaX[num][k]->Clone());
1606           hist->GetXaxis()->SetTitle(name);
1607           hist->GetYaxis()->SetTitle("<E_{HCAL}/(p-E_{ECAL})>");
1608           hist->GetYaxis()->SetRangeUser(0.4, 1.6);
1609           hist->Fit("pol0", "q");
1610           drawHist(hist, pad);
1611           if (save) {
1612             sprintf(name, "c_%s%d%d%d.gif", prefix_.c_str(), type, num, k);
1613             pad->Print(name);
1614           }
1615         }
1616       }
1617     }
1618   }
1619 }
1620
1621 template <class Hist>
1622 void CalibMonitor::drawHist(Hist *hist, TCanvas *pad) {
1623   hist->GetYaxis()->SetLabelOffset(0.005);
1624   hist->GetYaxis()->SetTitleOffset(1.20);
1625   hist->Draw();
1626   pad->Update();
1627   TPaveStats *st1 = (TPaveStats *)hist->GetListOfFunctions()->FindObject("stats");
1628   if (st1 != NULL) {
1629     st1->SetY1NDC(0.70);
1630     st1->SetY2NDC(0.90);
1631     st1->SetX1NDC(0.55);
1632     st1->SetX2NDC(0.90);
1633   }
1634   pad->Modified();
1635   pad->Update();
1636 }
1637
1638 void CalibMonitor::savePlot(const std::string &theName, bool append, bool all) {
1639   TFile *theFile(0);
1640   if (append) {
1641     theFile = new TFile(theName.c_str(), "UPDATE");
1642   } else {
1643     theFile = new TFile(theName.c_str(), "RECREATE");
1644   }
1645
1646   theFile->cd();
1647   for (unsigned int k = 0; k < ps_.size(); ++k) {
1648     for (unsigned int j = 0; j <= ietasL_.size(); ++j) {
1649       if ((all || k == kp50) && h_pp[k].size() > j && h_pp[k][j] != 0) {
1650         TH1D *hist = (TH1D *)h_pp[k][j]->Clone();
1651         hist->Write();
1652       }
1653     }
1654     if (plotType_ <= 1) {
1655       for (unsigned int i = 0; i < nvx_.size(); ++i) {
1656         if (h_etaX[k][i] != 0) {
1657           TProfile *hnew = (TProfile *)h_etaX[k][i]->Clone();
1658           hnew->Write();
1659         }
1660         if (h_nvxR[k].size() > i && h_nvxR[k][i] != 0) {
1661           TH1D *hist = (TH1D *)h_nvxR[k][i]->Clone();
1662           hist->Write();
1663         }
1664       }
1665     }
1666     for (unsigned int j = 0; j < etas_.size(); ++j) {
1667       if ((plotType_ <= 1) && (h_etaR[k][j] != 0)) {
1668         TH1D *hist = (TH1D *)h_etaR[k][j]->Clone();
1669         hist->Write();
1670       }
1671       if (h_etaF[k].size() > j && h_etaF[k][j] != 0 && (all || (k == kp50))) {
1672         TH1D *hist = (TH1D *)h_etaF[k][j]->Clone();
1673         hist->Write();
1674       }
1675     }
1676     if (plotType_ <= 1) {
1677       for (unsigned int j = 0; j < dl1_.size(); ++j) {
1678         if (h_dL1R[k][j] != 0) {
1679           TH1D *hist = (TH1D *)h_dL1R[k][j]->Clone();
1680           hist->Write();
1681         }
1682       }
1683     }
1684     for (unsigned int j = 0; j < ietasL_.size(); ++j) {
1685       if (h_etaB[k].size() > j && h_etaB[k][j] != 0 && (all || (k == kp50))) {
1686         TH1D *hist = (TH1D *)h_etaB[k][j]->Clone();
1687         hist->Write();
1688       }
1689     }
1690   }
1691   if (selRBX_) {
1692     for (unsigned int k = 0; k < 18; ++k) {
1693       if (h_rbx[k] != 0) {
1694         TH1D *h1 = (TH1D *)h_rbx[k]->Clone();
1695         h1->Write();
1696       }
1697     }
1698   }
1699   std::cout << "All done" << std::endl;
1700   theFile->Close();
1701 }
1702
1703 void CalibMonitor::correctEnergy(double &eHcal, const Long64_t &entry) {
1704   bool debug(false);
1705   double pmom = (useGen_ && (t_gentrackP > 0)) ? t_gentrackP : t_p;
1706   if ((ifDepth_ == 3) && (cFactor_ != nullptr)) {
1707     double cfac = cFactor_->getCorr(entry);
1708     eHcal *= cfac;
1709     if (debug)
1710       std::cout << "PU Factor for " << ifDepth_ << ":" << entry << " = " << cfac << ":" << eHcal << std::endl;
1711   } else if ((corrPU_ < 0) && (pmom > 0)) {
1712     double ediff = (t_eHcal30 - t_eHcal10);
1713     if (t_DetIds1 != 0 && t_DetIds3 != 0) {
1714       double Etot1(0), Etot3(0);
1715       // The masks are defined in DataFormats/HcalDetId/interface/HcalDetId.h
1716       for (unsigned int idet = 0; idet < (*t_DetIds1).size(); idet++) {
1717         unsigned int id = truncateId((*t_DetIds1)[idet], truncateFlag_, false);
1718         double cfac = corrFactor_->getCorr(id);
1719         if (cFactor_ != 0)
1720           cfac *= cFactor_->getCorr(t_Run, (*t_DetIds1)[idet]);
1721         double hitEn = cfac * (*t_HitEnergies1)[idet];
1722         Etot1 += hitEn;
1723       }
1724       for (unsigned int idet = 0; idet < (*t_DetIds3).size(); idet++) {
1725         unsigned int id = truncateId((*t_DetIds3)[idet], truncateFlag_, false);
1726         double cfac = corrFactor_->getCorr(id);
1727         if (cFactor_ != 0)
1728           cfac *= cFactor_->getCorr(t_Run, (*t_DetIds3)[idet]);
1729         double hitEn = cfac * (*t_HitEnergies3)[idet];
1730         Etot3 += hitEn;
1731       }
1732       ediff = (Etot3 - Etot1);
1733     }
1734     double fac = puFactor(-corrPU_, t_ieta, pmom, eHcal, ediff, false);
1735     if (debug) {
1736       double fac1 = puFactor(-corrPU_, t_ieta, pmom, eHcal, ediff, true);
1737       double fac2 = puFactor(2, t_ieta, pmom, eHcal, ediff, true);
1738       std::cout << "PU Factor for " << -corrPU_ << " = " << fac1 << "; for 2 = " << fac2 << std::endl;
1739     }
1740     eHcal *= fac;
1741   } else if (corrPU_ > 0) {
1742     eHcal = puFactorRho(corrPU_, t_ieta, t_rhoh, eHcal);
1743   }
1744 }
1745
1746 class GetEntries {
1747 public:
1748   TTree *fChain;   //!pointer to the analyzed TTree/TChain
1749   Int_t fCurrent;  //!current Tree number in a TChain
1750
1751   // Declaration of leaf types
1752   UInt_t t_RunNo;
1753   UInt_t t_EventNo;
1754   Int_t t_Tracks;
1755   Int_t t_TracksProp;
1756   Int_t t_TracksSaved;
1757   Int_t t_TracksLoose;
1758   Int_t t_TracksTight;
1759   Int_t t_allvertex;
1760   Bool_t t_TrigPass;
1761   Bool_t t_TrigPassSel;
1762   Bool_t t_L1Bit;
1763   std::vector<Bool_t> *t_hltbits;
1764   std::vector<int> *t_ietaAll;
1765   std::vector<int> *t_ietaGood;
1766   std::vector<int> *t_trackType;
1767
1768   // List of branches
1769   TBranch *b_t_RunNo;        //!
1770   TBranch *b_t_EventNo;      //!
1771   TBranch *b_t_Tracks;       //!
1772   TBranch *b_t_TracksProp;   //!
1773   TBranch *b_t_TracksSaved;  //!
1774   TBranch *b_t_TracksLoose;  //!
1775   TBranch *b_t_TracksTight;  //!
1776   TBranch *b_t_allvertex;    //!
1777   TBranch *b_t_TrigPass;     //!
1778   TBranch *b_t_TrigPassSel;  //!
1779   TBranch *b_t_L1Bit;        //!
1780   TBranch *b_t_hltbits;      //!
1781   TBranch *b_t_ietaAll;      //!
1782   TBranch *b_t_ietaGood;     //!
1783   TBranch *b_t_trackType;    //!
1784
1785   GetEntries(const std::string &fname,
1786              const std::string &dirname,
1787              const char *dupFileName,
1788              const unsigned int bit1,
1789              const unsigned int bit2);
1790   virtual ~GetEntries();
1791   virtual Int_t Cut(Long64_t entry);
1792   virtual Int_t GetEntry(Long64_t entry);
1793   virtual Long64_t LoadTree(Long64_t entry);
1794   virtual void Init(TTree *tree, const char *dupFileName);
1795   virtual void Loop();
1796   virtual Bool_t Notify();
1797   virtual void Show(Long64_t entry = -1);
1798
1799 private:
1800   unsigned int bit_[2];
1801   std::vector<Long64_t> entries_;
1802   TH1I *h_tk[3], *h_eta[4], *h_pvx[3];
1803   TH1D *h_eff[3];
1804 };
1805
1806 GetEntries::GetEntries(const std::string &fname,
1807                        const std::string &dirnm,
1808                        const char *dupFileName,
1809                        const unsigned int bit1,
1810                        const unsigned int bit2) {
1811   TFile *file = new TFile(fname.c_str());
1812   TDirectory *dir = (TDirectory *)file->FindObjectAny(dirnm.c_str());
1813   std::cout << fname << " file " << file << " " << dirnm << " " << dir << std::endl;
1814   TTree *tree = (TTree *)dir->Get("EventInfo");
1815   std::cout << "CalibTree " << tree << std::endl;
1816   bit_[0] = bit1;
1817   bit_[1] = bit2;
1818   Init(tree, dupFileName);
1819 }
1820
1821 GetEntries::~GetEntries() {
1822   if (!fChain)
1823     return;
1824   delete fChain->GetCurrentFile();
1825 }
1826
1827 Int_t GetEntries::GetEntry(Long64_t entry) {
1828   // Read contents of entry.
1829   if (!fChain)
1830     return 0;
1831   return fChain->GetEntry(entry);
1832 }
1833
1834 Long64_t GetEntries::LoadTree(Long64_t entry) {
1835   // Set the environment to read one entry
1836   if (!fChain)
1837     return -5;
1838   Long64_t centry = fChain->LoadTree(entry);
1839   if (centry < 0)
1840     return centry;
1841   if (!fChain->InheritsFrom(TChain::Class()))
1842     return centry;
1843   TChain *chain = (TChain *)fChain;
1844   if (chain->GetTreeNumber() != fCurrent) {
1845     fCurrent = chain->GetTreeNumber();
1846     Notify();
1847   }
1848   return centry;
1849 }
1850
1851 void GetEntries::Init(TTree *tree, const char *dupFileName) {
1852   // The Init() function is called when the selector needs to initialize
1853   // a new tree or chain. Typically here the branch addresses and branch
1854   // pointers of the tree will be set.
1855   // It is normally not necessary to make changes to the generated
1856   // code, but the routine can be extended by the user if needed.
1857   // Init() will be called many times when running on PROOF
1858   // (once per file to be processed).
1859
1860   // Set branch addresses and branch pointers
1861   // Set object pointer
1862   t_hltbits = 0;
1863   t_ietaAll = 0;
1864   t_ietaGood = 0;
1865   t_trackType = 0;
1866   t_L1Bit = false;
1867   fChain = tree;
1868   fCurrent = -1;
1869   if (!tree)
1870     return;
1871   fChain->SetMakeClass(1);
1872   fChain->SetBranchAddress("t_RunNo", &t_RunNo, &b_t_RunNo);
1873   fChain->SetBranchAddress("t_EventNo", &t_EventNo, &b_t_EventNo);
1874   fChain->SetBranchAddress("t_Tracks", &t_Tracks, &b_t_Tracks);
1875   fChain->SetBranchAddress("t_TracksProp", &t_TracksProp, &b_t_TracksProp);
1876   fChain->SetBranchAddress("t_TracksSaved", &t_TracksSaved, &b_t_TracksSaved);
1877   fChain->SetBranchAddress("t_TracksLoose", &t_TracksLoose, &b_t_TracksLoose);
1878   fChain->SetBranchAddress("t_TracksTight", &t_TracksTight, &b_t_TracksTight);
1879   fChain->SetBranchAddress("t_allvertex", &t_allvertex, &b_t_allvertex);
1880   fChain->SetBranchAddress("t_TrigPass", &t_TrigPass, &b_t_TrigPass);
1881   fChain->SetBranchAddress("t_TrigPassSel", &t_TrigPassSel, &b_t_TrigPassSel);
1882   fChain->SetBranchAddress("t_L1Bit", &t_L1Bit, &b_t_L1Bit);
1883   fChain->SetBranchAddress("t_hltbits", &t_hltbits, &b_t_hltbits);
1884   fChain->SetBranchAddress("t_ietaAll", &t_ietaAll, &b_t_ietaAll);
1885   fChain->SetBranchAddress("t_ietaGood", &t_ietaGood, &b_t_ietaGood);
1886   fChain->SetBranchAddress("t_trackType", &t_trackType, &b_t_trackType);
1887   Notify();
1888
1889   if (std::string(dupFileName) != "") {
1890     ifstream infile(dupFileName);
1891     if (!infile.is_open()) {
1892       std::cout << "Cannot open " << dupFileName << std::endl;
1893     } else {
1894       while (1) {
1895         Long64_t jentry;
1896         infile >> jentry;
1897         if (!infile.good())
1898           break;
1899         entries_.push_back(jentry);
1900       }
1901       infile.close();
1902       std::cout << "Reads a list of " << entries_.size() << " events from " << dupFileName << std::endl;
1903     }
1904   }
1905
1906   h_tk[0] = new TH1I("Track0", "# of tracks produced", 2000, 0, 2000);
1907   h_tk[1] = new TH1I("Track1", "# of tracks propagated", 2000, 0, 2000);
1908   h_tk[2] = new TH1I("Track2", "# of tracks saved in tree", 2000, 0, 2000);
1909   h_eta[0] = new TH1I("Eta0", "i#eta (All Tracks)", 60, -30, 30);
1910   h_eta[1] = new TH1I("Eta1", "i#eta (Good Tracks)", 60, -30, 30);
1911   h_eta[2] = new TH1I("Eta2", "i#eta (Loose Isolated Tracks)", 60, -30, 30);
1912   h_eta[3] = new TH1I("Eta3", "i#eta (Tight Isolated Tracks)", 60, -30, 30);
1913   h_eff[0] = new TH1D("Eff0", "i#eta (Selection Efficiency)", 60, -30, 30);
1914   h_eff[1] = new TH1D("Eff1", "i#eta (Loose Isolation Efficiency)", 60, -30, 30);
1915   h_eff[2] = new TH1D("Eff2", "i#eta (Tight Isolation Efficiency)", 60, -30, 30);
1916   h_pvx[0] = new TH1I("Pvx0", "Number of Good Vertex (All)", 100, 0, 100);
1917   h_pvx[1] = new TH1I("Pvx1", "Number of Good Vertex (Loose)", 100, 0, 100);
1918   h_pvx[2] = new TH1I("Pvx2", "Number of Good Vertex (Tight)", 100, 0, 100);
1919 }
1920
1921 Bool_t GetEntries::Notify() {
1922   // The Notify() function is called when a new file is opened. This
1923   // can be either for a new TTree in a TChain or when when a new TTree
1924   // is started when using PROOF. It is normally not necessary to make changes
1925   // to the generated code, but the routine can be extended by the
1926   // user if needed. The return value is currently not used.
1927
1928   return kTRUE;
1929 }
1930
1931 void GetEntries::Show(Long64_t entry) {
1932   // Print contents of entry.
1933   // If entry is not specified, print current entry
1934   if (!fChain)
1935     return;
1936   fChain->Show(entry);
1937 }
1938
1939 Int_t GetEntries::Cut(Long64_t) {
1940   // This function may be called from Loop.
1941   // returns  1 if entry is accepted.
1942   // returns -1 otherwise.
1943   return 1;
1944 }
1945
1946 void GetEntries::Loop() {
1947   //   In a ROOT session, you can do:
1948   //      Root > .L CalibMonitor.C+g
1949   //      Root > GetEntries t
1950   //      Root > t.GetEntry(12); // Fill t data members with entry number 12
1951   //      Root > t.Show();       // Show values of entry 12
1952   //      Root > t.Show(16);     // Read and show values of entry 16
1953   //      Root > t.Loop();       // Loop on all entries
1954   //
1955
1956   //     This is the loop skeleton where:
1957   //    jentry is the global entry number in the chain
1958   //    ientry is the entry number in the current Tree
1959   //  Note that the argument to GetEntry must be:
1960   //    jentry for TChain::GetEntry
1961   //    ientry for TTree::GetEntry and TBranch::GetEntry
1962   //
1963   //       To read only selected branches, Insert statements like:
1964   // METHOD1:
1965   //    fChain->SetBranchStatus("*",0);  // disable all branches
1966   //    fChain->SetBranchStatus("branchname",1);  // activate branchname
1967   // METHOD2: replace line
1968   //    fChain->GetEntry(jentry);       //read all branches
1969   //by  b_branchname->GetEntry(ientry); //read only this branch
1970   if (fChain == 0)
1971     return;
1972
1973   Long64_t nentries = fChain->GetEntriesFast();
1974   Long64_t nbytes = 0, nb = 0;
1975   unsigned int kount(0), duplicate(0), selected(0);
1976   int l1(0), hlt(0), loose(0), tight(0);
1977   int allHLT[3] = {0, 0, 0};
1978   int looseHLT[3] = {0, 0, 0};
1979   int tightHLT[3] = {0, 0, 0};
1980   for (Long64_t jentry = 0; jentry < nentries; jentry++) {
1981     Long64_t ientry = LoadTree(jentry);
1982     if (ientry < 0)
1983       break;
1984     nb = fChain->GetEntry(jentry);
1985     nbytes += nb;
1986     bool select = (std::find(entries_.begin(), entries_.end(), jentry) == entries_.end());
1987     if (!select) {
1988       ++duplicate;
1989       continue;
1990     }
1991     h_tk[0]->Fill(t_Tracks);
1992     h_tk[1]->Fill(t_TracksProp);
1993     h_tk[2]->Fill(t_TracksSaved);
1994     h_pvx[0]->Fill(t_allvertex);
1995     if (t_TracksLoose > 0)
1996       h_pvx[1]->Fill(t_allvertex);
1997     if (t_TracksTight > 0)
1998       h_pvx[2]->Fill(t_allvertex);
1999     if (t_L1Bit) {
2000       ++l1;
2001       if (t_TracksLoose > 0)
2002         ++loose;
2003       if (t_TracksTight > 0)
2004         ++tight;
2005       if (t_TrigPass)
2006         ++hlt;
2007     }
2008     if (t_TrigPass) {
2009       ++kount;
2010       if (t_TrigPassSel)
2011         ++selected;
2012     }
2013     bool passt[2] = {false, false}, pass(false);
2014     for (unsigned k = 0; k < t_hltbits->size(); ++k) {
2015       if ((*t_hltbits)[k] > 0) {
2016         pass = true;
2017         for (int i = 0; i < 2; ++i)
2018           if (k == bit_[i])
2019             passt[i] = true;
2020       }
2021     }
2022     if (pass) {
2023       ++allHLT[0];
2024       for (int i = 0; i < 2; ++i)
2025         if (passt[i])
2026           ++allHLT[i + 1];
2027       if (t_TracksLoose > 0) {
2028         ++looseHLT[0];
2029         for (int i = 0; i < 2; ++i)
2030           if (passt[i])
2031             ++looseHLT[i + 1];
2032       }
2033       if (t_TracksTight > 0) {
2034         ++tightHLT[0];
2035         for (int i = 0; i < 2; ++i)
2036           if (passt[i])
2037             ++tightHLT[i + 1];
2038       }
2039     }
2040     for (unsigned int k = 0; k < t_ietaAll->size(); ++k)
2041       h_eta[0]->Fill((*t_ietaAll)[k]);
2042     for (unsigned int k = 0; k < t_ietaGood->size(); ++k) {
2043       h_eta[1]->Fill((*t_ietaGood)[k]);
2044       if (t_trackType->size() > 0) {
2045         if ((*t_trackType)[k] > 0)
2046           h_eta[2]->Fill((*t_ietaGood)[k]);
2047         if ((*t_trackType)[k] > 1)
2048           h_eta[3]->Fill((*t_ietaGood)[k]);
2049       }
2050     }
2051   }
2052   double ymaxk(0);
2053   for (int i = 1; i <= h_eff[0]->GetNbinsX(); ++i) {
2054     double rat(0), drat(0);
2055     if (h_eta[0]->GetBinContent(i) > ymaxk)
2056       ymaxk = h_eta[0]->GetBinContent(i);
2057     if ((h_eta[1]->GetBinContent(i) > 0) && (h_eta[0]->GetBinContent(i) > 0)) {
2058       rat = h_eta[1]->GetBinContent(i) / h_eta[0]->GetBinContent(i);
2059       drat = rat * std::sqrt(pow((h_eta[1]->GetBinError(i) / h_eta[1]->GetBinContent(i)), 2) +
2060                              pow((h_eta[0]->GetBinError(i) / h_eta[0]->GetBinContent(i)), 2));
2061     }
2062     h_eff[0]->SetBinContent(i, rat);
2063     h_eff[0]->SetBinError(i, drat);
2064     if ((h_eta[1]->GetBinContent(i) > 0) && (h_eta[2]->GetBinContent(i) > 0)) {
2065       rat = h_eta[2]->GetBinContent(i) / h_eta[1]->GetBinContent(i);
2066       drat = rat * std::sqrt(pow((h_eta[2]->GetBinError(i) / h_eta[2]->GetBinContent(i)), 2) +
2067                              pow((h_eta[1]->GetBinError(i) / h_eta[1]->GetBinContent(i)), 2));
2068     } else {
2069       rat = drat = 0;
2070     }
2071     h_eff[1]->SetBinContent(i, rat);
2072     h_eff[1]->SetBinError(i, drat);
2073     if ((h_eta[1]->GetBinContent(i) > 0) && (h_eta[3]->GetBinContent(i) > 0)) {
2074       rat = h_eta[3]->GetBinContent(i) / h_eta[1]->GetBinContent(i);
2075       drat = rat * std::sqrt(pow((h_eta[3]->GetBinError(i) / h_eta[3]->GetBinContent(i)), 2) +
2076                              pow((h_eta[1]->GetBinError(i) / h_eta[1]->GetBinContent(i)), 2));
2077     } else {
2078       rat = drat = 0;
2079     }
2080     h_eff[1]->SetBinContent(i, rat);
2081     h_eff[1]->SetBinError(i, drat);
2082   }
2083   std::cout << "===== Remove " << duplicate << " events from " << nentries << "\n===== " << kount
2084             << " events passed trigger of which " << selected << " events get selected =====\n"
2085             << std::endl;
2086   std::cout << "===== " << l1 << " events passed L1 " << hlt << " events passed HLT and " << loose << ":" << tight
2087             << " events have at least 1 track candidate with loose:tight"
2088             << " isolation cut =====\n"
2089             << std::endl;
2090   for (int i = 0; i < 3; ++i) {
2091     char tbit[20];
2092     if (i == 0)
2093       sprintf(tbit, "Any");
2094     else
2095       sprintf(tbit, "%3d", bit_[i - 1]);
2096     std::cout << "Satisfying HLT trigger bit " << tbit << " Kount " << allHLT[i] << " Loose " << looseHLT[i]
2097               << " Tight " << tightHLT[i] << std::endl;
2098   }
2099
2100   gStyle->SetCanvasBorderMode(0);
2101   gStyle->SetCanvasColor(kWhite);
2102   gStyle->SetPadColor(kWhite);
2103   gStyle->SetFillColor(kWhite);
2104   gStyle->SetOptStat(1110);
2105   gStyle->SetOptTitle(0);
2106   int color[5] = {kBlack, kRed, kBlue, kMagenta, kCyan};
2107   int lines[5] = {1, 2, 3, 4, 5};
2108   TCanvas *pad1 = new TCanvas("c_track", "c_track", 500, 500);
2109   pad1->SetRightMargin(0.10);
2110   pad1->SetTopMargin(0.10);
2111   pad1->SetFillColor(kWhite);
2112   std::string titl1[3] = {"Reconstructed", "Propagated", "Saved"};
2113   TLegend *legend1 = new TLegend(0.11, 0.80, 0.50, 0.89);
2114   legend1->SetFillColor(kWhite);
2115   double ymax(0), xmax(0);
2116   for (int k = 0; k < 3; ++k) {
2117     int total(0), totaltk(0);
2118     for (int i = 1; i <= h_tk[k]->GetNbinsX(); ++i) {
2119       if (ymax < h_tk[k]->GetBinContent(i))
2120         ymax = h_tk[k]->GetBinContent(i);
2121       if (i > 1)
2122         total += (int)(h_tk[k]->GetBinContent(i));
2123       totaltk += (int)(h_tk[k]->GetBinContent(i)) * (i - 1);
2124       if (h_tk[k]->GetBinContent(i) > 0) {
2125         if (xmax < h_tk[k]->GetBinLowEdge(i) + h_tk[k]->GetBinWidth(i))
2126           xmax = h_tk[k]->GetBinLowEdge(i) + h_tk[k]->GetBinWidth(i);
2127       }
2128     }
2129     h_tk[k]->SetLineColor(color[k]);
2130     h_tk[k]->SetMarkerColor(color[k]);
2131     h_tk[k]->SetLineStyle(lines[k]);
2132     std::cout << h_tk[k]->GetTitle() << " Entries " << h_tk[k]->GetEntries() << " Events " << total << " Tracks "
2133               << totaltk << std::endl;
2134     legend1->AddEntry(h_tk[k], titl1[k].c_str(), "l");
2135   }
2136   int i1 = (int)(0.1 * xmax) + 1;
2137   xmax = 10.0 * i1;
2138   int i2 = (int)(0.01 * ymax) + 1;
2139
2140   ymax = 100.0 * i2;
2141   for (int k = 0; k < 3; ++k) {
2142     h_tk[k]->GetXaxis()->SetRangeUser(0, xmax);
2143     h_tk[k]->GetYaxis()->SetRangeUser(0.1, ymax);
2144     h_tk[k]->GetXaxis()->SetTitle("# Tracks");
2145     h_tk[k]->GetYaxis()->SetTitle("Events");
2146     h_tk[k]->GetYaxis()->SetLabelOffset(0.005);
2147     h_tk[k]->GetYaxis()->SetTitleOffset(1.20);
2148     if (k == 0)
2149       h_tk[k]->Draw("hist");
2150     else
2151       h_tk[k]->Draw("hist sames");
2152   }
2153   pad1->Update();
2154   pad1->SetLogy();
2155   ymax = 0.90;
2156   for (int k = 0; k < 3; ++k) {
2157     TPaveStats *st1 = (TPaveStats *)h_tk[k]->GetListOfFunctions()->FindObject("stats");
2158     if (st1 != NULL) {
2159       st1->SetLineColor(color[k]);
2160       st1->SetTextColor(color[k]);
2161       st1->SetY1NDC(ymax - 0.09);
2162       st1->SetY2NDC(ymax);
2163       st1->SetX1NDC(0.55);
2164       st1->SetX2NDC(0.90);
2165       ymax -= 0.09;
2166     }
2167   }
2168   pad1->Modified();
2169   pad1->Update();
2170   legend1->Draw("same");
2171   pad1->Update();
2172
2173   TCanvas *pad2 = new TCanvas("c_ieta", "c_ieta", 500, 500);
2174   pad2->SetRightMargin(0.10);
2175   pad2->SetTopMargin(0.10);
2176   pad2->SetFillColor(kWhite);
2177   pad2->SetLogy();
2178   std::string titl2[4] = {"All Tracks", "Selected Tracks", "Loose Isolation", "Tight Isolation"};
2179   TLegend *legend2 = new TLegend(0.11, 0.75, 0.50, 0.89);
2180   legend2->SetFillColor(kWhite);
2181   i2 = (int)(0.001 * ymaxk) + 1;
2182   ymax = 1000.0 * i2;
2183   for (int k = 0; k < 4; ++k) {
2184     h_eta[k]->GetYaxis()->SetRangeUser(1, ymax);
2185     h_eta[k]->SetLineColor(color[k]);
2186     h_eta[k]->SetMarkerColor(color[k]);
2187     h_eta[k]->SetLineStyle(lines[k]);
2188     h_eta[k]->GetXaxis()->SetTitle("i#eta");
2189     h_eta[k]->GetYaxis()->SetTitle("Tracks");
2190     h_eta[k]->GetYaxis()->SetLabelOffset(0.005);
2191     h_eta[k]->GetYaxis()->SetTitleOffset(1.20);
2192     legend2->AddEntry(h_eta[k], titl2[k].c_str(), "l");
2193     if (k == 0)
2194       h_eta[k]->Draw("hist");
2195     else
2196       h_eta[k]->Draw("hist sames");
2197   }
2198   pad2->Update();
2199   ymax = 0.90;
2200   //double ymin = 0.10;
2201   for (int k = 0; k < 4; ++k) {
2202     TPaveStats *st1 = (TPaveStats *)h_eta[k]->GetListOfFunctions()->FindObject("stats");
2203     if (st1 != NULL) {
2204       st1->SetLineColor(color[k]);
2205       st1->SetTextColor(color[k]);
2206       st1->SetY1NDC(ymax - 0.09);
2207       st1->SetY2NDC(ymax);
2208       st1->SetX1NDC(0.55);
2209       st1->SetX2NDC(0.90);
2210       ymax -= 0.09;
2211     }
2212   }
2213   pad2->Modified();
2214   pad2->Update();
2215   legend2->Draw("same");
2216   pad2->Update();
2217
2218   std::string titl3[3] = {"Selection", "Loose Isolation", "Tight Isolation"};
2219   TCanvas *pad3 = new TCanvas("c_effi", "c_effi", 500, 500);
2220   TLegend *legend3 = new TLegend(0.11, 0.785, 0.50, 0.89);
2221   pad3->SetRightMargin(0.10);
2222   pad3->SetTopMargin(0.10);
2223   pad3->SetFillColor(kWhite);
2224   pad3->SetLogy();
2225   for (int k = 0; k < 3; ++k) {
2226     h_eff[k]->SetStats(0);
2227     h_eff[k]->SetMarkerStyle(20);
2228     h_eff[k]->SetLineColor(color[k]);
2229     h_eff[k]->SetMarkerColor(color[k]);
2230     h_eff[k]->GetXaxis()->SetTitle("i#eta");
2231     h_eff[k]->GetYaxis()->SetTitle("Efficiency");
2232     h_eff[k]->GetYaxis()->SetLabelOffset(0.005);
2233     h_eff[k]->GetYaxis()->SetTitleOffset(1.20);
2234     if (k == 0)
2235       h_eff[k]->Draw("");
2236     else
2237       h_eff[k]->Draw("same");
2238     legend3->AddEntry(h_eff[k], titl3[k].c_str(), "l");
2239   }
2240   pad3->Modified();
2241   pad3->Update();
2242   legend3->Draw("same");
2243   pad3->Update();
2244
2245   std::string titl4[3] = {"All events", "Events with Loose Isolation", "Events with Tight Isolation"};
2246   TLegend *legend4 = new TLegend(0.11, 0.785, 0.50, 0.89);
2247   TCanvas *pad4 = new TCanvas("c_nvx", "c_nvx", 500, 500);
2248   pad4->SetRightMargin(0.10);
2249   pad4->SetTopMargin(0.10);
2250   pad4->SetFillColor(kWhite);
2251   pad4->SetLogy();
2252   for (int k = 0; k < 3; ++k) {
2253     h_pvx[k]->SetStats(1110);
2254     h_pvx[k]->SetMarkerStyle(20);
2255     h_pvx[k]->SetLineColor(color[k]);
2256     h_pvx[k]->SetMarkerColor(color[k]);
2257     h_pvx[k]->GetXaxis()->SetTitle("N_{PV}");
2258     h_pvx[k]->GetYaxis()->SetTitle("Events");
2259     h_pvx[k]->GetYaxis()->SetLabelOffset(0.005);
2260     h_pvx[k]->GetYaxis()->SetTitleOffset(1.20);
2261     if (k == 0)
2262       h_pvx[k]->Draw("");
2263     else
2264       h_pvx[k]->Draw("sames");
2265     legend4->AddEntry(h_pvx[k], titl4[k].c_str(), "l");
2266   }
2267   pad4->Update();
2268   ymax = 0.90;
2269   for (int k = 0; k < 3; ++k) {
2270     TPaveStats *st1 = (TPaveStats *)h_pvx[k]->GetListOfFunctions()->FindObject("stats");
2271     if (st1 != NULL) {
2272       st1->SetLineColor(color[k]);
2273       st1->SetTextColor(color[k]);
2274       st1->SetY1NDC(ymax - 0.09);
2275       st1->SetY2NDC(ymax);
2276       st1->SetX1NDC(0.55);
2277       st1->SetX2NDC(0.90);
2278       ymax -= 0.09;
2279     }
2280   }
2281   pad4->Modified();
2282   pad4->Update();
2283   legend4->Draw("same");
2284   pad4->Update();
2285 }
2286
2287 class CalibPlotProperties {
2288 public:
2289   TChain *fChain;  //!pointer to the analyzed TTree or TChain
2290   Int_t fCurrent;  //!current Tree number in a TChain
2291
2292   // Declaration of leaf types
2293   Int_t t_Run;
2294   Int_t t_Event;
2295   Int_t t_DataType;
2296   Int_t t_ieta;
2297   Int_t t_iphi;
2298   Double_t t_EventWeight;
2299   Int_t t_nVtx;
2300   Int_t t_nTrk;
2301   Int_t t_goodPV;
2302   Double_t t_l1pt;
2303   Double_t t_l1eta;
2304   Double_t t_l1phi;
2305   Double_t t_l3pt;
2306   Double_t t_l3eta;
2307   Double_t t_l3phi;
2308   Double_t t_p;
2309   Double_t t_pt;
2310   Double_t t_phi;
2311   Double_t t_mindR1;
2312   Double_t t_mindR2;
2313   Double_t t_eMipDR;
2314   Double_t t_eHcal;
2315   Double_t t_eHcal10;
2316   Double_t t_eHcal30;
2317   Double_t t_hmaxNearP;
2318   Double_t t_rhoh;
2319   Bool_t t_selectTk;
2320   Bool_t t_qltyFlag;
2321   Bool_t t_qltyMissFlag;
2322   Bool_t t_qltyPVFlag;
2323   Double_t t_gentrackP;
2324   std::vector<unsigned int> *t_DetIds;
2325   std::vector<double> *t_HitEnergies;
2326   std::vector<bool> *t_trgbits;
2327   std::vector<unsigned int> *t_DetIds1;
2328   std::vector<unsigned int> *t_DetIds3;
2329   std::vector<double> *t_HitEnergies1;
2330   std::vector<double> *t_HitEnergies3;
2331
2332   // List of branches
2333   TBranch *b_t_Run;           //!
2334   TBranch *b_t_Event;         //!
2335   TBranch *b_t_DataType;      //!
2336   TBranch *b_t_ieta;          //!
2337   TBranch *b_t_iphi;          //!
2338   TBranch *b_t_EventWeight;   //!
2339   TBranch *b_t_nVtx;          //!
2340   TBranch *b_t_nTrk;          //!
2341   TBranch *b_t_goodPV;        //!
2342   TBranch *b_t_l1pt;          //!
2343   TBranch *b_t_l1eta;         //!
2344   TBranch *b_t_l1phi;         //!
2345   TBranch *b_t_l3pt;          //!
2346   TBranch *b_t_l3eta;         //!
2347   TBranch *b_t_l3phi;         //!
2348   TBranch *b_t_p;             //!
2349   TBranch *b_t_pt;            //!
2350   TBranch *b_t_phi;           //!
2351   TBranch *b_t_mindR1;        //!
2352   TBranch *b_t_mindR2;        //!
2353   TBranch *b_t_eMipDR;        //!
2354   TBranch *b_t_eHcal;         //!
2355   TBranch *b_t_eHcal10;       //!
2356   TBranch *b_t_eHcal30;       //!
2357   TBranch *b_t_hmaxNearP;     //!
2358   TBranch *b_t_rhoh;          //!
2359   TBranch *b_t_selectTk;      //!
2360   TBranch *b_t_qltyFlag;      //!
2361   TBranch *b_t_qltyMissFlag;  //!
2362   TBranch *b_t_qltyPVFlag;    //!
2363   TBranch *b_t_gentrackP;     //!
2364   TBranch *b_t_DetIds;        //!
2365   TBranch *b_t_HitEnergies;   //!
2366   TBranch *b_t_trgbits;       //!
2367   TBranch *b_t_DetIds1;       //!
2368   TBranch *b_t_DetIds3;       //!
2369   TBranch *b_t_HitEnergies1;  //!
2370   TBranch *b_t_HitEnergies3;  //!
2371
2372   CalibPlotProperties(const char *fname,
2373                       const std::string &dirname,
2374                       const char *dupFileName,
2375                       const std::string &prefix = "",
2376                       const char *corrFileName = "",
2377                       const char *rcorFileName = "",
2378                       int puCorr = -2,
2379                       int flag = 1031,
2380                       bool dataMC = true,
2381                       int truncateFlag = 0,
2382                       bool useGen = false,
2383                       double scale = 1.0,
2384                       int useScale = 0,
2385                       int etalo = 0,
2386                       int etahi = 30,
2387                       int runlo = 0,
2388                       int runhi = 99999999,
2389                       int phimin = 1,
2390                       int phimax = 72,
2391                       int zside = 1,
2392                       int nvxlo = 0,
2393                       int nvxhi = 1000,
2394                       int rbx = 0,
2395                       bool exclude = false,
2396                       bool etamax = false);
2397   virtual ~CalibPlotProperties();
2398   virtual Int_t Cut(Long64_t entry);
2399   virtual Int_t GetEntry(Long64_t entry);
2400   virtual Long64_t LoadTree(Long64_t entry);
2401   virtual void Init(TChain *, const char *);
2402   virtual void Loop();
2403   virtual Bool_t Notify();
2404   virtual void Show(Long64_t entry = -1);
2405   bool goodTrack(double &eHcal, double &cut, const Long64_t &entry, bool debug);
2406   bool selectPhi(bool debug);
2407   void savePlot(const std::string &theName, bool append = true, bool all = false);
2408   void correctEnergy(double &ener, const Long64_t &entry);
2409
2410 private:
2411   static const unsigned int npbin = 6, kp50 = 3, ndepth = 7;
2412   CalibCorrFactor *corrFactor_;
2413   CalibCorr *cFactor_;
2414   CalibSelectRBX *cSelect_;
2415   const std::string fname_, dirnm_, prefix_, outFileName_;
2416   const int corrPU_, flag_;
2417   const bool dataMC_, useGen_;
2418   const int truncateFlag_;
2419   const int etalo_, etahi_;
2420   int runlo_, runhi_;
2421   const int phimin_, phimax_, zside_, nvxlo_, nvxhi_, rbx_;
2422   int ifDepth_;
2423   bool exclude_, corrE_, cutL1T_;
2424   bool includeRun_, getHist_;
2425   int flexibleSelect_;
2426   bool plotBasic_, plotEnergy_, plotHists_;
2427   double log2by18_;
2428   std::ofstream fileout_;
2429   std::vector<Long64_t> entries_;
2430   std::vector<std::pair<int, int> > events_;
2431   std::vector<double> etas_, ps_, dl1_;
2432   std::vector<int> nvx_, ietas_;
2433   TH1D *h_p[5], *h_eta[5], *h_nvtx;
2434   std::vector<TH1D *> h_eta0, h_eta1, h_eta2, h_eta3, h_eta4;
2435   std::vector<TH1D *> h_dL1, h_vtx;
2436   std::vector<TH1D *> h_etaEH[npbin], h_etaEp[npbin], h_etaEE[npbin];
2437   std::vector<TH1D *> h_mom, h_eEcal, h_eHcal;
2438   std::vector<TH1F *> h_bvlist, h_bvlist2, h_evlist, h_evlist2;
2439   std::vector<TH1F *> h_bvlist3, h_evlist3;
2440   TH2F *h_etaE;
2441 };
2442
2443 CalibPlotProperties::CalibPlotProperties(const char *fname,
2444                                          const std::string &dirnm,
2445                                          const char *dupFileName,
2446                                          const std::string &prefix,
2447                                          const char *corrFileName,
2448                                          const char *rcorFileName,
2449                                          int puCorr,
2450                                          int flag,
2451                                          bool dataMC,
2452                                          int truncate,
2453                                          bool useGen,
2454                                          double scl,
2455                                          int useScale,
2456                                          int etalo,
2457                                          int etahi,
2458                                          int runlo,
2459                                          int runhi,
2460                                          int phimin,
2461                                          int phimax,
2462                                          int zside,
2463                                          int nvxlo,
2464                                          int nvxhi,
2465                                          int rbx,
2466                                          bool exc,
2467                                          bool etam)
2468     : corrFactor_(nullptr),
2469       cFactor_(nullptr),
2470       cSelect_(nullptr),
2471       fname_(fname),
2472       dirnm_(dirnm),
2473       prefix_(prefix),
2474       corrPU_(puCorr),
2475       flag_(flag),
2476       dataMC_(dataMC),
2477       useGen_(useGen),
2478       truncateFlag_(truncate),
2479       etalo_(etalo),
2480       etahi_(etahi),
2481       runlo_(runlo),
2482       runhi_(runhi),
2483       phimin_(phimin),
2484       phimax_(phimax),
2485       zside_(zside),
2486       nvxlo_(nvxlo),
2487       nvxhi_(nvxhi),
2488       rbx_(rbx),
2489       exclude_(exc),
2490       includeRun_(true) {
2491   // if parameter tree is not specified (or zero), connect the file
2492   // used to generate this class and read the Tree
2493
2494   flexibleSelect_ = ((flag_ / 1) % 10);
2495   plotBasic_ = (((flag_ / 10) % 10) > 0);
2496   plotEnergy_ = (((flag_ / 10) % 10) > 0);
2497   int oneplace = ((flag_ / 1000) % 10);
2498   cutL1T_ = (oneplace % 2);
2499   bool marina = ((oneplace / 2) % 2);
2500   ifDepth_ = ((flag_ / 10000) % 10);
2501   plotHists_ = (((flag_ / 100000) % 10) > 0);
2502   log2by18_ = std::log(2.5) / 18.0;
2503   if (runlo_ < 0 || runhi_ < 0) {
2504     runlo_ = std::abs(runlo_);
2505     runhi_ = std::abs(runhi_);
2506     includeRun_ = false;
2507   }
2508   char treeName[400];
2509   sprintf(treeName, "%s/CalibTree", dirnm.c_str());
2510   TChain *chain = new TChain(treeName);
2511   std::cout << "Create a chain for " << treeName << " from " << fname << " flags " << flexibleSelect_ << "|"
2512             << plotBasic_ << "|"
2513             << "|" << plotEnergy_ << "|" << plotHists_ << "|" << corrPU_ << "\n cons " << log2by18_ << " eta range "
2514             << etalo_ << ":" << etahi_ << " run range " << runlo_ << ":" << runhi_ << " (inclusion flag " << includeRun_
2515             << ")\n Selection of RBX" << rbx << " Vertex Range " << nvxlo_ << ":" << nvxhi_
2516             << "\n corrFileName: " << corrFileName << " useScale " << useScale << ":" << scl << ":" << etam
2517             << "\n rcorFileName: " << rcorFileName << " flag " << ifDepth_ << std::endl;
2518   if (!fillChain(chain, fname)) {
2519     std::cout << "*****No valid tree chain can be obtained*****" << std::endl;
2520   } else {
2521     std::cout << "Proceed with a tree chain with " << chain->GetEntries() << " entries" << std::endl;
2522     Init(chain, dupFileName);
2523     corrFactor_ = new CalibCorrFactor(corrFileName, useScale, scl, etam, marina, false);
2524     if (std::string(rcorFileName) != "") {
2525       cFactor_ = new CalibCorr(rcorFileName, ifDepth_, false);
2526     } else {
2527       ifDepth_ = 0;
2528     }
2529     if (rbx != 0)
2530       cSelect_ = new CalibSelectRBX(rbx, false);
2531   }
2532 }
2533
2534 CalibPlotProperties::~CalibPlotProperties() {
2535   delete corrFactor_;
2536   delete cFactor_;
2537   delete cSelect_;
2538   if (!fChain)
2539     return;
2540   delete fChain->GetCurrentFile();
2541 }
2542
2543 Int_t CalibPlotProperties::GetEntry(Long64_t entry) {
2544   // Read contents of entry.
2545   if (!fChain)
2546     return 0;
2547   return fChain->GetEntry(entry);
2548 }
2549
2550 Long64_t CalibPlotProperties::LoadTree(Long64_t entry) {
2551   // Set the environment to read one entry
2552   if (!fChain)
2553     return -5;
2554   Long64_t centry = fChain->LoadTree(entry);
2555   if (centry < 0)
2556     return centry;
2557   if (!fChain->InheritsFrom(TChain::Class()))
2558     return centry;
2559   TChain *chain = (TChain *)fChain;
2560   if (chain->GetTreeNumber() != fCurrent) {
2561     fCurrent = chain->GetTreeNumber();
2562     Notify();
2563   }
2564   return centry;
2565 }
2566
2567 void CalibPlotProperties::Init(TChain *tree, const char *dupFileName) {
2568   // The Init() function is called when the selector needs to initialize
2569   // a new tree or chain. Typically here the branch addresses and branch
2570   // pointers of the tree will be set.
2571   // It is normally not necessary to make changes to the generated
2572   // code, but the routine can be extended by the user if needed.
2573   // Init() will be called many times when running on PROOF
2574   // (once per file to be processed).
2575
2576   // Set object pointer
2577   t_DetIds = 0;
2578   t_DetIds1 = 0;
2579   t_DetIds3 = 0;
2580   t_HitEnergies = 0;
2581   t_HitEnergies1 = 0;
2582   t_HitEnergies3 = 0;
2583   t_trgbits = 0;
2584   // Set branch addresses and branch pointers
2585   fChain = tree;
2586   fCurrent = -1;
2587   if (!tree)
2588     return;
2589   fChain->SetMakeClass(1);
2590
2591   fChain->SetBranchAddress("t_Run", &t_Run, &b_t_Run);
2592   fChain->SetBranchAddress("t_Event", &t_Event, &b_t_Event);
2593   fChain->SetBranchAddress("t_DataType", &t_DataType, &b_t_DataType);
2594   fChain->SetBranchAddress("t_ieta", &t_ieta, &b_t_ieta);
2595   fChain->SetBranchAddress("t_iphi", &t_iphi, &b_t_iphi);
2596   fChain->SetBranchAddress("t_EventWeight", &t_EventWeight, &b_t_EventWeight);
2597   fChain->SetBranchAddress("t_nVtx", &t_nVtx, &b_t_nVtx);
2598   fChain->SetBranchAddress("t_nTrk", &t_nTrk, &b_t_nTrk);
2599   fChain->SetBranchAddress("t_goodPV", &t_goodPV, &b_t_goodPV);
2600   fChain->SetBranchAddress("t_l1pt", &t_l1pt, &b_t_l1pt);
2601   fChain->SetBranchAddress("t_l1eta", &t_l1eta, &b_t_l1eta);
2602   fChain->SetBranchAddress("t_l1phi", &t_l1phi, &b_t_l1phi);
2603   fChain->SetBranchAddress("t_l3pt", &t_l3pt, &b_t_l3pt);
2604   fChain->SetBranchAddress("t_l3eta", &t_l3eta, &b_t_l3eta);
2605   fChain->SetBranchAddress("t_l3phi", &t_l3phi, &b_t_l3phi);
2606   fChain->SetBranchAddress("t_p", &t_p, &b_t_p);
2607   fChain->SetBranchAddress("t_pt", &t_pt, &b_t_pt);
2608   fChain->SetBranchAddress("t_phi", &t_phi, &b_t_phi);
2609   fChain->SetBranchAddress("t_mindR1", &t_mindR1, &b_t_mindR1);
2610   fChain->SetBranchAddress("t_mindR2", &t_mindR2, &b_t_mindR2);
2611   fChain->SetBranchAddress("t_eMipDR", &t_eMipDR, &b_t_eMipDR);
2612   fChain->SetBranchAddress("t_eHcal", &t_eHcal, &b_t_eHcal);
2613   fChain->SetBranchAddress("t_eHcal10", &t_eHcal10, &b_t_eHcal10);
2614   fChain->SetBranchAddress("t_eHcal30", &t_eHcal30, &b_t_eHcal30);
2615   fChain->SetBranchAddress("t_hmaxNearP", &t_hmaxNearP, &b_t_hmaxNearP);
2616   fChain->SetBranchAddress("t_rhoh", &t_rhoh, &b_t_rhoh);
2617   fChain->SetBranchAddress("t_selectTk", &t_selectTk, &b_t_selectTk);
2618   fChain->SetBranchAddress("t_qltyFlag", &t_qltyFlag, &b_t_qltyFlag);
2619   fChain->SetBranchAddress("t_qltyMissFlag", &t_qltyMissFlag, &b_t_qltyMissFlag);
2620   fChain->SetBranchAddress("t_qltyPVFlag", &t_qltyPVFlag, &b_t_qltyPVFlag);
2621   fChain->SetBranchAddress("t_gentrackP", &t_gentrackP, &b_t_gentrackP);
2622   fChain->SetBranchAddress("t_DetIds", &t_DetIds, &b_t_DetIds);
2623   fChain->SetBranchAddress("t_HitEnergies", &t_HitEnergies, &b_t_HitEnergies);
2624   fChain->SetBranchAddress("t_trgbits", &t_trgbits, &b_t_trgbits);
2625   fChain->SetBranchAddress("t_DetIds1", &t_DetIds1, &b_t_DetIds1);
2626   fChain->SetBranchAddress("t_DetIds3", &t_DetIds3, &b_t_DetIds3);
2627   fChain->SetBranchAddress("t_HitEnergies1", &t_HitEnergies1, &b_t_HitEnergies1);
2628   fChain->SetBranchAddress("t_HitEnergies3", &t_HitEnergies3, &b_t_HitEnergies3);
2629   Notify();
2630
2631   if (std::string(dupFileName) != "") {
2632     ifstream infil1(dupFileName);
2633     if (!infil1.is_open()) {
2634       std::cout << "Cannot open duplicate file " << dupFileName << std::endl;
2635     } else {
2636       while (1) {
2637         Long64_t jentry;
2638         infil1 >> jentry;
2639         if (!infil1.good())
2640           break;
2641         entries_.push_back(jentry);
2642       }
2643       infil1.close();
2644       std::cout << "Reads a list of " << entries_.size() << " events from " << dupFileName << std::endl;
2645     }
2646   } else {
2647     std::cout << "No duplicate events in the input file" << std::endl;
2648   }
2649
2650   int ipbin[npbin] = {10, 20, 30, 40, 60, 100};
2651   for (unsigned int i = 0; i < npbin; ++i)
2652     ps_.push_back((double)(ipbin[i]));
2653   int ietas[4] = {0, 13, 18, 23};
2654   for (int i = 0; i < 4; ++i)
2655     ietas_.push_back(ietas[i]);
2656
2657   char name[20], title[200];
2658   unsigned int kp = ps_.size() - 1;
2659   unsigned int kk(0);
2660   std::string titl[5] = {
2661       "all tracks", "good quality tracks", "selected tracks", "isolated good tracks", "tracks having MIP in ECAL"};
2662
2663   if (plotBasic_) {
2664     std::cout << "Book Basic Histos" << std::endl;
2665     for (int k = 0; k < 5; ++k) {
2666       sprintf(name, "%sp%d", prefix_.c_str(), k);
2667       sprintf(title, "Momentum for %s", titl[k].c_str());
2668       h_p[k] = new TH1D(name, title, 100, 10.0, 110.0);
2669       sprintf(name, "%seta%d", prefix_.c_str(), k);
2670       sprintf(title, "#eta for %s", titl[k].c_str());
2671       h_eta[k] = new TH1D(name, title, 60, -30.0, 30.0);
2672     }
2673     for (unsigned int k = 0; k < kp; ++k) {
2674       sprintf(name, "%seta0%d", prefix_.c_str(), k);
2675       sprintf(title, "#eta for %s (p = %d:%d GeV)", titl[0].c_str(), ipbin[k], ipbin[k + 1]);
2676       h_eta0.push_back(new TH1D(name, title, 60, -30.0, 30.0));
2677       kk = h_eta0.size() - 1;
2678       h_eta0[kk]->Sumw2();
2679       sprintf(name, "%seta1%d", prefix_.c_str(), k);
2680       sprintf(title, "#eta for %s (p = %d:%d GeV)", titl[1].c_str(), ipbin[k], ipbin[k + 1]);
2681       h_eta1.push_back(new TH1D(name, title, 60, -30.0, 30.0));
2682       kk = h_eta1.size() - 1;
2683       h_eta1[kk]->Sumw2();
2684       sprintf(name, "%seta2%d", prefix_.c_str(), k);
2685       sprintf(title, "#eta for %s (p = %d:%d GeV)", titl[2].c_str(), ipbin[k], ipbin[k + 1]);
2686       h_eta2.push_back(new TH1D(name, title, 60, -30.0, 30.0));
2687       kk = h_eta2.size() - 1;
2688       h_eta2[kk]->Sumw2();
2689       sprintf(name, "%seta3%d", prefix_.c_str(), k);
2690       sprintf(title, "#eta for %s (p = %d:%d GeV)", titl[3].c_str(), ipbin[k], ipbin[k + 1]);
2691       h_eta3.push_back(new TH1D(name, title, 60, -30.0, 30.0));
2692       kk = h_eta3.size() - 1;
2693       h_eta3[kk]->Sumw2();
2694       sprintf(name, "%seta4%d", prefix_.c_str(), k);
2695       sprintf(title, "#eta for %s (p = %d:%d GeV)", titl[4].c_str(), ipbin[k], ipbin[k + 1]);
2696       h_eta4.push_back(new TH1D(name, title, 60, -30.0, 30.0));
2697       kk = h_eta4.size() - 1;
2698       h_eta4[kk]->Sumw2();
2699       sprintf(name, "%sdl1%d", prefix_.c_str(), k);
2700       sprintf(title, "Distance from L1 (p = %d:%d GeV)", ipbin[k], ipbin[k + 1]);
2701       h_dL1.push_back(new TH1D(name, title, 160, 0.0, 8.0));
2702       kk = h_dL1.size() - 1;
2703       h_dL1[kk]->Sumw2();
2704       sprintf(name, "%svtx%d", prefix_.c_str(), k);
2705       sprintf(title, "N_{Vertex} (p = %d:%d GeV)", ipbin[k], ipbin[k + 1]);
2706       h_vtx.push_back(new TH1D(name, title, 100, 0.0, 100.0));
2707       kk = h_vtx.size() - 1;
2708       h_vtx[kk]->Sumw2();
2709     }
2710   }
2711
2712   if (plotEnergy_) {
2713     std::cout << "Make plots for good tracks" << std::endl;
2714     for (unsigned int k = 0; k < kp; ++k) {
2715       for (int j = etalo_; j <= etahi_ + 1; ++j) {
2716         sprintf(name, "%senergyH%d%d", prefix_.c_str(), k, j);
2717         if (j > etahi_)
2718           sprintf(title,
2719                   "HCAL energy for %s (p = %d:%d GeV |#eta| = %d:%d)",
2720                   titl[3].c_str(),
2721                   ipbin[k],
2722                   ipbin[k + 1],
2723                   etalo_,
2724                   etahi_);
2725         else
2726           sprintf(title, "HCAL energy for %s (p = %d:%d GeV |#eta| = %d)", titl[3].c_str(), ipbin[k], ipbin[k + 1], j);
2727         h_etaEH[k].push_back(new TH1D(name, title, 200, 0, 200));
2728         kk = h_etaEH[k].size() - 1;
2729         h_etaEH[k][kk]->Sumw2();
2730         sprintf(name, "%senergyP%d%d", prefix_.c_str(), k, j);
2731         if (j > etahi_)
2732           sprintf(title,
2733                   "momentum for %s (p = %d:%d GeV |#eta| = %d:%d)",
2734                   titl[3].c_str(),
2735                   ipbin[k],
2736                   ipbin[k + 1],
2737                   etalo_,
2738                   etahi_);
2739         else
2740           sprintf(title, "momentum for %s (p = %d:%d GeV |#eta| = %d)", titl[3].c_str(), ipbin[k], ipbin[k + 1], j);
2741         h_etaEp[k].push_back(new TH1D(name, title, 100, 0, 100));
2742         kk = h_etaEp[k].size() - 1;
2743         h_etaEp[k][kk]->Sumw2();
2744         sprintf(name, "%senergyE%d%d", prefix_.c_str(), k, j);
2745         if (j > etahi_)
2746           sprintf(title,
2747                   "ECAL energy for %s (p = %d:%d GeV |#eta| = %d:%d)",
2748                   titl[3].c_str(),
2749                   ipbin[k],
2750                   ipbin[k + 1],
2751                   etalo_,
2752                   etahi_);
2753         else
2754           sprintf(title, "ECAL energy for %s (p = %d:%d GeV |#eta| = %d)", titl[3].c_str(), ipbin[k], ipbin[k + 1], j);
2755         h_etaEE[k].push_back(new TH1D(name, title, 100, 0, 10));
2756         kk = h_etaEE[k].size() - 1;
2757         h_etaEE[k][kk]->Sumw2();
2758       }
2759     }
2760
2761     for (unsigned int j = 0; j < ietas_.size(); ++j) {
2762       sprintf(name, "%senergyH%d", prefix_.c_str(), j);
2763       if (j == 0)
2764         sprintf(title, "HCAL energy for %s (All)", titl[3].c_str());
2765       else
2766         sprintf(title, "HCAL energy for %s (|#eta| = %d:%d)", titl[3].c_str(), ietas_[j - 1], ietas_[j]);
2767       h_eHcal.push_back(new TH1D(name, title, 200, 0, 200));
2768       kk = h_eHcal.size() - 1;
2769       h_eHcal[kk]->Sumw2();
2770       sprintf(name, "%senergyP%d", prefix_.c_str(), j);
2771       if (j == 0)
2772         sprintf(title, "Track momentum for %s (All)", titl[3].c_str());
2773       else
2774         sprintf(title, "Track momentum for %s (|#eta| = %d:%d)", titl[3].c_str(), ietas_[j - 1], ietas_[j]);
2775       h_mom.push_back(new TH1D(name, title, 100, 0, 100));
2776       kk = h_mom.size() - 1;
2777       h_mom[kk]->Sumw2();
2778       sprintf(name, "%senergyE%d", prefix_.c_str(), j);
2779       if (j == 0)
2780         sprintf(title, "ECAL energy for %s (All)", titl[3].c_str());
2781       else
2782         sprintf(title, "ECAL energy for %s (|#eta| = %d:%d)", titl[3].c_str(), ietas_[j - 1], ietas_[j]);
2783       h_eEcal.push_back(new TH1D(name, title, 100, 0, 10));
2784       kk = h_eEcal.size() - 1;
2785       h_eEcal[kk]->Sumw2();
2786     }
2787   }
2788
2789   if (plotHists_) {
2790     h_nvtx = new TH1D("hnvtx", "Number of vertices", 10, 0, 100);
2791     h_nvtx->Sumw2();
2792     for (unsigned int i = 0; i < ndepth; i++) {
2793       sprintf(name, "b_edepth%d", i);
2794       sprintf(title, "Total RecHit energy in depth %d (Barrel)", i + 1);
2795       h_bvlist.push_back(new TH1F(name, title, 1000, 0, 100));
2796       h_bvlist[i]->Sumw2();
2797       sprintf(name, "b_recedepth%d", i);
2798       sprintf(title, "RecHit energy in depth %d (Barrel)", i + 1);
2799       h_bvlist2.push_back(new TH1F(name, title, 1000, 0, 100));
2800       h_bvlist2[i]->Sumw2();
2801       sprintf(name, "b_nrecdepth%d", i);
2802       sprintf(title, "#RecHits in depth %d (Barrel)", i + 1);
2803       h_bvlist3.push_back(new TH1F(name, title, 1000, 0, 100));
2804       h_bvlist3[i]->Sumw2();
2805       sprintf(name, "e_edepth%d", i);
2806       sprintf(title, "Total RecHit energy in depth %d (Endcap)", i + 1);
2807       h_evlist.push_back(new TH1F(name, title, 1000, 0, 100));
2808       h_evlist[i]->Sumw2();
2809       sprintf(name, "e_recedepth%d", i);
2810       sprintf(title, "RecHit energy in depth %d (Endcap)", i + 1);
2811       h_evlist2.push_back(new TH1F(name, title, 1000, 0, 100));
2812       h_evlist2[i]->Sumw2();
2813       sprintf(name, "e_nrecdepth%d", i);
2814       sprintf(title, "#RecHits in depth %d (Endcap)", i + 1);
2815       h_evlist3.push_back(new TH1F(name, title, 1000, 0, 100));
2816       h_evlist3[i]->Sumw2();
2817     }
2818     h_etaE = new TH2F("heta", "", 50, -25, 25, 100, 0, 100);
2819   }
2820 }
2821
2822 Bool_t CalibPlotProperties::Notify() {
2823   // The Notify() function is called when a new file is opened. This
2824   // can be either for a new TTree in a TChain or when when a new TTree
2825   // is started when using PROOF. It is normally not necessary to make changes
2826   // to the generated code, but the routine can be extended by the
2827   // user if needed. The return value is currently not used.
2828
2829   return kTRUE;
2830 }
2831
2832 void CalibPlotProperties::Show(Long64_t entry) {
2833   // Print contents of entry.
2834   // If entry is not specified, print current entry
2835   if (!fChain)
2836     return;
2837   fChain->Show(entry);
2838 }
2839
2840 Int_t CalibPlotProperties::Cut(Long64_t) {
2841   // This function may be called from Loop.
2842   // returns  1 if entry is accepted.
2843   // returns -1 otherwise.
2844   return 1;
2845 }
2846
2847 void CalibPlotProperties::Loop() {
2848   //   In a ROOT session, you can do:
2849   //      Root > .L CalibMonitor.C
2850   //      Root > CalibMonitor t
2851   //      Root > t.GetEntry(12); // Fill t data members with entry number 12
2852   //      Root > t.Show();       // Show values of entry 12
2853   //      Root > t.Show(16);     // Read and show values of entry 16
2854   //      Root > t.Loop();       // Loop on all entries
2855   //
2856
2857   //   This is the loop skeleton where:
2858   //      jentry is the global entry number in the chain
2859   //      ientry is the entry number in the current Tree
2860   //   Note that the argument to GetEntry must be:
2861   //      jentry for TChain::GetEntry
2862   //      ientry for TTree::GetEntry and TBranch::GetEntry
2863   //
2864   //       To read only selected branches, Insert statements like:
2865   // METHOD1:
2866   //    fChain->SetBranchStatus("*",0);  // disable all branches
2867   //    fChain->SetBranchStatus("branchname",1);  // activate branchname
2868   // METHOD2: replace line
2869   //    fChain->GetEntry(jentry);       //read all branches
2870   //by  b_branchname->GetEntry(ientry); //read only this branch
2871   const bool debug(false);
2872
2873   if (fChain == 0)
2874     return;
2875
2876   // Find list of duplicate events
2877   Long64_t nentries = fChain->GetEntriesFast();
2878   std::cout << "Total entries " << nentries << std::endl;
2879   Long64_t nbytes(0), nb(0);
2880   unsigned int duplicate(0), good(0), kount(0);
2881   double sel(0), selHB(0), selHE(0);
2882   for (Long64_t jentry = 0; jentry < nentries; jentry++) {
2883     Long64_t ientry = LoadTree(jentry);
2884     if (ientry < 0)
2885       break;
2886     nb = fChain->GetEntry(jentry);
2887     nbytes += nb;
2888     if (jentry % 1000000 == 0)
2889       std::cout << "Entry " << jentry << " Run " << t_Run << " Event " << t_Event << std::endl;
2890     bool select = (std::find(entries_.begin(), entries_.end(), jentry) == entries_.end());
2891     if (!select) {
2892       ++duplicate;
2893       if (debug)
2894         std::cout << "Duplicate event " << t_Run << " " << t_Event << " " << t_p << std::endl;
2895       continue;
2896     }
2897     bool selRun = (includeRun_ ? ((t_Run >= runlo_) && (t_Run <= runhi_)) : ((t_Run < runlo_) || (t_Run > runhi_)));
2898     select =
2899         (selRun && (fabs(t_ieta) >= etalo_) && (fabs(t_ieta) <= etahi_) && (t_nVtx >= nvxlo_) && (t_nVtx <= nvxhi_));
2900     if (!select) {
2901       if (debug)
2902         std::cout << "Run # " << t_Run << " out of range of " << runlo_ << ":" << runhi_ << " or ieta " << t_ieta
2903                   << " (" << etalo_ << ":" << etahi_ << ") or nvtx " << t_nVtx << " (" << nvxlo_ << ":" << nvxhi_
2904                   << ") out of range" << std::endl;
2905       continue;
2906     }
2907     if (cSelect_ != nullptr) {
2908       if (exclude_) {
2909         if (cSelect_->isItRBX(t_DetIds))
2910           continue;
2911       } else {
2912         if (!(cSelect_->isItRBX(t_ieta, t_iphi)))
2913           continue;
2914       }
2915     }
2916     select = (!cutL1T_ || (t_mindR1 >= 0.5));
2917     if (!select) {
2918       if (debug)
2919         std::cout << "Reject Run # " << t_Run << " Event # " << t_Event << " too close to L1 trigger " << t_mindR1
2920                   << std::endl;
2921       continue;
2922     }
2923     select = ((events_.size() == 0) ||
2924               (std::find(events_.begin(), events_.end(), std::pair<int, int>(t_Run, t_Event)) != events_.end()));
2925     if (!select) {
2926       if (debug)
2927         std::cout << "Reject Run # " << t_Run << " Event # " << t_Event << " not in the selection list" << std::endl;
2928       continue;
2929     }
2930
2931     if ((ifDepth_ == 3) && (cFactor_ != nullptr) && (cFactor_->absent(ientry)))
2932       continue;
2933     // if (Cut(ientry) < 0) continue;
2934     unsigned int kp = ps_.size();
2935     double pmom = (useGen_ && (t_gentrackP > 0)) ? t_gentrackP : t_p;
2936     for (unsigned int k = 1; k < ps_.size(); ++k) {
2937       if (pmom >= ps_[k - 1] && pmom < ps_[k]) {
2938         kp = k - 1;
2939         break;
2940       }
2941     }
2942     int jp1 = (((std::abs(t_ieta) >= etalo_) && (std::abs(t_ieta) <= etahi_)) ? (std::abs(t_ieta) - etalo_) : -1);
2943     int jp2 = (etahi_ - etalo_ + 1);
2944     unsigned int je1 = ietas_.size();
2945     for (unsigned int j = 1; j < ietas_.size(); ++j) {
2946       if (std::abs(t_ieta) > ietas_[j - 1] && std::abs(t_ieta) <= ietas_[j]) {
2947         je1 = j;
2948         break;
2949       }
2950     }
2951     int je2 = (je1 != ietas_.size()) ? 0 : -1;
2952     if (debug)
2953       std::cout << "Bin " << kp << ":" << je1 << ":" << je2 << ":" << jp1 << ":" << jp2 << std::endl;
2954     double cut = (pmom > 20) ? ((flexibleSelect_ == 0) ? 2.0 : 10.0) : 0.0;
2955     double rcut(-1000.0);
2956
2957     // Selection of good track and energy measured in Hcal
2958     double eHcal(t_eHcal);
2959     if (corrFactor_->doCorr() || (cFactor_ != 0)) {
2960       eHcal = 0;
2961       for (unsigned int k = 0; k < t_HitEnergies->size(); ++k) {
2962         // The masks are defined in DataFormats/HcalDetId/interface/HcalDetId.h
2963         double cfac(1.0);
2964         if (corrFactor_->doCorr()) {
2965           unsigned int id = truncateId((*t_DetIds)[k], truncateFlag_, false);
2966           cfac = corrFactor_->getCorr(id);
2967         }
2968         if (cFactor_ != 0)
2969           cfac *= cFactor_->getCorr(t_Run, (*t_DetIds)[k]);
2970         eHcal += (cfac * ((*t_HitEnergies)[k]));
2971         if (debug) {
2972           int subdet, zside, ieta, iphi, depth;
2973           unpackDetId((*t_DetIds)[k], subdet, zside, ieta, iphi, depth);
2974           std::cout << zside << ":" << ieta << ":" << depth << " Corr " << cfac << " " << (*t_HitEnergies)[k] << " Out "
2975                     << eHcal << std::endl;
2976         }
2977       }
2978     }
2979     bool goodTk = goodTrack(eHcal, cut, jentry, debug);
2980     bool selPhi = selectPhi(debug);
2981     double rat = (pmom > 0) ? (eHcal / (pmom - t_eMipDR)) : 1.0;
2982     if (debug)
2983       std::cout << "Entry " << jentry << " p|eHcal|ratio " << pmom << "|" << t_eHcal << "|" << eHcal << "|" << rat
2984                 << "|" << kp << " Cuts " << t_qltyFlag << "|" << t_selectTk << "|" << (t_hmaxNearP < cut) << "|"
2985                 << (t_eMipDR < 1.0) << "|" << goodTk << "|" << (rat > rcut) << " Select Phi " << selPhi << std::endl;
2986     if (plotBasic_) {
2987       h_p[0]->Fill(pmom, t_EventWeight);
2988       h_eta[0]->Fill(t_ieta, t_EventWeight);
2989       if (kp < ps_.size())
2990         h_eta0[kp]->Fill(t_ieta, t_EventWeight);
2991       if (t_qltyFlag) {
2992         h_p[1]->Fill(pmom, t_EventWeight);
2993         h_eta[1]->Fill(t_ieta, t_EventWeight);
2994         if (kp < ps_.size())
2995           h_eta1[kp]->Fill(t_ieta, t_EventWeight);
2996         if (t_selectTk) {
2997           h_p[2]->Fill(pmom, t_EventWeight);
2998           h_eta[2]->Fill(t_ieta, t_EventWeight);
2999           if (kp < ps_.size())
3000             h_eta2[kp]->Fill(t_ieta, t_EventWeight);
3001           if (t_hmaxNearP < cut) {
3002             h_p[3]->Fill(pmom, t_EventWeight);
3003             h_eta[3]->Fill(t_ieta, t_EventWeight);
3004             if (kp < ps_.size())
3005               h_eta3[kp]->Fill(t_ieta, t_EventWeight);
3006             if (t_eMipDR < 1.0) {
3007               h_p[4]->Fill(pmom, t_EventWeight);
3008               h_eta[4]->Fill(t_ieta, t_EventWeight);
3009               if (kp < ps_.size()) {
3010                 h_eta4[kp]->Fill(t_ieta, t_EventWeight);
3011                 h_dL1[kp]->Fill(t_mindR1, t_EventWeight);
3012                 h_vtx[kp]->Fill(t_goodPV, t_EventWeight);
3013               }
3014             }
3015           }
3016         }
3017       }
3018     }
3019
3020     if (goodTk && kp < ps_.size() && selPhi) {
3021       if (rat > rcut) {
3022         if (plotEnergy_) {
3023           if (jp1 >= 0) {
3024             h_etaEH[kp][jp1]->Fill(eHcal, t_EventWeight);
3025             h_etaEH[kp][jp2]->Fill(eHcal, t_EventWeight);
3026             h_etaEp[kp][jp1]->Fill(pmom, t_EventWeight);
3027             h_etaEp[kp][jp2]->Fill(pmom, t_EventWeight);
3028             h_etaEE[kp][jp1]->Fill(t_eMipDR, t_EventWeight);
3029             h_etaEE[kp][jp2]->Fill(t_eMipDR, t_EventWeight);
3030           }
3031           if (kp == kp50) {
3032             if (je1 != ietas_.size()) {
3033               h_eHcal[je1]->Fill(eHcal, t_EventWeight);
3034               h_eHcal[je2]->Fill(eHcal, t_EventWeight);
3035               h_mom[je1]->Fill(pmom, t_EventWeight);
3036               h_mom[je2]->Fill(pmom, t_EventWeight);
3037               h_eEcal[je1]->Fill(t_eMipDR, t_EventWeight);
3038               h_eEcal[je2]->Fill(t_eMipDR, t_EventWeight);
3039             }
3040           }
3041         }
3042
3043         if (plotHists_) {
3044           h_nvtx->Fill(t_nVtx);
3045           bool bad(false);
3046           for (unsigned int k = 0; k < t_HitEnergies->size(); ++k) {
3047             unsigned int id = truncateId((*t_DetIds)[k], truncateFlag_, false);
3048             double cfac = corrFactor_->getCorr(id);
3049             if (cFactor_ != 0)
3050               cfac *= cFactor_->getCorr(t_Run, (*t_DetIds)[k]);
3051             double ener = cfac * (*t_HitEnergies)[k];
3052             if (corrPU_)
3053               correctEnergy(ener, jentry);
3054             if (ener < 0.001)
3055               bad = true;
3056           }
3057           if ((!bad) && (std::fabs(rat - 1) < 0.15) && (kp == kp50) &&
3058               ((std::abs(t_ieta) < 15) || (std::abs(t_ieta) > 17))) {
3059             float weight = (dataMC_ ? t_EventWeight : t_EventWeight * puweight(t_nVtx));
3060             h_etaE->Fill(t_ieta, eHcal, weight);
3061             sel += weight;
3062             std::vector<float> bv(7, 0.0f), ev(7, 0.0f);
3063             std::vector<int> bnrec(7, 0), enrec(7, 0);
3064             double eb(0), ee(0);
3065             for (unsigned int k = 0; k < t_HitEnergies->size(); ++k) {
3066               unsigned int id = truncateId((*t_DetIds)[k], truncateFlag_, false);
3067               double cfac = corrFactor_->getCorr(id);
3068               if (cFactor_ != 0)
3069                 cfac *= cFactor_->getCorr(t_Run, (*t_DetIds)[k]);
3070               double ener = cfac * (*t_HitEnergies)[k];
3071               if (corrPU_)
3072                 correctEnergy(ener, jentry);
3073               unsigned int idx = (unsigned int)((*t_DetIds)[k]);
3074               int subdet, zside, ieta, iphi, depth;
3075               unpackDetId(idx, subdet, zside, ieta, iphi, depth);
3076               if (depth > 0 && depth <= (int)(ndepth)) {
3077                 if (subdet == 1) {
3078                   eb += ener;
3079                   bv[depth - 1] += ener;
3080                   h_bvlist2[depth - 1]->Fill(ener, weight);
3081                   ++bnrec[depth - 1];
3082                 } else if (subdet == 2) {
3083                   ee += ener;
3084                   ev[depth - 1] += ener;
3085                   h_evlist2[depth - 1]->Fill(ener, weight);
3086                   ++enrec[depth - 1];
3087                 }
3088               }
3089             }
3090             bool barrel = (eb > ee);
3091             if (barrel)
3092               selHB += weight;
3093             else
3094               selHE += weight;
3095             for (unsigned int i = 0; i < ndepth; i++) {
3096               if (barrel) {
3097                 h_bvlist[i]->Fill(bv[i], weight);
3098                 h_bvlist3[i]->Fill((bnrec[i] + 0.001), weight);
3099               } else {
3100                 h_evlist[i]->Fill(ev[i], weight);
3101                 h_evlist3[i]->Fill((enrec[i] + 0.001), weight);
3102               }
3103             }
3104           }
3105         }
3106       }
3107       good++;
3108     }
3109     ++kount;
3110   }
3111   std::cout << "Finds " << duplicate << " Duplicate events out of " << kount << " events in this file and " << good
3112             << " selected events" << std::endl;
3113   if (plotHists_)
3114     std::cout << "Number of weighted selected events " << sel << " HB " << selHB << " HE " << selHE << std::endl;
3115 }
3116
3117 bool CalibPlotProperties::goodTrack(double &eHcal, double &cuti, const Long64_t &entry, bool debug) {
3118   bool select(true);
3119   double cut(cuti);
3120   if (debug) {
3121     std::cout << "goodTrack input " << eHcal << ":" << cut;
3122   }
3123   if (flexibleSelect_ > 1) {
3124     double eta = (t_ieta > 0) ? t_ieta : -t_ieta;
3125     cut = 8.0 * exp(eta * log2by18_);
3126   }
3127   correctEnergy(eHcal, entry);
3128   select = ((t_qltyFlag) && (t_selectTk) && (t_hmaxNearP < cut) && (t_eMipDR < 100.0) && (eHcal > 0.001));
3129   if (debug) {
3130     std::cout << " output " << eHcal << ":" << cut << ":" << select << std::endl;
3131   }
3132   return select;
3133 }
3134
3135 bool CalibPlotProperties::selectPhi(bool debug) {
3136   bool select(true);
3137   if (phimin_ > 1 || phimax_ < 72) {
3138     double eTotal(0), eSelec(0);
3139     // The masks are defined in DataFormats/HcalDetId/interface/HcalDetId.h
3140     for (unsigned int k = 0; k < t_HitEnergies->size(); ++k) {
3141       int iphi = ((*t_DetIds)[k]) & (0x3FF);
3142       int zside = ((*t_DetIds)[k] & 0x80000) ? (1) : (-1);
3143       eTotal += ((*t_HitEnergies)[k]);
3144       if (iphi >= phimin_ && iphi <= phimax_ && zside == zside_)
3145         eSelec += ((*t_HitEnergies)[k]);
3146     }
3147     if (eSelec < 0.9 * eTotal)
3148       select = false;
3149     if (debug) {
3150       std::cout << "Etotal " << eTotal << " and ESelec " << eSelec << " (phi " << phimin_ << ":" << phimax_ << " z "
3151                 << zside_ << ") Selection " << select << std::endl;
3152     }
3153   }
3154   return select;
3155 }
3156
3157 void CalibPlotProperties::savePlot(const std::string &theName, bool append, bool all) {
3158   TFile *theFile(0);
3159   if (append) {
3160     theFile = new TFile(theName.c_str(), "UPDATE");
3161   } else {
3162     theFile = new TFile(theName.c_str(), "RECREATE");
3163   }
3164
3165   theFile->cd();
3166
3167   if (plotBasic_) {
3168     for (int k = 0; k < 5; ++k) {
3169       if (h_p[k] != 0) {
3170         TH1D *h1 = (TH1D *)h_p[k]->Clone();
3171         h1->Write();
3172       }
3173       if (h_eta[k] != 0) {
3174         TH1D *h2 = (TH1D *)h_eta[k]->Clone();
3175         h2->Write();
3176       }
3177     }
3178     for (unsigned int k = 0; k < h_eta0.size(); ++k) {
3179       if (h_eta0[k] != 0) {
3180         TH1D *h1 = (TH1D *)h_eta0[k]->Clone();
3181         h1->Write();
3182       }
3183       if (h_eta1[k] != 0) {
3184         TH1D *h2 = (TH1D *)h_eta1[k]->Clone();
3185         h2->Write();
3186       }
3187       if (h_eta2[k] != 0) {
3188         TH1D *h3 = (TH1D *)h_eta2[k]->Clone();
3189         h3->Write();
3190       }
3191       if (h_eta3[k] != 0) {
3192         TH1D *h4 = (TH1D *)h_eta3[k]->Clone();
3193         h4->Write();
3194       }
3195       if (h_eta4[k] != 0) {
3196         TH1D *h5 = (TH1D *)h_eta4[k]->Clone();
3197         h5->Write();
3198       }
3199       if (h_dL1[k] != 0) {
3200         TH1D *h6 = (TH1D *)h_dL1[k]->Clone();
3201         h6->Write();
3202       }
3203       if (h_vtx[k] != 0) {
3204         TH1D *h7 = (TH1D *)h_vtx[k]->Clone();
3205         h7->Write();
3206       }
3207     }
3208   }
3209
3210   if (plotEnergy_) {
3211     for (unsigned int k = 0; k < ps_.size() - 1; ++k) {
3212       for (unsigned int j = 0; j <= (unsigned int)(etahi_ - etalo_); ++j) {
3213         if (h_etaEH[k].size() > j && h_etaEH[k][j] != nullptr && (all || (k == kp50))) {
3214           TH1D *hist = (TH1D *)h_etaEH[k][j]->Clone();
3215           hist->Write();
3216         }
3217         if (h_etaEp[k].size() > j && h_etaEp[k][j] != nullptr && (all || (k == kp50))) {
3218           TH1D *hist = (TH1D *)h_etaEp[k][j]->Clone();
3219           hist->Write();
3220         }
3221         if (h_etaEE[k].size() > j && h_etaEE[k][j] != nullptr && (all || (k == kp50))) {
3222           TH1D *hist = (TH1D *)h_etaEE[k][j]->Clone();
3223           hist->Write();
3224         }
3225       }
3226     }
3227
3228     for (unsigned int j = 0; j < ietas_.size(); ++j) {
3229       if (h_eHcal.size() > j && (h_eHcal[j] != nullptr)) {
3230         TH1D *hist = (TH1D *)h_eHcal[j]->Clone();
3231         hist->Write();
3232       }
3233       if (h_mom.size() > j && (h_mom[j] != nullptr)) {
3234         TH1D *hist = (TH1D *)h_mom[j]->Clone();
3235         hist->Write();
3236       }
3237       if (h_eEcal.size() > j && (h_eEcal[j] != nullptr)) {
3238         TH1D *hist = (TH1D *)h_eEcal[j]->Clone();
3239         hist->Write();
3240       }
3241     }
3242   }
3243
3244   if (plotHists_) {
3245     if (h_nvtx != 0) {
3246       TH1D *h1 = (TH1D *)h_nvtx->Clone();
3247       h1->Write();
3248     }
3249     if (h_etaE != 0) {
3250       TH2D *h2 = (TH2D *)h_etaE->Clone();
3251       h2->Write();
3252     }
3253     for (unsigned int i = 0; i < ndepth; ++i) {
3254       if (h_bvlist[i] != 0) {
3255         TH1D *h1 = (TH1D *)h_bvlist[i]->Clone();
3256         h1->Write();
3257       }
3258       if (h_bvlist2[i] != 0) {
3259         TH1D *h2 = (TH1D *)h_bvlist2[i]->Clone();
3260         h2->Write();
3261       }
3262       if (h_bvlist3[i] != 0) {
3263         TH1D *h3 = (TH1D *)h_bvlist3[i]->Clone();
3264         h3->Write();
3265       }
3266       if (h_evlist[i] != 0) {
3267         TH1D *h4 = (TH1D *)h_evlist[i]->Clone();
3268         h4->Write();
3269       }
3270       if (h_evlist2[i] != 0) {
3271         TH1D *h5 = (TH1D *)h_evlist2[i]->Clone();
3272         h5->Write();
3273       }
3274       if (h_evlist3[i] != 0) {
3275         TH1D *h6 = (TH1D *)h_evlist3[i]->Clone();
3276         h6->Write();
3277       }
3278     }
3279   }
3280   std::cout << "All done" << std::endl;
3281   theFile->Close();
3282 }
3283
3284 void CalibPlotProperties::correctEnergy(double &eHcal, const Long64_t &entry) {
3285   double pmom = (useGen_ && (t_gentrackP > 0)) ? t_gentrackP : t_p;
3286   if ((ifDepth_ == 3) && (cFactor_ != nullptr)) {
3287     double cfac = cFactor_->getCorr(entry);
3288     eHcal *= cfac;
3289   } else if ((corrPU_ < 0) && (pmom > 0)) {
3290     double ediff = (t_eHcal30 - t_eHcal10);
3291     if (t_DetIds1 != 0 && t_DetIds3 != 0) {
3292       double Etot1(0), Etot3(0);
3293       // The masks are defined in DataFormats/HcalDetId/interface/HcalDetId.h
3294       for (unsigned int idet = 0; idet < (*t_DetIds1).size(); idet++) {
3295         unsigned int id = truncateId((*t_DetIds1)[idet], truncateFlag_, false);
3296         double cfac = corrFactor_->getCorr(id);
3297         if (cFactor_ != 0)
3298           cfac *= cFactor_->getCorr(t_Run, (*t_DetIds1)[idet]);
3299         double hitEn = cfac * (*t_HitEnergies1)[idet];
3300         Etot1 += hitEn;
3301       }
3302       for (unsigned int idet = 0; idet < (*t_DetIds3).size(); idet++) {
3303         unsigned int id = truncateId((*t_DetIds3)[idet], truncateFlag_, false);
3304         double cfac = corrFactor_->getCorr(id);
3305         if (cFactor_ != 0)
3306           cfac *= cFactor_->getCorr(t_Run, (*t_DetIds3)[idet]);
3307         double hitEn = cfac * (*t_HitEnergies3)[idet];
3308         Etot3 += hitEn;
3309       }
3310       ediff = (Etot3 - Etot1);
3311     }
3312     double fac = puFactor(-corrPU_, t_ieta, pmom, eHcal, ediff);
3313     eHcal *= fac;
3314   } else if (corrPU_ > 0) {
3315     eHcal = puFactorRho(corrPU_, t_ieta, t_rhoh, eHcal);
3316   }
3317 }