Back to home page

Project CMSSW displayed by LXR

 
 

    


File indexing completed on 2025-01-04 00:29:13

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