Back to home page

Project CMSSW displayed by LXR
 
 

    

File indexing completed on 2021-10-29 22:50:49 

0001 //////////////////////////////////////////////////////////////////////////////
0002 // Usage:
0003 // .L CalibPlotProperties.C+g
0004 //  CalibPlotProperties c1(fname, dirname, dupFileName, prefix, corrFileName,
0005 //                     rcorFileName, puCorr, flag, dataMC, truncateFlag,
0006 //                         useGen, scale, useScale, etalo, etahi, runlo, runhi,
0007 //                         phimin, phimax, zside, nvxlo, nvxhi, rbx, exclude,
0008 //                         etamax);
0009 //  c1.Loop();
0010 //  c1.savePlot(histFileName, append, all, debug);
0011 //
0012 //        This will prepare a set of histograms with properties of the tracks
0013 //        which can be displayed by the method in this file
0014 //
0015 //  PlotHist(histFileName, prefix, text, flagC, etalo, etahi, save)
0016 //
0017 //        This will plot the heistograms and save the canvases
0018 //
0019 //
0020 //   where:
0021 //
0022 //   fname   (std::string)     = file name of the input ROOT tree
0023 //                               or name of the file containing a list of
0024 //                               file names of input ROOT trees
0025 //   dirname (std::string)     = name of the directory where Tree resides
0026 //                               (use "HcalIsoTrkAnalyzer")
0027 //   dupFileName (char*)       = name of the file containing list of entries
0028 //                               of duplicate events
0029 //   prefix (std::string)      = String to be added to the name of histogram
0030 //                               (usually a 4 character string; default="")
0031 //   corrFileName (char*)      = name of the text file having the correction
0032 //                               factors to be used (default="", no corr.)
0033 //   rcorFileName (char*)      = name of the text file having the correction
0034 //                               factors as a function of run numbers or depth
0035 //                               to be used for raddam/depth dependent
0036 //                               correction  (default="", no corr.)
0037 //   puCorr (int)              = PU correction to be applied or not: 0 no
0038 //                               correction; < 0 use eDelta; > 0 rho dependent
0039 //                               correction (-2)
0040 //   flag (int)                = 6 digit integer (mlthdo) with control
0041 //                               information (m=0/1 for controlling creation
0042 //                               of depth depedendent histograms;
0043 //                               l=2/1/0 for type of rcorFileName (2 for overall
0044 //                               response corrections; 1 for depth dependence
0045 //                               corrections; 0 for raddam corrections);
0046 //                               t = bit information (lower bit set will
0047 //                               apply a cut on L1 closeness; and higher bit
0048 //                               set read correction file with Marina format);
0049 //                               h =0/1 flag to create energy histograms
0050 //                               d =0/1 flag to create basic set of histograms;
0051 //                               o =0/1/2 for tight / loose / flexible
0052 //                               selection). Default = 101111
0053 //   dataMC (bool)             = true/false for data/MC (default true)
0054 //   truncateFlag    (int)     = Flag to treat different depths differently (0)
0055 //                               both depths of ieta 15, 16 of HB as depth 1 (1)
0056 //                               all depths as depth 1 (2), all depths in HE
0057 //                               with values > 1 as depth 2 (3), all depths in
0058 //                               HB with values > 1 as depth 2 (4), all depths
0059 //                               in HB and HE with values > 1 as depth 2 (5)
0060 //                               (Default 0)
0061 //   useGen (bool)             = true/false to use generator level momentum
0062 //                               or reconstruction level momentum (def false)
0063 //   scale (double)            = energy scale if correction factor to be used
0064 //                               (default = 1.0)
0065 //   useScale (int)            = application of scale factor (0: nowehere,
0066 //                               1: barrel; 2: endcap, 3: everywhere)
0067 //                               barrel => |ieta| < 16; endcap => |ieta| > 15
0068 //   etalo/etahi (int,int)     = |eta| ranges (0:30)
0069 //   runlo  (int)              = lower value of run number to be included (+ve)
0070 //                               or excluded (-ve) (default 0)
0071 //   runhi  (int)              = higher value of run number to be included
0072 //                               (+ve) or excluded (-ve) (def 9999999)
0073 //   phimin          (int)     = minimum iphi value (1)
0074 //   phimax          (int)     = maximum iphi value (72)
0075 //   zside           (int)     = the side of the detector if phimin and phimax
0076 //                               differ from 1-72 (1)
0077 //   nvxlo           (int)     = minimum # of vertex in event to be used (0)
0078 //   nvxhi           (int)     = maximum # of vertex in event to be used (1000)
0079 //   rbx             (int)     = zside*(Subdet*100+RBX #) to be consdered (0)
0080 //   exclude         (bool)    = RBX specified by *rbx* to be exluded or only
0081 //                               considered (false)
0082 //   etamax          (bool)    = if set and if the corr-factor not found in the
0083 //                               corrFactor table, the corr-factor for the
0084 //                               corresponding zside, depth=1 and maximum ieta
0085 //                               in the table is taken (false)
0086 //
0087 //   histFileName (std::string)= name of the file containing saved histograms
0088 //   append (bool)             = true/false if the histogram file to be opened
0089 //                               in append/output mode
0090 //   all (bool)                = true/false if all histograms to be saved or
0091 //                               not (def false)
0092 //
0093 //   text  (string)            = string to be put in the title
0094 //   flagC (int)               = 3 digit integer (hdo) with control
0095 //                               information (h=0/1 for plottting the depth
0096 //                               depedendent histograms;
0097 //                               d =0/1 flag to plot energy histograms;
0098 //                               o =0/1 flag to plot basic set of histograms;
0099 //                               Default = 111
0100 //   save (int)                = flag to create or not save the plot in a file
0101 //                               (0 = no save, > 0 pdf file, < 0 hpg file)
0102 //////////////////////////////////////////////////////////////////////////////
0103 #include <TROOT.h>
0104 #include <TChain.h>
0105 #include <TFile.h>
0106 #include <TF1.h>
0107 #include <TH1D.h>
0108 #include <TH2F.h>
0109 #include <TProfile.h>
0110 #include <TFitResult.h>
0111 #include <TFitResultPtr.h>
0112 #include <TStyle.h>
0113 #include <TCanvas.h>
0114 #include <TLegend.h>
0115 #include <TPaveStats.h>
0116 #include <TPaveText.h>
0117 
0118 #include <algorithm>
0119 #include <iomanip>
0120 #include <iostream>
0121 #include <fstream>
0122 #include <map>
0123 #include <vector>
0124 #include <string>
0125 
0126 #include "CalibCorr.C"
0127 
0128 namespace CalibPlots {
0129   static const int npbin = 4;
0130   static const int netabin = 4;
0131   static const int ndepth = 7;
0132   static const int ntitles = 5;
0133   static const int npbin0 = (npbin + 1);
0134   int getP(int k) {
0135     const int ipbin[npbin0] = {20, 30, 40, 60, 100};
0136     return ((k >= 0 && k < npbin0) ? ipbin[k] : 0);
0137   }
0138   double getMomentum(int k) { return (double)(getP(k)); }
0139   int getEta(int k) {
0140     int ietas[netabin] = {0, 13, 18, 23};
0141     return ((k >= 0 && k < netabin) ? ietas[k] : -1);
0142   }
0143   std::string getTitle(int k) {
0144     std::string titl[ntitles] = {
0145         "all tracks", "good quality tracks", "selected tracks", "isolated good tracks", "tracks having MIP in ECAL"};
0146     return ((k >= 0 && k < ntitles) ? titl[k] : "");
0147   }
0148 }  // namespace CalibPlots
0149 
0150 class CalibPlotProperties {
0151 public:
0152   TChain *fChain;  //!pointer to the analyzed TTree or TChain
0153   Int_t fCurrent;  //!current Tree number in a TChain
0154 
0155   // Declaration of leaf types
0156   Int_t t_Run;
0157   Int_t t_Event;
0158   Int_t t_DataType;
0159   Int_t t_ieta;
0160   Int_t t_iphi;
0161   Double_t t_EventWeight;
0162   Int_t t_nVtx;
0163   Int_t t_nTrk;
0164   Int_t t_goodPV;
0165   Double_t t_l1pt;
0166   Double_t t_l1eta;
0167   Double_t t_l1phi;
0168   Double_t t_l3pt;
0169   Double_t t_l3eta;
0170   Double_t t_l3phi;
0171   Double_t t_p;
0172   Double_t t_pt;
0173   Double_t t_phi;
0174   Double_t t_mindR1;
0175   Double_t t_mindR2;
0176   Double_t t_eMipDR;
0177   Double_t t_eHcal;
0178   Double_t t_eHcal10;
0179   Double_t t_eHcal30;
0180   Double_t t_hmaxNearP;
0181   Double_t t_rhoh;
0182   Bool_t t_selectTk;
0183   Bool_t t_qltyFlag;
0184   Bool_t t_qltyMissFlag;
0185   Bool_t t_qltyPVFlag;
0186   Double_t t_gentrackP;
0187   std::vector<unsigned int> *t_DetIds;
0188   std::vector<double> *t_HitEnergies;
0189   std::vector<bool> *t_trgbits;
0190   std::vector<unsigned int> *t_DetIds1;
0191   std::vector<unsigned int> *t_DetIds3;
0192   std::vector<double> *t_HitEnergies1;
0193   std::vector<double> *t_HitEnergies3;
0194 
0195   // List of branches
0196   TBranch *b_t_Run;           //!
0197   TBranch *b_t_Event;         //!
0198   TBranch *b_t_DataType;      //!
0199   TBranch *b_t_ieta;          //!
0200   TBranch *b_t_iphi;          //!
0201   TBranch *b_t_EventWeight;   //!
0202   TBranch *b_t_nVtx;          //!
0203   TBranch *b_t_nTrk;          //!
0204   TBranch *b_t_goodPV;        //!
0205   TBranch *b_t_l1pt;          //!
0206   TBranch *b_t_l1eta;         //!
0207   TBranch *b_t_l1phi;         //!
0208   TBranch *b_t_l3pt;          //!
0209   TBranch *b_t_l3eta;         //!
0210   TBranch *b_t_l3phi;         //!
0211   TBranch *b_t_p;             //!
0212   TBranch *b_t_pt;            //!
0213   TBranch *b_t_phi;           //!
0214   TBranch *b_t_mindR1;        //!
0215   TBranch *b_t_mindR2;        //!
0216   TBranch *b_t_eMipDR;        //!
0217   TBranch *b_t_eHcal;         //!
0218   TBranch *b_t_eHcal10;       //!
0219   TBranch *b_t_eHcal30;       //!
0220   TBranch *b_t_hmaxNearP;     //!
0221   TBranch *b_t_rhoh;          //!
0222   TBranch *b_t_selectTk;      //!
0223   TBranch *b_t_qltyFlag;      //!
0224   TBranch *b_t_qltyMissFlag;  //!
0225   TBranch *b_t_qltyPVFlag;    //!
0226   TBranch *b_t_gentrackP;     //!
0227   TBranch *b_t_DetIds;        //!
0228   TBranch *b_t_HitEnergies;   //!
0229   TBranch *b_t_trgbits;       //!
0230   TBranch *b_t_DetIds1;       //!
0231   TBranch *b_t_DetIds3;       //!
0232   TBranch *b_t_HitEnergies1;  //!
0233   TBranch *b_t_HitEnergies3;  //!
0234 
0235   CalibPlotProperties(const char *fname,
0236                       const std::string &dirname,
0237                       const char *dupFileName,
0238                       const std::string &prefix = "",
0239                       const char *corrFileName = "",
0240                       const char *rcorFileName = "",
0241                       int puCorr = -1,
0242                       int flag = 101111,
0243                       bool dataMC = true,
0244                       int truncateFlag = 0,
0245                       bool useGen = false,
0246                       double scale = 1.0,
0247                       int useScale = 0,
0248                       int etalo = 0,
0249                       int etahi = 30,
0250                       int runlo = 0,
0251                       int runhi = 99999999,
0252                       int phimin = 1,
0253                       int phimax = 72,
0254                       int zside = 1,
0255                       int nvxlo = 0,
0256                       int nvxhi = 1000,
0257                       int rbx = 0,
0258                       bool exclude = false,
0259                       bool etamax = false);
0260   virtual ~CalibPlotProperties();
0261   virtual Int_t Cut(Long64_t entry);
0262   virtual Int_t GetEntry(Long64_t entry);
0263   virtual Long64_t LoadTree(Long64_t entry);
0264   virtual void Init(TChain *, const char *);
0265   virtual void Loop(Long64_t nentries = -1);
0266   virtual Bool_t Notify();
0267   virtual void Show(Long64_t entry = -1);
0268   bool goodTrack(double &eHcal, double &cut, bool debug);
0269   bool selectPhi(bool debug);
0270   void savePlot(const std::string &theName, bool append, bool all = false, bool debug = false);
0271   void correctEnergy(double &ener);
0272 
0273 private:
0274   static const int kp50 = 2;
0275   CalibCorrFactor *corrFactor_;
0276   CalibCorr *cFactor_;
0277   CalibSelectRBX *cSelect_;
0278   const std::string fname_, dirnm_, prefix_, outFileName_;
0279   const int corrPU_, flag_;
0280   const bool dataMC_, useGen_;
0281   const int truncateFlag_;
0282   const int etalo_, etahi_;
0283   int runlo_, runhi_;
0284   const int phimin_, phimax_, zside_, nvxlo_, nvxhi_, rbx_;
0285   bool exclude_, corrE_, cutL1T_;
0286   bool includeRun_, getHist_;
0287   int flexibleSelect_;
0288   bool plotBasic_, plotEnergy_, plotHists_;
0289   double log2by18_;
0290   std::ofstream fileout_;
0291   std::vector<Long64_t> entries_;
0292   std::vector<std::pair<int, int> > events_;
0293   TH1D *h_p[CalibPlots::ntitles];
0294   TH1D *h_eta[CalibPlots::ntitles], *h_nvtx;
0295   std::vector<TH1D *> h_eta0, h_eta1, h_eta2, h_eta3, h_eta4;
0296   std::vector<TH1D *> h_dL1, h_vtx;
0297   std::vector<TH1D *> h_etaEH[CalibPlots::npbin0];
0298   std::vector<TH1D *> h_etaEp[CalibPlots::npbin0];
0299   std::vector<TH1D *> h_etaEE[CalibPlots::npbin0];
0300   std::vector<TH1D *> h_etaEE0[CalibPlots::npbin0];
0301   std::vector<TH1D *> h_mom, h_eEcal, h_eHcal;
0302   std::vector<TH1F *> h_bvlist, h_bvlist2, h_evlist, h_evlist2;
0303   std::vector<TH1F *> h_bvlist3, h_evlist3;
0304   TH2F *h_etaE;
0305 };
0306 
0307 CalibPlotProperties::CalibPlotProperties(const char *fname,
0308                                          const std::string &dirnm,
0309                                          const char *dupFileName,
0310                                          const std::string &prefix,
0311                                          const char *corrFileName,
0312                                          const char *rcorFileName,
0313                                          int puCorr,
0314                                          int flag,
0315                                          bool dataMC,
0316                                          int truncate,
0317                                          bool useGen,
0318                                          double scl,
0319                                          int useScale,
0320                                          int etalo,
0321                                          int etahi,
0322                                          int runlo,
0323                                          int runhi,
0324                                          int phimin,
0325                                          int phimax,
0326                                          int zside,
0327                                          int nvxlo,
0328                                          int nvxhi,
0329                                          int rbx,
0330                                          bool exc,
0331                                          bool etam)
0332     : corrFactor_(nullptr),
0333       cFactor_(nullptr),
0334       cSelect_(nullptr),
0335       fname_(fname),
0336       dirnm_(dirnm),
0337       prefix_(prefix),
0338       corrPU_(puCorr),
0339       flag_(flag),
0340       dataMC_(dataMC),
0341       useGen_(useGen),
0342       truncateFlag_(truncate),
0343       etalo_(etalo),
0344       etahi_(etahi),
0345       runlo_(runlo),
0346       runhi_(runhi),
0347       phimin_(phimin),
0348       phimax_(phimax),
0349       zside_(zside),
0350       nvxlo_(nvxlo),
0351       nvxhi_(nvxhi),
0352       rbx_(rbx),
0353       exclude_(exc),
0354       includeRun_(true) {
0355   // if parameter tree is not specified (or zero), connect the file
0356   // used to generate this class and read the Tree
0357 
0358   flexibleSelect_ = ((flag_ / 1) % 10);
0359   plotBasic_ = (((flag_ / 10) % 10) > 0);
0360   plotEnergy_ = (((flag_ / 10) % 10) > 0);
0361   int oneplace = ((flag_ / 1000) % 10);
0362   cutL1T_ = (oneplace % 2);
0363   bool marina = ((oneplace / 2) % 2);
0364   bool ifDepth = (((flag_ / 10000) % 10) > 0);
0365   plotHists_ = (((flag_ / 100000) % 10) > 0);
0366   log2by18_ = std::log(2.5) / 18.0;
0367   if (runlo_ < 0 || runhi_ < 0) {
0368     runlo_ = std::abs(runlo_);
0369     runhi_ = std::abs(runhi_);
0370     includeRun_ = false;
0371   }
0372   char treeName[400];
0373   sprintf(treeName, "%s/CalibTree", dirnm.c_str());
0374   TChain *chain = new TChain(treeName);
0375   std::cout << "Create a chain for " << treeName << " from " << fname << " flags " << flexibleSelect_ << "|"
0376             << plotBasic_ << "|"
0377             << "|" << plotEnergy_ << "|" << plotHists_ << "|" << corrPU_ << " cons " << log2by18_ << " eta range "
0378             << etalo_ << ":" << etahi_ << " run range " << runlo_ << ":" << runhi_ << " (inclusion flag " << includeRun_
0379             << ") Vertex Range " << nvxlo_ << ":" << nvxhi_ << std::endl;
0380   corrFactor_ = new CalibCorrFactor(corrFileName, useScale, scl, etam, marina, false);
0381   if (!fillChain(chain, fname)) {
0382     std::cout << "*****No valid tree chain can be obtained*****" << std::endl;
0383   } else {
0384     std::cout << "Proceed with a tree chain with " << chain->GetEntries() << " entries" << std::endl;
0385     Init(chain, dupFileName);
0386     if (std::string(rcorFileName) != "")
0387       cFactor_ = new CalibCorr(rcorFileName, ifDepth, false);
0388     if (rbx != 0)
0389       cSelect_ = new CalibSelectRBX(rbx, false);
0390   }
0391 }
0392 
0393 CalibPlotProperties::~CalibPlotProperties() {
0394   delete corrFactor_;
0395   delete cFactor_;
0396   delete cSelect_;
0397   if (!fChain)
0398     return;
0399   delete fChain->GetCurrentFile();
0400 }
0401 
0402 Int_t CalibPlotProperties::GetEntry(Long64_t entry) {
0403   // Read contents of entry.
0404   if (!fChain)
0405     return 0;
0406   return fChain->GetEntry(entry);
0407 }
0408 
0409 Long64_t CalibPlotProperties::LoadTree(Long64_t entry) {
0410   // Set the environment to read one entry
0411   if (!fChain)
0412     return -5;
0413   Long64_t centry = fChain->LoadTree(entry);
0414   if (centry < 0)
0415     return centry;
0416   if (!fChain->InheritsFrom(TChain::Class()))
0417     return centry;
0418   TChain *chain = (TChain *)fChain;
0419   if (chain->GetTreeNumber() != fCurrent) {
0420     fCurrent = chain->GetTreeNumber();
0421     Notify();
0422   }
0423   return centry;
0424 }
0425 
0426 void CalibPlotProperties::Init(TChain *tree, const char *dupFileName) {
0427   // The Init() function is called when the selector needs to initialize
0428   // a new tree or chain. Typically here the branch addresses and branch
0429   // pointers of the tree will be set.
0430   // It is normally not necessary to make changes to the generated
0431   // code, but the routine can be extended by the user if needed.
0432   // Init() will be called many times when running on PROOF
0433   // (once per file to be processed).
0434 
0435   // Set object pointer
0436   t_DetIds = 0;
0437   t_DetIds1 = 0;
0438   t_DetIds3 = 0;
0439   t_HitEnergies = 0;
0440   t_HitEnergies1 = 0;
0441   t_HitEnergies3 = 0;
0442   t_trgbits = 0;
0443   // Set branch addresses and branch pointers
0444   fChain = tree;
0445   fCurrent = -1;
0446   if (!tree)
0447     return;
0448   fChain->SetMakeClass(1);
0449 
0450   fChain->SetBranchAddress("t_Run", &t_Run, &b_t_Run);
0451   fChain->SetBranchAddress("t_Event", &t_Event, &b_t_Event);
0452   fChain->SetBranchAddress("t_DataType", &t_DataType, &b_t_DataType);
0453   fChain->SetBranchAddress("t_ieta", &t_ieta, &b_t_ieta);
0454   fChain->SetBranchAddress("t_iphi", &t_iphi, &b_t_iphi);
0455   fChain->SetBranchAddress("t_EventWeight", &t_EventWeight, &b_t_EventWeight);
0456   fChain->SetBranchAddress("t_nVtx", &t_nVtx, &b_t_nVtx);
0457   fChain->SetBranchAddress("t_nTrk", &t_nTrk, &b_t_nTrk);
0458   fChain->SetBranchAddress("t_goodPV", &t_goodPV, &b_t_goodPV);
0459   fChain->SetBranchAddress("t_l1pt", &t_l1pt, &b_t_l1pt);
0460   fChain->SetBranchAddress("t_l1eta", &t_l1eta, &b_t_l1eta);
0461   fChain->SetBranchAddress("t_l1phi", &t_l1phi, &b_t_l1phi);
0462   fChain->SetBranchAddress("t_l3pt", &t_l3pt, &b_t_l3pt);
0463   fChain->SetBranchAddress("t_l3eta", &t_l3eta, &b_t_l3eta);
0464   fChain->SetBranchAddress("t_l3phi", &t_l3phi, &b_t_l3phi);
0465   fChain->SetBranchAddress("t_p", &t_p, &b_t_p);
0466   fChain->SetBranchAddress("t_pt", &t_pt, &b_t_pt);
0467   fChain->SetBranchAddress("t_phi", &t_phi, &b_t_phi);
0468   fChain->SetBranchAddress("t_mindR1", &t_mindR1, &b_t_mindR1);
0469   fChain->SetBranchAddress("t_mindR2", &t_mindR2, &b_t_mindR2);
0470   fChain->SetBranchAddress("t_eMipDR", &t_eMipDR, &b_t_eMipDR);
0471   fChain->SetBranchAddress("t_eHcal", &t_eHcal, &b_t_eHcal);
0472   fChain->SetBranchAddress("t_eHcal10", &t_eHcal10, &b_t_eHcal10);
0473   fChain->SetBranchAddress("t_eHcal30", &t_eHcal30, &b_t_eHcal30);
0474   fChain->SetBranchAddress("t_hmaxNearP", &t_hmaxNearP, &b_t_hmaxNearP);
0475   fChain->SetBranchAddress("t_rhoh", &t_rhoh, &b_t_rhoh);
0476   fChain->SetBranchAddress("t_selectTk", &t_selectTk, &b_t_selectTk);
0477   fChain->SetBranchAddress("t_qltyFlag", &t_qltyFlag, &b_t_qltyFlag);
0478   fChain->SetBranchAddress("t_qltyMissFlag", &t_qltyMissFlag, &b_t_qltyMissFlag);
0479   fChain->SetBranchAddress("t_qltyPVFlag", &t_qltyPVFlag, &b_t_qltyPVFlag);
0480   fChain->SetBranchAddress("t_gentrackP", &t_gentrackP, &b_t_gentrackP);
0481   fChain->SetBranchAddress("t_DetIds", &t_DetIds, &b_t_DetIds);
0482   fChain->SetBranchAddress("t_HitEnergies", &t_HitEnergies, &b_t_HitEnergies);
0483   fChain->SetBranchAddress("t_trgbits", &t_trgbits, &b_t_trgbits);
0484   fChain->SetBranchAddress("t_DetIds1", &t_DetIds1, &b_t_DetIds1);
0485   fChain->SetBranchAddress("t_DetIds3", &t_DetIds3, &b_t_DetIds3);
0486   fChain->SetBranchAddress("t_HitEnergies1", &t_HitEnergies1, &b_t_HitEnergies1);
0487   fChain->SetBranchAddress("t_HitEnergies3", &t_HitEnergies3, &b_t_HitEnergies3);
0488   Notify();
0489 
0490   if (std::string(dupFileName) != "") {
0491     ifstream infil1(dupFileName);
0492     if (!infil1.is_open()) {
0493       std::cout << "Cannot open duplicate file " << dupFileName << std::endl;
0494     } else {
0495       while (1) {
0496         Long64_t jentry;
0497         infil1 >> jentry;
0498         if (!infil1.good())
0499           break;
0500         entries_.push_back(jentry);
0501       }
0502       infil1.close();
0503       std::cout << "Reads a list of " << entries_.size() << " events from " << dupFileName << std::endl;
0504     }
0505   }
0506 
0507   char name[20], title[200];
0508   unsigned int kk(0);
0509 
0510   if (plotBasic_) {
0511     std::cout << "Book Basic Histos" << std::endl;
0512     for (int k = 0; k < CalibPlots::ntitles; ++k) {
0513       sprintf(name, "%sp%d", prefix_.c_str(), k);
0514       sprintf(title, "Momentum for %s", CalibPlots::getTitle(k).c_str());
0515       h_p[k] = new TH1D(name, title, 150, 0.0, 150.0);
0516       sprintf(name, "%seta%d", prefix_.c_str(), k);
0517       sprintf(title, "#eta for %s", CalibPlots::getTitle(k).c_str());
0518       h_eta[k] = new TH1D(name, title, 60, -30.0, 30.0);
0519     }
0520     for (int k = 0; k < CalibPlots::npbin; ++k) {
0521       sprintf(name, "%seta0%d", prefix_.c_str(), k);
0522       sprintf(title,
0523               "#eta for %s (p = %d:%d GeV)",
0524               CalibPlots::getTitle(0).c_str(),
0525               CalibPlots::getP(k),
0526               CalibPlots::getP(k + 1));
0527       h_eta0.push_back(new TH1D(name, title, 60, -30.0, 30.0));
0528       kk = h_eta0.size() - 1;
0529       h_eta0[kk]->Sumw2();
0530       sprintf(name, "%seta1%d", prefix_.c_str(), k);
0531       sprintf(title,
0532               "#eta for %s (p = %d:%d GeV)",
0533               CalibPlots::getTitle(1).c_str(),
0534               CalibPlots::getP(k),
0535               CalibPlots::getP(k + 1));
0536       h_eta1.push_back(new TH1D(name, title, 60, -30.0, 30.0));
0537       kk = h_eta1.size() - 1;
0538       h_eta1[kk]->Sumw2();
0539       sprintf(name, "%seta2%d", prefix_.c_str(), k);
0540       sprintf(title,
0541               "#eta for %s (p = %d:%d GeV)",
0542               CalibPlots::getTitle(2).c_str(),
0543               CalibPlots::getP(k),
0544               CalibPlots::getP(k + 1));
0545       h_eta2.push_back(new TH1D(name, title, 60, -30.0, 30.0));
0546       kk = h_eta2.size() - 1;
0547       h_eta2[kk]->Sumw2();
0548       sprintf(name, "%seta3%d", prefix_.c_str(), k);
0549       sprintf(title,
0550               "#eta for %s (p = %d:%d GeV)",
0551               CalibPlots::getTitle(3).c_str(),
0552               CalibPlots::getP(k),
0553               CalibPlots::getP(k + 1));
0554       h_eta3.push_back(new TH1D(name, title, 60, -30.0, 30.0));
0555       kk = h_eta3.size() - 1;
0556       h_eta3[kk]->Sumw2();
0557       sprintf(name, "%seta4%d", prefix_.c_str(), k);
0558       sprintf(title,
0559               "#eta for %s (p = %d:%d GeV)",
0560               CalibPlots::getTitle(4).c_str(),
0561               CalibPlots::getP(k),
0562               CalibPlots::getP(k + 1));
0563       h_eta4.push_back(new TH1D(name, title, 60, -30.0, 30.0));
0564       kk = h_eta4.size() - 1;
0565       h_eta4[kk]->Sumw2();
0566       sprintf(name, "%sdl1%d", prefix_.c_str(), k);
0567       sprintf(title, "Distance from L1 (p = %d:%d GeV)", CalibPlots::getP(k), CalibPlots::getP(k + 1));
0568       h_dL1.push_back(new TH1D(name, title, 160, 0.0, 8.0));
0569       kk = h_dL1.size() - 1;
0570       h_dL1[kk]->Sumw2();
0571       sprintf(name, "%svtx%d", prefix_.c_str(), k);
0572       sprintf(title, "N_{Vertex} (p = %d:%d GeV)", CalibPlots::getP(k), CalibPlots::getP(k + 1));
0573       h_vtx.push_back(new TH1D(name, title, 100, 0.0, 100.0));
0574       kk = h_vtx.size() - 1;
0575       h_vtx[kk]->Sumw2();
0576     }
0577   }
0578 
0579   if (plotEnergy_) {
0580     std::cout << "Make plots for good tracks" << std::endl;
0581     for (int k = 0; k < CalibPlots::npbin; ++k) {
0582       for (int j = etalo_; j <= etahi_ + 1; ++j) {
0583         sprintf(name, "%senergyH%d%d", prefix_.c_str(), k, j);
0584         if (j > etahi_)
0585           sprintf(title,
0586                   "HCAL energy for %s (p = %d:%d GeV |#eta| = %d:%d)",
0587                   CalibPlots::getTitle(3).c_str(),
0588                   CalibPlots::getP(k),
0589                   CalibPlots::getP(k + 1),
0590                   etalo_,
0591                   etahi_);
0592         else
0593           sprintf(title,
0594                   "HCAL energy for %s (p = %d:%d GeV |#eta| = %d)",
0595                   CalibPlots::getTitle(3).c_str(),
0596                   CalibPlots::getP(k),
0597                   CalibPlots::getP(k + 1),
0598                   j);
0599         h_etaEH[k].push_back(new TH1D(name, title, 200, 0, 200));
0600         kk = h_etaEH[k].size() - 1;
0601         h_etaEH[k][kk]->Sumw2();
0602         sprintf(name, "%senergyP%d%d", prefix_.c_str(), k, j);
0603         if (j > etahi_)
0604           sprintf(title,
0605                   "momentum for %s (p = %d:%d GeV |#eta| = %d:%d)",
0606                   CalibPlots::getTitle(3).c_str(),
0607                   CalibPlots::getP(k),
0608                   CalibPlots::getP(k + 1),
0609                   etalo_,
0610                   etahi_);
0611         else
0612           sprintf(title,
0613                   "momentum for %s (p = %d:%d GeV |#eta| = %d)",
0614                   CalibPlots::getTitle(3).c_str(),
0615                   CalibPlots::getP(k),
0616                   CalibPlots::getP(k + 1),
0617                   j);
0618         h_etaEp[k].push_back(new TH1D(name, title, 100, 0, 100));
0619         kk = h_etaEp[k].size() - 1;
0620         h_etaEp[k][kk]->Sumw2();
0621         sprintf(name, "%senergyE%d%d", prefix_.c_str(), k, j);
0622         if (j > etahi_)
0623           sprintf(title,
0624                   "ECAL energy for %s (p = %d:%d GeV |#eta| = %d:%d)",
0625                   CalibPlots::getTitle(3).c_str(),
0626                   CalibPlots::getP(k),
0627                   CalibPlots::getP(k + 1),
0628                   etalo_,
0629                   etahi_);
0630         else
0631           sprintf(title,
0632                   "ECAL energy for %s (p = %d:%d GeV |#eta| = %d)",
0633                   CalibPlots::getTitle(3).c_str(),
0634                   CalibPlots::getP(k),
0635                   CalibPlots::getP(k + 1),
0636                   j);
0637         h_etaEE[k].push_back(new TH1D(name, title, 100, 0, 10));
0638         kk = h_etaEE[k].size() - 1;
0639         h_etaEE[k][kk]->Sumw2();
0640         sprintf(name, "%senergyER%d%d", prefix_.c_str(), k, j);
0641         h_etaEE0[k].push_back(new TH1D(name, title, 100, 0, 1));
0642         kk = h_etaEE0[k].size() - 1;
0643         h_etaEE0[k][kk]->Sumw2();
0644       }
0645     }
0646 
0647     for (int j = 0; j < CalibPlots::netabin; ++j) {
0648       sprintf(name, "%senergyH%d", prefix_.c_str(), j);
0649       if (j == 0)
0650         sprintf(title, "HCAL energy for %s (All)", CalibPlots::getTitle(3).c_str());
0651       else
0652         sprintf(title,
0653                 "HCAL energy for %s (|#eta| = %d:%d)",
0654                 CalibPlots::getTitle(3).c_str(),
0655                 CalibPlots::getEta(j - 1),
0656                 CalibPlots::getEta(j));
0657       h_eHcal.push_back(new TH1D(name, title, 200, 0, 200));
0658       kk = h_eHcal.size() - 1;
0659       h_eHcal[kk]->Sumw2();
0660       sprintf(name, "%senergyP%d", prefix_.c_str(), j);
0661       if (j == 0)
0662         sprintf(title, "Track momentum for %s (All)", CalibPlots::getTitle(3).c_str());
0663       else
0664         sprintf(title,
0665                 "Track momentum for %s (|#eta| = %d:%d)",
0666                 CalibPlots::getTitle(3).c_str(),
0667                 CalibPlots::getEta(j - 1),
0668                 CalibPlots::getEta(j));
0669       h_mom.push_back(new TH1D(name, title, 100, 0, 100));
0670       kk = h_mom.size() - 1;
0671       h_mom[kk]->Sumw2();
0672       sprintf(name, "%senergyE%d", prefix_.c_str(), j);
0673       if (j == 0)
0674         sprintf(title, "ECAL energy for %s (All)", CalibPlots::getTitle(3).c_str());
0675       else
0676         sprintf(title,
0677                 "ECAL energy for %s (|#eta| = %d:%d)",
0678                 CalibPlots::getTitle(3).c_str(),
0679                 CalibPlots::getEta(j - 1),
0680                 CalibPlots::getEta(j));
0681       h_eEcal.push_back(new TH1D(name, title, 100, 0, 10));
0682       kk = h_eEcal.size() - 1;
0683       h_eEcal[kk]->Sumw2();
0684     }
0685   }
0686 
0687   if (plotHists_) {
0688     h_nvtx = new TH1D("hnvtx", "Number of vertices", 10, 0, 100);
0689     h_nvtx->Sumw2();
0690     for (int i = 0; i < CalibPlots::ndepth; i++) {
0691       sprintf(name, "b_edepth%d", i);
0692       sprintf(title, "Total RecHit energy in depth %d (Barrel)", i + 1);
0693       h_bvlist.push_back(new TH1F(name, title, 1000, 0, 100));
0694       h_bvlist[i]->Sumw2();
0695       sprintf(name, "b_recedepth%d", i);
0696       sprintf(title, "RecHit energy in depth %d (Barrel)", i + 1);
0697       h_bvlist2.push_back(new TH1F(name, title, 1000, 0, 100));
0698       h_bvlist2[i]->Sumw2();
0699       sprintf(name, "b_nrecdepth%d", i);
0700       sprintf(title, "#RecHits in depth %d (Barrel)", i + 1);
0701       h_bvlist3.push_back(new TH1F(name, title, 1000, 0, 100));
0702       h_bvlist3[i]->Sumw2();
0703       sprintf(name, "e_edepth%d", i);
0704       sprintf(title, "Total RecHit energy in depth %d (Endcap)", i + 1);
0705       h_evlist.push_back(new TH1F(name, title, 1000, 0, 100));
0706       h_evlist[i]->Sumw2();
0707       sprintf(name, "e_recedepth%d", i);
0708       sprintf(title, "RecHit energy in depth %d (Endcap)", i + 1);
0709       h_evlist2.push_back(new TH1F(name, title, 1000, 0, 100));
0710       h_evlist2[i]->Sumw2();
0711       sprintf(name, "e_nrecdepth%d", i);
0712       sprintf(title, "#RecHits in depth %d (Endcap)", i + 1);
0713       h_evlist3.push_back(new TH1F(name, title, 1000, 0, 100));
0714       h_evlist3[i]->Sumw2();
0715     }
0716     h_etaE = new TH2F("heta", "", 60, -30, 30, 100, 0, 100);
0717   }
0718 }
0719 
0720 Bool_t CalibPlotProperties::Notify() {
0721   // The Notify() function is called when a new file is opened. This
0722   // can be either for a new TTree in a TChain or when when a new TTree
0723   // is started when using PROOF. It is normally not necessary to make changes
0724   // to the generated code, but the routine can be extended by the
0725   // user if needed. The return value is currently not used.
0726 
0727   return kTRUE;
0728 }
0729 
0730 void CalibPlotProperties::Show(Long64_t entry) {
0731   // Print contents of entry.
0732   // If entry is not specified, print current entry
0733   if (!fChain)
0734     return;
0735   fChain->Show(entry);
0736 }
0737 
0738 Int_t CalibPlotProperties::Cut(Long64_t) {
0739   // This function may be called from Loop.
0740   // returns  1 if entry is accepted.
0741   // returns -1 otherwise.
0742   return 1;
0743 }
0744 
0745 void CalibPlotProperties::Loop(Long64_t nentries) {
0746   //   In a ROOT session, you can do:
0747   //      Root > .L CalibMonitor.C
0748   //      Root > CalibMonitor t
0749   //      Root > t.GetEntry(12); // Fill t data members with entry number 12
0750   //      Root > t.Show();       // Show values of entry 12
0751   //      Root > t.Show(16);     // Read and show values of entry 16
0752   //      Root > t.Loop();       // Loop on all entries
0753   //
0754 
0755   //   This is the loop skeleton where:
0756   //      jentry is the global entry number in the chain
0757   //      ientry is the entry number in the current Tree
0758   //   Note that the argument to GetEntry must be:
0759   //      jentry for TChain::GetEntry
0760   //      ientry for TTree::GetEntry and TBranch::GetEntry
0761   //
0762   //       To read only selected branches, Insert statements like:
0763   // METHOD1:
0764   //    fChain->SetBranchStatus("*",0);  // disable all branches
0765   //    fChain->SetBranchStatus("branchname",1);  // activate branchname
0766   // METHOD2: replace line
0767   //    fChain->GetEntry(jentry);       //read all branches
0768   //by  b_branchname->GetEntry(ientry); //read only this branch
0769   const bool debug(false);
0770 
0771   if (fChain == 0)
0772     return;
0773 
0774   // Find list of duplicate events
0775   if (nentries < 0)
0776     nentries = fChain->GetEntriesFast();
0777   std::cout << "Total entries " << nentries << ":" << fChain->GetEntriesFast() << std::endl;
0778   Long64_t nbytes(0), nb(0);
0779   unsigned int duplicate(0), good(0), kount(0);
0780   double sel(0), selHB(0), selHE(0);
0781   for (Long64_t jentry = 0; jentry < nentries; jentry++) {
0782     Long64_t ientry = LoadTree(jentry);
0783     if (ientry < 0)
0784       break;
0785     nb = fChain->GetEntry(jentry);
0786     nbytes += nb;
0787     if (jentry % 1000000 == 0)
0788       std::cout << "Entry " << jentry << " Run " << t_Run << " Event " << t_Event << std::endl;
0789     bool select = (std::find(entries_.begin(), entries_.end(), jentry) == entries_.end());
0790     if (!select) {
0791       ++duplicate;
0792       if (debug)
0793         std::cout << "Duplicate event " << t_Run << " " << t_Event << " " << t_p << std::endl;
0794       continue;
0795     }
0796     bool selRun = (includeRun_ ? ((t_Run >= runlo_) && (t_Run <= runhi_)) : ((t_Run < runlo_) || (t_Run > runhi_)));
0797     select =
0798         (selRun && (fabs(t_ieta) >= etalo_) && (fabs(t_ieta) <= etahi_) && (t_nVtx >= nvxlo_) && (t_nVtx <= nvxhi_));
0799     if (!select) {
0800       if (debug)
0801         std::cout << "Run # " << t_Run << " out of range of " << runlo_ << ":" << runhi_ << " or ieta " << t_ieta
0802                   << " (" << etalo_ << ":" << etahi_ << ") or nvtx " << t_nVtx << " (" << nvxlo_ << ":" << nvxhi_
0803                   << ") out of range" << std::endl;
0804       continue;
0805     }
0806     if (cSelect_ != nullptr) {
0807       if (exclude_) {
0808         if (cSelect_->isItRBX(t_DetIds))
0809           continue;
0810       } else {
0811         if (!(cSelect_->isItRBX(t_ieta, t_iphi)))
0812           continue;
0813       }
0814     }
0815     select = (!cutL1T_ || (t_mindR1 >= 0.5));
0816     if (!select) {
0817       if (debug)
0818         std::cout << "Reject Run # " << t_Run << " Event # " << t_Event << " too close to L1 trigger " << t_mindR1
0819                   << std::endl;
0820       continue;
0821     }
0822     select = ((events_.size() == 0) ||
0823               (std::find(events_.begin(), events_.end(), std::pair<int, int>(t_Run, t_Event)) != events_.end()));
0824     if (!select) {
0825       if (debug)
0826         std::cout << "Reject Run # " << t_Run << " Event # " << t_Event << " not in the selection list" << std::endl;
0827       continue;
0828     }
0829 
0830     // if (Cut(ientry) < 0) continue;
0831     int kp = CalibPlots::npbin0;
0832     double pmom = (useGen_ && (t_gentrackP > 0)) ? t_gentrackP : t_p;
0833     for (int k = 1; k < CalibPlots::npbin0; ++k) {
0834       if (pmom >= CalibPlots::getMomentum(k - 1) && pmom < CalibPlots::getMomentum(k)) {
0835         kp = k - 1;
0836         break;
0837       }
0838     }
0839     int jp1 = (((std::abs(t_ieta) >= etalo_) && (std::abs(t_ieta) <= etahi_)) ? (std::abs(t_ieta) - etalo_) : -1);
0840     int jp2 = (etahi_ - etalo_ + 1);
0841     int je1 = CalibPlots::netabin;
0842     for (int j = 1; j < CalibPlots::netabin; ++j) {
0843       if (std::abs(t_ieta) > CalibPlots::getEta(j - 1) && std::abs(t_ieta) <= CalibPlots::getEta(j)) {
0844         je1 = j;
0845         break;
0846       }
0847     }
0848     int je2 = (je1 != CalibPlots::netabin) ? 0 : -1;
0849     if (debug)
0850       std::cout << "Bin " << kp << ":" << je1 << ":" << je2 << ":" << jp1 << ":" << jp2 << std::endl;
0851     double cut = (pmom > 20) ? ((flexibleSelect_ == 0) ? 2.0 : 10.0) : 0.0;
0852     double rcut(-1000.0);
0853 
0854     // Selection of good track and energy measured in Hcal
0855     double eHcal(t_eHcal);
0856     if (corrFactor_->doCorr()) {
0857       eHcal = 0;
0858       for (unsigned int k = 0; k < t_HitEnergies->size(); ++k) {
0859         // The masks are defined in DataFormats/HcalDetId/interface/HcalDetId.h
0860         unsigned int id = truncateId((*t_DetIds)[k], truncateFlag_, false);
0861         double cfac = corrFactor_->getCorr(id);
0862         if (cFactor_ != 0)
0863           cfac *= cFactor_->getCorr(t_Run, (*t_DetIds)[k]);
0864         eHcal += (cfac * ((*t_HitEnergies)[k]));
0865         if (debug) {
0866           int subdet, zside, ieta, iphi, depth;
0867           unpackDetId(id, subdet, zside, ieta, iphi, depth);
0868           std::cout << zside << ":" << ieta << ":" << depth << " Corr " << cfac << " " << (*t_HitEnergies)[k] << " Out "
0869                     << eHcal << std::endl;
0870         }
0871       }
0872     }
0873     bool goodTk = goodTrack(eHcal, cut, debug);
0874     bool selPhi = selectPhi(debug);
0875     double rat = (pmom > 0) ? (eHcal / (pmom - t_eMipDR)) : 1.0;
0876     if (debug)
0877       std::cout << "Entry " << jentry << " p|eHcal|ratio " << pmom << "|" << t_eHcal << "|" << eHcal << "|" << rat
0878                 << "|" << kp << " Cuts " << t_qltyFlag << "|" << t_selectTk << "|" << (t_hmaxNearP < cut) << "|"
0879                 << (t_eMipDR < 1.0) << "|" << goodTk << "|" << (rat > rcut) << " Select Phi " << selPhi << std::endl;
0880     if (plotBasic_) {
0881       h_p[0]->Fill(pmom, t_EventWeight);
0882       h_eta[0]->Fill(t_ieta, t_EventWeight);
0883       if (kp < CalibPlots::npbin0)
0884         h_eta0[kp]->Fill(t_ieta, t_EventWeight);
0885       if (t_qltyFlag) {
0886         h_p[1]->Fill(pmom, t_EventWeight);
0887         h_eta[1]->Fill(t_ieta, t_EventWeight);
0888         if (kp < CalibPlots::npbin0)
0889           h_eta1[kp]->Fill(t_ieta, t_EventWeight);
0890         if (t_selectTk) {
0891           h_p[2]->Fill(pmom, t_EventWeight);
0892           h_eta[2]->Fill(t_ieta, t_EventWeight);
0893           if (kp < CalibPlots::npbin0)
0894             h_eta2[kp]->Fill(t_ieta, t_EventWeight);
0895           if (t_hmaxNearP < cut) {
0896             h_p[3]->Fill(pmom, t_EventWeight);
0897             h_eta[3]->Fill(t_ieta, t_EventWeight);
0898             if (kp < CalibPlots::npbin0)
0899               h_eta3[kp]->Fill(t_ieta, t_EventWeight);
0900             if (t_eMipDR < 1.0) {
0901               h_p[4]->Fill(pmom, t_EventWeight);
0902               h_eta[4]->Fill(t_ieta, t_EventWeight);
0903               if (kp < CalibPlots::npbin0) {
0904                 h_eta4[kp]->Fill(t_ieta, t_EventWeight);
0905                 h_dL1[kp]->Fill(t_mindR1, t_EventWeight);
0906                 h_vtx[kp]->Fill(t_goodPV, t_EventWeight);
0907               }
0908             }
0909           }
0910         }
0911       }
0912     }
0913 
0914     if (goodTk && kp != CalibPlots::npbin0 && selPhi) {
0915       if (rat > rcut) {
0916         if (plotEnergy_) {
0917           if (jp1 >= 0) {
0918             h_etaEH[kp][jp1]->Fill(eHcal, t_EventWeight);
0919             h_etaEH[kp][jp2]->Fill(eHcal, t_EventWeight);
0920             h_etaEp[kp][jp1]->Fill(pmom, t_EventWeight);
0921             h_etaEp[kp][jp2]->Fill(pmom, t_EventWeight);
0922             h_etaEE[kp][jp1]->Fill(t_eMipDR, t_EventWeight);
0923             h_etaEE[kp][jp2]->Fill(t_eMipDR, t_EventWeight);
0924             h_etaEE0[kp][jp1]->Fill(t_eMipDR, t_EventWeight);
0925             h_etaEE0[kp][jp2]->Fill(t_eMipDR, t_EventWeight);
0926           }
0927           if (kp == kp50) {
0928             if (je1 != CalibPlots::netabin) {
0929               h_eHcal[je1]->Fill(eHcal, t_EventWeight);
0930               h_eHcal[je2]->Fill(eHcal, t_EventWeight);
0931               h_mom[je1]->Fill(pmom, t_EventWeight);
0932               h_mom[je2]->Fill(pmom, t_EventWeight);
0933               h_eEcal[je1]->Fill(t_eMipDR, t_EventWeight);
0934               h_eEcal[je2]->Fill(t_eMipDR, t_EventWeight);
0935             }
0936           }
0937         }
0938 
0939         if (plotHists_) {
0940           h_nvtx->Fill(t_nVtx);
0941           if ((std::fabs(rat - 1) < 0.15) && (kp == kp50) && ((std::abs(t_ieta) < 15) || (std::abs(t_ieta) > 17))) {
0942             float weight = (dataMC_ ? t_EventWeight : t_EventWeight * puweight(t_nVtx));
0943             h_etaE->Fill(t_ieta, eHcal, weight);
0944             sel += weight;
0945             std::vector<float> bv(7, 0.0f), ev(7, 0.0f);
0946             std::vector<int> bnrec(7, 0), enrec(7, 0);
0947             double eb(0), ee(0);
0948             for (unsigned int k = 0; k < t_HitEnergies->size(); ++k) {
0949               unsigned int id = truncateId((*t_DetIds)[k], truncateFlag_, false);
0950               double cfac = corrFactor_->getCorr(id);
0951               if (cFactor_ != 0)
0952                 cfac *= cFactor_->getCorr(t_Run, (*t_DetIds)[k]);
0953               double ener = cfac * (*t_HitEnergies)[k];
0954               if (corrPU_)
0955                 correctEnergy(ener);
0956               unsigned int idx = (unsigned int)((*t_DetIds)[k]);
0957               int subdet, zside, ieta, iphi, depth;
0958               unpackDetId(idx, subdet, zside, ieta, iphi, depth);
0959               if (depth > 0 && depth <= CalibPlots::ndepth) {
0960                 if (subdet == 1) {
0961                   eb += ener;
0962                   bv[depth - 1] += ener;
0963                   h_bvlist2[depth - 1]->Fill(ener, weight);
0964                   ++bnrec[depth - 1];
0965                 } else if (subdet == 2) {
0966                   ee += ener;
0967                   ev[depth - 1] += ener;
0968                   h_evlist2[depth - 1]->Fill(ener, weight);
0969                   ++enrec[depth - 1];
0970                 }
0971               }
0972             }
0973             bool barrel = (eb > ee);
0974             if (barrel)
0975               selHB += weight;
0976             else
0977               selHE += weight;
0978             for (int i = 0; i < CalibPlots::ndepth; i++) {
0979               if (barrel) {
0980                 h_bvlist[i]->Fill(bv[i], weight);
0981                 h_bvlist3[i]->Fill((bnrec[i] + 0.001), weight);
0982               } else {
0983                 h_evlist[i]->Fill(ev[i], weight);
0984                 h_evlist3[i]->Fill((enrec[i] + 0.001), weight);
0985               }
0986             }
0987           }
0988         }
0989       }
0990       good++;
0991     }
0992     ++kount;
0993   }
0994   std::cout << "Finds " << duplicate << " Duplicate events out of " << kount << " events in this file and " << good
0995             << " selected events" << std::endl;
0996   if (plotHists_)
0997     std::cout << "Number of weighted selected events " << sel << " HB " << selHB << " HE " << selHE << std::endl;
0998 }
0999 
1000 bool CalibPlotProperties::goodTrack(double &eHcal, double &cuti, bool debug) {
1001   bool select(true);
1002   double cut(cuti);
1003   if (debug) {
1004     std::cout << "goodTrack input " << eHcal << ":" << cut;
1005   }
1006   if (flexibleSelect_ > 1) {
1007     double eta = (t_ieta > 0) ? t_ieta : -t_ieta;
1008     cut = 8.0 * exp(eta * log2by18_);
1009   }
1010   correctEnergy(eHcal);
1011   select = ((t_qltyFlag) && (t_selectTk) && (t_hmaxNearP < cut) && (t_eMipDR < 100.0));
1012   if (debug) {
1013     std::cout << " output " << eHcal << ":" << cut << ":" << select << std::endl;
1014   }
1015   return select;
1016 }
1017 
1018 bool CalibPlotProperties::selectPhi(bool debug) {
1019   bool select(true);
1020   if (phimin_ > 1 || phimax_ < 72) {
1021     double eTotal(0), eSelec(0);
1022     // The masks are defined in DataFormats/HcalDetId/interface/HcalDetId.h
1023     for (unsigned int k = 0; k < t_HitEnergies->size(); ++k) {
1024       int iphi = ((*t_DetIds)[k]) & (0x3FF);
1025       int zside = ((*t_DetIds)[k] & 0x80000) ? (1) : (-1);
1026       eTotal += ((*t_HitEnergies)[k]);
1027       if (iphi >= phimin_ && iphi <= phimax_ && zside == zside_)
1028         eSelec += ((*t_HitEnergies)[k]);
1029     }
1030     if (eSelec < 0.9 * eTotal)
1031       select = false;
1032     if (debug) {
1033       std::cout << "Etotal " << eTotal << " and ESelec " << eSelec << " (phi " << phimin_ << ":" << phimax_ << " z "
1034                 << zside_ << ") Selection " << select << std::endl;
1035     }
1036   }
1037   return select;
1038 }
1039 
1040 void CalibPlotProperties::savePlot(const std::string &theName, bool append, bool all, bool debug) {
1041   TFile *theFile(0);
1042   if (append) {
1043     theFile = new TFile(theName.c_str(), "UPDATE");
1044   } else {
1045     theFile = new TFile(theName.c_str(), "RECREATE");
1046   }
1047 
1048   theFile->cd();
1049 
1050   if (plotBasic_) {
1051     for (int k = 0; k < CalibPlots::ntitles; ++k) {
1052       if (debug)
1053         std::cout << "[" << k << "] p " << h_p[k] << " eta " << h_eta[k] << std::endl;
1054       h_p[k]->Write();
1055       h_eta[k]->Write();
1056     }
1057     for (int k = 0; k < CalibPlots::npbin; ++k) {
1058       if (debug)
1059         std::cout << "[" << k << "] eta " << h_eta0[k] << ":" << h_eta1[k] << ":" << h_eta2[k] << ":" << h_eta3[k]
1060                   << ":" << h_eta4[k] << " dl " << h_dL1[k] << " vtx " << h_vtx[k] << std::endl;
1061       if (h_eta0[k] != 0) {
1062         TH1D *h1 = (TH1D *)h_eta0[k]->Clone();
1063         h1->Write();
1064       }
1065       if (h_eta1[k] != 0) {
1066         TH1D *h2 = (TH1D *)h_eta1[k]->Clone();
1067         h2->Write();
1068       }
1069       if (h_eta2[k] != 0) {
1070         TH1D *h3 = (TH1D *)h_eta2[k]->Clone();
1071         h3->Write();
1072       }
1073       if (h_eta3[k] != 0) {
1074         TH1D *h4 = (TH1D *)h_eta3[k]->Clone();
1075         h4->Write();
1076       }
1077       if (h_eta4[k] != 0) {
1078         TH1D *h5 = (TH1D *)h_eta4[k]->Clone();
1079         h5->Write();
1080       }
1081       if (h_dL1[k] != 0) {
1082         TH1D *h6 = (TH1D *)h_dL1[k]->Clone();
1083         h6->Write();
1084       }
1085       if (h_vtx[k] != 0) {
1086         TH1D *h7 = (TH1D *)h_vtx[k]->Clone();
1087         h7->Write();
1088       }
1089     }
1090   }
1091 
1092   if (plotEnergy_) {
1093     for (int k = 0; k < CalibPlots::npbin; ++k) {
1094       if (debug)
1095         std::cout << "Energy[" << k << "] "
1096                   << " eta " << etalo_ << ":" << etahi_ << ":" << CalibPlots::netabin << " etaEH " << h_etaEH[k].size()
1097                   << " etaEp " << h_etaEp[k].size() << " etaEE " << h_etaEE[k].size() << " eHcal " << h_eHcal.size()
1098                   << " mom " << h_mom.size() << " eEcal " << h_eEcal.size() << std::endl;
1099       for (unsigned int j = 0; j <= (unsigned int)(etahi_ - etalo_); ++j) {
1100         if (h_etaEH[k].size() > j && h_etaEH[k][j] != nullptr && (all || (k == kp50))) {
1101           TH1D *hist = (TH1D *)h_etaEH[k][j]->Clone();
1102           hist->Write();
1103         }
1104         if (h_etaEp[k].size() > j && h_etaEp[k][j] != nullptr && (all || (k == kp50))) {
1105           TH1D *hist = (TH1D *)h_etaEp[k][j]->Clone();
1106           hist->Write();
1107         }
1108         if (h_etaEE[k].size() > j && h_etaEE[k][j] != nullptr && (all || (k == kp50))) {
1109           TH1D *hist = (TH1D *)h_etaEE[k][j]->Clone();
1110           hist->Write();
1111         }
1112         if (h_etaEE0[k].size() > j && h_etaEE0[k][j] != nullptr && (all || (k == kp50))) {
1113           TH1D *hist = (TH1D *)h_etaEE0[k][j]->Clone();
1114           hist->Write();
1115         }
1116       }
1117     }
1118 
1119     for (int j = 0; j < CalibPlots::netabin; ++j) {
1120       if (h_eHcal.size() > (unsigned int)(j) && (h_eHcal[j] != nullptr)) {
1121         TH1D *hist = (TH1D *)h_eHcal[j]->Clone();
1122         hist->Write();
1123       }
1124       if (h_mom.size() > (unsigned int)(j) && (h_mom[j] != nullptr)) {
1125         TH1D *hist = (TH1D *)h_mom[j]->Clone();
1126         hist->Write();
1127       }
1128       if (h_eEcal.size() > (unsigned int)(j) && (h_eEcal[j] != nullptr)) {
1129         TH1D *hist = (TH1D *)h_eEcal[j]->Clone();
1130         hist->Write();
1131       }
1132     }
1133   }
1134 
1135   if (plotHists_) {
1136     if (debug)
1137       std::cout << "nvtx " << h_nvtx << " etaE " << h_etaE << std::endl;
1138     h_nvtx->Write();
1139     h_etaE->Write();
1140     for (int i = 0; i < CalibPlots::ndepth; ++i) {
1141       if (debug)
1142         std::cout << "Depth[" << i << "] bv " << h_bvlist[i] << ":" << h_bvlist2[i] << ":" << h_bvlist3[i] << " ev "
1143                   << h_evlist[i] << ":" << h_evlist2[i] << ":" << h_evlist3[i] << std::endl;
1144       h_bvlist[i]->Write();
1145       h_bvlist2[i]->Write();
1146       h_bvlist3[i]->Write();
1147       h_evlist[i]->Write();
1148       h_evlist2[i]->Write();
1149       h_evlist3[i]->Write();
1150     }
1151   }
1152   std::cout << "All done" << std::endl;
1153   theFile->Close();
1154 }
1155 
1156 void CalibPlotProperties::correctEnergy(double &eHcal) {
1157   double pmom = (useGen_ && (t_gentrackP > 0)) ? t_gentrackP : t_p;
1158   if ((corrPU_ < 0) && (pmom > 0)) {
1159     double ediff = (t_eHcal30 - t_eHcal10);
1160     if (t_DetIds1 != 0 && t_DetIds3 != 0) {
1161       double Etot1(0), Etot3(0);
1162       // The masks are defined in DataFormats/HcalDetId/interface/HcalDetId.h
1163       for (unsigned int idet = 0; idet < (*t_DetIds1).size(); idet++) {
1164         unsigned int id = truncateId((*t_DetIds1)[idet], truncateFlag_, false);
1165         double cfac = corrFactor_->getCorr(id);
1166         if (cFactor_ != 0)
1167           cfac *= cFactor_->getCorr(t_Run, (*t_DetIds1)[idet]);
1168         double hitEn = cfac * (*t_HitEnergies1)[idet];
1169         Etot1 += hitEn;
1170       }
1171       for (unsigned int idet = 0; idet < (*t_DetIds3).size(); idet++) {
1172         unsigned int id = truncateId((*t_DetIds3)[idet], truncateFlag_, false);
1173         double cfac = corrFactor_->getCorr(id);
1174         if (cFactor_ != 0)
1175           cfac *= cFactor_->getCorr(t_Run, (*t_DetIds3)[idet]);
1176         double hitEn = cfac * (*t_HitEnergies3)[idet];
1177         Etot3 += hitEn;
1178       }
1179       ediff = (Etot3 - Etot1);
1180     }
1181     double fac = puFactor(-corrPU_, t_ieta, pmom, eHcal, ediff);
1182     eHcal *= fac;
1183   } else if (corrPU_ > 1) {
1184     eHcal = puFactorRho(corrPU_, t_ieta, t_rhoh, eHcal);
1185   }
1186 }
1187 
1188 void PlotThisHist(TH1D *hist, const std::string &text, int save) {
1189   char namep[120];
1190   sprintf(namep, "c_%s", hist->GetName());
1191   TCanvas *pad = new TCanvas(namep, namep, 700, 500);
1192   pad->SetRightMargin(0.10);
1193   pad->SetTopMargin(0.10);
1194   hist->GetXaxis()->SetTitleSize(0.04);
1195   hist->GetYaxis()->SetTitle("Tracks");
1196   hist->GetYaxis()->SetLabelOffset(0.005);
1197   hist->GetYaxis()->SetTitleSize(0.04);
1198   hist->GetYaxis()->SetLabelSize(0.035);
1199   hist->GetYaxis()->SetTitleOffset(1.10);
1200   hist->SetMarkerStyle(20);
1201   hist->SetMarkerColor(2);
1202   hist->SetLineColor(2);
1203   hist->Draw("Hist");
1204   pad->Modified();
1205   pad->Update();
1206   TPaveStats *st1 = (TPaveStats *)hist->GetListOfFunctions()->FindObject("stats");
1207   TPaveText *txt0 = new TPaveText(0.12, 0.91, 0.49, 0.96, "blNDC");
1208   txt0->SetFillColor(0);
1209   char txt[100];
1210   sprintf(txt, "CMS Simulation Preliminary");
1211   txt0->AddText(txt);
1212   txt0->Draw("same");
1213   TPaveText *txt1 = new TPaveText(0.51, 0.91, 0.90, 0.96, "blNDC");
1214   txt1->SetFillColor(0);
1215   sprintf(txt, "%s", text.c_str());
1216   txt1->AddText(txt);
1217   txt1->Draw("same");
1218   if (st1 != nullptr) {
1219     st1->SetY1NDC(0.70);
1220     st1->SetY2NDC(0.90);
1221     st1->SetX1NDC(0.65);
1222     st1->SetX2NDC(0.90);
1223   }
1224   pad->Update();
1225   if (save != 0) {
1226     if (save > 0)
1227       sprintf(namep, "%s.pdf", pad->GetName());
1228     else
1229       sprintf(namep, "%s.jpg", pad->GetName());
1230     pad->Print(namep);
1231   }
1232 }
1233 
1234 void PlotHist(const char *hisFileName,
1235               const std::string &prefix = "",
1236               const std::string &text = "",
1237               int flagC = 111,
1238               int etalo = 0,
1239               int etahi = 30,
1240               int save = 0) {
1241   gStyle->SetCanvasBorderMode(0);
1242   gStyle->SetCanvasColor(kWhite);
1243   gStyle->SetPadColor(kWhite);
1244   gStyle->SetFillColor(kWhite);
1245   gStyle->SetOptTitle(0);
1246   gStyle->SetOptStat(1110);
1247 
1248   bool plotBasic = (((flagC / 1) % 10) > 0);
1249   bool plotEnergy = (((flagC / 10) % 10) > 0);
1250   bool plotHists = (((flagC / 100) % 10) > 0);
1251 
1252   TFile *file = new TFile(hisFileName);
1253   char name[100], title[200];
1254   TH1D *hist;
1255   if ((file != nullptr) && plotBasic) {
1256     std::cout << "Plot Basic Histos" << std::endl;
1257     for (int k = 0; k < CalibPlots::ntitles; ++k) {
1258       sprintf(name, "%sp%d", prefix.c_str(), k);
1259       hist = (TH1D *)(file->FindObjectAny(name));
1260       if (hist != nullptr) {
1261         sprintf(title, "Momentum for %s (GeV)", CalibPlots::getTitle(k).c_str());
1262         hist->GetXaxis()->SetTitle(title);
1263         PlotThisHist(hist, text, save);
1264       }
1265       sprintf(name, "%seta%d", prefix.c_str(), k);
1266       hist = (TH1D *)(file->FindObjectAny(name));
1267       if (hist != nullptr) {
1268         sprintf(title, "#eta for %s", CalibPlots::getTitle(k).c_str());
1269         hist->GetXaxis()->SetTitle(title);
1270         PlotThisHist(hist, text, save);
1271       }
1272     }
1273     for (int k = 0; k < CalibPlots::npbin; ++k) {
1274       sprintf(name, "%seta0%d", prefix.c_str(), k);
1275       hist = (TH1D *)(file->FindObjectAny(name));
1276       if (hist != nullptr) {
1277         sprintf(title,
1278                 "#eta for %s (p = %d:%d GeV)",
1279                 CalibPlots::getTitle(0).c_str(),
1280                 CalibPlots::getP(k),
1281                 CalibPlots::getP(k + 1));
1282         hist->GetXaxis()->SetTitle(title);
1283         PlotThisHist(hist, text, save);
1284       }
1285       sprintf(name, "%seta1%d", prefix.c_str(), k);
1286       hist = (TH1D *)(file->FindObjectAny(name));
1287       if (hist != nullptr) {
1288         sprintf(title,
1289                 "#eta for %s (p = %d:%d GeV)",
1290                 CalibPlots::getTitle(1).c_str(),
1291                 CalibPlots::getP(k),
1292                 CalibPlots::getP(k + 1));
1293         hist->GetXaxis()->SetTitle(title);
1294         PlotThisHist(hist, text, save);
1295       }
1296       sprintf(name, "%seta2%d", prefix.c_str(), k);
1297       hist = (TH1D *)(file->FindObjectAny(name));
1298       if (hist != nullptr) {
1299         sprintf(title,
1300                 "#eta for %s (p = %d:%d GeV)",
1301                 CalibPlots::getTitle(2).c_str(),
1302                 CalibPlots::getP(k),
1303                 CalibPlots::getP(k + 1));
1304         hist->GetXaxis()->SetTitle(title);
1305         PlotThisHist(hist, text, save);
1306       }
1307       sprintf(name, "%seta3%d", prefix.c_str(), k);
1308       hist = (TH1D *)(file->FindObjectAny(name));
1309       if (hist != nullptr) {
1310         sprintf(title,
1311                 "#eta for %s (p = %d:%d GeV)",
1312                 CalibPlots::getTitle(3).c_str(),
1313                 CalibPlots::getP(k),
1314                 CalibPlots::getP(k + 1));
1315         hist->GetXaxis()->SetTitle(title);
1316         PlotThisHist(hist, text, save);
1317       }
1318       sprintf(name, "%seta4%d", prefix.c_str(), k);
1319       hist = (TH1D *)(file->FindObjectAny(name));
1320       if (hist != nullptr) {
1321         sprintf(title,
1322                 "#eta for %s (p = %d:%d GeV)",
1323                 CalibPlots::getTitle(4).c_str(),
1324                 CalibPlots::getP(k),
1325                 CalibPlots::getP(k + 1));
1326         hist->GetXaxis()->SetTitle(title);
1327         PlotThisHist(hist, text, save);
1328       }
1329       sprintf(name, "%sdl1%d", prefix.c_str(), k);
1330       hist = (TH1D *)(file->FindObjectAny(name));
1331       if (hist != nullptr) {
1332         sprintf(title, "Distance from L1 (p = %d:%d GeV)", CalibPlots::getP(k), CalibPlots::getP(k + 1));
1333         hist->GetXaxis()->SetTitle(title);
1334         PlotThisHist(hist, text, save);
1335       }
1336       sprintf(name, "%svtx%d", prefix.c_str(), k);
1337       hist = (TH1D *)(file->FindObjectAny(name));
1338       if (hist != nullptr) {
1339         sprintf(title, "N_{Vertex} (p = %d:%d GeV)", CalibPlots::getP(k), CalibPlots::getP(k + 1));
1340         hist->GetXaxis()->SetTitle(title);
1341         PlotThisHist(hist, text, save);
1342       }
1343     }
1344   }
1345 
1346   if ((file != nullptr) && plotEnergy) {
1347     std::cout << "Make plots for good tracks" << std::endl;
1348     for (int k = 0; k < CalibPlots::npbin; ++k) {
1349       for (int j = etalo; j <= etahi + 1; ++j) {
1350         sprintf(name, "%senergyH%d%d", prefix.c_str(), k, j);
1351         hist = (TH1D *)(file->FindObjectAny(name));
1352         if (hist != nullptr) {
1353           if (j > etahi)
1354             sprintf(title,
1355                     "HCAL energy for %s (p = %d:%d GeV |#eta| = %d:%d)",
1356                     CalibPlots::getTitle(3).c_str(),
1357                     CalibPlots::getP(k),
1358                     CalibPlots::getP(k + 1),
1359                     etalo,
1360                     etahi);
1361           else
1362             sprintf(title,
1363                     "HCAL energy for %s (p = %d:%d GeV |#eta| = %d)",
1364                     CalibPlots::getTitle(3).c_str(),
1365                     CalibPlots::getP(k),
1366                     CalibPlots::getP(k + 1),
1367                     j);
1368           hist->GetXaxis()->SetTitle(title);
1369           PlotThisHist(hist, text, save);
1370         }
1371         sprintf(name, "%senergyP%d%d", prefix.c_str(), k, j);
1372         hist = (TH1D *)(file->FindObjectAny(name));
1373         if (hist != nullptr) {
1374           if (j > etahi)
1375             sprintf(title,
1376                     "momentum for %s (p = %d:%d GeV |#eta| = %d:%d)",
1377                     CalibPlots::getTitle(3).c_str(),
1378                     CalibPlots::getP(k),
1379                     CalibPlots::getP(k + 1),
1380                     etalo,
1381                     etahi);
1382           else
1383             sprintf(title,
1384                     "momentum for %s (p = %d:%d GeV |#eta| = %d)",
1385                     CalibPlots::getTitle(3).c_str(),
1386                     CalibPlots::getP(k),
1387                     CalibPlots::getP(k + 1),
1388                     j);
1389           hist->GetXaxis()->SetTitle(title);
1390           PlotThisHist(hist, text, save);
1391         }
1392         sprintf(name, "%senergyE%d%d", prefix.c_str(), k, j);
1393         hist = (TH1D *)(file->FindObjectAny(name));
1394         if (hist != nullptr) {
1395           if (j > etahi)
1396             sprintf(title,
1397                     "ECAL energy for %s (p = %d:%d GeV |#eta| = %d:%d)",
1398                     CalibPlots::getTitle(3).c_str(),
1399                     CalibPlots::getP(k),
1400                     CalibPlots::getP(k + 1),
1401                     etalo,
1402                     etahi);
1403           else
1404             sprintf(title,
1405                     "ECAL energy for %s (p = %d:%d GeV |#eta| = %d)",
1406                     CalibPlots::getTitle(3).c_str(),
1407                     CalibPlots::getP(k),
1408                     CalibPlots::getP(k + 1),
1409                     j);
1410           hist->GetXaxis()->SetTitle(title);
1411           PlotThisHist(hist, text, save);
1412         }
1413         sprintf(name, "%senergyER%d%d", prefix.c_str(), k, j);
1414         hist = (TH1D *)(file->FindObjectAny(name));
1415         if (hist != nullptr) {
1416           std::cout << name << " Mean " << hist->GetMean() << " +- " << hist->GetMeanError() << " Entries "
1417                     << hist->GetEntries() << " RMS " << hist->GetRMS() << std::endl;
1418         }
1419       }
1420     }
1421 
1422     for (int j = 0; j < CalibPlots::netabin; ++j) {
1423       sprintf(name, "%senergyH%d", prefix.c_str(), j);
1424       hist = (TH1D *)(file->FindObjectAny(name));
1425       if (hist != nullptr) {
1426         if (j == 0)
1427           sprintf(title, "HCAL energy for %s (All)", CalibPlots::getTitle(3).c_str());
1428         else
1429           sprintf(title,
1430                   "HCAL energy for %s (|#eta| = %d:%d)",
1431                   CalibPlots::getTitle(3).c_str(),
1432                   CalibPlots::getEta(j - 1),
1433                   CalibPlots::getEta(j));
1434         hist->GetXaxis()->SetTitle(title);
1435         PlotThisHist(hist, text, save);
1436       }
1437       sprintf(name, "%senergyP%d", prefix.c_str(), j);
1438       hist = (TH1D *)(file->FindObjectAny(name));
1439       if (hist != nullptr) {
1440         if (j == 0)
1441           sprintf(title, "Track momentum for %s (All)", CalibPlots::getTitle(3).c_str());
1442         else
1443           sprintf(title,
1444                   "Track momentum for %s (|#eta| = %d:%d)",
1445                   CalibPlots::getTitle(3).c_str(),
1446                   CalibPlots::getEta(j - 1),
1447                   CalibPlots::getEta(j));
1448         hist->GetXaxis()->SetTitle(title);
1449         PlotThisHist(hist, text, save);
1450       }
1451       sprintf(name, "%senergyE%d", prefix.c_str(), j);
1452       hist = (TH1D *)(file->FindObjectAny(name));
1453       if (hist != nullptr) {
1454         if (j == 0)
1455           sprintf(title, "ECAL energy for %s (All)", CalibPlots::getTitle(3).c_str());
1456         else
1457           sprintf(title,
1458                   "ECAL energy for %s (|#eta| = %d:%d)",
1459                   CalibPlots::getTitle(3).c_str(),
1460                   CalibPlots::getEta(j - 1),
1461                   CalibPlots::getEta(j));
1462         hist->GetXaxis()->SetTitle(title);
1463         PlotThisHist(hist, text, save);
1464       }
1465     }
1466   }
1467 
1468   if (plotHists) {
1469     hist = (TH1D *)(file->FindObjectAny("hnvtx"));
1470     if (hist != nullptr) {
1471       hist->GetXaxis()->SetTitle("Number of vertices");
1472       PlotThisHist(hist, text, save);
1473     }
1474     for (int i = 0; i < CalibPlots::ndepth; i++) {
1475       sprintf(name, "b_edepth%d", i);
1476       hist = (TH1D *)(file->FindObjectAny(name));
1477       if (hist != nullptr) {
1478         sprintf(title, "Total RecHit energy in depth %d (Barrel)", i + 1);
1479         hist->GetXaxis()->SetTitle(title);
1480         PlotThisHist(hist, text, save);
1481       }
1482       sprintf(name, "b_recedepth%d", i);
1483       hist = (TH1D *)(file->FindObjectAny(name));
1484       if (hist != nullptr) {
1485         sprintf(title, "RecHit energy in depth %d (Barrel)", i + 1);
1486         hist->GetXaxis()->SetTitle(title);
1487         PlotThisHist(hist, text, save);
1488       }
1489       sprintf(name, "b_nrecdepth%d", i);
1490       hist = (TH1D *)(file->FindObjectAny(name));
1491       if (hist != nullptr) {
1492         sprintf(title, "#RecHits in depth %d (Barrel)", i + 1);
1493         hist->GetXaxis()->SetTitle(title);
1494         PlotThisHist(hist, text, save);
1495       }
1496       sprintf(name, "e_edepth%d", i);
1497       hist = (TH1D *)(file->FindObjectAny(name));
1498       if (hist != nullptr) {
1499         sprintf(title, "Total RecHit energy in depth %d (Endcap)", i + 1);
1500         hist->GetXaxis()->SetTitle(title);
1501         PlotThisHist(hist, text, save);
1502       }
1503       sprintf(name, "e_recedepth%d", i);
1504       hist = (TH1D *)(file->FindObjectAny(name));
1505       if (hist != nullptr) {
1506         sprintf(title, "RecHit energy in depth %d (Endcap)", i + 1);
1507         hist->GetXaxis()->SetTitle(title);
1508         PlotThisHist(hist, text, save);
1509       }
1510       sprintf(name, "e_nrecdepth%d", i);
1511       hist = (TH1D *)(file->FindObjectAny(name));
1512       if (hist != nullptr) {
1513         sprintf(title, "#RecHits in depth %d (Endcap)", i + 1);
1514         hist->GetXaxis()->SetTitle(title);
1515         PlotThisHist(hist, text, save);
1516       }
1517     }
1518     TH2F *h_etaE = (TH2F *)(file->FindObjectAny("heta"));
1519     if (h_etaE != nullptr) {
1520       h_etaE->GetXaxis()->SetTitle("i#eta");
1521       h_etaE->GetYaxis()->SetTitle("Energy (GeV)");
1522       char namep[120];
1523       sprintf(namep, "c_%s", h_etaE->GetName());
1524       TCanvas *pad = new TCanvas(namep, namep, 700, 700);
1525       pad->SetRightMargin(0.10);
1526       pad->SetTopMargin(0.10);
1527       h_etaE->GetXaxis()->SetTitleSize(0.04);
1528       h_etaE->GetYaxis()->SetLabelOffset(0.005);
1529       h_etaE->GetYaxis()->SetTitleSize(0.04);
1530       h_etaE->GetYaxis()->SetLabelSize(0.035);
1531       h_etaE->GetYaxis()->SetTitleOffset(1.10);
1532       h_etaE->SetMarkerStyle(20);
1533       h_etaE->SetMarkerColor(2);
1534       h_etaE->SetLineColor(2);
1535       h_etaE->Draw();
1536       pad->Update();
1537       TPaveStats *st1 = (TPaveStats *)h_etaE->GetListOfFunctions()->FindObject("stats");
1538       if (st1 != nullptr) {
1539         st1->SetY1NDC(0.70);
1540         st1->SetY2NDC(0.90);
1541         st1->SetX1NDC(0.65);
1542         st1->SetX2NDC(0.90);
1543       }
1544       if (save != 0) {
1545         if (save > 0)
1546           sprintf(namep, "%s.pdf", pad->GetName());
1547         else
1548           sprintf(namep, "%s.jpg", pad->GetName());
1549         pad->Print(namep);
1550       }
1551     }
1552   }
1553 }
[ Source navigation ] [ Diff markup ] [ Identifier search ] [ general search ]

This page was automatically generated by the 2.2.1 LXR engine.

The LXR team
Valid CSS 2.1! Valid HTML 4.01!