Back to home page

Project CMSSW displayed by LXR
 
 

    

File indexing completed on 2025-07-10 23:56:30 

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 //      badRunFile, writeDebugHisto, pmin, pmax, 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) Assign all depths > 4

0065 //                              as depth = 5; (9) Assign all depth = 1 as

0066 //                              depth = 2. The digit *d* is used if zside is

0067 //                              to be ignored (1) or not (0)

0068 //                              (Default 0)

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

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

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

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

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

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

0075 //                              from phi-symmetry; (5) use reults from several

0076 //                              phi-symmetry studies drive by run numeber.

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

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

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

0080 //                              selected; (3) list of run ranges and for each

0081 //                              range, ieta, depth where gain has changed.

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 //                              4: 2025 Begin of Year; 5: Derived from the file

0086 //                              PFCuts_IOV_362975.txt.

0087 //                              (Default 0)

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

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

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

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

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

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

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

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

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

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

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

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

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

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

0102 //                              For HEP17 it will be 217

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

0104 //                              be exluded or only considered (false)

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

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

0107 //                              and depth 1 (1) (default 1)

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

0109 //  badRunFile      (char *)  = Name of the file containing a list of runs

0110 //                              to be excluded

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

0112 //                              in o/p file (false)

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

0114 //                              estimating the correction factor

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

0116 //                              estimating the correction factor

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

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

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

0120 //

0121 //  doIt(inFileName, dupFileName)

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

0123 ////////////////////////////////////////////////////////////////////////////////

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

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

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

0262   TH1D *h_pbyE, *h_cvg;
0263   TProfile *h_Ebyp_bfr, *h_Ebyp_aftr;
0264 
0265 private:
0266   // Declaration of leaf types

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

0307   TBranch *b_t_Run;           //!

0308   TBranch *b_t_Event;         //!

0309   TBranch *b_t_DataType;      //!

0310   TBranch *b_t_ieta;          //!

0311   TBranch *b_t_iphi;          //!

0312   TBranch *b_t_EventWeight;   //!

0313   TBranch *b_t_nVtx;          //!

0314   TBranch *b_t_nTrk;          //!

0315   TBranch *b_t_goodPV;        //!

0316   TBranch *b_t_l1pt;          //!

0317   TBranch *b_t_l1eta;         //!

0318   TBranch *b_t_l1phi;         //!

0319   TBranch *b_t_l3pt;          //!

0320   TBranch *b_t_l3eta;         //!

0321   TBranch *b_t_l3phi;         //!

0322   TBranch *b_t_p;             //!

0323   TBranch *b_t_pt;            //!

0324   TBranch *b_t_phi;           //!

0325   TBranch *b_t_mindR1;        //!

0326   TBranch *b_t_mindR2;        //!

0327   TBranch *b_t_eMipDR;        //!

0328   TBranch *b_t_eHcal;         //!

0329   TBranch *b_t_eHcal10;       //!

0330   TBranch *b_t_eHcal30;       //!

0331   TBranch *b_t_hmaxNearP;     //!

0332   TBranch *b_t_rhoh;          //!

0333   TBranch *b_t_selectTk;      //!

0334   TBranch *b_t_qltyFlag;      //!

0335   TBranch *b_t_qltyMissFlag;  //!

0336   TBranch *b_t_qltyPVFlag;    //!

0337   TBranch *b_t_gentrackP;     //!

0338   TBranch *b_t_DetIds;        //!

0339   TBranch *b_t_HitEnergies;   //!

0340   TBranch *b_t_trgbits;       //!

0341   TBranch *b_t_DetIds1;       //!

0342   TBranch *b_t_DetIds3;       //!

0343   TBranch *b_t_HitEnergies1;  //!

0344   TBranch *b_t_HitEnergies3;  //!

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

