Back to home page

Project CMSSW displayed by LXR

 
 

    


File indexing completed on 2023-01-31 03:06:13

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