Back to home page

Project CMSSW displayed by LXR
 
 

    

File indexing completed on 2022-11-02 03:49:50 

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, rbx, exclude, higheta, fraction, writeDebugHisto,

0009 //      debug, nmax);

0010 //

0011 //  where:

0012 //

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

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

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

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

0017 //                              ("HcalIsoTrkAnalyzer")

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

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

0020 //                              ("Silver_out.root")

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

0022 //                              factors ("Silver_corr.txt")

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

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

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

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

0027 //                              to be used for raddam/depth dependent

0028 //                              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)     = Flag to treat different depths differently (0)

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

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

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

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

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

0059 //                              (Default 0)

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

0077 //                              HEP17 should correspond to 217

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

0079 //                              considered (false)

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

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

0082 //                              and depth 1 (1) (default 1)

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

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

0085 //                              in o/p file (false)

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

0087 //                              (false)

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

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

0090 //

0091 //  doIt(inFileName, dupFileName)

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

0093 ////////////////////////////////////////////////////////////////////////////////

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

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

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

0227   TH1D *h_pbyE, *h_cvg;
0228   TProfile *h_Ebyp_bfr, *h_Ebyp_aftr;
0229 
0230 private:
0231   // Declaration of leaf types

0232   Int_t t_Run;
0233   Int_t t_Event;
0234   Int_t t_DataType;
0235   Int_t t_ieta;
0236   Int_t t_iphi;
0237   Double_t t_EventWeight;
0238   Int_t t_nVtx;
0239   Int_t t_nTrk;
0240   Int_t t_goodPV;
0241   Double_t t_l1pt;
0242   Double_t t_l1eta;
0243   Double_t t_l1phi;
0244   Double_t t_l3pt;
0245   Double_t t_l3eta;
0246   Double_t t_l3phi;
0247   Double_t t_p;
0248   Double_t t_pt;
0249   Double_t t_phi;
0250   Double_t t_mindR1;
0251   Double_t t_mindR2;
0252   Double_t t_eMipDR;
0253   Double_t t_eHcal;
0254   Double_t t_eHcal10;
0255   Double_t t_eHcal30;
0256   Double_t t_hmaxNearP;
0257   Double_t t_rhoh;
0258   Bool_t t_selectTk;
0259   Bool_t t_qltyFlag;
0260   Bool_t t_qltyMissFlag;
0261   Bool_t t_qltyPVFlag;
0262   Double_t t_gentrackP;
0263   std::vector<unsigned int> *t_DetIds;
0264   std::vector<double> *t_HitEnergies;
0265   std::vector<bool> *t_trgbits;
0266   std::vector<unsigned int> *t_DetIds1;
0267   std::vector<unsigned int> *t_DetIds3;
0268   std::vector<double> *t_HitEnergies1;
0269   std::vector<double> *t_HitEnergies3;
0270 
0271   // List of branches

0272   TBranch *b_t_Run;           //!

0273   TBranch *b_t_Event;         //!

0274   TBranch *b_t_DataType;      //!

0275   TBranch *b_t_ieta;          //!

0276   TBranch *b_t_iphi;          //!

0277   TBranch *b_t_EventWeight;   //!

0278   TBranch *b_t_nVtx;          //!

0279   TBranch *b_t_nTrk;          //!

0280   TBranch *b_t_goodPV;        //!

0281   TBranch *b_t_l1pt;          //!

0282   TBranch *b_t_l1eta;         //!

0283   TBranch *b_t_l1phi;         //!

0284   TBranch *b_t_l3pt;          //!

0285   TBranch *b_t_l3eta;         //!

0286   TBranch *b_t_l3phi;         //!

0287   TBranch *b_t_p;             //!

0288   TBranch *b_t_pt;            //!

0289   TBranch *b_t_phi;           //!

0290   TBranch *b_t_mindR1;        //!

0291   TBranch *b_t_mindR2;        //!

0292   TBranch *b_t_eMipDR;        //!

0293   TBranch *b_t_eHcal;         //!

0294   TBranch *b_t_eHcal10;       //!

0295   TBranch *b_t_eHcal30;       //!

0296   TBranch *b_t_hmaxNearP;     //!

0297   TBranch *b_t_rhoh;          //!

0298   TBranch *b_t_selectTk;      //!

0299   TBranch *b_t_qltyFlag;      //!

0300   TBranch *b_t_qltyMissFlag;  //!

0301   TBranch *b_t_qltyPVFlag;    //!

0302   TBranch *b_t_gentrackP;     //!

0303   TBranch *b_t_DetIds;        //!

0304   TBranch *b_t_HitEnergies;   //!

0305   TBranch *b_t_trgbits;       //!

0306   TBranch *b_t_DetIds1;       //!

0307   TBranch *b_t_DetIds3;       //!

0308   TBranch *b_t_HitEnergies1;  //!

0309   TBranch *b_t_HitEnergies3;  //!

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