0662   if (!fChain)
0663     return 0;
0664   return fChain->GetEntry(entry);
0665 }
0666 
0667 Long64_t CalibTree::LoadTree(Long64_t entry) {
0668   // Set the environment to read one entry

0669   if (!fChain)
0670     return -5;
0671   Long64_t centry = fChain->LoadTree(entry);
0672   if (centry < 0)
0673     return centry;
0674   if (fChain->GetTreeNumber() != fCurrent) {
0675     fCurrent = fChain->GetTreeNumber();
0676     Notify();
0677   }
0678   return centry;
0679 }
0680 
0681 void CalibTree::Init(TChain *tree) {
0682   // The Init() function is called when the selector needs to initialize

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

0684   // pointers of the tree will be set.

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

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

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

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

0689 
0690   // Set object pointer

0691   t_DetIds = 0;
0692   t_HitEnergies = 0;
0693   t_trgbits = 0;
0694   t_DetIds1 = 0;
0695   t_DetIds3 = 0;
0696   t_HitEnergies1 = 0;
0697   t_HitEnergies3 = 0;
0698   // Set branch addresses and branch pointers

0699   fChain = tree;
0700   if (!tree)
0701     return;
0702   fCurrent = -1;
0703   fChain->SetMakeClass(1);
0704 
0705   fChain->SetBranchAddress("t_Run", &t_Run, &b_t_Run);
0706   fChain->SetBranchAddress("t_Event", &t_Event, &b_t_Event);
0707   fChain->SetBranchAddress("t_DataType", &t_DataType, &b_t_DataType);
0708   fChain->SetBranchAddress("t_ieta", &t_ieta, &b_t_ieta);
0709   fChain->SetBranchAddress("t_iphi", &t_iphi, &b_t_iphi);
0710   fChain->SetBranchAddress("t_EventWeight", &t_EventWeight, &b_t_EventWeight);
0711   fChain->SetBranchAddress("t_nVtx", &t_nVtx, &b_t_nVtx);
0712   fChain->SetBranchAddress("t_nTrk", &t_nTrk, &b_t_nTrk);
0713   fChain->SetBranchAddress("t_goodPV", &t_goodPV, &b_t_goodPV);
0714   fChain->SetBranchAddress("t_l1pt", &t_l1pt, &b_t_l1pt);
0715   fChain->SetBranchAddress("t_l1eta", &t_l1eta, &b_t_l1eta);
0716   fChain->SetBranchAddress("t_l1phi", &t_l1phi, &b_t_l1phi);
0717   fChain->SetBranchAddress("t_l3pt", &t_l3pt, &b_t_l3pt);
0718   fChain->SetBranchAddress("t_l3eta", &t_l3eta, &b_t_l3eta);
0719   fChain->SetBranchAddress("t_l3phi", &t_l3phi, &b_t_l3phi);
0720   fChain->SetBranchAddress("t_p", &t_p, &b_t_p);
0721   fChain->SetBranchAddress("t_pt", &t_pt, &b_t_pt);
0722   fChain->SetBranchAddress("t_phi", &t_phi, &b_t_phi);
0723   fChain->SetBranchAddress("t_mindR1", &t_mindR1, &b_t_mindR1);
0724   fChain->SetBranchAddress("t_mindR2", &t_mindR2, &b_t_mindR2);
0725   fChain->SetBranchAddress("t_eMipDR", &t_eMipDR, &b_t_eMipDR);
0726   fChain->SetBranchAddress("t_eHcal", &t_eHcal, &b_t_eHcal);
0727   fChain->SetBranchAddress("t_eHcal10", &t_eHcal10, &b_t_eHcal10);
0728   fChain->SetBranchAddress("t_eHcal30", &t_eHcal30, &b_t_eHcal30);
0729   fChain->SetBranchAddress("t_hmaxNearP", &t_hmaxNearP, &b_t_hmaxNearP);
0730   fChain->SetBranchAddress("t_rhoh", &t_rhoh, &b_t_rhoh);
0731   fChain->SetBranchAddress("t_selectTk", &t_selectTk, &b_t_selectTk);
0732   fChain->SetBranchAddress("t_qltyFlag", &t_qltyFlag, &b_t_qltyFlag);
0733   fChain->SetBranchAddress("t_qltyMissFlag", &t_qltyMissFlag, &b_t_qltyMissFlag);
0734   fChain->SetBranchAddress("t_qltyPVFlag", &t_qltyPVFlag, &b_t_qltyPVFlag);
0735   fChain->SetBranchAddress("t_gentrackP", &t_gentrackP, &b_t_gentrackP);
0736   fChain->SetBranchAddress("t_DetIds", &t_DetIds, &b_t_DetIds);
0737   fChain->SetBranchAddress("t_HitEnergies", &t_HitEnergies, &b_t_HitEnergies);
0738   fChain->SetBranchAddress("t_trgbits", &t_trgbits, &b_t_trgbits);
0739   fChain->SetBranchAddress("t_DetIds1", &t_DetIds1, &b_t_DetIds1);
0740   fChain->SetBranchAddress("t_DetIds3", &t_DetIds3, &b_t_DetIds3);
0741   fChain->SetBranchAddress("t_HitEnergies1", &t_HitEnergies1, &b_t_HitEnergies1);
0742   fChain->SetBranchAddress("t_HitEnergies3", &t_HitEnergies3, &b_t_HitEnergies3);
0743   Notify();
0744 }
0745 
0746 Bool_t CalibTree::Notify() {
0747   // The Notify() function is called when a new file is opened. This

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

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

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

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

0752 
0753   return kTRUE;
0754 }
0755 
0756 void CalibTree::Show(Long64_t entry) {
0757   // Print contents of entry.

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

0759   if (!fChain)
0760     return;
0761   fChain->Show(entry);
0762 }
0763 
0764 Int_t CalibTree::Cut(Long64_t) {
0765   // This function may be called from Loop.

0766   // returns  1 if entry is accepted.

0767   // returns -1 otherwise.

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

0870             bool okcell = (thrForm_ == 0) || ((*t_HitEnergies)[idet] > (cThr_->threshold((*t_DetIds)[idet])));
0871             if (okcell && selectPhi((*t_DetIds)[idet])) {
0872               unsigned int id = (*t_DetIds)[idet];
0873               unsigned int detid = truncateId(id, truncateFlag_, false);
0874               double hitEn = 0.0;
0875               if (debug) {
0876                 std::cout << "idet " << idet << " detid/hitenergy : " << std::hex << (*t_DetIds)[idet] << ":" << detid
0877                           << "/" << (*t_HitEnergies)[idet] << std::endl;
0878               }
0879               if (Cprev.find(detid) != Cprev.end())
0880                 hitEn = Cprev[detid].first * (*t_HitEnergies)[idet];
0881               else
0882                 hitEn = (*t_HitEnergies)[idet];
0883               if ((rcorForm_ != 3) && (rcorForm_ >= 0) && (cFactor_))
0884                 hitEn *= cFactor_->getCorr(t_Run, id);
0885               if ((cDuplicate_ != nullptr) && (cDuplicate_->doCorr(1)))
0886                 hitEn *= cDuplicate_->getWeight(id);
0887               if ((cDuplicate_ != nullptr) && (cDuplicate_->doCorr(3))) {
0888                 int subdet, zside, ieta, iphi, depth;
0889                 unpackDetId((*t_DetIds)[idet], subdet, zside, ieta, iphi, depth);
0890                 hitEn *= cDuplicate_->getCorr(t_Run, ieta, depth);
0891               }
0892               double Wi = evWt * hitEn / en.Etot;
0893               double Fac = (inverse) ? (en.ehcal / (pmom - t_eMipDR)) : ((pmom - t_eMipDR) / en.ehcal);
0894               double Fac2 = Wi * Fac * Fac;
0895               TH1D *hist(0);
0896               std::map<unsigned int, TH1D *>::const_iterator itr = histos_.find(detid);
0897               if (itr != histos_.end())
0898                 hist = itr->second;
0899               if (debug) {
0900                 std::cout << "Det Id " << std::hex << detid << std::dec << " " << hist << std::endl;
0901               }
0902               if (hist != 0)
0903                 hist->Fill(Fac, Wi);  //////histola

0904               Fac *= Wi;
0905               if (SumW.find(detid) != SumW.end()) {
0906                 Wi += SumW[detid].fact0;
0907                 Fac += SumW[detid].fact1;
0908                 Fac2 += SumW[detid].fact2;
0909                 int kount = SumW[detid].kount + 1;
0910                 SumW[detid] = myEntry(kount, Wi, Fac, Fac2);
0911                 nTrks[detid] += evWt;
0912               } else {
0913                 SumW.insert(std::pair<unsigned int, myEntry>(detid, myEntry(1, Wi, Fac, Fac2)));
0914                 nTrks.insert(std::pair<unsigned int, unsigned int>(detid, evWt));
0915               }
0916             }
0917           }
0918         }
0919       }
0920     }
0921   }
0922   if (debug)
0923     std::cout << "# of Good Tracks " << ntkgood << " out of " << nentries << std::endl;
0924   if (loop == 0) {
0925     h_pbyE->Write("h_pbyE");
0926     h_Ebyp_bfr->Write("h_Ebyp_bfr");
0927   }
0928   if (last) {
0929     h_Ebyp_aftr->Write("h_Ebyp_aftr");
0930   }
0931 
0932   std::map<unsigned int, std::pair<double, double> > cfactors;
0933   unsigned int kount(0), kountus(0);
0934   double sumfactor(0);
0935   for (std::map<unsigned int, TH1D *>::const_iterator itr = histos_.begin(); itr != histos_.end(); ++itr) {
0936     if (writeHisto) {
0937       //      std::pair<double, double> result_write = fitMean(itr->second, 0);

0938       (itr->second)->Write();
0939     }
0940     // The masks are defined in DataFormats/HcalDetId/interface/HcalDetId.h

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

1124       bool selRun = (includeRun_ ? ((t_Run >= runlo_) && (t_Run <= runhi_)) : ((t_Run < runlo_) || (t_Run > runhi_)));
1125       bool selTrack = ((ietaTrack <= 0) || (abs(t_ieta) <= ietaTrack));
1126       if (selRun && (t_nVtx >= nvxlo_) && (t_nVtx <= nvxhi_) && selTrack) {
1127         bool isItRBX(false);
1128         if (cSelect_ != nullptr) {
1129           bool temp = cSelect_->isItRBX(t_DetIds);
1130           if (exclude_)
1131             isItRBX = temp;
1132           else
1133             isItRBX = !(temp);
1134         }
1135         if ((cDuplicate_ != nullptr) && (cDuplicate_->doCorr(2)) && (!isItRBX))
1136           isItRBX = (cDuplicate_->select(t_ieta, t_iphi));
1137         ++kprint;
1138         if (!(isItRBX)) {
1139           for (unsigned int idet = 0; idet < (*t_DetIds).size(); idet++) {
1140             if (selectPhi((*t_DetIds)[idet])) {
1141               unsigned int detid = truncateId((*t_DetIds)[idet], truncateFlag_, debug);
1142               if (debug && (kprint <= 10)) {
1143                 std::cout << "DetId[" << idet << "] Original " << std::hex << (*t_DetIds)[idet] << " truncated "
1144                           << detid << std::dec;
1145               }
1146               if (std::find(detIds_.begin(), detIds_.end(), detid) == detIds_.end()) {
1147                 detIds_.push_back(detid);
1148                 if (debug && (kprint <= 10))
1149                   std::cout << " new";
1150               }
1151               if (debug && (kprint <= 10))
1152                 std::cout << std::endl;
1153             }
1154           }
1155           // Also look at the neighbouring cells if available

1156           if (t_DetIds3 != 0) {
1157             for (unsigned int idet = 0; idet < (*t_DetIds3).size(); idet++) {
1158               if (selectPhi((*t_DetIds3)[idet])) {
1159                 unsigned int detid = truncateId((*t_DetIds3)[idet], truncateFlag_, debug);
1160                 if (std::find(detIds_.begin(), detIds_.end(), detid) == detIds_.end()) {
1161                   detIds_.push_back(detid);
1162                 }
1163               }
1164             }
1165           }
1166         }
1167       }
1168     }
1169   }
1170   if (debug) {
1171     std::cout << "Total of " << detIds_.size() << " detIds" << std::endl;
1172     // The masks are defined in DataFormats/HcalDetId/interface/HcalDetId.h

1173     for (unsigned int k = 0; k < detIds_.size(); ++k) {
1174       int subdet, depth, zside, ieta, iphi;
1175       unpackDetId(detIds_[k], subdet, zside, ieta, iphi, depth);
1176       std::cout << "DetId[" << k << "] " << subdet << ":" << zside * ieta << ":" << depth << ":" << iphi << "  "
1177                 << std::hex << detIds_[k] << std::dec << std::endl;
1178     }
1179   }
1180 }
1181 
1182 void CalibTree::bookHistos(int loop, bool debug) {
1183   unsigned int k(0);
1184   for (std::map<unsigned int, TH1D *>::const_iterator itr = histos_.begin(); itr != histos_.end(); ++itr, ++k) {
1185     if (debug) {
1186       std::cout << "histos[" << k << "] " << std::hex << itr->first << std::dec << " " << itr->second;
1187       if (itr->second != 0)
1188         std::cout << " " << itr->second->GetTitle();
1189       std::cout << std::endl;
1190     }
1191     if (itr->second != 0)
1192       itr->second->Delete();
1193   }
1194 
1195   for (unsigned int k = 0; k < detIds_.size(); ++k) {
1196     char name[20], title[100];
1197     sprintf(name, "Hist%d_%d", detIds_[k], loop);
1198     int subdet, depth, zside, ieta, iphi;
1199     unpackDetId(detIds_[k], subdet, zside, ieta, iphi, depth);
1200     sprintf(title, "Correction for Subdet %d #eta %d depth %d (Loop %d)", subdet, zside * ieta, depth, loop);
1201     TH1D *hist = new TH1D(name, title, 100, 0.0, 5.0);
1202     hist->Sumw2();
1203     if (debug)
1204       std::cout << "Book Histo " << k << " " << title << std::endl;
1205     histos_[detIds_[k]] = hist;
1206   }
1207   std::cout << "Total of " << detIds_.size() << " detIds and " << histos_.size() << std::endl;
1208 }
1209 
1210 bool CalibTree::goodTrack() {
1211   bool ok(true);
1212   double cut(2.0);
1213   double pmom = (useGen_ && (t_gentrackP > 0)) ? t_gentrackP : t_p;
1214   if (sysmode_ == 1) {
1215     ok = ((t_qltyFlag) && (t_hmaxNearP < cut) && (t_eMipDR < 1.0) && (t_mindR1 > 1.0) && (pmom > pmin_) &&
1216           (pmom < pmax_));
1217   } else if (sysmode_ == 2) {
1218     ok = ((t_qltyFlag) && (t_qltyPVFlag) && (t_hmaxNearP < cut) && (t_eMipDR < 1.0) && (t_mindR1 > 1.0) &&
1219           (pmom > pmin_) && (pmom < pmax_));
1220   } else if (sysmode_ == 3) {
1221     ok = ((t_selectTk) && (t_hmaxNearP < cut) && (t_eMipDR < 1.0) && (t_mindR1 > 1.0) && (pmom > pmin_) &&
1222           (pmom < pmax_));
1223   } else if (sysmode_ == 4) {
1224     ok = ((t_selectTk) && (t_qltyMissFlag) && (t_hmaxNearP < 0.0) && (t_eMipDR < 1.0) && (t_mindR1 > 1.0) &&
1225           (pmom > pmin_) && (pmom < pmax_));
1226   } else if (sysmode_ == 5) {
1227     ok = ((t_selectTk) && (t_qltyMissFlag) && (t_hmaxNearP < cut) && (t_eMipDR < 0.5) && (t_mindR1 > 1.0) &&
1228           (pmom > pmin_) && (pmom < pmax_));
1229   } else if (sysmode_ == 6) {
1230     ok = ((t_selectTk) && (t_qltyMissFlag) && (t_hmaxNearP < cut) && (t_eMipDR < 2.0) && (t_mindR1 > 1.0) &&
1231           (pmom > pmin_) && (pmom < pmax_));
1232   } else if (sysmode_ == 7) {
1233     ok = ((t_selectTk) && (t_qltyMissFlag) && (t_hmaxNearP < cut) && (t_eMipDR < 1.0) && (t_mindR1 > 0.5) &&
1234           (pmom > pmin_) && (pmom < pmax_));
1235   } else {
1236     if (sysmode_ < 0) {
1237       double eta = (t_ieta > 0) ? t_ieta : -t_ieta;
1238       if (sysmode_ == -2)
1239         cut = 8.0 * exp(eta * log2by18_);
1240       else
1241         cut = 10.0;
1242     }
1243     ok = ((t_selectTk) && (t_qltyMissFlag) && (t_hmaxNearP < cut) && (t_eMipDR < 1.0) && (t_mindR1 > 1.0) &&
1244           (pmom > pmin_) && (pmom < pmax_));
1245   }
1246   return ok;
1247 }
1248 
1249 void CalibTree::writeCorrFactor(const char *corrFileName, int ietaMax) {
1250   std::ofstream myfile;
1251   myfile.open(corrFileName);
1252   if (!myfile.is_open()) {
1253     std::cout << "** ERROR: Can't open '" << corrFileName << std::endl;
1254   } else {
1255     myfile << "#" << std::setprecision(4) << std::setw(10) << "detId" << std::setw(10) << "ieta" << std::setw(10)
1256            << "depth" << std::setw(15) << "corrFactor" << std::endl;
1257     std::map<unsigned int, std::pair<double, double> >::const_iterator itr;
1258     for (itr = Cprev.begin(); itr != Cprev.end(); ++itr) {
1259       unsigned int detId = itr->first;
1260       int subdet, depth, zside, ieta, iphi;
1261       unpackDetId(detId, subdet, zside, ieta, iphi, depth);
1262       if (ieta <= ietaMax) {
1263         double corrf = ((itr->second.first > 0.1) && (itr->second.first < 4.0)) ? itr->second.first : 1.0;
1264         double dcorr = ((itr->second.first > 0.1) && (itr->second.first < 4.0)) ? itr->second.second : 0.0;
1265         myfile << std::setw(10) << std::hex << detId << std::setw(10) << std::dec << zside * ieta << std::setw(10)
1266                << depth << std::setw(10) << corrf << " " << std::setw(10) << dcorr << std::endl;
1267         std::cout << corrf << ",";
1268       }
1269     }
1270     myfile.close();
1271     std::cout << std::endl;
1272   }
1273 }
1274 
1275 bool CalibTree::selectPhi(unsigned int detId) {
1276   bool flag(true);
1277   // The masks are defined in DataFormats/HcalDetId/interface/HcalDetId.h

1278   if ((phimin_ > 1) || (phimax_ < 72) || (zside_ != 0)) {
1279     int subdet, depth, zside, ieta, iphi;
1280     unpackDetId(detId, subdet, zside, ieta, iphi, depth);
1281     if (phimin_ > 1 || phimax_ < 72) {
1282       if ((iphi < phimin_) || (iphi > phimax_))
1283         flag = false;
1284     }
1285     if (zside_ != 0) {
1286       if (zside != zside_)
1287         flag = false;
1288     }
1289   }
1290   return flag;
1291 }
1292 
1293 std::pair<double, double> CalibTree::fitMean(TH1D *hist, int mode) {
1294   std::pair<double, double> results = std::pair<double, double>(1, 0);
1295   if (hist != 0) {
1296     double mean = hist->GetMean(), rms = hist->GetRMS();
1297     double LowEdge(0.7), HighEdge(1.3);
1298     char option[20];
1299     if (mode == 1) {
1300       LowEdge = mean - 1.5 * rms;
1301       HighEdge = mean + 1.5 * rms;
1302       int nbin = hist->GetNbinsX();
1303       if (LowEdge < hist->GetBinLowEdge(1))
1304         LowEdge = hist->GetBinLowEdge(1);
1305       if (HighEdge > hist->GetBinLowEdge(nbin) + hist->GetBinWidth(nbin))
1306         HighEdge = hist->GetBinLowEdge(nbin) + hist->GetBinWidth(nbin);
1307     }
1308     if (hist->GetEntries() > 100)
1309       sprintf(option, "+QRLS");
1310     else
1311       sprintf(option, "+QRWLS");
1312     double value(mean);
1313     double error = rms / sqrt(hist->GetEntries());
1314     if (hist->GetEntries() > 20) {
1315       TFitResultPtr Fit = hist->Fit("gaus", option, "", LowEdge, HighEdge);
1316       value = Fit->Value(1);
1317       error = Fit->FitResult::Error(1);
1318       /*

1319       LowEdge  = value - 1.5*error;

1320       HighEdge = value + 1.5*error;

1321       value    = Fit->Value(1);

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

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

1324       */
1325     }
1326     results = std::pair<double, double>(value, error);
1327   }
1328   return results;
1329 }
1330 
1331 void CalibTree::makeplots(
1332     double rmin, double rmax, int ietaMax, bool useweight, double fraction, bool debug, Long64_t nmax) {
1333   if (fChain == 0)
1334     return;
1335   Long64_t nentryTot = fChain->GetEntriesFast();
1336   Long64_t nentries = (fraction > 0.01 && fraction < 0.99) ? (Long64_t)(fraction * nentryTot) : nentryTot;
1337   int32_t oddEven = (nmax == -2) ? 1 : ((nmax == -3) ? -1 : 0);
1338   if ((nentries > nmax) && (nmax > 0))
1339     nentries = nmax;
1340 
1341   // Book the histograms

1342   std::map<int, std::pair<TH1D *, TH1D *> > histos;
1343   for (int ieta = -ietaMax; ieta <= ietaMax; ++ieta) {
1344     char name[20], title[100];
1345     sprintf(name, "begin%d", ieta);
1346     if (ieta == 0)
1347       sprintf(title, "Ratio at start");
1348     else
1349       sprintf(title, "Ratio at start for i#eta=%d", ieta);
1350     TH1D *h1 = new TH1D(name, title, 50, rmin, rmax);
1351     h1->Sumw2();
1352     sprintf(name, "end%d", ieta);
1353     if (ieta == 0)
1354       sprintf(title, "Ratio at the end");
1355     else
1356       sprintf(title, "Ratio at the end for i#eta=%d", ieta);
1357     TH1D *h2 = new TH1D(name, title, 50, rmin, rmax);
1358     h2->Sumw2();
1359     histos[ieta] = std::pair<TH1D *, TH1D *>(h1, h2);
1360   }
1361   //Fill the histograms

1362   Long64_t nbytes(0), nb(0);
1363   for (Long64_t jentry = 0; jentry < nentries; jentry++) {
1364     Long64_t ientry = LoadTree(jentry);
1365     nb = fChain->GetEntry(jentry);
1366     nbytes += nb;
1367     if (ientry < 0)
1368       break;
1369     if (oddEven != 0) {
1370       if ((oddEven < 0) && (jentry % 2 == 0))
1371         continue;
1372       else if ((oddEven > 0) && (jentry % 2 != 0))
1373         continue;
1374     }
1375     bool select = ((cDuplicate_ != nullptr) && (cDuplicate_->doCorr(0))) ? (cDuplicate_->isDuplicate(jentry)) : true;
1376     if (!select)
1377       continue;
1378     if ((cDuplicate_ != nullptr) && (cDuplicate_->doCorr(2))) {
1379       select = !(cDuplicate_->select(t_ieta, t_iphi));
1380       if (!select)
1381         continue;
1382     }
1383     if (goodTrack()) {
1384       double pmom = (useGen_ && (t_gentrackP > 0)) ? t_gentrackP : t_p;
1385       CalibTree::energyCalor en1 = energyHcal(pmom, jentry, false);
1386       CalibTree::energyCalor en2 = energyHcal(pmom, jentry, true);
1387       if ((en1.ehcal > 0.001) && (en2.ehcal > 0.001)) {
1388         double evWt = (useweight) ? t_EventWeight : 1.0;
1389         double ratioi = en1.ehcal / (pmom - t_eMipDR);
1390         double ratiof = en2.ehcal / (pmom - t_eMipDR);
1391         if (t_ieta >= -ietaMax && t_ieta <= ietaMax && t_ieta != 0) {
1392           if (ratioi >= rmin && ratioi <= rmax) {
1393             histos[0].first->Fill(ratioi, evWt);
1394             histos[t_ieta].first->Fill(ratioi, evWt);
1395           }
1396           if (ratiof >= rmin && ratiof <= rmax) {
1397             histos[0].second->Fill(ratiof, evWt);
1398             histos[t_ieta].second->Fill(ratiof, evWt);
1399           }
1400         }
1401       }
1402     }
1403   }
1404 
1405   //Fit the histograms

1406   TH1D *hbef1 = new TH1D("Eta1Bf", "Mean vs i#eta", 2 * ietaMax, -ietaMax, ietaMax);
1407   TH1D *hbef2 = new TH1D("Eta2Bf", "Median vs i#eta", 2 * ietaMax, -ietaMax, ietaMax);
1408   TH1D *haft1 = new TH1D("Eta1Af", "Mean vs i#eta", 2 * ietaMax, -ietaMax, ietaMax);
1409   TH1D *haft2 = new TH1D("Eta2Af", "Median vs i#eta", 2 * ietaMax, -ietaMax, ietaMax);
1410   for (int ieta = -ietaMax; ieta <= ietaMax; ++ieta) {
1411     int bin = (ieta < 0) ? (ieta + ietaMax + 1) : (ieta + ietaMax);
1412     TH1D *h1 = histos[ieta].first;
1413     double mean1 = h1->GetMean();
1414     double err1 = h1->GetMeanError();
1415     std::pair<double, double> fit1 = fitMean(h1, 1);
1416     if (debug) {
1417       std::cout << ieta << " " << h1->GetName() << " " << mean1 << " +- " << err1 << " and " << fit1.first << " +- "
1418                 << fit1.second << std::endl;
1419     }
1420     if (ieta != 0) {
1421       hbef1->SetBinContent(bin, mean1);
1422       hbef1->SetBinError(bin, err1);
1423       hbef2->SetBinContent(bin, fit1.first);
1424       hbef2->SetBinError(bin, fit1.second);
1425     }
1426     h1->Write();
1427     TH1D *h2 = histos[ieta].second;
1428     double mean2 = h2->GetMean();
1429     double err2 = h2->GetMeanError();
1430     std::pair<double, double> fit2 = fitMean(h2, 1);
1431     if (debug) {
1432       std::cout << ieta << " " << h2->GetName() << " " << mean2 << " +- " << err2 << " and " << fit2.first << " +- "
1433                 << fit2.second << std::endl;
1434     }
1435     if (ieta != 0) {
1436       haft1->SetBinContent(bin, mean2);
1437       haft1->SetBinError(bin, err2);
1438       haft2->SetBinContent(bin, fit2.first);
1439       haft2->SetBinError(bin, fit2.second);
1440     }
1441     h2->Write();
1442   }
1443   fitPol0(hbef1, debug);
1444   fitPol0(hbef2, debug);
1445   fitPol0(haft1, debug);
1446   fitPol0(haft2, debug);
1447 }
1448 
1449 void CalibTree::fitPol0(TH1D *hist, bool debug) {
1450   hist->GetXaxis()->SetTitle("i#eta");
1451   hist->GetYaxis()->SetTitle("<E_{HCAL}/(p-E_{ECAL})>");
1452   hist->GetYaxis()->SetRangeUser(0.4, 1.6);
1453   TFitResultPtr Fit = hist->Fit("pol0", "+QRWLS");
1454   if (debug) {
1455     std::cout << "Fit to Pol0 to " << hist->GetTitle() << ": " << Fit->Value(0) << " +- " << Fit->FitResult::Error(0)
1456               << std::endl;
1457   }
1458   hist->Write();
1459 }
1460 
1461 void CalibTree::highEtaFactors(int ietaMax, bool debug) {
1462   std::map<unsigned int, std::pair<double, double> >::const_iterator itr;
1463   std::pair<double, double> cfacp, cfacn;
1464   cfacp = cfacn = std::pair<double, double>(1.0, 0.0);
1465   for (itr = Cprev.begin(); itr != Cprev.end(); ++itr) {
1466     unsigned int detid = itr->first;
1467     int subdet, depth, zside, ieta, iphi;
1468     unpackDetId(detid, subdet, zside, ieta, iphi, depth);
1469     if ((ieta == ietaMax) && (depth == 1)) {
1470       if (zside > 0)
1471         cfacp = itr->second;
1472       else
1473         cfacn = itr->second;
1474     }
1475   }
1476   if (debug) {
1477     std::cout << "Correction factor for (" << ietaMax << ",1) = " << cfacp.first << ":" << cfacp.second << " and ("
1478               << -ietaMax << ",1) = " << cfacn.first << ":" << cfacn.second << std::endl;
1479   }
1480   for (itr = Cprev.begin(); itr != Cprev.end(); ++itr) {
1481     unsigned int detid = itr->first;
1482     int subdet, depth, zside, ieta, iphi;
1483     unpackDetId(detid, subdet, zside, ieta, iphi, depth);
1484     if (ieta > ietaMax) {
1485       Cprev[detid] = (zside > 0) ? cfacp : cfacn;
1486       if (debug) {
1487         std::cout << "Set correction factor for (" << zside * ieta << "," << depth << ") = " << Cprev[detid].first
1488                   << ":" << Cprev[detid].second << std::endl;
1489       }
1490     }
1491   }
1492 }
1493 
1494 CalibTree::energyCalor CalibTree::energyHcal(double pmom, const Long64_t &entry, bool final) {
1495   double etot = t_eHcal;
1496   double etot2 = t_eHcal;
1497   double ediff = (t_eHcal30 - t_eHcal10);
1498   if (final) {
1499     etot = etot2 = 0;
1500     for (unsigned int idet = 0; idet < (*t_DetIds).size(); idet++) {
1501       // Apply thresholds if necessary

1502       bool okcell = (thrForm_ == 0) || ((*t_HitEnergies)[idet] > (cThr_->threshold((*t_DetIds)[idet])));
1503       if (okcell && selectPhi((*t_DetIds)[idet])) {
1504         unsigned int id = (*t_DetIds)[idet];
1505         double hitEn(0);
1506         unsigned int detid = truncateId(id, truncateFlag_, false);
1507         if (Cprev.find(detid) != Cprev.end())
1508           hitEn = Cprev[detid].first * (*t_HitEnergies)[idet];
1509         else
1510           hitEn = (*t_HitEnergies)[idet];
1511         if ((rcorForm_ != 3) && (rcorForm_ >= 0) && (cFactor_))
1512           hitEn *= cFactor_->getCorr(t_Run, id);
1513         if ((cDuplicate_ != nullptr) && (cDuplicate_->doCorr(1)))
1514           hitEn *= cDuplicate_->getWeight(id);
1515         if ((cDuplicate_ != nullptr) && (cDuplicate_->doCorr(3))) {
1516           int subdet, zside, ieta, iphi, depth;
1517           unpackDetId((*t_DetIds)[idet], subdet, zside, ieta, iphi, depth);
1518           hitEn *= cDuplicate_->getCorr(t_Run, ieta, depth);
1519         }
1520         etot += hitEn;
1521         etot2 += ((*t_HitEnergies)[idet]);
1522       }
1523     }
1524     // Now the outer cone

1525     double etot1(0), etot3(0);
1526     if (t_DetIds1 != 0 && t_DetIds3 != 0) {
1527       for (unsigned int idet = 0; idet < (*t_DetIds1).size(); idet++) {
1528         // Apply thresholds if necessary

1529         bool okcell = (thrForm_ == 0) || ((*t_HitEnergies1)[idet] > (cThr_->threshold((*t_DetIds1)[idet])));
1530         if (okcell && selectPhi((*t_DetIds1)[idet])) {
1531           unsigned int id = (*t_DetIds1)[idet];
1532           unsigned int detid = truncateId(id, truncateFlag_, false);
1533           double hitEn(0);
1534           if (Cprev.find(detid) != Cprev.end())
1535             hitEn = Cprev[detid].first * (*t_HitEnergies1)[idet];
1536           else
1537             hitEn = (*t_HitEnergies1)[idet];
1538           if ((rcorForm_ != 3) && (rcorForm_ >= 0) && (cFactor_))
1539             hitEn *= cFactor_->getCorr(t_Run, id);
1540           if ((cDuplicate_ != nullptr) && (cDuplicate_->doCorr(1)))
1541             hitEn *= cDuplicate_->getWeight(id);
1542           if ((cDuplicate_ != nullptr) && (cDuplicate_->doCorr(3))) {
1543             int subdet, zside, ieta, iphi, depth;
1544             unpackDetId((*t_DetIds1)[idet], subdet, zside, ieta, iphi, depth);
1545             hitEn *= cDuplicate_->getCorr(t_Run, ieta, depth);
1546           }
1547           etot1 += hitEn;
1548         }
1549       }
1550       for (unsigned int idet = 0; idet < (*t_DetIds3).size(); idet++) {
1551         // Apply thresholds if necessary

1552         bool okcell = (thrForm_ == 0) || ((*t_HitEnergies3)[idet] > (cThr_->threshold((*t_DetIds3)[idet])));
1553         if (okcell && selectPhi((*t_DetIds3)[idet])) {
1554           unsigned int id = (*t_DetIds3)[idet];
1555           unsigned int detid = truncateId(id, truncateFlag_, false);
1556           double hitEn(0);
1557           if (Cprev.find(detid) != Cprev.end())
1558             hitEn = Cprev[detid].first * (*t_HitEnergies3)[idet];
1559           else
1560             hitEn = (*t_HitEnergies3)[idet];
1561           if ((rcorForm_ != 3) && (rcorForm_ >= 0) && (cFactor_))
1562             hitEn *= cFactor_->getCorr(t_Run, id);
1563           if ((cDuplicate_ != nullptr) && (cDuplicate_->doCorr(1)))
1564             hitEn *= cDuplicate_->getWeight(id);
1565           if ((cDuplicate_ != nullptr) && (cDuplicate_->doCorr(3))) {
1566             int subdet, zside, ieta, iphi, depth;
1567             unpackDetId((*t_DetIds3)[idet], subdet, zside, ieta, iphi, depth);
1568             hitEn *= cDuplicate_->getCorr(t_Run, ieta, depth);
1569           }
1570           etot3 += hitEn;
1571         }
1572       }
1573     }
1574     ediff = etot3 - etot1;
1575   }
1576   // PU correction only for loose isolation cut

1577   double ehcal = (((rcorForm_ == 3) && (cFactor_ != nullptr))
1578                       ? (etot * cFactor_->getCorr(entry))
1579                       : ((puCorr_ == 0) ? etot
1580                                         : ((puCorr_ < 0) ? (etot * puFactor(-puCorr_, t_ieta, pmom, etot, ediff))
1581                                                          : puFactorRho(puCorr_, t_ieta, t_rhoh, etot))));
1582   return CalibTree::energyCalor(etot, etot2, ehcal);
1583 }
[ Source navigation ] [ Diff markup ] [ Identifier search ] [ general search ]

This page was automatically generated by the 2.2.1 LXR engine.

The LXR team
Valid CSS 2.1! Valid HTML 4.01!