Back to home page

Project CMSSW displayed by LXR

 
 

    


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

0001 ////////////////////////////////////////////////////////////////////////////////

0002 // Usage:

0003 // .L CalibTree.C+g

0004 //  Run(inFileName, dirName, treeName, outFileName, corrFileName, dupFileName,

0005 //      rcorFileName, useIter, useweight, useMean, nMin, inverse, ratMin,

0006 //      ratMax, ietaMax, ietaTrack, sysmode, puCorr, applyL1Cut, l1Cut,

0007 //      truncateFlag, maxIter, corForm, useGen, runlo, runhi, phimin, phimax,

0008 //      zside, nvxlo, nvxhi, rbxFile, exclude, higheta, fraction,

0009 //      writeDebugHisto, pmin, pmax, debug, nmax);

0010 //

0011 //  where:

0012 //

0013 //  inFileName  (const char*) = file name of the input ROOT tree or name of

0014 //                              the file containing a list of file names of

0015 //                              input ROOT trees ("Silver.root")

0016 //  dirName     (const char*) = name of the directory where the Tree resides

0017 //                              ("HcalIsoTrkAnalyzer")

0018 //  treeName    (const char*) = name of the Tree ("CalibTree")

0019 //  outFileName (const char*) = name of the output ROOT file

0020 //                              ("Silver_out.root")

0021 //  corrFileName(const char*) = name of the output text file with correction

0022 //                              factors ("Silver_corr.txt")

0023 //  dupFileName (const char*) = name of the file containing list of sequence

0024 //                              numbers of duplicate entry ("events_DXS2.txt")

0025 //  rcorFileName (const char*)= name of the text file having the correction

0026 //                              factors as a function of run numbers or depth

0027 //                              to be used for raddam/depth/pileup/phisym

0028 //                              dependent correction  (default="", no corr.)

0029 //  useIter         (bool)    = Flag to use iterative process or minimization

0030 //                              (true)

0031 //  useweight       (bool)    = Flag to use event weight (true)

0032 //  useMean         (bool)    = Flag to use Mean of Most probable value

0033 //                              (false -- use MPV)

0034 //  nMin            (int)     = Minmum entries for a given cell which will be

0035 //                              used in evaluating convergence criterion (0)

0036 //  inverse         (bool)    = Use the ratio E/p or p/E in determining the

0037 //                              coefficients (true -- use p/E)

0038 //  ratMin          (double)  = Lower  cut on E/p to select a track (0.25)

0039 //  ratMax          (double)  = Higher cut on E/p to select a track (3.0)

0040 //  ietaMax         (int)     = Maximum ieta value for which correcttion

0041 //                              factor is to be determined (25)

0042 //  ietaTrack       (int)     = Maximum track ieta (at HCAL) to be considered

0043 //                              for these; -1 means no check on track ieta (-1)

0044 //  sysmode         (int)     = systematic error study (-1 if default)

0045 //                              -1 loose, -2 flexible, > 0 for systematic

0046 //  puCorr          (int)     = PU correction to be applied or not: 0 no

0047 //                              correction; < 0 use eDelta; > 0 rho dependent

0048 //                              correction (-8: 2022 Mahi version)

0049 //  applyL1Cut      (int)     = Flag to see if closeness to L1 object to be

0050 //                              applied: 0 no check; 1 only to events with

0051 //                              datatype not equal to 1; 2 to all (1)

0052 //  l1Cut           (double)  = Cut value for the closeness parameter (0.5)

0053 //  truncateFlag    (int)     = A two digit flag (dr) with the default value 0.

0054 //                              The digit *r* is used to treat depth values:

0055 //                              (0) treat each depth independently; (1) all

0056 //                              depths of ieta 15, 16 of HB as depth 1; (2)

0057 //                              all depths in HB and HE as depth 1; (3) all

0058 //                              depths in HE with values > 1 as depth 2; (4)

0059 //                              all depths in HB with values > 1 as depth 2;

0060 //                              (5) all depths in HB and HE with values > 1

0061 //                              as depth 2.

0062 //                              The digit *d* is used if zside is to be

0063 //                              ignored (1) or not (0)

0064 //                              (Default 0)

0065 //  maxIter         (int)     = number of iterations (30)

0066 //  drForm          (int)     = type of threshold/dupFileName/rcorFileName (hdr)

0067 //                              For rccorFileName r: (0) for Raddam correction,

0068 //                              (1) for depth dependent corrections; (2) for

0069 //                              RespCorr corrections; (3) use machine learning

0070 //                              method for pileup correction; (4) use results

0071 //                              from phi-symmetry.

0072 //                              For dupFileName d: (0) contains list of

0073 //                              duplicate entries; (1) depth dependent weights;

0074 //                              (2) list of  (ieta, iphi) of channels to be

0075 //                              selected.

0076 //                              For threshold h: the format for threshold

0077 //                              application, 0: no threshold; 1: 2022 prompt

0078 //                              data; 2: 2022 reco data; 3: 2023 prompt data.

0079 //                              (Default 0)

0080 //  useGen          (bool)    = use generator level momentum information (false)

0081 //  runlo           (int)     = lower value of run number to be included (+ve)

0082 //                              or excluded (-ve) (default 0)

0083 //  runhi           (int)     = higher value of run number to be included

0084 //                              (+ve) or excluded (-ve) (def 9999999)

0085 //  phimin          (int)     = minimum iphi value (1)

0086 //  phimax          (int)     = maximum iphi value (72)

0087 //  zside           (int)     = the side of the detector (0)

0088 //                              (if 0 no selection on zside will be made)

0089 //  nvxlo           (int)     = minimum # of vertex in event to be used (0)

0090 //  nvxhi           (int)     = maximum # of vertex in event to be used (1000)

0091 //  rbxFile         (char *)  = Name of the file containing a list of RBX's

0092 //                              to be consdered (default = ""). RBX's are

0093 //                              specified by zside*(Subdet*100+RBX #).

0094 //                              For HEP17 it will be 217

0095 //  exclude         (bool)    = RBX specified by the contents in *rbxFile* to

0096 //                              be exluded or only considered (false)

0097 //  higheta         (int)     = take correction factors for |ieta| > ietamax

0098 //                              as 1 (0) or of ieta = ietamax with same sign

0099 //                              and depth 1 (1) (default 1)

0100 //  fraction        (double)  = fraction of events to be done (1.0)

0101 //  writeDebugHisto (bool)    = Flag to check writing intermediate histograms

0102 //                              in o/p file (false)

0103 //  pmin            (double)  = minimum momentum of tracks to be used in

0104 //                              estimating the correction factor

0105 //  pmax            (double)  = maximum momentum of tracks to be used in

0106 //                              estimating the correction factor

0107 //  debug           (bool)    = To produce more debug printing on screen

0108 //                              (false)

0109 //  nmax            (Long64_t)= maximum number of entries to be processed,

0110 //                               if -1, all entries to be processed (-1)

0111 //

0112 //  doIt(inFileName, dupFileName)

0113 //  calls Run 5 times reducing # of events by a factor of 2 in each case

0114 ////////////////////////////////////////////////////////////////////////////////

0115 
0116 #include <TStyle.h>
0117 #include <TCanvas.h>
0118 #include <TROOT.h>
0119 #include <TChain.h>
0120 #include <TFile.h>
0121 #include <TFitResult.h>
0122 #include <TFitResultPtr.h>
0123 #include <TTree.h>
0124 #include <TH1.h>
0125 #include <TGraph.h>
0126 #include <TProfile.h>
0127 #include <algorithm>
0128 #include <vector>
0129 #include <string>
0130 #include <iomanip>
0131 #include <iostream>
0132 #include <fstream>
0133 #include <sstream>
0134 
0135 void unpackDetId(unsigned int, int &, int &, int &, int &, int &);
0136 #include "CalibCorr.C"
0137 
0138 void Run(const char *inFileName = "Silver.root",
0139          const char *dirName = "HcalIsoTrkAnalyzer",
0140          const char *treeName = "CalibTree",
0141          const char *outFileName = "Silver_out.root",
0142          const char *corrFileName = "Silver_corr.txt",
0143          const char *dupFileName = "events_DXS2.txt",
0144          const char *rcorFileName = "",
0145          bool useIter = true,
0146          bool useweight = true,
0147          bool useMean = false,
0148          int nMin = 0,
0149          bool inverse = true,
0150          double ratMin = 0.25,
0151          double ratMax = 3.,
0152          int ietaMax = 25,
0153          int ietaTrack = -1,
0154          int sysmode = -1,
0155          int puCorr = -8,
0156          int applyL1Cut = 1,
0157          double l1Cut = 0.5,
0158          int truncateFlag = 0,
0159          int maxIter = 30,
0160          int drForm = 0,
0161          bool useGen = false,
0162          int runlo = 0,
0163          int runhi = 99999999,
0164          int phimin = 1,
0165          int phimax = 72,
0166          int zside = 0,
0167          int nvxlo = 0,
0168          int nvxhi = 1000,
0169          const char *rbxFile = "",
0170          bool exclude = true,
0171          int higheta = 1,
0172          double fraction = 1.0,
0173          bool writeDebugHisto = false,
0174          double pmin = 40.0,
0175          double pmax = 60.0,
0176          bool debug = false,
0177          Long64_t nmax = -1);
0178 
0179 // Fixed size dimensions of array or collections stored in the TTree if any.