0582   if (!fChain)
0583     return 0;
0584   return fChain->GetEntry(entry);
0585 }
0586 
0587 Long64_t CalibTree::LoadTree(Long64_t entry) {
0588   // Set the environment to read one entry

0589   if (!fChain)
0590     return -5;
0591   Long64_t centry = fChain->LoadTree(entry);
0592   if (centry < 0)
0593     return centry;
0594   if (fChain->GetTreeNumber() != fCurrent) {
0595     fCurrent = fChain->GetTreeNumber();
0596     Notify();
0597   }
0598   return centry;
0599 }
0600 
0601 void CalibTree::Init(TChain *tree, const char *dupFileName) {
0602   // The Init() function is called when the selector needs to initialize

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

0604   // pointers of the tree will be set.

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

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

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

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

0609 
0610   // Set object pointer

0611   t_DetIds = 0;
0612   t_HitEnergies = 0;
0613   t_trgbits = 0;
0614   t_DetIds1 = 0;
0615   t_DetIds3 = 0;
0616   t_HitEnergies1 = 0;
0617   t_HitEnergies3 = 0;
0618   // Set branch addresses and branch pointers

0619   fChain = tree;
0620   if (!tree)
0621     return;
0622   fCurrent = -1;
0623   fChain->SetMakeClass(1);
0624 
0625   fChain->SetBranchAddress("t_Run", &t_Run, &b_t_Run);
0626   fChain->SetBranchAddress("t_Event", &t_Event, &b_t_Event);
0627   fChain->SetBranchAddress("t_DataType", &t_DataType, &b_t_DataType);
0628   fChain->SetBranchAddress("t_ieta", &t_ieta, &b_t_ieta);
0629   fChain->SetBranchAddress("t_iphi", &t_iphi, &b_t_iphi);
0630   fChain->SetBranchAddress("t_EventWeight", &t_EventWeight, &b_t_EventWeight);
0631   fChain->SetBranchAddress("t_nVtx", &t_nVtx, &b_t_nVtx);
0632   fChain->SetBranchAddress("t_nTrk", &t_nTrk, &b_t_nTrk);
0633   fChain->SetBranchAddress("t_goodPV", &t_goodPV, &b_t_goodPV);
0634   fChain->SetBranchAddress("t_l1pt", &t_l1pt, &b_t_l1pt);
0635   fChain->SetBranchAddress("t_l1eta", &t_l1eta, &b_t_l1eta);
0636   fChain->SetBranchAddress("t_l1phi", &t_l1phi, &b_t_l1phi);
0637   fChain->SetBranchAddress("t_l3pt", &t_l3pt, &b_t_l3pt);
0638   fChain->SetBranchAddress("t_l3eta", &t_l3eta, &b_t_l3eta);
0639   fChain->SetBranchAddress("t_l3phi", &t_l3phi, &b_t_l3phi);
0640   fChain->SetBranchAddress("t_p", &t_p, &b_t_p);
0641   fChain->SetBranchAddress("t_pt", &t_pt, &b_t_pt);
0642   fChain->SetBranchAddress("t_phi", &t_phi, &b_t_phi);
0643   fChain->SetBranchAddress("t_mindR1", &t_mindR1, &b_t_mindR1);
0644   fChain->SetBranchAddress("t_mindR2", &t_mindR2, &b_t_mindR2);
0645   fChain->SetBranchAddress("t_eMipDR", &t_eMipDR, &b_t_eMipDR);
0646   fChain->SetBranchAddress("t_eHcal", &t_eHcal, &b_t_eHcal);
0647   fChain->SetBranchAddress("t_eHcal10", &t_eHcal10, &b_t_eHcal10);
0648   fChain->SetBranchAddress("t_eHcal30", &t_eHcal30, &b_t_eHcal30);
0649   fChain->SetBranchAddress("t_hmaxNearP", &t_hmaxNearP, &b_t_hmaxNearP);
0650   fChain->SetBranchAddress("t_rhoh", &t_rhoh, &b_t_rhoh);
0651   fChain->SetBranchAddress("t_selectTk", &t_selectTk, &b_t_selectTk);
0652   fChain->SetBranchAddress("t_qltyFlag", &t_qltyFlag, &b_t_qltyFlag);
0653   fChain->SetBranchAddress("t_qltyMissFlag", &t_qltyMissFlag, &b_t_qltyMissFlag);
0654   fChain->SetBranchAddress("t_qltyPVFlag", &t_qltyPVFlag, &b_t_qltyPVFlag);
0655   fChain->SetBranchAddress("t_gentrackP", &t_gentrackP, &b_t_gentrackP);
0656   fChain->SetBranchAddress("t_DetIds", &t_DetIds, &b_t_DetIds);
0657   fChain->SetBranchAddress("t_HitEnergies", &t_HitEnergies, &b_t_HitEnergies);
0658   fChain->SetBranchAddress("t_trgbits", &t_trgbits, &b_t_trgbits);
0659   fChain->SetBranchAddress("t_DetIds1", &t_DetIds1, &b_t_DetIds1);
0660   fChain->SetBranchAddress("t_DetIds3", &t_DetIds3, &b_t_DetIds3);
0661   fChain->SetBranchAddress("t_HitEnergies1", &t_HitEnergies1, &b_t_HitEnergies1);
0662   fChain->SetBranchAddress("t_HitEnergies3", &t_HitEnergies3, &b_t_HitEnergies3);
0663   Notify();
0664 
0665   if (std::string(dupFileName) != "") {
0666     std::ifstream infil1(dupFileName);
0667     if (!infil1.is_open()) {
0668       std::cout << "Cannot open " << dupFileName << std::endl;
0669     } else {
0670       while (1) {
0671         Long64_t jentry;
0672         infil1 >> jentry;
0673         if (!infil1.good())
0674           break;
0675         entries.push_back(jentry);
0676       }
0677       infil1.close();
0678       std::cout << "Reads a list of " << entries.size() << " events from " << dupFileName << std::endl;
0679     }
0680   } else {
0681     std::cout << "No duplicate events in the input file" << std::endl;
0682   }
0683 }
0684 
0685 Bool_t CalibTree::Notify() {
0686   // The Notify() function is called when a new file is opened. This

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

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

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

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

0691 
0692   return kTRUE;
0693 }
0694 
0695 void CalibTree::Show(Long64_t entry) {
0696   // Print contents of entry.

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

0698   if (!fChain)
0699     return;
0700   fChain->Show(entry);
0701 }
0702 
0703 Int_t CalibTree::Cut(Long64_t) {
0704   // This function may be called from Loop.

0705   // returns  1 if entry is accepted.

0706   // returns -1 otherwise.

0707   return 1;
0708 }
0709 
0710 Double_t CalibTree::Loop(int loop,
0711                          TFile *fout,
0712                          bool useweight,
0713                          int nMin,
0714                          bool inverse,
0715                          double rmin,
0716                          double rmax,
0717                          int ietaMax,
0718                          int ietaTrack,
0719                          int applyL1Cut,
0720                          double l1Cut,
0721                          bool last,
0722                          double fraction,
0723                          bool writeHisto,
0724                          bool debug,
0725                          Long64_t nmax) {
0726   Long64_t nbytes(0), nb(0);
0727   Long64_t nentryTot = fChain->GetEntriesFast();
0728   Long64_t nentries = (fraction > 0.01 && fraction < 0.99) ? (Long64_t)(fraction * nentryTot) : nentryTot;
0729   if ((nentries > nmax) && (nmax > 0))
0730     nentries = nmax;
0731 
0732   bookHistos(loop, debug);
0733   std::map<unsigned int, myEntry> SumW;
0734   std::map<unsigned int, double> nTrks;
0735 
0736   int ntkgood(0);
0737   for (Long64_t jentry = 0; jentry < nentries; jentry++) {
0738     Long64_t ientry = LoadTree(jentry);
0739     if (ientry < 0)
0740       break;
0741     nb = fChain->GetEntry(jentry);
0742     nbytes += nb;
0743     if (jentry % 1000000 == 0)
0744       std::cout << "Entry " << jentry << " Run " << t_Run << " Event " << t_Event << std::endl;
0745     if (std::find(entries.begin(), entries.end(), jentry) != entries.end())
0746       continue;
0747     bool selRun = (includeRun_ ? ((t_Run >= runlo_) && (t_Run <= runhi_)) : ((t_Run < runlo_) || (t_Run > runhi_)));
0748     if (!selRun)
0749       continue;
0750     if ((t_nVtx < nvxlo_) || (t_nVtx > nvxhi_))
0751       continue;
0752     if (cSelect_ != nullptr) {
0753       if (exclude_) {
0754         if (cSelect_->isItRBX(t_DetIds))
0755           continue;
0756       } else {
0757         if (!(cSelect_->isItRBX(t_ieta, t_iphi)))
0758           continue;
0759       }
0760     }
0761     bool selTrack = ((ietaTrack <= 0) || (abs(t_ieta) <= ietaTrack));
0762     if (!selTrack)
0763       continue;
0764     if ((rcorForm_ == 3) && (cFactor_ != nullptr) && (cFactor_->absent(ientry)))
0765       continue;
0766 
0767     if (debug) {
0768       std::cout << "***Entry (Track) Number : " << ientry << std::endl;
0769       std::cout << "p/eHCal/eMipDR/nDets : " << t_p << "/" << t_eHcal << "/" << t_eMipDR << "/" << (*t_DetIds).size()
0770                 << std::endl;
0771     }
0772     double pmom = (useGen_ && (t_gentrackP > 0)) ? t_gentrackP : t_p;
0773     if (goodTrack()) {
0774       ++ntkgood;
0775       CalibTree::energyCalor en = energyHcal(pmom, jentry, true);
0776       double evWt = (useweight) ? t_EventWeight : 1.0;
0777       if (en.ehcal > 0.001) {
0778         double pufac = (en.Etot > 0) ? (en.ehcal / en.Etot) : 1.0;
0779         double ratio = en.ehcal / (pmom - t_eMipDR);
0780         if (debug)
0781           std::cout << " Weights " << evWt << ":" << pufac << " Energy " << en.Etot2 << ":" << en.Etot << ":" << pmom
0782                     << ":" << t_eMipDR << ":" << t_eHcal << ":" << en.ehcal << " ratio " << ratio << std::endl;
0783         if (loop == 0) {
0784           h_pbyE->Fill(ratio, evWt);
0785           h_Ebyp_bfr->Fill(t_ieta, ratio, evWt);
0786         }
0787         if (last) {
0788           h_Ebyp_aftr->Fill(t_ieta, ratio, evWt);
0789         }
0790         bool l1c(true);
0791         if (applyL1Cut != 0)
0792           l1c = ((t_mindR1 >= l1Cut) || ((applyL1Cut == 1) && (t_DataType == 1)));
0793         if ((rmin >= 0 && ratio > rmin) && (rmax >= 0 && ratio < rmax) && l1c) {
0794           for (unsigned int idet = 0; idet < (*t_DetIds).size(); idet++) {
0795             if (selectPhi((*t_DetIds)[idet])) {
0796               unsigned int id = (*t_DetIds)[idet];
0797               unsigned int detid = truncateId(id, truncateFlag_, false);
0798               double hitEn = 0.0;
0799               if (debug) {
0800                 std::cout << "idet " << idet << " detid/hitenergy : " << std::hex << (*t_DetIds)[idet] << ":" << detid
0801                           << "/" << (*t_HitEnergies)[idet] << std::endl;
0802               }
0803               if (Cprev.find(detid) != Cprev.end())
0804                 hitEn = Cprev[detid].first * (*t_HitEnergies)[idet];
0805               else
0806                 hitEn = (*t_HitEnergies)[idet];
0807               if (cFactor_)
0808                 hitEn *= cFactor_->getCorr(t_Run, id);
0809               double Wi = evWt * hitEn / en.Etot;
0810               double Fac = (inverse) ? (en.ehcal / (pmom - t_eMipDR)) : ((pmom - t_eMipDR) / en.ehcal);
0811               double Fac2 = Wi * Fac * Fac;
0812               TH1D *hist(0);
0813               std::map<unsigned int, TH1D *>::const_iterator itr = histos_.find(detid);
0814               if (itr != histos_.end())
0815                 hist = itr->second;
0816               if (debug) {
0817                 std::cout << "Det Id " << std::hex << detid << std::dec << " " << hist << std::endl;
0818               }
0819               if (hist != 0)
0820                 hist->Fill(Fac, Wi);  //////histola

0821               Fac *= Wi;
0822               if (SumW.find(detid) != SumW.end()) {
0823                 Wi += SumW[detid].fact0;
0824                 Fac += SumW[detid].fact1;
0825                 Fac2 += SumW[detid].fact2;
0826                 int kount = SumW[detid].kount + 1;
0827                 SumW[detid] = myEntry(kount, Wi, Fac, Fac2);
0828                 nTrks[detid] += evWt;
0829               } else {
0830                 SumW.insert(std::pair<unsigned int, myEntry>(detid, myEntry(1, Wi, Fac, Fac2)));
0831                 nTrks.insert(std::pair<unsigned int, unsigned int>(detid, evWt));
0832               }
0833             }
0834           }
0835         }
0836       }
0837     }
0838   }
0839   if (debug)
0840     std::cout << "# of Good Tracks " << ntkgood << " out of " << nentries << std::endl;
0841   if (loop == 0) {
0842     h_pbyE->Write("h_pbyE");
0843     h_Ebyp_bfr->Write("h_Ebyp_bfr");
0844   }
0845   if (last) {
0846     h_Ebyp_aftr->Write("h_Ebyp_aftr");
0847   }
0848 
0849   std::map<unsigned int, std::pair<double, double> > cfactors;
0850   unsigned int kount(0), kountus(0);
0851   double sumfactor(0);
0852   for (std::map<unsigned int, TH1D *>::const_iterator itr = histos_.begin(); itr != histos_.end(); ++itr) {
0853     if (writeHisto) {
0854       //      std::pair<double, double> result_write = fitMean(itr->second, 0);

0855       (itr->second)->Write();
0856     }
0857     // The masks are defined in DataFormats/HcalDetId/interface/HcalDetId.h

0858     int subdet, depth, zside, ieta, iphi;
0859     unpackDetId(itr->first, subdet, zside, ieta, iphi, depth);
0860     if (debug) {
0861       std::cout << "DETID :" << subdet << "  IETA :" << ieta << " HIST ENTRIES :" << (itr->second)->GetEntries()
0862                 << std::endl;
0863     }
0864   }
0865 
0866   if (debug)
0867     std::cout << "Histos with " << histos_.size() << " entries\n";
0868   for (std::map<unsigned int, TH1D *>::const_iterator itr = histos_.begin(); itr != histos_.end(); ++itr, ++kount) {
0869     std::pair<double, double> result = fitMean(itr->second, 0);
0870     double factor = (inverse) ? (2. - result.first) : result.first;
0871     if (debug) {
0872       int subdet, depth, zside, ieta, iphi;
0873       unpackDetId(itr->first, subdet, zside, ieta, iphi, depth);
0874       std::cout << "DetId[" << kount << "] " << subdet << ":" << zside * ieta << ":" << depth << " Factor " << factor
0875                 << " +- " << result.second << std::endl;
0876     }
0877     if (!useMean_) {
0878       cfactors[itr->first] = std::pair<double, double>(factor, result.second);
0879       if (itr->second->GetEntries() > nMin) {
0880         kountus++;
0881         if (factor > 1)
0882           sumfactor += (1 - 1 / factor);
0883         else
0884           sumfactor += (1 - factor);
0885       }
0886     }
0887   }
0888 
0889   if (debug)
0890     std::cout << "SumW with " << SumW.size() << " entries\n";
0891   std::map<unsigned int, myEntry>::const_iterator SumWItr = SumW.begin();
0892   for (; SumWItr != SumW.end(); SumWItr++) {
0893     unsigned int detid = SumWItr->first;
0894     int subdet, depth, zside, ieta, iphi;
0895     unpackDetId(detid, subdet, zside, ieta, iphi, depth);
0896     if (debug) {
0897       std::cout << "Detid|kount|SumWi|SumFac|myId : " << subdet << ":" << zside * ieta << ":" << depth << " | "
0898                 << (SumWItr->second).kount << " | " << (SumWItr->second).fact0 << "|" << (SumWItr->second).fact1 << "|"
0899                 << (SumWItr->second).fact2 << std::endl;
0900     }
0901     double factor = (SumWItr->second).fact1 / (SumWItr->second).fact0;
0902     double dfac1 = ((SumWItr->second).fact2 / (SumWItr->second).fact0 - factor * factor);
0903     if (dfac1 < 0)
0904       dfac1 = 0;
0905     double dfac = sqrt(dfac1 / (SumWItr->second).kount);
0906     if (debug) {
0907       std::cout << "Factor " << factor << " " << dfac1 << " " << dfac << std::endl;
0908     }
0909     if (inverse)
0910       factor = 2. - factor;
0911     if (useMean_) {
0912       cfactors[detid] = std::pair<double, double>(factor, dfac);
0913       if ((SumWItr->second).kount > nMin) {
0914         kountus++;
0915         if (factor > 1)
0916           sumfactor += (1 - 1 / factor);
0917         else
0918           sumfactor += (1 - factor);
0919       }
0920     }
0921   }
0922 
0923   static const int maxch = 500;
0924   double dets[maxch], cfacs[maxch], wfacs[maxch], myId[maxch], nTrk[maxch];
0925   std::cout << "cafctors: " << cfactors.size() << ":" << maxch << std::endl;
0926   kount = 0;
0927   std::map<unsigned int, std::pair<double, double> >::const_iterator itr = cfactors.begin();
0928   const double factorMin(0.1);
0929   for (; itr != cfactors.end(); ++itr, ++kount) {
0930     unsigned int detid = itr->first;
0931     int subdet, depth, zside, ieta, iphi;
0932     unpackDetId(detid, subdet, zside, ieta, iphi, depth);
0933     double id = ieta * zside + 0.25 * (depth - 1);
0934     double factor = (itr->second).first;
0935     double dfac = (itr->second).second;
0936     if ((ieta > ietaMax) || (factor < factorMin)) {
0937       factor = 1;
0938       dfac = 0;
0939     }
0940     std::pair<double, double> cfac(factor, dfac);
0941     if (Cprev.find(detid) != Cprev.end()) {
0942       dfac /= factor;
0943       factor *= Cprev[detid].first;
0944       dfac *= factor;
0945       Cprev[detid] = std::pair<double, double>(factor, dfac);
0946       cfacs[kount] = factor;
0947     } else {
0948       Cprev[detid] = std::pair<double, double>(factor, dfac);
0949       cfacs[kount] = factor;
0950     }
0951     wfacs[kount] = factor;
0952     dets[kount] = detid;
0953     myId[kount] = id;
0954     nTrk[kount] = nTrks[detid];
0955   }
0956   if (higheta_ > 0)
0957     highEtaFactors(ietaMax, debug);
0958 
0959   std::cout << kountus << " detids out of " << kount << " have tracks > " << nMin << std::endl;
0960 
0961   char fname[50];
0962   fout->cd();
0963   TGraph *g_fac1 = new TGraph(kount, dets, cfacs);
0964   sprintf(fname, "Cfacs%d", loop);
0965   g_fac1->SetMarkerStyle(7);
0966   g_fac1->SetMarkerSize(5.0);
0967   g_fac1->Draw("AP");
0968   g_fac1->Write(fname);
0969   TGraph *g_fac2 = new TGraph(kount, dets, wfacs);
0970   sprintf(fname, "Wfacs%d", loop);
0971   g_fac2->SetMarkerStyle(7);
0972   g_fac2->SetMarkerSize(5.0);
0973   g_fac2->Draw("AP");
0974   g_fac2->Write(fname);
0975   TGraph *g_fac3 = new TGraph(kount, myId, cfacs);
0976   sprintf(fname, "CfacsVsMyId%d", loop);
0977   g_fac3->SetMarkerStyle(7);
0978   g_fac3->SetMarkerSize(5.0);
0979   g_fac3->Draw("AP");
0980   g_fac3->Write(fname);
0981   TGraph *g_fac4 = new TGraph(kount, myId, wfacs);
0982   sprintf(fname, "WfacsVsMyId%d", loop);
0983   g_fac4->SetMarkerStyle(7);
0984   g_fac4->SetMarkerSize(5.0);
0985   g_fac4->Draw("AP");
0986   g_fac4->Write(fname);
0987   TGraph *g_nTrk = new TGraph(kount, myId, nTrk);
0988   sprintf(fname, "nTrk");
0989   if (loop == 0) {
0990     g_nTrk->SetMarkerStyle(7);
0991     g_nTrk->SetMarkerSize(5.0);
0992     g_nTrk->Draw("AP");
0993     g_nTrk->Write(fname);
0994   }
0995   std::cout << "The new factors are :" << std::endl;
0996   std::map<unsigned int, std::pair<double, double> >::const_iterator CprevItr = Cprev.begin();
0997   unsigned int indx(0);
0998   for (; CprevItr != Cprev.end(); CprevItr++, indx++) {
0999     unsigned int detid = CprevItr->first;
1000     int subdet, depth, zside, ieta, iphi;
1001     unpackDetId(detid, subdet, zside, ieta, iphi, depth);
1002     std::cout << "DetId[" << indx << "] " << std::hex << detid << std::dec << "(" << ieta * zside << "," << depth
1003               << ") (nTrks:" << nTrks[detid] << ") : " << CprevItr->second.first << " +- " << CprevItr->second.second
1004               << std::endl;
1005   }
1006   double mean = (kountus > 0) ? (sumfactor / kountus) : 0;
1007   std::cout << "Mean deviation " << mean << " from 1 for " << kountus << " DetIds" << std::endl;
1008   h_cvg->SetBinContent(loop + 1, mean);
1009   if (last)
1010     h_cvg->Write("Cvg0");
1011   return mean;
1012 }
1013 
1014 void CalibTree::getDetId(double fraction, int ietaTrack, bool debug, Long64_t nmax) {
1015   if (fChain != 0) {
1016     Long64_t nbytes(0), nb(0), kprint(0);
1017     Long64_t nentryTot = fChain->GetEntriesFast();
1018     Long64_t nentries = (fraction > 0.01 && fraction < 0.99) ? (Long64_t)(fraction * nentryTot) : nentryTot;
1019     if ((nentries > nmax) && (nmax > 0))
1020       nentries = nmax;
1021 
1022     for (Long64_t jentry = 0; jentry < nentries; jentry++) {
1023       Long64_t ientry = LoadTree(jentry);
1024       if (ientry < 0)
1025         break;
1026       nb = fChain->GetEntry(jentry);
1027       nbytes += nb;
1028       if (jentry % 1000000 == 0)
1029         std::cout << "Entry " << jentry << " Run " << t_Run << " Event " << t_Event << std::endl;
1030       // Find DetIds contributing to the track

1031       bool selRun = (includeRun_ ? ((t_Run >= runlo_) && (t_Run <= runhi_)) : ((t_Run < runlo_) || (t_Run > runhi_)));
1032       bool selTrack = ((ietaTrack <= 0) || (abs(t_ieta) <= ietaTrack));
1033       if (selRun && (t_nVtx >= nvxlo_) && (t_nVtx <= nvxhi_) && selTrack) {
1034         bool isItRBX(false);
1035         if (cSelect_ != nullptr) {
1036           bool temp = cSelect_->isItRBX(t_DetIds);
1037           if (exclude_)
1038             isItRBX = temp;
1039           else
1040             isItRBX = !(temp);
1041         }
1042         ++kprint;
1043         if (!(isItRBX)) {
1044           for (unsigned int idet = 0; idet < (*t_DetIds).size(); idet++) {
1045             if (selectPhi((*t_DetIds)[idet])) {
1046               unsigned int detid = truncateId((*t_DetIds)[idet], truncateFlag_, debug);
1047               if (debug && (kprint <= 10)) {
1048                 std::cout << "DetId[" << idet << "] Original " << std::hex << (*t_DetIds)[idet] << " truncated "
1049                           << detid << std::dec;
1050               }
1051               if (std::find(detIds_.begin(), detIds_.end(), detid) == detIds_.end()) {
1052                 detIds_.push_back(detid);
1053                 if (debug && (kprint <= 10))
1054                   std::cout << " new";
1055               }
1056               if (debug && (kprint <= 10))
1057                 std::cout << std::endl;
1058             }
1059           }
1060           // Also look at the neighbouring cells if available

1061           if (t_DetIds3 != 0) {
1062             for (unsigned int idet = 0; idet < (*t_DetIds3).size(); idet++) {
1063               if (selectPhi((*t_DetIds3)[idet])) {
1064                 unsigned int detid = truncateId((*t_DetIds3)[idet], truncateFlag_, debug);
1065                 if (std::find(detIds_.begin(), detIds_.end(), detid) == detIds_.end()) {
1066                   detIds_.push_back(detid);
1067                 }
1068               }
1069             }
1070           }
1071         }
1072       }
1073     }
1074   }
1075   if (debug) {
1076     std::cout << "Total of " << detIds_.size() << " detIds" << std::endl;
1077     // The masks are defined in DataFormats/HcalDetId/interface/HcalDetId.h

1078     for (unsigned int k = 0; k < detIds_.size(); ++k) {
1079       int subdet, depth, zside, ieta, iphi;
1080       unpackDetId(detIds_[k], subdet, zside, ieta, iphi, depth);
1081       std::cout << "DetId[" << k << "] " << subdet << ":" << zside * ieta << ":" << depth << ":" << iphi << "  "
1082                 << std::hex << detIds_[k] << std::dec << std::endl;
1083     }
1084   }
1085 }
1086 
1087 void CalibTree::bookHistos(int loop, bool debug) {
1088   unsigned int k(0);
1089   for (std::map<unsigned int, TH1D *>::const_iterator itr = histos_.begin(); itr != histos_.end(); ++itr, ++k) {
1090     if (debug) {
1091       std::cout << "histos[" << k << "] " << std::hex << itr->first << std::dec << " " << itr->second;
1092       if (itr->second != 0)
1093         std::cout << " " << itr->second->GetTitle();
1094       std::cout << std::endl;
1095     }
1096     if (itr->second != 0)
1097       itr->second->Delete();
1098   }
1099 
1100   for (unsigned int k = 0; k < detIds_.size(); ++k) {
1101     char name[20], title[100];
1102     sprintf(name, "Hist%d_%d", detIds_[k], loop);
1103     int subdet, depth, zside, ieta, iphi;
1104     unpackDetId(detIds_[k], subdet, zside, ieta, iphi, depth);
1105     sprintf(title, "Correction for Subdet %d #eta %d depth %d (Loop %d)", subdet, zside * ieta, depth, loop);
1106     TH1D *hist = new TH1D(name, title, 100, 0.0, 5.0);
1107     hist->Sumw2();
1108     if (debug)
1109       std::cout << "Book Histo " << k << " " << title << std::endl;
1110     histos_[detIds_[k]] = hist;
1111   }
1112   std::cout << "Total of " << detIds_.size() << " detIds and " << histos_.size() << std::endl;
1113 }
1114 
1115 bool CalibTree::goodTrack() {
1116   bool ok(true);
1117   double cut(2.0);
1118   double pmom = (useGen_ && (t_gentrackP > 0)) ? t_gentrackP : t_p;
1119   if (sysmode_ == 1) {
1120     ok =
1121         ((t_qltyFlag) && (t_hmaxNearP < cut) && (t_eMipDR < 1.0) && (t_mindR1 > 1.0) && (pmom > 40.0) && (pmom < 60.0));
1122   } else if (sysmode_ == 2) {
1123     ok = ((t_qltyFlag) && (t_qltyPVFlag) && (t_hmaxNearP < cut) && (t_eMipDR < 1.0) && (t_mindR1 > 1.0) &&
1124           (pmom > 40.0) && (pmom < 60.0));
1125   } else if (sysmode_ == 3) {
1126     ok =
1127         ((t_selectTk) && (t_hmaxNearP < cut) && (t_eMipDR < 1.0) && (t_mindR1 > 1.0) && (pmom > 40.0) && (pmom < 60.0));
1128   } else if (sysmode_ == 4) {
1129     ok = ((t_selectTk) && (t_qltyMissFlag) && (t_hmaxNearP < 0.0) && (t_eMipDR < 1.0) && (t_mindR1 > 1.0) &&
1130           (pmom > 40.0) && (pmom < 60.0));
1131   } else if (sysmode_ == 5) {
1132     ok = ((t_selectTk) && (t_qltyMissFlag) && (t_hmaxNearP < cut) && (t_eMipDR < 0.5) && (t_mindR1 > 1.0) &&
1133           (pmom > 40.0) && (pmom < 60.0));
1134   } else if (sysmode_ == 6) {
1135     ok = ((t_selectTk) && (t_qltyMissFlag) && (t_hmaxNearP < cut) && (t_eMipDR < 2.0) && (t_mindR1 > 1.0) &&
1136           (pmom > 40.0) && (pmom < 60.0));
1137   } else if (sysmode_ == 7) {
1138     ok = ((t_selectTk) && (t_qltyMissFlag) && (t_hmaxNearP < cut) && (t_eMipDR < 1.0) && (t_mindR1 > 0.5) &&
1139           (pmom > 40.0) && (pmom < 60.0));
1140   } else {
1141     if (sysmode_ < 0) {
1142       double eta = (t_ieta > 0) ? t_ieta : -t_ieta;
1143       if (sysmode_ == -2)
1144         cut = 8.0 * exp(eta * log2by18_);
1145       else
1146         cut = 10.0;
1147     }
1148     ok = ((t_selectTk) && (t_qltyMissFlag) && (t_hmaxNearP < cut) && (t_eMipDR < 1.0) && (t_mindR1 > 1.0) &&
1149           (pmom > 40.0) && (pmom < 60.0));
1150   }
1151   return ok;
1152 }
1153 
1154 void CalibTree::writeCorrFactor(const char *corrFileName, int ietaMax) {
1155   std::ofstream myfile;
1156   myfile.open(corrFileName);
1157   if (!myfile.is_open()) {
1158     std::cout << "** ERROR: Can't open '" << corrFileName << std::endl;
1159   } else {
1160     myfile << "#" << std::setprecision(4) << std::setw(10) << "detId" << std::setw(10) << "ieta" << std::setw(10)
1161            << "depth" << std::setw(15) << "corrFactor" << std::endl;
1162     std::map<unsigned int, std::pair<double, double> >::const_iterator itr;
1163     for (itr = Cprev.begin(); itr != Cprev.end(); ++itr) {
1164       unsigned int detId = itr->first;
1165       int subdet, depth, zside, ieta, iphi;
1166       unpackDetId(detId, subdet, zside, ieta, iphi, depth);
1167       if (ieta <= ietaMax) {
1168         myfile << std::setw(10) << std::hex << detId << std::setw(10) << std::dec << zside * ieta << std::setw(10)
1169                << depth << std::setw(10) << itr->second.first << " " << std::setw(10) << itr->second.second
1170                << std::endl;
1171         std::cout << itr->second.first << ",";
1172       }
1173     }
1174     myfile.close();
1175     std::cout << std::endl;
1176   }
1177 }
1178 
1179 bool CalibTree::selectPhi(unsigned int detId) {
1180   bool flag(true);
1181   // The masks are defined in DataFormats/HcalDetId/interface/HcalDetId.h

1182   if ((phimin_ > 1) || (phimax_ < 72) || (zside_ != 0)) {
1183     int subdet, depth, zside, ieta, iphi;
1184     unpackDetId(detId, subdet, zside, ieta, iphi, depth);
1185     if (phimin_ > 1 || phimax_ < 72) {
1186       if ((iphi < phimin_) || (iphi > phimax_))
1187         flag = false;
1188     }
1189     if (zside_ != 0) {
1190       if (zside != zside_)
1191         flag = false;
1192     }
1193   }
1194   return flag;
1195 }
1196 
1197 std::pair<double, double> CalibTree::fitMean(TH1D *hist, int mode) {
1198   std::pair<double, double> results = std::pair<double, double>(1, 0);
1199   if (hist != 0) {
1200     double mean = hist->GetMean(), rms = hist->GetRMS();
1201     double LowEdge(0.7), HighEdge(1.3);
1202     char option[20];
1203     if (mode == 1) {
1204       LowEdge = mean - 1.5 * rms;
1205       HighEdge = mean + 1.5 * rms;
1206       int nbin = hist->GetNbinsX();
1207       if (LowEdge < hist->GetBinLowEdge(1))
1208         LowEdge = hist->GetBinLowEdge(1);
1209       if (HighEdge > hist->GetBinLowEdge(nbin) + hist->GetBinWidth(nbin))
1210         HighEdge = hist->GetBinLowEdge(nbin) + hist->GetBinWidth(nbin);
1211     }
1212     if (hist->GetEntries() > 100)
1213       sprintf(option, "+QRLS");
1214     else
1215       sprintf(option, "+QRWLS");
1216     double value(mean);
1217     double error = rms / sqrt(hist->GetEntries());
1218     if (hist->GetEntries() > 20) {
1219       TFitResultPtr Fit = hist->Fit("gaus", option, "", LowEdge, HighEdge);
1220       value = Fit->Value(1);
1221       error = Fit->FitResult::Error(1);
1222       /*

1223       LowEdge  = value - 1.5*error;

1224       HighEdge = value + 1.5*error;

1225       value    = Fit->Value(1);

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

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

1228       */
1229     }
1230     results = std::pair<double, double>(value, error);
1231   }
1232   return results;
1233 }
1234 
1235 void CalibTree::makeplots(
1236     double rmin, double rmax, int ietaMax, bool useweight, double fraction, bool debug, Long64_t nmax) {
1237   if (fChain == 0)
1238     return;
1239   Long64_t nentryTot = fChain->GetEntriesFast();
1240   Long64_t nentries = (fraction > 0.01 && fraction < 0.99) ? (Long64_t)(fraction * nentryTot) : nentryTot;
1241   if ((nentries > nmax) && (nmax > 0))
1242     nentries = nmax;
1243 
1244   // Book the histograms

1245   std::map<int, std::pair<TH1D *, TH1D *> > histos;
1246   for (int ieta = -ietaMax; ieta <= ietaMax; ++ieta) {
1247     char name[20], title[100];
1248     sprintf(name, "begin%d", ieta);
1249     if (ieta == 0)
1250       sprintf(title, "Ratio at start");
1251     else
1252       sprintf(title, "Ratio at start for i#eta=%d", ieta);
1253     TH1D *h1 = new TH1D(name, title, 50, rmin, rmax);
1254     h1->Sumw2();
1255     sprintf(name, "end%d", ieta);
1256     if (ieta == 0)
1257       sprintf(title, "Ratio at the end");
1258     else
1259       sprintf(title, "Ratio at the end for i#eta=%d", ieta);
1260     TH1D *h2 = new TH1D(name, title, 50, rmin, rmax);
1261     h2->Sumw2();
1262     histos[ieta] = std::pair<TH1D *, TH1D *>(h1, h2);
1263   }
1264   //Fill the histograms

1265   Long64_t nbytes(0), nb(0);
1266   for (Long64_t jentry = 0; jentry < nentries; jentry++) {
1267     Long64_t ientry = LoadTree(jentry);
1268     if (ientry < 0)
1269       break;
1270     if (std::find(entries.begin(), entries.end(), jentry) != entries.end())
1271       break;
1272     nb = fChain->GetEntry(jentry);
1273     nbytes += nb;
1274     if (goodTrack()) {
1275       double pmom = (useGen_ && (t_gentrackP > 0)) ? t_gentrackP : t_p;
1276       CalibTree::energyCalor en1 = energyHcal(pmom, jentry, false);
1277       CalibTree::energyCalor en2 = energyHcal(pmom, jentry, true);
1278       if ((en1.ehcal > 0.001) && (en2.ehcal > 0.001)) {
1279         double evWt = (useweight) ? t_EventWeight : 1.0;
1280         double ratioi = en1.ehcal / (pmom - t_eMipDR);
1281         double ratiof = en2.ehcal / (pmom - t_eMipDR);
1282         if (t_ieta >= -ietaMax && t_ieta <= ietaMax && t_ieta != 0) {
1283           if (ratioi >= rmin && ratioi <= rmax) {
1284             histos[0].first->Fill(ratioi, evWt);
1285             histos[t_ieta].first->Fill(ratioi, evWt);
1286           }
1287           if (ratiof >= rmin && ratiof <= rmax) {
1288             histos[0].second->Fill(ratiof, evWt);
1289             histos[t_ieta].second->Fill(ratiof, evWt);
1290           }
1291         }
1292       }
1293     }
1294   }
1295 
1296   //Fit the histograms

1297   TH1D *hbef1 = new TH1D("Eta1Bf", "Mean vs i#eta", 2 * ietaMax, -ietaMax, ietaMax);
1298   TH1D *hbef2 = new TH1D("Eta2Bf", "Median vs i#eta", 2 * ietaMax, -ietaMax, ietaMax);
1299   TH1D *haft1 = new TH1D("Eta1Af", "Mean vs i#eta", 2 * ietaMax, -ietaMax, ietaMax);
1300   TH1D *haft2 = new TH1D("Eta2Af", "Median vs i#eta", 2 * ietaMax, -ietaMax, ietaMax);
1301   for (int ieta = -ietaMax; ieta <= ietaMax; ++ieta) {
1302     int bin = (ieta < 0) ? (ieta + ietaMax + 1) : (ieta + ietaMax);
1303     TH1D *h1 = histos[ieta].first;
1304     double mean1 = h1->GetMean();
1305     double err1 = h1->GetMeanError();
1306     std::pair<double, double> fit1 = fitMean(h1, 1);
1307     if (debug) {
1308       std::cout << ieta << " " << h1->GetName() << " " << mean1 << " +- " << err1 << " and " << fit1.first << " +- "
1309                 << fit1.second << std::endl;
1310     }
1311     if (ieta != 0) {
1312       hbef1->SetBinContent(bin, mean1);
1313       hbef1->SetBinError(bin, err1);
1314       hbef2->SetBinContent(bin, fit1.first);
1315       hbef2->SetBinError(bin, fit1.second);
1316     }
1317     h1->Write();
1318     TH1D *h2 = histos[ieta].second;
1319     double mean2 = h2->GetMean();
1320     double err2 = h2->GetMeanError();
1321     std::pair<double, double> fit2 = fitMean(h2, 1);
1322     if (debug) {
1323       std::cout << ieta << " " << h2->GetName() << " " << mean2 << " +- " << err2 << " and " << fit2.first << " +- "
1324                 << fit2.second << std::endl;
1325     }
1326     if (ieta != 0) {
1327       haft1->SetBinContent(bin, mean2);
1328       haft1->SetBinError(bin, err2);
1329       haft2->SetBinContent(bin, fit2.first);
1330       haft2->SetBinError(bin, fit2.second);
1331     }
1332     h2->Write();
1333   }
1334   fitPol0(hbef1, debug);
1335   fitPol0(hbef2, debug);
1336   fitPol0(haft1, debug);
1337   fitPol0(haft2, debug);
1338 }
1339 
1340 void CalibTree::fitPol0(TH1D *hist, bool debug) {
1341   hist->GetXaxis()->SetTitle("i#eta");
1342   hist->GetYaxis()->SetTitle("<E_{HCAL}/(p-E_{ECAL})>");
1343   hist->GetYaxis()->SetRangeUser(0.4, 1.6);
1344   TFitResultPtr Fit = hist->Fit("pol0", "+QRWLS");
1345   if (debug) {
1346     std::cout << "Fit to Pol0 to " << hist->GetTitle() << ": " << Fit->Value(0) << " +- " << Fit->FitResult::Error(0)
1347               << std::endl;
1348   }
1349   hist->Write();
1350 }
1351 
1352 void CalibTree::highEtaFactors(int ietaMax, bool debug) {
1353   std::map<unsigned int, std::pair<double, double> >::const_iterator itr;
1354   std::pair<double, double> cfacp, cfacn;
1355   cfacp = cfacn = std::pair<double, double>(1.0, 0.0);
1356   for (itr = Cprev.begin(); itr != Cprev.end(); ++itr) {
1357     unsigned int detid = itr->first;
1358     int subdet, depth, zside, ieta, iphi;
1359     unpackDetId(detid, subdet, zside, ieta, iphi, depth);
1360     if ((ieta == ietaMax) && (depth == 1)) {
1361       if (zside > 0)
1362         cfacp = itr->second;
1363       else
1364         cfacn = itr->second;
1365     }
1366   }
1367   if (debug) {
1368     std::cout << "Correction factor for (" << ietaMax << ",1) = " << cfacp.first << ":" << cfacp.second << " and ("
1369               << -ietaMax << ",1) = " << cfacn.first << ":" << cfacn.second << std::endl;
1370   }
1371   for (itr = Cprev.begin(); itr != Cprev.end(); ++itr) {
1372     unsigned int detid = itr->first;
1373     int subdet, depth, zside, ieta, iphi;
1374     unpackDetId(detid, subdet, zside, ieta, iphi, depth);
1375     if (ieta > ietaMax) {
1376       Cprev[detid] = (zside > 0) ? cfacp : cfacn;
1377       if (debug) {
1378         std::cout << "Set correction factor for (" << zside * ieta << "," << depth << ") = " << Cprev[detid].first
1379                   << ":" << Cprev[detid].second << std::endl;
1380       }
1381     }
1382   }
1383 }
1384 
1385 CalibTree::energyCalor CalibTree::energyHcal(double pmom, const Long64_t &entry, bool final) {
1386   double etot = t_eHcal;
1387   double etot2 = t_eHcal;
1388   double ediff = (t_eHcal30 - t_eHcal10);
1389   if (final) {
1390     etot = etot2 = 0;
1391     for (unsigned int idet = 0; idet < (*t_DetIds).size(); idet++) {
1392       if (selectPhi((*t_DetIds)[idet])) {
1393         unsigned int id = (*t_DetIds)[idet];
1394         double hitEn(0);
1395         unsigned int detid = truncateId(id, truncateFlag_, false);
1396         if (Cprev.find(detid) != Cprev.end())
1397           hitEn = Cprev[detid].first * (*t_HitEnergies)[idet];
1398         else
1399           hitEn = (*t_HitEnergies)[idet];
1400         if (cFactor_)
1401           hitEn *= cFactor_->getCorr(t_Run, id);
1402         etot += hitEn;
1403         etot2 += ((*t_HitEnergies)[idet]);
1404       }
1405     }
1406     // Now the outer cone

1407     double etot1(0), etot3(0);
1408     if (t_DetIds1 != 0 && t_DetIds3 != 0) {
1409       for (unsigned int idet = 0; idet < (*t_DetIds1).size(); idet++) {
1410         if (selectPhi((*t_DetIds1)[idet])) {
1411           unsigned int id = (*t_DetIds1)[idet];
1412           unsigned int detid = truncateId(id, truncateFlag_, false);
1413           double hitEn(0);
1414           if (Cprev.find(detid) != Cprev.end())
1415             hitEn = Cprev[detid].first * (*t_HitEnergies1)[idet];
1416           else
1417             hitEn = (*t_HitEnergies1)[idet];
1418           if (cFactor_)
1419             hitEn *= cFactor_->getCorr(t_Run, id);
1420           etot1 += hitEn;
1421         }
1422       }
1423       for (unsigned int idet = 0; idet < (*t_DetIds3).size(); idet++) {
1424         if (selectPhi((*t_DetIds3)[idet])) {
1425           unsigned int id = (*t_DetIds3)[idet];
1426           unsigned int detid = truncateId(id, truncateFlag_, false);
1427           double hitEn(0);
1428           if (Cprev.find(detid) != Cprev.end())
1429             hitEn = Cprev[detid].first * (*t_HitEnergies3)[idet];
1430           else
1431             hitEn = (*t_HitEnergies3)[idet];
1432           if (cFactor_)
1433             hitEn *= cFactor_->getCorr(t_Run, id);
1434           etot3 += hitEn;
1435         }
1436       }
1437     }
1438     ediff = etot3 - etot1;
1439   }
1440   // PU correction only for loose isolation cut

1441   double ehcal = (((rcorForm_ == 3) && (cFactor_ != nullptr))
1442                       ? (etot * cFactor_->getCorr(entry))
1443                       : ((puCorr_ == 0) ? etot
1444                                         : ((puCorr_ < 0) ? (etot * puFactor(-puCorr_, t_ieta, pmom, etot, ediff))
1445                                                          : puFactorRho(puCorr_, t_ieta, t_rhoh, etot))));
1446   return CalibTree::energyCalor(etot, etot2, ehcal);
1447 }
[ 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!