Back to home page

Project CMSSW displayed by LXR

 
 

    


File indexing completed on 2023-10-25 09:34:30

0001 //////////////////////////////////////////////////////////////////////////////
0002 // Usage:
0003 // .L CalibSort.C+g (for the tree "CalibTree")
0004 //  CalibSort c1(fname, dirname, prefix, flag, mipCut);
0005 //  c1.Loop("events.txt");
0006 //  findDuplicate(infile, outfile, debug)
0007 //
0008 // .L CalibSort.C+g (for the tree "EventInfo")
0009 //  CalibSortEvent c2(fname, dirname, prefix, append);
0010 //  c2.Loop();
0011 //  findDuplicateEvent(infile, outfile, debug)
0012 //
0013 //  This will prepare a list of dupliate entries from comombined data sets
0014 //
0015 //   where:
0016 //
0017 //   fname    (const char*)    = file name of the input ROOT tree
0018 //                               or name of the file containing a list of
0019 //                               file names of input ROOT trees
0020 //   dirname  (std::string)    = name of the directory where Tree resides
0021 //                               (default "HcalIsoTrkAnalyzer")
0022 //   prefix  (std::string)     = String to be added to the name
0023 //                               (usually a 4 character string; default="")
0024 //   flag    (int)             = 2 digit integer (do) with specific control
0025 //                               information (d = 0/1 for debug off/on;
0026 //                               o = 0/1 for creating the output "events.txt"
0027 //                               file in append/output mode. Default = 0
0028 //   mipCut   (double)         = cut off on ECAL energy (default = 2.0)
0029 //
0030 //   infile   (std::string)    = name of the input file containing run, event,
0031 //                               information (created by CalibSort)
0032 //   outfile  (std::string)    = name of the file containing the entry numbers
0033 //                               of duplicate events
0034 //   debug    (bool)           = true/false for getting or not debug printouts
0035 //                               (default = false)
0036 //   append   (bool)           = flag used for opening o/p file. Default=false
0037 //
0038 //
0039 // .L CalibSort.C+g
0040 //  combine(fname1, fname2, fname0, outfile, debug)
0041 //
0042 //  Combines the 2 files created by isotrackApplyRegressor.py from
0043 //  Models 1 (endcap) and 2 (barrel)
0044 //
0045 //   fname1/2 (std::string)    = file names for endcap/barrel correction
0046 //                               factors
0047 //   fname0   (std::string)    = output combined file
0048 //   outfile  (std::string)    = root file name for storing histograms if
0049 //                               a vlaid name is given (default = "")
0050 //   debug    (int)            = verbosity flag (default = 1)
0051 //
0052 // .L CalibSort.C+g
0053 //  plotCombine(infile, flag, drawStatBox, save)
0054 //
0055 //   infile       (const char*)= name of input ROOT file
0056 //   flag         (int)        = indicates the source of the file.
0057 //                               -1 : created by the method combine
0058 //                               >=0 : created by CalibPlotCombine with
0059 //                                     the same meaning as in that method
0060 //   drawStatBox  (bool)       = flag to draw the statistics box (true)
0061 //   save         (bool)       = save the plots (true)
0062 //
0063 // .L CalibSort.C+g (for the tree "CalibTree")
0064 //  CalibPlotCombine c1(fname, corrFileName, dirname, flag);
0065 //  c1.Loop()
0066 //  c1.savePlot(histFileName,append)
0067 //
0068 //   fname        (const char*)= file name of the input ROOT tree
0069 //                               or name of the file containing a list of
0070 //                               file names of input ROOT trees
0071 //   corrFileName (const char*)= name of the text file having the correction
0072 //                               factors as a function of entry number
0073 //   dirname      (std::string)= name of the directory where Tree resides
0074 //                               ("HcalIsoTrkAnalyzer")
0075 //   flag         (int)        = two digit number "to". o: specifies if
0076 //                               charge isolation to be applied (1) or not
0077 //                               (0); t: if all tracks to be included (0)
0078 //                               or only with momentum 40-60 GeV (1)
0079 //                               Default (10)
0080 //   histFileName (std::string)= name of the histogram file
0081 //   append       (bool)       = true/false if the histogram file to be opened
0082 //                               in append/output mode
0083 //
0084 // .L CalibSort.C+g (for the o/p of isotrackRootTreeMaker.py)
0085 //  CalibFitPU c1(fname);
0086 //  c1.Loop(extractPUparams, fileName, dumpEntry);
0087 //
0088 //   fname        (const char*)= file name of the input ROOT tree which is
0089 //                               output of isotrackRootTreeMaker.py
0090 //   extractPUparams (bool)    = flag to see if extraction is needed
0091 //   filename     (std::string)= root name of all files to be created/read
0092 //                               (_par2d.txt, _parProf.txt, _parGraph.txt
0093 //                               will be names of files of parameters from
0094 //                               2D, profile, graphs; .root for storing all
0095 //                               histograms created)
0096 //   dumpEntry    (int)          Entries to be dumped.  Default (0)
0097 //
0098 // .L CalibSort.C+g (for merging information from PU and NoPU files)
0099 //  CalibMerge c1;
0100 //  c1.LoopNoPU(fileNameNoPU, dirnm);
0101 //  c1.LoopPU(fileNamePU, fname, dirnm);
0102 //  c1.close();
0103 //
0104 //   fileNameNoPU (const char*)= file name of the input ROOT tree for the
0105 //                               no PileUP Monte Carlo
0106 //   fileNamePU   (const char*)= file name of the input ROOT tree for the
0107 //                               PileUP Monte Carlo
0108 //   fname        (const char*)= file name of the output ROOT tree where the
0109 //                               combined information is written
0110 //   dirnm        (std::string) = name of the directory where Tree resides
0111 //                               (default "HcalIsoTrkAnalyzer")
0112 //
0113 // .L CalibSort.C+g
0114 //  combineML(inputFileList, outfile);
0115 //
0116 //  Combines the ML values otained from the analysis of muon analysis to
0117 //  determine depth dependent correction factors
0118 //
0119 //   inputFileList (const char*) = file containing filenames having the ML
0120 //                                 values for a given depth
0121 //   outfile       (const char*) = name of the output file where the
0122 //                                 depth dependent correction factors
0123 //                                 will be stored
0124 //   Example of a inputFileList:
0125 //      depth Name of the file
0126 //          1 ml_values_depth1.txt
0127 //          2 ml_values_depth2.txt
0128 //          3 ml_values_depth3.txt
0129 //          4 ml_values_depth4.txt
0130 //   where each file conatins
0131 //   (4 quantities per line with no tab separating each item):
0132 //   #ieta depth  ml    uncertainity-in-ml
0133 //
0134 // .L CalibSort.C+g
0135 //  calCorrCombine(inFileCorr, inFileDepth, outFileCorr, truncateFlag, etaMax);
0136 //
0137 //  Combines the correction factors obtained from CalibTree and the depth
0138 //  dependent correction factors used in the CalibTree command
0139 //
0140 //   inputFileCorr  (const char*) = file containing correction factors obtained
0141 //                                  from the CalibTree
0142 //   inputFileDepth (const char*) = file containing depth dependent correction
0143 //                                  factors used in the CalibTree step
0144 //   outfileCorr    (const char*) = name of the output file where the combined
0145 //                                  correction factors
0146 //   truncateFlag   (int)         = Truncate flag used in the CalibTree step
0147 //   etaMax         (int)         = Maximum eta value
0148 //
0149 //////////////////////////////////////////////////////////////////////////////
0150 
0151 #include <TCanvas.h>
0152 #include <TChain.h>
0153 #include <TF1.h>
0154 #include <TFile.h>
0155 #include <TFitResult.h>
0156 #include <TFitResultPtr.h>
0157 #include <TGraph.h>
0158 #include <TH1D.h>
0159 #include <TH2D.h>
0160 #include <TLegend.h>
0161 #include <TPaveStats.h>
0162 #include <TPaveText.h>
0163 #include <TProfile.h>
0164 #include <TROOT.h>
0165 #include <TStyle.h>
0166 
0167 #include <algorithm>
0168 #include <fstream>
0169 #include <iomanip>
0170 #include <iostream>
0171 #include <map>
0172 #include <sstream>
0173 #include <string>
0174 #include <vector>
0175 
0176 #include "CalibCorr.C"
0177 
0178 struct record {
0179   record(int ser = 0, int ent = 0, int r = 0, int ev = 0, int ie = 0, double p = 0)
0180       : serial_(ser), entry_(ent), run_(r), event_(ev), ieta_(ie), p_(p) {}
0181 
0182   int serial_, entry_, run_, event_, ieta_;
0183   double p_;
0184 };
0185 
0186 struct recordLess {
0187   bool operator()(const record &a, const record &b) {
0188     return ((a.run_ < b.run_) || ((a.run_ == b.run_) && (a.event_ < b.event_)) ||
0189             ((a.run_ == b.run_) && (a.event_ == b.event_) && (a.ieta_ < b.ieta_)));
0190   }
0191 };
0192 
0193 struct recordEvent {
0194   recordEvent(unsigned int ser = 0, unsigned int ent = 0, unsigned int r = 0, unsigned int ev = 0)
0195       : serial_(ser), entry_(ent), run_(r), event_(ev) {}
0196 
0197   unsigned int serial_, entry_, run_, event_;
0198 };
0199 
0200 struct recordEventLess {
0201   bool operator()(const recordEvent &a, const recordEvent &b) {
0202     return ((a.run_ < b.run_) || ((a.run_ == b.run_) && (a.event_ < b.event_)));
0203   }
0204 };
0205 
0206 struct listML {
0207   listML(int dep = 0, double ml = 0, double dml = 0) : depth_(dep), ml_(ml), dml_(dml) {}
0208 
0209   int depth_;
0210   double ml_, dml_;
0211 };
0212 
0213 class CalibSort {
0214 public:
0215   CalibSort(const char *fname,
0216             std::string dirname = "HcalIsoTrkAnalyzer",
0217             std::string prefix = "",
0218             bool allEvent = false,
0219             int flag = 0,
0220             double mipCut = 2.0);
0221   virtual ~CalibSort();
0222   virtual Int_t Cut(Long64_t entry);
0223   virtual Int_t GetEntry(Long64_t entry);
0224   virtual Long64_t LoadTree(Long64_t entry);
0225   virtual void Init(TChain *);
0226   virtual void Loop(const char *);
0227   virtual Bool_t Notify();
0228   virtual void Show(Long64_t entry = -1);
0229 
0230 private:
0231   TChain *fChain;  //!pointer to the analyzed TTree or TChain
0232   Int_t fCurrent;  //!current Tree number in a TChain
0233 
0234   // Declaration of leaf types
0235   Int_t t_Run;
0236   Int_t t_Event;
0237   Int_t t_DataType;
0238   Int_t t_ieta;
0239   Int_t t_iphi;
0240   Double_t t_EventWeight;
0241   Int_t t_nVtx;
0242   Int_t t_nTrk;
0243   Int_t t_goodPV;
0244   Double_t t_l1pt;
0245   Double_t t_l1eta;
0246   Double_t t_l1phi;
0247   Double_t t_l3pt;
0248   Double_t t_l3eta;
0249   Double_t t_l3phi;
0250   Double_t t_p;
0251   Double_t t_pt;
0252   Double_t t_phi;
0253   Double_t t_mindR1;
0254   Double_t t_mindR2;
0255   Double_t t_eMipDR;
0256   Double_t t_eHcal;
0257   Double_t t_eHcal10;
0258   Double_t t_eHcal30;
0259   Double_t t_hmaxNearP;
0260   Double_t t_rhoh;
0261   Bool_t t_selectTk;
0262   Bool_t t_qltyFlag;
0263   Bool_t t_qltyMissFlag;
0264   Bool_t t_qltyPVFlag;
0265   Double_t t_gentrackP;
0266   std::vector<unsigned int> *t_DetIds;
0267   std::vector<double> *t_HitEnergies;
0268   std::vector<bool> *t_trgbits;
0269   std::vector<unsigned int> *t_DetIds1;
0270   std::vector<unsigned int> *t_DetIds3;
0271   std::vector<double> *t_HitEnergies1;
0272   std::vector<double> *t_HitEnergies3;
0273 
0274   // List of branches
0275   TBranch *b_t_Run;           //!
0276   TBranch *b_t_Event;         //!
0277   TBranch *b_t_DataType;      //!
0278   TBranch *b_t_ieta;          //!
0279   TBranch *b_t_iphi;          //!
0280   TBranch *b_t_EventWeight;   //!
0281   TBranch *b_t_nVtx;          //!
0282   TBranch *b_t_nTrk;          //!
0283   TBranch *b_t_goodPV;        //!
0284   TBranch *b_t_l1pt;          //!
0285   TBranch *b_t_l1eta;         //!
0286   TBranch *b_t_l1phi;         //!
0287   TBranch *b_t_l3pt;          //!
0288   TBranch *b_t_l3eta;         //!
0289   TBranch *b_t_l3phi;         //!
0290   TBranch *b_t_p;             //!
0291   TBranch *b_t_pt;            //!
0292   TBranch *b_t_phi;           //!
0293   TBranch *b_t_mindR1;        //!
0294   TBranch *b_t_mindR2;        //!
0295   TBranch *b_t_eMipDR;        //!
0296   TBranch *b_t_eHcal;         //!
0297   TBranch *b_t_eHcal10;       //!
0298   TBranch *b_t_eHcal30;       //!
0299   TBranch *b_t_hmaxNearP;     //!
0300   TBranch *b_t_rhoh;          //!
0301   TBranch *b_t_selectTk;      //!
0302   TBranch *b_t_qltyFlag;      //!
0303   TBranch *b_t_qltyMissFlag;  //!
0304   TBranch *b_t_qltyPVFlag;    //!
0305   TBranch *b_t_gentrackP;     //!
0306   TBranch *b_t_DetIds;        //!
0307   TBranch *b_t_HitEnergies;   //!
0308   TBranch *b_t_trgbits;       //!
0309   TBranch *b_t_DetIds1;       //!
0310   TBranch *b_t_DetIds3;       //!
0311   TBranch *b_t_HitEnergies1;  //!
0312   TBranch *b_t_HitEnergies3;  //!
0313 
0314   std::string fname_, dirnm_, prefix_;
0315   bool allEvent_;
0316   int flag_;
0317   double mipCut_;
0318 };
0319 
0320 CalibSort::CalibSort(const char *fname, std::string dirnm, std::string prefix, bool allEvent, int flag, double mipCut)
0321     : fname_(fname), dirnm_(dirnm), prefix_(prefix), allEvent_(allEvent), flag_(flag), mipCut_(mipCut) {
0322   char treeName[400];
0323   sprintf(treeName, "%s/CalibTree", dirnm.c_str());
0324   TChain *chain = new TChain(treeName);
0325   std::cout << "Create a chain for " << treeName << " from " << fname << std::endl;
0326   if (!fillChain(chain, fname)) {
0327     std::cout << "*****No valid tree chain can be obtained*****" << std::endl;
0328   } else {
0329     std::cout << "Proceed with a tree chain with " << chain->GetEntries() << " entries" << std::endl;
0330     Init(chain);
0331   }
0332 }
0333 
0334 CalibSort::~CalibSort() {
0335   if (!fChain)
0336     return;
0337   delete fChain->GetCurrentFile();
0338 }
0339 
0340 Int_t CalibSort::GetEntry(Long64_t entry) {
0341   // Read contents of entry.
0342   if (!fChain)
0343     return 0;
0344   return fChain->GetEntry(entry);
0345 }
0346 
0347 Long64_t CalibSort::LoadTree(Long64_t entry) {
0348   // Set the environment to read one entry
0349   if (!fChain)
0350     return -5;
0351   Long64_t centry = fChain->LoadTree(entry);
0352   if (centry < 0)
0353     return centry;
0354   if (!fChain->InheritsFrom(TChain::Class()))
0355     return centry;
0356   TChain *chain = (TChain *)fChain;
0357   if (chain->GetTreeNumber() != fCurrent) {
0358     fCurrent = chain->GetTreeNumber();
0359     Notify();
0360   }
0361   return centry;
0362 }
0363 
0364 void CalibSort::Init(TChain *tree) {
0365   // The Init() function is called when the selector needs to initialize
0366   // a new tree or chain. Typically here the branch addresses and branch
0367   // pointers of the tree will be set.
0368   // It is normally not necessary to make changes to the generated
0369   // code, but the routine can be extended by the user if needed.
0370   // Init() will be called many times when running on PROOF
0371   // (once per file to be processed).
0372 
0373   // Set object pointer
0374   t_DetIds = 0;
0375   t_DetIds1 = 0;
0376   t_DetIds3 = 0;
0377   t_HitEnergies = 0;
0378   t_HitEnergies1 = 0;
0379   t_HitEnergies3 = 0;
0380   t_trgbits = 0;
0381   // Set branch addresses and branch pointers
0382   fChain = tree;
0383   fCurrent = -1;
0384   if (!tree)
0385     return;
0386   fChain->SetMakeClass(1);
0387 
0388   fChain->SetBranchAddress("t_Run", &t_Run, &b_t_Run);
0389   fChain->SetBranchAddress("t_Event", &t_Event, &b_t_Event);
0390   fChain->SetBranchAddress("t_DataType", &t_DataType, &b_t_DataType);
0391   fChain->SetBranchAddress("t_ieta", &t_ieta, &b_t_ieta);
0392   fChain->SetBranchAddress("t_iphi", &t_iphi, &b_t_iphi);
0393   fChain->SetBranchAddress("t_EventWeight", &t_EventWeight, &b_t_EventWeight);
0394   fChain->SetBranchAddress("t_nVtx", &t_nVtx, &b_t_nVtx);
0395   fChain->SetBranchAddress("t_nTrk", &t_nTrk, &b_t_nTrk);
0396   fChain->SetBranchAddress("t_goodPV", &t_goodPV, &b_t_goodPV);
0397   fChain->SetBranchAddress("t_l1pt", &t_l1pt, &b_t_l1pt);
0398   fChain->SetBranchAddress("t_l1eta", &t_l1eta, &b_t_l1eta);
0399   fChain->SetBranchAddress("t_l1phi", &t_l1phi, &b_t_l1phi);
0400   fChain->SetBranchAddress("t_l3pt", &t_l3pt, &b_t_l3pt);
0401   fChain->SetBranchAddress("t_l3eta", &t_l3eta, &b_t_l3eta);
0402   fChain->SetBranchAddress("t_l3phi", &t_l3phi, &b_t_l3phi);
0403   fChain->SetBranchAddress("t_p", &t_p, &b_t_p);
0404   fChain->SetBranchAddress("t_pt", &t_pt, &b_t_pt);
0405   fChain->SetBranchAddress("t_phi", &t_phi, &b_t_phi);
0406   fChain->SetBranchAddress("t_mindR1", &t_mindR1, &b_t_mindR1);
0407   fChain->SetBranchAddress("t_mindR2", &t_mindR2, &b_t_mindR2);
0408   fChain->SetBranchAddress("t_eMipDR", &t_eMipDR, &b_t_eMipDR);
0409   fChain->SetBranchAddress("t_eHcal", &t_eHcal, &b_t_eHcal);
0410   fChain->SetBranchAddress("t_eHcal10", &t_eHcal10, &b_t_eHcal10);
0411   fChain->SetBranchAddress("t_eHcal30", &t_eHcal30, &b_t_eHcal30);
0412   fChain->SetBranchAddress("t_hmaxNearP", &t_hmaxNearP, &b_t_hmaxNearP);
0413   fChain->SetBranchAddress("t_rhoh", &t_rhoh, &b_t_rhoh);
0414   fChain->SetBranchAddress("t_selectTk", &t_selectTk, &b_t_selectTk);
0415   fChain->SetBranchAddress("t_qltyFlag", &t_qltyFlag, &b_t_qltyFlag);
0416   fChain->SetBranchAddress("t_qltyMissFlag", &t_qltyMissFlag, &b_t_qltyMissFlag);
0417   fChain->SetBranchAddress("t_qltyPVFlag", &t_qltyPVFlag, &b_t_qltyPVFlag);
0418   fChain->SetBranchAddress("t_gentrackP", &t_gentrackP, &b_t_gentrackP);
0419   fChain->SetBranchAddress("t_DetIds", &t_DetIds, &b_t_DetIds);
0420   fChain->SetBranchAddress("t_HitEnergies", &t_HitEnergies, &b_t_HitEnergies);
0421   fChain->SetBranchAddress("t_trgbits", &t_trgbits, &b_t_trgbits);
0422   fChain->SetBranchAddress("t_DetIds1", &t_DetIds1, &b_t_DetIds1);
0423   fChain->SetBranchAddress("t_DetIds3", &t_DetIds3, &b_t_DetIds3);
0424   fChain->SetBranchAddress("t_HitEnergies1", &t_HitEnergies1, &b_t_HitEnergies1);
0425   fChain->SetBranchAddress("t_HitEnergies3", &t_HitEnergies3, &b_t_HitEnergies3);
0426 
0427   Notify();
0428 }
0429 
0430 Bool_t CalibSort::Notify() {
0431   // The Notify() function is called when a new file is opened. This
0432   // can be either for a new TTree in a TChain or when when a new TTree
0433   // is started when using PROOF. It is normally not necessary to make changes
0434   // to the generated code, but the routine can be extended by the
0435   // user if needed. The return value is currently not used.
0436 
0437   return kTRUE;
0438 }
0439 
0440 void CalibSort::Show(Long64_t entry) {
0441   // Print contents of entry.
0442   // If entry is not specified, print current entry
0443   if (!fChain)
0444     return;
0445   fChain->Show(entry);
0446 }
0447 
0448 Int_t CalibSort::Cut(Long64_t) {
0449   // This function may be called from Loop.
0450   // returns  1 if entry is accepted.
0451   // returns -1 otherwise.
0452   return 1;
0453 }
0454 
0455 void CalibSort::Loop(const char *outFile) {
0456   //   In a ROOT session, you can do:
0457   //      Root > .L CalibSort.C
0458   //      Root > CalibSort t
0459   //      Root > t.GetEntry(12); // Fill t data members with entry number 12
0460   //      Root > t.Show();       // Show values of entry 12
0461   //      Root > t.Show(16);     // Read and show values of entry 16
0462   //      Root > t.Loop();       // Loop on all entries
0463   //
0464 
0465   //   This is the loop skeleton where:
0466   //      jentry is the global entry number in the chain
0467   //      ientry is the entry number in the current Tree
0468   //   Note that the argument to GetEntry must be:
0469   //      jentry for TChain::GetEntry
0470   //      ientry for TTree::GetEntry and TBranch::GetEntry
0471   //
0472   //       To read only selected branches, Insert statements like:
0473   // METHOD1:
0474   //    fChain->SetBranchStatus("*",0);  // disable all branches
0475   //    fChain->SetBranchStatus("branchname",1);  // activate branchname
0476   // METHOD2: replace line
0477   //    fChain->GetEntry(jentry);       //read all branches
0478   //by  b_branchname->GetEntry(ientry); //read only this branch
0479   if (fChain == 0)
0480     return;
0481 
0482   std::ofstream fileout;
0483   if ((flag_ % 10) == 1) {
0484     fileout.open(outFile, std::ofstream::out);
0485     std::cout << "Opens " << outFile << " in output mode" << std::endl;
0486   } else {
0487     fileout.open(outFile, std::ofstream::app);
0488     std::cout << "Opens " << outFile << " in append mode" << std::endl;
0489   }
0490   if (!allEvent_)
0491     fileout << "Input file: " << fname_ << " Directory: " << dirnm_ << " Prefix: " << prefix_ << std::endl;
0492   Int_t runLow(99999999), runHigh(0);
0493   Long64_t nbytes(0), nb(0), good(0);
0494   Long64_t nentries = fChain->GetEntriesFast();
0495   for (Long64_t jentry = 0; jentry < nentries; jentry++) {
0496     Long64_t ientry = LoadTree(jentry);
0497     if (ientry < 0)
0498       break;
0499     nb = fChain->GetEntry(jentry);
0500     nbytes += nb;
0501     if (t_Run > 200000 && t_Run < 800000) {
0502       if (t_Run < runLow)
0503         runLow = t_Run;
0504       if (t_Run > runHigh)
0505         runHigh = t_Run;
0506     }
0507     double cut = (t_p > 20) ? 10.0 : 0.0;
0508     if ((flag_ / 10) % 10 > 0)
0509       std::cout << "Entry " << jentry << " p " << t_p << " Cuts " << t_qltyFlag << "|" << t_selectTk << "|"
0510                 << (t_hmaxNearP < cut) << "|" << (t_eMipDR < mipCut_) << std::endl;
0511     if ((t_qltyFlag && t_selectTk && (t_hmaxNearP < cut) && (t_eMipDR < mipCut_)) || allEvent_) {
0512       good++;
0513       fileout << good << " " << jentry << " " << t_Run << " " << t_Event << " " << t_ieta << " " << t_p << std::endl;
0514     }
0515   }
0516   fileout.close();
0517   std::cout << "Writes " << good << " events in the file " << outFile << " from " << nentries
0518             << " entries in run range " << runLow << ":" << runHigh << std::endl;
0519 }
0520 
0521 class CalibSortEvent {
0522 public:
0523   TChain *fChain;  //!pointer to the analyzed TTree/TChain
0524   Int_t fCurrent;  //!current Tree number in a TChain
0525 
0526   // Declaration of leaf types
0527   UInt_t t_RunNo;
0528   UInt_t t_EventNo;
0529   Int_t t_Tracks;
0530   Int_t t_TracksProp;
0531   Int_t t_TracksSaved;
0532   Int_t t_TracksLoose;
0533   Int_t t_TracksTight;
0534   Int_t t_allvertex;
0535   Bool_t t_TrigPass;
0536   Bool_t t_TrigPassSel;
0537   Bool_t t_L1Bit;
0538   std::vector<Bool_t> *t_hltbits;
0539   std::vector<int> *t_ietaAll;
0540   std::vector<int> *t_ietaGood;
0541   std::vector<int> *t_trackType;
0542 
0543   // List of branches
0544   TBranch *b_t_RunNo;        //!
0545   TBranch *b_t_EventNo;      //!
0546   TBranch *b_t_Tracks;       //!
0547   TBranch *b_t_TracksProp;   //!
0548   TBranch *b_t_TracksSaved;  //!
0549   TBranch *b_t_TracksLoose;  //!
0550   TBranch *b_t_TracksTight;  //!
0551   TBranch *b_t_allvertex;    //!
0552   TBranch *b_t_TrigPass;     //!
0553   TBranch *b_t_TrigPassSel;  //!
0554   TBranch *b_t_L1Bit;        //!
0555   TBranch *b_t_hltbits;      //!
0556   TBranch *b_t_ietaAll;      //!
0557   TBranch *b_t_ietaGood;     //!
0558   TBranch *b_t_trackType;    //!
0559 
0560   CalibSortEvent(const char *fname, std::string dirname, std::string prefix = "", bool append = false);
0561   virtual ~CalibSortEvent();
0562   virtual Int_t Cut(Long64_t entry);
0563   virtual Int_t GetEntry(Long64_t entry);
0564   virtual Long64_t LoadTree(Long64_t entry);
0565   virtual void Init(TChain *tree);
0566   virtual void Loop();
0567   virtual Bool_t Notify();
0568   virtual void Show(Long64_t entry = -1);
0569 
0570 private:
0571   std::string fname_, dirnm_, prefix_;
0572   bool append_;
0573 };
0574 
0575 CalibSortEvent::CalibSortEvent(const char *fname, std::string dirnm, std::string prefix, bool append)
0576     : fname_(fname), dirnm_(dirnm), prefix_(prefix), append_(append) {
0577   char treeName[400];
0578   sprintf(treeName, "%s/EventInfo", dirnm.c_str());
0579   TChain *chain = new TChain(treeName);
0580   std::cout << "Create a chain for " << treeName << " from " << fname << std::endl;
0581   if (!fillChain(chain, fname)) {
0582     std::cout << "*****No valid tree chain can be obtained*****" << std::endl;
0583   } else {
0584     std::cout << "Proceed with a tree chain with " << chain->GetEntries() << " entries" << std::endl;
0585     Init(chain);
0586   }
0587 }
0588 
0589 CalibSortEvent::~CalibSortEvent() {
0590   if (!fChain)
0591     return;
0592   delete fChain->GetCurrentFile();
0593 }
0594 
0595 Int_t CalibSortEvent::GetEntry(Long64_t entry) {
0596   // Read contents of entry.
0597   if (!fChain)
0598     return 0;
0599   return fChain->GetEntry(entry);
0600 }
0601 
0602 Long64_t CalibSortEvent::LoadTree(Long64_t entry) {
0603   // Set the environment to read one entry
0604   if (!fChain)
0605     return -5;
0606   Long64_t centry = fChain->LoadTree(entry);
0607   if (centry < 0)
0608     return centry;
0609   if (!fChain->InheritsFrom(TChain::Class()))
0610     return centry;
0611   TChain *chain = (TChain *)fChain;
0612   if (chain->GetTreeNumber() != fCurrent) {
0613     fCurrent = chain->GetTreeNumber();
0614     Notify();
0615   }
0616   return centry;
0617 }
0618 
0619 void CalibSortEvent::Init(TChain *tree) {
0620   // The Init() function is called when the selector needs to initialize
0621   // a new tree or chain. Typically here the branch addresses and branch
0622   // pointers of the tree will be set.
0623   // It is normally not necessary to make changes to the generated
0624   // code, but the routine can be extended by the user if needed.
0625   // Init() will be called many times when running on PROOF
0626   // (once per file to be processed).
0627 
0628   // Set branch addresses and branch pointers
0629   // Set object pointer
0630   t_hltbits = 0;
0631   t_ietaAll = 0;
0632   t_ietaGood = 0;
0633   t_trackType = 0;
0634   t_L1Bit = false;
0635   fChain = tree;
0636   fCurrent = -1;
0637   if (!tree)
0638     return;
0639   fChain->SetMakeClass(1);
0640   fChain->SetBranchAddress("t_RunNo", &t_RunNo, &b_t_RunNo);
0641   fChain->SetBranchAddress("t_EventNo", &t_EventNo, &b_t_EventNo);
0642   fChain->SetBranchAddress("t_Tracks", &t_Tracks, &b_t_Tracks);
0643   fChain->SetBranchAddress("t_TracksProp", &t_TracksProp, &b_t_TracksProp);
0644   fChain->SetBranchAddress("t_TracksSaved", &t_TracksSaved, &b_t_TracksSaved);
0645   fChain->SetBranchAddress("t_TracksLoose", &t_TracksLoose, &b_t_TracksLoose);
0646   fChain->SetBranchAddress("t_TracksTight", &t_TracksTight, &b_t_TracksTight);
0647   fChain->SetBranchAddress("t_allvertex", &t_allvertex, &b_t_allvertex);
0648   fChain->SetBranchAddress("t_TrigPass", &t_TrigPass, &b_t_TrigPass);
0649   fChain->SetBranchAddress("t_TrigPassSel", &t_TrigPassSel, &b_t_TrigPassSel);
0650   fChain->SetBranchAddress("t_L1Bit", &t_L1Bit, &b_t_L1Bit);
0651   fChain->SetBranchAddress("t_hltbits", &t_hltbits, &b_t_hltbits);
0652   fChain->SetBranchAddress("t_ietaAll", &t_ietaAll, &b_t_ietaAll);
0653   fChain->SetBranchAddress("t_ietaGood", &t_ietaGood, &b_t_ietaGood);
0654   fChain->SetBranchAddress("t_trackType", &t_trackType, &b_t_trackType);
0655   Notify();
0656 }
0657 
0658 Bool_t CalibSortEvent::Notify() {
0659   // The Notify() function is called when a new file is opened. This
0660   // can be either for a new TTree in a TChain or when when a new TTree
0661   // is started when using PROOF. It is normally not necessary to make changes
0662   // to the generated code, but the routine can be extended by the
0663   // user if needed. The return value is currently not used.
0664 
0665   return kTRUE;
0666 }
0667 
0668 void CalibSortEvent::Show(Long64_t entry) {
0669   // Print contents of entry.
0670   // If entry is not specified, print current entry
0671   if (!fChain)
0672     return;
0673   fChain->Show(entry);
0674 }
0675 
0676 Int_t CalibSortEvent::Cut(Long64_t) {
0677   // This function may be called from Loop.
0678   // returns  1 if entry is accepted.
0679   // returns -1 otherwise.
0680   return 1;
0681 }
0682 
0683 void CalibSortEvent::Loop() {
0684   //   In a ROOT session, you can do:
0685   //      Root > .L CalibMonitor.C+g
0686   //      Root > CalibSortEvent t
0687   //      Root > t.GetEntry(12); // Fill t data members with entry number 12
0688   //      Root > t.Show();       // Show values of entry 12
0689   //      Root > t.Show(16);     // Read and show values of entry 16
0690   //      Root > t.Loop();       // Loop on all entries
0691   //
0692 
0693   //     This is the loop skeleton where:
0694   //    jentry is the global entry number in the chain
0695   //    ientry is the entry number in the current Tree
0696   //  Note that the argument to GetEntry must be:
0697   //    jentry for TChain::GetEntry
0698   //    ientry for TTree::GetEntry and TBranch::GetEntry
0699   //
0700   //       To read only selected branches, Insert statements like:
0701   // METHOD1:
0702   //    fChain->SetBranchStatus("*",0);  // disable all branches
0703   //    fChain->SetBranchStatus("branchname",1);  // activate branchname
0704   // METHOD2: replace line
0705   //    fChain->GetEntry(jentry);       //read all branches
0706   //by  b_branchname->GetEntry(ientry); //read only this branch
0707   if (fChain == 0)
0708     return;
0709 
0710   std::ofstream fileout;
0711   if (!append_) {
0712     fileout.open("runevents.txt", std::ofstream::out);
0713     std::cout << "Opens runevents.txt in output mode" << std::endl;
0714   } else {
0715     fileout.open("runevents.txt", std::ofstream::app);
0716     std::cout << "Opens runevents.txt in append mode" << std::endl;
0717   }
0718   fileout << "Input file: " << fname_ << " Directory: " << dirnm_ << " Prefix: " << prefix_ << std::endl;
0719   UInt_t runLow(99999999), runHigh(0);
0720   Long64_t nbytes(0), nb(0), good(0);
0721   Long64_t nentries = fChain->GetEntriesFast();
0722   for (Long64_t jentry = 0; jentry < nentries; jentry++) {
0723     Long64_t ientry = LoadTree(jentry);
0724     if (ientry < 0)
0725       break;
0726     nb = fChain->GetEntry(jentry);
0727     nbytes += nb;
0728     if (t_RunNo > 200000 && t_RunNo < 800000) {
0729       if (t_RunNo < runLow)
0730         runLow = t_RunNo;
0731       if (t_RunNo > runHigh)
0732         runHigh = t_RunNo;
0733     }
0734     good++;
0735     fileout << good << " " << jentry << " " << t_RunNo << " " << t_EventNo << std::endl;
0736   }
0737   fileout.close();
0738   std::cout << "Writes " << good << " events in the file events.txt from " << nentries << " entries in run range "
0739             << runLow << ":" << runHigh << std::endl;
0740 }
0741 
0742 void readRecords(std::string fname, std::vector<record> &records, bool debug) {
0743   records.clear();
0744   ifstream infile(fname.c_str());
0745   if (!infile.is_open()) {
0746     std::cout << "Cannot open " << fname << std::endl;
0747   } else {
0748     while (1) {
0749       int ser, ent, r, ev, ie;
0750       double p;
0751       infile >> ser >> ent >> r >> ev >> ie >> p;
0752       if (!infile.good())
0753         break;
0754       record rec(ser, ent, r, ev, ie, p);
0755       records.push_back(rec);
0756     }
0757     infile.close();
0758   }
0759   std::cout << "Reads " << records.size() << " records from " << fname << std::endl;
0760   if (debug) {
0761     for (unsigned int k = 0; k < records.size(); ++k) {
0762       if (k % 100 == 0)
0763         std::cout << "[" << records[k].serial_ << ":" << records[k].entry_ << "] " << records[k].run_ << ":"
0764                   << records[k].event_ << " " << records[k].ieta_ << " " << records[k].p_ << std::endl;
0765     }
0766   }
0767 }
0768 
0769 void readMap(std::string fname, std::map<std::pair<int, int>, int> &records, bool debug) {
0770   records.clear();
0771   ifstream infile(fname.c_str());
0772   if (!infile.is_open()) {
0773     std::cout << "Cannot open " << fname << std::endl;
0774   } else {
0775     while (1) {
0776       int ser, ent, r, ev, ie;
0777       double p;
0778       infile >> ser >> ent >> r >> ev >> ie >> p;
0779       if (!infile.good())
0780         break;
0781       std::pair<int, int> key(r, ev);
0782       if (records.find(key) == records.end())
0783         records[key] = ent;
0784     }
0785     infile.close();
0786   }
0787   std::cout << "Reads " << records.size() << " records from " << fname << std::endl;
0788   if (debug) {
0789     unsigned k(0);
0790     for (std::map<std::pair<int, int>, int>::iterator itr = records.begin(); itr != records.end(); ++itr, ++k) {
0791       if (k % 100 == 0)
0792         std::cout << "[" << k << "] " << itr->second << ":" << (itr->first).first << ":" << (itr->first).second << "\n";
0793     }
0794   }
0795 }
0796 
0797 void sort(std::vector<record> &records, bool debug) {
0798   // Use std::sort
0799   std::sort(records.begin(), records.end(), recordLess());
0800   if (debug) {
0801     for (unsigned int k = 0; k < records.size(); ++k) {
0802       std::cout << "[" << k << ":" << records[k].serial_ << ":" << records[k].entry_ << "] " << records[k].run_ << ":"
0803                 << records[k].event_ << " " << records[k].ieta_ << " " << records[k].p_ << std::endl;
0804     }
0805   }
0806 }
0807 
0808 void duplicate(std::string fname, std::vector<record> &records, bool debug) {
0809   std::ofstream file;
0810   file.open(fname.c_str(), std::ofstream::out);
0811   std::cout << "List of entry names of duplicate events" << std::endl;
0812   int duplicate(0), dupl40(0);
0813   for (unsigned int k = 1; k < records.size(); ++k) {
0814     if ((records[k].run_ == records[k - 1].run_) && (records[k].event_ == records[k - 1].event_) &&
0815         (records[k].ieta_ == records[k - 1].ieta_) && (fabs(records[k].p_ - records[k - 1].p_) < 0.0001)) {
0816       // This is a duplicate event - reject the one with larger serial #
0817       if (records[k].entry_ < records[k - 1].entry_) {
0818         record swap = records[k - 1];
0819         records[k - 1] = records[k];
0820         records[k] = swap;
0821       }
0822       if (debug) {
0823         std::cout << "Serial " << records[k - 1].serial_ << ":" << records[k].serial_ << " Entry "
0824                   << records[k - 1].entry_ << ":" << records[k].entry_ << " Run " << records[k - 1].run_ << ":"
0825                   << records[k].run_ << " Event " << records[k - 1].event_ << " " << records[k].event_ << " Eta "
0826                   << records[k - 1].ieta_ << " " << records[k].ieta_ << " p " << records[k - 1].p_ << ":"
0827                   << records[k].p_ << std::endl;
0828       }
0829       file << records[k].entry_ << std::endl;
0830       duplicate++;
0831       if (records[k].p_ >= 40.0 && records[k].p_ <= 60.0)
0832         dupl40++;
0833     }
0834   }
0835   file.close();
0836   std::cout << "Total # of duplcate events " << duplicate << " (" << dupl40 << " with p 40:60)" << std::endl;
0837 }
0838 
0839 void findDuplicate(std::string infile, std::string outfile, bool debug = false) {
0840   std::vector<record> records;
0841   readRecords(infile, records, debug);
0842   sort(records, debug);
0843   duplicate(outfile, records, debug);
0844 }
0845 
0846 void findCommon(std::string infile1, std::string infile2, std::string infile3, std::string outfile, bool debug = false) {
0847   std::map<std::pair<int, int>, int> map1, map2, map3;
0848   readMap(infile1, map1, debug);
0849   readMap(infile2, map2, debug);
0850   readMap(infile3, map3, debug);
0851   bool check3 = (map3.size() > 0);
0852   std::ofstream file;
0853   file.open(outfile.c_str(), std::ofstream::out);
0854   unsigned int k(0), good(0);
0855   for (std::map<std::pair<int, int>, int>::iterator itr = map1.begin(); itr != map1.end(); ++itr, ++k) {
0856     std::pair<int, int> key = itr->first;
0857     bool ok = (map2.find(key) != map2.end());
0858     if (ok && check3)
0859       ok = (map3.find(key) != map3.end());
0860     if (debug && k % 100 == 0)
0861       std::cout << "[" << k << "] Run " << key.first << " Event " << key.second << " Flag " << ok << std::endl;
0862     if (ok) {
0863       ++good;
0864       file << key.first << "   " << key.second << std::endl;
0865     }
0866   }
0867   file.close();
0868   std::cout << "Total # of common events " << good << " written to o/p file " << outfile << std::endl;
0869 }
0870 
0871 void readRecordEvents(std::string fname, std::vector<recordEvent> &records, bool debug) {
0872   records.clear();
0873   ifstream infile(fname.c_str());
0874   if (!infile.is_open()) {
0875     std::cout << "Cannot open " << fname << std::endl;
0876   } else {
0877     while (1) {
0878       unsigned int ser, ent, r, ev;
0879       infile >> ser >> ent >> r >> ev;
0880       if (!infile.good())
0881         break;
0882       recordEvent rec(ser, ent, r, ev);
0883       records.push_back(rec);
0884     }
0885     infile.close();
0886   }
0887   std::cout << "Reads " << records.size() << " records from " << fname << std::endl;
0888   if (debug) {
0889     for (unsigned int k = 0; k < records.size(); ++k) {
0890       if (k % 100 == 0)
0891         std::cout << "[" << records[k].serial_ << ":" << records[k].entry_ << "] " << records[k].run_ << ":"
0892                   << records[k].event_ << std::endl;
0893     }
0894   }
0895 }
0896 
0897 void sortEvent(std::vector<recordEvent> &records, bool debug) {
0898   // Use std::sort
0899   std::sort(records.begin(), records.end(), recordEventLess());
0900   if (debug) {
0901     for (unsigned int k = 0; k < records.size(); ++k) {
0902       std::cout << "[" << k << ":" << records[k].serial_ << ":" << records[k].entry_ << "] " << records[k].run_ << ":"
0903                 << records[k].event_ << std::endl;
0904     }
0905   }
0906 }
0907 
0908 void duplicateEvent(std::string fname, std::vector<recordEvent> &records, bool debug) {
0909   std::ofstream file;
0910   file.open(fname.c_str(), std::ofstream::out);
0911   std::cout << "List of entry names of duplicate events" << std::endl;
0912   int duplicate(0);
0913   for (unsigned int k = 1; k < records.size(); ++k) {
0914     if ((records[k].run_ == records[k - 1].run_) && (records[k].event_ == records[k - 1].event_)) {
0915       // This is a duplicate event - reject the one with larger serial #
0916       if (records[k].entry_ < records[k - 1].entry_) {
0917         recordEvent swap = records[k - 1];
0918         records[k - 1] = records[k];
0919         records[k] = swap;
0920       }
0921       if (debug) {
0922         std::cout << "Serial " << records[k - 1].serial_ << ":" << records[k].serial_ << " Entry "
0923                   << records[k - 1].entry_ << ":" << records[k].entry_ << " Run " << records[k - 1].run_ << ":"
0924                   << records[k].run_ << " Event " << records[k - 1].event_ << " " << records[k].event_ << std::endl;
0925       }
0926       file << records[k].entry_ << std::endl;
0927       duplicate++;
0928     }
0929   }
0930   file.close();
0931   std::cout << "Total # of duplcate events " << duplicate << std::endl;
0932 }
0933 
0934 void findDuplicateEvent(std::string infile, std::string outfile, bool debug = false) {
0935   std::vector<recordEvent> records;
0936   readRecordEvents(infile, records, debug);
0937   sortEvent(records, debug);
0938   duplicateEvent(outfile, records, debug);
0939 }
0940 
0941 void combine(const std::string fname1,
0942              const std::string fname2,
0943              const std::string fname,
0944              const std::string outfile = "",
0945              int debug = 0) {
0946   ifstream infile1(fname1.c_str());
0947   ifstream infile2(fname2.c_str());
0948   if ((!infile1.is_open()) || (!infile2.is_open())) {
0949     std::cout << "Cannot open " << fname1 << " or " << fname2 << std::endl;
0950   } else {
0951     std::ofstream file;
0952     file.open(fname.c_str(), std::ofstream::out);
0953     TH1D *hist0 = new TH1D("W0", "Correction Factor (All)", 100, 0.0, 10.0);
0954     TH1D *hist1 = new TH1D("W1", "Correction Factor (Barrel)", 100, 0.0, 10.0);
0955     TH1D *hist2 = new TH1D("W2", "Correction Factor (Endcap)", 100, 0.0, 10.0);
0956     TProfile *prof = new TProfile("P", "Correction vs i#eta", 60, -30, 30, 0, 100, "");
0957     unsigned int kount(0), kout(0);
0958     double wmaxb(0), wmaxe(0);
0959     while (1) {
0960       Long64_t entry1, entry2;
0961       int ieta1, ieta2;
0962       double wt1, wt2;
0963       infile1 >> entry1 >> ieta1 >> wt1;
0964       infile2 >> entry2 >> ieta2 >> wt2;
0965       if ((!infile1.good()) || (!infile2.good()))
0966         break;
0967       ++kount;
0968       if (debug > 1) {
0969         std::cout << kount << " Enrty No. " << entry1 << ":" << entry2 << " eta " << ieta1 << ":" << ieta2 << " wt "
0970                   << wt1 << ":" << wt2 << std::endl;
0971       }
0972       if (entry1 == entry2) {
0973         int ieta = fabs(ieta1);
0974         double wt = (ieta >= 16) ? wt1 : wt2;
0975         file << std::setw(8) << entry1 << " " << std::setw(12) << std::setprecision(8) << wt << std::endl;
0976         hist0->Fill(wt);
0977         if (ieta >= 16) {
0978           hist2->Fill(wt);
0979           if (wt > wmaxe)
0980             wmaxe = wt;
0981         } else {
0982           hist1->Fill(wt);
0983           if (wt > wmaxb)
0984             wmaxb = wt;
0985         }
0986         prof->Fill(ieta1, wt);
0987         ++kout;
0988       } else if (debug > 0) {
0989         std::cout << kount << " Entry " << entry1 << ":" << entry2 << " eta " << ieta1 << ":" << ieta2 << " wt " << wt1
0990                   << ":" << wt2 << " mismatch" << std::endl;
0991       }
0992     }
0993     infile1.close();
0994     infile2.close();
0995     file.close();
0996     std::cout << "Writes " << kout << " entries to " << fname << " from " << kount << " events from " << fname1
0997               << " and " << fname2 << " maximum correction factor " << wmaxb << " (Barrel) " << wmaxe << " (Endcap)"
0998               << std::endl;
0999     if (outfile != "") {
1000       TFile *theFile = new TFile(outfile.c_str(), "RECREATE");
1001       theFile->cd();
1002       hist0->Write();
1003       hist1->Write();
1004       hist2->Write();
1005       prof->Write();
1006       theFile->Close();
1007     }
1008   }
1009 }
1010 
1011 void plotCombine(const char *infile, int flag = -1, bool drawStatBox = true, bool save = true) {
1012   int flag1 = (flag >= 0 ? (flag / 10) % 10 : 0);
1013   int flag2 = (flag >= 0 ? (flag % 10) : 0);
1014   std::string name[4] = {"W0", "W1", "W2", "P"};
1015   std::string xtitl[4] = {
1016       "Correction Factor (All)", "Correction Factor (Barrel)", "Correction Factor (Endcap)", "i#eta"};
1017   std::string ytitl[4] = {"Track", "Track", "Track", "Correction Factor"};
1018   char title[100];
1019   std::string mom[2] = {"all momentum", "p = 40:60 GeV"};
1020   std::string iso[2] = {"loose", "tight"};
1021   if (flag < 0) {
1022     sprintf(title, "All tracks");
1023   } else {
1024     sprintf(title, "Good tracks of %s with %s isolation", mom[flag1].c_str(), iso[flag2].c_str());
1025   }
1026 
1027   gStyle->SetCanvasBorderMode(0);
1028   gStyle->SetCanvasColor(kWhite);
1029   gStyle->SetPadColor(kWhite);
1030   gStyle->SetFillColor(kWhite);
1031   gStyle->SetOptTitle(0);
1032   gStyle->SetOptStat(111110);
1033   gStyle->SetOptFit(0);
1034   TFile *file = new TFile(infile);
1035   for (int k = 0; k < 4; ++k) {
1036     char nameh[40], namep[50];
1037     if (flag >= 0)
1038       sprintf(nameh, "%s%d%d", name[k].c_str(), flag1, flag2);
1039     else
1040       sprintf(nameh, "%s", name[k].c_str());
1041     sprintf(namep, "c_%s", nameh);
1042     TCanvas *pad(nullptr);
1043     if (k == 3) {
1044       TProfile *hist = (TProfile *)file->FindObjectAny(nameh);
1045       if (hist != nullptr) {
1046         pad = new TCanvas(namep, namep, 700, 500);
1047         pad->SetRightMargin(0.10);
1048         pad->SetTopMargin(0.10);
1049         hist->GetXaxis()->SetTitleSize(0.04);
1050         hist->GetXaxis()->SetTitle(xtitl[k].c_str());
1051         hist->GetYaxis()->SetTitle(ytitl[k].c_str());
1052         hist->GetYaxis()->SetLabelOffset(0.005);
1053         hist->GetYaxis()->SetTitleSize(0.04);
1054         hist->GetYaxis()->SetLabelSize(0.035);
1055         hist->GetYaxis()->SetTitleOffset(1.10);
1056         hist->SetMarkerStyle(20);
1057         hist->Draw();
1058         pad->Update();
1059         TPaveStats *st1 = (TPaveStats *)hist->GetListOfFunctions()->FindObject("stats");
1060         if (st1 != nullptr && drawStatBox) {
1061           st1->SetY1NDC(0.70);
1062           st1->SetY2NDC(0.90);
1063           st1->SetX1NDC(0.65);
1064           st1->SetX2NDC(0.90);
1065         }
1066       }
1067     } else {
1068       TH1D *hist = (TH1D *)file->FindObjectAny(nameh);
1069       if (hist != nullptr) {
1070         pad = new TCanvas(namep, namep, 700, 500);
1071         pad->SetRightMargin(0.10);
1072         pad->SetTopMargin(0.10);
1073         pad->SetLogy();
1074         hist->GetXaxis()->SetTitleSize(0.04);
1075         hist->GetXaxis()->SetTitle(xtitl[k].c_str());
1076         hist->GetYaxis()->SetTitle(ytitl[k].c_str());
1077         hist->GetYaxis()->SetLabelOffset(0.005);
1078         hist->GetYaxis()->SetTitleSize(0.04);
1079         hist->GetYaxis()->SetLabelSize(0.035);
1080         hist->GetYaxis()->SetTitleOffset(1.10);
1081         hist->SetMarkerStyle(20);
1082         hist->Draw();
1083         pad->Update();
1084         TPaveStats *st1 = (TPaveStats *)hist->GetListOfFunctions()->FindObject("stats");
1085         if (st1 != nullptr && drawStatBox) {
1086           st1->SetY1NDC(0.70);
1087           st1->SetY2NDC(0.90);
1088           st1->SetX1NDC(0.65);
1089           st1->SetX2NDC(0.90);
1090         }
1091       }
1092       TPaveText *txt0 = new TPaveText(0.10, 0.91, 0.90, 0.96, "blNDC");
1093       txt0->SetFillColor(0);
1094       txt0->AddText(title);
1095       txt0->Draw("same");
1096       pad->Update();
1097     }
1098     if (pad != nullptr) {
1099       pad->Modified();
1100       pad->Update();
1101       if (save) {
1102         sprintf(namep, "%s.pdf", pad->GetName());
1103         pad->Print(namep);
1104       }
1105     }
1106   }
1107 }
1108 
1109 class CalibPlotCombine {
1110 public:
1111   TChain *fChain;  //!pointer to the analyzed TTree or TChain
1112   Int_t fCurrent;  //!current Tree number in a TChain
1113 
1114   // Declaration of leaf types
1115   Int_t t_Run;
1116   Int_t t_Event;
1117   Int_t t_DataType;
1118   Int_t t_ieta;
1119   Int_t t_iphi;
1120   Double_t t_EventWeight;
1121   Int_t t_nVtx;
1122   Int_t t_nTrk;
1123   Int_t t_goodPV;
1124   Double_t t_l1pt;
1125   Double_t t_l1eta;
1126   Double_t t_l1phi;
1127   Double_t t_l3pt;
1128   Double_t t_l3eta;
1129   Double_t t_l3phi;
1130   Double_t t_p;
1131   Double_t t_pt;
1132   Double_t t_phi;
1133   Double_t t_mindR1;
1134   Double_t t_mindR2;
1135   Double_t t_eMipDR;
1136   Double_t t_eHcal;
1137   Double_t t_eHcal10;
1138   Double_t t_eHcal30;
1139   Double_t t_hmaxNearP;
1140   Double_t t_rhoh;
1141   Bool_t t_selectTk;
1142   Bool_t t_qltyFlag;
1143   Bool_t t_qltyMissFlag;
1144   Bool_t t_qltyPVFlag;
1145   Double_t t_gentrackP;
1146   std::vector<unsigned int> *t_DetIds;
1147   std::vector<double> *t_HitEnergies;
1148   std::vector<bool> *t_trgbits;
1149   std::vector<unsigned int> *t_DetIds1;
1150   std::vector<unsigned int> *t_DetIds3;
1151   std::vector<double> *t_HitEnergies1;
1152   std::vector<double> *t_HitEnergies3;
1153 
1154   // List of branches
1155   TBranch *b_t_Run;           //!
1156   TBranch *b_t_Event;         //!
1157   TBranch *b_t_DataType;      //!
1158   TBranch *b_t_ieta;          //!
1159   TBranch *b_t_iphi;          //!
1160   TBranch *b_t_EventWeight;   //!
1161   TBranch *b_t_nVtx;          //!
1162   TBranch *b_t_nTrk;          //!
1163   TBranch *b_t_goodPV;        //!
1164   TBranch *b_t_l1pt;          //!
1165   TBranch *b_t_l1eta;         //!
1166   TBranch *b_t_l1phi;         //!
1167   TBranch *b_t_l3pt;          //!
1168   TBranch *b_t_l3eta;         //!
1169   TBranch *b_t_l3phi;         //!
1170   TBranch *b_t_p;             //!
1171   TBranch *b_t_pt;            //!
1172   TBranch *b_t_phi;           //!
1173   TBranch *b_t_mindR1;        //!
1174   TBranch *b_t_mindR2;        //!
1175   TBranch *b_t_eMipDR;        //!
1176   TBranch *b_t_eHcal;         //!
1177   TBranch *b_t_eHcal10;       //!
1178   TBranch *b_t_eHcal30;       //!
1179   TBranch *b_t_hmaxNearP;     //!
1180   TBranch *b_t_rhoh;          //!
1181   TBranch *b_t_selectTk;      //!
1182   TBranch *b_t_qltyFlag;      //!
1183   TBranch *b_t_qltyMissFlag;  //!
1184   TBranch *b_t_qltyPVFlag;    //!
1185   TBranch *b_t_gentrackP;     //!
1186   TBranch *b_t_DetIds;        //!
1187   TBranch *b_t_HitEnergies;   //!
1188   TBranch *b_t_trgbits;       //!
1189   TBranch *b_t_DetIds1;       //!
1190   TBranch *b_t_DetIds3;       //!
1191   TBranch *b_t_HitEnergies1;  //!
1192   TBranch *b_t_HitEnergies3;  //!
1193 
1194   CalibPlotCombine(const char *fname,
1195                    const char *corrFileName,
1196                    const std::string dirname = "HcalIsoTrkAnalyzer",
1197                    int flag = 10);
1198   virtual ~CalibPlotCombine();
1199   virtual Int_t Cut(Long64_t entry);
1200   virtual Int_t GetEntry(Long64_t entry);
1201   virtual Long64_t LoadTree(Long64_t entry);
1202   virtual void Init(TChain *);
1203   virtual void Loop();
1204   virtual Bool_t Notify();
1205   virtual void Show(Long64_t entry = -1);
1206   bool goodTrack(double eHcal, double cut);
1207   void savePlot(const std::string &theName, bool append = false);
1208 
1209 private:
1210   CalibCorr *cFactor_;
1211   int flag_, ifDepth_;
1212   bool selectP_, tightSelect_;
1213   TH1D *hist0_, *hist1_, *hist2_;
1214   TProfile *prof_;
1215 };
1216 
1217 CalibPlotCombine::CalibPlotCombine(const char *fname, const char *corrFileName, const std::string dirnm, int flag)
1218     : cFactor_(nullptr), flag_(flag), ifDepth_(3) {
1219   selectP_ = (((flag_ / 10) % 10) > 0);
1220   tightSelect_ = ((flag_ % 10) > 0);
1221   flag_ = 0;
1222   if (selectP_)
1223     flag_ += 10;
1224   if (tightSelect_)
1225     ++flag_;
1226   std::cout << "Selection on momentum range: " << selectP_ << std::endl
1227             << "Tight isolation flag:        " << tightSelect_ << std::endl
1228             << "Flag:                        " << flag_ << std::endl;
1229   char treeName[400];
1230   sprintf(treeName, "%s/CalibTree", dirnm.c_str());
1231   TChain *chain = new TChain(treeName);
1232   std::cout << "Create a chain for " << treeName << " from " << fname << std::endl;
1233   if (!fillChain(chain, fname)) {
1234     std::cout << "*****No valid tree chain can be obtained*****" << std::endl;
1235   } else {
1236     std::cout << "Proceed with a tree chain with " << chain->GetEntries() << " entries" << std::endl;
1237     Init(chain);
1238   }
1239   cFactor_ = new CalibCorr(corrFileName, ifDepth_, false);
1240 }
1241 
1242 CalibPlotCombine::~CalibPlotCombine() {
1243   delete cFactor_;
1244   if (!fChain)
1245     return;
1246   delete fChain->GetCurrentFile();
1247 }
1248 
1249 Int_t CalibPlotCombine::GetEntry(Long64_t entry) {
1250   // Read contents of entry.
1251   if (!fChain)
1252     return 0;
1253   return fChain->GetEntry(entry);
1254 }
1255 
1256 Long64_t CalibPlotCombine::LoadTree(Long64_t entry) {
1257   // Set the environment to read one entry
1258   if (!fChain)
1259     return -5;
1260   Long64_t centry = fChain->LoadTree(entry);
1261   if (centry < 0)
1262     return centry;
1263   if (!fChain->InheritsFrom(TChain::Class()))
1264     return centry;
1265   TChain *chain = (TChain *)fChain;
1266   if (chain->GetTreeNumber() != fCurrent) {
1267     fCurrent = chain->GetTreeNumber();
1268     Notify();
1269   }
1270   return centry;
1271 }
1272 
1273 void CalibPlotCombine::Init(TChain *tree) {
1274   // The Init() function is called when the selector needs to initialize
1275   // a new tree or chain. Typically here the branch addresses and branch
1276   // pointers of the tree will be set.
1277   // It is normally not necessary to make changes to the generated
1278   // code, but the routine can be extended by the user if needed.
1279   // Init() will be called many times when running on PROOF
1280   // (once per file to be processed).
1281 
1282   // Set object pointer
1283   t_DetIds = 0;
1284   t_DetIds1 = 0;
1285   t_DetIds3 = 0;
1286   t_HitEnergies = 0;
1287   t_HitEnergies1 = 0;
1288   t_HitEnergies3 = 0;
1289   t_trgbits = 0;
1290   // Set branch addresses and branch pointers
1291   fChain = tree;
1292   fCurrent = -1;
1293   if (!tree)
1294     return;
1295   fChain->SetMakeClass(1);
1296 
1297   fChain->SetBranchAddress("t_Run", &t_Run, &b_t_Run);
1298   fChain->SetBranchAddress("t_Event", &t_Event, &b_t_Event);
1299   fChain->SetBranchAddress("t_DataType", &t_DataType, &b_t_DataType);
1300   fChain->SetBranchAddress("t_ieta", &t_ieta, &b_t_ieta);
1301   fChain->SetBranchAddress("t_iphi", &t_iphi, &b_t_iphi);
1302   fChain->SetBranchAddress("t_EventWeight", &t_EventWeight, &b_t_EventWeight);
1303   fChain->SetBranchAddress("t_nVtx", &t_nVtx, &b_t_nVtx);
1304   fChain->SetBranchAddress("t_nTrk", &t_nTrk, &b_t_nTrk);
1305   fChain->SetBranchAddress("t_goodPV", &t_goodPV, &b_t_goodPV);
1306   fChain->SetBranchAddress("t_l1pt", &t_l1pt, &b_t_l1pt);
1307   fChain->SetBranchAddress("t_l1eta", &t_l1eta, &b_t_l1eta);
1308   fChain->SetBranchAddress("t_l1phi", &t_l1phi, &b_t_l1phi);
1309   fChain->SetBranchAddress("t_l3pt", &t_l3pt, &b_t_l3pt);
1310   fChain->SetBranchAddress("t_l3eta", &t_l3eta, &b_t_l3eta);
1311   fChain->SetBranchAddress("t_l3phi", &t_l3phi, &b_t_l3phi);
1312   fChain->SetBranchAddress("t_p", &t_p, &b_t_p);
1313   fChain->SetBranchAddress("t_pt", &t_pt, &b_t_pt);
1314   fChain->SetBranchAddress("t_phi", &t_phi, &b_t_phi);
1315   fChain->SetBranchAddress("t_mindR1", &t_mindR1, &b_t_mindR1);
1316   fChain->SetBranchAddress("t_mindR2", &t_mindR2, &b_t_mindR2);
1317   fChain->SetBranchAddress("t_eMipDR", &t_eMipDR, &b_t_eMipDR);
1318   fChain->SetBranchAddress("t_eHcal", &t_eHcal, &b_t_eHcal);
1319   fChain->SetBranchAddress("t_eHcal10", &t_eHcal10, &b_t_eHcal10);
1320   fChain->SetBranchAddress("t_eHcal30", &t_eHcal30, &b_t_eHcal30);
1321   fChain->SetBranchAddress("t_hmaxNearP", &t_hmaxNearP, &b_t_hmaxNearP);
1322   fChain->SetBranchAddress("t_rhoh", &t_rhoh, &b_t_rhoh);
1323   fChain->SetBranchAddress("t_selectTk", &t_selectTk, &b_t_selectTk);
1324   fChain->SetBranchAddress("t_qltyFlag", &t_qltyFlag, &b_t_qltyFlag);
1325   fChain->SetBranchAddress("t_qltyMissFlag", &t_qltyMissFlag, &b_t_qltyMissFlag);
1326   fChain->SetBranchAddress("t_qltyPVFlag", &t_qltyPVFlag, &b_t_qltyPVFlag);
1327   fChain->SetBranchAddress("t_gentrackP", &t_gentrackP, &b_t_gentrackP);
1328   fChain->SetBranchAddress("t_DetIds", &t_DetIds, &b_t_DetIds);
1329   fChain->SetBranchAddress("t_HitEnergies", &t_HitEnergies, &b_t_HitEnergies);
1330   fChain->SetBranchAddress("t_trgbits", &t_trgbits, &b_t_trgbits);
1331   fChain->SetBranchAddress("t_DetIds1", &t_DetIds1, &b_t_DetIds1);
1332   fChain->SetBranchAddress("t_DetIds3", &t_DetIds3, &b_t_DetIds3);
1333   fChain->SetBranchAddress("t_HitEnergies1", &t_HitEnergies1, &b_t_HitEnergies1);
1334   fChain->SetBranchAddress("t_HitEnergies3", &t_HitEnergies3, &b_t_HitEnergies3);
1335   Notify();
1336 
1337   // Book the histograms
1338   char name[20], title[100];
1339   std::string mom[2] = {"all momentum", "p = 40:60 GeV"};
1340   std::string iso[2] = {"loose", "tight"};
1341   std::string rng[3] = {"All", "Barrel", "Endcap"};
1342   int flag1 = (flag_ / 10) % 10;
1343   int flag2 = (flag_ % 10);
1344   sprintf(name, "W0%d%d", flag1, flag2);
1345   sprintf(title,
1346           "Correction Factor (%s) for tracks with %s and %s isolation",
1347           rng[0].c_str(),
1348           mom[flag1].c_str(),
1349           iso[flag2].c_str());
1350   hist0_ = new TH1D(name, title, 100, 0.0, 10.0);
1351   sprintf(name, "W1%d%d", flag1, flag2);
1352   sprintf(title,
1353           "Correction Factor (%s) for tracks with %s and %s isolation",
1354           rng[1].c_str(),
1355           mom[flag1].c_str(),
1356           iso[flag2].c_str());
1357   hist1_ = new TH1D(name, title, 100, 0.0, 10.0);
1358   sprintf(name, "W2%d%d", flag1, flag2);
1359   sprintf(title,
1360           "Correction Factor (%s) for tracks with %s and %s isolation",
1361           rng[2].c_str(),
1362           mom[flag1].c_str(),
1363           iso[flag2].c_str());
1364   hist2_ = new TH1D(name, title, 100, 0.0, 10.0);
1365   sprintf(name, "P%d%d", flag1, flag2);
1366   sprintf(
1367       title, "Correction Factor vs i#eta for tracks with %s and %s isolation", mom[flag1].c_str(), iso[flag2].c_str());
1368   prof_ = new TProfile(name, title, 60, -30, 30, 0, 100, "");
1369 }
1370 
1371 Bool_t CalibPlotCombine::Notify() {
1372   // The Notify() function is called when a new file is opened. This
1373   // can be either for a new TTree in a TChain or when when a new TTree
1374   // is started when using PROOF. It is normally not necessary to make changes
1375   // to the generated code, but the routine can be extended by the
1376   // user if needed. The return value is currently not used.
1377 
1378   return kTRUE;
1379 }
1380 
1381 void CalibPlotCombine::Show(Long64_t entry) {
1382   // Print contents of entry.
1383   // If entry is not specified, print current entry
1384   if (!fChain)
1385     return;
1386   fChain->Show(entry);
1387 }
1388 
1389 Int_t CalibPlotCombine::Cut(Long64_t) {
1390   // This function may be called from Loop.
1391   // returns  1 if entry is accepted.
1392   // returns -1 otherwise.
1393   return 1;
1394 }
1395 
1396 void CalibPlotCombine::Loop() {
1397   //   In a ROOT session, you can do:
1398   //      Root > .L CalibPlotCombine.C
1399   //      Root > CalibPlotCombine t
1400   //      Root > t.GetEntry(12); // Fill t data members with entry number 12
1401   //      Root > t.Show();       // Show values of entry 12
1402   //      Root > t.Show(16);     // Read and show values of entry 16
1403   //      Root > t.Loop();       // Loop on all entries
1404   //
1405 
1406   //   This is the loop skeleton where:
1407   //      jentry is the global entry number in the chain
1408   //      ientry is the entry number in the current Tree
1409   //   Note that the argument to GetEntry must be:
1410   //      jentry for TChain::GetEntry
1411   //      ientry for TTree::GetEntry and TBranch::GetEntry
1412   //
1413   //       To read only selected branches, Insert statements like:
1414   // METHOD1:
1415   //    fChain->SetBranchStatus("*",0);  // disable all branches
1416   //    fChain->SetBranchStatus("branchname",1);  // activate branchname
1417   // METHOD2: replace line
1418   //    fChain->GetEntry(jentry);       //read all branches
1419   //by  b_branchname->GetEntry(ientry); //read only this branch
1420 
1421   if (fChain == 0)
1422     return;
1423   Long64_t nentries = fChain->GetEntriesFast();
1424   std::cout << "Total entries " << nentries << std::endl;
1425   Long64_t nbytes(0), nb(0);
1426   double wminb(99999999), wmaxb(0), wmine(99999999), wmaxe(0);
1427   for (Long64_t jentry = 0; jentry < nentries; jentry++) {
1428     Long64_t ientry = LoadTree(jentry);
1429     if (ientry < 0)
1430       break;
1431     nb = fChain->GetEntry(jentry);
1432     nbytes += nb;
1433     if (jentry % 1000000 == 0)
1434       std::cout << "Entry " << jentry << " Run " << t_Run << " Event " << t_Event << std::endl;
1435     double cut = (t_p > 20) ? (tightSelect_ ? 2.0 : 10.0) : 0.0;
1436     double eHcal(t_eHcal);
1437     if (goodTrack(eHcal, cut)) {
1438       if ((!selectP_) || ((t_p > 40.0) && (t_p < 60.0))) {
1439         int ieta = fabs(t_ieta);
1440         double wt = cFactor_->getTrueCorr(ientry);
1441         hist0_->Fill(wt);
1442         if (ieta >= 16) {
1443           hist2_->Fill(wt);
1444           wmaxe = std::max(wmaxe, wt);
1445           wmine = std::min(wmine, wt);
1446         } else {
1447           hist1_->Fill(wt);
1448           wmaxb = std::max(wmaxb, wt);
1449           wminb = std::min(wminb, wt);
1450         }
1451         prof_->Fill(t_ieta, wt);
1452       }
1453     }
1454   }
1455   std::cout << "Minimum and maximum correction factors " << wminb << ":" << wmaxb << " (Barrel) " << wmine << ":"
1456             << wmaxe << " (Endcap)" << std::endl;
1457 }
1458 
1459 bool CalibPlotCombine::goodTrack(double eHcal, double cut) {
1460   return ((t_qltyFlag) && t_selectTk && (t_hmaxNearP < cut) && (t_eMipDR < 1.0) && (eHcal > 0.001));
1461 }
1462 
1463 void CalibPlotCombine::savePlot(const std::string &theName, bool append) {
1464   TFile *theFile = (append ? (new TFile(theName.c_str(), "UPDATE")) : (new TFile(theName.c_str(), "RECREATE")));
1465   theFile->cd();
1466   TH1D *hist;
1467   hist = (TH1D *)(hist0_->Clone());
1468   hist->Write();
1469   hist = (TH1D *)(hist1_->Clone());
1470   hist->Write();
1471   hist = (TH1D *)(hist2_->Clone());
1472   hist->Write();
1473   TProfile *prof = (TProfile *)(prof_->Clone());
1474   prof->Write();
1475   std::cout << "All done" << std::endl;
1476   theFile->Close();
1477 }
1478 
1479 class CalibFitPU {
1480 private:
1481   TTree *fChain;   //!pointer to the analyzed TTree or TChain
1482   Int_t fCurrent;  //!current Tree number in a TChain
1483 
1484   // Fixed size dimensions of array or collections stored in the TTree if any.
1485 
1486   // Declaration of leaf types
1487   Int_t t_Event;
1488   Double_t t_p_PU;
1489   Double_t t_eHcal_PU;
1490   Double_t t_delta_PU;
1491   Double_t t_p_NoPU;
1492   Double_t t_eHcal_noPU;
1493   Double_t t_delta_NoPU;
1494   Int_t t_ieta;
1495 
1496   // List of branches
1497   TBranch *b_t_Event;
1498   TBranch *b_t_p_PU;
1499   TBranch *b_t_eHcal_PU;
1500   TBranch *b_t_delta_PU;
1501   TBranch *b_t_p_NoPU;
1502   TBranch *b_t_eHcal_noPU;
1503   TBranch *b_t_delta_NoPU;
1504   TBranch *b_t_ieta;
1505 
1506 public:
1507   CalibFitPU(const char *fname = "isotrackRelval.root");
1508   virtual ~CalibFitPU();
1509   virtual Int_t Cut(Long64_t entry);
1510   virtual Int_t GetEntry(Long64_t entry);
1511   virtual Long64_t LoadTree(Long64_t entry);
1512   virtual void Init(TTree *tree);
1513   virtual void Loop(bool extract_PU_parameters, std::string fileName, int dumpmax = 0);
1514   virtual Bool_t Notify();
1515   virtual void Show(Long64_t entry = -1);
1516 };
1517 
1518 CalibFitPU::CalibFitPU(const char *fname) : fChain(0) {
1519   // if parameter tree is not specified (or zero), connect the file
1520   // used to generate this class and read the Tree.
1521   TFile *f = new TFile(fname);
1522   TTree *tree = new TTree();
1523   f->GetObject("tree", tree);
1524   std::cout << "Find tree Tree in " << tree << " from " << fname << std::endl;
1525   Init(tree);
1526 }
1527 
1528 CalibFitPU::~CalibFitPU() {
1529   if (!fChain)
1530     return;
1531   delete fChain->GetCurrentFile();
1532 }
1533 
1534 Int_t CalibFitPU::GetEntry(Long64_t entry) {
1535   // Read contents of entry.
1536   if (!fChain)
1537     return 0;
1538   return fChain->GetEntry(entry);
1539 }
1540 
1541 Long64_t CalibFitPU::LoadTree(Long64_t entry) {
1542   // Set the environment to read one entry
1543   if (!fChain)
1544     return -5;
1545   Long64_t centry = fChain->LoadTree(entry);
1546   if (centry < 0)
1547     return centry;
1548   if (fChain->GetTreeNumber() != fCurrent) {
1549     fCurrent = fChain->GetTreeNumber();
1550     Notify();
1551   }
1552   return centry;
1553 }
1554 
1555 void CalibFitPU::Init(TTree *tree) {
1556   // The Init() function is called when the selector needs to initialize
1557   // a new tree or chain. Typically here the branch addresses and branch
1558   // pointers of the tree will be set.
1559   // It is normally not necessary to make changes to the generated
1560   // code, but the routine can be extended by the user if needed.
1561   // Init() will be called many times when running on PROOF
1562   // (once per file to be processed).
1563 
1564   // Set branch addresses and branch pointers
1565   if (!tree)
1566     return;
1567   fChain = tree;
1568   fCurrent = -1;
1569   fChain->SetMakeClass(1);
1570 
1571   fChain->SetBranchAddress("t_Event", &t_Event, &b_t_Event);
1572   fChain->SetBranchAddress("t_p_PU", &t_p_PU, &b_t_p_PU);
1573   fChain->SetBranchAddress("t_eHcal_PU", &t_eHcal_PU, &b_t_eHcal_PU);
1574   fChain->SetBranchAddress("t_delta_PU", &t_delta_PU, &b_t_delta_PU);
1575   fChain->SetBranchAddress("t_p_NoPU", &t_p_NoPU, &b_t_p_NoPU);
1576   fChain->SetBranchAddress("t_eHcal_noPU", &t_eHcal_noPU, &b_t_eHcal_noPU);
1577   fChain->SetBranchAddress("t_delta_NoPU", &t_delta_NoPU, &b_t_delta_NoPU);
1578   fChain->SetBranchAddress("t_ieta", &t_ieta, &b_t_ieta);
1579   Notify();
1580 }
1581 
1582 Bool_t CalibFitPU::Notify() {
1583   // The Notify() function is called when a new file is opened. This
1584   // can be either for a new TTree in a TChain or when when a new TTree
1585   // is started when using PROOF. It is normally not necessary to make changes
1586   // to the generated code, but the routine can be extended by the
1587   // user if needed. The return value is currently not used.
1588 
1589   return kTRUE;
1590 }
1591 
1592 void CalibFitPU::Show(Long64_t entry) {
1593   // Print contents of entry.
1594   // If entry is not specified, print current entry
1595   if (!fChain)
1596     return;
1597   fChain->Show(entry);
1598 }
1599 
1600 Int_t CalibFitPU::Cut(Long64_t) {
1601   // This function may be called from Loop.
1602   // returns  1 if entry is accepted.
1603   // returns -1 otherwise.
1604   return 1;
1605 }
1606 
1607 void CalibFitPU::Loop(bool extract_PU_parameters, std::string fileName, int dumpmax) {
1608   //   In a ROOT session, you can do:
1609   //      root> .L CalibFitPU.C
1610   //      root> CalibFitPU t
1611   //      root> t.GetEntry(12); // Fill t data members with entry number 12
1612   //      root> t.Show();       // Show values of entry 12
1613   //      root> t.Show(16);     // Read and show values of entry 16
1614   //      root> t.Loop();       // Loop on all entries
1615   //
1616 
1617   //     This is the loop skeleton where:
1618   //    jentry is the global entry number in the chain
1619   //    ientry is the entry number in the current Tree
1620   //  Note that the argument to GetEntry must be:
1621   //    jentry for TChain::GetEntry
1622   //    ientry for TTree::GetEntry and TBranch::GetEntry
1623   //
1624   //       To read only selected branches, Insert statements like:
1625   // METHOD1:
1626   //    fChain->SetBranchStatus("*",0);  // disable all branches
1627   //    fChain->SetBranchStatus("branchname",1);  // activate branchname
1628   // METHOD2: replace line
1629   //    fChain->GetEntry(jentry);       //read all branches
1630   //by  b_branchname->GetEntry(ientry); //read only this branch
1631   if (fChain == 0)
1632     return;
1633 
1634   Long64_t nentries = fChain->GetEntriesFast();
1635   Long64_t nbytes = 0, nb = 0;
1636 
1637   char filename[100];
1638 
1639   gStyle->SetCanvasBorderMode(0);
1640   gStyle->SetCanvasColor(kWhite);
1641   gStyle->SetPadColor(kWhite);
1642   gStyle->SetFillColor(kWhite);
1643 
1644   const int n = 7;
1645   int ieta_grid[n] = {7, 16, 25, 26, 27, 28, 29};
1646   double a00[n], a10[n], a20[n], a01[n], a11[n], a21[n], a02[n], a12[n], a22[n];
1647   const int nbin1 = 15;
1648   double bins1[nbin1 + 1] = {0.0, 0.1, 0.2, 0.3, 0.4, 0.6, 0.8, 1.0, 1.5, 2.0, 2.5, 3.0, 3.5, 4.0, 6.0, 9.0};
1649   const int nbin2 = 7;
1650   double bins2[nbin1 + 1] = {1.0, 1.5, 2.0, 2.5, 3.0, 3.5, 4.0, 5.0};
1651   int nbins[n] = {4, 5, 12, nbin2, nbin2, nbin2, nbin2};
1652   char name[20], title[100];
1653 
1654   std::vector<TH2D *> vec_h2;
1655   std::vector<TGraph *> vec_gr;
1656   std::vector<TProfile *> vec_hp;
1657 
1658   if (extract_PU_parameters) {
1659     for (int k = 0; k < n; k++) {
1660       sprintf(name, "h2_ieta%d", k);
1661       int ieta1 = (k == 0) ? 1 : ieta_grid[k - 1];
1662       sprintf(title, "PU Energy vs #Delta/p (i#eta = %d:%d)", ieta1, (ieta_grid[k] - 1));
1663       vec_h2.push_back(new TH2D(name, title, 100, 0, 10, 100, 0, 2));
1664       vec_gr.push_back(new TGraph());
1665       sprintf(name, "hp_ieta%d", k);
1666       if (k < 3)
1667         vec_hp.push_back(new TProfile(name, title, nbins[k], bins1));
1668       else
1669         vec_hp.push_back(new TProfile(name, title, nbins[k], bins2));
1670     }
1671 
1672     int points[7] = {0, 0, 0, 0, 0, 0, 0};
1673     //=======================================Starting of event Loop=======================================================
1674     std::cout << "Get " << nentries << " events from file\n";
1675 
1676     for (Long64_t jentry = 0; jentry < nentries; jentry++) {
1677       Long64_t ientry = LoadTree(jentry);
1678       if (ientry < 0)
1679         break;
1680       nb = fChain->GetEntry(jentry);
1681       nbytes += nb;
1682 
1683       double deltaOvP = t_delta_PU / t_p_PU;
1684       double diffEpuEnopuOvP = t_eHcal_noPU / t_eHcal_PU;
1685       if (jentry < dumpmax)
1686         std::cout << "Entry " << jentry << " ieta " << t_ieta << " p " << t_delta_PU << ":" << t_p_PU << ":" << deltaOvP
1687                   << " ehcal " << t_eHcal_noPU << ":" << t_eHcal_PU << ":" << diffEpuEnopuOvP << std::endl;
1688       for (int k = 0; k < n; k++) {
1689         if (std::abs(t_ieta) < ieta_grid[k]) {
1690           points[k]++;
1691           vec_h2[k]->Fill(deltaOvP, diffEpuEnopuOvP);
1692           vec_gr[k]->SetPoint(points[k], deltaOvP, diffEpuEnopuOvP);
1693           vec_hp[k]->Fill(deltaOvP, diffEpuEnopuOvP);
1694           break;
1695         }
1696       }
1697 
1698     }  //End of Event Loop to extract PU correction parameters
1699 
1700     for (int k = 0; k < n; k++)
1701       std::cout << "Bin " << k << " for 2D Hist " << vec_h2[k] << ":" << vec_h2[k]->GetEntries() << " Graph "
1702                 << vec_gr[k] << ":" << points[k] << " Profile " << vec_hp[k] << ":" << vec_hp[k]->GetEntries()
1703                 << std::endl;
1704 
1705     std::ofstream myfile0, myfile1, myfile2;
1706     sprintf(filename, "%s_par2d.txt", fileName.c_str());
1707     myfile0.open(filename);
1708     sprintf(filename, "%s_parProf.txt", fileName.c_str());
1709     myfile1.open(filename);
1710     sprintf(filename, "%s_parGraph.txt", fileName.c_str());
1711     myfile2.open(filename);
1712 
1713     char namepng[20];
1714     for (int k = 0; k < n; k++) {
1715       gStyle->SetOptStat(1100);
1716       gStyle->SetOptFit(1);
1717 
1718       TF1 *f1 = ((k < 2) ? (new TF1("f1", "[0]+[1]*x", 0, 5)) : (new TF1("f1", "[0]+[1]*x+[2]*x*x", 0, 5)));
1719       sprintf(name, "c_ieta%d2D", k);
1720       TCanvas *pad1 = new TCanvas(name, name, 500, 500);
1721       pad1->SetLeftMargin(0.10);
1722       pad1->SetRightMargin(0.10);
1723       pad1->SetTopMargin(0.10);
1724       vec_h2[k]->GetXaxis()->SetLabelSize(0.035);
1725       vec_h2[k]->GetYaxis()->SetLabelSize(0.035);
1726       vec_h2[k]->GetXaxis()->SetTitle("#Delta/p");
1727       vec_h2[k]->GetYaxis()->SetTitle("E_{NoPU} / E_{PU}");
1728       vec_h2[k]->Draw();
1729       vec_h2[k]->Fit(f1);
1730 
1731       a00[k] = f1->GetParameter(0);
1732       a10[k] = f1->GetParameter(1);
1733       a20[k] = (k < 2) ? 0 : f1->GetParameter(2);
1734       myfile0 << k << "\t" << a00[k] << "\t" << a10[k] << "\t" << a20[k] << "\n";
1735       pad1->Update();
1736       TPaveStats *st1 = (TPaveStats *)vec_h2[k]->GetListOfFunctions()->FindObject("stats");
1737       if (st1 != nullptr) {
1738         st1->SetY1NDC(0.70);
1739         st1->SetY2NDC(0.90);
1740         st1->SetX1NDC(0.65);
1741         st1->SetX2NDC(0.90);
1742       }
1743       pad1->Update();
1744       sprintf(namepng, "%s.png", pad1->GetName());
1745       pad1->Print(namepng);
1746 
1747       TF1 *f2 = ((k < 2) ? (new TF1("f2", "[0]+[1]*x", 0, 5)) : (new TF1("f2", "[0]+[1]*x+[2]*x*x", 0, 5)));
1748       sprintf(name, "c_ieta%dPr", k);
1749       TCanvas *pad2 = new TCanvas(name, name, 500, 500);
1750       pad2->SetLeftMargin(0.10);
1751       pad2->SetRightMargin(0.10);
1752       pad2->SetTopMargin(0.10);
1753       vec_hp[k]->GetXaxis()->SetLabelSize(0.035);
1754       vec_hp[k]->GetYaxis()->SetLabelSize(0.035);
1755       vec_hp[k]->GetXaxis()->SetTitle("#Delta/p");
1756       vec_hp[k]->GetYaxis()->SetTitle("E_{NoPU} / E_{PU}");
1757       vec_hp[k]->Draw();
1758       vec_hp[k]->Fit(f2);
1759 
1760       a01[k] = f2->GetParameter(0);
1761       a11[k] = f2->GetParameter(1);
1762       a21[k] = (k < 2) ? 0 : f2->GetParameter(2);
1763       myfile1 << k << "\t" << a01[k] << "\t" << a11[k] << "\t" << a21[k] << "\n";
1764       pad2->Update();
1765       TPaveStats *st2 = (TPaveStats *)vec_hp[k]->GetListOfFunctions()->FindObject("stats");
1766       if (st2 != nullptr) {
1767         st2->SetY1NDC(0.70);
1768         st2->SetY2NDC(0.90);
1769         st2->SetX1NDC(0.65);
1770         st2->SetX2NDC(0.90);
1771       }
1772       pad2->Update();
1773       sprintf(namepng, "%s.png", pad2->GetName());
1774       pad2->Print(namepng);
1775 
1776       TF1 *f3 = ((k < 2) ? (new TF1("f3", "[0]+[1]*x", 0, 5)) : (new TF1("f3", "[0]+[1]*x+[2]*x*x", 0, 5)));
1777       sprintf(name, "c_ieta%dGr", k);
1778       TCanvas *pad3 = new TCanvas(name, name, 500, 500);
1779       pad3->SetLeftMargin(0.10);
1780       pad3->SetRightMargin(0.10);
1781       pad3->SetTopMargin(0.10);
1782       gStyle->SetOptFit(1111);
1783       vec_gr[k]->GetXaxis()->SetLabelSize(0.035);
1784       vec_gr[k]->GetYaxis()->SetLabelSize(0.035);
1785       vec_gr[k]->GetXaxis()->SetTitle("#Delta/p");
1786       vec_gr[k]->GetYaxis()->SetTitle("E_{PU} - E_{NoPU}/p");
1787       vec_gr[k]->Fit(f3, "R");
1788       vec_gr[k]->Draw("Ap");
1789       f3->Draw("same");
1790 
1791       a02[k] = f3->GetParameter(0);
1792       a12[k] = f3->GetParameter(1);
1793       a22[k] = (k < 2) ? 0 : f3->GetParameter(2);
1794       myfile2 << k << "\t" << a02[k] << "\t" << a12[k] << "\t" << a22[k] << "\n";
1795       pad3->Update();
1796       TPaveStats *st3 = (TPaveStats *)vec_gr[k]->GetListOfFunctions()->FindObject("stats");
1797       if (st3 != nullptr) {
1798         st3->SetY1NDC(0.70);
1799         st3->SetY2NDC(0.90);
1800         st3->SetX1NDC(0.65);
1801         st3->SetX2NDC(0.90);
1802       }
1803       pad3->Update();
1804       pad3->Modified();
1805       sprintf(namepng, "%s.png", pad3->GetName());
1806       pad3->Print(namepng);
1807     }
1808   } else {
1809     std::string line;
1810     double number;
1811     sprintf(filename, "%s_par2d.txt", fileName.c_str());
1812     std::ifstream myfile0(filename);
1813     if (myfile0.is_open()) {
1814       int iii = 0;
1815       while (getline(myfile0, line)) {
1816         std::istringstream iss2(line);
1817         int ii = 0;
1818         while (iss2 >> number) {
1819           if (ii == 0)
1820             a00[iii] = number;
1821           else if (ii == 1)
1822             a10[iii] = number;
1823           else if (ii == 2)
1824             a20[iii] = number;
1825           ++ii;
1826         }
1827         ++iii;
1828       }
1829     }
1830     sprintf(filename, "%s_parProf.txt", fileName.c_str());
1831     std::ifstream myfile1(filename);
1832     if (myfile1.is_open()) {
1833       int iii = 0;
1834       while (getline(myfile1, line)) {
1835         std::istringstream iss2(line);
1836         int ii = 0;
1837         while (iss2 >> number) {
1838           if (ii == 0)
1839             a01[iii] = number;
1840           else if (ii == 1)
1841             a11[iii] = number;
1842           else if (ii == 2)
1843             a21[iii] = number;
1844           ++ii;
1845         }
1846         ++iii;
1847       }
1848     }
1849     sprintf(filename, "%s_parGraph.txt", fileName.c_str());
1850     std::ifstream myfile2(filename);
1851     if (myfile2.is_open()) {
1852       int iii = 0;
1853       while (getline(myfile2, line)) {
1854         std::istringstream iss2(line);
1855         int ii = 0;
1856         while (iss2 >> number) {
1857           if (ii == 0)
1858             a02[iii] = number;
1859           else if (ii == 1)
1860             a12[iii] = number;
1861           else if (ii == 2)
1862             a22[iii] = number;
1863           ++ii;
1864         }
1865         ++iii;
1866       }
1867     }
1868   }
1869 
1870   std::cout << "\nParameter Values:\n";
1871   for (int i = 0; i < n; i++) {
1872     std::cout << "[" << i << "] (" << a00[i] << ", " << a10[i] << ", " << a20[i] << ")  (" << a01[i] << ", " << a11[i]
1873               << ", " << a21[i] << ")  (" << a02[i] << ", " << a12[i] << ", " << a22[i] << ")\n";
1874   }
1875   std::cout << "\n\n";
1876   std::vector<TH1F *> vec_EovP_UnCorr_PU, vec_EovP_UnCorr_NoPU;
1877   std::vector<TH1F *> vec_EovP_Corr0_PU, vec_EovP_Corr0_NoPU;
1878   std::vector<TH1F *> vec_EovP_Corr1_PU, vec_EovP_Corr1_NoPU;
1879   std::vector<TH1F *> vec_EovP_Corr2_PU, vec_EovP_Corr2_NoPU;
1880   TH1F *h_current;
1881 
1882   for (int k = 0; k < (2 * 28 + 1); k++) {
1883     if (k != 28) {
1884       sprintf(name, "EovP_ieta%dUnPU", k - 28);
1885       sprintf(title, "E/p (Uncorrected PU) for i#eta = %d", k - 28);
1886       h_current = new TH1F(name, title, 100, 0, 5);
1887       h_current->GetXaxis()->SetTitle(title);
1888       h_current->GetYaxis()->SetTitle("Tracks");
1889       vec_EovP_UnCorr_PU.push_back(h_current);
1890       sprintf(name, "EovP_ieta%dUnNoPU", k - 28);
1891       sprintf(title, "E/p (Uncorrected No PU) for i#eta = %d", k - 28);
1892       h_current = new TH1F(name, title, 100, 0, 5);
1893       h_current->GetXaxis()->SetTitle(title);
1894       h_current->GetYaxis()->SetTitle("Tracks");
1895       vec_EovP_UnCorr_NoPU.push_back(h_current);
1896       sprintf(name, "EovP_ieta%dCor0NoPU", k - 28);
1897       sprintf(title, "E/p (Corrected using 2D No PU) for i#eta = %d", k - 28);
1898       h_current = new TH1F(name, title, 100, 0, 5);
1899       h_current->GetXaxis()->SetTitle(title);
1900       h_current->GetYaxis()->SetTitle("Tracks");
1901       vec_EovP_Corr0_NoPU.push_back(h_current);
1902       sprintf(name, "EovP_ieta%dCor0PU", k - 28);
1903       sprintf(title, "E/p (Corrected using 2D PU) for i#eta = %d", k - 28);
1904       h_current = new TH1F(name, title, 100, 0, 5);
1905       h_current->GetXaxis()->SetTitle(title);
1906       h_current->GetYaxis()->SetTitle("Tracks");
1907       vec_EovP_Corr0_PU.push_back(h_current);
1908       sprintf(name, "EovP_ieta%dCor1NoPU", k - 28);
1909       sprintf(title, "E/p (Corrected using profile No PU) for i#eta = %d", k - 28);
1910       h_current = new TH1F(name, title, 100, 0, 5);
1911       h_current->GetXaxis()->SetTitle(title);
1912       h_current->GetYaxis()->SetTitle("Tracks");
1913       vec_EovP_Corr1_NoPU.push_back(h_current);
1914       sprintf(name, "EovP_ieta%dCor1PU", k - 28);
1915       sprintf(title, "E/p (Corrected using profile PU) for i#eta = %d", k - 28);
1916       h_current = new TH1F(name, title, 100, 0, 5);
1917       h_current->GetXaxis()->SetTitle(title);
1918       h_current->GetYaxis()->SetTitle("Tracks");
1919       vec_EovP_Corr1_PU.push_back(h_current);
1920       sprintf(name, "EovP_ieta%dCor2NoPU", k - 28);
1921       sprintf(title, "E/p (Corrected using graph No PU) for i#eta = %d", k - 28);
1922       h_current = new TH1F(name, title, 100, 0, 5);
1923       h_current->GetXaxis()->SetTitle(title);
1924       h_current->GetYaxis()->SetTitle("Tracks");
1925       vec_EovP_Corr2_NoPU.push_back(h_current);
1926       sprintf(name, "EovP_ieta%dCor2PU", k - 28);
1927       sprintf(title, "E/p (Corrected using graph PU) for i#eta = %d", k - 28);
1928       h_current = new TH1F(name, title, 100, 0, 5);
1929       h_current->GetXaxis()->SetTitle(title);
1930       h_current->GetYaxis()->SetTitle("Tracks");
1931       vec_EovP_Corr2_PU.push_back(h_current);
1932     }
1933   }
1934   std::cout << "Book " << (8 * vec_EovP_UnCorr_PU.size()) << " histograms\n";
1935 
1936   //==============================================================================================================================================
1937 
1938   nbytes = nb = 0;
1939   std::cout << nentries << " entries in the root file" << std::endl;
1940   for (Long64_t jentry = 0; jentry < nentries; jentry++) {
1941     Long64_t ientry = LoadTree(jentry);
1942     if (ientry < 0)
1943       break;
1944     nb = fChain->GetEntry(jentry);
1945     nbytes += nb;
1946 
1947     int i1 = n;
1948     for (int k = 0; k < n; k++) {
1949       if (std::abs(t_ieta) < ieta_grid[k]) {
1950         i1 = k;
1951         break;
1952       }
1953     }
1954     if ((t_ieta == 0) || (i1 >= n))
1955       continue;
1956     int i2 = (t_ieta < 0) ? (t_ieta + 28) : (t_ieta + 27);
1957 
1958     double EpuOvP = t_eHcal_PU / t_p_PU;
1959     double EnopuOvP = t_eHcal_noPU / t_p_PU;
1960     double deltaOvP = t_delta_PU / t_p_PU;
1961     double deltaNopuOvP = t_delta_NoPU / t_p_PU;
1962 
1963     vec_EovP_UnCorr_PU[i2]->Fill(EpuOvP);
1964     vec_EovP_UnCorr_NoPU[i2]->Fill(EnopuOvP);
1965 
1966     double c0p =
1967         (((i1 < 3) || (deltaOvP > 1.0)) ? (a00[i1] + a10[i1] * deltaOvP + a20[i1] * deltaOvP * deltaOvP) : 1.0);
1968     double c1p =
1969         (((i1 < 3) || (deltaOvP > 1.0)) ? (a01[i1] + a11[i1] * deltaOvP + a21[i1] * deltaOvP * deltaOvP) : 1.0);
1970     double c2p =
1971         (((i1 < 3) || (deltaOvP > 1.0)) ? (a02[i1] + a12[i1] * deltaOvP + a22[i1] * deltaOvP * deltaOvP) : 1.0);
1972     double c0np =
1973         (((i1 < 3) || (deltaNopuOvP > 1.0)) ? (a00[i1] + a10[i1] * deltaNopuOvP + a12[i1] * deltaNopuOvP * deltaNopuOvP)
1974                                             : 1.0);
1975     double c1np =
1976         (((i1 < 3) || (deltaNopuOvP > 1.0)) ? (a01[i1] + a11[i1] * deltaNopuOvP + a21[i1] * deltaNopuOvP * deltaNopuOvP)
1977                                             : 1.0);
1978     double c2np =
1979         (((i1 < 3) || (deltaNopuOvP > 1.0)) ? (a02[i1] + a12[i1] * deltaNopuOvP + a22[i1] * deltaNopuOvP * deltaNopuOvP)
1980                                             : 1.0);
1981 
1982     vec_EovP_Corr0_PU[i2]->Fill(EpuOvP * c0p);
1983     vec_EovP_Corr0_NoPU[i2]->Fill(EnopuOvP * c0np);
1984     vec_EovP_Corr1_PU[i2]->Fill(EpuOvP * c1p);
1985     vec_EovP_Corr1_NoPU[i2]->Fill(EnopuOvP * c1np);
1986     vec_EovP_Corr2_PU[i2]->Fill(EpuOvP * c2p);
1987     vec_EovP_Corr2_NoPU[i2]->Fill(EnopuOvP * c2np);
1988   }
1989 
1990   sprintf(filename, "%s.root", fileName.c_str());
1991   TFile *f1 = new TFile(filename, "RECREATE");
1992   f1->cd();
1993   for (unsigned int k = 0; k < vec_EovP_UnCorr_PU.size(); k++) {
1994     vec_EovP_UnCorr_PU[k]->Write();
1995     vec_EovP_UnCorr_NoPU[k]->Write();
1996     vec_EovP_Corr0_PU[k]->Write();
1997     vec_EovP_Corr0_NoPU[k]->Write();
1998     vec_EovP_Corr1_PU[k]->Write();
1999     vec_EovP_Corr1_NoPU[k]->Write();
2000     vec_EovP_Corr2_PU[k]->Write();
2001     vec_EovP_Corr2_NoPU[k]->Write();
2002   }
2003   for (unsigned int k = 0; k < vec_hp.size(); ++k) {
2004     vec_h2[k]->Write();
2005     vec_hp[k]->Write();
2006     vec_gr[k]->Write();
2007   }
2008 }
2009 
2010 class CalibMerge {
2011 public:
2012   struct isotk {
2013     isotk(int pv = 0, double rho = 0, double pp = 0, double eH = 0, double eH10 = 0, double eH30 = 0, double del = 0)
2014         : goodPV(pv), rhoh(rho), p(pp), eHcal(eH), eHcal10(eH10), eHcal30(eH30), delta(del) {}
2015 
2016     int goodPV;
2017     double rhoh, p, eHcal, eHcal10, eHcal30, delta;
2018   };
2019 
2020   CalibMerge();
2021   ~CalibMerge();
2022   Int_t Cut(Long64_t entry);
2023   Int_t GetEntry(Long64_t entry);
2024   Long64_t LoadTree(Long64_t entry);
2025   Bool_t Notify();
2026   void Show(Long64_t entry = -1);
2027   void LoopNoPU(const char *fnameNoPU, std::string dirnm = "HcalIsoTrkAnalyzer");
2028   void LoopPU(const char *fnamePU, const char *fnameOut, std::string dirnm = "HcalIsoTrkAnalyzer");
2029   bool Init(const char *fname, const std::string &dirnm);
2030   void bookTree(const char *fname);
2031   void close();
2032 
2033 private:
2034   TChain *fChain;  //!pointer to the analyzed TTree or TChain
2035   Int_t fCurrent;  //!current Tree number in a TChain
2036 
2037   std::map<std::pair<int, int>, std::vector<isotk> > isotks_;
2038 
2039   // Declaration of leaf types
2040   Int_t t_Run;
2041   Int_t t_Event;
2042   Int_t t_DataType;
2043   Int_t t_ieta;
2044   Int_t t_iphi;
2045   Double_t t_EventWeight;
2046   Int_t t_nVtx;
2047   Int_t t_nTrk;
2048   Int_t t_goodPV;
2049   Double_t t_l1pt;
2050   Double_t t_l1eta;
2051   Double_t t_l1phi;
2052   Double_t t_l3pt;
2053   Double_t t_l3eta;
2054   Double_t t_l3phi;
2055   Double_t t_p;
2056   Double_t t_pt;
2057   Double_t t_phi;
2058   Double_t t_mindR1;
2059   Double_t t_mindR2;
2060   Double_t t_eMipDR;
2061   Double_t t_eHcal;
2062   Double_t t_eHcal10;
2063   Double_t t_eHcal30;
2064   Double_t t_hmaxNearP;
2065   Double_t t_rhoh;
2066   Bool_t t_selectTk;
2067   Bool_t t_qltyFlag;
2068   Bool_t t_qltyMissFlag;
2069   Bool_t t_qltyPVFlag;
2070   Double_t t_gentrackP;
2071   std::vector<unsigned int> *t_DetIds;
2072   std::vector<double> *t_HitEnergies;
2073   std::vector<bool> *t_trgbits;
2074   std::vector<unsigned int> *t_DetIds1;
2075   std::vector<unsigned int> *t_DetIds3;
2076   std::vector<double> *t_HitEnergies1;
2077   std::vector<double> *t_HitEnergies3;
2078 
2079   // List of branches
2080   TBranch *b_t_Run;           //!
2081   TBranch *b_t_Event;         //!
2082   TBranch *b_t_DataType;      //!
2083   TBranch *b_t_ieta;          //!
2084   TBranch *b_t_iphi;          //!
2085   TBranch *b_t_EventWeight;   //!
2086   TBranch *b_t_nVtx;          //!
2087   TBranch *b_t_nTrk;          //!
2088   TBranch *b_t_goodPV;        //!
2089   TBranch *b_t_l1pt;          //!
2090   TBranch *b_t_l1eta;         //!
2091   TBranch *b_t_l1phi;         //!
2092   TBranch *b_t_l3pt;          //!
2093   TBranch *b_t_l3eta;         //!
2094   TBranch *b_t_l3phi;         //!
2095   TBranch *b_t_p;             //!
2096   TBranch *b_t_pt;            //!
2097   TBranch *b_t_phi;           //!
2098   TBranch *b_t_mindR1;        //!
2099   TBranch *b_t_mindR2;        //!
2100   TBranch *b_t_eMipDR;        //!
2101   TBranch *b_t_eHcal;         //!
2102   TBranch *b_t_eHcal10;       //!
2103   TBranch *b_t_eHcal30;       //!
2104   TBranch *b_t_hmaxNearP;     //!
2105   TBranch *b_t_rhoh;          //!
2106   TBranch *b_t_selectTk;      //!
2107   TBranch *b_t_qltyFlag;      //!
2108   TBranch *b_t_qltyMissFlag;  //!
2109   TBranch *b_t_qltyPVFlag;    //!
2110   TBranch *b_t_gentrackP;     //!
2111   TBranch *b_t_DetIds;        //!
2112   TBranch *b_t_HitEnergies;   //!
2113   TBranch *b_t_trgbits;       //!
2114   TBranch *b_t_DetIds1;       //!
2115   TBranch *b_t_DetIds3;       //!
2116   TBranch *b_t_HitEnergies1;  //!
2117   TBranch *b_t_HitEnergies3;  //!
2118 
2119   TFile *outputFile_;
2120   TTree *outtree_;
2121   Int_t event_;
2122   Int_t ieta_;
2123   Int_t nvtx_;
2124   Double_t rhoh_;
2125   Double_t p_NoPU_;
2126   Double_t p_PU_;
2127   Double_t eHcal_NoPU_;
2128   Double_t eHcal_PU_;
2129   Double_t eHcal10_NoPU_;
2130   Double_t eHcal10_PU_;
2131   Double_t eHcal30_NoPU_;
2132   Double_t eHcal30_PU_;
2133   Double_t delta_NoPU_;
2134   Double_t delta_PU_;
2135 };
2136 
2137 CalibMerge::CalibMerge() : outputFile_(nullptr), outtree_(nullptr) {}
2138 
2139 CalibMerge::~CalibMerge() { close(); }
2140 
2141 Int_t CalibMerge::GetEntry(Long64_t entry) {
2142   // Read contents of entry.
2143   if (!fChain)
2144     return 0;
2145   return fChain->GetEntry(entry);
2146 }
2147 
2148 Long64_t CalibMerge::LoadTree(Long64_t entry) {
2149   // Set the environment to read one entry
2150   if (!fChain)
2151     return -5;
2152   Long64_t centry = fChain->LoadTree(entry);
2153   if (centry < 0)
2154     return centry;
2155   if (fChain->GetTreeNumber() != fCurrent) {
2156     fCurrent = fChain->GetTreeNumber();
2157     Notify();
2158   }
2159   return centry;
2160 }
2161 
2162 Bool_t CalibMerge::Notify() {
2163   // The Notify() function is called when a new file is opened. This
2164   // can be either for a new TTree in a TChain or when when a new TTree
2165   // is started when using PROOF. It is normally not necessary to make changes
2166   // to the generated code, but the routine can be extended by the
2167   // user if needed. The return value is currently not used.
2168 
2169   return kTRUE;
2170 }
2171 
2172 void CalibMerge::Show(Long64_t entry) {
2173   // Print contents of entry.
2174   // If entry is not specified, print current entry
2175   if (!fChain)
2176     return;
2177   fChain->Show(entry);
2178 }
2179 
2180 Int_t CalibMerge::Cut(Long64_t) {
2181   // This function may be called from Loop.
2182   // returns  1 if entry is accepted.
2183   // returns -1 otherwise.
2184   return 1;
2185 }
2186 
2187 void CalibMerge::LoopNoPU(const char *fname, std::string dirnm) {
2188   bool ok = Init(fname, dirnm);
2189   std::cout << "Opens no PU file " << fname << " with flag " << ok << std::endl;
2190   isotks_.clear();
2191   if (ok) {
2192     Long64_t nentries = fChain->GetEntriesFast();
2193     Long64_t nbytes = 0, nb = 0;
2194     for (Long64_t jentry = 0; jentry < nentries; jentry++) {
2195       Long64_t ientry = LoadTree(jentry);
2196       if (ientry < 0)
2197         break;
2198       nb = fChain->GetEntry(jentry);
2199       nbytes += nb;
2200       std::pair<int, int> key(t_Event, t_ieta);
2201       isotk tk(t_goodPV, t_rhoh, t_p, t_eHcal, t_eHcal10, t_eHcal30, (t_eHcal30 - t_eHcal10));
2202       std::map<std::pair<int, int>, std::vector<isotk> >::iterator itr = isotks_.find(key);
2203       if (itr == isotks_.end()) {
2204         std::vector<isotk> v;
2205         v.push_back(tk);
2206         isotks_[key] = v;
2207       } else {
2208         (itr->second).push_back(tk);
2209       }
2210     }
2211     std::cout << "Finds " << isotks_.size() << " tracks from " << nentries << " entries with " << nbytes << " bytes"
2212               << std::endl;
2213   }
2214 }
2215 
2216 void CalibMerge::LoopPU(const char *fname, const char *fout, std::string dirnm) {
2217   bool ok = Init(fname, dirnm);
2218   std::cout << "Opens PU file " << fname << " with flag " << ok << std::endl;
2219   if (ok && (isotks_.size() > 0)) {
2220     bookTree(fout);
2221     Long64_t nentries = fChain->GetEntriesFast();
2222     Long64_t nbytes(0), nb(0), merge(0);
2223     for (Long64_t jentry = 0; jentry < nentries; jentry++) {
2224       Long64_t ientry = LoadTree(jentry);
2225       if (ientry < 0)
2226         break;
2227       nb = fChain->GetEntry(jentry);
2228       nbytes += nb;
2229       std::pair<int, int> key(t_Event, t_ieta);
2230       std::map<std::pair<int, int>, std::vector<isotk> >::iterator itr = isotks_.find(key);
2231       if (itr != isotks_.end()) {
2232         for (unsigned int k = 0; (itr->second).size(); ++k) {
2233           if (std::abs(t_p - (itr->second)[k].p) < 0.01) {
2234             event_ = (itr->first).first;
2235             ieta_ = (itr->first).second;
2236             nvtx_ = (itr->second)[k].goodPV;
2237             rhoh_ = (itr->second)[k].rhoh;
2238             p_NoPU_ = (itr->second)[k].p;
2239             p_PU_ = t_p;
2240             eHcal_NoPU_ = (itr->second)[k].eHcal;
2241             eHcal_PU_ = t_eHcal;
2242             eHcal10_NoPU_ = (itr->second)[k].eHcal10;
2243             eHcal10_PU_ = t_eHcal10;
2244             eHcal30_NoPU_ = (itr->second)[k].eHcal30;
2245             eHcal30_PU_ = t_eHcal30;
2246             delta_NoPU_ = (itr->second)[k].delta;
2247             delta_PU_ = (t_eHcal30 - t_eHcal10);
2248             outtree_->Fill();
2249             ++merge;
2250             break;
2251           }
2252         }
2253       }
2254     }
2255     std::cout << "Use " << isotks_.size() << ":" << nentries << " from noPU|PU files to produce " << merge
2256               << " merged tracks";
2257   }
2258 }
2259 
2260 bool CalibMerge::Init(const char *fname, const std::string &dirnm) {
2261   char treeName[400];
2262   sprintf(treeName, "%s/CalibTree", dirnm.c_str());
2263   TChain *fChain = new TChain(treeName);
2264   if (!fillChain(fChain, fname)) {
2265     std::cout << "Creation of the chain for " << treeName << " from " << fname << " failed " << std::endl;
2266     return false;
2267   } else {
2268     std::cout << "Proceed with a tree chain with " << fChain->GetEntries() << " entries" << std::endl;
2269     fCurrent = -1;
2270     fChain->SetMakeClass(1);
2271     // Set object pointer
2272     t_DetIds = 0;
2273     t_DetIds1 = 0;
2274     t_DetIds3 = 0;
2275     t_HitEnergies = 0;
2276     t_HitEnergies1 = 0;
2277     t_HitEnergies3 = 0;
2278     t_trgbits = 0;
2279     fChain->SetBranchAddress("t_Run", &t_Run, &b_t_Run);
2280     fChain->SetBranchAddress("t_Event", &t_Event, &b_t_Event);
2281     fChain->SetBranchAddress("t_DataType", &t_DataType, &b_t_DataType);
2282     fChain->SetBranchAddress("t_ieta", &t_ieta, &b_t_ieta);
2283     fChain->SetBranchAddress("t_iphi", &t_iphi, &b_t_iphi);
2284     fChain->SetBranchAddress("t_EventWeight", &t_EventWeight, &b_t_EventWeight);
2285     fChain->SetBranchAddress("t_nVtx", &t_nVtx, &b_t_nVtx);
2286     fChain->SetBranchAddress("t_nTrk", &t_nTrk, &b_t_nTrk);
2287     fChain->SetBranchAddress("t_goodPV", &t_goodPV, &b_t_goodPV);
2288     fChain->SetBranchAddress("t_l1pt", &t_l1pt, &b_t_l1pt);
2289     fChain->SetBranchAddress("t_l1eta", &t_l1eta, &b_t_l1eta);
2290     fChain->SetBranchAddress("t_l1phi", &t_l1phi, &b_t_l1phi);
2291     fChain->SetBranchAddress("t_l3pt", &t_l3pt, &b_t_l3pt);
2292     fChain->SetBranchAddress("t_l3eta", &t_l3eta, &b_t_l3eta);
2293     fChain->SetBranchAddress("t_l3phi", &t_l3phi, &b_t_l3phi);
2294     fChain->SetBranchAddress("t_p", &t_p, &b_t_p);
2295     fChain->SetBranchAddress("t_pt", &t_pt, &b_t_pt);
2296     fChain->SetBranchAddress("t_phi", &t_phi, &b_t_phi);
2297     fChain->SetBranchAddress("t_mindR1", &t_mindR1, &b_t_mindR1);
2298     fChain->SetBranchAddress("t_mindR2", &t_mindR2, &b_t_mindR2);
2299     fChain->SetBranchAddress("t_eMipDR", &t_eMipDR, &b_t_eMipDR);
2300     fChain->SetBranchAddress("t_eHcal", &t_eHcal, &b_t_eHcal);
2301     fChain->SetBranchAddress("t_eHcal10", &t_eHcal10, &b_t_eHcal10);
2302     fChain->SetBranchAddress("t_eHcal30", &t_eHcal30, &b_t_eHcal30);
2303     fChain->SetBranchAddress("t_hmaxNearP", &t_hmaxNearP, &b_t_hmaxNearP);
2304     fChain->SetBranchAddress("t_rhoh", &t_rhoh, &b_t_rhoh);
2305     fChain->SetBranchAddress("t_selectTk", &t_selectTk, &b_t_selectTk);
2306     fChain->SetBranchAddress("t_qltyFlag", &t_qltyFlag, &b_t_qltyFlag);
2307     fChain->SetBranchAddress("t_qltyMissFlag", &t_qltyMissFlag, &b_t_qltyMissFlag);
2308     fChain->SetBranchAddress("t_qltyPVFlag", &t_qltyPVFlag, &b_t_qltyPVFlag);
2309     fChain->SetBranchAddress("t_gentrackP", &t_gentrackP, &b_t_gentrackP);
2310     fChain->SetBranchAddress("t_DetIds", &t_DetIds, &b_t_DetIds);
2311     fChain->SetBranchAddress("t_HitEnergies", &t_HitEnergies, &b_t_HitEnergies);
2312     fChain->SetBranchAddress("t_trgbits", &t_trgbits, &b_t_trgbits);
2313     fChain->SetBranchAddress("t_DetIds1", &t_DetIds1, &b_t_DetIds1);
2314     fChain->SetBranchAddress("t_DetIds3", &t_DetIds3, &b_t_DetIds3);
2315     fChain->SetBranchAddress("t_HitEnergies1", &t_HitEnergies1, &b_t_HitEnergies1);
2316     fChain->SetBranchAddress("t_HitEnergies3", &t_HitEnergies3, &b_t_HitEnergies3);
2317     Notify();
2318     return true;
2319   }
2320 }
2321 
2322 void CalibMerge::bookTree(const char *fname) {
2323   std::cout << "BookTree in " << fname << std::endl;
2324   outputFile_ = TFile::Open(fname, "RECREATE");
2325   outputFile_->cd();
2326   outtree_ = new TTree("Tree", "Tree");
2327   outtree_->Branch("t_Event", &event_);
2328   outtree_->Branch("t_ieta", &ieta_);
2329   outtree_->Branch("t_nvtx", &nvtx_);
2330   outtree_->Branch("t_rhoh", &rhoh_);
2331   outtree_->Branch("t_p_noPU", &p_NoPU_);
2332   outtree_->Branch("t_p_PU", &p_PU_);
2333   outtree_->Branch("t_eHcal_noPU", &eHcal_NoPU_);
2334   outtree_->Branch("t_eHcal_PU", &eHcal_PU_);
2335   outtree_->Branch("t_eHcal10_noPU", &eHcal10_NoPU_);
2336   outtree_->Branch("t_eHcal10_PU", &eHcal10_PU_);
2337   outtree_->Branch("t_eHcal30_noPU", &eHcal30_NoPU_);
2338   outtree_->Branch("t_eHcal30_PU", &eHcal30_PU_);
2339   outtree_->Branch("t_delta_noPU", &delta_NoPU_);
2340   outtree_->Branch("t_delta_PU", &delta_PU_);
2341 }
2342 
2343 void CalibMerge::close() {
2344   if (outputFile_ != nullptr) {
2345     outputFile_->cd();
2346     std::cout << "file yet to be Written" << std::endl;
2347     outtree_->Write();
2348     std::cout << "file Written" << std::endl;
2349     outputFile_->Close();
2350   }
2351   outputFile_ = nullptr;
2352   std::cout << "now doing return" << std::endl;
2353 }
2354 
2355 void combineML(const char *inputFileList, const char *outfile) {
2356   std::map<int, std::string> files;
2357   std::ifstream infile(inputFileList);
2358   if (!infile.is_open()) {
2359     std::cout << "** ERROR: Can't open '" << inputFileList << "' for input" << std::endl;
2360   } else {
2361     while (1) {
2362       int depth;
2363       std::string fname;
2364       infile >> depth >> fname;
2365       if (!infile.good())
2366         break;
2367       files[depth] = fname;
2368     }
2369     infile.close();
2370   }
2371   std::cout << "Gets a list of " << files.size() << " file names from " << inputFileList << std::endl;
2372   if (files.size() > 0) {
2373     std::map<int, std::vector<listML> > mlList;
2374     for (std::map<int, std::string>::const_iterator itr = files.begin(); itr != files.end(); ++itr) {
2375       int depth = itr->first;
2376       std::string fname = itr->second;
2377       std::ifstream fInput(fname.c_str());
2378       if (!fInput.good()) {
2379         std::cout << "Cannot open file " << fname << std::endl;
2380       } else {
2381         char buffer[1024];
2382         unsigned int all(0), good1(0), good2(0);
2383         while (fInput.getline(buffer, 1024)) {
2384           ++all;
2385           if (buffer[0] == '#')
2386             continue;  //ignore comment
2387           std::vector<std::string> items = splitString(std::string(buffer));
2388           if (items.size() != 4) {
2389             std::cout << "Ignore  line: " << buffer << std::endl;
2390           } else {
2391             ++good1;
2392             int depth0 = std::atoi(items[1].c_str());
2393             if (depth0 == depth) {
2394               ++good2;
2395               int ieta = std::atoi(items[0].c_str());
2396               double ml = std::atoi(items[2].c_str());
2397               double dml = std::atoi(items[3].c_str());
2398               listML l0(depth, ml, dml);
2399               if (mlList.find(ieta) == mlList.end()) {
2400                 std::vector<listML> l0v;
2401                 mlList[ieta] = l0v;
2402               }
2403               (mlList[ieta]).push_back(l0);
2404             }
2405           }
2406         }
2407         fInput.close();
2408         std::cout << "Reads total of " << all << " and " << good1 << ":" << good2 << " good records for depth " << depth
2409                   << std::endl;
2410       }
2411     }
2412     if (mlList.size() > 0) {
2413       std::ofstream fout(outfile);
2414       for (std::map<int, std::vector<listML> >::const_iterator itr = mlList.begin(); itr != mlList.end(); ++itr) {
2415         int ieta = itr->first;
2416         std::vector<listML> l0v = itr->second;
2417         double den(0), dden(0);
2418         for (unsigned int k = 0; k < l0v.size(); ++k) {
2419           if (l0v[k].depth_ == 2) {
2420             den = l0v[k].ml_;
2421             dden = l0v[k].dml_;
2422           }
2423         }
2424         if (den > 0) {
2425           fout << std::setw(4) << ieta << "   " << l0v.size();
2426           for (unsigned int k = 0; k < l0v.size(); ++k) {
2427             double ml = den / l0v[k].ml_;
2428             double dml = l0v[k].dml_ * den / (l0v[k].ml_ * l0v[k].ml_);
2429             fout << "  " << l0v[k].depth_ << "  " << std::setw(6) << ml << "  " << std::setw(6) << dml;
2430           }
2431           fout << std::endl;
2432         }
2433       }
2434       fout.close();
2435     }
2436   }
2437 }
2438 
2439 void calCorrCombine(
2440     const char *inFileCorr, const char *inFileDepth, const char *outFileCorr, int truncateFlag, int etaMax) {
2441   const int neta = 58;
2442   int ietas[neta] = {1,   2,   3,   4,   5,   6,   7,   8,   9,   10,  11,  12,  13,  14,  15,  16,  17,  18, 19,  20,
2443                      21,  22,  23,  24,  25,  26,  27,  28,  29,  -1,  -2,  -3,  -4,  -5,  -6,  -7,  -8,  -9, -10, -11,
2444                      -12, -13, -14, -15, -16, -17, -18, -19, -20, -21, -22, -23, -24, -25, -26, -27, -28, -29};
2445   int depthMin[neta] = {1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 2, 2, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1,
2446                         1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 2, 2, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1};
2447   int depthMax[neta] = {4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 3, 5, 6, 6, 6, 6, 6, 6, 6, 7, 7, 7, 3,
2448                         4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 3, 5, 6, 6, 6, 6, 6, 6, 6, 7, 7, 7, 3};
2449 
2450   // Read the correction factors
2451   std::map<unsigned int, std::pair<double, double> > cfactors;
2452   if (std::string(inFileCorr) != "") {
2453     std::ifstream fInput(inFileCorr);
2454     if (!fInput.good()) {
2455       std::cout << "Cannot open file " << inFileCorr << std::endl;
2456     } else {
2457       char buffer[1024];
2458       unsigned int all(0), good(0);
2459       while (fInput.getline(buffer, 1024)) {
2460         ++all;
2461         if (buffer[0] == '#')
2462           continue;  //ignore comment
2463         std::vector<std::string> items = splitString(std::string(buffer));
2464         if (items.size() != 5) {
2465           std::cout << "Ignore  line: " << buffer << std::endl;
2466         } else {
2467           ++good;
2468           std::istringstream converter(items[0].c_str());
2469           unsigned int det;
2470           converter >> std::hex >> det;
2471           double corrf = std::atof(items[3].c_str());
2472           double dcorr = std::atof(items[4].c_str());
2473           cfactors[det] = std::pair<double, double>(corrf, dcorr);
2474         }
2475       }
2476       fInput.close();
2477       std::cout << "Reads total of " << all << " and " << good << " good records" << std::endl;
2478     }
2479   }
2480 
2481   // Read in the depth dependent correction factors
2482   std::map<int, std::vector<double> > weights;
2483   if (std::string(inFileDepth) != "") {
2484     std::ifstream infile(inFileDepth);
2485     if (!infile.is_open()) {
2486       std::cout << "Cannot open duplicate file " << inFileDepth << std::endl;
2487     } else {
2488       unsigned int all(0), good(0);
2489       char buffer[1024];
2490       while (infile.getline(buffer, 1024)) {
2491         ++all;
2492         std::string bufferString(buffer);
2493         if (bufferString.substr(0, 1) == "#") {
2494           continue;  //ignore other comments
2495         } else {
2496           std::vector<std::string> items = splitString(bufferString);
2497           if (items.size() < 3) {
2498             std::cout << "Ignore  line: " << buffer << " Size " << items.size() << std::endl;
2499           } else {
2500             ++good;
2501             int ieta = std::atoi(items[0].c_str());
2502             std::vector<double> weight;
2503             for (unsigned int k = 1; k < items.size(); ++k) {
2504               double corrf = std::atof(items[k].c_str());
2505               weight.push_back(corrf);
2506             }
2507             weights[ieta] = weight;
2508           }
2509         }
2510       }
2511       infile.close();
2512       std::cout << "Reads total of " << all << " and " << good << " good records of depth dependent factors from "
2513                 << inFileDepth << std::endl;
2514     }
2515   }
2516 
2517   // Now combine the two information
2518   std::map<unsigned int, std::pair<double, double> > cfacFinal;
2519   std::ofstream outfile(outFileCorr);
2520   outfile << "#" << std::setprecision(4) << std::setw(10) << "detId" << std::setw(10) << "ieta" << std::setw(10)
2521           << "depth" << std::setw(15) << "corrFactor" << std::endl;
2522   for (int k = 0; k < neta; ++k) {
2523     int ieta = ietas[k];
2524     for (int depth = depthMin[k]; depth <= depthMax[k]; ++depth) {
2525       int subdet = ifHB(ieta, depth) ? 1 : 2;
2526       int d = truncateDepth(ieta, depth, truncateFlag);
2527       unsigned int key = repackId(subdet, ieta, 0, d);
2528       double c1(1), dc1(0), c2(1);
2529       if (cfactors.find(key) != cfactors.end()) {
2530         c1 = cfactors[key].first;
2531         dc1 = cfactors[key].second;
2532       }
2533       if (weights.find(ieta) != weights.end()) {
2534         if (static_cast<unsigned int>(depth) <= weights[ieta].size())
2535           c2 = (weights[ieta])[depth - 1];
2536       }
2537       double cf = c1 * c2;
2538       double dcf = dc1 * c2;
2539       if (std::abs(ieta) > etaMax) {
2540         int ieta0 = (ieta > 0) ? etaMax : -etaMax;
2541         key = repackId(subdet, ieta0, 0, depth);
2542         if (cfacFinal.find(key) != cfacFinal.end()) {
2543           cf = cfacFinal[key].first;
2544           dcf = cfacFinal[key].second;
2545         }
2546       }
2547       key = repackId(subdet, ieta, 0, depth);
2548       cfacFinal[key] = std::pair<double, double>(cf, dcf);
2549       if (std::abs(ieta) <= etaMax)
2550         outfile << std::setw(10) << std::hex << key << std::setw(10) << std::dec << ieta << std::setw(10) << depth
2551                 << std::setw(10) << cf << " " << std::setw(10) << dcf << std::endl;
2552     }
2553   }
2554   outfile.close();
2555 }