Back to home page

Project CMSSW displayed by LXR
 
 

    

File indexing completed on 2021-04-09 02:50:18 

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

0002 // Usage:

0003 // .L CalibTree.C+g

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

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

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

0007 //      maxIter, corForm, useGen, runlo, runhi, phimin, phimax, zside, nvxlo,

0008 //      nvxhi, rbx, exclude, higheta, fraction, writeDebugHisto, debug);

0009 //

0010 //  where:

0011 //

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

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

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

0015 //  dirName     (std::string) = name of the directory where the Tree resides

0016 //                              ("HcalIsoTrkAnalyzer")

0017 //  treeName    (std::string) = name of the Tree ("CalibTree")

0018 //  outFileName (std::string) = name of the output ROOT file

0019 //                              ("Silver_out.root")

0020 //  corrFileName(std::string) = name of the output text file with correction

0021 //                              factors ("Silver_corr.txt")

0022 //  dupFileName (std::string) = name of the file containing list of sequence

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

0024 //  rcorFileName (std::string)= name of the text file having the correction

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

0026 //                              to be used for raddam/depth dependent

0027 //                              correction  (default="", no corr.)

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

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

0030 //                              (false -- use MPV)

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

0032 //                              used in evaluating convergence criterion (0)

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

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

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

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

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

0038 //                              factor is to be determined (25)

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

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

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

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

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

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

0045 //                              correction (-5: 2018 version)

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

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

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

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

0050 //  truncateFlag    (int)     = Flag to treat different depths differently (0)

0051 //                              both depths of ieta 15, 16 of HB as depth 1 (1)

0052 //                              all depths as depth 1 (2), all depths in HE

0053 //                              with values > 1 as depth 2 (3), all depths in

0054 //                              HB with values > 1 as depth 2 (4), all depths

0055 //                              in HB and HE with values > 1 as depth 2 (5)

0056 //                              (Default 0)

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

0058 //  rcorForm        (int)     = type of rcorFileName: (0) for Raddam correction,

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

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

0061 //                              method for pileup correction. (Default 0)

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

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

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

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

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

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

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

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

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

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

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

0073 //  rbx             (int)     = zside*(Subdet*100+RBX #) to be consdered (0)

0074 //                              HEP17 should correspond to 217

0075 //  exclude         (bool)    = RBX specified by *rbx* to be exluded or only

0076 //                              considered (false)

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

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

0079 //                              and depth 1 (1) (default 1)

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

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

0082 //                              in o/p file (false)

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

0084 //                              (false)

0085 //

0086 //  doIt(inFileName, dupFileName)

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

0088 ////////////////////////////////////////////////////////////////////////////////

0089 
0090 #include <TStyle.h>
0091 #include <TCanvas.h>
0092 #include <TROOT.h>
0093 #include <TChain.h>
0094 #include <TFile.h>
0095 #include <TFitResult.h>
0096 #include <TFitResultPtr.h>
0097 #include <TTree.h>
0098 #include <TH1.h>
0099 #include <TGraph.h>
0100 #include <TProfile.h>
0101 #include <algorithm>
0102 #include <vector>
0103 #include <string>
0104 #include <iomanip>
0105 #include <iostream>
0106 #include <fstream>
0107 #include <sstream>
0108 
0109 #include "CalibCorr.C"
0110 
0111 void Run(const char *inFileName = "Silver.root",
0112          const char *dirName = "HcalIsoTrkAnalyzer",
0113          const char *treeName = "CalibTree",
0114          const char *outFileName = "Silver_out.root",
0115          const char *corrFileName = "Silver_corr.txt",
0116          const char *dupFileName = "events_DXS2.txt",
0117          const char *rcorFileName = "",
0118          bool useweight = true,
0119          bool useMean = false,
0120          int nMin = 0,
0121          bool inverse = true,
0122          double ratMin = 0.25,
0123          double ratMax = 3.,
0124          int ietaMax = 25,
0125          int ietaTrack = -1,
0126          int sysmode = -1,
0127          int puCorr = -5,
0128          int applyL1Cut = 1,
0129          double l1Cut = 0.5,
0130          int truncateFlag = 0,
0131          int maxIter = 30,
0132          int rcorForm = 0,
0133          bool useGen = false,
0134          int runlo = 0,
0135          int runhi = 99999999,
0136          int phimin = 1,
0137          int phimax = 72,
0138          int zside = 0,
0139          int nvxlo = 0,
0140          int nvxhi = 1000,
0141          int rbx = 0,
0142          bool exclude = true,
0143          int higheta = 1,
0144          double fraction = 1.0,
0145          bool writeDebugHisto = false,
0146          bool debug = false);
0147 
0148 // Fixed size dimensions of array or collections stored in the TTree if any.

