Back to home page

Project CMSSW displayed by LXR

 
 

    


File indexing completed on 2024-04-06 11:58:54

0001 ///////////////////////////////////////////////////////////////////////////////
0002 //
0003 //  HBHEMuonOfflineAnalyzer h1(tree, outfile, rcorFile, flag, mode, maxDHB,
0004 //                             maxDHE, cutMu, cutP, nevMax, over, runLo,
0005 //                             runHi, etaMin, etaMax, debug);
0006 //  HBHEMuonOfflineAnalyzer h1(infile, outfile, rcorFile, flag, mode, maxDHB,
0007 //                             maxDHE, cutMu, cutP, nevMax, over, runLo,
0008 //                             runHi, etaMin, etaMax, debug);
0009 //   h1.Loop()
0010 //
0011 //   tree       TTree*       Pointer to the tree chain
0012 //   infile     const char*  Name of the input file
0013 //   outfile    const char*  Name of the output file
0014 //                           (dyll_PU20_25_output_10.root)
0015 //   rcorFile   consr char*  name of the text file having the correction factors
0016 //                           as a function of run numbers to be used for raddam
0017 //                           correction (default="", no corr.)
0018 //   flag       int          Flag of 2 digits ("to"): to where "o" decides if
0019 //                           corrected (1) or default (0) energy to be used;
0020 //                           "t" decides if all depths to be merged (1) or not
0021 //                           (0) (default is 0)
0022 //   mode       int          Geometry file used 0:(defined by maxDHB/HE);
0023 //                           1 (Run 1; valid till 2016); 2 (Run 2; 2018);
0024 //                           3 (Run 3; post LS2); 4 (2017 Plan 1);
0025 //                           5 (Run 4; post LS3); default (3)
0026 //   maxDHB     int          Maximum number of depths for HB (4)
0027 //   maxDHE     int          Maximum number of depths for HE (7)
0028 //   cutMu      int          Selection of muon type:
0029 //                           (tight:0; medium:1; loose:2) default (0)
0030 //   cutP       float        Minimum muon momentum; default (10)
0031 //   nevMax     int          Maximum # oe entries to be processed; -1 means
0032 //                           all entries (-1)
0033 //   over       int          Override some of the selection
0034 //                           (0: not to override; 1: override) default (0)
0035 //   runLO      int          Minimum run number (1)
0036 //   runHI      int          Maximum run number (99999999)
0037 //   etaMin     int          Minimum (absolute) eta value (1)
0038 //   etaMax     int          Maximum (absolute) eta value (29)
0039 //
0040 ///////////////////////////////////////////////////////////////////////////////
0041 
0042 #include <algorithm>
0043 #include <cmath>
0044 #include <fstream>
0045 #include <iostream>
0046 #include <sstream>
0047 #include <vector>
0048 #include <string>
0049 #include <TCanvas.h>
0050 #include <TChain.h>
0051 #include <TFile.h>
0052 #include <TH1.h>
0053 #include <TH1D.h>
0054 #include <TH2.h>
0055 #include <TH2D.h>
0056 #include <TProfile.h>
0057 #include <TROOT.h>
0058 #include <TStyle.h>
0059 #include <TString.h>
0060 #include <TTree.h>
0061 
0062 class HBHEMuonOfflineAnalyzer {
0063 public:
0064   TChain *fChain;  //!pointer to the analyzed TTree/TChain
0065   Int_t fCurrent;  //!current Tree number in a TChain
0066 
0067   // Fixed size dimensions of array or collections stored in the TTree if any.
0068   // Declaration of leaf types
0069 
0070   UInt_t Event_No;
0071   UInt_t Run_No;
0072   UInt_t LumiNumber;
0073   UInt_t BXNumber;
0074   UInt_t GoodVertex;
0075   bool PF_Muon;
0076   bool Global_Muon;
0077   bool Tracker_muon;
0078   bool MuonIsTight;
0079   bool MuonIsMedium;
0080   double pt_of_muon;
0081   double eta_of_muon;
0082   double phi_of_muon;
0083   double energy_of_muon;
0084   double p_of_muon;
0085   float muon_trkKink;
0086   float muon_chi2LocalPosition;
0087   float muon_segComp;
0088   int TrackerLayer;
0089   int NumPixelLayers;
0090   int InnerTrackPixelHits;
0091   bool innerTrack;
0092   double chiTracker;
0093   double DxyTracker;
0094   double DzTracker;
0095   double innerTrackpt;
0096   double innerTracketa;
0097   double innerTrackphi;
0098   double tight_validFraction;
0099   bool OuterTrack;
0100   double OuterTrackPt;
0101   double OuterTrackEta;
0102   double OuterTrackPhi;
0103   double OuterTrackChi;
0104   int OuterTrackHits;
0105   int OuterTrackRHits;
0106   bool GlobalTrack;
0107   double GlobalTrckPt;
0108   double GlobalTrckEta;
0109   double GlobalTrckPhi;
0110   int Global_Muon_Hits;
0111   int MatchedStations;
0112   double GlobTrack_Chi;
0113   double Tight_LongitudinalImpactparameter;
0114   double Tight_TransImpactparameter;
0115   double IsolationR04;
0116   double IsolationR03;
0117   double ecal_3into3;
0118   double hcal_3into3;
0119   double tracker_3into3;
0120   bool matchedId;
0121   bool hcal_cellHot;
0122   double ecal_3x3;
0123   double hcal_1x1;
0124   unsigned int ecal_detID;
0125   unsigned int hcal_detID;
0126   unsigned int ehcal_detID;
0127   int hcal_ieta;
0128   int hcal_iphi;
0129   double hcal_edepth1;
0130   double hcal_activeL1;
0131   double hcal_edepthHot1;
0132   double hcal_activeHotL1;
0133   double hcal_edepthCorrect1;
0134   double hcal_edepthHotCorrect1;
0135   double hcal_cdepthHot1;
0136   double hcal_cdepthHotBG1;
0137   bool hcal_depthMatch1;
0138   bool hcal_depthMatchHot1;
0139   double hcal_edepth2;
0140   double hcal_activeL2;
0141   double hcal_edepthHot2;
0142   double hcal_activeHotL2;
0143   double hcal_edepthCorrect2;
0144   double hcal_edepthHotCorrect2;
0145   double hcal_cdepthHot2;
0146   double hcal_cdepthHotBG2;
0147   bool hcal_depthMatch2;
0148   bool hcal_depthMatchHot2;
0149   double hcal_edepth3;
0150   double hcal_activeL3;
0151   double hcal_edepthHot3;
0152   double hcal_activeHotL3;
0153   double hcal_edepthCorrect3;
0154   double hcal_edepthHotCorrect3;
0155   double hcal_cdepthHot3;
0156   double hcal_cdepthHotBG3;
0157   bool hcal_depthMatch3;
0158   bool hcal_depthMatchHot3;
0159   double hcal_edepth4;
0160   double hcal_activeL4;
0161   double hcal_edepthHot4;
0162   double hcal_activeHotL4;
0163   double hcal_edepthCorrect4;
0164   double hcal_edepthHotCorrect4;
0165   double hcal_cdepthHot4;
0166   double hcal_cdepthHotBG4;
0167   bool hcal_depthMatch4;
0168   bool hcal_depthMatchHot4;
0169   double hcal_edepth5;
0170   double hcal_activeL5;
0171   double hcal_edepthHot5;
0172   double hcal_activeHotL5;
0173   double hcal_edepthCorrect5;
0174   double hcal_edepthHotCorrect5;
0175   double hcal_cdepthHot5;
0176   double hcal_cdepthHotBG5;
0177   bool hcal_depthMatch5;
0178   bool hcal_depthMatchHot5;
0179   double hcal_edepth6;
0180   double hcal_activeL6;
0181   double hcal_edepthHot6;
0182   double hcal_activeHotL6;
0183   double hcal_edepthCorrect6;
0184   double hcal_edepthHotCorrect6;
0185   double hcal_cdepthHot6;
0186   double hcal_cdepthHotBG6;
0187   bool hcal_depthMatch6;
0188   bool hcal_depthMatchHot6;
0189   double hcal_edepth7;
0190   double hcal_activeL7;
0191   double hcal_edepthHot7;
0192   double hcal_activeHotL7;
0193   double hcal_edepthCorrect7;
0194   double hcal_edepthHotCorrect7;
0195   double hcal_cdepthHot7;
0196   double hcal_cdepthHotBG7;
0197   bool hcal_depthMatch7;
0198   bool hcal_depthMatchHot7;
0199   double activeLength;
0200   double activeLengthHot;
0201   std::vector<int> *hltresults;
0202   std::vector<std::string> *all_triggers;
0203 
0204   // List of branches
0205   TBranch *b_Event_No;                           //!
0206   TBranch *b_Run_No;                             //!
0207   TBranch *b_LumiNumber;                         //!
0208   TBranch *b_BXNumber;                           //!
0209   TBranch *b_GoodVertex;                         //!
0210   TBranch *b_PF_Muon;                            //!
0211   TBranch *b_Global_Muon;                        //!
0212   TBranch *b_Tracker_muon;                       //!
0213   TBranch *b_MuonIsTight;                        //!
0214   TBranch *b_MuonIsMedium;                       //!
0215   TBranch *b_pt_of_muon;                         //!
0216   TBranch *b_eta_of_muon;                        //!
0217   TBranch *b_phi_of_muon;                        //!
0218   TBranch *b_energy_of_muon;                     //!
0219   TBranch *b_p_of_muon;                          //!
0220   TBranch *b_muon_trkKink;                       //!
0221   TBranch *b_muon_chi2LocalPosition;             //!
0222   TBranch *b_muon_segComp;                       //!
0223   TBranch *b_TrackerLayer;                       //!
0224   TBranch *b_NumPixelLayers;                     //!
0225   TBranch *b_InnerTrackPixelHits;                //!
0226   TBranch *b_innerTrack;                         //!
0227   TBranch *b_chiTracker;                         //!
0228   TBranch *b_DxyTracker;                         //!
0229   TBranch *b_DzTracker;                          //!
0230   TBranch *b_innerTrackpt;                       //!
0231   TBranch *b_innerTracketa;                      //!
0232   TBranch *b_innerTrackphi;                      //!
0233   TBranch *b_tight_validFraction;                //!
0234   TBranch *b_OuterTrack;                         //!
0235   TBranch *b_OuterTrackPt;                       //!
0236   TBranch *b_OuterTrackEta;                      //!
0237   TBranch *b_OuterTrackPhi;                      //!
0238   TBranch *b_OuterTrackChi;                      //!
0239   TBranch *b_OuterTrackHits;                     //!
0240   TBranch *b_OuterTrackRHits;                    //!
0241   TBranch *b_GlobalTrack;                        //!
0242   TBranch *b_GlobalTrckPt;                       //!
0243   TBranch *b_GlobalTrckEta;                      //!
0244   TBranch *b_GlobalTrckPhi;                      //!
0245   TBranch *b_Global_Muon_Hits;                   //!
0246   TBranch *b_MatchedStations;                    //!
0247   TBranch *b_GlobTrack_Chi;                      //!
0248   TBranch *b_Tight_LongitudinalImpactparameter;  //!
0249   TBranch *b_Tight_TransImpactparameter;         //!
0250   TBranch *b_IsolationR04;                       //!
0251   TBranch *b_IsolationR03;                       //!
0252   TBranch *b_ecal_3into3;                        //!
0253   TBranch *b_hcal_3into3;                        //!
0254   TBranch *b_tracker_3into3;                     //!
0255   TBranch *b_matchedId;                          //!
0256   TBranch *b_hcal_cellHot;                       //!
0257   TBranch *b_ecal_3x3;                           //!
0258   TBranch *b_hcal_1x1;                           //!
0259   TBranch *b_hcal_detID;                         //!
0260   TBranch *b_ecal_detID;                         //!
0261   TBranch *b_ehcal_detID;                        //!
0262   TBranch *b_hcal_ieta;                          //!
0263   TBranch *b_hcal_iphi;                          //!
0264   TBranch *b_hcal_edepth1;                       //!
0265   TBranch *b_hcal_activeL1;                      //!
0266   TBranch *b_hcal_edepthHot1;                    //!
0267   TBranch *b_hcal_activeHotL1;                   //!
0268   TBranch *b_hcal_edepthCorrect1;                //!
0269   TBranch *b_hcal_edepthHotCorrect1;             //!
0270   TBranch *b_hcal_cdepthHot1;                    //!
0271   TBranch *b_hcal_cdepthHotBG1;                  //!
0272   TBranch *b_hcal_depthMatch1;                   //!
0273   TBranch *b_hcal_depthMatchHot1;                //!
0274   TBranch *b_hcal_edepth2;                       //!
0275   TBranch *b_hcal_activeL2;                      //!
0276   TBranch *b_hcal_edepthHot2;                    //!
0277   TBranch *b_hcal_activeHotL2;                   //!
0278   TBranch *b_hcal_edepthCorrect2;                //!
0279   TBranch *b_hcal_edepthHotCorrect2;             //!
0280   TBranch *b_hcal_cdepthHot2;                    //!
0281   TBranch *b_hcal_cdepthHotBG2;                  //!
0282   TBranch *b_hcal_depthMatch2;                   //!
0283   TBranch *b_hcal_depthMatchHot2;                //!
0284   TBranch *b_hcal_edepth3;                       //!
0285   TBranch *b_hcal_activeL3;                      //!
0286   TBranch *b_hcal_edepthHot3;                    //!
0287   TBranch *b_hcal_activeHotL3;                   //!
0288   TBranch *b_hcal_edepthCorrect3;                //!
0289   TBranch *b_hcal_edepthHotCorrect3;             //!
0290   TBranch *b_hcal_cdepthHot3;                    //!
0291   TBranch *b_hcal_cdepthHotBG3;                  //!
0292   TBranch *b_hcal_depthMatch3;                   //!
0293   TBranch *b_hcal_depthMatchHot3;                //!
0294   TBranch *b_hcal_edepth4;                       //!
0295   TBranch *b_hcal_activeL4;                      //!
0296   TBranch *b_hcal_edepthHot4;                    //!
0297   TBranch *b_hcal_activeHotL4;                   //!
0298   TBranch *b_hcal_edepthCorrect4;                //!
0299   TBranch *b_hcal_edepthHotCorrect4;             //!
0300   TBranch *b_hcal_cdepthHot4;                    //!
0301   TBranch *b_hcal_cdepthHotBG4;                  //!
0302   TBranch *b_hcal_depthMatch4;                   //!
0303   TBranch *b_hcal_depthMatchHot4;                //!
0304   TBranch *b_hcal_edepth5;                       //!
0305   TBranch *b_hcal_activeL5;                      //!
0306   TBranch *b_hcal_edepthHot5;                    //!
0307   TBranch *b_hcal_activeHotL5;                   //!
0308   TBranch *b_hcal_edepthCorrect5;                //!
0309   TBranch *b_hcal_edepthHotCorrect5;             //!
0310   TBranch *b_hcal_cdepthHot5;                    //!
0311   TBranch *b_hcal_cdepthHotBG5;                  //!
0312   TBranch *b_hcal_depthMatch5;                   //!
0313   TBranch *b_hcal_depthMatchHot5;                //!
0314   TBranch *b_hcal_edepth6;                       //!
0315   TBranch *b_hcal_activeL6;                      //!
0316   TBranch *b_hcal_edepthHot6;                    //!
0317   TBranch *b_hcal_activeHotL6;                   //!
0318   TBranch *b_hcal_edepthCorrect6;                //!
0319   TBranch *b_hcal_edepthHotCorrect6;             //!
0320   TBranch *b_hcal_cdepthHot6;                    //!
0321   TBranch *b_hcal_cdepthHotBG6;                  //!
0322   TBranch *b_hcal_depthMatch6;                   //!
0323   TBranch *b_hcal_depthMatchHot6;                //!
0324   TBranch *b_hcal_edepth7;                       //!
0325   TBranch *b_hcal_activeL7;                      //!
0326   TBranch *b_hcal_edepthHot7;                    //!
0327   TBranch *b_hcal_activeHotL7;                   //!
0328   TBranch *b_hcal_edepthCorrect7;                //!
0329   TBranch *b_hcal_edepthHotCorrect7;             //!
0330   TBranch *b_hcal_cdepthHot7;                    //!
0331   TBranch *b_hcal_cdepthHotBG7;                  //!
0332   TBranch *b_hcal_depthMatch7;                   //!
0333   TBranch *b_hcal_depthMatchHot7;                //!
0334   TBranch *b_activeLength;                       //!
0335   TBranch *b_activeLengthHot;                    //!
0336   TBranch *b_hltresults;                         //!
0337   TBranch *b_all_triggers;                       //!
0338 
0339   HBHEMuonOfflineAnalyzer(TChain *tree = 0,
0340                           const char *outfile = "dyll_PU20_25_output_10.root",
0341                           const char *rcorFileName = "",
0342                           int flag = 0,
0343                           int mode = 3,
0344                           int maxDHB = 4,
0345                           int maxDHE = 7,
0346                           int cutMu = 0,
0347                           float cutP = 5,
0348                           int nevMax = -1,
0349                           int over = 0,
0350                           int runLo = 1,
0351                           int runHi = 99999999,
0352                           int etaMin = 1,
0353                           int etaMax = 29,
0354                           bool debug = false);
0355   HBHEMuonOfflineAnalyzer(const char *infile,
0356                           const char *outfile = "dyll_PU20_25_output_10.root",
0357                           const char *rcorFileName = "",
0358                           int flag = 0,
0359                           int mode = 3,
0360                           int maxDHB = 4,
0361                           int maxDHE = 7,
0362                           int cutMu = 0,
0363                           float cutP = 5,
0364                           int nevMax = -1,
0365                           int over = 0,
0366                           int runLo = 1,
0367                           int runHi = 99999999,
0368                           int etaMin = 1,
0369                           int etaMax = 29,
0370                           bool debug = false);
0371   // mode of LHC is kept 1 for 2017 scenario as no change in depth segmentation
0372   // mode of LHC is 3 for 2021
0373   virtual ~HBHEMuonOfflineAnalyzer();
0374 
0375   virtual Int_t Cut(Long64_t entry);
0376   virtual Int_t GetEntry(Long64_t entry);
0377   virtual Long64_t LoadTree(Long64_t entry);
0378   virtual void Init(TChain *tree,
0379                     const char *rcorFileName,
0380                     int flag,
0381                     int mode,
0382                     int maxDHB,
0383                     int maxDHE,
0384                     int runLo,
0385                     int runHi,
0386                     int etaMin,
0387                     int etaMax);
0388   virtual void Loop();
0389   virtual Bool_t Notify();
0390   virtual void Show(Long64_t entry = -1);
0391 
0392   bool fillChain(TChain *chain, const char *inputFileList);
0393   bool readCorr(const char *rcorFileName);
0394   void bookHistograms(const char *);
0395   bool getEnergy(int dep, double &enb, double &enu, double &enh, double &enc, double &chgS, double &chgB, double &actL);
0396   void writeHistograms();
0397   bool looseMuon();
0398   bool tightMuon();
0399   bool softMuon();
0400   bool mediumMuon2016();
0401   void etaPhiHcal(unsigned int detId, int &eta, int &phi, int &depth);
0402   void etaPhiEcal(unsigned int detId, int &type, int &zside, int &etaX, int &phiY, int &plane, int &strip);
0403   void calculateP(double pt, double eta, double &pM);
0404   void close();
0405   int nDepthBins(int ieta, int iphi);
0406   int nPhiBins(int ieta);
0407   float getCorr(int run, unsigned int id);
0408   std::vector<std::string> splitString(const std::string &);
0409   unsigned int getDetIdHBHE(int ieta, int iphi, int depth);
0410   unsigned int getDetId(int subdet, int ieta, int iphi, int depth);
0411   unsigned int correctDetId(const unsigned int &detId);
0412   void unpackDetId(unsigned int detId, int &subdet, int &zside, int &ieta, int &iphi, int &depth);
0413 
0414 private:
0415   static const int maxDep_ = 7;
0416   static const int maxEta_ = 29;
0417   static const int maxPhi_ = 72;
0418   //3x16x72x2 + 5x4x72x2 + 5x9x36x2
0419   static const int maxHist_ = 20000;  //13032;
0420   static const unsigned int nmax_ = 10;
0421   int nCut_;
0422   const double cutP_;
0423   const int nevMax_;
0424   const bool over_, debug_;
0425   int modeLHC_, maxDepthHB_, maxDepthHE_, maxDepth_;
0426   int runLo_, runHi_, etaMin_, etaMax_;
0427   bool cFactor_, useCorrect_, mergeDepth_;
0428   int nHist, nDepths[maxEta_], nDepthsPhi[maxEta_];
0429   int indxEta[maxEta_][maxDep_][maxPhi_];
0430   TFile *output_file;
0431   std::map<unsigned int, float> corrFac_[nmax_];
0432   std::vector<int> runlow_;
0433 
0434   TTree *outtree_;
0435   int t_ieta, t_iphi, t_nvtx;
0436   double t_p, t_ediff;
0437   std::vector<double> t_ene, t_enec, t_actln, t_charge;
0438   std::vector<int> t_depth;
0439 
0440   TH1D *h_evtype, *h_Pt_Muon, *h_Eta_Muon, *h_Phi_Muon, *h_P_Muon;
0441   TH1D *h_PF_Muon, *h_GlobTrack_Chi, *h_Global_Muon_Hits;
0442   TH1D *h_MatchedStations, *h_Tight_TransImpactparameter;
0443   TH1D *h_Tight_LongitudinalImpactparameter, *h_InnerTrackPixelHits;
0444   TH1D *h_TrackerLayer, *h_IsolationR04, *h_Global_Muon;
0445   TH1D *h_LongImpactParameter, *h_LongImpactParameterBin1, *h_LongImpactParameterBin2;
0446 
0447   TH1D *h_TransImpactParameter, *h_TransImpactParameterBin1, *h_TransImpactParameterBin2;
0448   TH1D *h_Hot_MuonEnergy_hcal_ClosestCell[maxHist_], *h_Hot_MuonEnergy_hcal_HotCell[maxHist_],
0449       *h_Hot_MuonEnergy_hcal_HotCell_VsActiveLength[maxHist_], *h_HotCell_MuonEnergy_phi[maxHist_],
0450       *h_active_length_Fill[maxHist_], *h_p_muon_ineta[maxHist_], *h_charge_signal[maxHist_], *h_charge_bg[maxHist_];
0451   TH2D *h_2D_Bin1, *h_2D_Bin2;
0452   TH1D *h_ecal_energy, *h_hcal_energy, *h_3x3_ecal, *h_1x1_hcal;
0453   TH1D *h_MuonHittingEcal, *h_HotCell, *h_MuonEnergy_hcal[maxHist_];
0454   TH1D *h_Hot_MuonEnergy_hcal[maxHist_];
0455   TH2D *hcal_ietaVsEnergy;
0456   TProfile *h_EtaX_hcal, *h_PhiY_hcal, *h_EtaX_ecal, *h_PhiY_ecal;
0457   TProfile *h_Eta_ecal, *h_Phi_ecal;
0458   TProfile *h_MuonEnergy_eta[maxDep_], *h_MuonEnergy_phi[maxDep_], *h_MuonEnergy_muon_eta[maxDep_];
0459   TProfile *h_Hot_MuonEnergy_eta[maxDep_], *h_Hot_MuonEnergy_phi[maxDep_], *h_Hot_MuonEnergy_muon_eta[maxDep_];
0460   TProfile *h_IsoHot_MuonEnergy_eta[maxDep_], *h_IsoHot_MuonEnergy_phi[maxDep_], *h_IsoHot_MuonEnergy_muon_eta[maxDep_];
0461   TProfile *h_IsoWithoutHot_MuonEnergy_eta[maxDep_], *h_IsoWithoutHot_MuonEnergy_phi[maxDep_],
0462       *h_IsoWithoutHot_MuonEnergy_muon_eta[maxDep_];
0463   TProfile *h_HotWithoutIso_MuonEnergy_eta[maxDep_], *h_HotWithoutIso_MuonEnergy_phi[maxDep_],
0464       *h_HotWithoutIso_MuonEnergy_muon_eta[maxDep_];
0465 };
0466 
0467 HBHEMuonOfflineAnalyzer::HBHEMuonOfflineAnalyzer(TChain *tree,
0468                                                  const char *outFileName,
0469                                                  const char *rcorFileName,
0470                                                  int flag,
0471                                                  int mode,
0472                                                  int maxDHB,
0473                                                  int maxDHE,
0474                                                  int cutMu,
0475                                                  float cutP,
0476                                                  int nevMax,
0477                                                  int over,
0478                                                  int runLo,
0479                                                  int runHi,
0480                                                  int etaMin,
0481                                                  int etaMax,
0482                                                  bool deb)
0483     : nCut_(cutMu), cutP_(cutP), nevMax_(nevMax), over_(over == 1), debug_(deb), cFactor_(false) {
0484   if ((nCut_ < 0) || (nCut_ > 2))
0485     nCut_ = 0;
0486   Init(tree, rcorFileName, flag, mode, maxDHB, maxDHE, runLo, runHi, etaMin, etaMax);
0487 
0488   //Now book histograms
0489   bookHistograms(outFileName);
0490 }
0491 
0492 HBHEMuonOfflineAnalyzer::HBHEMuonOfflineAnalyzer(const char *infile,
0493                                                  const char *outFileName,
0494                                                  const char *rcorFileName,
0495                                                  int flag,
0496                                                  int mode,
0497                                                  int maxDHB,
0498                                                  int maxDHE,
0499                                                  int cutMu,
0500                                                  float cutP,
0501                                                  int nevMax,
0502                                                  int over,
0503                                                  int runLo,
0504                                                  int runHi,
0505                                                  int etaMin,
0506                                                  int etaMax,
0507                                                  bool deb)
0508     : nCut_(cutMu), cutP_(cutP), nevMax_(nevMax), over_(over == 1), debug_(deb), cFactor_(false) {
0509   if ((nCut_ < 0) || (nCut_ > 2))
0510     nCut_ = 0;
0511   TChain *chain = new TChain("hcalHBHEMuon/TREE");
0512   if (!fillChain(chain, infile)) {
0513     std::cout << "*****No valid tree chain can be obtained*****" << std::endl;
0514   } else {
0515     std::cout << "Proceed with a tree chain with " << chain->GetEntries() << " entries" << std::endl;
0516     Init(chain, rcorFileName, flag, mode, maxDHB, maxDHE, runLo, runHi, etaMin, etaMax);
0517 
0518     //Now book histograms
0519     bookHistograms(outFileName);
0520   }
0521 }
0522 
0523 HBHEMuonOfflineAnalyzer::~HBHEMuonOfflineAnalyzer() {
0524   if (!fChain)
0525     return;
0526   delete fChain->GetCurrentFile();
0527 }
0528 
0529 Int_t HBHEMuonOfflineAnalyzer::Cut(Long64_t) {
0530   // This function may be called from Loop.
0531   // returns  1 if entry is accepted.
0532   // returns -1 otherwise.
0533   return 1;
0534 }
0535 
0536 Int_t HBHEMuonOfflineAnalyzer::GetEntry(Long64_t entry) {
0537   // Read contents of entry.
0538   if (!fChain)
0539     return 0;
0540   return fChain->GetEntry(entry);
0541 }
0542 
0543 Long64_t HBHEMuonOfflineAnalyzer::LoadTree(Long64_t entry) {
0544   // Set the environment to read one entry
0545   if (!fChain)
0546     return -5;
0547   Long64_t centry = fChain->LoadTree(entry);
0548   if (centry < 0)
0549     return centry;
0550   if (fChain->GetTreeNumber() != fCurrent) {
0551     fCurrent = fChain->GetTreeNumber();
0552     Notify();
0553   }
0554   return centry;
0555 }
0556 
0557 void HBHEMuonOfflineAnalyzer::Init(TChain *tree,
0558                                    const char *rcorFileName,
0559                                    int flag,
0560                                    int mode,
0561                                    int maxDHB,
0562                                    int maxDHE,
0563                                    int runLo,
0564                                    int runHi,
0565                                    int etaMin,
0566                                    int etaMax) {
0567   modeLHC_ = mode;
0568   maxDepthHB_ = maxDHB;
0569   maxDepthHE_ = maxDHE;
0570   maxDepth_ = (maxDepthHB_ > maxDepthHE_) ? maxDepthHB_ : maxDepthHE_;
0571   runLo_ = runLo;
0572   runHi_ = runHi;
0573   etaMin_ = (etaMin > 0) ? etaMin : 1;
0574   etaMax_ = (etaMax <= 29) ? etaMax : 29;
0575   if (etaMax_ <= etaMin_) {
0576     if (etaMax_ == 29)
0577       etaMin_ = etaMax_ - 1;
0578     else
0579       etaMax_ = etaMin_ + 1;
0580   }
0581   useCorrect_ = ((flag % 10) > 0);
0582   mergeDepth_ = (((flag / 10) % 10) > 0);
0583   if (std::string(rcorFileName) != "")
0584     cFactor_ = readCorr(rcorFileName);
0585 
0586   // The Init() function is called when the selector needs to initialize
0587   // a new tree or chain. Typically here the branch addresses and branch
0588   // pointers of the tree will be set.
0589   // It is normally not necessary to make changes to the generated
0590   // code, but the routine can be extended by the user if needed.
0591   // Init() will be called many times when running on PROOF
0592   // (once per file to be processed).
0593 
0594   // Set object pointer
0595 
0596   PF_Muon = 0;
0597   Global_Muon = 0;
0598   Tracker_muon = 0;
0599   pt_of_muon = 0;
0600   eta_of_muon = 0;
0601   phi_of_muon = 0;
0602   energy_of_muon = 0;
0603   p_of_muon = 0;
0604   muon_trkKink = 0;
0605   muon_chi2LocalPosition = 0;
0606   muon_segComp = 0;
0607   TrackerLayer = 0;
0608   NumPixelLayers = 0;
0609   InnerTrackPixelHits = 0;
0610   innerTrack = 0;
0611   chiTracker = 0;
0612   DxyTracker = 0;
0613   DzTracker = 0;
0614   innerTrackpt = 0;
0615   innerTracketa = 0;
0616   innerTrackphi = 0;
0617   tight_validFraction = 0;
0618   OuterTrack = 0;
0619   OuterTrackPt = 0;
0620   OuterTrackEta = 0;
0621   OuterTrackPhi = 0;
0622   OuterTrackHits = 0;
0623   OuterTrackRHits = 0;
0624   OuterTrackChi = 0;
0625   GlobalTrack = 0;
0626   GlobalTrckPt = 0;
0627   GlobalTrckEta = 0;
0628   GlobalTrckPhi = 0;
0629   Global_Muon_Hits = 0;
0630   MatchedStations = 0;
0631   GlobTrack_Chi = 0;
0632   Tight_LongitudinalImpactparameter = 0;
0633   Tight_TransImpactparameter = 0;
0634   IsolationR04 = 0;
0635   IsolationR03 = 0;
0636   ecal_3into3 = 0;
0637   hcal_3into3 = 0;
0638   tracker_3into3 = 0;
0639   matchedId = 0;
0640   hcal_cellHot = 0;
0641   ecal_3x3 = 0;
0642   hcal_1x1 = 0;
0643   ecal_detID = 0;
0644   hcal_detID = 0;
0645   ehcal_detID = 0;
0646   hcal_edepth1 = 0;
0647   hcal_activeL1 = 0;
0648   hcal_edepthHot1 = 0;
0649   hcal_activeHotL1 = 0;
0650   hcal_edepthCorrect1 = 0;
0651   hcal_edepthHotCorrect1 = 0;
0652   hcal_cdepthHot1 = 0;
0653   hcal_cdepthHotBG1 = 0;
0654   hcal_depthMatch1 = 0;
0655   hcal_depthMatchHot1 = 0;
0656   hcal_edepth2 = 0;
0657   hcal_activeL2 = 0;
0658   hcal_edepthHot2 = 0;
0659   hcal_activeHotL2 = 0;
0660   hcal_edepthCorrect2 = 0;
0661   hcal_edepthHotCorrect2 = 0;
0662   hcal_cdepthHot2 = 0;
0663   hcal_cdepthHotBG2 = 0;
0664   hcal_depthMatch2 = 0;
0665   hcal_depthMatchHot2 = 0;
0666   hcal_edepth3 = 0;
0667   hcal_activeL3 = 0;
0668   hcal_edepthHot3 = 0;
0669   hcal_activeHotL3 = 0;
0670   hcal_edepthCorrect3 = 0;
0671   hcal_edepthHotCorrect3 = 0;
0672   hcal_cdepthHot3 = 0;
0673   hcal_cdepthHotBG3 = 0;
0674   hcal_depthMatch3 = 0;
0675   hcal_depthMatchHot3 = 0;
0676   hcal_edepth4 = 0;
0677   hcal_activeL4 = 0;
0678   hcal_edepthHot4 = 0;
0679   hcal_activeHotL4 = 0;
0680   hcal_edepthCorrect4 = 0;
0681   hcal_edepthHotCorrect4 = 0;
0682   hcal_cdepthHot4 = 0;
0683   hcal_cdepthHotBG4 = 0;
0684   hcal_depthMatch4 = 0;
0685   hcal_depthMatchHot4 = 0;
0686   hcal_edepth5 = 0;
0687   hcal_activeL5 = 0;
0688   hcal_edepthHot5 = 0;
0689   hcal_activeHotL5 = 0;
0690   hcal_edepthCorrect5 = 0;
0691   hcal_edepthHotCorrect5 = 0;
0692   hcal_cdepthHot5 = 0;
0693   hcal_cdepthHotBG5 = 0;
0694   hcal_depthMatch5 = 0;
0695   hcal_depthMatchHot5 = 0;
0696   hcal_edepth6 = 0;
0697   hcal_activeL6 = 0;
0698   hcal_edepthHot6 = 0;
0699   hcal_activeHotL6 = 0;
0700   hcal_edepthCorrect6 = 0;
0701   hcal_edepthHotCorrect6 = 0;
0702   hcal_cdepthHot6 = 0;
0703   hcal_cdepthHotBG6 = 0;
0704   hcal_depthMatch6 = 0;
0705   hcal_depthMatchHot6 = 0;
0706   hcal_edepth7 = 0;
0707   hcal_activeL7 = 0;
0708   hcal_edepthHot7 = 0;
0709   hcal_activeHotL7 = 0;
0710   hcal_edepthCorrect7 = 0;
0711   hcal_edepthHotCorrect7 = 0;
0712   hcal_cdepthHot7 = 0;
0713   hcal_cdepthHotBG7 = 0;
0714   hcal_depthMatch7 = 0;
0715   hcal_depthMatchHot7 = 0;
0716   activeLength = 0;
0717   activeLengthHot = 0;
0718   hltresults = 0;
0719   all_triggers = 0;
0720   // Set branch addresses and branch pointers
0721   if (!tree)
0722     return;
0723   fChain = tree;
0724   fCurrent = -1;
0725   fChain->SetMakeClass(1);
0726 
0727   fChain->SetBranchAddress("Event_No", &Event_No, &b_Event_No);
0728   fChain->SetBranchAddress("Run_No", &Run_No, &b_Run_No);
0729   fChain->SetBranchAddress("LumiNumber", &LumiNumber, &b_LumiNumber);
0730   fChain->SetBranchAddress("BXNumber", &BXNumber, &b_BXNumber);
0731   fChain->SetBranchAddress("GoodVertex", &GoodVertex, &b_GoodVertex);
0732   fChain->SetBranchAddress("PF_Muon", &PF_Muon, &b_PF_Muon);
0733   fChain->SetBranchAddress("Global_Muon", &Global_Muon, &b_Global_Muon);
0734   fChain->SetBranchAddress("Tracker_muon", &Tracker_muon, &b_Tracker_muon);
0735   fChain->SetBranchAddress("MuonIsTight", &MuonIsTight, &b_MuonIsTight);
0736   fChain->SetBranchAddress("MuonIsMedium", &MuonIsMedium, &b_MuonIsMedium);
0737   fChain->SetBranchAddress("pt_of_muon", &pt_of_muon, &b_pt_of_muon);
0738   fChain->SetBranchAddress("eta_of_muon", &eta_of_muon, &b_eta_of_muon);
0739   fChain->SetBranchAddress("phi_of_muon", &phi_of_muon, &b_phi_of_muon);
0740   fChain->SetBranchAddress("energy_of_muon", &energy_of_muon, &b_energy_of_muon);
0741   fChain->SetBranchAddress("p_of_muon", &p_of_muon, &b_p_of_muon);
0742   fChain->SetBranchAddress("muon_trkKink", &muon_trkKink, &b_muon_trkKink);
0743   fChain->SetBranchAddress("muon_chi2LocalPosition", &muon_chi2LocalPosition, &b_muon_chi2LocalPosition);
0744   fChain->SetBranchAddress("muon_segComp", &muon_segComp, &b_muon_segComp);
0745   fChain->SetBranchAddress("TrackerLayer", &TrackerLayer, &b_TrackerLayer);
0746   fChain->SetBranchAddress("NumPixelLayers", &NumPixelLayers, &b_NumPixelLayers);
0747   fChain->SetBranchAddress("InnerTrackPixelHits", &InnerTrackPixelHits, &b_InnerTrackPixelHits);
0748   fChain->SetBranchAddress("innerTrack", &innerTrack, &b_innerTrack);
0749   fChain->SetBranchAddress("chiTracker", &chiTracker, &b_chiTracker);
0750   fChain->SetBranchAddress("DxyTracker", &DxyTracker, &b_DxyTracker);
0751   fChain->SetBranchAddress("DzTracker", &DzTracker, &b_DzTracker);
0752   fChain->SetBranchAddress("innerTrackpt", &innerTrackpt, &b_innerTrackpt);
0753   fChain->SetBranchAddress("innerTracketa", &innerTracketa, &b_innerTracketa);
0754   fChain->SetBranchAddress("innerTrackphi", &innerTrackphi, &b_innerTrackphi);
0755   fChain->SetBranchAddress("tight_validFraction", &tight_validFraction, &b_tight_validFraction);
0756   fChain->SetBranchAddress("OuterTrack", &OuterTrack, &b_OuterTrack);
0757   fChain->SetBranchAddress("OuterTrackPt", &OuterTrackPt, &b_OuterTrackPt);
0758   fChain->SetBranchAddress("OuterTrackEta", &OuterTrackEta, &b_OuterTrackEta);
0759   fChain->SetBranchAddress("OuterTrackPhi", &OuterTrackPhi, &b_OuterTrackPhi);
0760   fChain->SetBranchAddress("OuterTrackChi", &OuterTrackChi, &b_OuterTrackChi);
0761   fChain->SetBranchAddress("OuterTrackHits", &OuterTrackHits, &b_OuterTrackHits);
0762   fChain->SetBranchAddress("OuterTrackRHits", &OuterTrackRHits, &b_OuterTrackRHits);
0763   fChain->SetBranchAddress("GlobalTrack", &GlobalTrack, &b_GlobalTrack);
0764   fChain->SetBranchAddress("GlobalTrckPt", &GlobalTrckPt, &b_GlobalTrckPt);
0765   fChain->SetBranchAddress("GlobalTrckEta", &GlobalTrckEta, &b_GlobalTrckEta);
0766   fChain->SetBranchAddress("GlobalTrckPhi", &GlobalTrckPhi, &b_GlobalTrckPhi);
0767   fChain->SetBranchAddress("Global_Muon_Hits", &Global_Muon_Hits, &b_Global_Muon_Hits);
0768   fChain->SetBranchAddress("MatchedStations", &MatchedStations, &b_MatchedStations);
0769   fChain->SetBranchAddress("GlobTrack_Chi", &GlobTrack_Chi, &b_GlobTrack_Chi);
0770   fChain->SetBranchAddress(
0771       "Tight_LongitudinalImpactparameter", &Tight_LongitudinalImpactparameter, &b_Tight_LongitudinalImpactparameter);
0772   fChain->SetBranchAddress("Tight_TransImpactparameter", &Tight_TransImpactparameter, &b_Tight_TransImpactparameter);
0773   fChain->SetBranchAddress("IsolationR04", &IsolationR04, &b_IsolationR04);
0774   fChain->SetBranchAddress("IsolationR03", &IsolationR03, &b_IsolationR03);
0775   fChain->SetBranchAddress("ecal_3into3", &ecal_3into3, &b_ecal_3into3);
0776   fChain->SetBranchAddress("hcal_3into3", &hcal_3into3, &b_hcal_3into3);
0777   fChain->SetBranchAddress("tracker_3into3", &tracker_3into3, &b_tracker_3into3);
0778   fChain->SetBranchAddress("matchedId", &matchedId, &b_matchedId);
0779   fChain->SetBranchAddress("hcal_cellHot", &hcal_cellHot, &b_hcal_cellHot);
0780   fChain->SetBranchAddress("ecal_3x3", &ecal_3x3, &b_ecal_3x3);
0781   fChain->SetBranchAddress("hcal_1x1", &hcal_1x1, &b_hcal_1x1);
0782   fChain->SetBranchAddress("ecal_detID", &ecal_detID, &b_ecal_detID);
0783   fChain->SetBranchAddress("hcal_detID", &hcal_detID, &b_hcal_detID);
0784   fChain->SetBranchAddress("ehcal_detID", &ehcal_detID, &b_ehcal_detID);
0785   fChain->SetBranchAddress("hcal_ieta", &hcal_ieta, &b_hcal_ieta);
0786   fChain->SetBranchAddress("hcal_iphi", &hcal_iphi, &b_hcal_iphi);
0787   fChain->SetBranchAddress("hcal_edepth1", &hcal_edepth1, &b_hcal_edepth1);
0788   fChain->SetBranchAddress("hcal_activeL1", &hcal_activeL1, &b_hcal_activeL1);
0789   fChain->SetBranchAddress("hcal_edepthHot1", &hcal_edepthHot1, &b_hcal_edepthHot1);
0790   fChain->SetBranchAddress("hcal_activeHotL1", &hcal_activeHotL1, &b_hcal_activeHotL1);
0791   fChain->SetBranchAddress("hcal_edepthCorrect1", &hcal_edepthCorrect1, &b_hcal_edepthCorrect1);
0792   fChain->SetBranchAddress("hcal_edepthHotCorrect1", &hcal_edepthHotCorrect1, &b_hcal_edepthHotCorrect1);
0793   fChain->SetBranchAddress("hcal_cdepthHot1", &hcal_cdepthHot1, &b_hcal_cdepthHot1);
0794   fChain->SetBranchAddress("hcal_cdepthHotBG1", &hcal_cdepthHotBG1, &b_hcal_cdepthHotBG1);
0795   fChain->SetBranchAddress("hcal_depthMatch1", &hcal_depthMatch1, &b_hcal_depthMatch1);
0796   fChain->SetBranchAddress("hcal_depthMatchHot1", &hcal_depthMatchHot1, &b_hcal_depthMatchHot1);
0797   fChain->SetBranchAddress("hcal_edepth2", &hcal_edepth2, &b_hcal_edepth2);
0798   fChain->SetBranchAddress("hcal_activeL2", &hcal_activeL2, &b_hcal_activeL2);
0799   fChain->SetBranchAddress("hcal_edepthHot2", &hcal_edepthHot2, &b_hcal_edepthHot2);
0800   fChain->SetBranchAddress("hcal_activeHotL2", &hcal_activeHotL2, &b_hcal_activeHotL2);
0801   fChain->SetBranchAddress("hcal_edepthCorrect2", &hcal_edepthCorrect2, &b_hcal_edepthCorrect2);
0802   fChain->SetBranchAddress("hcal_edepthHotCorrect2", &hcal_edepthHotCorrect2, &b_hcal_edepthHotCorrect2);
0803   fChain->SetBranchAddress("hcal_cdepthHot2", &hcal_cdepthHot2, &b_hcal_cdepthHot2);
0804   fChain->SetBranchAddress("hcal_cdepthHotBG2", &hcal_cdepthHotBG2, &b_hcal_cdepthHotBG2);
0805   fChain->SetBranchAddress("hcal_depthMatch2", &hcal_depthMatch2, &b_hcal_depthMatch2);
0806   fChain->SetBranchAddress("hcal_depthMatchHot2", &hcal_depthMatchHot2, &b_hcal_depthMatchHot2);
0807   fChain->SetBranchAddress("hcal_edepth3", &hcal_edepth3, &b_hcal_edepth3);
0808   fChain->SetBranchAddress("hcal_activeL3", &hcal_activeL3, &b_hcal_activeL3);
0809   fChain->SetBranchAddress("hcal_edepthHot3", &hcal_edepthHot3, &b_hcal_edepthHot3);
0810   fChain->SetBranchAddress("hcal_activeHotL3", &hcal_activeHotL3, &b_hcal_activeHotL3);
0811   fChain->SetBranchAddress("hcal_edepthCorrect3", &hcal_edepthCorrect3, &b_hcal_edepthCorrect3);
0812   fChain->SetBranchAddress("hcal_edepthHotCorrect3", &hcal_edepthHotCorrect3, &b_hcal_edepthHotCorrect3);
0813   fChain->SetBranchAddress("hcal_cdepthHot3", &hcal_cdepthHot3, &b_hcal_cdepthHot3);
0814   fChain->SetBranchAddress("hcal_cdepthHotBG3", &hcal_cdepthHotBG3, &b_hcal_cdepthHotBG3);
0815   fChain->SetBranchAddress("hcal_depthMatch3", &hcal_depthMatch3, &b_hcal_depthMatch3);
0816   fChain->SetBranchAddress("hcal_depthMatchHot3", &hcal_depthMatchHot3, &b_hcal_depthMatchHot3);
0817   fChain->SetBranchAddress("hcal_edepth4", &hcal_edepth4, &b_hcal_edepth4);
0818   fChain->SetBranchAddress("hcal_activeL4", &hcal_activeL4, &b_hcal_activeL4);
0819   fChain->SetBranchAddress("hcal_edepthHot4", &hcal_edepthHot4, &b_hcal_edepthHot4);
0820   fChain->SetBranchAddress("hcal_activeHotL4", &hcal_activeHotL4, &b_hcal_activeHotL4);
0821   fChain->SetBranchAddress("hcal_edepthCorrect4", &hcal_edepthCorrect4, &b_hcal_edepthCorrect4);
0822   fChain->SetBranchAddress("hcal_edepthHotCorrect4", &hcal_edepthHotCorrect4, &b_hcal_edepthHotCorrect4);
0823   fChain->SetBranchAddress("hcal_cdepthHot4", &hcal_cdepthHot4, &b_hcal_cdepthHot4);
0824   fChain->SetBranchAddress("hcal_cdepthHotBG4", &hcal_cdepthHotBG4, &b_hcal_cdepthHotBG4);
0825   fChain->SetBranchAddress("hcal_depthMatch4", &hcal_depthMatch4, &b_hcal_depthMatch4);
0826   fChain->SetBranchAddress("hcal_depthMatchHot4", &hcal_depthMatchHot4, &b_hcal_depthMatchHot4);
0827   fChain->SetBranchAddress("hcal_edepth5", &hcal_edepth5, &b_hcal_edepth5);
0828   fChain->SetBranchAddress("hcal_activeL5", &hcal_activeL5, &b_hcal_activeL5);
0829   fChain->SetBranchAddress("hcal_edepthHot5", &hcal_edepthHot5, &b_hcal_edepthHot5);
0830   fChain->SetBranchAddress("hcal_activeHotL5", &hcal_activeHotL5, &b_hcal_activeHotL5);
0831   fChain->SetBranchAddress("hcal_edepthCorrect5", &hcal_edepthCorrect5, &b_hcal_edepthCorrect5);
0832   fChain->SetBranchAddress("hcal_edepthHotCorrect5", &hcal_edepthHotCorrect5, &b_hcal_edepthHotCorrect5);
0833   fChain->SetBranchAddress("hcal_cdepthHot5", &hcal_cdepthHot5, &b_hcal_cdepthHot5);
0834   fChain->SetBranchAddress("hcal_cdepthHotBG5", &hcal_cdepthHotBG5, &b_hcal_cdepthHotBG5);
0835   fChain->SetBranchAddress("hcal_depthMatch5", &hcal_depthMatch5, &b_hcal_depthMatch5);
0836   fChain->SetBranchAddress("hcal_depthMatchHot5", &hcal_depthMatchHot5, &b_hcal_depthMatchHot5);
0837   fChain->SetBranchAddress("hcal_edepth6", &hcal_edepth6, &b_hcal_edepth6);
0838   fChain->SetBranchAddress("hcal_activeL6", &hcal_activeL6, &b_hcal_activeL6);
0839   fChain->SetBranchAddress("hcal_edepthHot6", &hcal_edepthHot6, &b_hcal_edepthHot6);
0840   fChain->SetBranchAddress("hcal_activeHotL6", &hcal_activeHotL6, &b_hcal_activeHotL6);
0841   fChain->SetBranchAddress("hcal_edepthCorrect6", &hcal_edepthCorrect6, &b_hcal_edepthCorrect6);
0842   fChain->SetBranchAddress("hcal_edepthHotCorrect6", &hcal_edepthHotCorrect6, &b_hcal_edepthHotCorrect6);
0843   fChain->SetBranchAddress("hcal_cdepthHot6", &hcal_cdepthHot6, &b_hcal_cdepthHot6);
0844   fChain->SetBranchAddress("hcal_cdepthHotBG6", &hcal_cdepthHotBG6, &b_hcal_cdepthHotBG6);
0845   fChain->SetBranchAddress("hcal_depthMatch6", &hcal_depthMatch6, &b_hcal_depthMatch6);
0846   fChain->SetBranchAddress("hcal_depthMatchHot6", &hcal_depthMatchHot6, &b_hcal_depthMatchHot6);
0847   fChain->SetBranchAddress("hcal_edepth7", &hcal_edepth7, &b_hcal_edepth7);
0848   fChain->SetBranchAddress("hcal_activeL7", &hcal_activeL7, &b_hcal_activeL7);
0849   fChain->SetBranchAddress("hcal_edepthHot7", &hcal_edepthHot7, &b_hcal_edepthHot7);
0850   fChain->SetBranchAddress("hcal_activeHotL7", &hcal_activeHotL7, &b_hcal_activeHotL7);
0851   fChain->SetBranchAddress("hcal_edepthCorrect7", &hcal_edepthCorrect7, &b_hcal_edepthCorrect7);
0852   fChain->SetBranchAddress("hcal_edepthHotCorrect7", &hcal_edepthHotCorrect7, &b_hcal_edepthHotCorrect7);
0853   fChain->SetBranchAddress("hcal_cdepthHot7", &hcal_cdepthHot7, &b_hcal_cdepthHot7);
0854   fChain->SetBranchAddress("hcal_cdepthHotBG7", &hcal_cdepthHotBG7, &b_hcal_cdepthHotBG7);
0855   fChain->SetBranchAddress("hcal_depthMatch7", &hcal_depthMatch7, &b_hcal_depthMatch7);
0856   fChain->SetBranchAddress("hcal_depthMatchHot7", &hcal_depthMatchHot7, &b_hcal_depthMatchHot7);
0857   fChain->SetBranchAddress("activeLength", &activeLength, &b_activeLength);
0858   fChain->SetBranchAddress("activeLengthHot", &activeLengthHot, &b_activeLengthHot);
0859   fChain->SetBranchAddress("hltresults", &hltresults, &b_hltresults);
0860   fChain->SetBranchAddress("all_triggers", &all_triggers, &b_all_triggers);
0861 
0862   Notify();
0863 }
0864 
0865 void HBHEMuonOfflineAnalyzer::Loop() {
0866   //declarations
0867   if (fChain == 0)
0868     return;
0869 
0870   Long64_t nentries = fChain->GetEntriesFast();
0871   if (debug_)
0872     std::cout << "nevent = " << nentries << std::endl;
0873   if (nevMax_ > 0)
0874     nentries = nevMax_;
0875 
0876   Long64_t nbytes = 0, nb = 0, nsel1 = 0, nsel2 = 0;
0877   Long64_t nstep1 = 0, nstep2 = 0, nstep3 = 0, nstep4 = 0;
0878 
0879   for (Long64_t jentry = 0; jentry < nentries; jentry++) {
0880     Long64_t ientry = LoadTree(jentry);
0881     if (ientry < 0)
0882       break;
0883     nb = fChain->GetEntry(jentry);
0884     nbytes += nb;
0885     if ((int)(Run_No) < runLo_ || (int)(Run_No) > runHi_)
0886       continue;
0887     ++nstep1;
0888     if (debug_)
0889       std::cout << "Run " << Run_No << " Event " << Event_No << " Muon pt " << pt_of_muon << std::endl;
0890 
0891     bool loose(false), soft(false), tight(false), pcut(false), ptcut(false);
0892     t_ene.clear();
0893     t_enec.clear();
0894     t_charge.clear();
0895     t_actln.clear();
0896     t_depth.clear();
0897 
0898     if (debug_)
0899       std::cout << "ecal_det_id " << ecal_detID << std::endl;
0900 
0901     int typeEcal, etaXEcal, phiYEcal, zsideEcal, planeEcal, stripEcal;
0902     etaPhiEcal(ecal_detID, typeEcal, zsideEcal, etaXEcal, phiYEcal, planeEcal, stripEcal);
0903     double etaEcal = (etaXEcal - 0.5) * zsideEcal;
0904     double phiEcal = phiYEcal - 0.5;
0905 
0906     if (debug_)
0907       std::cout << "hcal_det_id " << std::hex << hcal_detID << std::dec;
0908 
0909     int etaHcal, phiHcal, depthHcal;
0910     etaPhiHcal(hcal_detID, etaHcal, phiHcal, depthHcal);
0911 
0912     int eta = (etaHcal > 0) ? etaHcal - 1 : -etaHcal - 1;
0913     double etaXHcal = (etaHcal > 0) ? etaHcal - 0.5 : etaHcal + 0.5;
0914     int nDepth = nDepthBins(eta + 1, phiHcal);
0915     int nPhi = nPhiBins(eta + 1);
0916     int PHI = (nPhi > 36) ? (phiHcal - 1) : (phiHcal - 1) / 2;
0917     double phiYHcal = (phiHcal - 0.5);
0918     t_ieta = etaHcal;
0919     t_iphi = PHI;
0920     t_p = p_of_muon;
0921     t_ediff = hcal_3into3 - hcal_1x1;
0922     t_nvtx = GoodVertex;
0923     if (p_of_muon > cutP_)
0924       pcut = true;
0925     if (pt_of_muon > cutP_)
0926       ptcut = true;
0927     if (looseMuon())
0928       loose = true;
0929     if (softMuon())
0930       soft = true;
0931     if (tightMuon())
0932       tight = true;
0933 
0934     if (debug_)
0935       std::cout << " etaHcal " << etaHcal << ":" << etaXHcal << " phiHcal " << phiHcal << ":" << phiYHcal << ":" << PHI
0936                 << " Depth " << nDepth << " Muon Pt " << pt_of_muon << " Isol " << IsolationR04 << std::endl;
0937 
0938     int cut(nCut_);
0939     bool select(false);
0940     if (cut == 0)
0941       select = tightMuon();
0942     else if (cut == 1)
0943       select = softMuon();
0944     else
0945       select = looseMuon();
0946 
0947     if (select)
0948       ++nstep2;
0949     if (select && ((eta + 1) >= etaMin_) && ((eta + 1) <= etaMax_)) {
0950       ++nstep3;
0951       h_Pt_Muon->Fill(pt_of_muon);
0952       h_Eta_Muon->Fill(eta_of_muon);
0953       h_Phi_Muon->Fill(phi_of_muon);
0954       h_PF_Muon->Fill(PF_Muon);
0955       h_GlobTrack_Chi->Fill(GlobTrack_Chi);
0956       h_Global_Muon_Hits->Fill(Global_Muon_Hits);
0957       h_MatchedStations->Fill(MatchedStations);
0958       h_Tight_TransImpactparameter->Fill(Tight_TransImpactparameter);
0959       h_Tight_LongitudinalImpactparameter->Fill(Tight_LongitudinalImpactparameter);
0960       h_InnerTrackPixelHits->Fill(InnerTrackPixelHits);
0961       h_TrackerLayer->Fill(TrackerLayer);
0962       h_IsolationR04->Fill(IsolationR04);
0963       h_Global_Muon->Fill(Global_Muon);
0964 
0965       h_TransImpactParameter->Fill(Tight_TransImpactparameter);
0966       h_LongImpactParameter->Fill(Tight_LongitudinalImpactparameter);
0967 
0968       //in Phi Bins
0969       if (((phi_of_muon) >= -1.5) || ((phi_of_muon) <= 0.5)) {
0970         h_TransImpactParameterBin1->Fill(Tight_TransImpactparameter);
0971         h_LongImpactParameterBin1->Fill(Tight_LongitudinalImpactparameter);
0972         h_2D_Bin1->Fill(Tight_TransImpactparameter, Tight_LongitudinalImpactparameter);
0973       }
0974 
0975       if ((phi_of_muon > 0.5) || (phi_of_muon < -1.5)) {
0976         h_TransImpactParameterBin2->Fill(Tight_TransImpactparameter);
0977         h_LongImpactParameterBin2->Fill(Tight_LongitudinalImpactparameter);
0978         h_2D_Bin2->Fill(Tight_TransImpactparameter, Tight_LongitudinalImpactparameter);
0979       }
0980 
0981       h_ecal_energy->Fill(ecal_3into3);
0982       h_3x3_ecal->Fill(ecal_3x3);
0983       h_Eta_ecal->Fill(eta_of_muon, ecal_3x3);
0984       h_Phi_ecal->Fill(phi_of_muon, ecal_3x3);
0985       h_MuonHittingEcal->Fill(typeEcal);
0986       if (typeEcal == 1) {
0987         h_EtaX_ecal->Fill(etaEcal, ecal_3x3);
0988         h_PhiY_ecal->Fill(phiEcal, ecal_3x3);
0989       }
0990 
0991       h_hcal_energy->Fill(hcal_3into3);
0992       h_1x1_hcal->Fill(hcal_1x1);
0993       h_EtaX_hcal->Fill(etaXHcal, hcal_1x1);
0994       h_PhiY_hcal->Fill(phiYHcal, hcal_1x1);
0995       h_HotCell->Fill(hcal_cellHot);
0996       if (mergeDepth_) {
0997         double en1(0), en2(0), actLTot(0), chargeS(0), chargeBG(0);
0998         double enh(0), enc(0);
0999         for (int dep = 0; dep < nDepth; ++dep) {
1000           double enb(0), enu(0), eh0(0), ec0(0), chgS(0), chgB(0), actL(0);
1001           getEnergy(dep, enb, enu, eh0, ec0, chgS, chgB, actL);
1002           en1 += ((useCorrect_) ? enu : enb);
1003           en2 += ((useCorrect_) ? ec0 : eh0);
1004           enh += (eh0);
1005           enc += (ec0);
1006           actLTot += (actL);
1007           chargeS += (chgS);
1008           chargeBG += (chgB);
1009         }
1010         int ind = (etaHcal > 0) ? indxEta[eta][0][PHI] : 1 + indxEta[eta][0][PHI];
1011         if (debug_)  // || eta==15 || eta==17)
1012           std::cout << "Matched Id " << matchedId << " Hot " << hcal_cellHot << " eta " << etaHcal << ":" << eta
1013                     << " phi " << phiHcal << ":" << PHI << " Index " << ind << " E " << en1 << ":" << en2 << ":" << enh
1014                     << ":" << enc << " L " << actLTot << " Charge " << chargeS << ":" << chargeBG << std::endl;
1015         if (!(matchedId) && !(over_))
1016           continue;
1017         if ((hcal_cellHot == 1) || over_) {
1018           if (actLTot > 0) {
1019             h_Hot_MuonEnergy_hcal_HotCell[ind]->Fill(en2);
1020             h_Hot_MuonEnergy_hcal_HotCell_VsActiveLength[ind]->Fill(en2 / actLTot);
1021             h_active_length_Fill[ind]->Fill(actLTot);
1022             h_p_muon_ineta[ind]->Fill(p_of_muon);
1023             h_charge_signal[ind]->Fill(chargeS);
1024             h_charge_bg[ind]->Fill(chargeBG);
1025 
1026             t_ene.push_back(enh);
1027             t_enec.push_back(enc);
1028             t_charge.push_back(chargeS);
1029             t_actln.push_back(actLTot);
1030             t_depth.push_back(0);
1031 
1032             outtree_->Fill();
1033             ++nsel1;
1034           }
1035         }
1036       } else {
1037         bool fillTree(false);
1038         ++nstep4;
1039         for (int dep = 0; dep < nDepth; ++dep) {
1040           if (debug_)
1041             std::cout << "dep:" << dep << std::endl;
1042 
1043           double actL(0), chargeS(-9999), chargeBG(-9999);
1044           double enh(-9999), enc(-9999), enb(0), enu(0);
1045           bool ok1 = getEnergy(dep, enb, enu, enh, enc, chargeS, chargeBG, actL);
1046           double en1 = ((useCorrect_) ? enu : enb);
1047           double en2 = ((useCorrect_) ? enc : enh);
1048           if (debug_)
1049             std::cout << "Hello in " << dep + 1 << " " << en1 << ":" << en2 << ":" << actL << std::endl;
1050 
1051           bool ok2 = ok1;
1052 
1053           if (debug_)
1054             std::cout << "Before Index " << ok1 << ":" << ok2 << std::endl;
1055 
1056           int ind = (etaHcal > 0) ? indxEta[eta][dep][PHI] : 1 + indxEta[eta][dep][PHI];
1057           if (debug_)  // || eta==15 || eta==17)
1058             std::cout << "Matched Id " << matchedId << " Hot " << hcal_cellHot << " eta " << etaHcal << ":" << eta
1059                       << " phi " << phiHcal << ":" << PHI << " depth " << dep << " Index " << ind << " E " << en1 << ":"
1060                       << en2 << ":" << enh << ":" << enc << " L " << actL << " Charge " << chargeS << ":" << chargeBG
1061                       << std::endl;
1062           if (debug_)
1063             std::cout << "matchedId " << matchedId << " Over " << over_ << " OK " << ok1 << ":" << ok2 << " cellHot "
1064                       << (hcal_cellHot == 1) << std::endl;
1065           if (!(matchedId) && !(over_))
1066             continue;
1067           if (ok1 || over_) {
1068             if (debug_)
1069               std::cout << "enter ok1" << std::endl;
1070 
1071             if ((hcal_cellHot == 1) || (over_)) {
1072               if (actL > 0) {
1073                 h_Hot_MuonEnergy_hcal_HotCell[ind]->Fill(en2);
1074                 h_Hot_MuonEnergy_hcal_HotCell_VsActiveLength[ind]->Fill(en2 / actL);
1075                 h_active_length_Fill[ind]->Fill(actL);
1076                 h_p_muon_ineta[ind]->Fill(p_of_muon);
1077                 h_charge_signal[ind]->Fill(chargeS);
1078                 h_charge_bg[ind]->Fill(chargeBG);
1079                 t_ene.push_back(enh);
1080                 t_enec.push_back(enc);
1081                 t_charge.push_back(chargeS);
1082                 t_actln.push_back(actL);
1083                 // added depth vector AmanKalsi
1084                 t_depth.push_back(dep);
1085                 fillTree = true;
1086               } else {
1087                 t_ene.push_back(-999.0);
1088                 t_enec.push_back(-999.0);
1089                 t_charge.push_back(-999.0);
1090                 t_actln.push_back(-999.0);
1091                 t_depth.push_back(-999.0);
1092               }
1093               if (debug_)
1094                 std::cout << "enter hot cell" << std::endl;
1095             }
1096           }
1097 
1098           if (ok2) {
1099             if (debug_)
1100               std::cout << "enter ok2" << std::endl;
1101             if (hcal_cellHot != 1) {
1102             }
1103           }
1104 
1105           if (debug_)
1106             std::cout << "ETA \t" << eta << "DEPTH \t" << dep << std::endl;
1107         }
1108         if (fillTree) {
1109           outtree_->Fill();
1110           ++nsel2;
1111         }
1112       }
1113     }
1114     int evtype(0);
1115     if (pcut)
1116       evtype += 1;
1117     if (ptcut)
1118       evtype += 2;
1119     if (loose)
1120       evtype += 4;
1121     if (soft)
1122       evtype += 8;
1123     if (tight)
1124       evtype += 16;
1125     h_evtype->Fill(evtype);
1126   }
1127   std::cout << "Number of events in the output root tree: " << nsel1 << ":" << nsel2 << ":" << nstep1 << ":" << nstep2
1128             << ":" << nstep3 << ":" << nstep4 << std::endl;
1129   close();
1130 }
1131 
1132 Bool_t HBHEMuonOfflineAnalyzer::Notify() {
1133   // The Notify() function is called when a new file is opened. This
1134   // can be either for a new TTree in a TChain or when when a new TTree
1135   // is started when using PROOF. It is normally not necessary to make changes
1136   // to the generated code, but the routine can be extended by the
1137   // user if needed. The return value is currently not used.
1138 
1139   return kTRUE;
1140 }
1141 
1142 void HBHEMuonOfflineAnalyzer::Show(Long64_t entry) {
1143   // Print contents of entry.
1144   // If entry is not specified, print current entry
1145   if (!fChain)
1146     return;
1147   fChain->Show(entry);
1148 }
1149 
1150 bool HBHEMuonOfflineAnalyzer::fillChain(TChain *chain, const char *inputFileList) {
1151   std::string fname(inputFileList);
1152   if (fname.substr(fname.size() - 5, 5) == ".root") {
1153     chain->Add(fname.c_str());
1154   } else {
1155     std::ifstream infile(inputFileList);
1156     if (!infile.is_open()) {
1157       std::cout << "** ERROR: Can't open '" << inputFileList << "' for input" << std::endl;
1158       return false;
1159     }
1160     while (1) {
1161       infile >> fname;
1162       if (!infile.good())
1163         break;
1164       chain->Add(fname.c_str());
1165     }
1166     infile.close();
1167   }
1168   std::cout << "No. of Entries in this tree : " << chain->GetEntries() << std::endl;
1169   return true;
1170 }
1171 
1172 bool HBHEMuonOfflineAnalyzer::readCorr(const char *infile) {
1173   std::ifstream fInput(infile);
1174   unsigned int ncorr(0), all(0), good(0);
1175   if (!fInput.good()) {
1176     std::cout << "Cannot open file " << infile << std::endl;
1177   } else {
1178     char buffer[1024];
1179     while (fInput.getline(buffer, 1024)) {
1180       ++all;
1181       std::string bufferString(buffer);
1182       if (bufferString.substr(0, 5) == "#IOVs") {
1183         std::vector<std::string> items = splitString(bufferString.substr(6));
1184         ncorr = items.size() - 1;
1185         for (unsigned int n = 0; n < ncorr; ++n) {
1186           int run = std::atoi(items[n].c_str());
1187           runlow_.push_back(run);
1188         }
1189         std::cout << ncorr << ":" << runlow_.size() << " Run ranges" << std::endl;
1190         for (unsigned int n = 0; n < runlow_.size(); ++n)
1191           std::cout << " [" << n << "] " << runlow_[n];
1192         std::cout << std::endl;
1193       } else if (buffer[0] == '#') {
1194         continue;  //ignore other comments
1195       } else {
1196         std::vector<std::string> items = splitString(bufferString);
1197         if (items.size() != ncorr + 3) {
1198           std::cout << "Ignore  line: " << buffer << std::endl;
1199         } else {
1200           ++good;
1201           int ieta = std::atoi(items[0].c_str());
1202           int iphi = std::atoi(items[1].c_str());
1203           int depth = std::atoi(items[2].c_str());
1204           unsigned int id = getDetIdHBHE(ieta, iphi, depth);
1205           for (unsigned int n = 0; n < ncorr; ++n) {
1206             float corrf = std::atof(items[n + 3].c_str());
1207             if (n < nmax_)
1208               corrFac_[n][id] = corrf;
1209           }
1210           if (debug_) {
1211             std::cout << "ID " << std::hex << id << std::dec << ":" << id << " (eta " << ieta << " phi " << iphi
1212                       << " depth " << depth << ")";
1213             for (unsigned int n = 0; n < ncorr; ++n)
1214               std::cout << " " << corrFac_[n][id];
1215             std::cout << std::endl;
1216           }
1217         }
1218       }
1219     }
1220     fInput.close();
1221     std::cout << "Reads total of " << all << " and " << good << " good records" << std::endl;
1222   }
1223   return (good > 0);
1224 }
1225 
1226 void HBHEMuonOfflineAnalyzer::bookHistograms(const char *fname) {
1227   std::cout << "BookHistograms" << std::endl;
1228   output_file = TFile::Open(fname, "RECREATE");
1229   output_file->cd();
1230   outtree_ = new TTree("Lep_Tree", "Lep_Tree");
1231   outtree_->Branch("t_ieta", &t_ieta);
1232   outtree_->Branch("t_iphi", &t_iphi);
1233   outtree_->Branch("t_nvtx", &t_nvtx);
1234   outtree_->Branch("t_p", &t_p);
1235   outtree_->Branch("t_ediff", &t_ediff);
1236   outtree_->Branch("t_ene", &t_ene);
1237   outtree_->Branch("t_enec", &t_enec);
1238   outtree_->Branch("t_charge", &t_charge);
1239   outtree_->Branch("t_actln", &t_actln);
1240   outtree_->Branch("t_depth", &t_depth);
1241 
1242   std::string type[] = {"tight", "soft", "loose"};
1243   char name[128], title[500];
1244 
1245   nHist = 0;
1246   for (int eta = etaMin_; eta <= etaMax_; ++eta) {
1247     int nDepth = nDepthBins(eta, -1);
1248     int nPhi = nPhiBins(eta);
1249     for (int depth = 0; depth < nDepth; depth++) {
1250       for (int PHI = 0; PHI < nPhi; ++PHI) {
1251         indxEta[eta - 1][depth][PHI] = nHist;
1252         nHist += 2;
1253       }
1254     }
1255   }
1256   if (nHist >= maxHist_) {
1257     std::cout << "Problem here " << nHist << ":" << maxHist_ << std::endl;
1258   }
1259   //    TDirectory *d_output_file[nCut_][29];
1260   //output_file->cd();
1261   h_evtype = new TH1D("EvType", "Event Type", 100, 0, 100);
1262   int i(nCut_);
1263   sprintf(name, "h_Pt_Muon_%s", type[i].c_str());
1264   sprintf(title, "p_{T} of %s muons (GeV)", type[i].c_str());
1265   h_Pt_Muon = new TH1D(name, title, 100, 0, 200);
1266 
1267   sprintf(name, "h_Eta_Muon_%s", type[i].c_str());
1268   sprintf(title, "#eta of %s muons", type[i].c_str());
1269   h_Eta_Muon = new TH1D(name, title, 50, -2.5, 2.5);
1270 
1271   sprintf(name, "h_Phi_Muon_%s", type[i].c_str());
1272   sprintf(title, "#phi of %s muons", type[i].c_str());
1273   h_Phi_Muon = new TH1D(name, title, 100, -3.1415926, 3.1415926);
1274 
1275   sprintf(name, "h_P_Muon_%s", type[i].c_str());
1276   sprintf(title, "p of %s muons (GeV)", type[i].c_str());
1277   h_P_Muon = new TH1D(name, title, 100, 0, 200);
1278 
1279   sprintf(name, "h_PF_Muon_%s", type[i].c_str());
1280   sprintf(title, "PF %s muons (GeV)", type[i].c_str());
1281   h_PF_Muon = new TH1D(name, title, 2, 0, 2);
1282 
1283   sprintf(name, "h_Global_Muon_Chi2_%s", type[i].c_str());
1284   sprintf(title, "Chi2 Global %s muons (GeV)", type[i].c_str());
1285   h_GlobTrack_Chi = new TH1D(name, title, 15, 0, 15);
1286 
1287   sprintf(name, "h_Global_Muon_Hits_%s", type[i].c_str());
1288   sprintf(title, "Global Hits %s muons (GeV)", type[i].c_str());
1289   h_Global_Muon_Hits = new TH1D(name, title, 10, 0, 10);
1290 
1291   sprintf(name, "h_Matched_Stations_%s", type[i].c_str());
1292   sprintf(title, "Matched Stations %s muons (GeV)", type[i].c_str());
1293   h_MatchedStations = new TH1D(name, title, 10, 0, 10);
1294 
1295   sprintf(name, "h_Transverse_ImpactParameter_%s", type[i].c_str());
1296   sprintf(title, "Transverse_ImpactParameter of %s muons (GeV)", type[i].c_str());
1297   h_Tight_TransImpactparameter = new TH1D(name, title, 50, 0, 10);
1298 
1299   sprintf(name, "h_Longitudinal_ImpactParameter_%s", type[i].c_str());
1300   sprintf(title, "Longitudinal_ImpactParameter of %s muons (GeV)", type[i].c_str());
1301   h_Tight_LongitudinalImpactparameter = new TH1D(name, title, 20, 0, 10);
1302 
1303   sprintf(name, "h_InnerTrack_PixelHits_%s", type[i].c_str());
1304   sprintf(title, "InnerTrack_PixelHits of %s muons (GeV)", type[i].c_str());
1305   h_InnerTrackPixelHits = new TH1D(name, title, 20, 0, 20);
1306 
1307   sprintf(name, "h_TrackLayers_%s", type[i].c_str());
1308   sprintf(title, "No. of Tracker Layers of %s muons (GeV)", type[i].c_str());
1309   h_TrackerLayer = new TH1D(name, title, 20, 0, 20);
1310 
1311   sprintf(name, "h_IsolationR04_%s", type[i].c_str());
1312   sprintf(title, "IsolationR04 %s muons (GeV)", type[i].c_str());
1313   h_IsolationR04 = new TH1D(name, title, 45, 0, 5);
1314 
1315   sprintf(name, "h_Global_Muon_%s", type[i].c_str());
1316   sprintf(title, "Global %s muons (GeV)", type[i].c_str());
1317   h_Global_Muon = new TH1D(name, title, 2, 0, 2);
1318 
1319   sprintf(name, "h_TransImpactParameter_%s", type[i].c_str());
1320   sprintf(title, "TransImpactParameter of %s muons (GeV)", type[i].c_str());
1321   h_TransImpactParameter = new TH1D(name, title, 100, 0, 0.5);
1322 
1323   sprintf(name, "h_TransImpactParameterBin1_%s", type[i].c_str());
1324   sprintf(title, "TransImpactParameter of %s muons (GeV) in -1.5 <= #phi <= 0.5", type[i].c_str());
1325   h_TransImpactParameterBin1 = new TH1D(name, title, 100, 0, 0.5);
1326 
1327   sprintf(name, "h_TransImpactParameterBin2_%s", type[i].c_str());
1328   sprintf(title, "TransImpactParameter of %s muons (GeV) in #phi> 0.5 and #phi< -1.5 ", type[i].c_str());
1329   h_TransImpactParameterBin2 = new TH1D(name, title, 100, 0, 0.5);
1330 
1331   sprintf(name, "h_LongImpactParameter_%s", type[i].c_str());
1332   sprintf(title, "LongImpactParameter of %s muons (GeV)", type[i].c_str());
1333   h_LongImpactParameter = new TH1D(name, title, 100, 0, 30);
1334 
1335   sprintf(name, "h_LongImpactParameterBin1_%s", type[i].c_str());
1336   sprintf(title, "LongImpactParameter of %s muons (GeV) in -1.5 <= #phi <= 0.5", type[i].c_str());
1337   h_LongImpactParameterBin1 = new TH1D(name, title, 100, 0, 30);
1338 
1339   sprintf(name, "h_LongImpactParameterBin2_%s", type[i].c_str());
1340   sprintf(title, "LongImpactParameter of %s muons (GeV) in #phi> 0.5 and #phi< -1.5 ", type[i].c_str());
1341   h_LongImpactParameterBin2 = new TH1D(name, title, 100, 0, 30);
1342 
1343   sprintf(name, "h_2D_Bin1_%s", type[i].c_str());
1344   sprintf(title, "Trans/Long ImpactParameter of %s muons (GeV) in -1.5 <= #phi< 0.5 ", type[i].c_str());
1345   h_2D_Bin1 = new TH2D(name, title, 100, 0, 0.5, 100, 0, 30);
1346 
1347   sprintf(name, "h_2D_Bin2_%s", type[i].c_str());
1348   sprintf(title, "Trans/Long ImpactParameter of %s muons (GeV) in #phi> 0.5 and #phi< -1.5 ", type[i].c_str());
1349   h_2D_Bin2 = new TH2D(name, title, 100, 0, 0.5, 100, 0, 30);
1350 
1351   sprintf(name, "h_ecal_energy_%s", type[i].c_str());
1352   sprintf(title, "ECAL energy for %s muons", type[i].c_str());
1353   h_ecal_energy = new TH1D(name, title, 1000, -10.0, 90.0);
1354 
1355   sprintf(name, "h_hcal_energy_%s", type[i].c_str());
1356   sprintf(title, "HCAL energy for %s muons", type[i].c_str());
1357   h_hcal_energy = new TH1D(name, title, 500, -10.0, 90.0);
1358 
1359   sprintf(name, "h_3x3_ecal_%s", type[i].c_str());
1360   sprintf(title, "ECAL energy in 3x3 for %s muons", type[i].c_str());
1361   h_3x3_ecal = new TH1D(name, title, 1000, -10.0, 90.0);
1362 
1363   sprintf(name, "h_1x1_hcal_%s", type[i].c_str());
1364   sprintf(title, "HCAL energy in 1x1 for %s muons", type[i].c_str());
1365   h_1x1_hcal = new TH1D(name, title, 500, -10.0, 90.0);
1366 
1367   sprintf(name, "h_EtaX_hcal_%s", type[i].c_str());
1368   sprintf(title, "HCAL energy as a function of i#eta for %s muons", type[i].c_str());
1369   h_EtaX_hcal = new TProfile(name, title, 60, -30.0, 30.0);
1370 
1371   sprintf(name, "h_PhiY_hcal_%s", type[i].c_str());
1372   sprintf(title, "HCAL energy as a function of i#phi for %s muons", type[i].c_str());
1373   h_PhiY_hcal = new TProfile(name, title, 72, 0, 72);
1374 
1375   sprintf(name, "h_EtaX_ecal_%s", type[i].c_str());
1376   sprintf(title, "EB energy as a function of i#eta for %s muons", type[i].c_str());
1377   h_EtaX_ecal = new TProfile(name, title, 170, -85.0, 85.0);
1378 
1379   sprintf(name, "h_PhiY_ecal_%s", type[i].c_str());
1380   sprintf(title, "EB energy as a function of i#phi for %s muons", type[i].c_str());
1381   h_PhiY_ecal = new TProfile(name, title, 360, 0, 360);
1382 
1383   sprintf(name, "h_Eta_ecal_%s", type[i].c_str());
1384   sprintf(title, "ECAL energy as a function of #eta for %s muons", type[i].c_str());
1385   h_Eta_ecal = new TProfile(name, title, 100, -2.5, 2.5);
1386 
1387   sprintf(name, "h_Phi_ecal_%s", type[i].c_str());
1388   sprintf(title, "ECAL energy as a function of #phi for %s muons", type[i].c_str());
1389   h_Phi_ecal = new TProfile(name, title, 100, -3.1415926, 3.1415926);
1390 
1391   sprintf(name, "h_MuonHittingEcal_%s", type[i].c_str());
1392   sprintf(title, "%s muons hitting ECAL", type[i].c_str());
1393   h_MuonHittingEcal = new TH1D(name, title, 100, 0, 5.0);
1394 
1395   sprintf(name, "h_HotCell_%s", type[i].c_str());
1396   sprintf(title, "Hot cell for %s muons", type[i].c_str());
1397   h_HotCell = new TH1D(name, title, 100, 0, 2);
1398 
1399   //        output_file->cd();
1400   for (int eta = etaMin_; eta <= etaMax_; ++eta) {
1401     int nDepth = nDepthBins(eta, -1);
1402     int nPhi = nPhiBins(eta);
1403     //sprintf(name, "Dir_muon_type_%s_ieta%d",type[i].c_str(), eta);
1404     //d_output_file[i][eta]= output_file->mkdir(name);
1405     //output_file->cd(name);
1406     //d_output_file[i][eta]->cd();
1407     for (int depth = 0; depth < nDepth; ++depth) {
1408       for (int PHI = 0; PHI < nPhi; ++PHI) {
1409         int PHI0 = (nPhi == 72) ? PHI + 1 : 2 * PHI + 1;
1410         int ih = indxEta[eta - 1][depth][PHI];
1411         if (debug_)
1412           std::cout << "eta:" << eta << " depth:" << depth << " PHI:" << PHI << ":" << PHI0 << " ih:" << ih
1413                     << std::endl;
1414 
1415         sprintf(name, "h_Hot_MuonEnergy_hc_%d_%d_%d_%s_HotCell", eta, (depth + 1), PHI0, type[i].c_str());
1416         sprintf(title,
1417                 "HCAL energy in hot tower (i#eta=%d, depth=%d, i#phi = %d) for extrapolated %s muons (Hot Cell)",
1418                 eta,
1419                 (depth + 1),
1420                 PHI0,
1421                 type[i].c_str());
1422         h_Hot_MuonEnergy_hcal_HotCell[ih] = new TH1D(name, title, 4000, 0.0, 10.0);
1423         h_Hot_MuonEnergy_hcal_HotCell[ih]->Sumw2();
1424 
1425         sprintf(
1426             name, "h_Hot_MuonEnergy_hc_%d_%d_%d_%s_HotCell_ByActiveLength", eta, (depth + 1), PHI0, type[i].c_str());
1427         sprintf(title,
1428                 "HCAL energy in hot tower (i#eta=%d, depth=%d, i#phi = %d) for extrapolated %s muons (Hot Cell) "
1429                 "divided by Active Length",
1430                 eta,
1431                 (depth + 1),
1432                 PHI0,
1433                 type[i].c_str());
1434         h_Hot_MuonEnergy_hcal_HotCell_VsActiveLength[ih] = new TH1D(name, title, 4000, 0.0, 10.0);
1435         h_Hot_MuonEnergy_hcal_HotCell_VsActiveLength[ih]->Sumw2();
1436 
1437         sprintf(name, "h_active_length_Fill_%d_%d_%d_%s", eta, (depth + 1), PHI0, type[i].c_str());
1438         sprintf(title, "active_length%d_%d_%d_%s", eta, (depth + 1), PHI0, type[i].c_str());
1439         h_active_length_Fill[ih] = new TH1D(name, title, 20, 0, 20);
1440         h_active_length_Fill[ih]->Sumw2();
1441 
1442         sprintf(name, "h_p_muon_in_%d_%d_%d_%s", eta, (depth + 1), PHI0, type[i].c_str());
1443         sprintf(title, "p_muon_in%d_%d_%d_%s", eta, (depth + 1), PHI0, type[i].c_str());
1444         h_p_muon_ineta[ih] = new TH1D(name, title, 500, 0, 500);
1445         h_p_muon_ineta[ih]->Sumw2();
1446 
1447         sprintf(name, "h_charge_signal_in_%d_%d_%d_%s", eta, (depth + 1), PHI0, type[i].c_str());
1448         sprintf(name, "charge_signal_in_%d_%d_%d_%s", eta, (depth + 1), PHI0, type[i].c_str());
1449         h_charge_signal[ih] = new TH1D(name, title, 500, 0, 500);
1450         h_charge_signal[ih]->Sumw2();
1451 
1452         sprintf(name, "h_charge_bg_in_%d_%d_%d_%s", eta, (depth + 1), PHI0, type[i].c_str());
1453         sprintf(name, "charge_bg_in_%d_%d_%d_%s", eta, (depth + 1), PHI0, type[i].c_str());
1454         h_charge_bg[ih] = new TH1D(name, title, 500, 0, 500);
1455         h_charge_bg[ih]->Sumw2();
1456 
1457         ih++;
1458 
1459         sprintf(name, "h_Hot_MuonEnergy_hc_%d_%d_%d_%s_HotCell", -eta, (depth + 1), PHI0, type[i].c_str());
1460         sprintf(title,
1461                 "HCAL energy in hot tower (i#eta=%d, depth=%d, i#phi = %d) for extrapolated %s muons (Hot Cell)",
1462                 -eta,
1463                 (depth + 1),
1464                 PHI0,
1465                 type[i].c_str());
1466         h_Hot_MuonEnergy_hcal_HotCell[ih] = new TH1D(name, title, 4000, 0.0, 10.0);
1467         h_Hot_MuonEnergy_hcal_HotCell[ih]->Sumw2();
1468 
1469         sprintf(
1470             name, "h_Hot_MuonEnergy_hc_%d_%d_%d_%s_HotCell_ByActiveLength", -eta, (depth + 1), PHI0, type[i].c_str());
1471         sprintf(title,
1472                 "HCAL energy in hot tower (i#eta=%d, depth=%d, i#phi=%d) for extrapolated %s muons (Hot Cell) "
1473                 "divided by Active Length",
1474                 -eta,
1475                 (depth + 1),
1476                 PHI0,
1477                 type[i].c_str());
1478         h_Hot_MuonEnergy_hcal_HotCell_VsActiveLength[ih] = new TH1D(name, title, 4000, 0.0, 10.0);
1479         h_Hot_MuonEnergy_hcal_HotCell_VsActiveLength[ih]->Sumw2();
1480 
1481         sprintf(name, "h_active_length_Fill_%d_%d_%d_%s", -eta, (depth + 1), PHI0, type[i].c_str());
1482         sprintf(title, "active_length%d_%d_%d_%s", -eta, (depth + 1), PHI0, type[i].c_str());
1483         h_active_length_Fill[ih] = new TH1D(name, title, 20, 0, 20);
1484         h_active_length_Fill[ih]->Sumw2();
1485 
1486         sprintf(name, "h_p_muon_in_%d_%d_%d_%s", -eta, (depth + 1), PHI0, type[i].c_str());
1487         sprintf(title, "p_muon_in%d_%d_%d_%s", -eta, (depth + 1), PHI0, type[i].c_str());
1488         h_p_muon_ineta[ih] = new TH1D(name, title, 500, 0, 500);
1489         h_p_muon_ineta[ih]->Sumw2();
1490 
1491         sprintf(name, "h_charge_signal_in_%d_%d_%d_%s", -eta, (depth + 1), PHI0, type[i].c_str());
1492         sprintf(name, "charge_signal_in_%d_%d_%d_%s", -eta, (depth + 1), PHI0, type[i].c_str());
1493         h_charge_signal[ih] = new TH1D(name, title, 500, 0, 500);
1494         h_charge_signal[ih]->Sumw2();
1495 
1496         sprintf(name, "h_charge_bg_in_%d_%d_%d_%s", -eta, (depth + 1), PHI0, type[i].c_str());
1497         sprintf(name, "charge_bg_in_%d_%d_%d_%s", -eta, (depth + 1), PHI0, type[i].c_str());
1498         h_charge_bg[ih] = new TH1D(name, title, 500, 0, 500);
1499         h_charge_bg[ih]->Sumw2();
1500       }
1501     }
1502     //output_file->cd();
1503   }
1504   //output_file->cd();
1505 }
1506 
1507 bool HBHEMuonOfflineAnalyzer::getEnergy(
1508     int dep, double &enb, double &enu, double &enh, double &enc, double &chgS, double &chgB, double &actL) {
1509   double cfac(1.0);
1510   bool flag(true);
1511   if (cFactor_) {
1512     int ieta = hcal_ieta;
1513     int iphi = hcal_iphi;
1514     unsigned int detId = getDetIdHBHE(ieta, iphi, dep + 1);
1515     cfac = getCorr(Run_No, detId);
1516   }
1517   if (dep == 0) {
1518     enb = cfac * hcal_edepth1;
1519     enu = cfac * hcal_edepthCorrect1;
1520     enh = cfac * hcal_edepthHot1;
1521     enc = cfac * hcal_edepthHotCorrect1;
1522     chgS = hcal_cdepthHot1;
1523     actL = hcal_activeHotL1;
1524     chgB = hcal_cdepthHotBG1;
1525   } else if (dep == 1) {
1526     enb = cfac * hcal_edepth2;
1527     enu = cfac * hcal_edepthCorrect2;
1528     enh = cfac * hcal_edepthHot2;
1529     enc = cfac * hcal_edepthHotCorrect2;
1530     chgS = hcal_cdepthHot2;
1531     actL = hcal_activeHotL2;
1532     chgB = hcal_cdepthHotBG2;
1533   } else if (dep == 2) {
1534     enb = cfac * hcal_edepth3;
1535     enu = cfac * hcal_edepthCorrect3;
1536     enh = cfac * hcal_edepthHot3;
1537     enc = cfac * hcal_edepthHotCorrect3;
1538     chgS = hcal_cdepthHot3;
1539     actL = hcal_activeHotL3;
1540     chgB = hcal_cdepthHotBG3;
1541   } else if (dep == 3) {
1542     enb = cfac * hcal_edepth4;
1543     enu = cfac * hcal_edepthCorrect4;
1544     enh = cfac * hcal_edepthHot4;
1545     enc = cfac * hcal_edepthHotCorrect4;
1546     chgS = hcal_cdepthHot4;
1547     actL = hcal_activeHotL4;
1548     chgB = hcal_cdepthHotBG4;
1549   } else if (dep == 4) {
1550     enb = cfac * hcal_edepth5;
1551     enu = cfac * hcal_edepthCorrect5;
1552     enh = cfac * hcal_edepthHot5;
1553     enc = cfac * hcal_edepthHotCorrect5;
1554     chgS = hcal_cdepthHot5;
1555     actL = hcal_activeHotL5;
1556     chgB = hcal_cdepthHotBG5;
1557   } else if (dep == 5) {
1558     if (dep <= maxDepth_) {
1559       enb = cfac * hcal_edepth6;
1560       enu = cfac * hcal_edepthCorrect6;
1561       enh = cfac * hcal_edepthHot6;
1562       enc = cfac * hcal_edepthHotCorrect6;
1563       chgS = hcal_cdepthHot6;
1564       actL = hcal_activeHotL6;
1565       chgB = hcal_cdepthHotBG6;
1566     } else {
1567       enb = enu = enh = enc = chgS = actL = chgB = 0;
1568       flag = false;
1569     }
1570   } else if (dep == 6) {
1571     if (dep <= maxDepth_) {
1572       enb = cfac * hcal_edepth7;
1573       enu = cfac * hcal_edepthCorrect7;
1574       enh = cfac * hcal_edepthHot7;
1575       enc = cfac * hcal_edepthHotCorrect7;
1576       chgS = hcal_cdepthHot7;
1577       actL = hcal_activeHotL7;
1578       chgB = hcal_cdepthHotBG7;
1579     } else {
1580       enb = enu = enh = enc = chgS = actL = chgB = 0;
1581       flag = false;
1582     }
1583   } else {
1584     enb = enu = enh = enc = chgS = actL = chgB = 0;
1585     flag = false;
1586   }
1587   return flag;
1588 }
1589 
1590 bool HBHEMuonOfflineAnalyzer::looseMuon() {
1591   if (pt_of_muon > 20.) {
1592     if (mediumMuon2016()) {
1593       if (IsolationR04 < 0.25) {
1594         return true;
1595       }
1596     }
1597   }
1598   return false;
1599 }
1600 
1601 bool HBHEMuonOfflineAnalyzer::softMuon() {
1602   if (pt_of_muon > 20.) {
1603     if (mediumMuon2016()) {
1604       if (IsolationR03 < 0.10) {
1605         return true;
1606       }
1607     }
1608   }
1609   return false;
1610 }
1611 
1612 bool HBHEMuonOfflineAnalyzer::tightMuon() {
1613   if (pt_of_muon > 20.) {
1614     if (mediumMuon2016()) {
1615       if (IsolationR04 < 0.15) {
1616         return true;
1617       }
1618     }
1619   }
1620   return false;
1621 }
1622 
1623 bool HBHEMuonOfflineAnalyzer::mediumMuon2016() {
1624   bool medium16 = (((PF_Muon) && (Global_Muon || Tracker_muon)) && (tight_validFraction > 0.49));
1625   if (!medium16)
1626     return medium16;
1627 
1628   bool goodGlob = (Global_Muon && GlobTrack_Chi < 3 && muon_chi2LocalPosition < 12 && muon_trkKink < 20);
1629   medium16 = muon_segComp > (goodGlob ? 0.303 : 0.451);
1630   return medium16;
1631 }
1632 
1633 void HBHEMuonOfflineAnalyzer::etaPhiHcal(unsigned int detId, int &eta, int &phi, int &depth) {
1634   int zside, etaAbs;
1635   if ((detId & 0x1000000) == 0) {
1636     zside = (detId & 0x2000) ? (1) : (-1);
1637     etaAbs = (detId >> 7) & 0x3F;
1638     phi = detId & 0x7F;
1639     depth = (detId >> 14) & 0x1F;
1640   } else {
1641     zside = (detId & 0x80000) ? (1) : (-1);
1642     etaAbs = (detId >> 10) & 0x1FF;
1643     phi = detId & 0x3FF;
1644     depth = (detId >> 20) & 0xF;
1645   }
1646   eta = etaAbs * zside;
1647 }
1648 
1649 void HBHEMuonOfflineAnalyzer::etaPhiEcal(
1650     unsigned int detId, int &type, int &zside, int &etaX, int &phiY, int &plane, int &strip) {
1651   type = ((detId >> 25) & 0x7);
1652   // std::cout<<"type"<<type<<std::endl;
1653   plane = strip = 0;
1654   if (type == 1) {
1655     //Ecal Barrel
1656     zside = (detId & 0x10000) ? (1) : (-1);
1657     etaX = (detId >> 9) & 0x7F;
1658     phiY = detId & 0x1FF;
1659   } else if (type == 2) {
1660     zside = (detId & 0x4000) ? (1) : (-1);
1661     etaX = (detId >> 7) & 0x7F;
1662     phiY = (detId & 0x7F);
1663   } else if (type == 3) {
1664     zside = (detId & 0x80000) ? (1) : (-1);
1665     etaX = (detId >> 6) & 0x3F;
1666     /** get the sensor iy */
1667     phiY = (detId >> 12) & 0x3F;
1668     /** get the strip */
1669     plane = ((detId >> 18) & 0x1) + 1;
1670     strip = detId & 0x3F;
1671   } else {
1672     zside = etaX = phiY = 0;
1673   }
1674 }
1675 
1676 void HBHEMuonOfflineAnalyzer::calculateP(double pt, double eta, double &pM) {
1677   pM = (pt * cos(2 * (1 / atan(exp(eta)))));
1678 }
1679 
1680 void HBHEMuonOfflineAnalyzer::close() {
1681   output_file->cd();
1682   std::cout << "file yet to be Written" << std::endl;
1683   outtree_->Write();
1684   writeHistograms();
1685   std::cout << "file Written" << std::endl;
1686   output_file->Close();
1687   std::cout << "now doing return" << std::endl;
1688 }
1689 
1690 void HBHEMuonOfflineAnalyzer::writeHistograms() {
1691   //output_file->cd();
1692   std::string type[] = {"tight", "soft", "loose"};
1693   char name[128];
1694 
1695   std::cout << "WriteHistograms" << std::endl;
1696   nHist = 0;
1697   for (int eta = etaMin_; eta <= etaMax_; ++eta) {
1698     int nDepth = nDepthBins(eta, -1);
1699     int nPhi = nPhiBins(eta);
1700     if (debug_)
1701       std::cout << "Eta:" << eta << " nDepths " << nDepth << " nPhis " << nPhi << std::endl;
1702     for (int depth = 0; depth < nDepth; ++depth) {
1703       if (debug_)
1704         std::cout << "Eta:" << eta << "Depth:" << depth << std::endl;
1705       for (int PHI = 0; PHI < nPhi; ++PHI) {
1706         indxEta[eta - 1][depth][PHI] = nHist;
1707         nHist += 2;
1708       }
1709     }
1710   }
1711 
1712   TDirectory *d_output_file[29];
1713   h_evtype->Write();
1714   //output_file->cd();
1715   int i(nCut_);
1716   h_Pt_Muon->Write();
1717   h_Eta_Muon->Write();
1718   h_Phi_Muon->Write();
1719 
1720   h_P_Muon->Write();
1721   h_PF_Muon->Write();
1722   h_GlobTrack_Chi->Write();
1723   h_Global_Muon_Hits->Write();
1724 
1725   h_MatchedStations->Write();
1726   h_Tight_TransImpactparameter->Write();
1727   h_Tight_LongitudinalImpactparameter->Write();
1728 
1729   h_InnerTrackPixelHits->Write();
1730   h_TrackerLayer->Write();
1731   h_IsolationR04->Write();
1732 
1733   h_Global_Muon->Write();
1734   h_TransImpactParameter->Write();
1735 
1736   h_TransImpactParameterBin1->Write();
1737   h_TransImpactParameterBin2->Write();
1738 
1739   h_LongImpactParameter->Write();
1740   h_LongImpactParameterBin1->Write();
1741   h_LongImpactParameterBin2->Write();
1742 
1743   h_ecal_energy->Write();
1744   h_hcal_energy->Write();
1745 
1746   h_3x3_ecal->Write();
1747   h_1x1_hcal->Write();
1748 
1749   h_EtaX_hcal->Write();
1750   h_PhiY_hcal->Write();
1751 
1752   h_EtaX_ecal->Write();
1753   h_PhiY_ecal->Write();
1754 
1755   h_Eta_ecal->Write();
1756   h_Phi_ecal->Write();
1757   h_MuonHittingEcal->Write();
1758   h_HotCell->Write();
1759 
1760   output_file->cd();
1761   for (int eta = etaMin_; eta <= etaMax_; ++eta) {
1762     int nDepth = nDepthBins(eta, -1);
1763     int nPhi = nPhiBins(eta);
1764     sprintf(name, "Dir_muon_type_%s_ieta%d", type[i].c_str(), eta);
1765     d_output_file[eta - 1] = output_file->mkdir(name);
1766     //output_file->cd(name);
1767     d_output_file[eta - 1]->cd();
1768     for (int depth = 0; depth < nDepth; ++depth) {
1769       for (int PHI = 0; PHI < nPhi; ++PHI) {
1770         //  std::cout<<"eta:"<<eta<<"depth:"<<depth<<"PHI:"<<PHI<<std::endl;
1771         int ih = indxEta[eta - 1][depth][PHI];
1772         //  std::cout<<"ih:"<<ih<<std::endl;
1773         h_Hot_MuonEnergy_hcal_HotCell[ih]->Write();
1774 
1775         h_Hot_MuonEnergy_hcal_HotCell_VsActiveLength[ih]->Write();
1776         h_active_length_Fill[ih]->Write();
1777         h_p_muon_ineta[ih]->Write();
1778         h_charge_signal[ih]->Write();
1779         h_charge_bg[ih]->Write();
1780         ih++;
1781         h_Hot_MuonEnergy_hcal_HotCell[ih]->Write();
1782         h_Hot_MuonEnergy_hcal_HotCell_VsActiveLength[ih]->Write();
1783         h_active_length_Fill[ih]->Write();
1784         h_p_muon_ineta[ih]->Write();
1785         h_charge_signal[ih]->Write();
1786         h_charge_bg[ih]->Write();
1787       }
1788     }
1789     output_file->cd();
1790   }
1791   output_file->cd();
1792 }
1793 
1794 int HBHEMuonOfflineAnalyzer::nDepthBins(int eta, int phi) {
1795   // Run 1 scenario
1796   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};
1797   // Run 2 scenario from 2018
1798   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};
1799   // Run 3 scenario
1800   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};
1801   // Run 4 scenario
1802   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};
1803   // for a test scenario with multi depth segmentation considered during Run 1
1804   //    int  nDepth[29]={3,3,3,3,3,3,3,3,3,3,3,3,3,3,3,3,5,5,5,5,5,5,5,5,5,5,5,5,5};
1805   // modeLHC_ = 0 --> nbin defined maxDepthHB/HE
1806   //          = 1 -->      corresponds to Run 1 (valid till 2016)
1807   //          = 2 -->      corresponds to Run 2 (2018 geometry)
1808   //          = 3 -->      corresponds to Run 3 (post LS2)
1809   //          = 4 -->      corresponds to 2017 (Plan 1)
1810   //          = 5 -->      corresponds to Run 4 (post LS3)
1811   int nbin(0);
1812   if (modeLHC_ == 0) {
1813     if (eta <= 15) {
1814       nbin = maxDepthHB_;
1815     } else if (eta == 16) {
1816       nbin = 4;
1817     } else {
1818       nbin = maxDepthHE_;
1819     }
1820   } else if (modeLHC_ == 1) {
1821     nbin = nDepthR1[eta - 1];
1822   } else if (modeLHC_ == 2) {
1823     nbin = nDepthR2[eta - 1];
1824   } else if (modeLHC_ == 3) {
1825     nbin = nDepthR3[eta - 1];
1826   } else if (modeLHC_ == 4) {
1827     if (phi > 0) {
1828       if (eta >= 16 && phi >= 63 && phi <= 66) {
1829         nbin = nDepthR2[eta - 1];
1830       } else {
1831         nbin = nDepthR1[eta - 1];
1832       }
1833     } else {
1834       if (eta >= 16) {
1835         nbin = (nDepthR2[eta - 1] > nDepthR1[eta - 1]) ? nDepthR2[eta - 1] : nDepthR1[eta - 1];
1836       } else {
1837         nbin = nDepthR1[eta - 1];
1838       }
1839     }
1840   } else {
1841     if (eta > 0 && eta < 30) {
1842       nbin = nDepthR4[eta - 1];
1843     } else {
1844       nbin = nDepthR4[28];
1845     }
1846   }
1847   return nbin;
1848 }
1849 
1850 int HBHEMuonOfflineAnalyzer::nPhiBins(int eta) {
1851   int nphi = (eta <= 20) ? 72 : 36;
1852   if (modeLHC_ == 5 && eta > 16)
1853     nphi = 360;
1854   return nphi;
1855 }
1856 
1857 float HBHEMuonOfflineAnalyzer::getCorr(int run, unsigned int id) {
1858   float cfac(1.0);
1859   int ip(-1);
1860   for (unsigned int k = 0; k < runlow_.size(); ++k) {
1861     unsigned int i = runlow_.size() - k - 1;
1862     if (run >= runlow_[i]) {
1863       ip = (int)(i);
1864       break;
1865     }
1866   }
1867   if (debug_) {
1868     std::cout << "Run " << run << " Perdiod " << ip << std::endl;
1869   }
1870   unsigned idx = correctDetId(id);
1871   if (ip >= 0) {
1872     std::map<unsigned int, float>::iterator itr = corrFac_[ip].find(idx);
1873     if (itr != corrFac_[ip].end())
1874       cfac = itr->second;
1875   }
1876   if (debug_) {
1877     int subdet, zside, ieta, iphi, depth;
1878     unpackDetId(idx, subdet, zside, ieta, iphi, depth);
1879     std::cout << "ID " << std::hex << id << std::dec << " (Sub " << subdet << " eta " << zside * ieta << " phi " << iphi
1880               << " depth " << depth << ")  Factor " << cfac << std::endl;
1881   }
1882   return cfac;
1883 }
1884 
1885 std::vector<std::string> HBHEMuonOfflineAnalyzer::splitString(const std::string &fLine) {
1886   std::vector<std::string> result;
1887   int start = 0;
1888   bool empty = true;
1889   for (unsigned i = 0; i <= fLine.size(); i++) {
1890     if (fLine[i] == ' ' || i == fLine.size()) {
1891       if (!empty) {
1892         std::string item(fLine, start, i - start);
1893         result.push_back(item);
1894         empty = true;
1895       }
1896       start = i + 1;
1897     } else {
1898       if (empty)
1899         empty = false;
1900     }
1901   }
1902   return result;
1903 }
1904 
1905 unsigned int HBHEMuonOfflineAnalyzer::getDetIdHBHE(int ieta, int iphi, int depth) {
1906   int eta = std::abs(ieta);
1907   int subdet(0);
1908   if (eta > 16)
1909     subdet = 2;
1910   else if (eta == 16 && depth > 2)
1911     subdet = 2;
1912   else
1913     subdet = 1;
1914   return getDetId(subdet, ieta, iphi, depth);
1915 }
1916 
1917 unsigned int HBHEMuonOfflineAnalyzer::getDetId(int subdet, int ieta, int iphi, int depth) {
1918   // All numbers used here are described as masks/offsets in
1919   // DataFormats/HcalDetId/interface/HcalDetId.h
1920   unsigned int id = ((4 << 28) | ((subdet & 0x7) << 25));
1921   id |= ((0x1000000) | ((depth & 0xF) << 20) | ((ieta > 0) ? (0x80000 | (ieta << 10)) : ((-ieta) << 10)) |
1922          (iphi & 0x3FF));
1923   return id;
1924 }
1925 
1926 unsigned int HBHEMuonOfflineAnalyzer::correctDetId(const unsigned int &detId) {
1927   int subdet, ieta, zside, depth, iphi;
1928   unpackDetId(detId, subdet, zside, ieta, iphi, depth);
1929   if (subdet == 0) {
1930     if (ieta > 16)
1931       subdet = 2;
1932     else if (ieta == 16 && depth > 2)
1933       subdet = 2;
1934     else
1935       subdet = 1;
1936   }
1937   unsigned int id = getDetId(subdet, ieta * zside, iphi, depth);
1938   if ((id != detId) && debug_) {
1939     std::cout << "Correct Id " << std::hex << detId << " to " << id << std::dec << "(Sub " << subdet << " eta "
1940               << ieta * zside << " phi " << iphi << " depth " << depth << ")" << std::endl;
1941   }
1942   return id;
1943 }
1944 
1945 void HBHEMuonOfflineAnalyzer::unpackDetId(
1946     unsigned int detId, int &subdet, int &zside, int &ieta, int &iphi, int &depth) {
1947   // The maskings are defined in DataFormats/DetId/interface/DetId.h
1948   //                      and in DataFormats/HcalDetId/interface/HcalDetId.h
1949   // The macro does not invoke the classes there and use them
1950   subdet = ((detId >> 25) & (0x7));
1951   if ((detId & 0x1000000) == 0) {
1952     ieta = ((detId >> 7) & 0x3F);
1953     zside = (detId & 0x2000) ? (1) : (-1);
1954     depth = ((detId >> 14) & 0x1F);
1955     iphi = (detId & 0x3F);
1956   } else {
1957     ieta = ((detId >> 10) & 0x1FF);
1958     zside = (detId & 0x80000) ? (1) : (-1);
1959     depth = ((detId >> 20) & 0xF);
1960     iphi = (detId & 0x3FF);
1961   }
1962 }