Back to home page

Project CMSSW displayed by LXR

 
 

    


File indexing completed on 2024-06-07 02:29:09

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) ignore

0058 //                              depth index in HE (depth index set to 1); (4)

0059 //                              ignore depth index in HB (depth index set 1);

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

0061 //                              as depth 2; (6) for depth = 1 and 2, depth =

0062 //                              1, else depth = 2; (7) in case of HB, depths

0063 //                              1 and 2 are set to 1, else depth =2; for HE

0064 //                              ignore depth index; (8) in case of HE, depths

0065 //                              1 and 2 are set to 1, else depth =2; for HB

0066 //                              ignore depth index; (9) Ignore depth index for

0067 //                              depth > 1 in HB and all depth index for HE.

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

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

0070 //                              (Default 0)

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

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

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

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

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

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

0077 //                              from phi-symmetry.

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

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

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

0081 //                              selected.

0082 //                              For threshold h: the format for threshold

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

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

0085 //                              (Default 0)

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

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

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

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

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

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

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

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

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

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

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

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

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

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

0100 //                              For HEP17 it will be 217

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

0102 //                              be exluded or only considered (false)

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

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

0105 //                              and depth 1 (1) (default 1)

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

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

0108 //                              in o/p file (false)

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

0110 //                              estimating the correction factor

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

0112 //                              estimating the correction factor

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

0114 //                              (false)

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

0116 //                              if -1, all entries to be processed; -2 take

0117 //                              all odd entries; -3 take all even entries (-1)

0118 //

0119 //  doIt(inFileName, dupFileName)

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

0121 ////////////////////////////////////////////////////////////////////////////////

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

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

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

0259   TH1D *h_pbyE, *h_cvg;
0260   TProfile *h_Ebyp_bfr, *h_Ebyp_aftr;
0261 
0262 private:
0263   // Declaration of leaf types

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

0304   TBranch *b_t_Run;           //!

0305   TBranch *b_t_Event;         //!

0306   TBranch *b_t_DataType;      //!

0307   TBranch *b_t_ieta;          //!

0308   TBranch *b_t_iphi;          //!

0309   TBranch *b_t_EventWeight;   //!

0310   TBranch *b_t_nVtx;          //!

0311   TBranch *b_t_nTrk;          //!

0312   TBranch *b_t_goodPV;        //!

0313   TBranch *b_t_l1pt;          //!

0314   TBranch *b_t_l1eta;         //!

0315   TBranch *b_t_l1phi;         //!

0316   TBranch *b_t_l3pt;          //!

0317   TBranch *b_t_l3eta;         //!

0318   TBranch *b_t_l3phi;         //!

0319   TBranch *b_t_p;             //!

0320   TBranch *b_t_pt;            //!

0321   TBranch *b_t_phi;           //!

0322   TBranch *b_t_mindR1;        //!

0323   TBranch *b_t_mindR2;        //!

0324   TBranch *b_t_eMipDR;        //!

0325   TBranch *b_t_eHcal;         //!

0326   TBranch *b_t_eHcal10;       //!

0327   TBranch *b_t_eHcal30;       //!

0328   TBranch *b_t_hmaxNearP;     //!

0329   TBranch *b_t_rhoh;          //!

0330   TBranch *b_t_selectTk;      //!

0331   TBranch *b_t_qltyFlag;      //!

0332   TBranch *b_t_qltyMissFlag;  //!

0333   TBranch *b_t_qltyPVFlag;    //!

0334   TBranch *b_t_gentrackP;     //!

0335   TBranch *b_t_DetIds;        //!

0336   TBranch *b_t_HitEnergies;   //!

0337   TBranch *b_t_trgbits;       //!

0338   TBranch *b_t_DetIds1;       //!

0339   TBranch *b_t_DetIds3;       //!

0340   TBranch *b_t_HitEnergies1;  //!

0341   TBranch *b_t_HitEnergies3;  //!

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

