Back to home page

Project CMSSW displayed by LXR

 
 

    


File indexing completed on 2021-02-14 12:48:12

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