0149 
0150 class CalibTree {
0151 public:
0152   struct myEntry {
0153     myEntry(int k = 0, double f0 = 0, double f1 = 0, double f2 = 0) : kount(k), fact0(f0), fact1(f1), fact2(f2) {}
0154     int kount;
0155     double fact0, fact1, fact2;
0156   };
0157 
0158   struct energyCalor {
0159     energyCalor(double e1 = 0, double e2 = 0, double e3 = 0) : Etot(e1), Etot2(e2), ehcal(e3) {}
0160     double Etot, Etot2, ehcal;
0161   };
0162 
0163   CalibTree(const char *dupFileName,
0164             const char *rcorFileName,
0165             int truncateFlag,
0166             bool useMean,
0167             int runlo,
0168             int runhi,
0169             int phimin,
0170             int phimax,
0171             int zside,
0172             int nvxlo,
0173             int nvxhi,
0174             int sysmode,
0175             int rbx,
0176             int puCorr,
0177             int rcorForm,
0178             bool useGen,
0179             bool exclude,
0180             int higheta,
0181             TChain *tree);
0182   virtual ~CalibTree();
0183   virtual Int_t Cut(Long64_t entry);
0184   virtual Int_t GetEntry(Long64_t entry);
0185   virtual Long64_t LoadTree(Long64_t entry);
0186   virtual void Init(TChain *tree, const char *dupFileName);
0187   virtual Double_t Loop(int k,
0188                         TFile *fout,
0189                         bool useweight,
0190                         int nMin,
0191                         bool inverse,
0192                         double rMin,
0193                         double rMax,
0194                         int ietaMax,
0195                         int ietaTrack,
0196                         int applyL1Cut,
0197                         double l1Cut,
0198                         bool last,
0199                         double fraction,
0200                         bool writeHisto,
0201                         bool debug);
0202   virtual Bool_t Notify();
0203   virtual void Show(Long64_t entry = -1);
0204   bool goodTrack();
0205   void writeCorrFactor(const char *corrFileName, int ietaMax);
0206   bool selectPhi(unsigned int detId);
0207   std::pair<double, double> fitMean(TH1D *, int);
0208   void makeplots(double rmin, double rmax, int ietaMax, bool useWeight, double fraction, bool debug);
0209   void fitPol0(TH1D *hist, bool debug);
0210   void highEtaFactors(int ietaMax, bool debug);
0211   energyCalor energyHcal(double pmom, const Long64_t &entry, bool final);
0212 
0213   TChain *fChain;  //!pointer to the analyzed TTree or TChain

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

0215   TH1D *h_pbyE, *h_cvg;
0216   TProfile *h_Ebyp_bfr, *h_Ebyp_aftr;
0217 
0218 private:
0219   // Declaration of leaf types

0220   Int_t t_Run;
0221   Int_t t_Event;
0222   Int_t t_DataType;
0223   Int_t t_ieta;
0224   Int_t t_iphi;
0225   Double_t t_EventWeight;
0226   Int_t t_nVtx;
0227   Int_t t_nTrk;
0228   Int_t t_goodPV;
0229   Double_t t_l1pt;
0230   Double_t t_l1eta;
0231   Double_t t_l1phi;
0232   Double_t t_l3pt;
0233   Double_t t_l3eta;
0234   Double_t t_l3phi;
0235   Double_t t_p;
0236   Double_t t_pt;
0237   Double_t t_phi;
0238   Double_t t_mindR1;
0239   Double_t t_mindR2;
0240   Double_t t_eMipDR;
0241   Double_t t_eHcal;
0242   Double_t t_eHcal10;
0243   Double_t t_eHcal30;
0244   Double_t t_hmaxNearP;
0245   Double_t t_rhoh;
0246   Bool_t t_selectTk;
0247   Bool_t t_qltyFlag;
0248   Bool_t t_qltyMissFlag;
0249   Bool_t t_qltyPVFlag;
0250   Double_t t_gentrackP;
0251   std::vector<unsigned int> *t_DetIds;
0252   std::vector<double> *t_HitEnergies;
0253   std::vector<bool> *t_trgbits;
0254   std::vector<unsigned int> *t_DetIds1;
0255   std::vector<unsigned int> *t_DetIds3;
0256   std::vector<double> *t_HitEnergies1;
0257   std::vector<double> *t_HitEnergies3;
0258 
0259   // List of branches

0260   TBranch *b_t_Run;           //!

0261   TBranch *b_t_Event;         //!

0262   TBranch *b_t_DataType;      //!

0263   TBranch *b_t_ieta;          //!

0264   TBranch *b_t_iphi;          //!

0265   TBranch *b_t_EventWeight;   //!

0266   TBranch *b_t_nVtx;          //!

0267   TBranch *b_t_nTrk;          //!

0268   TBranch *b_t_goodPV;        //!

0269   TBranch *b_t_l1pt;          //!

0270   TBranch *b_t_l1eta;         //!

0271   TBranch *b_t_l1phi;         //!

0272   TBranch *b_t_l3pt;          //!

0273   TBranch *b_t_l3eta;         //!

0274   TBranch *b_t_l3phi;         //!

0275   TBranch *b_t_p;             //!

0276   TBranch *b_t_pt;            //!

0277   TBranch *b_t_phi;           //!

0278   TBranch *b_t_mindR1;        //!

0279   TBranch *b_t_mindR2;        //!

0280   TBranch *b_t_eMipDR;        //!

0281   TBranch *b_t_eHcal;         //!

0282   TBranch *b_t_eHcal10;       //!

0283   TBranch *b_t_eHcal30;       //!

0284   TBranch *b_t_hmaxNearP;     //!

0285   TBranch *b_t_rhoh;          //!

0286   TBranch *b_t_selectTk;      //!

0287   TBranch *b_t_qltyFlag;      //!

0288   TBranch *b_t_qltyMissFlag;  //!

0289   TBranch *b_t_qltyPVFlag;    //!

0290   TBranch *b_t_gentrackP;     //!

0291   TBranch *b_t_DetIds;        //!

0292   TBranch *b_t_HitEnergies;   //!

0293   TBranch *b_t_trgbits;       //!

0294   TBranch *b_t_DetIds1;       //!

0295   TBranch *b_t_DetIds3;       //!

0296   TBranch *b_t_HitEnergies1;  //!

0297   TBranch *b_t_HitEnergies3;  //!

0298 
0299   CalibCorr *cFactor_;
0300   CalibSelectRBX *cSelect_;
0301   const int truncateFlag_;
0302   const bool useMean_;
0303   int runlo_, runhi_;
0304   const int phimin_, phimax_;
0305   const int zside_, nvxlo_, nvxhi_;
0306   const int sysmode_, rbx_, puCorr_;
0307   int rcorForm_;
0308   const bool useGen_, exclude_;
0309   const int higheta_;
0310   bool includeRun_;
0311   double log2by18_, eHcalDelta_;
0312   std::vector<Long64_t> entries;
0313   std::vector<unsigned int> detIds;
0314   std::map<unsigned int, TH1D *> histos;
0315   std::map<unsigned int, std::pair<double, double> > Cprev;
0316 };
0317 
0318 void doIt(const char *infile, const char *dup) {
0319   char outf1[100], outf2[100];
0320   double lumt(1.0), fac(0.5);
0321   for (int k = 0; k < 5; ++k) {
0322     sprintf(outf1, "%s_%d.root", infile, k);
0323     sprintf(outf2, "%s_%d.txt", infile, k);
0324     double lumi = (k == 0) ? -1 : lumt;
0325     lumt *= fac;
0326     Run(infile,
0327         "HcalIsoTrkAnalyzer",
0328         "CalibTree",
0329         outf1,
0330         outf2,
0331         dup,
0332         "",
0333         true,
0334         false,
0335         0,
0336         true,
0337         0.25,
0338         3.0,
0339         25,
0340         -1,
0341         -1,
0342         -5,
0343         1,
0344         0.5,
0345         0,
0346         30,
0347         0,
0348         false,
0349         0,
0350         99999999,
0351         1,
0352         72,
0353         0,
0354         1000,
0355         0,
0356         true,
0357         1,
0358         lumi,
0359         false,
0360         false);
0361   }
0362 }
0363 
0364 void Run(const char *inFileName,
0365          const char *dirName,
0366          const char *treeName,
0367          const char *outFileName,
0368          const char *corrFileName,
0369          const char *dupFileName,
0370          const char *rcorFileName,
0371          bool useweight,
0372          bool useMean,
0373          int nMin,
0374          bool inverse,
0375          double ratMin,
0376          double ratMax,
0377          int ietaMax,
0378          int ietaTrack,
0379          int sysmode,
0380          int puCorr,
0381          int applyL1Cut,
0382          double l1Cut,
0383          int truncateFlag,
0384          int maxIter,
0385          int rcorForm,
0386          bool useGen,
0387          int runlo,
0388          int runhi,
0389          int phimin,
0390          int phimax,
0391          int zside,
0392          int nvxlo,
0393          int nvxhi,
0394          int rbx,
0395          bool exclude,
0396          int higheta,
0397          double fraction,
0398          bool writeHisto,
0399          bool debug) {
0400   char name[500];
0401   sprintf(name, "%s/%s", dirName, treeName);
0402   TChain *chain = new TChain(name);
0403   std::cout << "Create a chain for " << name << " from " << inFileName << std::endl;
0404   if (!fillChain(chain, inFileName)) {
0405     std::cout << "*****No valid tree chain can be obtained*****" << std::endl;
0406   } else {
0407     std::cout << "Proceed with a tree chain with " << chain->GetEntries() << " entries" << std::endl;
0408     Long64_t nentryTot = chain->GetEntries();
0409     Long64_t nentries = (fraction > 0.01 && fraction < 0.99) ? (Long64_t)(fraction * nentryTot) : nentryTot;
0410     static const int maxIterMax = 100;
0411     if (maxIter > maxIterMax)
0412       maxIter = maxIterMax;
0413     std::cout << "Tree " << name << " " << chain << " in directory " << dirName << " from file " << inFileName
0414               << " with nentries (tracks): " << nentries << std::endl;
0415     unsigned int k(0), kmax(maxIter);
0416     CalibTree t(dupFileName,
0417                 rcorFileName,
0418                 truncateFlag,
0419                 useMean,
0420                 runlo,
0421                 runhi,
0422                 phimin,
0423                 phimax,
0424                 zside,
0425                 nvxlo,
0426                 nvxhi,
0427                 sysmode,
0428                 rbx,
0429                 puCorr,
0430                 rcorForm,
0431                 useGen,
0432                 exclude,
0433                 higheta,
0434                 chain);
0435     t.h_pbyE = new TH1D("pbyE", "pbyE", 100, -1.0, 9.0);
0436     t.h_Ebyp_bfr = new TProfile("Ebyp_bfr", "Ebyp_bfr", 60, -30, 30, 0, 10);
0437     t.h_Ebyp_aftr = new TProfile("Ebyp_aftr", "Ebyp_aftr", 60, -30, 30, 0, 10);
0438     t.h_cvg = new TH1D("Cvg0", "Convergence", kmax, 0, kmax);
0439     t.h_cvg->SetMarkerStyle(7);
0440     t.h_cvg->SetMarkerSize(5.0);
0441 
0442     TFile *fout = new TFile(outFileName, "RECREATE");
0443     std::cout << "Output file: " << outFileName << " opened in recreate mode" << std::endl;
0444     fout->cd();
0445 
0446     double cvgs[maxIterMax], itrs[maxIterMax];
0447     for (; k <= kmax; ++k) {
0448       std::cout << "Calling Loop() " << k << "th time" << std::endl;
0449       double cvg = t.Loop(k,
0450                           fout,
0451                           useweight,
0452                           nMin,
0453                           inverse,
0454                           ratMin,
0455                           ratMax,
0456                           ietaMax,
0457                           ietaTrack,
0458                           applyL1Cut,
0459                           l1Cut,
0460                           k == kmax,
0461                           fraction,
0462                           writeHisto,
0463                           debug);
0464       itrs[k] = k;
0465       cvgs[k] = cvg;
0466       if (cvg < 0.00001)
0467         break;
0468     }
0469 
0470     t.writeCorrFactor(corrFileName, ietaMax);
0471 
0472     fout->cd();
0473     TGraph *g_cvg;
0474     g_cvg = new TGraph(k, itrs, cvgs);
0475     g_cvg->SetMarkerStyle(7);
0476     g_cvg->SetMarkerSize(5.0);
0477     g_cvg->Draw("AP");
0478     g_cvg->Write("Cvg");
0479     std::cout << "Finish looping after " << k << " iterations" << std::endl;
0480     t.makeplots(ratMin, ratMax, ietaMax, useweight, fraction, debug);
0481     fout->Close();
0482   }
0483 }
0484 
0485 CalibTree::CalibTree(const char *dupFileName,
0486                      const char *rcorFileName,
0487                      int flag,
0488                      bool useMean,
0489                      int runlo,
0490                      int runhi,
0491                      int phimin,
0492                      int phimax,
0493                      int zside,
0494                      int nvxlo,
0495                      int nvxhi,
0496                      int mode,
0497                      int rbx,
0498                      int pu,
0499                      int rForm,
0500                      bool gen,
0501                      bool excl,
0502                      int heta,
0503                      TChain *tree)
0504     : fChain(nullptr),
0505       cFactor_(nullptr),
0506       cSelect_(nullptr),
0507       truncateFlag_(flag),
0508       useMean_(useMean),
0509       runlo_(runlo),
0510       runhi_(runhi),
0511       phimin_(phimin),
0512       phimax_(phimax),
0513       zside_(zside),
0514       nvxlo_(nvxlo),
0515       nvxhi_(nvxhi),
0516       sysmode_(mode),
0517       rbx_(rbx),
0518       puCorr_(pu),
0519       rcorForm_(rForm),
0520       useGen_(gen),
0521       exclude_(excl),
0522       higheta_(heta),
0523       includeRun_(true) {
0524   if (runlo_ < 0 || runhi_ < 0) {
0525     runlo_ = std::abs(runlo_);
0526     runhi_ = std::abs(runhi_);
0527     includeRun_ = false;
0528   }
0529   log2by18_ = std::log(2.5) / 18.0;
0530   eHcalDelta_ = 0;
0531   std::cout << "Initialize CalibTree with TruncateFlag " << truncateFlag_ << " UseMean " << useMean_ << " Run Range "
0532             << runlo_ << ":" << runhi_ << " Phi Range " << phimin_ << ":" << phimax_ << ":" << zside_
0533             << " Vertex Range " << nvxlo_ << ":" << nvxhi_ << " Mode " << sysmode_ << " PU " << puCorr_
0534             << " rcorFormat " << rcorForm_ << " Gen " << useGen_ << " High Eta " << higheta_ << std::endl;
0535   std::cout << "Duplicate events read from " << dupFileName << " RadDam Corrections read from " << rcorFileName
0536             << " Treat RBX " << rbx_ << " with exclusion mode " << exclude_ << std::endl;
0537   Init(tree, dupFileName);
0538   if (std::string(rcorFileName) != "") {
0539     cFactor_ = new CalibCorr(rcorFileName, rcorForm_, false);
0540   } else {
0541     rcorForm_ = 0;
0542   }
0543   if (rbx != 0)
0544     cSelect_ = new CalibSelectRBX(rbx);
0545 }
0546 
0547 CalibTree::~CalibTree() {
0548   delete cFactor_;
0549   delete cSelect_;
0550   if (!fChain)
0551     return;
0552   delete fChain->GetCurrentFile();
0553 }
0554 
0555 Int_t CalibTree::GetEntry(Long64_t entry) {
0556   // Read contents of entry.

0557   if (!fChain)
0558     return 0;
0559   return fChain->GetEntry(entry);
0560 }
0561 
0562 Long64_t CalibTree::LoadTree(Long64_t entry) {
0563   // Set the environment to read one entry

0564   if (!fChain)
0565     return -5;
0566   Long64_t centry = fChain->LoadTree(entry);
0567   if (centry < 0)
0568     return centry;
0569   if (fChain->GetTreeNumber() != fCurrent) {
0570     fCurrent = fChain->GetTreeNumber();
0571     Notify();
0572   }
0573   return centry;
0574 }
0575 
0576 void CalibTree::Init(TChain *tree, const char *dupFileName) {
0577   // The Init() function is called when the selector needs to initialize

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

0579   // pointers of the tree will be set.

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

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

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

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

0584 
0585   // Set object pointer

0586   t_DetIds = 0;
0587   t_HitEnergies = 0;
0588   t_trgbits = 0;
0589   t_DetIds1 = 0;
0590   t_DetIds3 = 0;
0591   t_HitEnergies1 = 0;
0592   t_HitEnergies3 = 0;
0593   // Set branch addresses and branch pointers

0594   fChain = tree;
0595   if (!tree)
0596     return;
0597   fCurrent = -1;
0598   fChain->SetMakeClass(1);
0599 
0600   fChain->SetBranchAddress("t_Run", &t_Run, &b_t_Run);
0601   fChain->SetBranchAddress("t_Event", &t_Event, &b_t_Event);
0602   fChain->SetBranchAddress("t_DataType", &t_DataType, &b_t_DataType);
0603   fChain->SetBranchAddress("t_ieta", &t_ieta, &b_t_ieta);
0604   fChain->SetBranchAddress("t_iphi", &t_iphi, &b_t_iphi);
0605   fChain->SetBranchAddress("t_EventWeight", &t_EventWeight, &b_t_EventWeight);
0606   fChain->SetBranchAddress("t_nVtx", &t_nVtx, &b_t_nVtx);
0607   fChain->SetBranchAddress("t_nTrk", &t_nTrk, &b_t_nTrk);
0608   fChain->SetBranchAddress("t_goodPV", &t_goodPV, &b_t_goodPV);
0609   fChain->SetBranchAddress("t_l1pt", &t_l1pt, &b_t_l1pt);
0610   fChain->SetBranchAddress("t_l1eta", &t_l1eta, &b_t_l1eta);
0611   fChain->SetBranchAddress("t_l1phi", &t_l1phi, &b_t_l1phi);
0612   fChain->SetBranchAddress("t_l3pt", &t_l3pt, &b_t_l3pt);
0613   fChain->SetBranchAddress("t_l3eta", &t_l3eta, &b_t_l3eta);
0614   fChain->SetBranchAddress("t_l3phi", &t_l3phi, &b_t_l3phi);
0615   fChain->SetBranchAddress("t_p", &t_p, &b_t_p);
0616   fChain->SetBranchAddress("t_pt", &t_pt, &b_t_pt);
0617   fChain->SetBranchAddress("t_phi", &t_phi, &b_t_phi);
0618   fChain->SetBranchAddress("t_mindR1", &t_mindR1, &b_t_mindR1);
0619   fChain->SetBranchAddress("t_mindR2", &t_mindR2, &b_t_mindR2);
0620   fChain->SetBranchAddress("t_eMipDR", &t_eMipDR, &b_t_eMipDR);
0621   fChain->SetBranchAddress("t_eHcal", &t_eHcal, &b_t_eHcal);
0622   fChain->SetBranchAddress("t_eHcal10", &t_eHcal10, &b_t_eHcal10);
0623   fChain->SetBranchAddress("t_eHcal30", &t_eHcal30, &b_t_eHcal30);
0624   fChain->SetBranchAddress("t_hmaxNearP", &t_hmaxNearP, &b_t_hmaxNearP);
0625   fChain->SetBranchAddress("t_rhoh", &t_rhoh, &b_t_rhoh);
0626   fChain->SetBranchAddress("t_selectTk", &t_selectTk, &b_t_selectTk);
0627   fChain->SetBranchAddress("t_qltyFlag", &t_qltyFlag, &b_t_qltyFlag);
0628   fChain->SetBranchAddress("t_qltyMissFlag", &t_qltyMissFlag, &b_t_qltyMissFlag);
0629   fChain->SetBranchAddress("t_qltyPVFlag", &t_qltyPVFlag, &b_t_qltyPVFlag);
0630   fChain->SetBranchAddress("t_gentrackP", &t_gentrackP, &b_t_gentrackP);
0631   fChain->SetBranchAddress("t_DetIds", &t_DetIds, &b_t_DetIds);
0632   fChain->SetBranchAddress("t_HitEnergies", &t_HitEnergies, &b_t_HitEnergies);
0633   fChain->SetBranchAddress("t_trgbits", &t_trgbits, &b_t_trgbits);
0634   fChain->SetBranchAddress("t_DetIds1", &t_DetIds1, &b_t_DetIds1);
0635   fChain->SetBranchAddress("t_DetIds3", &t_DetIds3, &b_t_DetIds3);
0636   fChain->SetBranchAddress("t_HitEnergies1", &t_HitEnergies1, &b_t_HitEnergies1);
0637   fChain->SetBranchAddress("t_HitEnergies3", &t_HitEnergies3, &b_t_HitEnergies3);
0638   Notify();
0639 
0640   if (std::string(dupFileName) != "") {
0641     ifstream infil1(dupFileName);
0642     if (!infil1.is_open()) {
0643       std::cout << "Cannot open " << dupFileName << std::endl;
0644     } else {
0645       while (1) {
0646         Long64_t jentry;
0647         infil1 >> jentry;
0648         if (!infil1.good())
0649           break;
0650         entries.push_back(jentry);
0651       }
0652       infil1.close();
0653       std::cout << "Reads a list of " << entries.size() << " events from " << dupFileName << std::endl;
0654     }
0655   } else {
0656     std::cout << "No duplicate events in the input file" << std::endl;
0657   }
0658 }
0659 
0660 Bool_t CalibTree::Notify() {
0661   // The Notify() function is called when a new file is opened. This

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

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

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

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

0666 
0667   return kTRUE;
0668 }
0669 
0670 void CalibTree::Show(Long64_t entry) {
0671   // Print contents of entry.

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

0673   if (!fChain)
0674     return;
0675   fChain->Show(entry);
0676 }
0677 
0678 Int_t CalibTree::Cut(Long64_t) {
0679   // This function may be called from Loop.

0680   // returns  1 if entry is accepted.

0681   // returns -1 otherwise.

0682   return 1;
0683 }
0684 
0685 Double_t CalibTree::Loop(int loop,
0686                          TFile *fout,
0687                          bool useweight,
0688                          int nMin,
0689                          bool inverse,
0690                          double rmin,
0691                          double rmax,
0692                          int ietaMax,
0693                          int ietaTrack,
0694                          int applyL1Cut,
0695                          double l1Cut,
0696                          bool last,
0697                          double fraction,
0698                          bool writeHisto,
0699                          bool debug) {
0700   if (fChain == 0)
0701     return 0;
0702   Long64_t nbytes(0), nb(0), kprint(0);
0703   Long64_t nentryTot = fChain->GetEntriesFast();
0704   Long64_t nentries = (fraction > 0.01 && fraction < 0.99) ? (Long64_t)(fraction * nentryTot) : nentryTot;
0705   if (detIds.size() == 0) {
0706     for (Long64_t jentry = 0; jentry < nentries; jentry++) {
0707       Long64_t ientry = LoadTree(jentry);
0708       if (ientry < 0)
0709         break;
0710       nb = fChain->GetEntry(jentry);
0711       nbytes += nb;
0712       if (jentry % 1000000 == 0)
0713         std::cout << "Entry " << jentry << " Run " << t_Run << " Event " << t_Event << std::endl;
0714       // Find DetIds contributing to the track

0715       bool selRun = (includeRun_ ? ((t_Run >= runlo_) && (t_Run <= runhi_)) : ((t_Run < runlo_) || (t_Run > runhi_)));
0716       bool selTrack = ((ietaTrack <= 0) || (abs(t_ieta) <= ietaTrack));
0717       if (selRun && (t_nVtx >= nvxlo_) && (t_nVtx <= nvxhi_) && selTrack) {
0718         bool isItRBX(false);
0719         if (cSelect_ != nullptr) {
0720           bool temp = cSelect_->isItRBX(t_DetIds);
0721           if (exclude_)
0722             isItRBX = temp;
0723           else
0724             isItRBX = !(temp);
0725         }
0726         ++kprint;
0727         if (!(isItRBX)) {
0728           for (unsigned int idet = 0; idet < (*t_DetIds).size(); idet++) {
0729             if (selectPhi((*t_DetIds)[idet])) {
0730               unsigned int detid = truncateId((*t_DetIds)[idet], truncateFlag_, debug);
0731               if (debug && (kprint <= 10)) {
0732                 std::cout << "DetId[" << idet << "] Original " << std::hex << (*t_DetIds)[idet] << " truncated "
0733                           << detid << std::dec;
0734               }
0735               if (std::find(detIds.begin(), detIds.end(), detid) == detIds.end()) {
0736                 detIds.push_back(detid);
0737                 if (debug && (kprint <= 10))
0738                   std::cout << " new";
0739               }
0740               if (debug && (kprint <= 10))
0741                 std::cout << std::endl;
0742             }
0743           }
0744           // Also look at the neighbouring cells if available

0745           if (t_DetIds3 != 0) {
0746             for (unsigned int idet = 0; idet < (*t_DetIds3).size(); idet++) {
0747               if (selectPhi((*t_DetIds3)[idet])) {
0748                 unsigned int detid = truncateId((*t_DetIds3)[idet], truncateFlag_, debug);
0749                 if (std::find(detIds.begin(), detIds.end(), detid) == detIds.end()) {
0750                   detIds.push_back(detid);
0751                 }
0752               }
0753             }
0754           }
0755         }
0756       }
0757     }
0758   }
0759   if (debug) {
0760     std::cout << "Total of " << detIds.size() << " detIds and " << histos.size() << " histos found" << std::endl;
0761     // The masks are defined in DataFormats/HcalDetId/interface/HcalDetId.h

0762     for (unsigned int k = 0; k < detIds.size(); ++k) {
0763       int subdet, depth, zside, ieta, iphi;
0764       unpackDetId(detIds[k], subdet, zside, ieta, iphi, depth);
0765       std::cout << "DetId[" << k << "] " << subdet << ":" << zside * ieta << ":" << depth << ":" << iphi << "  "
0766                 << std::hex << detIds[k] << std::dec << std::endl;
0767     }
0768   }
0769   unsigned int k(0);
0770   for (std::map<unsigned int, TH1D *>::const_iterator itr = histos.begin(); itr != histos.end(); ++itr, ++k) {
0771     if (debug) {
0772       std::cout << "histos[" << k << "] " << std::hex << itr->first << std::dec << " " << itr->second;
0773       if (itr->second != 0)
0774         std::cout << " " << itr->second->GetTitle();
0775       std::cout << std::endl;
0776     }
0777     if (itr->second != 0)
0778       itr->second->Delete();
0779   }
0780 
0781   for (unsigned int k = 0; k < detIds.size(); ++k) {
0782     char name[20], title[100];
0783     sprintf(name, "Hist%d_%d", detIds[k], loop);
0784     int subdet, depth, zside, ieta, iphi;
0785     unpackDetId(detIds[k], subdet, zside, ieta, iphi, depth);
0786     sprintf(title, "Correction for Subdet %d #eta %d depth %d (Loop %d)", subdet, zside * ieta, depth, loop);
0787     TH1D *hist = new TH1D(name, title, 100, 0.0, 5.0);
0788     hist->Sumw2();
0789     if (debug)
0790       std::cout << "Book Histo " << k << " " << title << std::endl;
0791     histos[detIds[k]] = hist;
0792   }
0793   std::cout << "Total of " << detIds.size() << " detIds and " << histos.size() << " found in " << nentries << std::endl;
0794 
0795   nbytes = nb = 0;
0796   std::map<unsigned int, myEntry> SumW;
0797   std::map<unsigned int, double> nTrks;
0798 
0799   int ntkgood(0);
0800   for (Long64_t jentry = 0; jentry < nentries; jentry++) {
0801     Long64_t ientry = LoadTree(jentry);
0802     if (ientry < 0)
0803       break;
0804     nb = fChain->GetEntry(jentry);
0805     nbytes += nb;
0806     if (jentry % 1000000 == 0)
0807       std::cout << "Entry " << jentry << " Run " << t_Run << " Event " << t_Event << std::endl;
0808     if (std::find(entries.begin(), entries.end(), jentry) != entries.end())
0809       continue;
0810     bool selRun = (includeRun_ ? ((t_Run >= runlo_) && (t_Run <= runhi_)) : ((t_Run < runlo_) || (t_Run > runhi_)));
0811     if (!selRun)
0812       continue;
0813     if ((t_nVtx < nvxlo_) || (t_nVtx > nvxhi_))
0814       continue;
0815     if (cSelect_ != nullptr) {
0816       if (exclude_) {
0817         if (cSelect_->isItRBX(t_DetIds))
0818           continue;
0819       } else {
0820         if (!(cSelect_->isItRBX(t_ieta, t_iphi)))
0821           continue;
0822       }
0823     }
0824     bool selTrack = ((ietaTrack <= 0) || (abs(t_ieta) <= ietaTrack));
0825     if (!selTrack)
0826       continue;
0827     if ((rcorForm_ == 3) && (cFactor_ != nullptr) && (cFactor_->absent(ientry)))
0828       continue;
0829 
0830     if (debug) {
0831       std::cout << "***Entry (Track) Number : " << ientry << std::endl;
0832       std::cout << "p/eHCal/eMipDR/nDets : " << t_p << "/" << t_eHcal << "/" << t_eMipDR << "/" << (*t_DetIds).size()
0833                 << std::endl;
0834     }
0835     double pmom = (useGen_ && (t_gentrackP > 0)) ? t_gentrackP : t_p;
0836     if (goodTrack()) {
0837       ++ntkgood;
0838       CalibTree::energyCalor en = energyHcal(pmom, jentry, true);
0839       double evWt = (useweight) ? t_EventWeight : 1.0;
0840       if (en.ehcal > 0.001) {
0841         double pufac = (en.Etot > 0) ? (en.ehcal / en.Etot) : 1.0;
0842         double ratio = en.ehcal / (pmom - t_eMipDR);
0843         if (debug)
0844           std::cout << " Weights " << evWt << ":" << pufac << " Energy " << en.Etot2 << ":" << en.Etot << ":" << pmom
0845                     << ":" << t_eMipDR << ":" << t_eHcal << ":" << en.ehcal << " ratio " << ratio << std::endl;
0846         if (loop == 0) {
0847           h_pbyE->Fill(ratio, evWt);
0848           h_Ebyp_bfr->Fill(t_ieta, ratio, evWt);
0849         }
0850         if (last) {
0851           h_Ebyp_aftr->Fill(t_ieta, ratio, evWt);
0852         }
0853         bool l1c(true);
0854         if (applyL1Cut != 0)
0855           l1c = ((t_mindR1 >= l1Cut) || ((applyL1Cut == 1) && (t_DataType == 1)));
0856         if ((rmin >= 0 && ratio > rmin) && (rmax >= 0 && ratio < rmax) && l1c) {
0857           for (unsigned int idet = 0; idet < (*t_DetIds).size(); idet++) {
0858             if (selectPhi((*t_DetIds)[idet])) {
0859               unsigned int id = (*t_DetIds)[idet];
0860               unsigned int detid = truncateId(id, truncateFlag_, false);
0861               double hitEn = 0.0;
0862               if (debug) {
0863                 std::cout << "idet " << idet << " detid/hitenergy : " << std::hex << (*t_DetIds)[idet] << ":" << detid
0864                           << "/" << (*t_HitEnergies)[idet] << std::endl;
0865               }
0866               if (Cprev.find(detid) != Cprev.end())
0867                 hitEn = Cprev[detid].first * (*t_HitEnergies)[idet];
0868               else
0869                 hitEn = (*t_HitEnergies)[idet];
0870               if (cFactor_)
0871                 hitEn *= cFactor_->getCorr(t_Run, id);
0872               double Wi = evWt * hitEn / en.Etot;
0873               double Fac = (inverse) ? (en.ehcal / (pmom - t_eMipDR)) : ((pmom - t_eMipDR) / en.ehcal);
0874               double Fac2 = Wi * Fac * Fac;
0875               TH1D *hist(0);
0876               std::map<unsigned int, TH1D *>::const_iterator itr = histos.find(detid);
0877               if (itr != histos.end())
0878                 hist = itr->second;
0879               if (debug) {
0880                 std::cout << "Det Id " << std::hex << detid << std::dec << " " << hist << std::endl;
0881               }
0882               if (hist != 0)
0883                 hist->Fill(Fac, Wi);  //////histola

0884               Fac *= Wi;
0885               if (SumW.find(detid) != SumW.end()) {
0886                 Wi += SumW[detid].fact0;
0887                 Fac += SumW[detid].fact1;
0888                 Fac2 += SumW[detid].fact2;
0889                 int kount = SumW[detid].kount + 1;
0890                 SumW[detid] = myEntry(kount, Wi, Fac, Fac2);
0891                 nTrks[detid] += evWt;
0892               } else {
0893                 SumW.insert(std::pair<unsigned int, myEntry>(detid, myEntry(1, Wi, Fac, Fac2)));
0894                 nTrks.insert(std::pair<unsigned int, unsigned int>(detid, evWt));
0895               }
0896             }
0897           }
0898         }
0899       }
0900     }
0901   }
0902   if (debug)
0903     std::cout << "# of Good Tracks " << ntkgood << " out of " << nentries << std::endl;
0904   if (loop == 0) {
0905     h_pbyE->Write("h_pbyE");
0906     h_Ebyp_bfr->Write("h_Ebyp_bfr");
0907   }
0908   if (last) {
0909     h_Ebyp_aftr->Write("h_Ebyp_aftr");
0910   }
0911 
0912   std::map<unsigned int, std::pair<double, double> > cfactors;
0913   unsigned int kount(0), kountus(0);
0914   double sumfactor(0);
0915   for (std::map<unsigned int, TH1D *>::const_iterator itr = histos.begin(); itr != histos.end(); ++itr) {
0916     if (writeHisto) {
0917       std::pair<double, double> result_write = fitMean(itr->second, 0);
0918       (itr->second)->Write();
0919     }
0920     // The masks are defined in DataFormats/HcalDetId/interface/HcalDetId.h

0921     int subdet, depth, zside, ieta, iphi;
0922     unpackDetId(itr->first, subdet, zside, ieta, iphi, depth);
0923     if (debug) {
0924       std::cout << "DETID :" << subdet << "  IETA :" << ieta << " HIST ENTRIES :" << (itr->second)->GetEntries()
0925                 << std::endl;
0926     }
0927   }
0928 
0929   if (debug)
0930     std::cout << "Histos with " << histos.size() << " entries\n";
0931   for (std::map<unsigned int, TH1D *>::const_iterator itr = histos.begin(); itr != histos.end(); ++itr, ++kount) {
0932     std::pair<double, double> result = fitMean(itr->second, 0);
0933     double factor = (inverse) ? (2. - result.first) : result.first;
0934     if (debug) {
0935       int subdet, depth, zside, ieta, iphi;
0936       unpackDetId(itr->first, subdet, zside, ieta, iphi, depth);
0937       std::cout << "DetId[" << kount << "] " << subdet << ":" << zside * ieta << ":" << depth << " Factor " << factor
0938                 << " +- " << result.second << std::endl;
0939     }
0940     if (!useMean_) {
0941       cfactors[itr->first] = std::pair<double, double>(factor, result.second);
0942       if (itr->second->GetEntries() > nMin) {
0943         kountus++;
0944         if (factor > 1)
0945           sumfactor += (1 - 1 / factor);
0946         else
0947           sumfactor += (1 - factor);
0948       }
0949     }
0950   }
0951 
0952   if (debug)
0953     std::cout << "SumW with " << SumW.size() << " entries\n";
0954   std::map<unsigned int, myEntry>::const_iterator SumWItr = SumW.begin();
0955   for (; SumWItr != SumW.end(); SumWItr++) {
0956     unsigned int detid = SumWItr->first;
0957     int subdet, depth, zside, ieta, iphi;
0958     unpackDetId(detid, subdet, zside, ieta, iphi, depth);
0959     if (debug) {
0960       std::cout << "Detid|kount|SumWi|SumFac|myId : " << subdet << ":" << zside * ieta << ":" << depth << " | "
0961                 << (SumWItr->second).kount << " | " << (SumWItr->second).fact0 << "|" << (SumWItr->second).fact1 << "|"
0962                 << (SumWItr->second).fact2 << std::endl;
0963     }
0964     double factor = (SumWItr->second).fact1 / (SumWItr->second).fact0;
0965     double dfac1 = ((SumWItr->second).fact2 / (SumWItr->second).fact0 - factor * factor);
0966     if (dfac1 < 0)
0967       dfac1 = 0;
0968     double dfac = sqrt(dfac1 / (SumWItr->second).kount);
0969     if (debug) {
0970       std::cout << "Factor " << factor << " " << dfac1 << " " << dfac << std::endl;
0971     }
0972     if (inverse)
0973       factor = 2. - factor;
0974     if (useMean_) {
0975       cfactors[detid] = std::pair<double, double>(factor, dfac);
0976       if ((SumWItr->second).kount > nMin) {
0977         kountus++;
0978         if (factor > 1)
0979           sumfactor += (1 - 1 / factor);
0980         else
0981           sumfactor += (1 - factor);
0982       }
0983     }
0984   }
0985 
0986   static const int maxch = 500;
0987   double dets[maxch], cfacs[maxch], wfacs[maxch], myId[maxch], nTrk[maxch];
0988   std::cout << "cafctors: " << cfactors.size() << ":" << maxch << std::endl;
0989   kount = 0;
0990   std::map<unsigned int, std::pair<double, double> >::const_iterator itr = cfactors.begin();
0991   for (; itr != cfactors.end(); ++itr, ++kount) {
0992     unsigned int detid = itr->first;
0993     int subdet, depth, zside, ieta, iphi;
0994     unpackDetId(detid, subdet, zside, ieta, iphi, depth);
0995     double id = ieta * zside + 0.25 * (depth - 1);
0996     double factor = (itr->second).first;
0997     double dfac = (itr->second).second;
0998     if (ieta > ietaMax) {
0999       factor = 1;
1000       dfac = 0;
1001     }
1002     std::pair<double, double> cfac(factor, dfac);
1003     if (Cprev.find(detid) != Cprev.end()) {
1004       dfac /= factor;
1005       factor *= Cprev[detid].first;
1006       dfac *= factor;
1007       Cprev[detid] = std::pair<double, double>(factor, dfac);
1008       cfacs[kount] = factor;
1009     } else {
1010       Cprev[detid] = std::pair<double, double>(factor, dfac);
1011       cfacs[kount] = factor;
1012     }
1013     wfacs[kount] = factor;
1014     dets[kount] = detid;
1015     myId[kount] = id;
1016     nTrk[kount] = nTrks[detid];
1017   }
1018   if (higheta_ > 0)
1019     highEtaFactors(ietaMax, debug);
1020 
1021   std::cout << kountus << " detids out of " << kount << " have tracks > " << nMin << std::endl;
1022 
1023   char fname[50];
1024   fout->cd();
1025   TGraph *g_fac1 = new TGraph(kount, dets, cfacs);
1026   sprintf(fname, "Cfacs%d", loop);
1027   g_fac1->SetMarkerStyle(7);
1028   g_fac1->SetMarkerSize(5.0);
1029   g_fac1->Draw("AP");
1030   g_fac1->Write(fname);
1031   TGraph *g_fac2 = new TGraph(kount, dets, wfacs);
1032   sprintf(fname, "Wfacs%d", loop);
1033   g_fac2->SetMarkerStyle(7);
1034   g_fac2->SetMarkerSize(5.0);
1035   g_fac2->Draw("AP");
1036   g_fac2->Write(fname);
1037   TGraph *g_fac3 = new TGraph(kount, myId, cfacs);
1038   sprintf(fname, "CfacsVsMyId%d", loop);
1039   g_fac3->SetMarkerStyle(7);
1040   g_fac3->SetMarkerSize(5.0);
1041   g_fac3->Draw("AP");
1042   g_fac3->Write(fname);
1043   TGraph *g_fac4 = new TGraph(kount, myId, wfacs);
1044   sprintf(fname, "WfacsVsMyId%d", loop);
1045   g_fac4->SetMarkerStyle(7);
1046   g_fac4->SetMarkerSize(5.0);
1047   g_fac4->Draw("AP");
1048   g_fac4->Write(fname);
1049   TGraph *g_nTrk = new TGraph(kount, myId, nTrk);
1050   sprintf(fname, "nTrk");
1051   if (loop == 0) {
1052     g_nTrk->SetMarkerStyle(7);
1053     g_nTrk->SetMarkerSize(5.0);
1054     g_nTrk->Draw("AP");
1055     g_nTrk->Write(fname);
1056   }
1057   std::cout << "The new factors are :" << std::endl;
1058   std::map<unsigned int, std::pair<double, double> >::const_iterator CprevItr = Cprev.begin();
1059   unsigned int indx(0);
1060   for (; CprevItr != Cprev.end(); CprevItr++, indx++) {
1061     unsigned int detid = CprevItr->first;
1062     int subdet, depth, zside, ieta, iphi;
1063     unpackDetId(detid, subdet, zside, ieta, iphi, depth);
1064     std::cout << "DetId[" << indx << "] " << std::hex << detid << std::dec << "(" << ieta * zside << "," << depth
1065               << ") (nTrks:" << nTrks[detid] << ") : " << CprevItr->second.first << " +- " << CprevItr->second.second
1066               << std::endl;
1067   }
1068   double mean = (kountus > 0) ? (sumfactor / kountus) : 0;
1069   std::cout << "Mean deviation " << mean << " from 1 for " << kountus << " DetIds" << std::endl;
1070   h_cvg->SetBinContent(loop + 1, mean);
1071   if (last)
1072     h_cvg->Write("Cvg0");
1073   return mean;
1074 }
1075 
1076 bool CalibTree::goodTrack() {
1077   bool ok(true);
1078   double cut(2.0);
1079   double pmom = (useGen_ && (t_gentrackP > 0)) ? t_gentrackP : t_p;
1080   if (sysmode_ == 1) {
1081     ok =
1082         ((t_qltyFlag) && (t_hmaxNearP < cut) && (t_eMipDR < 1.0) && (t_mindR1 > 1.0) && (pmom > 40.0) && (pmom < 60.0));
1083   } else if (sysmode_ == 2) {
1084     ok = ((t_qltyFlag) && (t_qltyPVFlag) && (t_hmaxNearP < cut) && (t_eMipDR < 1.0) && (t_mindR1 > 1.0) &&
1085           (pmom > 40.0) && (pmom < 60.0));
1086   } else if (sysmode_ == 3) {
1087     ok =
1088         ((t_selectTk) && (t_hmaxNearP < cut) && (t_eMipDR < 1.0) && (t_mindR1 > 1.0) && (pmom > 40.0) && (pmom < 60.0));
1089   } else if (sysmode_ == 4) {
1090     ok = ((t_selectTk) && (t_qltyMissFlag) && (t_hmaxNearP < 0.0) && (t_eMipDR < 1.0) && (t_mindR1 > 1.0) &&
1091           (pmom > 40.0) && (pmom < 60.0));
1092   } else if (sysmode_ == 5) {
1093     ok = ((t_selectTk) && (t_qltyMissFlag) && (t_hmaxNearP < cut) && (t_eMipDR < 0.5) && (t_mindR1 > 1.0) &&
1094           (pmom > 40.0) && (pmom < 60.0));
1095   } else if (sysmode_ == 6) {
1096     ok = ((t_selectTk) && (t_qltyMissFlag) && (t_hmaxNearP < cut) && (t_eMipDR < 2.0) && (t_mindR1 > 1.0) &&
1097           (pmom > 40.0) && (pmom < 60.0));
1098   } else if (sysmode_ == 7) {
1099     ok = ((t_selectTk) && (t_qltyMissFlag) && (t_hmaxNearP < cut) && (t_eMipDR < 1.0) && (t_mindR1 > 0.5) &&
1100           (pmom > 40.0) && (pmom < 60.0));
1101   } else {
1102     if (sysmode_ < 0) {
1103       double eta = (t_ieta > 0) ? t_ieta : -t_ieta;
1104       if (sysmode_ == -2)
1105         cut = 8.0 * exp(eta * log2by18_);
1106       else
1107         cut = 10.0;
1108     }
1109     ok = ((t_selectTk) && (t_qltyMissFlag) && (t_hmaxNearP < cut) && (t_eMipDR < 1.0) && (t_mindR1 > 1.0) &&
1110           (pmom > 40.0) && (pmom < 60.0));
1111   }
1112   return ok;
1113 }
1114 
1115 void CalibTree::writeCorrFactor(const char *corrFileName, int ietaMax) {
1116   ofstream myfile;
1117   myfile.open(corrFileName);
1118   if (!myfile.is_open()) {
1119     std::cout << "** ERROR: Can't open '" << corrFileName << std::endl;
1120   } else {
1121     myfile << "#" << std::setprecision(4) << std::setw(10) << "detId" << std::setw(10) << "ieta" << std::setw(10)
1122            << "depth" << std::setw(15) << "corrFactor" << std::endl;
1123     std::map<unsigned int, std::pair<double, double> >::const_iterator itr;
1124     for (itr = Cprev.begin(); itr != Cprev.end(); ++itr) {
1125       unsigned int detId = itr->first;
1126       int subdet, depth, zside, ieta, iphi;
1127       unpackDetId(detId, subdet, zside, ieta, iphi, depth);
1128       if (ieta <= ietaMax) {
1129         myfile << std::setw(10) << std::hex << detId << std::setw(10) << std::dec << zside * ieta << std::setw(10)
1130                << depth << std::setw(10) << itr->second.first << " " << std::setw(10) << itr->second.second
1131                << std::endl;
1132         std::cout << itr->second.first << ",";
1133       }
1134     }
1135     myfile.close();
1136     std::cout << std::endl;
1137   }
1138 }
1139 
1140 bool CalibTree::selectPhi(unsigned int detId) {
1141   bool flag(true);
1142   // The masks are defined in DataFormats/HcalDetId/interface/HcalDetId.h

1143   if ((phimin_ > 1) || (phimax_ < 72) || (zside_ != 0)) {
1144     int subdet, depth, zside, ieta, iphi;
1145     unpackDetId(detId, subdet, zside, ieta, iphi, depth);
1146     if (phimin_ > 1 || phimax_ < 72) {
1147       if ((iphi < phimin_) || (iphi > phimax_))
1148         flag = false;
1149     }
1150     if (zside_ != 0) {
1151       if (zside != zside_)
1152         flag = false;
1153     }
1154   }
1155   return flag;
1156 }
1157 
1158 std::pair<double, double> CalibTree::fitMean(TH1D *hist, int mode) {
1159   std::pair<double, double> results = std::pair<double, double>(1, 0);
1160   if (hist != 0) {
1161     double mean = hist->GetMean(), rms = hist->GetRMS();
1162     double LowEdge(0.7), HighEdge(1.3);
1163     char option[20];
1164     if (mode == 1) {
1165       LowEdge = mean - 1.5 * rms;
1166       HighEdge = mean + 1.5 * rms;
1167       int nbin = hist->GetNbinsX();
1168       if (LowEdge < hist->GetBinLowEdge(1))
1169         LowEdge = hist->GetBinLowEdge(1);
1170       if (HighEdge > hist->GetBinLowEdge(nbin) + hist->GetBinWidth(nbin))
1171         HighEdge = hist->GetBinLowEdge(nbin) + hist->GetBinWidth(nbin);
1172     }
1173     if (hist->GetEntries() > 100)
1174       sprintf(option, "+QRLS");
1175     else
1176       sprintf(option, "+QRWLS");
1177     double value(mean);
1178     double error = rms / sqrt(hist->GetEntries());
1179     if (hist->GetEntries() > 20) {
1180       TFitResultPtr Fit = hist->Fit("gaus", option, "", LowEdge, HighEdge);
1181       value = Fit->Value(1);
1182       error = Fit->FitResult::Error(1);
1183       /*

1184       LowEdge  = value - 1.5*error;

1185       HighEdge = value + 1.5*error;

1186       value    = Fit->Value(1);

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

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

1189       */
1190     }
1191     results = std::pair<double, double>(value, error);
1192   }
1193   return results;
1194 }
1195 
1196 void CalibTree::makeplots(double rmin, double rmax, int ietaMax, bool useweight, double fraction, bool debug) {
1197   if (fChain == 0)
1198     return;
1199   Long64_t nentryTot = fChain->GetEntriesFast();
1200   Long64_t nentries = (fraction > 0.01 && fraction < 0.99) ? (Long64_t)(fraction * nentryTot) : nentryTot;
1201 
1202   // Book the histograms

1203   std::map<int, std::pair<TH1D *, TH1D *> > histos;
1204   for (int ieta = -ietaMax; ieta <= ietaMax; ++ieta) {
1205     char name[20], title[100];
1206     sprintf(name, "begin%d", ieta);
1207     if (ieta == 0)
1208       sprintf(title, "Ratio at start");
1209     else
1210       sprintf(title, "Ratio at start for i#eta=%d", ieta);
1211     TH1D *h1 = new TH1D(name, title, 50, rmin, rmax);
1212     h1->Sumw2();
1213     sprintf(name, "end%d", ieta);
1214     if (ieta == 0)
1215       sprintf(title, "Ratio at the end");
1216     else
1217       sprintf(title, "Ratio at the end for i#eta=%d", ieta);
1218     TH1D *h2 = new TH1D(name, title, 50, rmin, rmax);
1219     h2->Sumw2();
1220     histos[ieta] = std::pair<TH1D *, TH1D *>(h1, h2);
1221   }
1222   //Fill the histograms

1223   Long64_t nbytes(0), nb(0);
1224   for (Long64_t jentry = 0; jentry < nentries; jentry++) {
1225     Long64_t ientry = LoadTree(jentry);
1226     if (ientry < 0)
1227       break;
1228     if (std::find(entries.begin(), entries.end(), jentry) != entries.end())
1229       break;
1230     nb = fChain->GetEntry(jentry);
1231     nbytes += nb;
1232     if (goodTrack()) {
1233       double pmom = (useGen_ && (t_gentrackP > 0)) ? t_gentrackP : t_p;
1234       CalibTree::energyCalor en1 = energyHcal(pmom, jentry, false);
1235       CalibTree::energyCalor en2 = energyHcal(pmom, jentry, true);
1236       if ((en1.ehcal > 0.001) && (en2.ehcal > 0.001)) {
1237         double evWt = (useweight) ? t_EventWeight : 1.0;
1238         double ratioi = en1.ehcal / (pmom - t_eMipDR);
1239         double ratiof = en2.ehcal / (pmom - t_eMipDR);
1240         if (t_ieta >= -ietaMax && t_ieta <= ietaMax && t_ieta != 0) {
1241           if (ratioi >= rmin && ratioi <= rmax) {
1242             histos[0].first->Fill(ratioi, evWt);
1243             histos[t_ieta].first->Fill(ratioi, evWt);
1244           }
1245           if (ratiof >= rmin && ratiof <= rmax) {
1246             histos[0].second->Fill(ratiof, evWt);
1247             histos[t_ieta].second->Fill(ratiof, evWt);
1248           }
1249         }
1250       }
1251     }
1252   }
1253 
1254   //Fit the histograms

1255   TH1D *hbef1 = new TH1D("Eta1Bf", "Mean vs i#eta", 2 * ietaMax, -ietaMax, ietaMax);
1256   TH1D *hbef2 = new TH1D("Eta2Bf", "Median vs i#eta", 2 * ietaMax, -ietaMax, ietaMax);
1257   TH1D *haft1 = new TH1D("Eta1Af", "Mean vs i#eta", 2 * ietaMax, -ietaMax, ietaMax);
1258   TH1D *haft2 = new TH1D("Eta2Af", "Median vs i#eta", 2 * ietaMax, -ietaMax, ietaMax);
1259   for (int ieta = -ietaMax; ieta <= ietaMax; ++ieta) {
1260     int bin = (ieta < 0) ? (ieta + ietaMax + 1) : (ieta + ietaMax);
1261     TH1D *h1 = histos[ieta].first;
1262     double mean1 = h1->GetMean();
1263     double err1 = h1->GetMeanError();
1264     std::pair<double, double> fit1 = fitMean(h1, 1);
1265     if (debug) {
1266       std::cout << ieta << " " << h1->GetName() << " " << mean1 << " +- " << err1 << " and " << fit1.first << " +- "
1267                 << fit1.second << std::endl;
1268     }
1269     if (ieta != 0) {
1270       hbef1->SetBinContent(bin, mean1);
1271       hbef1->SetBinError(bin, err1);
1272       hbef2->SetBinContent(bin, fit1.first);
1273       hbef2->SetBinError(bin, fit1.second);
1274     }
1275     h1->Write();
1276     TH1D *h2 = histos[ieta].second;
1277     double mean2 = h2->GetMean();
1278     double err2 = h2->GetMeanError();
1279     std::pair<double, double> fit2 = fitMean(h2, 1);
1280     if (debug) {
1281       std::cout << ieta << " " << h2->GetName() << " " << mean2 << " +- " << err2 << " and " << fit2.first << " +- "
1282                 << fit2.second << std::endl;
1283     }
1284     if (ieta != 0) {
1285       haft1->SetBinContent(bin, mean2);
1286       haft1->SetBinError(bin, err2);
1287       haft2->SetBinContent(bin, fit2.first);
1288       haft2->SetBinError(bin, fit2.second);
1289     }
1290     h2->Write();
1291   }
1292   fitPol0(hbef1, debug);
1293   fitPol0(hbef2, debug);
1294   fitPol0(haft1, debug);
1295   fitPol0(haft2, debug);
1296 }
1297 
1298 void CalibTree::fitPol0(TH1D *hist, bool debug) {
1299   hist->GetXaxis()->SetTitle("i#eta");
1300   hist->GetYaxis()->SetTitle("<E_{HCAL}/(p-E_{ECAL})>");
1301   hist->GetYaxis()->SetRangeUser(0.4, 1.6);
1302   TFitResultPtr Fit = hist->Fit("pol0", "+QRWLS");
1303   if (debug) {
1304     std::cout << "Fit to Pol0 to " << hist->GetTitle() << ": " << Fit->Value(0) << " +- " << Fit->FitResult::Error(0)
1305               << std::endl;
1306   }
1307   hist->Write();
1308 }
1309 
1310 void CalibTree::highEtaFactors(int ietaMax, bool debug) {
1311   std::map<unsigned int, std::pair<double, double> >::const_iterator itr;
1312   std::pair<double, double> cfacp, cfacn;
1313   cfacp = cfacn = std::pair<double, double>(1.0, 0.0);
1314   for (itr = Cprev.begin(); itr != Cprev.end(); ++itr) {
1315     unsigned int detid = itr->first;
1316     int subdet, depth, zside, ieta, iphi;
1317     unpackDetId(detid, subdet, zside, ieta, iphi, depth);
1318     if ((ieta == ietaMax) && (depth == 1)) {
1319       if (zside > 0)
1320         cfacp = itr->second;
1321       else
1322         cfacn = itr->second;
1323     }
1324   }
1325   if (debug) {
1326     std::cout << "Correction factor for (" << ietaMax << ",1) = " << cfacp.first << ":" << cfacp.second << " and ("
1327               << -ietaMax << ",1) = " << cfacn.first << ":" << cfacn.second << std::endl;
1328   }
1329   for (itr = Cprev.begin(); itr != Cprev.end(); ++itr) {
1330     unsigned int detid = itr->first;
1331     int subdet, depth, zside, ieta, iphi;
1332     unpackDetId(detid, subdet, zside, ieta, iphi, depth);
1333     if (ieta > ietaMax) {
1334       Cprev[detid] = (zside > 0) ? cfacp : cfacn;
1335       if (debug) {
1336         std::cout << "Set correction factor for (" << zside * ieta << "," << depth << ") = " << Cprev[detid].first
1337                   << ":" << Cprev[detid].second << std::endl;
1338       }
1339     }
1340   }
1341 }
1342 
1343 CalibTree::energyCalor CalibTree::energyHcal(double pmom, const Long64_t &entry, bool final) {
1344   double etot = t_eHcal;
1345   double etot2 = t_eHcal;
1346   double ediff = (t_eHcal30 - t_eHcal10);
1347   if (final) {
1348     etot = etot2 = 0;
1349     for (unsigned int idet = 0; idet < (*t_DetIds).size(); idet++) {
1350       if (selectPhi((*t_DetIds)[idet])) {
1351         unsigned int id = (*t_DetIds)[idet];
1352         double hitEn(0);
1353         unsigned int detid = truncateId(id, truncateFlag_, false);
1354         if (Cprev.find(detid) != Cprev.end())
1355           hitEn = Cprev[detid].first * (*t_HitEnergies)[idet];
1356         else
1357           hitEn = (*t_HitEnergies)[idet];
1358         if (cFactor_)
1359           hitEn *= cFactor_->getCorr(t_Run, id);
1360         etot += hitEn;
1361         etot2 += ((*t_HitEnergies)[idet]);
1362       }
1363     }
1364     // Now the outer cone

1365     double etot1(0), etot3(0);
1366     if (t_DetIds1 != 0 && t_DetIds3 != 0) {
1367       for (unsigned int idet = 0; idet < (*t_DetIds1).size(); idet++) {
1368         if (selectPhi((*t_DetIds1)[idet])) {
1369           unsigned int id = (*t_DetIds1)[idet];
1370           unsigned int detid = truncateId(id, truncateFlag_, false);
1371           double hitEn(0);
1372           if (Cprev.find(detid) != Cprev.end())
1373             hitEn = Cprev[detid].first * (*t_HitEnergies1)[idet];
1374           else
1375             hitEn = (*t_HitEnergies1)[idet];
1376           if (cFactor_)
1377             hitEn *= cFactor_->getCorr(t_Run, id);
1378           etot1 += hitEn;
1379         }
1380       }
1381       for (unsigned int idet = 0; idet < (*t_DetIds3).size(); idet++) {
1382         if (selectPhi((*t_DetIds3)[idet])) {
1383           unsigned int id = (*t_DetIds3)[idet];
1384           unsigned int detid = truncateId(id, truncateFlag_, false);
1385           double hitEn(0);
1386           if (Cprev.find(detid) != Cprev.end())
1387             hitEn = Cprev[detid].first * (*t_HitEnergies3)[idet];
1388           else
1389             hitEn = (*t_HitEnergies3)[idet];
1390           if (cFactor_)
1391             hitEn *= cFactor_->getCorr(t_Run, id);
1392           etot3 += hitEn;
1393         }
1394       }
1395     }
1396     ediff = etot3 - etot1;
1397   }
1398   // PU correction only for loose isolation cut

1399   double ehcal = (((rcorForm_ == 3) && (cFactor_ != nullptr))
1400                       ? (etot * cFactor_->getCorr(entry))
1401                       : ((puCorr_ == 0) ? etot
1402                                         : ((puCorr_ < 0) ? (etot * puFactor(-puCorr_, t_ieta, pmom, etot, ediff))
1403                                                          : puFactorRho(puCorr_, t_ieta, t_rhoh, etot))));
1404   return CalibTree::energyCalor(etot, etot2, ehcal);
1405 }
[ 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!