Back to home page

Project CMSSW displayed by LXR

 
 

    


File indexing completed on 2021-02-14 12:48:07

0001 //////////////////////////////////////////////////////////////////////////////
0002 // This class analyzes the "Lep Tree" created by HBHEMuonOfflineAnalyzer
0003 // It has two constructors using either a pointer to the tree chain or
0004 // the file name where the tree resides
0005 // There are 2 additional arguments:
0006 //      mode (a packed integer) with bits specifying some action
0007 //           Bit 0 : 0 will produce set of histograms of energy deposit,
0008 //                     active length, active length corrected energy;
0009 //                   1 will produce plots of p, nvx and scatter plots for
0010 //                     each ieta
0011 //               1 : 0 ignores the information on number of vertex
0012 //                   1 groups ranges of # of vertex 0:15:20:25:30:100
0013 //               2 : 0 ignores the information on track momentum
0014 //                   1 separate plots for certain momentum range
0015 //                     (the range depends on ieta)
0016 //              3-4: 0 separate plot for each depth
0017 //                   1 sums up all depths
0018 //                   2 collapse depth
0019 //              5-6: 0 no check on iphi
0020 //                   1 separate plot for each iphi
0021 //                   2 separate plot for each RBX
0022 //                   3 exclude the RBX specified by bits 6-10
0023 //             7-11:   RBX # to be excluded (maximum 5 bits needed for RBX)
0024 //               12: 0 varying ranges of p depending on |ieta|
0025 //                   1 constant p ranges
0026 //            13-15: 0 no cut on ediff; 1-4 cuts at 5, 10, 15, 20 GeV
0027 //      modeLHC (integer) specifies the detector condition
0028 //              0      Run1   detector (till 2016)
0029 //              1      Plan36 detector (2018)
0030 //              2      Phase1 detector (Run3)
0031 //              3      Plan1  detector (2017)
0032 //              4      Phase2 detector
0033 //
0034 //   AnalyzeLepTree a1(tree, mode, modeLHC);
0035 //   AnalyzeLepTree a1(fname, mode, modeLHC);
0036 //   a1.Loop();
0037 //   a1.writeHisto(outfile);
0038 //   a1.plotHisto(drawStatBox,type,save);
0039 //   a1.writeMeanError(outfileMean);
0040 //
0041 //      outfile     (char*)  Name of the file where histograms to be written
0042 //      outfileMean (char*)  Name of the file where means with errors to be
0043 //                           written when it is run with mode bit 0 set to 1
0044 //      drawStatBox (bool)   If Statbox to be drawn or not
0045 //      type        (int)    Each bit says what plots to be drawn
0046 //                           If bit 0 of "mode" is 1
0047 //                           (0: momentum for each ieta;
0048 //                            1: number of vertex for each ieta;
0049 //                            2: 2D plot for nvx vs p for each ieta;
0050 //                            3: number of vertex for each ieta & depth;
0051 //                            4: momentum for each ieta & depth
0052 //                            5: energy in the outer ring)
0053 //                           If bit 0 of "mode" is 0
0054 //                           plots for all or specific depth & phi if
0055 //                           "depth" as of (type/16)&15 is 0 or the specified
0056 //                           value and "phi" as of (type/256)&127 is 0 or
0057 //                           the specified value. Bit 0 set plots the energy
0058 //                           distribution; 1 set plots active length corrected
0059 //                           energy; 2 set plots charge; 3 set plots active
0060 //                           length corrected charge distributions
0061 //      save        (bool)   If plots to be saved as pdf file or not
0062 ///////////////////////////////////////////////////////////////////////////////
0063 
0064 #include <TCanvas.h>
0065 #include <TChain.h>
0066 #include <TColor.h>
0067 #include <TFile.h>
0068 #include <TH2.h>
0069 #include <TPaveStats.h>
0070 #include <TPaveText.h>
0071 #include <TROOT.h>
0072 #include <TStyle.h>
0073 
0074 #include <cmath>
0075 #include <iomanip>
0076 #include <iostream>
0077 #include <fstream>
0078 #include <sstream>
0079 #include <map>
0080 #include <string>
0081 #include <vector>
0082 
0083 class AnalyzeLepTree {
0084 public:
0085   AnalyzeLepTree(TChain* tree, int mode = 0, int modeLHC = 1);
0086   AnalyzeLepTree(const char* fname, int mode = 0, int modeLHC = 1);
0087   virtual ~AnalyzeLepTree();
0088   virtual Int_t Cut(Long64_t entry);
0089   virtual Int_t GetEntry(Long64_t entry);
0090   virtual Long64_t LoadTree(Long64_t entry);
0091   virtual void Init(TChain* tree);
0092   virtual void Loop(Long64_t nmax = -1, bool debug = false);
0093   virtual Bool_t Notify();
0094   virtual void Show(Long64_t entry = -1);
0095   void writeHisto(const char* outfile);
0096   void writeMeanError(const char* outfile);
0097   std::vector<TCanvas*> plotHisto(bool drawStatBox, int type, bool save = false);
0098 
0099 private:
0100   bool fillChain(TChain* chain, const char* fname);
0101   void bookHisto();
0102   void plotHisto(std::map<unsigned int, TH1D*> hists, std::vector<TCanvas*>& cvs, bool save);
0103   void plotHisto(std::map<unsigned int, TH1D*> hists, int phi0, int depth0, std::vector<TCanvas*>& cvs, bool save);
0104   TCanvas* plotHisto(TH1D* hist);
0105   void plot2DHisto(std::map<unsigned int, TH2D*> hists, std::vector<TCanvas*>& cvs, bool save);
0106   int getCollapsedDepth(int eta, int phi, int depth);
0107   int getRBX(int eta);
0108   int getPBin(int eta);
0109   int getVxBin();
0110   int getDepthBin(int depth);
0111   int getPhiBin(int eta);
0112   void makeVxBins(int modeLHC);
0113   int nDepthBins(int eta, int phi, int modeLHC);
0114   int nPhiBins(int eta);
0115   int nPBins(int eta);
0116   int nVxBins();
0117   unsigned int packID(int zside, int eta, int phi, int depth, int nvx, int ipbin);
0118   void unpackID(unsigned int id, int& zside, int& eta, int& phi, int& depth, int& nvx, int& ipbin);
0119   void getBins(int type, int eta, int phi, int depth, int& nbins, double& xmax);
0120 
0121 private:
0122   TChain* fChain;  //!pointer to the analyzed TTree or TChain
0123   Int_t fCurrent;  //!current Tree number in a TChain
0124 
0125   // Declaration of leaf types
0126   Int_t t_ieta;
0127   Int_t t_iphi;
0128   Int_t t_nvtx;
0129   Double_t t_p;
0130   Double_t t_ediff;
0131   std::vector<double>* t_ene;
0132   std::vector<double>* t_enec;
0133   std::vector<double>* t_charge;
0134   std::vector<double>* t_actln;
0135   std::vector<int>* t_depth;
0136 
0137   // List of branches
0138   TBranch* b_t_ieta;    //!
0139   TBranch* b_t_iphi;    //!
0140   TBranch* b_t_nvtx;    //!
0141   TBranch* b_t_p;       //!
0142   TBranch* b_t_ediff;   //!
0143   TBranch* b_t_ene;     //!
0144   TBranch* b_t_enec;    //!
0145   TBranch* b_t_charge;  //!
0146   TBranch* b_t_actln;   //!
0147   TBranch* b_t_depth;   //!
0148 
0149   static const int etamax_ = 26, npbin_ = 9, nvbin_ = 6;
0150   static const bool debug_ = false;
0151   int mode_, modeLHC_, exRBX_, kphi_, kdepth_;
0152   std::vector<int> npvbin_, iprange_;
0153   std::vector<double> prange_[5];
0154   double cutEdiff_;
0155   std::map<unsigned int, TH1D*> h_p_, h_nv_;
0156   std::map<unsigned int, TH2D*> h_pnv_;
0157   std::map<unsigned int, TH1D*> h_p2_, h_nv2_;
0158   std::map<unsigned int, TH1D*> h_Energy_, h_Ecorr_, h_Charge_, h_Chcorr_;
0159   std::map<unsigned int, TH1D*> h_EnergyC_, h_EcorrC_;
0160   std::map<unsigned int, TH1D*> h_ediff_, h_ediff_nvtx_;
0161 };
0162 
0163 AnalyzeLepTree::AnalyzeLepTree(TChain* tree, int mode1, int mode2) : mode_(mode1), modeLHC_(mode2) {
0164   std::cout << "Proceed with a tree chain with " << tree->GetEntries() << " entries" << std::endl;
0165   Init(tree);
0166 }
0167 
0168 AnalyzeLepTree::AnalyzeLepTree(const char* fname, int mode1, int mode2) : mode_(mode1), modeLHC_(mode2) {
0169   TChain* chain = new TChain("Lep_Tree");
0170   if (!fillChain(chain, fname)) {
0171     std::cout << "*****No valid tree chain can be obtained*****" << std::endl;
0172   } else {
0173     std::cout << "Proceed with a tree chain with " << chain->GetEntries() << " entries" << std::endl;
0174     Init(chain);
0175   }
0176 }
0177 
0178 AnalyzeLepTree::~AnalyzeLepTree() {
0179   if (!fChain)
0180     return;
0181   delete fChain->GetCurrentFile();
0182 }
0183 
0184 Int_t AnalyzeLepTree::GetEntry(Long64_t entry) {
0185   // Read contents of entry.
0186   if (!fChain)
0187     return 0;
0188   return fChain->GetEntry(entry);
0189 }
0190 
0191 Long64_t AnalyzeLepTree::LoadTree(Long64_t entry) {
0192   // Set the environment to read one entry
0193   if (!fChain)
0194     return -5;
0195   Long64_t centry = fChain->LoadTree(entry);
0196   if (centry < 0)
0197     return centry;
0198   if (!fChain->InheritsFrom(TChain::Class()))
0199     return centry;
0200   TChain* chain = (TChain*)fChain;
0201   if (chain->GetTreeNumber() != fCurrent) {
0202     fCurrent = chain->GetTreeNumber();
0203     Notify();
0204   }
0205   return centry;
0206 }
0207 
0208 void AnalyzeLepTree::Init(TChain* tree) {
0209   // The Init() function is called when the selector needs to initialize
0210   // a new tree or chain. Typically here the branch addresses and branch
0211   // pointers of the tree will be set.
0212   // It is normally not necessary to make changes to the generated
0213   // code, but the routine can be extended by the user if needed.
0214   // Init() will be called many times when running on PROOF
0215   // (once per file to be processed).
0216 
0217   makeVxBins(modeLHC_);
0218   exRBX_ = (mode_ / 128) % 32;
0219   kphi_ = (mode_ / 32) % 4;
0220   kdepth_ = (mode_ / 8) % 4;
0221   if ((mode_ % 2) == 0)
0222     std::cout << "Produce set of histograms of energy, "
0223               << " active length, active length corrected "
0224               << "energy for 3 types" << std::endl;
0225   else
0226     std::cout << "Produce plots of p, nvx and scatter plots "
0227               << "for each ieta" << std::endl;
0228   if (((mode_ / 2) % 2) == 0) {
0229     std::cout << "Ignore the information on number of vertex iformation" << std::endl;
0230   } else {
0231     std::cout << "Group ranges of # of vertex ";
0232     for (unsigned int k = 0; k < npvbin_.size(); ++k)
0233       std::cout << npvbin_[k] << ":";
0234     std::cout << std::endl;
0235   }
0236   if (((mode_ / 4) % 2) == 0)
0237     std::cout << "Ignore the information on track "
0238               << "momentum" << std::endl;
0239   else
0240     std::cout << "Separate plots for certain momentum "
0241               << "range (the range depends on ieta)\n";
0242   if (kdepth_ == 0)
0243     std::cout << "Generate separate plot for each depth" << std::endl;
0244   else if (kdepth_ == 1)
0245     std::cout << "Sums up all depths for plots\n";
0246   else
0247     std::cout << "Collapse depths to Run 1 scenario\n";
0248   if (kphi_ == 0)
0249     std::cout << "Make no check on iphi" << std::endl;
0250   else if (kphi_ == 1)
0251     std::cout << "Make separate plot for each iphi\n";
0252   else if (kphi_ == 2)
0253     std::cout << "Make separate plot for each RBX\n";
0254   else
0255     std::cout << "Exclude the RBX " << exRBX_ << std::endl;
0256   if (modeLHC_ == 0)
0257     std::cout << "This is Run1 detector (till 2016)\n";
0258   else if (modeLHC_ == 1)
0259     std::cout << "This is Plan36 detector (2018)\n";
0260   else if (modeLHC_ == 2)
0261     std::cout << "This is Phase1 detector (after 2021)\n";
0262   else if (modeLHC_ == 3)
0263     std::cout << "This is Plan1  detector (2017)\n";
0264   else
0265     std::cout << "This is Phase2 detector (after 2026)\n";
0266   static const double cuts[8] = {200, 5, 10, 15, 20, 25, 30, 40};
0267   int cutE = (mode_ / 4096) % 8;
0268   cutEdiff_ = cuts[cutE];
0269   std::cout << "Cut off for energy in the 8 neighbouring towers " << cutEdiff_ << std::endl;
0270 
0271   // Set object pointer
0272   t_ene = 0;
0273   t_enec = 0;
0274   t_charge = 0;
0275   t_actln = 0;
0276   t_depth = 0;
0277   fChain = tree;
0278   // Set branch addresses and branch pointers
0279   if (!tree)
0280     return;
0281   fCurrent = -1;
0282   fChain->SetMakeClass(1);
0283 
0284   fChain->SetBranchAddress("t_ieta", &t_ieta, &b_t_ieta);
0285   fChain->SetBranchAddress("t_iphi", &t_iphi, &b_t_iphi);
0286   fChain->SetBranchAddress("t_nvtx", &t_nvtx, &b_t_nvtx);
0287   fChain->SetBranchAddress("t_p", &t_p, &b_t_p);
0288   fChain->SetBranchAddress("t_ediff", &t_ediff, &b_t_ediff);
0289   fChain->SetBranchAddress("t_ene", &t_ene, &b_t_ene);
0290   fChain->SetBranchAddress("t_enec", &t_enec, &b_t_enec);
0291   fChain->SetBranchAddress("t_charge", &t_charge, &b_t_charge);
0292   fChain->SetBranchAddress("t_actln", &t_actln, &b_t_actln);
0293   fChain->SetBranchAddress("t_depth", &t_depth, &b_t_depth);
0294   Notify();
0295 
0296   t_ediff = 0;
0297 
0298   bookHisto();
0299 }
0300 
0301 Bool_t AnalyzeLepTree::Notify() {
0302   // The Notify() function is called when a new file is opened. This
0303   // can be either for a new TTree in a TChain or when when a new TTree
0304   // is started when using PROOF. It is normally not necessary to make changes
0305   // to the generated code, but the routine can be extended by the
0306   // user if needed. The return value is currently not used.
0307   return kTRUE;
0308 }
0309 
0310 void AnalyzeLepTree::Show(Long64_t entry) {
0311   // Print contents of entry.
0312   // If entry is not specified, print current entry
0313   if (!fChain)
0314     return;
0315   fChain->Show(entry);
0316 }
0317 
0318 Int_t AnalyzeLepTree::Cut(Long64_t) {
0319   // This function may be called from Loop.
0320   // returns  1 if entry is accepted.
0321   // returns -1 otherwise.
0322   return 1;
0323 }
0324 
0325 void AnalyzeLepTree::Loop(Long64_t nmax, bool debug) {
0326   //   In a ROOT session, you can do:
0327   //      Root > .L AnalyzeLepTree.C
0328   //      Root > AnalyzeLepTree t
0329   //      Root > t.GetEntry(12); // Fill t data members with entry number 12
0330   //      Root > t.Show();       // Show values of entry 12
0331   //      Root > t.Show(16);     // Read and show values of entry 16
0332   //      Root > t.Loop();       // Loop on all entries
0333   //
0334 
0335   //     This is the loop skeleton where:
0336   //    jentry is the global entry number in the chain
0337   //    ientry is the entry number in the current Tree
0338   //  Note that the argument to GetEntry must be:
0339   //    jentry for TChain::GetEntry
0340   //    ientry for TTree::GetEntry and TBranch::GetEntry
0341   //
0342   //       To read only selected branches, Insert statements like:
0343   // METHOD1:
0344   //    fChain->SetBranchStatus("*",0);  // disable all branches
0345   //    fChain->SetBranchStatus("branchname",1);  // activate branchname
0346   // METHOD2: replace line
0347   //    fChain->GetEntry(jentry);       //read all branches
0348   //by  b_branchname->GetEntry(ientry); //read only this branch
0349   if (fChain == 0)
0350     return;
0351 
0352   Long64_t nentries = fChain->GetEntriesFast();
0353   std::cout << "Number of entries: " << nentries << ":" << nmax << std::endl;
0354   if (nmax > 0 && nmax < nentries)
0355     nentries = nmax;
0356   const double ethr = 0.00001;  // Threshold of 10 keV
0357 
0358   Long64_t nbytes = 0, nb = 0;
0359   int32_t n15(0), n16(0);
0360   for (Long64_t jentry = 0; jentry < nentries; jentry++) {
0361     Long64_t ientry = LoadTree(jentry);
0362     if ((jentry % 1000000 == 0) || debug)
0363       std::cout << "Entry " << jentry << ":" << ientry << std::endl;
0364     if (ientry < 0)
0365       break;
0366     nb = fChain->GetEntry(jentry);
0367     nbytes += nb;
0368     // if (Cut(ientry) < 0) continue;
0369     int zside = (t_ieta > 0) ? 1 : -1;
0370     int eta = (t_ieta > 0) ? t_ieta : -t_ieta;
0371     int phi = getPhiBin(eta);
0372     int pbin = getPBin(eta);
0373     int vbin = getVxBin();
0374     if ((mode_ / 1) % 2 == 1) {
0375       unsigned int id0 = packID(zside, eta, 1, 1, 1, 1);
0376       std::map<unsigned int, TH1D*>::iterator it1 = h_p_.find(id0);
0377       if (it1 != h_p_.end())
0378         (it1->second)->Fill(t_p);
0379       std::map<unsigned int, TH1D*>::iterator it2 = h_nv_.find(id0);
0380       if (it2 != h_nv_.end())
0381         (it2->second)->Fill(t_nvtx);
0382       std::map<unsigned int, TH2D*>::iterator it3 = h_pnv_.find(id0);
0383       if (it3 != h_pnv_.end())
0384         (it3->second)->Fill(t_nvtx, t_p);
0385       unsigned int id1 = packID(zside, eta, 1, 1, 1, pbin);
0386       std::map<unsigned int, TH1D*>::iterator it4 = h_p2_.find(id1);
0387       if (it4 != h_p2_.end())
0388         (it4->second)->Fill(t_p);
0389       unsigned int id2 = packID(zside, eta, 1, 1, vbin, 1);
0390       std::map<unsigned int, TH1D*>::iterator it5 = h_nv2_.find(id2);
0391       if (it5 != h_nv2_.end())
0392         (it5->second)->Fill(t_nvtx);
0393       std::map<unsigned int, TH1D*>::iterator it6 = h_ediff_.find(id0);
0394       if (it6 != h_ediff_.end())
0395         (it6->second)->Fill(t_ediff);
0396       std::map<unsigned int, TH1D*>::iterator it7 = h_ediff_nvtx_.find(id2);
0397       if (it7 != h_ediff_nvtx_.end())
0398         (it7->second)->Fill(t_ediff);
0399     } else {
0400       if (phi > 0 && pbin >= 0 && vbin >= 0) {
0401         if (kdepth_ == 0) {
0402           for (unsigned int k = 0; k < t_depth->size(); ++k) {
0403             if (eta == 15)
0404               ++n15;
0405             else if (eta == 16)
0406               ++n16;
0407             int depth = (*t_depth)[k];
0408             unsigned int id = packID(zside, eta, phi, depth + 1, vbin, pbin);
0409             double ene = (*t_ene)[k];
0410             double enec = (*t_enec)[k];
0411             double charge = (*t_charge)[k];
0412             double actl = (*t_actln)[k];
0413             if (ene > ethr && actl > 0 && charge > 0 && t_ediff < cutEdiff_) {
0414               std::map<unsigned int, TH1D*>::iterator it1 = h_Energy_.find(id);
0415               if (it1 != h_Energy_.end())
0416                 (it1->second)->Fill(ene);
0417               std::map<unsigned int, TH1D*>::iterator it2 = h_Ecorr_.find(id);
0418               if (it2 != h_Ecorr_.end())
0419                 (it2->second)->Fill(ene / actl);
0420               std::map<unsigned int, TH1D*>::iterator it3 = h_EnergyC_.find(id);
0421               if (it3 != h_EnergyC_.end())
0422                 (it3->second)->Fill(enec);
0423               std::map<unsigned int, TH1D*>::iterator it4 = h_EcorrC_.find(id);
0424               if (it4 != h_EcorrC_.end())
0425                 (it4->second)->Fill(enec / actl);
0426               std::map<unsigned int, TH1D*>::iterator it5 = h_Charge_.find(id);
0427               if (it5 != h_Charge_.end())
0428                 (it5->second)->Fill(charge);
0429               std::map<unsigned int, TH1D*>::iterator it6 = h_Chcorr_.find(id);
0430               if (it6 != h_Chcorr_.end())
0431                 (it6->second)->Fill(charge / actl);
0432               if (debug_) {
0433                 //      if ((eta>20 && (t_iphi > 35)) || (t_iphi > 71)) std::cout << zside << ":" << eta << ":" << phi << ":" << t_iphi << ":" << depth+1 << ":" << vbin << ":" << pbin << " ID " << std::hex << id << std::dec << " Flags " <<  (it1 != h_Energy_.end()) << ":" << (it2 != h_Ecorr_.end()) << ":" <<  (it3 != h_EnergyC_.end()) << ":" << (it4 != h_EcorrC_.end()) << ":" << (it5 != h_Charge_.end()) << ":" << (it6 != h_Chcorr_.end()) << " E " << ene << " C " << charge << " L " << actl << std::endl;
0434                 if ((it1 == h_Energy_.end()) || (it2 == h_Ecorr_.end()) || (it3 == h_EnergyC_.end()) ||
0435                     (it4 == h_EcorrC_.end()) || (it5 == h_Charge_.end()) || (it6 == h_Chcorr_.end()))
0436                   std::cout << zside << ":" << eta << ":" << phi << ":" << t_iphi << ":" << depth + 1 << ":" << vbin
0437                             << ":" << pbin << " ID " << std::hex << id << std::dec << " Flags "
0438                             << (it1 != h_Energy_.end()) << ":" << (it2 != h_Ecorr_.end()) << ":"
0439                             << (it3 != h_Charge_.end()) << ":" << (it4 != h_Chcorr_.end()) << " E " << ene << " C "
0440                             << charge << " L " << actl << std::endl;
0441               }
0442             }
0443           }
0444         } else if (kdepth_ == 1) {
0445           double ene[2], enec[2], actl[2], charge[2];
0446           for (unsigned int k = 0; k < 2; ++k) {
0447             ene[k] = enec[k] = actl[k] = charge[k] = 0;
0448           }
0449           for (unsigned int k = 0; k < t_depth->size(); ++k) {
0450             if ((*t_ene)[k] > 0 && (*t_actln)[k] > 0) {
0451               int dep = (*t_depth)[k];
0452               int depth = ((eta != 16) ? 0 : ((dep > 1) ? 1 : 0));
0453               ene[depth] += (*t_ene)[k];
0454               enec[depth] += (*t_enec)[k];
0455               charge[depth] += (*t_charge)[k];
0456               actl[depth] += (*t_actln)[k];
0457             }
0458           }
0459           int nDepth = (eta == 16) ? 2 : 1;
0460           for (int k = 0; k < nDepth; ++k) {
0461             if (ene[k] > ethr && actl[k] > 0 && charge[k] > 0 && t_ediff < cutEdiff_) {
0462               if (eta == 15)
0463                 ++n15;
0464               else if (eta == 16)
0465                 ++n16;
0466               int depth = k + 1;
0467               unsigned int id = packID(zside, eta, phi, depth, vbin, pbin);
0468               std::map<unsigned int, TH1D*>::iterator it1 = h_Energy_.find(id);
0469               if (it1 != h_Energy_.end())
0470                 (it1->second)->Fill(ene[k]);
0471               std::map<unsigned int, TH1D*>::iterator it2 = h_Ecorr_.find(id);
0472               if (it2 != h_Ecorr_.end())
0473                 (it2->second)->Fill(ene[k] / actl[k]);
0474               std::map<unsigned int, TH1D*>::iterator it3 = h_EnergyC_.find(id);
0475               if (it3 != h_EnergyC_.end())
0476                 (it3->second)->Fill(enec[k]);
0477               std::map<unsigned int, TH1D*>::iterator it4 = h_EcorrC_.find(id);
0478               if (it4 != h_EcorrC_.end())
0479                 (it4->second)->Fill(enec[k] / actl[k]);
0480               std::map<unsigned int, TH1D*>::iterator it5 = h_Charge_.find(id);
0481               if (it5 != h_Charge_.end())
0482                 (it5->second)->Fill(charge[k]);
0483               std::map<unsigned int, TH1D*>::iterator it6 = h_Chcorr_.find(id);
0484               if (it6 != h_Chcorr_.end())
0485                 (it6->second)->Fill(charge[k] / actl[k]);
0486               if (((eta == 15) || (eta == 16)) && debug_)
0487                 std::cout << zside << ":" << eta << ":" << phi << ":" << t_iphi << ":" << depth << ":" << vbin << ":"
0488                           << pbin << " ID " << std::hex << id << std::dec << " Flags " << (it1 != h_Energy_.end())
0489                           << ":" << (it2 != h_Ecorr_.end()) << ":" << (it3 != h_Charge_.end()) << ":"
0490                           << (it4 != h_Chcorr_.end()) << " E " << ene << " C " << charge << " L " << actl << std::endl;
0491             }
0492           }
0493         } else {
0494           double ene[3], enec[3], actl[3], charge[3];
0495           for (unsigned int k = 0; k < 3; ++k) {
0496             ene[k] = enec[k] = actl[k] = charge[k] = 0;
0497           }
0498           for (unsigned int k = 0; k < t_depth->size(); ++k) {
0499             if ((*t_ene)[k] > 0 && (*t_actln)[k] > 0) {
0500               int dep = (*t_depth)[k];
0501               int depth = getCollapsedDepth(zside * eta, phi, dep) - 1;
0502               ene[depth] += (*t_ene)[k];
0503               enec[depth] += (*t_enec)[k];
0504               charge[depth] += (*t_charge)[k];
0505               actl[depth] += (*t_actln)[k];
0506             }
0507           }
0508           for (int k = 0; k < nDepthBins(eta, phi, 0); ++k) {
0509             if (ene[k] > ethr && actl[k] > 0 && charge[k] > 0 && t_ediff < cutEdiff_) {
0510               if (eta == 15)
0511                 ++n15;
0512               else if (eta == 16)
0513                 ++n16;
0514               int depth = k + 1;
0515               unsigned int id = packID(zside, eta, phi, depth, vbin, pbin);
0516               std::map<unsigned int, TH1D*>::iterator it1 = h_Energy_.find(id);
0517               if (it1 != h_Energy_.end())
0518                 (it1->second)->Fill(ene[k]);
0519               std::map<unsigned int, TH1D*>::iterator it2 = h_Ecorr_.find(id);
0520               if (it2 != h_Ecorr_.end())
0521                 (it2->second)->Fill(ene[k] / actl[k]);
0522               std::map<unsigned int, TH1D*>::iterator it3 = h_EnergyC_.find(id);
0523               if (it3 != h_EnergyC_.end())
0524                 (it3->second)->Fill(enec[k]);
0525               std::map<unsigned int, TH1D*>::iterator it4 = h_EcorrC_.find(id);
0526               if (it4 != h_EcorrC_.end())
0527                 (it4->second)->Fill(enec[k] / actl[k]);
0528               std::map<unsigned int, TH1D*>::iterator it5 = h_Charge_.find(id);
0529               if (it5 != h_Charge_.end())
0530                 (it5->second)->Fill(charge[k]);
0531               std::map<unsigned int, TH1D*>::iterator it6 = h_Chcorr_.find(id);
0532               if (it6 != h_Chcorr_.end())
0533                 (it6->second)->Fill(charge[k] / actl[k]);
0534               if (((eta == 15) || (eta == 16)) && debug_)
0535                 std::cout << zside << ":" << eta << ":" << phi << ":" << t_iphi << ":" << depth << ":" << vbin << ":"
0536                           << pbin << " ID " << std::hex << id << std::dec << " Flags " << (it1 != h_Energy_.end())
0537                           << ":" << (it2 != h_Ecorr_.end()) << ":" << (it3 != h_Charge_.end()) << ":"
0538                           << (it4 != h_Chcorr_.end()) << " E " << ene[k] << " C " << charge[k] << " L " << actl[k]
0539                           << std::endl;
0540             }
0541           }
0542         }
0543       }
0544     }
0545   }
0546   std::cout << "Number of events with eta15: " << n15 << " and eta16: " << n16 << std::endl;
0547 }
0548 
0549 bool AnalyzeLepTree::fillChain(TChain* chain, const char* inputFileList) {
0550   int kount(0);
0551   std::string fname(inputFileList);
0552   if (fname.substr(fname.size() - 5, 5) == ".root") {
0553     chain->Add(fname.c_str());
0554   } else {
0555     ifstream infile(inputFileList);
0556     if (!infile.is_open()) {
0557       std::cout << "** ERROR: Can't open '" << inputFileList << "' for input" << std::endl;
0558       return false;
0559     }
0560     while (1) {
0561       infile >> fname;
0562       if (!infile.good())
0563         break;
0564       chain->Add(fname.c_str());
0565       ++kount;
0566     }
0567     infile.close();
0568   }
0569   std::cout << "Adds " << kount << " files in the chain from " << fname << std::endl;
0570   return true;
0571 }
0572 
0573 void AnalyzeLepTree::bookHisto() {
0574   for (int i = 0; i < 5; ++i)
0575     prange_[i].clear();
0576   int ipbrng[30] = {0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 1, 1, 1, 1, 2, 2, 2, 2, 2, 2, 3, 3, 3, 4, 4, 4, 4};
0577   for (int i = 0; i < 30; ++i)
0578     iprange_.push_back(ipbrng[i]);
0579   double prange0[npbin_] = {0, 30, 45, 55, 75, 100, 125, 150, 500};
0580   double prange1[npbin_] = {0, 50, 75, 100, 125, 150, 200, 300, 500};
0581   double prange2[npbin_] = {0, 60, 75, 100, 125, 150, 200, 300, 500};
0582   double prange3[npbin_] = {0, 100, 125, 150, 175, 200, 300, 400, 500};
0583   double prange4[npbin_] = {0, 125, 150, 175, 200, 250, 300, 400, 500};
0584   double prangeX[npbin_] = {125, 150, 200, 250, 300, 350, 400, 500};
0585   for (int i = 0; i < npbin_; ++i) {
0586     if ((mode_ / 4096) % 2 == 0) {
0587       prange_[0].push_back(prange0[i]);
0588       prange_[1].push_back(prange1[i]);
0589       prange_[2].push_back(prange2[i]);
0590       prange_[3].push_back(prange3[i]);
0591       prange_[4].push_back(prange4[i]);
0592     } else {
0593       prange_[0].push_back(prangeX[i]);
0594       prange_[1].push_back(prangeX[i]);
0595       prange_[2].push_back(prangeX[i]);
0596       prange_[3].push_back(prangeX[i]);
0597       prange_[4].push_back(prangeX[i]);
0598     }
0599   }
0600   if (debug_) {
0601     std::cout << "Eta range " << -etamax_ << ":" << etamax_ << " # of vtx bins " << nVxBins() << std::endl;
0602     if ((mode_ / 1) % 2 == 0) {
0603       for (int ieta = -etamax_; ieta <= etamax_; ++ieta) {
0604         int eta = (ieta > 0) ? ieta : -ieta;
0605         if (eta != 0) {
0606           int ndepth = ((kdepth_ == 0) ? nDepthBins(eta, 63, modeLHC_)
0607                                        : ((kdepth_ != 1) ? nDepthBins(eta, 63, 0) : (eta == 16) ? 2 : 1));
0608           std::cout << "Eta " << ieta << " with " << nPhiBins(eta) << " phi bins " << ndepth << " maximum depths and "
0609                     << nPBins(eta) << " p bins" << std::endl;
0610         }
0611       }
0612     }
0613   }
0614 
0615   char name[100], title[200];
0616   unsigned int book1(0), book2(0);
0617   if ((mode_ / 1) % 2 == 1) {
0618     h_p_.clear();
0619     h_nv_.clear();
0620     h_pnv_.clear();
0621     h_nv2_.clear();
0622     h_p2_.clear();
0623     h_ediff_.clear();
0624     h_ediff_nvtx_.clear();
0625     for (int ieta = -etamax_; ieta <= etamax_; ++ieta) {
0626       if (ieta != 0) {
0627         int zside = (ieta > 0) ? 1 : -1;
0628         int eta = (ieta > 0) ? ieta : -ieta;
0629         unsigned int id0 = packID(zside, eta, 1, 1, 1, 1);
0630         sprintf(name, "peta%d", ieta);
0631         sprintf(title, "Momentum for i#eta = %d (GeV)", ieta);
0632         h_p_[id0] = new TH1D(name, title, 100, 0.0, 500.0);
0633         ++book1;
0634         sprintf(name, "Ediff_eta%d", ieta);
0635         sprintf(title, "Energy difference for i#eta = %d (GeV)", ieta);
0636         h_ediff_[id0] = new TH1D(name, title, 1000, -500.0, 500.0);
0637         ++book1;
0638         sprintf(name, "nveta%d", ieta);
0639         sprintf(title, "Number of Vertex for i#eta = %d", ieta);
0640         h_nv_[id0] = new TH1D(name, title, 100, 0.0, 100.0);
0641         ++book1;
0642         sprintf(name, "pnveta%d", ieta);
0643         sprintf(title, "i#eta = %d", ieta);
0644         TH2D* h2 = new TH2D(name, title, 100, 0.0, 100.0, 100, 0.0, 500.0);
0645         ++book2;
0646         h2->GetXaxis()->SetTitle("Number of Vertex");
0647         h2->GetYaxis()->SetTitle("Momentum (GeV)");
0648         h_pnv_[id0] = h2;
0649         ++book1;
0650         char etas[10];
0651         sprintf(etas, "i#eta=%d", ieta);
0652         char name[100], title[200];
0653         for (int pbin = 0; pbin < nPBins(eta); ++pbin) {
0654           char ps[20];
0655           if ((mode_ / 4) % 2 == 1) {
0656             int np = (eta >= 0 && eta < (int)(iprange_.size())) ? iprange_[eta] - 1 : iprange_[0];
0657             sprintf(ps, "p=%d:%d", (int)prange_[np][pbin], (int)prange_[np][pbin + 1]);
0658           };
0659           unsigned int id = packID(zside, eta, 1, 1, 1, pbin);
0660           sprintf(name, "pvx%d111%d", ieta, pbin);
0661           sprintf(title, "Momentum for %s %s", etas, ps);
0662           h_p2_[id] = new TH1D(name, title, 100, 0.0, 500.0);
0663           h_p2_[id]->Sumw2();
0664           ++book1;
0665         }
0666         for (int vbin = 0; vbin < nVxBins(); ++vbin) {
0667           char vtx[12];
0668           if ((mode_ / 2) % 2 == 1) {
0669             sprintf(vtx, "N_{vtx}=%d:%d", npvbin_[vbin], npvbin_[vbin + 1]);
0670           } else {
0671             sprintf(vtx, "all N_{vtx}");
0672           }
0673           unsigned int id = packID(zside, eta, 1, 1, vbin, 1);
0674           sprintf(name, "nvx%d11%d1", ieta, vbin);
0675           sprintf(title, "Number of vertex for %s %s", etas, vtx);
0676           h_nv2_[id] = new TH1D(name, title, 100, 0.0, 100.0);
0677           h_nv2_[id]->Sumw2();
0678           ++book1;
0679           sprintf(name, "Ediff_nvx%d11%d1", ieta, vbin);
0680           sprintf(title, "Energy difference for %s %s", etas, vtx);
0681           h_ediff_nvtx_[id] = new TH1D(name, title, 1000, -500.0, 500.0);
0682           h_ediff_nvtx_[id]->Sumw2();
0683           ++book1;
0684         }
0685       }
0686     }
0687   } else {
0688     h_Energy_.clear();
0689     h_Ecorr_.clear();
0690     h_Charge_.clear();
0691     h_Chcorr_.clear();
0692     h_EnergyC_.clear();
0693     h_EcorrC_.clear();
0694     for (int ieta = -etamax_; ieta <= etamax_; ++ieta) {
0695       if (ieta != 0) {
0696         int zside = (ieta > 0) ? 1 : -1;
0697         int eta = (ieta > 0) ? ieta : -ieta;
0698         char etas[20];
0699         sprintf(etas, "i#eta=%d", ieta);
0700         for (int iphi = 0; iphi < nPhiBins(eta); ++iphi) {
0701           char phis[20];
0702           int phi(1), phi0(63);
0703           if (kphi_ == 1) {
0704             phi0 = phi = (eta <= 20) ? iphi + 1 : 2 * iphi + 1;
0705             sprintf(phis, "i#phi=%d", phi);
0706           } else if (kphi_ == 2) {
0707             phi0 = 4 * iphi + 1;
0708             phi = iphi + 1;
0709             sprintf(phis, "RBX=%d", iphi + 1);
0710           } else if (kphi_ == 3) {
0711             sprintf(phis, "All except RBX %d", exRBX_);
0712           } else {
0713             sprintf(phis, "All i#phi");
0714           }
0715           int ndepth = ((kdepth_ == 0) ? nDepthBins(eta, phi0, modeLHC_)
0716                                        : ((kdepth_ != 1) ? nDepthBins(eta, phi0, 0) : (eta == 16) ? 2 : 1));
0717           for (int depth = 0; depth < ndepth; ++depth) {
0718             char deps[20];
0719             if (kdepth_ == 1) {
0720               if (depth == 0)
0721                 sprintf(deps, "all depths");
0722               else
0723                 sprintf(deps, "all endcap depths");
0724             } else {
0725               sprintf(deps, "Depth=%d", depth + 1);
0726             }
0727             for (int pbin = 0; pbin < nPBins(eta); ++pbin) {
0728               char ps[30];
0729               if ((mode_ / 4) % 2 == 1) {
0730                 int np = (eta >= 0 && eta < (int)(iprange_.size())) ? iprange_[eta] - 1 : iprange_[0];
0731                 sprintf(ps, "p=%d:%d", (int)prange_[np][pbin], (int)prange_[np][pbin + 1]);
0732               } else {
0733                 sprintf(ps, "all p");
0734               };
0735               for (int vbin = 0; vbin < nVxBins(); ++vbin) {
0736                 int nbin(4000);
0737                 double xmax(10.0);
0738                 char vtx[20];
0739                 if ((mode_ / 2) % 2 == 1) {
0740                   sprintf(vtx, "N_{vtx}=%d:%d", npvbin_[vbin], npvbin_[vbin + 1]);
0741                 } else {
0742                   sprintf(vtx, "all N_{vtx}");
0743                 }
0744                 unsigned int id = packID(zside, eta, phi, depth + 1, vbin, pbin);
0745                 char name[100], title[200];
0746                 sprintf(name, "EdepE%dF%dD%dV%dP%d", ieta, phi, depth, vbin, pbin);
0747                 sprintf(title, "Deposited energy for %s %s %s %s %s (GeV)", etas, phis, deps, ps, vtx);
0748                 getBins(0, ieta, phi0, depth + 1, nbin, xmax);
0749                 h_Energy_[id] = new TH1D(name, title, nbin, 0.0, xmax);
0750                 ++book1;
0751                 sprintf(name, "EcorE%dF%dD%dV%dP%d", ieta, phi, depth, vbin, pbin);
0752                 sprintf(title, "Active length corrected energy for %s %s %s %s %s (GeV/cm)", etas, phis, deps, ps, vtx);
0753                 getBins(1, ieta, phi0, depth + 1, nbin, xmax);
0754                 h_Ecorr_[id] = new TH1D(name, title, nbin, 0.0, xmax);
0755                 ++book1;
0756                 sprintf(name, "EdepCE%dF%dD%dV%dP%d", ieta, phi, depth, vbin, pbin);
0757                 sprintf(
0758                     title, "Response Corrected deposited energy for %s %s %s %s %s (GeV)", etas, phis, deps, ps, vtx);
0759                 getBins(2, ieta, phi0, depth + 1, nbin, xmax);
0760                 h_EnergyC_[id] = new TH1D(name, title, nbin, 0.0, xmax);
0761                 ++book1;
0762                 sprintf(name, "EcorCE%dF%dD%dV%dP%d", ieta, phi, depth, vbin, pbin);
0763                 sprintf(title,
0764                         "Response and active length corrected energy for %s %s %s %s %s (GeV/cm)",
0765                         etas,
0766                         phis,
0767                         deps,
0768                         ps,
0769                         vtx);
0770                 getBins(3, ieta, phi0, depth + 1, nbin, xmax);
0771                 h_EcorrC_[id] = new TH1D(name, title, nbin, 0.0, xmax);
0772                 ++book1;
0773                 sprintf(name, "ChrgE%dF%dD%dV%dP%d", ieta, phi, depth, vbin, pbin);
0774                 sprintf(title, "Measured charge for %s %s %s %s %s (cm)", etas, phis, deps, ps, vtx);
0775                 getBins(4, ieta, phi0, depth + 1, nbin, xmax);
0776                 h_Charge_[id] = new TH1D(name, title, nbin, 0.0, xmax);
0777                 ++book1;
0778                 sprintf(name, "ChcorE%dF%dD%dV%dP%d", ieta, phi, depth, vbin, pbin);
0779                 sprintf(title, "Active length corrected charge for %s %s %s %s %s (cm)", etas, phis, deps, ps, vtx);
0780                 getBins(5, ieta, phi0, depth + 1, nbin, xmax);
0781                 h_Chcorr_[id] = new TH1D(name, title, nbin, 0.0, xmax);
0782                 ++book1;
0783               }
0784             }
0785           }
0786         }
0787       }
0788     }
0789   }
0790   std::cout << "Booked " << book1 << " 1D and " << book2 << " 2D Histos\n";
0791 }
0792 
0793 void AnalyzeLepTree::writeHisto(const char* outfile) {
0794   TFile* output_file = TFile::Open(outfile, "RECREATE");
0795   output_file->cd();
0796   if ((mode_ / 1) % 2 == 1) {
0797     for (std::map<unsigned int, TH1D*>::const_iterator itr = h_p_.begin(); itr != h_p_.end(); ++itr)
0798       (itr->second)->Write();
0799     for (std::map<unsigned int, TH1D*>::const_iterator itr = h_nv_.begin(); itr != h_nv_.end(); ++itr)
0800       (itr->second)->Write();
0801     for (std::map<unsigned int, TH2D*>::const_iterator itr = h_pnv_.begin(); itr != h_pnv_.end(); ++itr)
0802       (itr->second)->Write();
0803     for (std::map<unsigned int, TH1D*>::const_iterator itr = h_p2_.begin(); itr != h_p2_.end(); ++itr)
0804       (itr->second)->Write();
0805     for (std::map<unsigned int, TH1D*>::const_iterator itr = h_nv2_.begin(); itr != h_nv2_.end(); ++itr)
0806       (itr->second)->Write();
0807     for (std::map<unsigned int, TH1D*>::const_iterator itr = h_ediff_.begin(); itr != h_ediff_.end(); ++itr)
0808       (itr->second)->Write();
0809     for (std::map<unsigned int, TH1D*>::const_iterator itr = h_ediff_nvtx_.begin(); itr != h_ediff_nvtx_.end(); ++itr)
0810       (itr->second)->Write();
0811   } else {
0812     for (int ieta = -etamax_; ieta <= etamax_; ++ieta) {
0813       if (ieta != 0) {
0814         char dirname[50];
0815         sprintf(dirname, "DieMuonEta%d", ieta);
0816         TDirectory* d_output = output_file->mkdir(dirname);
0817         d_output->cd();
0818         int zside = (ieta > 0) ? 1 : -1;
0819         int eta = (ieta > 0) ? ieta : -ieta;
0820         for (int iphi = 0; iphi < nPhiBins(eta); ++iphi) {
0821           int phi(1), phi0(1);
0822           if (kphi_ == 1) {
0823             phi0 = phi = (eta <= 20) ? iphi + 1 : 2 * iphi + 1;
0824           } else if (kphi_ == 2) {
0825             phi0 = 4 * iphi + 1;
0826             phi = iphi + 1;
0827           };
0828           int ndepth = ((kdepth_ == 0) ? nDepthBins(eta, phi0, modeLHC_)
0829                                        : ((kdepth_ != 1) ? nDepthBins(eta, phi0, 0) : (eta == 16) ? 2 : 1));
0830           for (int depth = 0; depth < ndepth; ++depth) {
0831             for (int pbin = 0; pbin < nPBins(eta); ++pbin) {
0832               for (int vbin = 0; vbin < nVxBins(); ++vbin) {
0833                 unsigned int id = packID(zside, eta, phi, depth + 1, vbin, pbin);
0834                 std::map<unsigned int, TH1D*>::const_iterator itr;
0835                 itr = h_Energy_.find(id);
0836                 if (itr != h_Energy_.end())
0837                   (itr->second)->Write();
0838                 itr = h_Ecorr_.find(id);
0839                 if (itr != h_Ecorr_.end())
0840                   (itr->second)->Write();
0841                 itr = h_EnergyC_.find(id);
0842                 if (itr != h_EnergyC_.end())
0843                   (itr->second)->Write();
0844                 itr = h_EcorrC_.find(id);
0845                 if (itr != h_EcorrC_.end())
0846                   (itr->second)->Write();
0847                 itr = h_Charge_.find(id);
0848                 if (itr != h_Charge_.end())
0849                   (itr->second)->Write();
0850                 itr = h_Chcorr_.find(id);
0851                 if (itr != h_Chcorr_.end())
0852                   (itr->second)->Write();
0853               }
0854             }
0855           }
0856         }
0857       }
0858       output_file->cd();
0859     }
0860   }
0861 }
0862 
0863 void AnalyzeLepTree::writeMeanError(const char* outfile) {
0864   if ((mode_ / 1) % 2 == 1) {
0865     std::ofstream fOutput(outfile);
0866     for (int vbin = 0; vbin < nVxBins(); ++vbin) {
0867       for (int ieta = -etamax_; ieta <= etamax_; ++ieta) {
0868         if (ieta != 0) {
0869           int zside = (ieta > 0) ? 1 : -1;
0870           int eta = (ieta > 0) ? ieta : -ieta;
0871           unsigned int id = packID(zside, eta, 1, 1, vbin, 1);
0872           char name[100];
0873           sprintf(name, "nvx%d11%d1", ieta, vbin);
0874           std::map<unsigned int, TH1D*>::iterator itr = h_nv2_.find(id);
0875           if (itr != h_nv2_.end()) {
0876             double mean = (itr->second)->GetMean();
0877             double error = (itr->second)->GetMeanError();
0878             char vtx[12];
0879             if ((mode_ / 2) % 2 == 1) {
0880               sprintf(vtx, "Nvtx=%3d:%3d", npvbin_[vbin], npvbin_[vbin + 1]);
0881             } else {
0882               sprintf(vtx, "all Nvtx");
0883             }
0884             fOutput << "Eta " << std::setw(3) << ieta << " " << vtx << " " << mean << " +- " << error << std::endl;
0885           }
0886         }
0887       }
0888     }
0889     fOutput.close();
0890   }
0891 }
0892 
0893 std::vector<TCanvas*> AnalyzeLepTree::plotHisto(bool drawStatBox, int type, bool save) {
0894   std::vector<TCanvas*> cvs;
0895   gStyle->SetCanvasBorderMode(0);
0896   gStyle->SetCanvasColor(kWhite);
0897   gStyle->SetPadColor(kWhite);
0898   gStyle->SetFillColor(kWhite);
0899   gStyle->SetOptTitle(0);
0900   if (drawStatBox) {
0901     gStyle->SetOptStat(111110);
0902     gStyle->SetOptFit(1);
0903   } else {
0904     gStyle->SetOptStat(0);
0905     gStyle->SetOptFit(0);
0906   }
0907 
0908   if ((mode_ / 1) % 2 == 1) {
0909     if (type % 2 > 0)
0910       plotHisto(h_p_, cvs, save);
0911     if ((type / 2) % 2 > 0)
0912       plotHisto(h_nv_, cvs, save);
0913     if ((type / 4) % 2 > 0)
0914       plot2DHisto(h_pnv_, cvs, save);
0915     if ((type / 8) % 2 > 0)
0916       plotHisto(h_nv2_, cvs, save);
0917     if ((type / 16) % 2 > 0)
0918       plotHisto(h_p2_, cvs, save);
0919     if ((type / 32) % 2 > 0)
0920       plotHisto(h_ediff_, cvs, save);
0921     if ((type / 32) % 2 > 0)
0922       plotHisto(h_ediff_nvtx_, cvs, save);
0923   } else {
0924     int depth0 = (type / 16) & 15;
0925     int phi0 = (type / 256) & 127;
0926     bool doEn = ((type / 1) % 2 > 0);
0927     bool doEnL = ((type / 2) % 2 > 0);
0928     bool doChg = ((type / 4) % 2 > 0);
0929     bool doChL = ((type / 8) % 2 > 0);
0930     if (doEn)
0931       plotHisto(h_Energy_, phi0, depth0, cvs, save);
0932     if (doEn)
0933       plotHisto(h_EnergyC_, phi0, depth0, cvs, save);
0934     if (doEnL)
0935       plotHisto(h_Ecorr_, phi0, depth0, cvs, save);
0936     if (doEnL)
0937       plotHisto(h_EcorrC_, phi0, depth0, cvs, save);
0938     if (doChg)
0939       plotHisto(h_Charge_, phi0, depth0, cvs, save);
0940     if (doChL)
0941       plotHisto(h_Chcorr_, phi0, depth0, cvs, save);
0942   }
0943   return cvs;
0944 }
0945 
0946 void AnalyzeLepTree::plotHisto(std::map<unsigned int, TH1D*> hists, std::vector<TCanvas*>& cvs, bool save) {
0947   for (std::map<unsigned int, TH1D*>::const_iterator itr = hists.begin(); itr != hists.end(); ++itr) {
0948     TH1D* hist = (itr->second);
0949     if (hist != 0) {
0950       TCanvas* pad = plotHisto(hist);
0951       cvs.push_back(pad);
0952       if (save) {
0953         char name[100];
0954         sprintf(name, "c_%s.pdf", pad->GetName());
0955         pad->Print(name);
0956       }
0957     }
0958   }
0959 }
0960 
0961 void AnalyzeLepTree::plotHisto(
0962     std::map<unsigned int, TH1D*> hists, int phi0, int depth0, std::vector<TCanvas*>& cvs, bool save) {
0963   for (std::map<unsigned int, TH1D*>::const_iterator itr = hists.begin(); itr != hists.end(); ++itr) {
0964     int zside, eta, phi, depth, pbin, vbin;
0965     unpackID(itr->first, zside, eta, phi, depth, vbin, pbin);
0966     TH1D* hist = itr->second;
0967     if (((depth == depth0) || (depth0 == 0)) && ((phi == phi0) || (phi0 == 0)) && (hist != 0)) {
0968       TCanvas* pad = plotHisto(hist);
0969       cvs.push_back(pad);
0970       if (save) {
0971         char name[100];
0972         sprintf(name, "c_%s.pdf", pad->GetName());
0973         pad->Print(name);
0974       }
0975     }
0976   }
0977 }
0978 
0979 TCanvas* AnalyzeLepTree::plotHisto(TH1D* hist) {
0980   TCanvas* pad = new TCanvas(hist->GetName(), hist->GetName(), 700, 500);
0981   pad->SetRightMargin(0.10);
0982   pad->SetTopMargin(0.10);
0983   hist->GetXaxis()->SetTitleSize(0.04);
0984   hist->GetXaxis()->SetTitle(hist->GetTitle());
0985   hist->GetYaxis()->SetTitle("Tracks");
0986   hist->GetYaxis()->SetLabelOffset(0.005);
0987   hist->GetYaxis()->SetTitleSize(0.04);
0988   hist->GetYaxis()->SetLabelSize(0.035);
0989   hist->GetYaxis()->SetTitleOffset(1.30);
0990   hist->SetMarkerStyle(20);
0991   hist->SetLineWidth(2);
0992   hist->Draw();
0993   pad->Update();
0994   TPaveStats* st1 = (TPaveStats*)hist->GetListOfFunctions()->FindObject("stats");
0995   if (st1 != NULL) {
0996     st1->SetY1NDC(0.70);
0997     st1->SetY2NDC(0.90);
0998     st1->SetX1NDC(0.65);
0999     st1->SetX2NDC(0.90);
1000   }
1001   pad->Modified();
1002   pad->Update();
1003   return pad;
1004 }
1005 
1006 void AnalyzeLepTree::plot2DHisto(std::map<unsigned int, TH2D*> hists, std::vector<TCanvas*>& cvs, bool save) {
1007   char name[100];
1008   for (std::map<unsigned int, TH2D*>::const_iterator itr = hists.begin(); itr != hists.end(); ++itr) {
1009     TH2D* hist = (itr->second);
1010     if (hist != 0) {
1011       TCanvas* pad = new TCanvas(hist->GetName(), hist->GetName(), 700, 700);
1012       pad->SetRightMargin(0.10);
1013       pad->SetTopMargin(0.10);
1014       hist->GetXaxis()->SetTitleSize(0.04);
1015       hist->GetYaxis()->SetLabelOffset(0.005);
1016       hist->GetYaxis()->SetTitleSize(0.04);
1017       hist->GetYaxis()->SetLabelSize(0.035);
1018       hist->GetYaxis()->SetTitleOffset(1.30);
1019       hist->Draw("COLZ");
1020       pad->Update();
1021       TPaveStats* st1 = (TPaveStats*)hist->GetListOfFunctions()->FindObject("stats");
1022       if (st1 != NULL) {
1023         st1->SetY1NDC(0.65);
1024         st1->SetY2NDC(0.90);
1025         st1->SetX1NDC(0.65);
1026         st1->SetX2NDC(0.90);
1027       }
1028       pad->Modified();
1029       pad->Update();
1030       cvs.push_back(pad);
1031       if (save) {
1032         sprintf(name, "c_%s.pdf", pad->GetName());
1033         pad->Print(name);
1034       }
1035     }
1036   }
1037 }
1038 
1039 int AnalyzeLepTree::getCollapsedDepth(int eta, int phi, int dep) {
1040   int depth = dep + 1;
1041   int ieta = (eta > 0) ? eta : -eta;
1042   if (ieta <= 14 || ieta == 17) {
1043     depth = 1;
1044   } else if (ieta == 15) {
1045     if (modeLHC_ > 3) {
1046       if (dep > 3)
1047         depth = 2;
1048       else
1049         depth = 1;
1050     }
1051   } else if (ieta == 16) {
1052     if (modeLHC_ == 0 || (modeLHC_ == 3 && (phi < 63 || phi > 66 || eta < 0))) {
1053     } else {
1054       if (dep > 2)
1055         depth = 3;
1056     }
1057   } else if (ieta < 26) {
1058     if (modeLHC_ == 0 || (modeLHC_ == 3 && (phi < 63 || phi > 66 || eta < 0))) {
1059     } else {
1060       if (dep < 3)
1061         depth = 1;
1062       else
1063         depth = 2;
1064     }
1065   } else if (ieta == 26) {
1066     if (modeLHC_ == 0 || (modeLHC_ == 3 && (phi < 63 || phi > 66 || eta < 0))) {
1067     } else {
1068       if (dep < 4)
1069         depth = 1;
1070       else
1071         depth = 2;
1072     }
1073   } else {
1074     if (modeLHC_ == 0 || (modeLHC_ == 3 && (phi < 63 || phi > 66 || eta < 0))) {
1075     } else {
1076       if (dep < 3)
1077         depth = 1;
1078       else if (dep == 3)
1079         depth = 2;
1080       else
1081         depth = 3;
1082     }
1083   }
1084   return depth;
1085 }
1086 
1087 int AnalyzeLepTree::getRBX(int eta) {
1088   int rbx(1);
1089   int phi = (eta > 20) ? (2 * t_iphi + 1) : (t_iphi + 1);
1090   if (phi > 2 && phi < 71)
1091     rbx = (phi + 5) / 4;
1092   return rbx;
1093 }
1094 
1095 int AnalyzeLepTree::getPBin(int eta) {
1096   int bin(0);
1097   if ((mode_ / 4) % 2 == 1) {
1098     int np = (eta >= 0 && eta < (int)(iprange_.size())) ? iprange_[eta] : iprange_[0];
1099     for (unsigned int k = 1; k < prange_[np].size(); ++k) {
1100       if (t_p < prange_[np][k])
1101         break;
1102       bin = k;
1103     }
1104   }
1105   return bin;
1106 }
1107 
1108 int AnalyzeLepTree::getVxBin() {
1109   int bin(0);
1110   if ((mode_ / 2) % 2 == 1) {
1111     for (unsigned int k = 1; k < npvbin_.size(); ++k) {
1112       if (t_nvtx < npvbin_[k])
1113         break;
1114       bin = k;
1115     }
1116   }
1117   return bin;
1118 }
1119 
1120 int AnalyzeLepTree::getDepthBin(int depth) {
1121   int bin = (kdepth_ == 0) ? depth : 1;
1122   return bin;
1123 }
1124 
1125 int AnalyzeLepTree::getPhiBin(int eta) {
1126   int bin(1);
1127   if (kphi_ == 1) {
1128     bin = (eta > 20) ? (2 * t_iphi + 1) : (t_iphi + 1);
1129   } else if (kphi_ == 2) {
1130     bin = getRBX(eta);
1131   } else if (kphi_ == 3) {
1132     if (exRBX_ == getRBX(eta))
1133       bin = -1;
1134   }
1135   return bin;
1136 }
1137 
1138 void AnalyzeLepTree::makeVxBins(int modeLHC) {
1139   int npvbin0[nvbin_] = {0, 15, 20, 25, 30, 100};
1140   int npvbin1[nvbin_] = {0, 20, 25, 30, 35, 100};
1141   int npvbin2[nvbin_] = {0, 25, 30, 35, 40, 100};
1142   int npvbin3[nvbin_] = {0, 30, 40, 50, 70, 200};
1143   npvbin_.clear();
1144   for (int i = 0; i < nvbin_; ++i) {
1145     if (modeLHC <= 0 || modeLHC == 3)
1146       npvbin_.push_back(npvbin0[i]);
1147     else if (modeLHC == 1)
1148       npvbin_.push_back(npvbin1[i]);
1149     else if (modeLHC == 2)
1150       npvbin_.push_back(npvbin2[i]);
1151     else
1152       npvbin_.push_back(npvbin3[i]);
1153   }
1154 }
1155 
1156 int AnalyzeLepTree::nDepthBins(int eta, int phi, int modeLHC) {
1157   // Run 1 scenario
1158   int nDepthR1[29] = {1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 2, 3, 1, 2, 2, 2, 2, 2, 2, 2, 2, 2, 3, 3, 2};
1159   // Run 2 scenario from 2018
1160   int nDepthR2[29] = {1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 2, 4, 3, 5, 6, 6, 6, 6, 6, 6, 6, 7, 7, 7, 3};
1161   // Run 3 scenario
1162   int nDepthR3[29] = {4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 3, 5, 6, 6, 6, 6, 6, 6, 6, 7, 7, 7, 3};
1163   // Run 4 scenario
1164   int nDepthR4[29] = {4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 7, 7, 7, 7, 7, 7, 7, 7, 7, 7, 7, 7, 7};
1165   // modeLHC = 0 -->      corresponds to Run 1 (valid till 2016)
1166   //         = 1 -->      corresponds to Run 2 (2018 geometry)
1167   //         = 2 -->      corresponds to Run 3 (post LS2)
1168   //         = 3 -->      corresponds to 2017  (Plan 1)
1169   //         = 4 -->      corresponds to Run4  (post LS3)
1170   int nbin(0);
1171   if (modeLHC == 0) {
1172     nbin = nDepthR1[eta - 1];
1173   } else if (modeLHC == 1) {
1174     nbin = nDepthR2[eta - 1];
1175   } else if (modeLHC == 2) {
1176     nbin = nDepthR3[eta - 1];
1177   } else if (modeLHC == 3) {
1178     if (phi > 0) {
1179       if (eta >= 16 && phi >= 63 && phi <= 66) {
1180         nbin = nDepthR2[eta - 1];
1181       } else {
1182         nbin = nDepthR1[eta - 1];
1183       }
1184     } else {
1185       if (eta >= 16) {
1186         nbin = (nDepthR2[eta - 1] > nDepthR1[eta - 1]) ? nDepthR2[eta - 1] : nDepthR1[eta - 1];
1187       } else {
1188         nbin = nDepthR1[eta - 1];
1189       }
1190     }
1191   } else {
1192     if (eta > 0 && eta < 30) {
1193       nbin = nDepthR4[eta - 1];
1194     } else {
1195       nbin = nDepthR4[28];
1196     }
1197   }
1198   return nbin;
1199 }
1200 
1201 int AnalyzeLepTree::nPhiBins(int eta) {
1202   int nphi = (eta <= 20) ? 72 : 36;
1203   if (modeLHC_ == 4 && eta > 16)
1204     nphi = 360;
1205   if (kphi_ == 0)
1206     nphi = 1;
1207   else if (kphi_ == 2)
1208     nphi = 18;
1209   else if (kphi_ == 3)
1210     nphi = 1;
1211   return nphi;
1212 }
1213 
1214 int AnalyzeLepTree::nPBins(int eta) {
1215   int bin(1);
1216   if ((mode_ / 4) % 2 == 1) {
1217     int np = (eta >= 0 && eta < (int)(iprange_.size())) ? iprange_[eta] - 1 : iprange_[0];
1218     bin = (int)(prange_[np].size()) - 1;
1219   }
1220   return bin;
1221 }
1222 
1223 int AnalyzeLepTree::nVxBins() {
1224   int nvx = ((mode_ / 2) % 2 == 1) ? ((int)(npvbin_.size()) - 1) : 1;
1225   return nvx;
1226 }
1227 
1228 unsigned int AnalyzeLepTree::packID(int zside, int eta, int phi, int depth, int nvx, int ipbin) {
1229   unsigned int id = (zside > 0) ? 1 : 0;
1230   id |= (((nvx & 7) << 19) | ((ipbin & 7) << 16) | ((depth & 7) << 13) | ((eta & 31) << 8) | ((phi & 127) << 1));
1231   return id;
1232 }
1233 
1234 void AnalyzeLepTree::unpackID(unsigned int id, int& zside, int& eta, int& phi, int& depth, int& nvx, int& ipbin) {
1235   zside = (id % 2 == 0) ? -1 : 1;
1236   phi = (id >> 1) & 127;
1237   eta = (id >> 8) & 31;
1238   depth = (id >> 13) & 7;
1239   ipbin = (id >> 16) & 7;
1240   nvx = (id >> 19) & 7;
1241 }
1242 
1243 void AnalyzeLepTree::getBins(int type, int ieta, int phi, int depth, int& nbin, double& xmax) {
1244   int eta = (ieta >= 0) ? ieta : -ieta;
1245   bool barrel = (eta < 16) || ((eta == 16) && (depth <= 2));
1246   bool rbx17 = (phi >= 63) && (phi <= 66) && (ieta >= 16) && (!barrel);
1247   nbin = 5000;
1248   xmax = 10.0;
1249   if (type >= 4) {
1250     if ((modeLHC_ == 0) || (((modeLHC_ == 1) || (modeLHC_ == 3)) && barrel) || ((modeLHC_ == 3) && (!rbx17))) {
1251       // HPD Channels
1252       xmax = 50.0;
1253     } else {
1254       // SiPM Channels
1255       xmax = 10000.0;
1256       nbin = 10000;
1257     }
1258   }
1259 }