0180 
0181 class CalibTree {
0182 public:
0183   struct myEntry {
0184     myEntry(int k = 0, double f0 = 0, double f1 = 0, double f2 = 0) : kount(k), fact0(f0), fact1(f1), fact2(f2) {}
0185     int kount;
0186     double fact0, fact1, fact2;
0187   };
0188 
0189   struct energyCalor {
0190     energyCalor(double e1 = 0, double e2 = 0, double e3 = 0) : Etot(e1), Etot2(e2), ehcal(e3) {}
0191     double Etot, Etot2, ehcal;
0192   };
0193 
0194   CalibTree(const char *dupFileName,
0195             const char *rcorFileName,
0196             int truncateFlag,
0197             bool useIter,
0198             bool useMean,
0199             int runlo,
0200             int runhi,
0201             int phimin,
0202             int phimax,
0203             int zside,
0204             int nvxlo,
0205             int nvxhi,
0206             int sysmode,
0207             const char *rbxFile,
0208             int puCorr,
0209             int drForm,
0210             bool useGen,
0211             bool exclude,
0212             int higheta,
0213             double pmin,
0214             double pmax,
0215             TChain *tree);
0216   virtual ~CalibTree();
0217   virtual Int_t Cut(Long64_t entry);
0218   virtual Int_t GetEntry(Long64_t entry);
0219   virtual Long64_t LoadTree(Long64_t entry);
0220   virtual void Init(TChain *tree);
0221   virtual Double_t Loop(int k,
0222                         TFile *fout,
0223                         bool useweight,
0224                         int nMin,
0225                         bool inverse,
0226                         double rMin,
0227                         double rMax,
0228                         int ietaMax,
0229                         int ietaTrack,
0230                         int applyL1Cut,
0231                         double l1Cut,
0232                         bool last,
0233                         double fraction,
0234                         bool writeHisto,
0235                         bool debug,
0236                         Long64_t nmax);
0237   virtual Bool_t Notify();
0238   virtual void Show(Long64_t entry = -1);
0239   void getDetId(double fraction, int ietaTrack, bool debug, Long64_t nmax);
0240   void bookHistos(int loop, bool debug);
0241   bool goodTrack();
0242   void writeCorrFactor(const char *corrFileName, int ietaMax);
0243   bool selectPhi(unsigned int detId);
0244   std::pair<double, double> fitMean(TH1D *, int);
0245   void makeplots(double rmin, double rmax, int ietaMax, bool useWeight, double fraction, bool debug, Long64_t nmax);
0246   void fitPol0(TH1D *hist, bool debug);
0247   void highEtaFactors(int ietaMax, bool debug);
0248   energyCalor energyHcal(double pmom, const Long64_t &entry, bool final);
0249 
0250   TChain *fChain;  //!pointer to the analyzed TTree or TChain

0251   Int_t fCurrent;  //!current Tree number in a TChain

0252   TH1D *h_pbyE, *h_cvg;
0253   TProfile *h_Ebyp_bfr, *h_Ebyp_aftr;
0254 
0255 private:
0256   // Declaration of leaf types

0257   Int_t t_Run;
0258   Int_t t_Event;
0259   Int_t t_DataType;
0260   Int_t t_ieta;
0261   Int_t t_iphi;
0262   Double_t t_EventWeight;
0263   Int_t t_nVtx;
0264   Int_t t_nTrk;
0265   Int_t t_goodPV;
0266   Double_t t_l1pt;
0267   Double_t t_l1eta;
0268   Double_t t_l1phi;
0269   Double_t t_l3pt;
0270   Double_t t_l3eta;
0271   Double_t t_l3phi;
0272   Double_t t_p;
0273   Double_t t_pt;
0274   Double_t t_phi;
0275   Double_t t_mindR1;
0276   Double_t t_mindR2;
0277   Double_t t_eMipDR;
0278   Double_t t_eHcal;
0279   Double_t t_eHcal10;
0280   Double_t t_eHcal30;
0281   Double_t t_hmaxNearP;
0282   Double_t t_rhoh;
0283   Bool_t t_selectTk;
0284   Bool_t t_qltyFlag;
0285   Bool_t t_qltyMissFlag;
0286   Bool_t t_qltyPVFlag;
0287   Double_t t_gentrackP;
0288   std::vector<unsigned int> *t_DetIds;
0289   std::vector<double> *t_HitEnergies;
0290   std::vector<bool> *t_trgbits;
0291   std::vector<unsigned int> *t_DetIds1;
0292   std::vector<unsigned int> *t_DetIds3;
0293   std::vector<double> *t_HitEnergies1;
0294   std::vector<double> *t_HitEnergies3;
0295 
0296   // List of branches

0297   TBranch *b_t_Run;           //!

0298   TBranch *b_t_Event;         //!

0299   TBranch *b_t_DataType;      //!

0300   TBranch *b_t_ieta;          //!

0301   TBranch *b_t_iphi;          //!

0302   TBranch *b_t_EventWeight;   //!

0303   TBranch *b_t_nVtx;          //!

0304   TBranch *b_t_nTrk;          //!

0305   TBranch *b_t_goodPV;        //!

0306   TBranch *b_t_l1pt;          //!

0307   TBranch *b_t_l1eta;         //!

0308   TBranch *b_t_l1phi;         //!

0309   TBranch *b_t_l3pt;          //!

0310   TBranch *b_t_l3eta;         //!

0311   TBranch *b_t_l3phi;         //!

0312   TBranch *b_t_p;             //!

0313   TBranch *b_t_pt;            //!

0314   TBranch *b_t_phi;           //!

0315   TBranch *b_t_mindR1;        //!

0316   TBranch *b_t_mindR2;        //!

0317   TBranch *b_t_eMipDR;        //!

0318   TBranch *b_t_eHcal;         //!

0319   TBranch *b_t_eHcal10;       //!

0320   TBranch *b_t_eHcal30;       //!

0321   TBranch *b_t_hmaxNearP;     //!

0322   TBranch *b_t_rhoh;          //!

0323   TBranch *b_t_selectTk;      //!

0324   TBranch *b_t_qltyFlag;      //!

0325   TBranch *b_t_qltyMissFlag;  //!

0326   TBranch *b_t_qltyPVFlag;    //!

0327   TBranch *b_t_gentrackP;     //!

0328   TBranch *b_t_DetIds;        //!

0329   TBranch *b_t_HitEnergies;   //!

0330   TBranch *b_t_trgbits;       //!

0331   TBranch *b_t_DetIds1;       //!

0332   TBranch *b_t_DetIds3;       //!

0333   TBranch *b_t_HitEnergies1;  //!

0334   TBranch *b_t_HitEnergies3;  //!

0335 
0336   CalibCorr *cFactor_;
0337   CalibSelectRBX *cSelect_;
0338   CalibDuplicate *cDuplicate_;
0339   const int truncateFlag_;
0340   const bool useIter_;
0341   const bool useMean_;
0342   int runlo_, runhi_;
0343   const int phimin_, phimax_;
0344   const int zside_, nvxlo_, nvxhi_;
0345   const int sysmode_;
0346   const char *rbxFile_;
0347   const int puCorr_, drForm_;
0348   int rcorForm_, duplicate_;
0349   const bool useGen_, exclude_;
0350   const int higheta_;
0351   const double pmin_, pmax_;
0352   bool includeRun_;
0353   double log2by18_, eHcalDelta_;
0354   int thrForm_;
0355   std::vector<unsigned int> detIds_;
0356   std::map<unsigned int, TH1D *> histos_;
0357   std::map<unsigned int, std::pair<double, double> > Cprev;
0358 };
0359 
0360 void doIt(const char *infile, const char *dup) {
0361   char outf1[100], outf2[100];
0362   double lumt(1.0), fac(0.5);
0363   for (int k = 0; k < 5; ++k) {
0364     sprintf(outf1, "%s_%d.root", infile, k);
0365     sprintf(outf2, "%s_%d.txt", infile, k);
0366     double lumi = (k == 0) ? -1 : lumt;
0367     lumt *= fac;
0368     Run(infile,
0369         "HcalIsoTrkAnalyzer",
0370         "CalibTree",
0371         outf1,
0372         outf2,
0373         dup,
0374         "",
0375         true,
0376         true,
0377         false,
0378         0,
0379         true,
0380         0.25,
0381         3.0,
0382         25,
0383         -1,
0384         -1,
0385         -5,
0386         1,
0387         0.5,
0388         0,
0389         30,
0390         0,
0391         false,
0392         0,
0393         99999999,
0394         1,
0395         72,
0396         0,
0397         0,
0398         1000,
0399         "",
0400         true,
0401         1,
0402         lumi,
0403         false,
0404         40.0,
0405         60.0,
0406         false,
0407         -1);
0408   }
0409 }
0410 
0411 void Run(const char *inFileName,
0412          const char *dirName,
0413          const char *treeName,
0414          const char *outFileName,
0415          const char *corrFileName,
0416          const char *dupFileName,
0417          const char *rcorFileName,
0418          bool useIter,
0419          bool useweight,
0420          bool useMean,
0421          int nMin,
0422          bool inverse,
0423          double ratMin,
0424          double ratMax,
0425          int ietaMax,
0426          int ietaTrack,
0427          int sysmode,
0428          int puCorr,
0429          int applyL1Cut,
0430          double l1Cut,
0431          int truncateFlag,
0432          int maxIter,
0433          int drForm,
0434          bool useGen,
0435          int runlo,
0436          int runhi,
0437          int phimin,
0438          int phimax,
0439          int zside,
0440          int nvxlo,
0441          int nvxhi,
0442          const char *rbxFile,
0443          bool exclude,
0444          int higheta,
0445          double fraction,
0446          bool writeHisto,
0447          double pmin,
0448          double pmax,
0449          bool debug,
0450          Long64_t nmax) {
0451   char name[500];
0452   sprintf(name, "%s/%s", dirName, treeName);
0453   TChain *chain = new TChain(name);
0454   std::cout << "Create a chain for " << name << " from " << inFileName << std::endl;
0455   if (!fillChain(chain, inFileName)) {
0456     std::cout << "*****No valid tree chain can be obtained*****" << std::endl;
0457   } else {
0458     std::cout << "Proceed with a tree chain with " << chain->GetEntries() << " entries" << std::endl;
0459     Long64_t nentryTot = chain->GetEntries();
0460     Long64_t nentries = (fraction > 0.01 && fraction < 0.99) ? (Long64_t)(fraction * nentryTot) : nentryTot;
0461     if ((nentries > nmax) && (nmax > 0))
0462       nentries = nmax;
0463     static const int maxIterMax = 100;
0464     if (maxIter > maxIterMax)
0465       maxIter = maxIterMax;
0466     std::cout << "Tree " << name << " " << chain << " in directory " << dirName << " from file " << inFileName
0467               << " with nentries (tracks): " << nentries << std::endl;
0468     unsigned int k(0), kmax(maxIter);
0469     CalibTree t(dupFileName,
0470                 rcorFileName,
0471                 truncateFlag,
0472                 useIter,
0473                 useMean,
0474                 runlo,
0475                 runhi,
0476                 phimin,
0477                 phimax,
0478                 zside,
0479                 nvxlo,
0480                 nvxhi,
0481                 sysmode,
0482                 rbxFile,
0483                 puCorr,
0484                 drForm,
0485                 useGen,
0486                 exclude,
0487                 higheta,
0488                 pmin,
0489                 pmax,
0490                 chain);
0491     t.h_pbyE = new TH1D("pbyE", "pbyE", 100, -1.0, 9.0);
0492     t.h_Ebyp_bfr = new TProfile("Ebyp_bfr", "Ebyp_bfr", 60, -30, 30, 0, 10);
0493     t.h_Ebyp_aftr = new TProfile("Ebyp_aftr", "Ebyp_aftr", 60, -30, 30, 0, 10);
0494     t.h_cvg = new TH1D("Cvg0", "Convergence", kmax, 0, kmax);
0495     t.h_cvg->SetMarkerStyle(7);
0496     t.h_cvg->SetMarkerSize(5.0);
0497 
0498     TFile *fout = new TFile(outFileName, "RECREATE");
0499     std::cout << "Output file: " << outFileName << " opened in recreate mode" << std::endl;
0500     fout->cd();
0501 
0502     double cvgs[maxIterMax], itrs[maxIterMax];
0503     t.getDetId(fraction, ietaTrack, debug, nmax);
0504 
0505     for (; k <= kmax; ++k) {
0506       std::cout << "Calling Loop() " << k << "th time" << std::endl;
0507       double cvg = t.Loop(k,
0508                           fout,
0509                           useweight,
0510                           nMin,
0511                           inverse,
0512                           ratMin,
0513                           ratMax,
0514                           ietaMax,
0515                           ietaTrack,
0516                           applyL1Cut,
0517                           l1Cut,
0518                           k == kmax,
0519                           fraction,
0520                           writeHisto,
0521                           debug,
0522                           nmax);
0523       itrs[k] = k;
0524       cvgs[k] = cvg;
0525       if (cvg < 0.00001)
0526         break;
0527     }
0528 
0529     t.writeCorrFactor(corrFileName, ietaMax);
0530 
0531     fout->cd();
0532     TGraph *g_cvg;
0533     g_cvg = new TGraph(k, itrs, cvgs);
0534     g_cvg->SetMarkerStyle(7);
0535     g_cvg->SetMarkerSize(5.0);
0536     g_cvg->Draw("AP");
0537     g_cvg->Write("Cvg");
0538     std::cout << "Finish looping after " << k << " iterations" << std::endl;
0539     t.makeplots(ratMin, ratMax, ietaMax, useweight, fraction, debug, nmax);
0540     fout->Close();
0541   }
0542 }
0543 
0544 CalibTree::CalibTree(const char *dupFileName,
0545                      const char *rcorFileName,
0546                      int flag,
0547                      bool useIter,
0548                      bool useMean,
0549                      int runlo,
0550                      int runhi,
0551                      int phimin,
0552                      int phimax,
0553                      int zside,
0554                      int nvxlo,
0555                      int nvxhi,
0556                      int mode,
0557                      const char *rbxFile,
0558                      int pu,
0559                      int drForm,
0560                      bool gen,
0561                      bool excl,
0562                      int heta,
0563                      double pmin,
0564                      double pmax,
0565                      TChain *tree)
0566     : fChain(nullptr),
0567       cFactor_(nullptr),
0568       cSelect_(nullptr),
0569       cDuplicate_(nullptr),
0570       truncateFlag_(flag),
0571       useIter_(useIter),
0572       useMean_(useMean),
0573       runlo_(runlo),
0574       runhi_(runhi),
0575       phimin_(phimin),
0576       phimax_(phimax),
0577       zside_(zside),
0578       nvxlo_(nvxlo),
0579       nvxhi_(nvxhi),
0580       sysmode_(mode),
0581       rbxFile_(rbxFile),
0582       puCorr_(pu),
0583       drForm_(drForm),
0584       useGen_(gen),
0585       exclude_(excl),
0586       higheta_(heta),
0587       pmin_(pmin),
0588       pmax_(pmax),
0589       includeRun_(true) {
0590   if (runlo_ < 0 || runhi_ < 0) {
0591     runlo_ = std::abs(runlo_);
0592     runhi_ = std::abs(runhi_);
0593     includeRun_ = false;
0594   }
0595   log2by18_ = std::log(2.5) / 18.0;
0596   duplicate_ = (drForm_ / 10) % 10;
0597   rcorForm_ = (drForm_ % 10);
0598   thrForm_ = (drForm_ / 100);
0599   eHcalDelta_ = 0;
0600   std::cout << "Initialize CalibTree with TruncateFlag " << truncateFlag_ << " UseMean " << useMean_ << " Run Range "
0601             << runlo_ << ":" << runhi_ << " Phi Range " << phimin_ << ":" << phimax_ << ":" << zside_
0602             << " Vertex Range " << nvxlo_ << ":" << nvxhi_ << " Mode " << sysmode_ << " PU " << puCorr_ << " Gen "
0603             << useGen_ << " High Eta " << higheta_ << " Threshold Flag " << thrForm_ << std::endl;
0604   std::cout << "Duplicate events read from " << dupFileName << " duplicateFormat " << duplicate_
0605             << " RadDam Corrections read from " << rcorFileName << " rcorFormat " << rcorForm_ << " Treat RBX "
0606             << rbxFile_ << " with exclusion mode " << exclude_ << std::endl;
0607   Init(tree);
0608   if (std::string(dupFileName) != "")
0609     cDuplicate_ = new CalibDuplicate(dupFileName, duplicate_, false);
0610   if (std::string(rcorFileName) != "") {
0611     cFactor_ = new CalibCorr(rcorFileName, rcorForm_, false);
0612     if (cFactor_->absent())
0613       rcorForm_ = -1;
0614   } else {
0615     rcorForm_ = -1;
0616   }
0617   if (std::string(rbxFile) != "")
0618     cSelect_ = new CalibSelectRBX(rbxFile, false);
0619 }
0620 
0621 CalibTree::~CalibTree() {
0622   delete cFactor_;
0623   delete cSelect_;
0624   delete cDuplicate_;
0625   if (!fChain)
0626     return;
0627   delete fChain->GetCurrentFile();
0628 }
0629 
0630 Int_t CalibTree::GetEntry(Long64_t entry) {
0631   // Read contents of entry.

0632   if (!fChain)
0633     return 0;
0634   return fChain->GetEntry(entry);
0635 }
0636 
0637 Long64_t CalibTree::LoadTree(Long64_t entry) {
0638   // Set the environment to read one entry

0639   if (!fChain)
0640     return -5;
0641   Long64_t centry = fChain->LoadTree(entry);
0642   if (centry < 0)
0643     return centry;
0644   if (fChain->GetTreeNumber() != fCurrent) {
0645     fCurrent = fChain->GetTreeNumber();
0646     Notify();
0647   }
0648   return centry;
0649 }
0650 
0651 void CalibTree::Init(TChain *tree) {
0652   // The Init() function is called when the selector needs to initialize

0653   // a new tree or chain. Typically here the branch addresses and branch

0654   // pointers of the tree will be set.

0655   // It is normally not necessary to make changes to the generated

0656   // code, but the routine can be extended by the user if needed.

0657   // Init() will be called many times when running on PROOF

0658   // (once per file to be processed).

0659 
0660   // Set object pointer

0661   t_DetIds = 0;
0662   t_HitEnergies = 0;
0663   t_trgbits = 0;
0664   t_DetIds1 = 0;
0665   t_DetIds3 = 0;
0666   t_HitEnergies1 = 0;
0667   t_HitEnergies3 = 0;
0668   // Set branch addresses and branch pointers

0669   fChain = tree;
0670   if (!tree)
0671     return;
0672   fCurrent = -1;
0673   fChain->SetMakeClass(1);
0674 
0675   fChain->SetBranchAddress("t_Run", &t_Run, &b_t_Run);
0676   fChain->SetBranchAddress("t_Event", &t_Event, &b_t_Event);
0677   fChain->SetBranchAddress("t_DataType", &t_DataType, &b_t_DataType);
0678   fChain->SetBranchAddress("t_ieta", &t_ieta, &b_t_ieta);
0679   fChain->SetBranchAddress("t_iphi", &t_iphi, &b_t_iphi);
0680   fChain->SetBranchAddress("t_EventWeight", &t_EventWeight, &b_t_EventWeight);
0681   fChain->SetBranchAddress("t_nVtx", &t_nVtx, &b_t_nVtx);
0682   fChain->SetBranchAddress("t_nTrk", &t_nTrk, &b_t_nTrk);
0683   fChain->SetBranchAddress("t_goodPV", &t_goodPV, &b_t_goodPV);
0684   fChain->SetBranchAddress("t_l1pt", &t_l1pt, &b_t_l1pt);
0685   fChain->SetBranchAddress("t_l1eta", &t_l1eta, &b_t_l1eta);
0686   fChain->SetBranchAddress("t_l1phi", &t_l1phi, &b_t_l1phi);
0687   fChain->SetBranchAddress("t_l3pt", &t_l3pt, &b_t_l3pt);
0688   fChain->SetBranchAddress("t_l3eta", &t_l3eta, &b_t_l3eta);
0689   fChain->SetBranchAddress("t_l3phi", &t_l3phi, &b_t_l3phi);
0690   fChain->SetBranchAddress("t_p", &t_p, &b_t_p);
0691   fChain->SetBranchAddress("t_pt", &t_pt, &b_t_pt);
0692   fChain->SetBranchAddress("t_phi", &t_phi, &b_t_phi);
0693   fChain->SetBranchAddress("t_mindR1", &t_mindR1, &b_t_mindR1);
0694   fChain->SetBranchAddress("t_mindR2", &t_mindR2, &b_t_mindR2);
0695   fChain->SetBranchAddress("t_eMipDR", &t_eMipDR, &b_t_eMipDR);
0696   fChain->SetBranchAddress("t_eHcal", &t_eHcal, &b_t_eHcal);
0697   fChain->SetBranchAddress("t_eHcal10", &t_eHcal10, &b_t_eHcal10);
0698   fChain->SetBranchAddress("t_eHcal30", &t_eHcal30, &b_t_eHcal30);
0699   fChain->SetBranchAddress("t_hmaxNearP", &t_hmaxNearP, &b_t_hmaxNearP);
0700   fChain->SetBranchAddress("t_rhoh", &t_rhoh, &b_t_rhoh);
0701   fChain->SetBranchAddress("t_selectTk", &t_selectTk, &b_t_selectTk);
0702   fChain->SetBranchAddress("t_qltyFlag", &t_qltyFlag, &b_t_qltyFlag);
0703   fChain->SetBranchAddress("t_qltyMissFlag", &t_qltyMissFlag, &b_t_qltyMissFlag);
0704   fChain->SetBranchAddress("t_qltyPVFlag", &t_qltyPVFlag, &b_t_qltyPVFlag);
0705   fChain->SetBranchAddress("t_gentrackP", &t_gentrackP, &b_t_gentrackP);
0706   fChain->SetBranchAddress("t_DetIds", &t_DetIds, &b_t_DetIds);
0707   fChain->SetBranchAddress("t_HitEnergies", &t_HitEnergies, &b_t_HitEnergies);
0708   fChain->SetBranchAddress("t_trgbits", &t_trgbits, &b_t_trgbits);
0709   fChain->SetBranchAddress("t_DetIds1", &t_DetIds1, &b_t_DetIds1);
0710   fChain->SetBranchAddress("t_DetIds3", &t_DetIds3, &b_t_DetIds3);
0711   fChain->SetBranchAddress("t_HitEnergies1", &t_HitEnergies1, &b_t_HitEnergies1);
0712   fChain->SetBranchAddress("t_HitEnergies3", &t_HitEnergies3, &b_t_HitEnergies3);
0713   Notify();
0714 }
0715 
0716 Bool_t CalibTree::Notify() {
0717   // The Notify() function is called when a new file is opened. This

0718   // can be either for a new TTree in a TChain or when when a new TTree

0719   // is started when using PROOF. It is normally not necessary to make changes

0720   // to the generated code, but the routine can be extended by the

0721   // user if needed. The return value is currently not used.

0722 
0723   return kTRUE;
0724 }
0725 
0726 void CalibTree::Show(Long64_t entry) {
0727   // Print contents of entry.

0728   // If entry is not specified, print current entry

0729   if (!fChain)
0730     return;
0731   fChain->Show(entry);
0732 }
0733 
0734 Int_t CalibTree::Cut(Long64_t) {
0735   // This function may be called from Loop.

0736   // returns  1 if entry is accepted.

0737   // returns -1 otherwise.

0738   return 1;
0739 }
0740 
0741 Double_t CalibTree::Loop(int loop,
0742                          TFile *fout,
0743                          bool useweight,
0744                          int nMin,
0745                          bool inverse,
0746                          double rmin,
0747                          double rmax,
0748                          int ietaMax,
0749                          int ietaTrack,
0750                          int applyL1Cut,
0751                          double l1Cut,
0752                          bool last,
0753                          double fraction,
0754                          bool writeHisto,
0755                          bool debug,
0756                          Long64_t nmax) {
0757   Long64_t nbytes(0), nb(0);
0758   Long64_t nentryTot = fChain->GetEntriesFast();
0759   Long64_t nentries = (fraction > 0.01 && fraction < 0.99) ? (Long64_t)(fraction * nentryTot) : nentryTot;
0760   if ((nentries > nmax) && (nmax > 0))
0761     nentries = nmax;
0762 
0763   bookHistos(loop, debug);
0764   std::map<unsigned int, myEntry> SumW;
0765   std::map<unsigned int, double> nTrks;
0766 
0767   int ntkgood(0);
0768   for (Long64_t jentry = 0; jentry < nentries; jentry++) {
0769     Long64_t ientry = LoadTree(jentry);
0770     if (ientry < 0)
0771       break;
0772     nb = fChain->GetEntry(jentry);
0773     nbytes += nb;
0774     if (jentry % 1000000 == 0)
0775       std::cout << "Entry " << jentry << " Run " << t_Run << " Event " << t_Event << std::endl;
0776     bool select = ((cDuplicate_ != nullptr) && (duplicate_ == 0)) ? (cDuplicate_->isDuplicate(jentry)) : true;
0777     if (!select)
0778       continue;
0779     bool selRun = (includeRun_ ? ((t_Run >= runlo_) && (t_Run <= runhi_)) : ((t_Run < runlo_) || (t_Run > runhi_)));
0780     if (!selRun)
0781       continue;
0782     if ((t_nVtx < nvxlo_) || (t_nVtx > nvxhi_))
0783       continue;
0784     if (cSelect_ != nullptr) {
0785       if (exclude_) {
0786         if (cSelect_->isItRBX(t_DetIds))
0787           continue;
0788       } else {
0789         if (!(cSelect_->isItRBX(t_ieta, t_iphi)))
0790           continue;
0791       }
0792     }
0793     if (cDuplicate_ != nullptr) {
0794       if (cDuplicate_->select(t_ieta, t_iphi))
0795         continue;
0796     }
0797     bool selTrack = ((ietaTrack <= 0) || (abs(t_ieta) <= ietaTrack));
0798     if (!selTrack)
0799       continue;
0800     if ((rcorForm_ == 3) && (cFactor_ != nullptr) && (cFactor_->absent(ientry)))
0801       continue;
0802 
0803     if (debug) {
0804       std::cout << "***Entry (Track) Number : " << ientry << std::endl;
0805       std::cout << "p/eHCal/eMipDR/nDets : " << t_p << "/" << t_eHcal << "/" << t_eMipDR << "/" << (*t_DetIds).size()
0806                 << std::endl;
0807     }
0808     double pmom = (useGen_ && (t_gentrackP > 0)) ? t_gentrackP : t_p;
0809     if (goodTrack()) {
0810       ++ntkgood;
0811       CalibTree::energyCalor en = energyHcal(pmom, jentry, true);
0812       double evWt = (useweight) ? t_EventWeight : 1.0;
0813       if (en.ehcal > 0.001) {
0814         double pufac = (en.Etot > 0) ? (en.ehcal / en.Etot) : 1.0;
0815         double ratio = en.ehcal / (pmom - t_eMipDR);
0816         if (debug)
0817           std::cout << " Weights " << evWt << ":" << pufac << " Energy " << en.Etot2 << ":" << en.Etot << ":" << pmom
0818                     << ":" << t_eMipDR << ":" << t_eHcal << ":" << en.ehcal << " ratio " << ratio << std::endl;
0819         if (loop == 0) {
0820           h_pbyE->Fill(ratio, evWt);
0821           h_Ebyp_bfr->Fill(t_ieta, ratio, evWt);
0822         }
0823         if (last) {
0824           h_Ebyp_aftr->Fill(t_ieta, ratio, evWt);
0825         }
0826         bool l1c(true);
0827         if (applyL1Cut != 0)
0828           l1c = ((t_mindR1 >= l1Cut) || ((applyL1Cut == 1) && (t_DataType == 1)));
0829         if ((rmin >= 0 && ratio > rmin) && (rmax >= 0 && ratio < rmax) && l1c) {
0830           for (unsigned int idet = 0; idet < (*t_DetIds).size(); idet++) {
0831             // Apply thresholds if necessary

0832             bool okcell = (thrForm_ == 0) || ((*t_HitEnergies)[idet] > threshold((*t_DetIds)[idet], thrForm_));
0833             if (okcell && selectPhi((*t_DetIds)[idet])) {
0834               unsigned int id = (*t_DetIds)[idet];
0835               unsigned int detid = truncateId(id, truncateFlag_, false);
0836               double hitEn = 0.0;
0837               if (debug) {
0838                 std::cout << "idet " << idet << " detid/hitenergy : " << std::hex << (*t_DetIds)[idet] << ":" << detid
0839                           << "/" << (*t_HitEnergies)[idet] << std::endl;
0840               }
0841               if (Cprev.find(detid) != Cprev.end())
0842                 hitEn = Cprev[detid].first * (*t_HitEnergies)[idet];
0843               else
0844                 hitEn = (*t_HitEnergies)[idet];
0845               if ((rcorForm_ != 3) && (rcorForm_ >= 0) && (cFactor_))
0846                 hitEn *= cFactor_->getCorr(t_Run, id);
0847               if ((cDuplicate_ != nullptr) && (cDuplicate_->doCorr()))
0848                 hitEn *= cDuplicate_->getWeight(id);
0849               double Wi = evWt * hitEn / en.Etot;
0850               double Fac = (inverse) ? (en.ehcal / (pmom - t_eMipDR)) : ((pmom - t_eMipDR) / en.ehcal);
0851               double Fac2 = Wi * Fac * Fac;
0852               TH1D *hist(0);
0853               std::map<unsigned int, TH1D *>::const_iterator itr = histos_.find(detid);
0854               if (itr != histos_.end())
0855                 hist = itr->second;
0856               if (debug) {
0857                 std::cout << "Det Id " << std::hex << detid << std::dec << " " << hist << std::endl;
0858               }
0859               if (hist != 0)
0860                 hist->Fill(Fac, Wi);  //////histola

0861               Fac *= Wi;
0862               if (SumW.find(detid) != SumW.end()) {
0863                 Wi += SumW[detid].fact0;
0864                 Fac += SumW[detid].fact1;
0865                 Fac2 += SumW[detid].fact2;
0866                 int kount = SumW[detid].kount + 1;
0867                 SumW[detid] = myEntry(kount, Wi, Fac, Fac2);
0868                 nTrks[detid] += evWt;
0869               } else {
0870                 SumW.insert(std::pair<unsigned int, myEntry>(detid, myEntry(1, Wi, Fac, Fac2)));
0871                 nTrks.insert(std::pair<unsigned int, unsigned int>(detid, evWt));
0872               }
0873             }
0874           }
0875         }
0876       }
0877     }
0878   }
0879   if (debug)
0880     std::cout << "# of Good Tracks " << ntkgood << " out of " << nentries << std::endl;
0881   if (loop == 0) {
0882     h_pbyE->Write("h_pbyE");
0883     h_Ebyp_bfr->Write("h_Ebyp_bfr");
0884   }
0885   if (last) {
0886     h_Ebyp_aftr->Write("h_Ebyp_aftr");
0887   }
0888 
0889   std::map<unsigned int, std::pair<double, double> > cfactors;
0890   unsigned int kount(0), kountus(0);
0891   double sumfactor(0);
0892   for (std::map<unsigned int, TH1D *>::const_iterator itr = histos_.begin(); itr != histos_.end(); ++itr) {
0893     if (writeHisto) {
0894       //      std::pair<double, double> result_write = fitMean(itr->second, 0);

0895       (itr->second)->Write();
0896     }
0897     // The masks are defined in DataFormats/HcalDetId/interface/HcalDetId.h

0898     int subdet, depth, zside, ieta, iphi;
0899     unpackDetId(itr->first, subdet, zside, ieta, iphi, depth);
0900     if (debug) {
0901       std::cout << "DETID :" << subdet << "  IETA :" << ieta << " HIST ENTRIES :" << (itr->second)->GetEntries()
0902                 << std::endl;
0903     }
0904   }
0905 
0906   if (debug)
0907     std::cout << "Histos with " << histos_.size() << " entries\n";
0908   for (std::map<unsigned int, TH1D *>::const_iterator itr = histos_.begin(); itr != histos_.end(); ++itr, ++kount) {
0909     std::pair<double, double> result = fitMean(itr->second, 0);
0910     double factor = (inverse) ? (2. - result.first) : result.first;
0911     if (debug) {
0912       int subdet, depth, zside, ieta, iphi;
0913       unpackDetId(itr->first, subdet, zside, ieta, iphi, depth);
0914       std::cout << "DetId[" << kount << "] " << subdet << ":" << zside * ieta << ":" << depth << " Factor " << factor
0915                 << " +- " << result.second << std::endl;
0916     }
0917     if (!useMean_) {
0918       cfactors[itr->first] = std::pair<double, double>(factor, result.second);
0919       if (itr->second->GetEntries() > nMin) {
0920         kountus++;
0921         if (factor > 1)
0922           sumfactor += (1 - 1 / factor);
0923         else
0924           sumfactor += (1 - factor);
0925       }
0926     }
0927   }
0928 
0929   if (debug)
0930     std::cout << "SumW with " << SumW.size() << " entries\n";
0931   std::map<unsigned int, myEntry>::const_iterator SumWItr = SumW.begin();
0932   for (; SumWItr != SumW.end(); SumWItr++) {
0933     unsigned int detid = SumWItr->first;
0934     int subdet, depth, zside, ieta, iphi;
0935     unpackDetId(detid, subdet, zside, ieta, iphi, depth);
0936     if (debug) {
0937       std::cout << "Detid|kount|SumWi|SumFac|myId : " << subdet << ":" << zside * ieta << ":" << depth << " | "
0938                 << (SumWItr->second).kount << " | " << (SumWItr->second).fact0 << "|" << (SumWItr->second).fact1 << "|"
0939                 << (SumWItr->second).fact2 << std::endl;
0940     }
0941     double factor = (SumWItr->second).fact1 / (SumWItr->second).fact0;
0942     double dfac1 = ((SumWItr->second).fact2 / (SumWItr->second).fact0 - factor * factor);
0943     if (dfac1 < 0)
0944       dfac1 = 0;
0945     double dfac = sqrt(dfac1 / (SumWItr->second).kount);
0946     if (debug) {
0947       std::cout << "Factor " << factor << " " << dfac1 << " " << dfac << std::endl;
0948     }
0949     if (inverse)
0950       factor = 2. - factor;
0951     if (useMean_) {
0952       cfactors[detid] = std::pair<double, double>(factor, dfac);
0953       if ((SumWItr->second).kount > nMin) {
0954         kountus++;
0955         if (factor > 1)
0956           sumfactor += (1 - 1 / factor);
0957         else
0958           sumfactor += (1 - factor);
0959       }
0960     }
0961   }
0962 
0963   static const int maxch = 500;
0964   double dets[maxch], cfacs[maxch], wfacs[maxch], myId[maxch], nTrk[maxch];
0965   std::cout << "cafctors: " << cfactors.size() << ":" << maxch << std::endl;
0966   kount = 0;
0967   std::map<unsigned int, std::pair<double, double> >::const_iterator itr = cfactors.begin();
0968   const double factorMin(0.1);
0969   for (; itr != cfactors.end(); ++itr, ++kount) {
0970     unsigned int detid = itr->first;
0971     int subdet, depth, zside, ieta, iphi;
0972     unpackDetId(detid, subdet, zside, ieta, iphi, depth);
0973     double id = ieta * zside + 0.25 * (depth - 1);
0974     double factor = (itr->second).first;
0975     double dfac = (itr->second).second;
0976     if ((ieta > ietaMax) || (factor < factorMin)) {
0977       factor = 1;
0978       dfac = 0;
0979     }
0980     std::pair<double, double> cfac(factor, dfac);
0981     if (Cprev.find(detid) != Cprev.end()) {
0982       dfac /= factor;
0983       factor *= Cprev[detid].first;
0984       dfac *= factor;
0985       Cprev[detid] = std::pair<double, double>(factor, dfac);
0986       cfacs[kount] = factor;
0987     } else {
0988       Cprev[detid] = std::pair<double, double>(factor, dfac);
0989       cfacs[kount] = factor;
0990     }
0991     wfacs[kount] = factor;
0992     dets[kount] = detid;
0993     myId[kount] = id;
0994     nTrk[kount] = nTrks[detid];
0995   }
0996   if (higheta_ > 0)
0997     highEtaFactors(ietaMax, debug);
0998 
0999   std::cout << kountus << " detids out of " << kount << " have tracks > " << nMin << std::endl;
1000 
1001   char fname[50];
1002   fout->cd();
1003   TGraph *g_fac1 = new TGraph(kount, dets, cfacs);
1004   sprintf(fname, "Cfacs%d", loop);
1005   g_fac1->SetMarkerStyle(7);
1006   g_fac1->SetMarkerSize(5.0);
1007   g_fac1->Draw("AP");
1008   g_fac1->Write(fname);
1009   TGraph *g_fac2 = new TGraph(kount, dets, wfacs);
1010   sprintf(fname, "Wfacs%d", loop);
1011   g_fac2->SetMarkerStyle(7);
1012   g_fac2->SetMarkerSize(5.0);
1013   g_fac2->Draw("AP");
1014   g_fac2->Write(fname);
1015   TGraph *g_fac3 = new TGraph(kount, myId, cfacs);
1016   sprintf(fname, "CfacsVsMyId%d", loop);
1017   g_fac3->SetMarkerStyle(7);
1018   g_fac3->SetMarkerSize(5.0);
1019   g_fac3->Draw("AP");
1020   g_fac3->Write(fname);
1021   TGraph *g_fac4 = new TGraph(kount, myId, wfacs);
1022   sprintf(fname, "WfacsVsMyId%d", loop);
1023   g_fac4->SetMarkerStyle(7);
1024   g_fac4->SetMarkerSize(5.0);
1025   g_fac4->Draw("AP");
1026   g_fac4->Write(fname);
1027   TGraph *g_nTrk = new TGraph(kount, myId, nTrk);
1028   sprintf(fname, "nTrk");
1029   if (loop == 0) {
1030     g_nTrk->SetMarkerStyle(7);
1031     g_nTrk->SetMarkerSize(5.0);
1032     g_nTrk->Draw("AP");
1033     g_nTrk->Write(fname);
1034   }
1035   std::cout << "The new factors are :" << std::endl;
1036   std::map<unsigned int, std::pair<double, double> >::const_iterator CprevItr = Cprev.begin();
1037   unsigned int indx(0);
1038   for (; CprevItr != Cprev.end(); CprevItr++, indx++) {
1039     unsigned int detid = CprevItr->first;
1040     int subdet, depth, zside, ieta, iphi;
1041     unpackDetId(detid, subdet, zside, ieta, iphi, depth);
1042     std::cout << "DetId[" << indx << "] " << std::hex << detid << std::dec << "(" << ieta * zside << "," << depth
1043               << ") (nTrks:" << nTrks[detid] << ") : " << CprevItr->second.first << " +- " << CprevItr->second.second
1044               << std::endl;
1045   }
1046   double mean = (kountus > 0) ? (sumfactor / kountus) : 0;
1047   std::cout << "Mean deviation " << mean << " from 1 for " << kountus << " DetIds" << std::endl;
1048   h_cvg->SetBinContent(loop + 1, mean);
1049   if (last)
1050     h_cvg->Write("Cvg0");
1051   return mean;
1052 }
1053 
1054 void CalibTree::getDetId(double fraction, int ietaTrack, bool debug, Long64_t nmax) {
1055   if (fChain != 0) {
1056     Long64_t nbytes(0), nb(0), kprint(0);
1057     Long64_t nentryTot = fChain->GetEntriesFast();
1058     Long64_t nentries = (fraction > 0.01 && fraction < 0.99) ? (Long64_t)(fraction * nentryTot) : nentryTot;
1059     if ((nentries > nmax) && (nmax > 0))
1060       nentries = nmax;
1061 
1062     for (Long64_t jentry = 0; jentry < nentries; jentry++) {
1063       Long64_t ientry = LoadTree(jentry);
1064       if (ientry < 0)
1065         break;
1066       nb = fChain->GetEntry(jentry);
1067       nbytes += nb;
1068       if (jentry % 1000000 == 0)
1069         std::cout << "Entry " << jentry << " Run " << t_Run << " Event " << t_Event << std::endl;
1070       bool select = ((cDuplicate_ != nullptr) && (duplicate_ == 0)) ? (cDuplicate_->isDuplicate(jentry)) : true;
1071       if (!select)
1072         continue;
1073       // Find DetIds contributing to the track

1074       bool selRun = (includeRun_ ? ((t_Run >= runlo_) && (t_Run <= runhi_)) : ((t_Run < runlo_) || (t_Run > runhi_)));
1075       bool selTrack = ((ietaTrack <= 0) || (abs(t_ieta) <= ietaTrack));
1076       if (selRun && (t_nVtx >= nvxlo_) && (t_nVtx <= nvxhi_) && selTrack) {
1077         bool isItRBX(false);
1078         if (cSelect_ != nullptr) {
1079           bool temp = cSelect_->isItRBX(t_DetIds);
1080           if (exclude_)
1081             isItRBX = temp;
1082           else
1083             isItRBX = !(temp);
1084         }
1085         if ((cDuplicate_ != nullptr) && (!isItRBX))
1086           isItRBX = (cDuplicate_->select(t_ieta, t_iphi));
1087         ++kprint;
1088         if (!(isItRBX)) {
1089           for (unsigned int idet = 0; idet < (*t_DetIds).size(); idet++) {
1090             if (selectPhi((*t_DetIds)[idet])) {
1091               unsigned int detid = truncateId((*t_DetIds)[idet], truncateFlag_, debug);
1092               if (debug && (kprint <= 10)) {
1093                 std::cout << "DetId[" << idet << "] Original " << std::hex << (*t_DetIds)[idet] << " truncated "
1094                           << detid << std::dec;
1095               }
1096               if (std::find(detIds_.begin(), detIds_.end(), detid) == detIds_.end()) {
1097                 detIds_.push_back(detid);
1098                 if (debug && (kprint <= 10))
1099                   std::cout << " new";
1100               }
1101               if (debug && (kprint <= 10))
1102                 std::cout << std::endl;
1103             }
1104           }
1105           // Also look at the neighbouring cells if available

1106           if (t_DetIds3 != 0) {
1107             for (unsigned int idet = 0; idet < (*t_DetIds3).size(); idet++) {
1108               if (selectPhi((*t_DetIds3)[idet])) {
1109                 unsigned int detid = truncateId((*t_DetIds3)[idet], truncateFlag_, debug);
1110                 if (std::find(detIds_.begin(), detIds_.end(), detid) == detIds_.end()) {
1111                   detIds_.push_back(detid);
1112                 }
1113               }
1114             }
1115           }
1116         }
1117       }
1118     }
1119   }
1120   if (debug) {
1121     std::cout << "Total of " << detIds_.size() << " detIds" << std::endl;
1122     // The masks are defined in DataFormats/HcalDetId/interface/HcalDetId.h

1123     for (unsigned int k = 0; k < detIds_.size(); ++k) {
1124       int subdet, depth, zside, ieta, iphi;
1125       unpackDetId(detIds_[k], subdet, zside, ieta, iphi, depth);
1126       std::cout << "DetId[" << k << "] " << subdet << ":" << zside * ieta << ":" << depth << ":" << iphi << "  "
1127                 << std::hex << detIds_[k] << std::dec << std::endl;
1128     }
1129   }
1130 }
1131 
1132 void CalibTree::bookHistos(int loop, bool debug) {
1133   unsigned int k(0);
1134   for (std::map<unsigned int, TH1D *>::const_iterator itr = histos_.begin(); itr != histos_.end(); ++itr, ++k) {
1135     if (debug) {
1136       std::cout << "histos[" << k << "] " << std::hex << itr->first << std::dec << " " << itr->second;
1137       if (itr->second != 0)
1138         std::cout << " " << itr->second->GetTitle();
1139       std::cout << std::endl;
1140     }
1141     if (itr->second != 0)
1142       itr->second->Delete();
1143   }
1144 
1145   for (unsigned int k = 0; k < detIds_.size(); ++k) {
1146     char name[20], title[100];
1147     sprintf(name, "Hist%d_%d", detIds_[k], loop);
1148     int subdet, depth, zside, ieta, iphi;
1149     unpackDetId(detIds_[k], subdet, zside, ieta, iphi, depth);
1150     sprintf(title, "Correction for Subdet %d #eta %d depth %d (Loop %d)", subdet, zside * ieta, depth, loop);
1151     TH1D *hist = new TH1D(name, title, 100, 0.0, 5.0);
1152     hist->Sumw2();
1153     if (debug)
1154       std::cout << "Book Histo " << k << " " << title << std::endl;
1155     histos_[detIds_[k]] = hist;
1156   }
1157   std::cout << "Total of " << detIds_.size() << " detIds and " << histos_.size() << std::endl;
1158 }
1159 
1160 bool CalibTree::goodTrack() {
1161   bool ok(true);
1162   double cut(2.0);
1163   double pmom = (useGen_ && (t_gentrackP > 0)) ? t_gentrackP : t_p;
1164   if (sysmode_ == 1) {
1165     ok = ((t_qltyFlag) && (t_hmaxNearP < cut) && (t_eMipDR < 1.0) && (t_mindR1 > 1.0) && (pmom > pmin_) &&
1166           (pmom < pmax_));
1167   } else if (sysmode_ == 2) {
1168     ok = ((t_qltyFlag) && (t_qltyPVFlag) && (t_hmaxNearP < cut) && (t_eMipDR < 1.0) && (t_mindR1 > 1.0) &&
1169           (pmom > pmin_) && (pmom < pmax_));
1170   } else if (sysmode_ == 3) {
1171     ok = ((t_selectTk) && (t_hmaxNearP < cut) && (t_eMipDR < 1.0) && (t_mindR1 > 1.0) && (pmom > pmin_) &&
1172           (pmom < pmax_));
1173   } else if (sysmode_ == 4) {
1174     ok = ((t_selectTk) && (t_qltyMissFlag) && (t_hmaxNearP < 0.0) && (t_eMipDR < 1.0) && (t_mindR1 > 1.0) &&
1175           (pmom > pmin_) && (pmom < pmax_));
1176   } else if (sysmode_ == 5) {
1177     ok = ((t_selectTk) && (t_qltyMissFlag) && (t_hmaxNearP < cut) && (t_eMipDR < 0.5) && (t_mindR1 > 1.0) &&
1178           (pmom > pmin_) && (pmom < pmax_));
1179   } else if (sysmode_ == 6) {
1180     ok = ((t_selectTk) && (t_qltyMissFlag) && (t_hmaxNearP < cut) && (t_eMipDR < 2.0) && (t_mindR1 > 1.0) &&
1181           (pmom > pmin_) && (pmom < pmax_));
1182   } else if (sysmode_ == 7) {
1183     ok = ((t_selectTk) && (t_qltyMissFlag) && (t_hmaxNearP < cut) && (t_eMipDR < 1.0) && (t_mindR1 > 0.5) &&
1184           (pmom > pmin_) && (pmom < pmax_));
1185   } else {
1186     if (sysmode_ < 0) {
1187       double eta = (t_ieta > 0) ? t_ieta : -t_ieta;
1188       if (sysmode_ == -2)
1189         cut = 8.0 * exp(eta * log2by18_);
1190       else
1191         cut = 10.0;
1192     }
1193     ok = ((t_selectTk) && (t_qltyMissFlag) && (t_hmaxNearP < cut) && (t_eMipDR < 1.0) && (t_mindR1 > 1.0) &&
1194           (pmom > pmin_) && (pmom < pmax_));
1195   }
1196   return ok;
1197 }
1198 
1199 void CalibTree::writeCorrFactor(const char *corrFileName, int ietaMax) {
1200   std::ofstream myfile;
1201   myfile.open(corrFileName);
1202   if (!myfile.is_open()) {
1203     std::cout << "** ERROR: Can't open '" << corrFileName << std::endl;
1204   } else {
1205     myfile << "#" << std::setprecision(4) << std::setw(10) << "detId" << std::setw(10) << "ieta" << std::setw(10)
1206            << "depth" << std::setw(15) << "corrFactor" << std::endl;
1207     std::map<unsigned int, std::pair<double, double> >::const_iterator itr;
1208     for (itr = Cprev.begin(); itr != Cprev.end(); ++itr) {
1209       unsigned int detId = itr->first;
1210       int subdet, depth, zside, ieta, iphi;
1211       unpackDetId(detId, subdet, zside, ieta, iphi, depth);
1212       if (ieta <= ietaMax) {
1213         myfile << std::setw(10) << std::hex << detId << std::setw(10) << std::dec << zside * ieta << std::setw(10)
1214                << depth << std::setw(10) << itr->second.first << " " << std::setw(10) << itr->second.second
1215                << std::endl;
1216         std::cout << itr->second.first << ",";
1217       }
1218     }
1219     myfile.close();
1220     std::cout << std::endl;
1221   }
1222 }
1223 
1224 bool CalibTree::selectPhi(unsigned int detId) {
1225   bool flag(true);
1226   // The masks are defined in DataFormats/HcalDetId/interface/HcalDetId.h

1227   if ((phimin_ > 1) || (phimax_ < 72) || (zside_ != 0)) {
1228     int subdet, depth, zside, ieta, iphi;
1229     unpackDetId(detId, subdet, zside, ieta, iphi, depth);
1230     if (phimin_ > 1 || phimax_ < 72) {
1231       if ((iphi < phimin_) || (iphi > phimax_))
1232         flag = false;
1233     }
1234     if (zside_ != 0) {
1235       if (zside != zside_)
1236         flag = false;
1237     }
1238   }
1239   return flag;
1240 }
1241 
1242 std::pair<double, double> CalibTree::fitMean(TH1D *hist, int mode) {
1243   std::pair<double, double> results = std::pair<double, double>(1, 0);
1244   if (hist != 0) {
1245     double mean = hist->GetMean(), rms = hist->GetRMS();
1246     double LowEdge(0.7), HighEdge(1.3);
1247     char option[20];
1248     if (mode == 1) {
1249       LowEdge = mean - 1.5 * rms;
1250       HighEdge = mean + 1.5 * rms;
1251       int nbin = hist->GetNbinsX();
1252       if (LowEdge < hist->GetBinLowEdge(1))
1253         LowEdge = hist->GetBinLowEdge(1);
1254       if (HighEdge > hist->GetBinLowEdge(nbin) + hist->GetBinWidth(nbin))
1255         HighEdge = hist->GetBinLowEdge(nbin) + hist->GetBinWidth(nbin);
1256     }
1257     if (hist->GetEntries() > 100)
1258       sprintf(option, "+QRLS");
1259     else
1260       sprintf(option, "+QRWLS");
1261     double value(mean);
1262     double error = rms / sqrt(hist->GetEntries());
1263     if (hist->GetEntries() > 20) {
1264       TFitResultPtr Fit = hist->Fit("gaus", option, "", LowEdge, HighEdge);
1265       value = Fit->Value(1);
1266       error = Fit->FitResult::Error(1);
1267       /*

1268       LowEdge  = value - 1.5*error;

1269       HighEdge = value + 1.5*error;

1270       value    = Fit->Value(1);

1271       error    = Fit->FitResult::Error(1); 

1272       Fit = hist->Fit("gaus",option,"",LowEdge,HighEdge);

1273       */
1274     }
1275     results = std::pair<double, double>(value, error);
1276   }
1277   return results;
1278 }
1279 
1280 void CalibTree::makeplots(
1281     double rmin, double rmax, int ietaMax, bool useweight, double fraction, bool debug, Long64_t nmax) {
1282   if (fChain == 0)
1283     return;
1284   Long64_t nentryTot = fChain->GetEntriesFast();
1285   Long64_t nentries = (fraction > 0.01 && fraction < 0.99) ? (Long64_t)(fraction * nentryTot) : nentryTot;
1286   if ((nentries > nmax) && (nmax > 0))
1287     nentries = nmax;
1288 
1289   // Book the histograms

1290   std::map<int, std::pair<TH1D *, TH1D *> > histos;
1291   for (int ieta = -ietaMax; ieta <= ietaMax; ++ieta) {
1292     char name[20], title[100];
1293     sprintf(name, "begin%d", ieta);
1294     if (ieta == 0)
1295       sprintf(title, "Ratio at start");
1296     else
1297       sprintf(title, "Ratio at start for i#eta=%d", ieta);
1298     TH1D *h1 = new TH1D(name, title, 50, rmin, rmax);
1299     h1->Sumw2();
1300     sprintf(name, "end%d", ieta);
1301     if (ieta == 0)
1302       sprintf(title, "Ratio at the end");
1303     else
1304       sprintf(title, "Ratio at the end for i#eta=%d", ieta);
1305     TH1D *h2 = new TH1D(name, title, 50, rmin, rmax);
1306     h2->Sumw2();
1307     histos[ieta] = std::pair<TH1D *, TH1D *>(h1, h2);
1308   }
1309   //Fill the histograms

1310   Long64_t nbytes(0), nb(0);
1311   for (Long64_t jentry = 0; jentry < nentries; jentry++) {
1312     Long64_t ientry = LoadTree(jentry);
1313     nb = fChain->GetEntry(jentry);
1314     nbytes += nb;
1315     if (ientry < 0)
1316       break;
1317     bool select = ((cDuplicate_ != nullptr) && (duplicate_ == 0)) ? (cDuplicate_->isDuplicate(jentry)) : true;
1318     if (!select)
1319       continue;
1320     if (cDuplicate_ != nullptr) {
1321       select = !(cDuplicate_->select(t_ieta, t_iphi));
1322       if (!select)
1323         continue;
1324     }
1325     if (goodTrack()) {
1326       double pmom = (useGen_ && (t_gentrackP > 0)) ? t_gentrackP : t_p;
1327       CalibTree::energyCalor en1 = energyHcal(pmom, jentry, false);
1328       CalibTree::energyCalor en2 = energyHcal(pmom, jentry, true);
1329       if ((en1.ehcal > 0.001) && (en2.ehcal > 0.001)) {
1330         double evWt = (useweight) ? t_EventWeight : 1.0;
1331         double ratioi = en1.ehcal / (pmom - t_eMipDR);
1332         double ratiof = en2.ehcal / (pmom - t_eMipDR);
1333         if (t_ieta >= -ietaMax && t_ieta <= ietaMax && t_ieta != 0) {
1334           if (ratioi >= rmin && ratioi <= rmax) {
1335             histos[0].first->Fill(ratioi, evWt);
1336             histos[t_ieta].first->Fill(ratioi, evWt);
1337           }
1338           if (ratiof >= rmin && ratiof <= rmax) {
1339             histos[0].second->Fill(ratiof, evWt);
1340             histos[t_ieta].second->Fill(ratiof, evWt);
1341           }
1342         }
1343       }
1344     }
1345   }
1346 
1347   //Fit the histograms

1348   TH1D *hbef1 = new TH1D("Eta1Bf", "Mean vs i#eta", 2 * ietaMax, -ietaMax, ietaMax);
1349   TH1D *hbef2 = new TH1D("Eta2Bf", "Median vs i#eta", 2 * ietaMax, -ietaMax, ietaMax);
1350   TH1D *haft1 = new TH1D("Eta1Af", "Mean vs i#eta", 2 * ietaMax, -ietaMax, ietaMax);
1351   TH1D *haft2 = new TH1D("Eta2Af", "Median vs i#eta", 2 * ietaMax, -ietaMax, ietaMax);
1352   for (int ieta = -ietaMax; ieta <= ietaMax; ++ieta) {
1353     int bin = (ieta < 0) ? (ieta + ietaMax + 1) : (ieta + ietaMax);
1354     TH1D *h1 = histos[ieta].first;
1355     double mean1 = h1->GetMean();
1356     double err1 = h1->GetMeanError();
1357     std::pair<double, double> fit1 = fitMean(h1, 1);
1358     if (debug) {
1359       std::cout << ieta << " " << h1->GetName() << " " << mean1 << " +- " << err1 << " and " << fit1.first << " +- "
1360                 << fit1.second << std::endl;
1361     }
1362     if (ieta != 0) {
1363       hbef1->SetBinContent(bin, mean1);
1364       hbef1->SetBinError(bin, err1);
1365       hbef2->SetBinContent(bin, fit1.first);
1366       hbef2->SetBinError(bin, fit1.second);
1367     }
1368     h1->Write();
1369     TH1D *h2 = histos[ieta].second;
1370     double mean2 = h2->GetMean();
1371     double err2 = h2->GetMeanError();
1372     std::pair<double, double> fit2 = fitMean(h2, 1);
1373     if (debug) {
1374       std::cout << ieta << " " << h2->GetName() << " " << mean2 << " +- " << err2 << " and " << fit2.first << " +- "
1375                 << fit2.second << std::endl;
1376     }
1377     if (ieta != 0) {
1378       haft1->SetBinContent(bin, mean2);
1379       haft1->SetBinError(bin, err2);
1380       haft2->SetBinContent(bin, fit2.first);
1381       haft2->SetBinError(bin, fit2.second);
1382     }
1383     h2->Write();
1384   }
1385   fitPol0(hbef1, debug);
1386   fitPol0(hbef2, debug);
1387   fitPol0(haft1, debug);
1388   fitPol0(haft2, debug);
1389 }
1390 
1391 void CalibTree::fitPol0(TH1D *hist, bool debug) {
1392   hist->GetXaxis()->SetTitle("i#eta");
1393   hist->GetYaxis()->SetTitle("<E_{HCAL}/(p-E_{ECAL})>");
1394   hist->GetYaxis()->SetRangeUser(0.4, 1.6);
1395   TFitResultPtr Fit = hist->Fit("pol0", "+QRWLS");
1396   if (debug) {
1397     std::cout << "Fit to Pol0 to " << hist->GetTitle() << ": " << Fit->Value(0) << " +- " << Fit->FitResult::Error(0)
1398               << std::endl;
1399   }
1400   hist->Write();
1401 }
1402 
1403 void CalibTree::highEtaFactors(int ietaMax, bool debug) {
1404   std::map<unsigned int, std::pair<double, double> >::const_iterator itr;
1405   std::pair<double, double> cfacp, cfacn;
1406   cfacp = cfacn = std::pair<double, double>(1.0, 0.0);
1407   for (itr = Cprev.begin(); itr != Cprev.end(); ++itr) {
1408     unsigned int detid = itr->first;
1409     int subdet, depth, zside, ieta, iphi;
1410     unpackDetId(detid, subdet, zside, ieta, iphi, depth);
1411     if ((ieta == ietaMax) && (depth == 1)) {
1412       if (zside > 0)
1413         cfacp = itr->second;
1414       else
1415         cfacn = itr->second;
1416     }
1417   }
1418   if (debug) {
1419     std::cout << "Correction factor for (" << ietaMax << ",1) = " << cfacp.first << ":" << cfacp.second << " and ("
1420               << -ietaMax << ",1) = " << cfacn.first << ":" << cfacn.second << std::endl;
1421   }
1422   for (itr = Cprev.begin(); itr != Cprev.end(); ++itr) {
1423     unsigned int detid = itr->first;
1424     int subdet, depth, zside, ieta, iphi;
1425     unpackDetId(detid, subdet, zside, ieta, iphi, depth);
1426     if (ieta > ietaMax) {
1427       Cprev[detid] = (zside > 0) ? cfacp : cfacn;
1428       if (debug) {
1429         std::cout << "Set correction factor for (" << zside * ieta << "," << depth << ") = " << Cprev[detid].first
1430                   << ":" << Cprev[detid].second << std::endl;
1431       }
1432     }
1433   }
1434 }
1435 
1436 CalibTree::energyCalor CalibTree::energyHcal(double pmom, const Long64_t &entry, bool final) {
1437   double etot = t_eHcal;
1438   double etot2 = t_eHcal;
1439   double ediff = (t_eHcal30 - t_eHcal10);
1440   if (final) {
1441     etot = etot2 = 0;
1442     for (unsigned int idet = 0; idet < (*t_DetIds).size(); idet++) {
1443       // Apply thresholds if necessary

1444       bool okcell = (thrForm_ == 0) || ((*t_HitEnergies)[idet] > threshold((*t_DetIds)[idet], thrForm_));
1445       if (okcell && selectPhi((*t_DetIds)[idet])) {
1446         unsigned int id = (*t_DetIds)[idet];
1447         double hitEn(0);
1448         unsigned int detid = truncateId(id, truncateFlag_, false);
1449         if (Cprev.find(detid) != Cprev.end())
1450           hitEn = Cprev[detid].first * (*t_HitEnergies)[idet];
1451         else
1452           hitEn = (*t_HitEnergies)[idet];
1453         if ((rcorForm_ != 3) && (rcorForm_ >= 0) && (cFactor_))
1454           hitEn *= cFactor_->getCorr(t_Run, id);
1455         if ((cDuplicate_ != nullptr) && (cDuplicate_->doCorr()))
1456           hitEn *= cDuplicate_->getWeight(id);
1457         etot += hitEn;
1458         etot2 += ((*t_HitEnergies)[idet]);
1459       }
1460     }
1461     // Now the outer cone

1462     double etot1(0), etot3(0);
1463     if (t_DetIds1 != 0 && t_DetIds3 != 0) {
1464       for (unsigned int idet = 0; idet < (*t_DetIds1).size(); idet++) {
1465         // Apply thresholds if necessary

1466         bool okcell = (thrForm_ == 0) || ((*t_HitEnergies1)[idet] > threshold((*t_DetIds1)[idet], thrForm_));
1467         if (okcell && selectPhi((*t_DetIds1)[idet])) {
1468           unsigned int id = (*t_DetIds1)[idet];
1469           unsigned int detid = truncateId(id, truncateFlag_, false);
1470           double hitEn(0);
1471           if (Cprev.find(detid) != Cprev.end())
1472             hitEn = Cprev[detid].first * (*t_HitEnergies1)[idet];
1473           else
1474             hitEn = (*t_HitEnergies1)[idet];
1475           if ((rcorForm_ != 3) && (rcorForm_ >= 0) && (cFactor_))
1476             hitEn *= cFactor_->getCorr(t_Run, id);
1477           if ((cDuplicate_ != nullptr) && (cDuplicate_->doCorr()))
1478             hitEn *= cDuplicate_->getWeight(id);
1479           etot1 += hitEn;
1480         }
1481       }
1482       for (unsigned int idet = 0; idet < (*t_DetIds3).size(); idet++) {
1483         // Apply thresholds if necessary

1484         bool okcell = (thrForm_ == 0) || ((*t_HitEnergies3)[idet] > threshold((*t_DetIds3)[idet], thrForm_));
1485         if (okcell && selectPhi((*t_DetIds3)[idet])) {
1486           unsigned int id = (*t_DetIds3)[idet];
1487           unsigned int detid = truncateId(id, truncateFlag_, false);
1488           double hitEn(0);
1489           if (Cprev.find(detid) != Cprev.end())
1490             hitEn = Cprev[detid].first * (*t_HitEnergies3)[idet];
1491           else
1492             hitEn = (*t_HitEnergies3)[idet];
1493           if ((rcorForm_ != 3) && (rcorForm_ >= 0) && (cFactor_))
1494             hitEn *= cFactor_->getCorr(t_Run, id);
1495           if ((cDuplicate_ != nullptr) && (cDuplicate_->doCorr()))
1496             hitEn *= cDuplicate_->getWeight(id);
1497           etot3 += hitEn;
1498         }
1499       }
1500     }
1501     ediff = etot3 - etot1;
1502   }
1503   // PU correction only for loose isolation cut

1504   double ehcal = (((rcorForm_ == 3) && (cFactor_ != nullptr))
1505                       ? (etot * cFactor_->getCorr(entry))
1506                       : ((puCorr_ == 0) ? etot
1507                                         : ((puCorr_ < 0) ? (etot * puFactor(-puCorr_, t_ieta, pmom, etot, ediff))
1508                                                          : puFactorRho(puCorr_, t_ieta, t_rhoh, etot))));
1509   return CalibTree::energyCalor(etot, etot2, ehcal);
1510 }