0639   if (!fChain)
0640     return 0;
0641   return fChain->GetEntry(entry);
0642 }
0643 
0644 Long64_t CalibTree::LoadTree(Long64_t entry) {
0645   // Set the environment to read one entry

0646   if (!fChain)
0647     return -5;
0648   Long64_t centry = fChain->LoadTree(entry);
0649   if (centry < 0)
0650     return centry;
0651   if (fChain->GetTreeNumber() != fCurrent) {
0652     fCurrent = fChain->GetTreeNumber();
0653     Notify();
0654   }
0655   return centry;
0656 }
0657 
0658 void CalibTree::Init(TChain *tree) {
0659   // The Init() function is called when the selector needs to initialize

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

0661   // pointers of the tree will be set.

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

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

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

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

0666 
0667   // Set object pointer

0668   t_DetIds = 0;
0669   t_HitEnergies = 0;
0670   t_trgbits = 0;
0671   t_DetIds1 = 0;
0672   t_DetIds3 = 0;
0673   t_HitEnergies1 = 0;
0674   t_HitEnergies3 = 0;
0675   // Set branch addresses and branch pointers

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

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

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

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

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

0729 
0730   return kTRUE;
0731 }
0732 
0733 void CalibTree::Show(Long64_t entry) {
0734   // Print contents of entry.

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

0736   if (!fChain)
0737     return;
0738   fChain->Show(entry);
0739 }
0740 
0741 Int_t CalibTree::Cut(Long64_t) {
0742   // This function may be called from Loop.

0743   // returns  1 if entry is accepted.

0744   // returns -1 otherwise.

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

0846             bool okcell = (thrForm_ == 0) || ((*t_HitEnergies)[idet] > threshold((*t_DetIds)[idet], thrForm_));
0847             if (okcell && selectPhi((*t_DetIds)[idet])) {
0848               unsigned int id = (*t_DetIds)[idet];
0849               unsigned int detid = truncateId(id, truncateFlag_, false);
0850               double hitEn = 0.0;
0851               if (debug) {
0852                 std::cout << "idet " << idet << " detid/hitenergy : " << std::hex << (*t_DetIds)[idet] << ":" << detid
0853                           << "/" << (*t_HitEnergies)[idet] << std::endl;
0854               }
0855               if (Cprev.find(detid) != Cprev.end())
0856                 hitEn = Cprev[detid].first * (*t_HitEnergies)[idet];
0857               else
0858                 hitEn = (*t_HitEnergies)[idet];
0859               if ((rcorForm_ != 3) && (rcorForm_ >= 0) && (cFactor_))
0860                 hitEn *= cFactor_->getCorr(t_Run, id);
0861               if ((cDuplicate_ != nullptr) && (cDuplicate_->doCorr()))
0862                 hitEn *= cDuplicate_->getWeight(id);
0863               double Wi = evWt * hitEn / en.Etot;
0864               double Fac = (inverse) ? (en.ehcal / (pmom - t_eMipDR)) : ((pmom - t_eMipDR) / en.ehcal);
0865               double Fac2 = Wi * Fac * Fac;
0866               TH1D *hist(0);
0867               std::map<unsigned int, TH1D *>::const_iterator itr = histos_.find(detid);
0868               if (itr != histos_.end())
0869                 hist = itr->second;
0870               if (debug) {
0871                 std::cout << "Det Id " << std::hex << detid << std::dec << " " << hist << std::endl;
0872               }
0873               if (hist != 0)
0874                 hist->Fill(Fac, Wi);  //////histola

0875               Fac *= Wi;
0876               if (SumW.find(detid) != SumW.end()) {
0877                 Wi += SumW[detid].fact0;
0878                 Fac += SumW[detid].fact1;
0879                 Fac2 += SumW[detid].fact2;
0880                 int kount = SumW[detid].kount + 1;
0881                 SumW[detid] = myEntry(kount, Wi, Fac, Fac2);
0882                 nTrks[detid] += evWt;
0883               } else {
0884                 SumW.insert(std::pair<unsigned int, myEntry>(detid, myEntry(1, Wi, Fac, Fac2)));
0885                 nTrks.insert(std::pair<unsigned int, unsigned int>(detid, evWt));
0886               }
0887             }
0888           }
0889         }
0890       }
0891     }
0892   }
0893   if (debug)
0894     std::cout << "# of Good Tracks " << ntkgood << " out of " << nentries << std::endl;
0895   if (loop == 0) {
0896     h_pbyE->Write("h_pbyE");
0897     h_Ebyp_bfr->Write("h_Ebyp_bfr");
0898   }
0899   if (last) {
0900     h_Ebyp_aftr->Write("h_Ebyp_aftr");
0901   }
0902 
0903   std::map<unsigned int, std::pair<double, double> > cfactors;
0904   unsigned int kount(0), kountus(0);
0905   double sumfactor(0);
0906   for (std::map<unsigned int, TH1D *>::const_iterator itr = histos_.begin(); itr != histos_.end(); ++itr) {
0907     if (writeHisto) {
0908       //      std::pair<double, double> result_write = fitMean(itr->second, 0);

0909       (itr->second)->Write();
0910     }
0911     // The masks are defined in DataFormats/HcalDetId/interface/HcalDetId.h

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

1095       bool selRun = (includeRun_ ? ((t_Run >= runlo_) && (t_Run <= runhi_)) : ((t_Run < runlo_) || (t_Run > runhi_)));
1096       bool selTrack = ((ietaTrack <= 0) || (abs(t_ieta) <= ietaTrack));
1097       if (selRun && (t_nVtx >= nvxlo_) && (t_nVtx <= nvxhi_) && selTrack) {
1098         bool isItRBX(false);
1099         if (cSelect_ != nullptr) {
1100           bool temp = cSelect_->isItRBX(t_DetIds);
1101           if (exclude_)
1102             isItRBX = temp;
1103           else
1104             isItRBX = !(temp);
1105         }
1106         if ((cDuplicate_ != nullptr) && (!isItRBX))
1107           isItRBX = (cDuplicate_->select(t_ieta, t_iphi));
1108         ++kprint;
1109         if (!(isItRBX)) {
1110           for (unsigned int idet = 0; idet < (*t_DetIds).size(); idet++) {
1111             if (selectPhi((*t_DetIds)[idet])) {
1112               unsigned int detid = truncateId((*t_DetIds)[idet], truncateFlag_, debug);
1113               if (debug && (kprint <= 10)) {
1114                 std::cout << "DetId[" << idet << "] Original " << std::hex << (*t_DetIds)[idet] << " truncated "
1115                           << detid << std::dec;
1116               }
1117               if (std::find(detIds_.begin(), detIds_.end(), detid) == detIds_.end()) {
1118                 detIds_.push_back(detid);
1119                 if (debug && (kprint <= 10))
1120                   std::cout << " new";
1121               }
1122               if (debug && (kprint <= 10))
1123                 std::cout << std::endl;
1124             }
1125           }
1126           // Also look at the neighbouring cells if available

1127           if (t_DetIds3 != 0) {
1128             for (unsigned int idet = 0; idet < (*t_DetIds3).size(); idet++) {
1129               if (selectPhi((*t_DetIds3)[idet])) {
1130                 unsigned int detid = truncateId((*t_DetIds3)[idet], truncateFlag_, debug);
1131                 if (std::find(detIds_.begin(), detIds_.end(), detid) == detIds_.end()) {
1132                   detIds_.push_back(detid);
1133                 }
1134               }
1135             }
1136           }
1137         }
1138       }
1139     }
1140   }
1141   if (debug) {
1142     std::cout << "Total of " << detIds_.size() << " detIds" << std::endl;
1143     // The masks are defined in DataFormats/HcalDetId/interface/HcalDetId.h

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

1248   if ((phimin_ > 1) || (phimax_ < 72) || (zside_ != 0)) {
1249     int subdet, depth, zside, ieta, iphi;
1250     unpackDetId(detId, subdet, zside, ieta, iphi, depth);
1251     if (phimin_ > 1 || phimax_ < 72) {
1252       if ((iphi < phimin_) || (iphi > phimax_))
1253         flag = false;
1254     }
1255     if (zside_ != 0) {
1256       if (zside != zside_)
1257         flag = false;
1258     }
1259   }
1260   return flag;
1261 }
1262 
1263 std::pair<double, double> CalibTree::fitMean(TH1D *hist, int mode) {
1264   std::pair<double, double> results = std::pair<double, double>(1, 0);
1265   if (hist != 0) {
1266     double mean = hist->GetMean(), rms = hist->GetRMS();
1267     double LowEdge(0.7), HighEdge(1.3);
1268     char option[20];
1269     if (mode == 1) {
1270       LowEdge = mean - 1.5 * rms;
1271       HighEdge = mean + 1.5 * rms;
1272       int nbin = hist->GetNbinsX();
1273       if (LowEdge < hist->GetBinLowEdge(1))
1274         LowEdge = hist->GetBinLowEdge(1);
1275       if (HighEdge > hist->GetBinLowEdge(nbin) + hist->GetBinWidth(nbin))
1276         HighEdge = hist->GetBinLowEdge(nbin) + hist->GetBinWidth(nbin);
1277     }
1278     if (hist->GetEntries() > 100)
1279       sprintf(option, "+QRLS");
1280     else
1281       sprintf(option, "+QRWLS");
1282     double value(mean);
1283     double error = rms / sqrt(hist->GetEntries());
1284     if (hist->GetEntries() > 20) {
1285       TFitResultPtr Fit = hist->Fit("gaus", option, "", LowEdge, HighEdge);
1286       value = Fit->Value(1);
1287       error = Fit->FitResult::Error(1);
1288       /*

1289       LowEdge  = value - 1.5*error;

1290       HighEdge = value + 1.5*error;

1291       value    = Fit->Value(1);

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

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

1294       */
1295     }
1296     results = std::pair<double, double>(value, error);
1297   }
1298   return results;
1299 }
1300 
1301 void CalibTree::makeplots(
1302     double rmin, double rmax, int ietaMax, bool useweight, double fraction, bool debug, Long64_t nmax) {
1303   if (fChain == 0)
1304     return;
1305   Long64_t nentryTot = fChain->GetEntriesFast();
1306   Long64_t nentries = (fraction > 0.01 && fraction < 0.99) ? (Long64_t)(fraction * nentryTot) : nentryTot;
1307   int32_t oddEven = (nmax == -2) ? 1 : ((nmax == -3) ? -1 : 0);
1308   if ((nentries > nmax) && (nmax > 0))
1309     nentries = nmax;
1310 
1311   // Book the histograms

1312   std::map<int, std::pair<TH1D *, TH1D *> > histos;
1313   for (int ieta = -ietaMax; ieta <= ietaMax; ++ieta) {
1314     char name[20], title[100];
1315     sprintf(name, "begin%d", ieta);
1316     if (ieta == 0)
1317       sprintf(title, "Ratio at start");
1318     else
1319       sprintf(title, "Ratio at start for i#eta=%d", ieta);
1320     TH1D *h1 = new TH1D(name, title, 50, rmin, rmax);
1321     h1->Sumw2();
1322     sprintf(name, "end%d", ieta);
1323     if (ieta == 0)
1324       sprintf(title, "Ratio at the end");
1325     else
1326       sprintf(title, "Ratio at the end for i#eta=%d", ieta);
1327     TH1D *h2 = new TH1D(name, title, 50, rmin, rmax);
1328     h2->Sumw2();
1329     histos[ieta] = std::pair<TH1D *, TH1D *>(h1, h2);
1330   }
1331   //Fill the histograms

1332   Long64_t nbytes(0), nb(0);
1333   for (Long64_t jentry = 0; jentry < nentries; jentry++) {
1334     Long64_t ientry = LoadTree(jentry);
1335     nb = fChain->GetEntry(jentry);
1336     nbytes += nb;
1337     if (ientry < 0)
1338       break;
1339     if (oddEven != 0) {
1340       if ((oddEven < 0) && (jentry % 2 == 0))
1341         continue;
1342       else if ((oddEven > 0) && (jentry % 2 != 0))
1343         continue;
1344     }
1345     bool select = ((cDuplicate_ != nullptr) && (duplicate_ == 0)) ? (cDuplicate_->isDuplicate(jentry)) : true;
1346     if (!select)
1347       continue;
1348     if (cDuplicate_ != nullptr) {
1349       select = !(cDuplicate_->select(t_ieta, t_iphi));
1350       if (!select)
1351         continue;
1352     }
1353     if (goodTrack()) {
1354       double pmom = (useGen_ && (t_gentrackP > 0)) ? t_gentrackP : t_p;
1355       CalibTree::energyCalor en1 = energyHcal(pmom, jentry, false);
1356       CalibTree::energyCalor en2 = energyHcal(pmom, jentry, true);
1357       if ((en1.ehcal > 0.001) && (en2.ehcal > 0.001)) {
1358         double evWt = (useweight) ? t_EventWeight : 1.0;
1359         double ratioi = en1.ehcal / (pmom - t_eMipDR);
1360         double ratiof = en2.ehcal / (pmom - t_eMipDR);
1361         if (t_ieta >= -ietaMax && t_ieta <= ietaMax && t_ieta != 0) {
1362           if (ratioi >= rmin && ratioi <= rmax) {
1363             histos[0].first->Fill(ratioi, evWt);
1364             histos[t_ieta].first->Fill(ratioi, evWt);
1365           }
1366           if (ratiof >= rmin && ratiof <= rmax) {
1367             histos[0].second->Fill(ratiof, evWt);
1368             histos[t_ieta].second->Fill(ratiof, evWt);
1369           }
1370         }
1371       }
1372     }
1373   }
1374 
1375   //Fit the histograms

1376   TH1D *hbef1 = new TH1D("Eta1Bf", "Mean vs i#eta", 2 * ietaMax, -ietaMax, ietaMax);
1377   TH1D *hbef2 = new TH1D("Eta2Bf", "Median vs i#eta", 2 * ietaMax, -ietaMax, ietaMax);
1378   TH1D *haft1 = new TH1D("Eta1Af", "Mean vs i#eta", 2 * ietaMax, -ietaMax, ietaMax);
1379   TH1D *haft2 = new TH1D("Eta2Af", "Median vs i#eta", 2 * ietaMax, -ietaMax, ietaMax);
1380   for (int ieta = -ietaMax; ieta <= ietaMax; ++ieta) {
1381     int bin = (ieta < 0) ? (ieta + ietaMax + 1) : (ieta + ietaMax);
1382     TH1D *h1 = histos[ieta].first;
1383     double mean1 = h1->GetMean();
1384     double err1 = h1->GetMeanError();
1385     std::pair<double, double> fit1 = fitMean(h1, 1);
1386     if (debug) {
1387       std::cout << ieta << " " << h1->GetName() << " " << mean1 << " +- " << err1 << " and " << fit1.first << " +- "
1388                 << fit1.second << std::endl;
1389     }
1390     if (ieta != 0) {
1391       hbef1->SetBinContent(bin, mean1);
1392       hbef1->SetBinError(bin, err1);
1393       hbef2->SetBinContent(bin, fit1.first);
1394       hbef2->SetBinError(bin, fit1.second);
1395     }
1396     h1->Write();
1397     TH1D *h2 = histos[ieta].second;
1398     double mean2 = h2->GetMean();
1399     double err2 = h2->GetMeanError();
1400     std::pair<double, double> fit2 = fitMean(h2, 1);
1401     if (debug) {
1402       std::cout << ieta << " " << h2->GetName() << " " << mean2 << " +- " << err2 << " and " << fit2.first << " +- "
1403                 << fit2.second << std::endl;
1404     }
1405     if (ieta != 0) {
1406       haft1->SetBinContent(bin, mean2);
1407       haft1->SetBinError(bin, err2);
1408       haft2->SetBinContent(bin, fit2.first);
1409       haft2->SetBinError(bin, fit2.second);
1410     }
1411     h2->Write();
1412   }
1413   fitPol0(hbef1, debug);
1414   fitPol0(hbef2, debug);
1415   fitPol0(haft1, debug);
1416   fitPol0(haft2, debug);
1417 }
1418 
1419 void CalibTree::fitPol0(TH1D *hist, bool debug) {
1420   hist->GetXaxis()->SetTitle("i#eta");
1421   hist->GetYaxis()->SetTitle("<E_{HCAL}/(p-E_{ECAL})>");
1422   hist->GetYaxis()->SetRangeUser(0.4, 1.6);
1423   TFitResultPtr Fit = hist->Fit("pol0", "+QRWLS");
1424   if (debug) {
1425     std::cout << "Fit to Pol0 to " << hist->GetTitle() << ": " << Fit->Value(0) << " +- " << Fit->FitResult::Error(0)
1426               << std::endl;
1427   }
1428   hist->Write();
1429 }
1430 
1431 void CalibTree::highEtaFactors(int ietaMax, bool debug) {
1432   std::map<unsigned int, std::pair<double, double> >::const_iterator itr;
1433   std::pair<double, double> cfacp, cfacn;
1434   cfacp = cfacn = std::pair<double, double>(1.0, 0.0);
1435   for (itr = Cprev.begin(); itr != Cprev.end(); ++itr) {
1436     unsigned int detid = itr->first;
1437     int subdet, depth, zside, ieta, iphi;
1438     unpackDetId(detid, subdet, zside, ieta, iphi, depth);
1439     if ((ieta == ietaMax) && (depth == 1)) {
1440       if (zside > 0)
1441         cfacp = itr->second;
1442       else
1443         cfacn = itr->second;
1444     }
1445   }
1446   if (debug) {
1447     std::cout << "Correction factor for (" << ietaMax << ",1) = " << cfacp.first << ":" << cfacp.second << " and ("
1448               << -ietaMax << ",1) = " << cfacn.first << ":" << cfacn.second << std::endl;
1449   }
1450   for (itr = Cprev.begin(); itr != Cprev.end(); ++itr) {
1451     unsigned int detid = itr->first;
1452     int subdet, depth, zside, ieta, iphi;
1453     unpackDetId(detid, subdet, zside, ieta, iphi, depth);
1454     if (ieta > ietaMax) {
1455       Cprev[detid] = (zside > 0) ? cfacp : cfacn;
1456       if (debug) {
1457         std::cout << "Set correction factor for (" << zside * ieta << "," << depth << ") = " << Cprev[detid].first
1458                   << ":" << Cprev[detid].second << std::endl;
1459       }
1460     }
1461   }
1462 }
1463 
1464 CalibTree::energyCalor CalibTree::energyHcal(double pmom, const Long64_t &entry, bool final) {
1465   double etot = t_eHcal;
1466   double etot2 = t_eHcal;
1467   double ediff = (t_eHcal30 - t_eHcal10);
1468   if (final) {
1469     etot = etot2 = 0;
1470     for (unsigned int idet = 0; idet < (*t_DetIds).size(); idet++) {
1471       // Apply thresholds if necessary

1472       bool okcell = (thrForm_ == 0) || ((*t_HitEnergies)[idet] > threshold((*t_DetIds)[idet], thrForm_));
1473       if (okcell && selectPhi((*t_DetIds)[idet])) {
1474         unsigned int id = (*t_DetIds)[idet];
1475         double hitEn(0);
1476         unsigned int detid = truncateId(id, truncateFlag_, false);
1477         if (Cprev.find(detid) != Cprev.end())
1478           hitEn = Cprev[detid].first * (*t_HitEnergies)[idet];
1479         else
1480           hitEn = (*t_HitEnergies)[idet];
1481         if ((rcorForm_ != 3) && (rcorForm_ >= 0) && (cFactor_))
1482           hitEn *= cFactor_->getCorr(t_Run, id);
1483         if ((cDuplicate_ != nullptr) && (cDuplicate_->doCorr()))
1484           hitEn *= cDuplicate_->getWeight(id);
1485         etot += hitEn;
1486         etot2 += ((*t_HitEnergies)[idet]);
1487       }
1488     }
1489     // Now the outer cone

1490     double etot1(0), etot3(0);
1491     if (t_DetIds1 != 0 && t_DetIds3 != 0) {
1492       for (unsigned int idet = 0; idet < (*t_DetIds1).size(); idet++) {
1493         // Apply thresholds if necessary

1494         bool okcell = (thrForm_ == 0) || ((*t_HitEnergies1)[idet] > threshold((*t_DetIds1)[idet], thrForm_));
1495         if (okcell && selectPhi((*t_DetIds1)[idet])) {
1496           unsigned int id = (*t_DetIds1)[idet];
1497           unsigned int detid = truncateId(id, truncateFlag_, false);
1498           double hitEn(0);
1499           if (Cprev.find(detid) != Cprev.end())
1500             hitEn = Cprev[detid].first * (*t_HitEnergies1)[idet];
1501           else
1502             hitEn = (*t_HitEnergies1)[idet];
1503           if ((rcorForm_ != 3) && (rcorForm_ >= 0) && (cFactor_))
1504             hitEn *= cFactor_->getCorr(t_Run, id);
1505           if ((cDuplicate_ != nullptr) && (cDuplicate_->doCorr()))
1506             hitEn *= cDuplicate_->getWeight(id);
1507           etot1 += hitEn;
1508         }
1509       }
1510       for (unsigned int idet = 0; idet < (*t_DetIds3).size(); idet++) {
1511         // Apply thresholds if necessary

1512         bool okcell = (thrForm_ == 0) || ((*t_HitEnergies3)[idet] > threshold((*t_DetIds3)[idet], thrForm_));
1513         if (okcell && selectPhi((*t_DetIds3)[idet])) {
1514           unsigned int id = (*t_DetIds3)[idet];
1515           unsigned int detid = truncateId(id, truncateFlag_, false);
1516           double hitEn(0);
1517           if (Cprev.find(detid) != Cprev.end())
1518             hitEn = Cprev[detid].first * (*t_HitEnergies3)[idet];
1519           else
1520             hitEn = (*t_HitEnergies3)[idet];
1521           if ((rcorForm_ != 3) && (rcorForm_ >= 0) && (cFactor_))
1522             hitEn *= cFactor_->getCorr(t_Run, id);
1523           if ((cDuplicate_ != nullptr) && (cDuplicate_->doCorr()))
1524             hitEn *= cDuplicate_->getWeight(id);
1525           etot3 += hitEn;
1526         }
1527       }
1528     }
1529     ediff = etot3 - etot1;
1530   }
1531   // PU correction only for loose isolation cut

1532   double ehcal = (((rcorForm_ == 3) && (cFactor_ != nullptr))
1533                       ? (etot * cFactor_->getCorr(entry))
1534                       : ((puCorr_ == 0) ? etot
1535                                         : ((puCorr_ < 0) ? (etot * puFactor(-puCorr_, t_ieta, pmom, etot, ediff))
1536                                                          : puFactorRho(puCorr_, t_ieta, t_rhoh, etot))));
1537   return CalibTree::energyCalor(etot, etot2, ehcal);
1538 }