Back to home page

Project CMSSW displayed by LXR

 
 

    


File indexing completed on 2022-09-08 03:21:21

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