Back to home page

Project CMSSW displayed by LXR
 
 

    

File indexing completed on 2023-03-17 10:42:55 

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

This page was automatically generated by the 2.2.1 LXR engine.

The LXR team
Valid CSS 2.1! Valid HTML 4.01!