IsolatedParticles/macros/IsoTrkOfflineAnalyzer.C

0001 //////////////////////////////////////////////////////////
0002 // This class has been automatically generated on
0003 // Thu Jul 30 20:56:15 2015 by ROOT version 5.26/00
0004 // from TTree CalibTree/CalibTree
0005 // found on file: output_all.root
0006 //////////////////////////////////////////////////////////
0007
0008 #include <TROOT.h>
0009 #include <TChain.h>
0010 #include <TFile.h>
0011 #include <TF1.h>
0012 #include <TH1D.h>
0013 #include <TProfile.h>
0014 #include <TFitResult.h>
0015 #include <TFitResultPtr.h>
0016 #include <TStyle.h>
0017 #include <TCanvas.h>
0018 #include <TLegend.h>
0019 #include <TPaveStats.h>
0020 #include <TPaveText.h>
0021 #include <vector>
0022 #include <string>
0023 #include <iomanip>
0024 #include <iostream>
0025 #include <fstream>
0026
0027 class IsoTrkOfflineAnalyzer {
0028 public :
0029   TTree          *fChain;   //!pointer to the analyzed TTree or TChain
0030   Int_t           fCurrent; //!current Tree number in a TChain
0031
0032   // Declaration of leaf types
0033   Int_t                      t_Run;
0034   Int_t                      t_Event;
0035   Int_t                      t_ieta;
0036   Int_t                      t_goodPV;
0037   Double_t                   t_EventWeight;
0038   Double_t                   t_l1pt;
0039   Double_t                   t_l1eta;
0040   Double_t                   t_l1phi;
0041   Double_t                   t_l3pt;
0042   Double_t                   t_l3eta;
0043   Double_t                   t_l3phi;
0044   Double_t                   t_p;
0045   Double_t                   t_mindR1;
0046   Double_t                   t_mindR2;
0047   Double_t                   t_eMipDR;
0048   Double_t                   t_eHcal;
0049   Double_t                   t_hmaxNearP;
0050   Bool_t                     t_selectTk;
0051   Bool_t                     t_qltyFlag;
0052   Bool_t                     t_qltyMissFlag;
0053   Bool_t                     t_qltyPVFlag;
0054   std::vector<unsigned int> *t_DetIds;
0055   std::vector<double>       *t_HitEnergies;
0056   std::vector<bool>         *t_trgbits;
0057
0058   // List of branches
0059   TBranch                   *b_t_Run;           //!
0060   TBranch                   *b_t_Event;         //!
0061   TBranch                   *b_t_ieta;          //!
0062   TBranch                   *b_t_goodPV;        //!
0063   TBranch                   *b_t_EventWeight;   //!
0064   TBranch                   *b_t_l1pt;          //!
0065   TBranch                   *b_t_l1eta;         //!
0066   TBranch                   *b_t_l1phi;         //!
0067   TBranch                   *b_t_l3pt;          //!
0068   TBranch                   *b_t_l3eta;         //!
0069   TBranch                   *b_t_l3phi;         //!
0070   TBranch                   *b_t_p;             //!
0071   TBranch                   *b_t_mindR1;        //!
0072   TBranch                   *b_t_mindR2;        //!
0073   TBranch                   *b_t_eMipDR;        //!
0074   TBranch                   *b_t_eHcal;         //!
0075   TBranch                   *b_t_hmaxNearP;     //!
0076   TBranch                   *b_t_selectTk;      //!
0077   TBranch                   *b_t_qltyFlag;      //!
0078   TBranch                   *b_t_qltyMissFlag;  //!
0079   TBranch                   *b_t_qltyPVFlag;    //!
0080   TBranch                   *b_t_DetIds;        //!
0081   TBranch                   *b_t_HitEnergies;   //!
0082   TBranch                   *b_t_trgbits;       //!
0083   struct record {
0084     record() {
0085       serial_ = entry_ = run_ = event_ = ieta_ = p_ = 0;
0086     };
0087     record(int ser, Long64_t ent, int r, int ev, int ie, double p) :
0088       serial_(ser), entry_(ent), run_(r), event_(ev), ieta_(ie), p_(p) {};
0089     int      serial_;
0090     Long64_t entry_;
0091     int      run_, event_, ieta_;
0092     double   p_;
0093   };
0094
0095   IsoTrkOfflineAnalyzer(std::string fname, std::string dirname, std::string prefix="", int flag=0, bool datMC=true, std::string listname="runEventFile.txt");
0096   virtual ~IsoTrkOfflineAnalyzer();
0097   virtual Int_t              Cut(Long64_t entry);
0098   virtual Int_t              GetEntry(Long64_t entry);
0099   virtual Long64_t           LoadTree(Long64_t entry);
0100   virtual void               Init(TTree *tree);
0101   virtual void               Loop();
0102   virtual Bool_t             Notify();
0103   virtual void               Show(Long64_t entry = -1);
0104   std::vector<Long64_t>      findDuplicate (std::vector<IsoTrkOfflineAnalyzer::record>&);
0105   void                       PlotHist(int type, int num, bool save=false);
0106   template<class Hist> void  DrawHist(Hist*, TCanvas*);
0107   void                       SavePlot(std::string theName, bool append, bool all=false);
0108 private:
0109
0110   static const unsigned int npbin=10, kp50=7;
0111   std::string               fname_, dirnm_, prefix_, listName_;
0112   int                       flag_;
0113   bool                      dataMC_, plotStandard_;
0114   std::vector<double>       etas_, ps_, dl1_;
0115   std::vector<int>          nvx_;
0116   TH1D                     *h_p[5], *h_eta[5];
0117   std::vector<TH1D*>        h_eta0, h_eta1, h_eta2, h_eta3, h_eta4;
0118   std::vector<TH1D*>        h_dL1,  h_vtx, h_etaF[npbin];
0119   std::vector<TProfile*>    h_etaX[npbin];
0120   std::vector<TH1D*>        h_etaR[npbin], h_nvxR[npbin], h_dL1R[npbin];
0121 };
0122
0123 IsoTrkOfflineAnalyzer::IsoTrkOfflineAnalyzer(std::string fname, std::string dirnm, std::string prefix, int flag, bool datMC, std::string listname) : fname_(fname), dirnm_(dirnm), prefix_(prefix), listName_(listname), flag_(flag), dataMC_(datMC) {
0124   // if parameter tree is not specified (or zero), connect the file
0125   // used to generate this class and read the Tree
0126   plotStandard_    = (((flag_/10000)%10) == 0);
0127   TFile      *file = new TFile(fname.c_str());
0128   TDirectory *dir  = (TDirectory*)file->FindObjectAny(dirnm.c_str());
0129   std::cout << fname << " file " << file << " " << dirnm << " " << dir << std::endl;
0130   TTree      *tree = (TTree*)dir->Get("CalibTree");
0131   std::cout << "CalibTree " << tree << std::endl;
0132   Init(tree);
0133 }
0134
0135 IsoTrkOfflineAnalyzer::~IsoTrkOfflineAnalyzer() {
0136   if (!fChain) return;
0137   delete fChain->GetCurrentFile();
0138 }
0139
0140 Int_t IsoTrkOfflineAnalyzer::GetEntry(Long64_t entry) {
0141   // Read contents of entry.
0142   if (!fChain) return 0;
0143   return fChain->GetEntry(entry);
0144 }
0145
0146 Long64_t IsoTrkOfflineAnalyzer::LoadTree(Long64_t entry) {
0147   // Set the environment to read one entry
0148   if (!fChain) return -5;
0149   Long64_t centry = fChain->LoadTree(entry);
0150   if (centry < 0) return centry;
0151   if (!fChain->InheritsFrom(TChain::Class()))  return centry;
0152   TChain *chain = (TChain*)fChain;
0153   if (chain->GetTreeNumber() != fCurrent) {
0154     fCurrent = chain->GetTreeNumber();
0155     Notify();
0156   }
0157   return centry;
0158 }
0159
0160 void IsoTrkOfflineAnalyzer::Init(TTree *tree) {
0161   // The Init() function is called when the selector needs to initialize
0162   // a new tree or chain. Typically here the branch addresses and branch
0163   // pointers of the tree will be set.
0164   // It is normally not necessary to make changes to the generated
0165   // code, but the routine can be extended by the user if needed.
0166   // Init() will be called many times when running on PROOF
0167   // (once per file to be processed).
0168
0169   // Set object pointer
0170   t_DetIds      = 0;
0171   t_HitEnergies = 0;
0172   t_trgbits     = 0;
0173   // Set branch addresses and branch pointers
0174   if (!tree) return;
0175   fChain = tree;
0176   fCurrent = -1;
0177   fChain->SetMakeClass(1);
0178
0179   fChain->SetBranchAddress("t_Run", &t_Run, &b_t_Run);
0180   fChain->SetBranchAddress("t_Event", &t_Event, &b_t_Event);
0181   fChain->SetBranchAddress("t_ieta", &t_ieta, &b_t_ieta);
0182   fChain->SetBranchAddress("t_goodPV", &t_goodPV, &b_t_goodPV);
0183   fChain->SetBranchAddress("t_EventWeight", &t_EventWeight, &b_t_EventWeight);
0184   fChain->SetBranchAddress("t_l1pt", &t_l1pt, &b_t_l1pt);
0185   fChain->SetBranchAddress("t_l1eta", &t_l1eta, &b_t_l1eta);
0186   fChain->SetBranchAddress("t_l1phi", &t_l1phi, &b_t_l1phi);
0187   fChain->SetBranchAddress("t_l3pt", &t_l3pt, &b_t_l3pt);
0188   fChain->SetBranchAddress("t_l3eta", &t_l3eta, &b_t_l3eta);
0189   fChain->SetBranchAddress("t_l3phi", &t_l3phi, &b_t_l3phi);
0190   fChain->SetBranchAddress("t_p", &t_p, &b_t_p);
0191   fChain->SetBranchAddress("t_mindR1", &t_mindR1, &b_t_mindR1);
0192   fChain->SetBranchAddress("t_mindR2", &t_mindR2, &b_t_mindR2);
0193   fChain->SetBranchAddress("t_eMipDR", &t_eMipDR, &b_t_eMipDR);
0194   fChain->SetBranchAddress("t_eHcal", &t_eHcal, &b_t_eHcal);
0195   fChain->SetBranchAddress("t_hmaxNearP", &t_hmaxNearP, &b_t_hmaxNearP);
0196   fChain->SetBranchAddress("t_selectTk", &t_selectTk, &b_t_selectTk);
0197   fChain->SetBranchAddress("t_qltyFlag", &t_qltyFlag, &b_t_qltyFlag);
0198   fChain->SetBranchAddress("t_qltyMissFlag", &t_qltyMissFlag, &b_t_qltyMissFlag);
0199   fChain->SetBranchAddress("t_qltyPVFlag", &t_qltyPVFlag, &b_t_qltyPVFlag);
0200   fChain->SetBranchAddress("t_DetIds", &t_DetIds, &b_t_DetIds);
0201   fChain->SetBranchAddress("t_HitEnergies", &t_HitEnergies, &b_t_HitEnergies);
0202   fChain->SetBranchAddress("t_trgbits", &t_trgbits, &b_t_trgbits);
0203   Notify();
0204
0205   double xbins[9] = {-21.0, -16.0, -12.0, -6.0, 0.0, 6.0, 12.0, 16.0, 21.0};
0206   double xbina[43]= {-21.5,-20.5,-19.5,-18.5,-17.5,-16.5,-15.5,-14.5,-13.5,
0207              -12.5,-11.5,-10.5,-9.5,-8.5,-7.5,-6.5,-5.5,-4.5,-3.5,
0208              -2.5,-1.5,0.0,1.5,2.5,3.5,4.5,5.5,6.5,7.5,8.5,9.5,10.5,
0209              11.5,12.5,13.5,14.5,15.5,16.5,17.5,18.5,19.5,20.5,21.5};
0210   if (plotStandard_) {
0211     for (int i=0; i<9; ++i) etas_.push_back(xbins[i]);
0212   } else {
0213     for (int i=0; i<43; ++i) etas_.push_back(xbina[i]);
0214   }
0215   int ipbin[10] = {2, 4, 7, 10, 15, 20, 30, 40, 60, 100};
0216   for (int i=0; i<10; ++i) ps_.push_back((double)(ipbin[i]));
0217   int npvtx[6]  = {0, 7,10, 13, 16,100};
0218   for (int i=0; i<6; ++i)  nvx_.push_back(npvtx[i]);
0219   double dl1s[9]= {0, 0.10, 0.20, 0.50, 1.0, 2.0, 2.5, 3.0, 10.0};
0220
0221   char        name[20], title[100];
0222   std::string titl[5] = {"All tracks", "Good quality tracks", "Selected tracks",
0223              "Tracks with charge isolation", "Tracks MIP in ECAL"};
0224   for (int i=0; i<9; ++i)  dl1_.push_back(dl1s[i]);
0225   if (plotStandard_) {
0226     std::cout << "Book Histos for Standard\n";
0227     for (int k=0; k<5; ++k) {
0228       sprintf (name, "%sp%d", prefix_.c_str(), k);
0229       sprintf (title,"%s", titl[k].c_str());
0230       h_p[k] = new TH1D(name, title, 100, 10.0, 110.0);
0231       sprintf (name, "%seta%d", prefix_.c_str(), k);
0232       sprintf (title,"%s", titl[k].c_str());
0233       h_eta[k] = new TH1D(name, title, 60, -30.0, 30.0);
0234     }
0235     unsigned int kp = (ps_.size()-1);
0236     for (unsigned int k=0; k<kp; ++k) {
0237       sprintf (name, "%seta0%d", prefix_.c_str(), k);
0238       sprintf (title,"%s (p = %d:%d GeV)",titl[0].c_str(),ipbin[k],ipbin[k+1]);
0239       h_eta0.push_back(new TH1D(name, title, 60, -30.0, 30.0));
0240       sprintf (name, "%seta1%d", prefix_.c_str(), k);
0241       sprintf (title,"%s (p = %d:%d GeV)",titl[1].c_str(),ipbin[k],ipbin[k+1]);
0242       h_eta1.push_back(new TH1D(name, title, 60, -30.0, 30.0));
0243       sprintf (name, "%seta2%d", prefix_.c_str(), k);
0244       sprintf (title,"%s (p = %d:%d GeV)",titl[2].c_str(),ipbin[k],ipbin[k+1]);
0245       h_eta2.push_back(new TH1D(name, title, 60, -30.0, 30.0));
0246       sprintf (name, "%seta3%d", prefix_.c_str(), k);
0247       sprintf (title,"%s (p = %d:%d GeV)",titl[3].c_str(),ipbin[k],ipbin[k+1]);
0248       h_eta3.push_back(new TH1D(name, title, 60, -30.0, 30.0));
0249       sprintf (name, "%seta4%d", prefix_.c_str(), k);
0250       sprintf (title,"%s (p = %d:%d GeV)",titl[4].c_str(),ipbin[k],ipbin[k+1]);
0251       h_eta4.push_back(new TH1D(name, title, 60, -30.0, 30.0));
0252       sprintf (name, "%sdl1%d", prefix_.c_str(), k);
0253       sprintf (title,"Distance from L1 (p = %d:%d GeV)",ipbin[k],ipbin[k+1]);
0254       h_dL1.push_back(new TH1D(name, title, 160, 0.0, 8.0));
0255       sprintf (name, "%svtx%d", prefix_.c_str(), k);
0256       sprintf (title,"N_{Vertex} (p = %d:%d GeV)",ipbin[k],ipbin[k+1]);
0257       h_vtx.push_back(new TH1D(name, title, 100, 0.0, 100.0));
0258       for (unsigned int i=0; i<nvx_.size(); ++i) {
0259     sprintf (name, "%setaX%d%d", prefix_.c_str(), k, i);
0260     if (i == 0) {
0261       sprintf (title,"%s (p = %d:%d GeV all vertices)",titl[4].c_str(),ipbin[k],ipbin[k+1]);
0262     } else {
0263       sprintf (title,"%s (p = %d:%d GeV # Vtx %d:%d)",titl[4].c_str(),ipbin[k],ipbin[k+1],nvx_[i-1],nvx_[i]);
0264     }
0265     h_etaX[k].push_back(new TProfile(name, title, 8, xbins));
0266     unsigned int kk = h_etaX[k].size()-1;
0267     h_etaX[k][kk]->Sumw2();
0268     sprintf (name, "%snvxR%d%d", prefix_.c_str(), k, i);
0269     if (i == 0) {
0270       sprintf (title,"E/p for %s (p = %d:%d GeV all vertices)",titl[4].c_str(),ipbin[k],ipbin[k+1]);
0271     } else {
0272       sprintf (title,"E/p for %s (p = %d:%d GeV # Vtx %d:%d)",titl[4].c_str(),ipbin[k],ipbin[k+1],nvx_[i-1],nvx_[i]);
0273     }
0274     h_nvxR[k].push_back(new TH1D(name,title,100,0.,5.));
0275     kk = h_nvxR[k].size()-1;
0276     h_nvxR[k][kk]->Sumw2();
0277       }
0278       for (unsigned int j=0; j<etas_.size(); ++j) {
0279     sprintf (name, "%sratio%d%d", prefix_.c_str(), k, j);
0280     if (j == 0) {
0281       sprintf (title,"E/p for %s (p = %d:%d GeV)",titl[4].c_str(),ipbin[k],ipbin[k+1]);
0282     } else {
0283       sprintf (title,"E/p for %s (p = %d:%d GeV #eta %4.1f:%4.1f)",titl[4].c_str(),ipbin[k],ipbin[k+1],etas_[j-1],etas_[j]);
0284     }
0285     h_etaF[k].push_back(new TH1D(name,title,100,0.,5.));
0286     unsigned int kk = h_etaF[k].size()-1;
0287     h_etaF[k][kk]->Sumw2();
0288     sprintf (name, "%setaR%d%d", prefix_.c_str(), k, j);
0289     h_etaR[k].push_back(new TH1D(name,title,100,0.,5.));
0290     kk = h_etaR[k].size()-1;
0291     h_etaR[k][kk]->Sumw2();
0292       }
0293       for (unsigned int j=0; j<dl1_.size(); ++j) {
0294     sprintf (name, "%sdl1R%d%d", prefix_.c_str(), k, j);
0295     if (j == 0) {
0296       sprintf (title,"E/p for %s (p = %d:%d GeV All d_{L1})",titl[4].c_str(),ipbin[k],ipbin[k+1]);
0297     } else {
0298       sprintf (title,"E/p for %s (p = %d:%d GeV d_{L1} %4.2f:%4.2f)",titl[4].c_str(),ipbin[k],ipbin[k+1],dl1_[j-1],dl1_[j]);
0299     }
0300     h_dL1R[k].push_back(new TH1D(name,title,100,0.,5.));
0301     unsigned int kk = h_dL1R[k].size()-1;
0302     h_dL1R[k][kk]->Sumw2();
0303       }
0304     }
0305     for (unsigned int i=0; i<nvx_.size(); ++i) {
0306       sprintf (name, "%setaX%d%d", prefix_.c_str(), kp, i);
0307       if (i == 0) {
0308     sprintf (title,"%s (All Momentum all vertices)",titl[4].c_str());
0309       } else {
0310     sprintf (title,"%s (All Momentum # Vtx %d:%d)",titl[4].c_str(),nvx_[i-1],nvx_[i]);
0311       }
0312       h_etaX[npbin-1].push_back(new TProfile(name, title, 8, xbins));
0313       unsigned int kk = h_etaX[npbin-1].size()-1;
0314       h_etaX[npbin-1][kk]->Sumw2();
0315       sprintf (name, "%snvxR%d%d", prefix_.c_str(), kp, i);
0316       if (i == 0) {
0317     sprintf (title,"E/p for %s (All Momentum all vertices)",titl[4].c_str());
0318       } else {
0319     sprintf (title,"E/p for %s (All Momentum # Vtx %d:%d)",titl[4].c_str(),nvx_[i-1],nvx_[i]);
0320       }
0321       h_nvxR[npbin-1].push_back(new TH1D(name,title,200,0.,10.));
0322       kk = h_nvxR[npbin-1].size()-1;
0323       h_nvxR[npbin-1][kk]->Sumw2();
0324     }
0325     for (unsigned int j=0; j<etas_.size(); ++j) {
0326       sprintf (name, "%sratio%d%d", prefix_.c_str(), kp, j);
0327       if (j == 0) {
0328     sprintf (title,"E/p for %s (All momentum)",titl[4].c_str());
0329       } else {
0330     sprintf (title,"E/p for %s (All momentum #eta %4.1f:%4.1f)",titl[4].c_str(),etas_[j-1],etas_[j]);
0331       }
0332       h_etaF[npbin-1].push_back(new TH1D(name,title,200,0.,10.));
0333       unsigned int kk = h_etaF[npbin-1].size()-1;
0334       h_etaF[npbin-1][kk]->Sumw2();
0335       sprintf (name, "%setaR%d%d", prefix_.c_str(), kp, j);
0336       h_etaR[npbin-1].push_back(new TH1D(name,title,200,0.,10.));
0337       kk = h_etaR[npbin-1].size()-1;
0338       h_etaR[npbin-1][kk]->Sumw2();
0339     }
0340     for (unsigned int j=0; j<dl1_.size(); ++j) {
0341       sprintf (name, "%sdl1R%d%d", prefix_.c_str(), kp, j);
0342       if (j == 0) {
0343     sprintf (title,"E/p for %s (All momentum)",titl[4].c_str());
0344       } else {
0345     sprintf (title,"E/p for %s (All momentum d_{L1} %4.2f:%4.2f)",titl[4].c_str(),dl1_[j-1],dl1_[j]);
0346       }
0347       h_dL1R[npbin-1].push_back(new TH1D(name,title,200,0.,10.));
0348       unsigned int kk = h_dL1R[npbin-1].size()-1;
0349       h_dL1R[npbin-1][kk]->Sumw2();
0350     }
0351   } else {
0352     std::cout << "Book Histos for Non-Standard " << etas_.size() << ":" << kp50 << "\n";
0353     for (unsigned int j=0; j<etas_.size(); ++j) {
0354       sprintf (name, "%sratio%d%d", prefix_.c_str(), kp50, j);
0355       if (j == 0) {
0356     sprintf (title,"E/p for %s (p = %d:%d GeV)",titl[4].c_str(),ipbin[kp50],ipbin[kp50+1]);
0357       } else {
0358     sprintf (title,"E/p for %s (p = %d:%d GeV #eta %4.1f:%4.1f)",titl[4].c_str(),ipbin[kp50],ipbin[kp50+1],etas_[j-1],etas_[j]);
0359       }
0360       h_etaF[kp50].push_back(new TH1D(name,title,100,0.,5.));
0361       unsigned int kk = h_etaF[kp50].size()-1;
0362       h_etaF[kp50][kk]->Sumw2();
0363     }
0364   }
0365 }
0366
0367 Bool_t IsoTrkOfflineAnalyzer::Notify() {
0368   // The Notify() function is called when a new file is opened. This
0369   // can be either for a new TTree in a TChain or when when a new TTree
0370   // is started when using PROOF. It is normally not necessary to make changes
0371   // to the generated code, but the routine can be extended by the
0372   // user if needed. The return value is currently not used.
0373
0374   return kTRUE;
0375 }
0376
0377 void IsoTrkOfflineAnalyzer::Show(Long64_t entry) {
0378   // Print contents of entry.
0379   // If entry is not specified, print current entry
0380   if (!fChain) return;
0381   fChain->Show(entry);
0382 }
0383
0384 Int_t IsoTrkOfflineAnalyzer::Cut(Long64_t) {
0385   // This function may be called from Loop.
0386   // returns  1 if entry is accepted.
0387   // returns -1 otherwise.
0388   return 1;
0389 }
0390
0391 void IsoTrkOfflineAnalyzer::Loop() {
0392   //   In a ROOT session, you can do:
0393   //      Root > .L IsoTrkOfflineAnalyzer.C
0394   //      Root > IsoTrkOfflineAnalyzer t
0395   //      Root > t.GetEntry(12); // Fill t data members with entry number 12
0396   //      Root > t.Show();       // Show values of entry 12
0397   //      Root > t.Show(16);     // Read and show values of entry 16
0398   //      Root > t.Loop();       // Loop on all entries
0399   //
0400
0401   //   This is the loop skeleton where:
0402   //      jentry is the global entry number in the chain
0403   //      ientry is the entry number in the current Tree
0404   //   Note that the argument to GetEntry must be:
0405   //      jentry for TChain::GetEntry
0406   //      ientry for TTree::GetEntry and TBranch::GetEntry
0407   //
0408   //       To read only selected branches, Insert statements like:
0409   // METHOD1:
0410   //    fChain->SetBranchStatus("*",0);  // disable all branches
0411   //    fChain->SetBranchStatus("branchname",1);  // activate branchname
0412   // METHOD2: replace line
0413   //    fChain->GetEntry(jentry);       //read all branches
0414   //by  b_branchname->GetEntry(ientry); //read only this branch
0415   if (fChain == 0) return;
0416
0417   // Get list of run, event #'s if flag_>0
0418   std::vector<int> runs, events;
0419   runs.clear(); events.clear();
0420   std::vector<double> p_s;
0421   if ((flag_%10)>0) {
0422     ifstream infile(listName_.c_str());
0423     if (!infile.is_open()) {
0424       std::cout << "Cannot open " << listName_ << std::endl;
0425     } else {
0426       while (1) {
0427     int run, event;
0428     infile >> run >> event;
0429     if (!infile.good()) break;
0430     runs.push_back(run); events.push_back(event);
0431       }
0432       infile.close();
0433       std::cout << "Reads a list of " << runs.size() << " events from " << listName_ << std::endl;
0434       for (unsigned int k=0; k<runs.size(); ++k) std::cout << "[" << k << "] " << runs[k] << " " << events[k] << std::endl;
0435     }
0436   }
0437   std::ofstream file;
0438   if (((flag_/10)%10)>0) {
0439     file.open(listName_.c_str(), std::ofstream::out);
0440     std::cout << "Opens " << listName_ << " file in output mode\n";
0441   }
0442   std::ofstream fileout;
0443   if (((flag_/1000)%10)>0) {
0444     if (((flag_/1000)%10)==2) {
0445       fileout.open("events.txt", std::ofstream::out);
0446       std::cout << "Opens events.txt in output mode\n";
0447     } else {
0448       fileout.open("events.txt", std::ofstream::app);
0449       std::cout << "Opens events.txt in append mode\n";
0450     }
0451     fileout << "Input file: " << fname_ << " Directory: " << dirnm_
0452         << " Prefix: " << prefix_ << "\n";
0453   }
0454
0455   // Find list of duplicate events
0456   Long64_t nentries = fChain->GetEntriesFast();
0457   std::cout << "Total entries " << nentries << std::endl;
0458   Long64_t nbytes(0), nb(0);
0459   std::vector<Long64_t> entries;
0460   if (((flag_/100)%10)==0) {
0461     std::vector<IsoTrkOfflineAnalyzer::record> records;
0462     int                                        kount(0);
0463     for (Long64_t jentry=0; jentry<nentries;jentry++) {
0464       Long64_t ientry = LoadTree(jentry);
0465       if (ientry < 0) break;
0466       nb = fChain->GetEntry(jentry);   nbytes += nb;
0467       double cut = (t_p > 20) ? 2.0 : 0.0;
0468       double rcut= (t_p > 20) ? 0.25: 0.1;
0469       if (t_qltyFlag) {
0470     if (t_selectTk) {
0471       if (t_hmaxNearP < cut) {
0472         if (t_eMipDR < 1.0) {
0473           double rat = (t_p > 0) ? (t_eHcal/(t_p-t_eMipDR)) : 1.0;
0474           if (rat > rcut) {
0475         kount++;
0476         IsoTrkOfflineAnalyzer::record rec(kount,jentry,t_Run,t_Event,t_ieta,t_p);
0477         records.push_back(rec);
0478           }
0479         }
0480       }
0481     }
0482       }
0483     }
0484     entries = findDuplicate (records);
0485   } else if (((flag_/100)%10)==1) {
0486     char filename[100];
0487     sprintf (filename, "events_%s.txt", prefix_.c_str());
0488     ifstream infile(filename);
0489     if (!infile.is_open()) {
0490       std::cout << "Cannot open " << filename << std::endl;
0491     } else {
0492       while (1) {
0493     Long64_t jentry;
0494     infile >> jentry;
0495     if (!infile.good()) break;
0496     entries.push_back(jentry);
0497       }
0498       infile.close();
0499       std::cout << "Reads a list of " << entries.size() << " events from "
0500         << filename << std::endl;
0501     }
0502   }
0503
0504   nbytes = nb = 0;
0505   unsigned int  kount(0), duplicate(0), good(0);
0506   for (Long64_t jentry=0; jentry<nentries;jentry++) {
0507     Long64_t ientry = LoadTree(jentry);
0508     if (ientry < 0) break;
0509     nb = fChain->GetEntry(jentry);   nbytes += nb;
0510     bool select(true);
0511     if (runs.size() > 0) {
0512       select = false;
0513       for (unsigned int k=0; k<runs.size(); ++k) {
0514     if (t_Run == runs[k] && t_Event == events[k]) {
0515       select = true;
0516       break;
0517     }
0518       }
0519       if (!select) std::cout << "Unwanted event " << t_Run << " " << t_Event
0520                  << std::endl;
0521     }
0522     if ((entries.size() > 0) && (((flag_/100)%10)<=1)) {
0523       for (unsigned int k=0; k<entries.size(); ++k) {
0524     if (jentry == entries[k]) {
0525 //    std::cout << "Duplicate event " << t_Run << " " << t_Event << " " << t_p << std::endl;
0526       select = false;
0527       duplicate++;
0528       break;
0529     }
0530       }
0531     }
0532     if (((flag_/10)%10)>0) {
0533       if (t_p >= 20.0) {
0534     file << t_Run << " " << t_Event << "\n";
0535     ++kount;
0536       }
0537     }
0538     if (!select) continue;
0539
0540     // if (Cut(ientry) < 0) continue;
0541     int kp(-1), jp(-1);
0542     for (unsigned int k=1; k<ps_.size(); ++k ) {
0543       if (t_p >= ps_[k-1] && t_p < ps_[k]) {
0544     kp = k - 1; break;
0545       }
0546     }
0547     unsigned int kp1 = ps_.size() - 1;
0548     unsigned int kv1 = 0;
0549     unsigned int kv  = nvx_.size() - 1;
0550     for (unsigned int k=1; k<nvx_.size(); ++k ) {
0551       if (t_goodPV >= nvx_[k-1] && t_goodPV < nvx_[k]) {
0552     kv = k; break;
0553       }
0554     }
0555     unsigned int kd1 =  0;
0556     unsigned int kd  = dl1_.size() - 1;
0557     for (unsigned int k=1; k<dl1_.size(); ++k ) {
0558       if (t_mindR1 >= dl1_[k-1] && t_mindR1 < dl1_[k]) {
0559     kd = k; break;
0560       }
0561     }
0562     double eta = (t_ieta > 0) ? ((double)(t_ieta)-0.001) : ((double)(t_ieta)+0.001);
0563     for (unsigned int j=1; j<etas_.size(); ++j) {
0564       if (eta > etas_[j-1] && eta < etas_[j]) {
0565     jp = j; break;
0566       }
0567     }
0568 //    std::cout << "Bin " << kp << ":" << kp1 << ":" << kv << ":" << kv1 << ":" << kd << ":" << kd1 << ":" << jp << std::endl;
0569     if (plotStandard_) {
0570       h_p[0]->Fill(t_p,t_EventWeight);
0571       h_eta[0]->Fill(t_ieta,t_EventWeight);
0572       if (kp >= 0) h_eta0[kp]->Fill(t_ieta,t_EventWeight);
0573     }
0574     double rat = (t_p > 0) ? (t_eHcal/(t_p-t_eMipDR)) : 1.0;
0575     double cut = (t_p > 20) ? 2.0 : 0.0;
0576     double rcut= (t_p > 20) ? 0.25: 0.1;
0577 //    std::cout << "Entry " << jentry << " p|ratio " << t_p << "|" << rat << "|" << kp << "|" << kv << "|" << jp << " Cuts " << t_qltyFlag << "|" << t_selectTk << "|" << (t_hmaxNearP < cut) << "|" << (t_eMipDR < 1.0) << std::endl;
0578     if (t_qltyFlag) {
0579       if (plotStandard_) {
0580     h_p[1]->Fill(t_p,t_EventWeight);
0581     h_eta[1]->Fill(t_ieta,t_EventWeight);
0582     if (kp >= 0) h_eta1[kp]->Fill(t_ieta,t_EventWeight);
0583       }
0584       if (t_selectTk) {
0585     if (plotStandard_) {
0586       h_p[2]->Fill(t_p,t_EventWeight);
0587       h_eta[2]->Fill(t_ieta,t_EventWeight);
0588       if (kp >= 0) h_eta2[kp]->Fill(t_ieta,t_EventWeight);
0589     }
0590     if (t_hmaxNearP < cut) {
0591       if (plotStandard_) {
0592         h_p[3]->Fill(t_p,t_EventWeight);
0593         h_eta[3]->Fill(t_ieta,t_EventWeight);
0594         if (kp >= 0) h_eta3[kp]->Fill(t_ieta,t_EventWeight);
0595       }
0596       if (t_eMipDR < 1.0) {
0597         if (plotStandard_) {
0598           h_p[4]->Fill(t_p,t_EventWeight);
0599           h_eta[4]->Fill(t_ieta,t_EventWeight);
0600         }
0601         if (kp >= 0) {
0602           if (plotStandard_) {
0603         h_eta4[kp]->Fill(t_ieta,t_EventWeight);
0604         h_dL1[kp]->Fill(t_mindR1,t_EventWeight);
0605         h_vtx[kp]->Fill(t_goodPV,t_EventWeight);
0606           }
0607           if (rat > rcut) {
0608         if (plotStandard_) {
0609           h_etaX[kp][kv]->Fill(eta,rat,t_EventWeight);
0610           h_etaX[kp][kv1]->Fill(eta,rat,t_EventWeight);
0611           h_nvxR[kp][kv]->Fill(rat,t_EventWeight);
0612           h_nvxR[kp][kv1]->Fill(rat,t_EventWeight);
0613           h_dL1R[kp][kd]->Fill(rat,t_EventWeight);
0614           h_dL1R[kp][kd1]->Fill(rat,t_EventWeight);
0615           if (jp > 0) h_etaR[kp][jp]->Fill(rat,t_EventWeight);
0616           h_etaR[kp][0]->Fill(rat,t_EventWeight);
0617         }
0618         if ((!dataMC_) || (t_mindR1 > 0.5)) {
0619           if (plotStandard_) {
0620             if (jp > 0) h_etaF[kp][jp]->Fill(rat,t_EventWeight);
0621             h_etaF[kp][0]->Fill(rat,t_EventWeight);
0622           }
0623         }
0624         if ((!plotStandard_) && (kp == (int)(kp50))) {
0625 //        std::cout << "kp " << kp << h_etaF[kp].size() << std::endl;
0626           if ((!dataMC_) || (t_mindR1>0.5) || (((flag_/10000)%10)==2)) {
0627             if (jp > 0) h_etaF[kp][jp]->Fill(rat,t_EventWeight);
0628             h_etaF[kp][0]->Fill(rat,t_EventWeight);
0629           }
0630         }
0631           }
0632         }
0633         if (t_p > 20.0 && rat > rcut) {
0634           if (plotStandard_) {
0635         h_etaX[kp1][kv]->Fill(eta,rat,t_EventWeight);
0636         h_etaX[kp1][kv1]->Fill(eta,rat,t_EventWeight);
0637         h_nvxR[kp1][kv]->Fill(rat,t_EventWeight);
0638         h_nvxR[kp1][kv1]->Fill(rat,t_EventWeight);
0639         h_dL1R[kp1][kd]->Fill(rat,t_EventWeight);
0640         h_dL1R[kp1][kd1]->Fill(rat,t_EventWeight);
0641         if (jp > 0) h_etaR[kp1][jp]->Fill(rat,t_EventWeight);
0642         h_etaR[kp1][0]->Fill(rat,t_EventWeight);
0643         if ((!dataMC_) || (t_mindR1 > 0.5)) {
0644           if (jp > 0) h_etaF[kp1][jp]->Fill(rat,t_EventWeight);
0645           h_etaF[kp1][0]->Fill(rat,t_EventWeight);
0646         }
0647           }
0648           if (((flag_/1000)%10)!=0) {
0649         good++;
0650         fileout << good << " " << jentry << " " << t_Run  << " "
0651             << t_Event << " " << t_ieta << " " << t_p << std::endl;
0652           }
0653         }
0654       }
0655     }
0656       }
0657     }
0658   }
0659   if (((flag_/10)%10)>0) {
0660     file.close();
0661     std::cout << "Writes " << kount << " records in " << listName_ << std::endl;
0662   }
0663   if (((flag_/1000)%10)>0) {
0664     fileout.close();
0665     std::cout << "Writes " << good << " events in the file events.txt\n";
0666   }
0667   if (((flag_/100)%10)<=1) {
0668     std::cout << "Finds " << duplicate << " Duplicate events in this file\n";
0669   }
0670 }
0671
0672 std::vector<Long64_t>  IsoTrkOfflineAnalyzer::findDuplicate (std::vector<IsoTrkOfflineAnalyzer::record>& records){
0673   // First sort by run number
0674   for (int c = 0 ; c < ((int)(records.size())-1); c++) {
0675     for (int d = 0; d < ((int)(records.size())-c-1); d++) {
0676       if (records[d].run_ > records[d+1].run_) {
0677         record swap  = records[d];
0678         records[d]   = records[d+1];
0679         records[d+1] = swap;
0680       }
0681     }
0682   }
0683   // Then sort by event number
0684   for (int c = 0 ; c < ((int)(records.size())-1); c++) {
0685     for (int d = 0; d < ((int)(records.size())-c-1); d++) {
0686       if ((records[d].run_ == records[d+1].run_) &&
0687       (records[d].event_ > records[d+1].event_)) {
0688         record swap  = records[d];
0689         records[d]   = records[d+1];
0690         records[d+1] = swap;
0691       }
0692     }
0693   }
0694   // Finally by ieta
0695   for (int c = 0 ; c < ((int)(records.size())-1); c++) {
0696     for (int d = 0; d < ((int)(records.size())-c-1); d++) {
0697       if ((records[d].run_ == records[d+1].run_) &&
0698       (records[d].event_ == records[d+1].event_) &&
0699       (records[d].ieta_ > records[d+1].ieta_)) {
0700         record swap  = records[d];
0701         records[d]   = records[d+1];
0702         records[d+1] = swap;
0703       }
0704     }
0705   }
0706   // First sort by run number
0707   for (int c = 0 ; c < ((int)(records.size())-1); c++) {
0708     for (int d = 0; d < ((int)(records.size())-c-1); d++) {
0709       if (records[d].run_ > records[d+1].run_) {
0710         record swap  = records[d];
0711         records[d]   = records[d+1];
0712         records[d+1] = swap;
0713       }
0714     }
0715   }
0716   // Then sort by event number
0717   for (int c = 0 ; c < ((int)(records.size())-1); c++) {
0718     for (int d = 0; d < ((int)(records.size())-c-1); d++) {
0719       if ((records[d].run_ == records[d+1].run_) &&
0720       (records[d].event_ > records[d+1].event_)) {
0721         record swap  = records[d];
0722         records[d]   = records[d+1];
0723         records[d+1] = swap;
0724       }
0725     }
0726   }
0727   // Finally by ieta
0728   for (int c = 0 ; c < ((int)(records.size())-1); c++) {
0729     for (int d = 0; d < ((int)(records.size())-c-1); d++) {
0730       if ((records[d].run_ == records[d+1].run_) &&
0731       (records[d].event_ == records[d+1].event_) &&
0732       (records[d].ieta_ > records[d+1].ieta_)) {
0733         record swap  = records[d];
0734         records[d]   = records[d+1];
0735         records[d+1] = swap;
0736       }
0737     }
0738   }
0739   // Find duplicate events
0740   std::vector<Long64_t> entries;
0741   for (unsigned int k=1; k<records.size(); ++k) {
0742     if ((records[k].run_ == records[k-1].run_) &&
0743     (records[k].event_ == records[k-1].event_) &&
0744     (records[k].ieta_ == records[k-1].ieta_) &&
0745     (fabs(records[k].p_-records[k-1].p_) < 0.0001)) {
0746       // This is a duplicate event
0747       entries.push_back(records[k].entry_);
0748     }
0749   }
0750   return entries;
0751 }
0752
0753 void IsoTrkOfflineAnalyzer::PlotHist(int itype, int inum, bool save) {
0754   gStyle->SetCanvasBorderMode(0); gStyle->SetCanvasColor(kWhite);
0755   gStyle->SetPadColor(kWhite);    gStyle->SetFillColor(kWhite);
0756   gStyle->SetOptStat(111110);     gStyle->SetOptFit(1);
0757   char name[100];
0758   int itmin = (itype >=0 && itype < 14) ? itype : 0;
0759   int itmax = (itype >=0 && itype < 14) ? itype : 13;
0760   std::string types[14] = {"p","i#eta","i#eta","i#eta","i#eta","i#eta",
0761                "i#eta","i#eta","E_{HCAL}/(p-E_{ECAL})",
0762                "E_{HCAL}/(p-E_{ECAL})","E_{HCAL}/(p-E_{ECAL})",
0763                "E_{HCAL}/(p-E_{ECAL})","dR_{L1}","Vertex"};
0764   int nmax[14] = {5, 5, npbin-1, npbin-1, npbin-1, npbin-1, npbin-1,
0765           npbin, npbin, npbin, npbin, npbin, npbin-1,npbin-1};
0766   for (int type=itmin; type<=itmax; ++type) {
0767     int inmin = (inum >=0 && inum < nmax[type]) ? inum : 0;
0768     int inmax = (inum >=0 && inum < nmax[type]) ? inum : nmax[type]-1;
0769     int kmax  = 1;
0770     if      (type == 8)  kmax = (int)(etas_.size());
0771     else if (type == 7)  kmax = (int)(nvx_.size());
0772     else if (type == 9)  kmax = (int)(etas_.size());
0773     else if (type == 10) kmax = (int)(nvx_.size());
0774     else if (type == 11) kmax = (int)(dl1_.size());
0775     for (int num=inmin; num<=inmax; ++num) {
0776       for (int k=0; k<kmax; ++k) {
0777     sprintf (name, "c_%d%d%d", type, num, k);
0778     TCanvas *pad = new TCanvas(name, name, 700, 500);
0779     pad->SetRightMargin(0.10);
0780     pad->SetTopMargin(0.10);
0781     sprintf (name, "%s", types[type].c_str());
0782     if (type != 7) {
0783       TH1D* hist(0);
0784       if      (type == 0)  hist = (TH1D*)(h_p[num]->Clone());
0785       else if (type == 1)  hist = (TH1D*)(h_eta[num]->Clone());
0786       else if (type == 2)  hist = (TH1D*)(h_eta0[num]->Clone());
0787       else if (type == 3)  hist = (TH1D*)(h_eta1[num]->Clone());
0788       else if (type == 4)  hist = (TH1D*)(h_eta2[num]->Clone());
0789       else if (type == 5)  hist = (TH1D*)(h_eta3[num]->Clone());
0790       else if (type == 6)  hist = (TH1D*)(h_eta4[num]->Clone());
0791       else if (type == 8)  hist = (TH1D*)(h_etaR[num][k]->Clone());
0792       else if (type == 9)  hist = (TH1D*)(h_etaF[num][k]->Clone());
0793       else if (type == 10) hist = (TH1D*)(h_nvxR[num][k]->Clone());
0794       else if (type == 11) hist = (TH1D*)(h_dL1R[num][k]->Clone());
0795       else if (type == 12) hist = (TH1D*)(h_dL1[num]->Clone());
0796       else                 hist = (TH1D*)(h_vtx[num]->Clone());
0797       hist->GetXaxis()->SetTitle(name);
0798       hist->GetYaxis()->SetTitle("Tracks");
0799       DrawHist(hist,pad);
0800       if (save) {
0801         sprintf (name, "c_%s%d%d%d.gif", prefix_.c_str(), type,num,k);
0802         pad->Print(name);
0803       }
0804     } else {
0805       TProfile* hist = (TProfile*)(h_etaX[num][k]->Clone());
0806       hist->GetXaxis()->SetTitle(name);
0807       hist->GetYaxis()->SetTitle("<E_{HCAL}/(p-E_{ECAL})>");
0808       hist->GetYaxis()->SetRangeUser(0.4,1.6);
0809       hist->Fit("pol0","q");
0810       DrawHist(hist,pad);
0811       if (save) {
0812         sprintf (name, "c_%s%d%d%d.gif", prefix_.c_str(), type,num,k);
0813         pad->Print(name);
0814       }
0815     }
0816       }
0817     }
0818   }
0819 }
0820
0821 template<class Hist> void IsoTrkOfflineAnalyzer::DrawHist(Hist* hist, TCanvas* pad) {
0822   hist->GetYaxis()->SetLabelOffset(0.005);
0823   hist->GetYaxis()->SetTitleOffset(1.20);
0824   hist->Draw();
0825   pad->Update();
0826   TPaveStats* st1 = (TPaveStats*)hist->GetListOfFunctions()->FindObject("stats");
0827   if (st1 != NULL) {
0828     st1->SetY1NDC(0.70); st1->SetY2NDC(0.90);
0829     st1->SetX1NDC(0.55); st1->SetX2NDC(0.90);
0830   }
0831   pad->Modified();
0832   pad->Update();
0833 }
0834
0835 void IsoTrkOfflineAnalyzer::SavePlot(std::string theName, bool append, bool all) {
0836
0837   TFile* theFile(0);
0838   if (append) {
0839     theFile = new TFile(theName.c_str(), "UPDATE");
0840   } else {
0841     theFile = new TFile(theName.c_str(), "RECREATE");
0842   }
0843
0844   theFile->cd();
0845   for (unsigned int k=0; k<ps_.size(); ++k) {
0846     if (plotStandard_) {
0847       for (unsigned int i=0; i<nvx_.size(); ++i) {
0848     if (h_etaX[k][i] != 0) {
0849       TProfile* hnew = (TProfile*)h_etaX[k][i]->Clone();
0850       hnew->Write();
0851     }
0852     if (h_nvxR[k].size() > i && h_nvxR[k][i] != 0) {
0853       TH1D* hist = (TH1D*)h_nvxR[k][i]->Clone();
0854       hist->Write();
0855     }
0856       }
0857     }
0858     for (unsigned int j=0; j<etas_.size(); ++j) {
0859       if (plotStandard_ && h_etaR[k][j] != 0) {
0860     TH1D* hist = (TH1D*)h_etaR[k][j]->Clone();
0861     hist->Write();
0862       }
0863       if (h_etaF[k].size() > j && h_etaF[k][j] != 0) {
0864     TH1D* hist = (TH1D*)h_etaF[k][j]->Clone();
0865     hist->Write();
0866       }
0867     }
0868     if (plotStandard_) {
0869       for (unsigned int j=0; j<dl1_.size(); ++j) {
0870     if (h_dL1R[k][j] != 0) {
0871       TH1D* hist = (TH1D*)h_dL1R[k][j]->Clone();
0872       hist->Write();
0873     }
0874       }
0875     }
0876     if (all && plotStandard_ && ((k+1) < ps_.size())) {
0877       if (h_eta0[k] != 0) {TH1D* h1 = (TH1D*)h_eta0[k]->Clone(); h1->Write();}
0878       if (h_eta1[k] != 0) {TH1D* h2 = (TH1D*)h_eta1[k]->Clone(); h2->Write();}
0879       if (h_eta2[k] != 0) {TH1D* h3 = (TH1D*)h_eta2[k]->Clone(); h3->Write();}
0880       if (h_eta3[k] != 0) {TH1D* h4 = (TH1D*)h_eta3[k]->Clone(); h4->Write();}
0881       if (h_eta4[k] != 0) {TH1D* h5 = (TH1D*)h_eta4[k]->Clone(); h5->Write();}
0882       if (h_dL1[k] != 0)  {TH1D* h6 = (TH1D*)h_dL1[k]->Clone();  h6->Write();}
0883       if (h_vtx[k] != 0)  {TH1D* h7 = (TH1D*)h_vtx[k]->Clone();  h7->Write();}
0884     }
0885   }
0886   std::cout << "All done\n";
0887   theFile->Close();
0888 }
0889
0890 class GetEntries {
0891 public :
0892   TTree                     *fChain;   //!pointer to the analyzed TTree/TChain
0893   Int_t                      fCurrent; //!current Tree number in a TChain
0894
0895   // Declaration of leaf types
0896   Int_t                      t_Tracks;
0897   Int_t                      t_TracksProp;
0898   Int_t                      t_TracksSaved;
0899   std::vector<int>          *t_ietaAll;
0900   std::vector<int>          *t_ietaGood;
0901
0902   // List of branches
0903   TBranch                   *b_t_Tracks;        //!
0904   TBranch                   *b_t_TracksProp;    //!
0905   TBranch                   *b_t_TracksSaved;   //!
0906   TBranch                   *b_t_ietaAll;       //!
0907   TBranch                   *b_t_ietaGood;      //!
0908
0909   GetEntries(std::string fname, std::string dirname, bool ifOld=true);
0910   virtual ~GetEntries();
0911   virtual Int_t    Cut(Long64_t entry);
0912   virtual Int_t    GetEntry(Long64_t entry);
0913   virtual Long64_t LoadTree(Long64_t entry);
0914   virtual void     Init(TTree *tree);
0915   virtual void     Loop();
0916   virtual Bool_t   Notify();
0917   virtual void     Show(Long64_t entry = -1);
0918
0919 private:
0920   bool             ifOld_;
0921   TH1I            *h_tk[3], *h_eta[2];
0922   TH1D            *h_eff;
0923 };
0924
0925 GetEntries::GetEntries(std::string fname, std::string dirnm, bool ifOld) : ifOld_(ifOld) {
0926
0927   TFile      *file = new TFile(fname.c_str());
0928   TDirectory *dir  = (TDirectory*)file->FindObjectAny(dirnm.c_str());
0929   std::cout << fname << " file " << file << " " << dirnm << " " << dir << std::endl;
0930   TTree      *tree = (TTree*)dir->Get("EventInfo");
0931   std::cout << "CalibTree " << tree << std::endl;
0932   Init(tree);
0933 }
0934
0935 GetEntries::~GetEntries() {
0936   if (!fChain) return;
0937   delete fChain->GetCurrentFile();
0938 }
0939
0940 Int_t GetEntries::GetEntry(Long64_t entry) {
0941   // Read contents of entry.
0942   if (!fChain) return 0;
0943   return fChain->GetEntry(entry);
0944 }
0945
0946 Long64_t GetEntries::LoadTree(Long64_t entry) {
0947   // Set the environment to read one entry
0948   if (!fChain) return -5;
0949   Long64_t centry = fChain->LoadTree(entry);
0950   if (centry < 0) return centry;
0951   if (!fChain->InheritsFrom(TChain::Class()))  return centry;
0952   TChain *chain = (TChain*)fChain;
0953   if (chain->GetTreeNumber() != fCurrent) {
0954     fCurrent = chain->GetTreeNumber();
0955     Notify();
0956   }
0957   return centry;
0958 }
0959
0960 void GetEntries::Init(TTree *tree) {
0961   // The Init() function is called when the selector needs to initialize
0962   // a new tree or chain. Typically here the branch addresses and branch
0963   // pointers of the tree will be set.
0964   // It is normally not necessary to make changes to the generated
0965   // code, but the routine can be extended by the user if needed.
0966   // Init() will be called many times when running on PROOF
0967   // (once per file to be processed).
0968
0969   // Set branch addresses and branch pointers
0970   if (!tree) return;
0971   fChain = tree;
0972   fCurrent = -1;
0973   fChain->SetMakeClass(1);
0974   fChain->SetBranchAddress("t_Tracks",      &t_Tracks,      &b_t_Tracks);
0975   fChain->SetBranchAddress("t_TracksProp",  &t_TracksProp,  &b_t_TracksProp);
0976   fChain->SetBranchAddress("t_TracksSaved", &t_TracksSaved, &b_t_TracksSaved);
0977   if (!ifOld_) {
0978     fChain->SetBranchAddress("t_ietaAll",     &t_ietaAll,     &b_t_ietaAll);
0979     fChain->SetBranchAddress("t_ietaGood",    &t_ietaGood,    &b_t_ietaGood);
0980   }
0981   Notify();
0982
0983   h_tk[0] = new TH1I("Track0", "# of tracks produced",      2000, 0, 2000);
0984   h_tk[1] = new TH1I("Track1", "# of tracks propagated",    2000, 0, 2000);
0985   h_tk[2] = new TH1I("Track2", "# of tracks saved in tree", 2000, 0, 2000);
0986   h_eta[0] = new TH1I("Eta0", "i#eta (All Tracks)",           60, -30, 30);
0987   h_eta[1] = new TH1I("Eta1", "i#eta (Good Tracks)",          60, -30, 30);
0988   h_eff    = new TH1D("Eta2", "i#eta (Selection Efficiency)", 60, -30, 30);
0989 }
0990
0991 Bool_t GetEntries::Notify() {
0992   // The Notify() function is called when a new file is opened. This
0993   // can be either for a new TTree in a TChain or when when a new TTree
0994   // is started when using PROOF. It is normally not necessary to make changes
0995   // to the generated code, but the routine can be extended by the
0996   // user if needed. The return value is currently not used.
0997
0998   return kTRUE;
0999 }
1000
1001 void GetEntries::Show(Long64_t entry) {
1002   // Print contents of entry.
1003   // If entry is not specified, print current entry
1004   if (!fChain) return;
1005   fChain->Show(entry);
1006 }
1007
1008 Int_t GetEntries::Cut(Long64_t) {
1009   // This function may be called from Loop.
1010   // returns  1 if entry is accepted.
1011   // returns -1 otherwise.
1012   return 1;
1013 }
1014
1015 void GetEntries::Loop() {
1016   //   In a ROOT session, you can do:
1017   //      Root > .L IsoTrkOfflineAnalyzer.C+g
1018   //      Root > GetEntries t
1019   //      Root > t.GetEntry(12); // Fill t data members with entry number 12
1020   //      Root > t.Show();       // Show values of entry 12
1021   //      Root > t.Show(16);     // Read and show values of entry 16
1022   //      Root > t.Loop();       // Loop on all entries
1023   //
1024
1025   //     This is the loop skeleton where:
1026   //    jentry is the global entry number in the chain
1027   //    ientry is the entry number in the current Tree
1028   //  Note that the argument to GetEntry must be:
1029   //    jentry for TChain::GetEntry
1030   //    ientry for TTree::GetEntry and TBranch::GetEntry
1031   //
1032   //       To read only selected branches, Insert statements like:
1033   // METHOD1:
1034   //    fChain->SetBranchStatus("*",0);  // disable all branches
1035   //    fChain->SetBranchStatus("branchname",1);  // activate branchname
1036   // METHOD2: replace line
1037   //    fChain->GetEntry(jentry);       //read all branches
1038   //by  b_branchname->GetEntry(ientry); //read only this branch
1039   if (fChain == 0) return;
1040
1041   Long64_t nentries = fChain->GetEntriesFast();
1042
1043   Long64_t nbytes = 0, nb = 0;
1044   for (Long64_t jentry=0; jentry<nentries;jentry++) {
1045     Long64_t ientry = LoadTree(jentry);
1046     if (ientry < 0) break;
1047     nb = fChain->GetEntry(jentry);   nbytes += nb;
1048     // if (Cut(ientry) < 0) continue;
1049     h_tk[0]->Fill(t_Tracks);
1050     h_tk[1]->Fill(t_TracksProp);
1051     h_tk[2]->Fill(t_TracksSaved);
1052     if (!ifOld_) {
1053       for (unsigned int k=0; k<t_ietaAll->size(); ++k)
1054     h_eta[0]->Fill((*t_ietaAll)[k]);
1055       for (unsigned int k=0; k<t_ietaGood->size(); ++k)
1056     h_eta[1]->Fill((*t_ietaGood)[k]);
1057     }
1058   }
1059   double ymaxk(0);
1060   if (!ifOld_) {
1061     for (int i=1; i<=h_eff->GetNbinsX(); ++i) {
1062       double rat(0), drat(0);
1063       if (h_eta[0]->GetBinContent(i) > ymaxk) ymaxk = h_eta[0]->GetBinContent(i);
1064       if (h_eta[1]->GetBinContent(i) > 0) {
1065     rat = h_eta[1]->GetBinContent(i)/h_eta[0]->GetBinContent(i);
1066     drat= rat*std::sqrt(pow((h_eta[1]->GetBinError(i)/h_eta[1]->GetBinContent(i)),2) +
1067                 pow((h_eta[0]->GetBinError(i)/h_eta[0]->GetBinContent(i)),2));
1068       }
1069       h_eff->SetBinContent(i,rat);
1070       h_eff->SetBinError(i,drat);
1071     }
1072   }
1073   gStyle->SetCanvasBorderMode(0); gStyle->SetCanvasColor(kWhite);
1074   gStyle->SetPadColor(kWhite);    gStyle->SetFillColor(kWhite);
1075   gStyle->SetOptStat(1110);       gStyle->SetOptTitle(0);
1076   int color[3] = {kBlack, kRed, kBlue};
1077   int lines[3] = {1, 2, 3};
1078   TCanvas *pad1 = new TCanvas("c_track", "c_track", 500, 500);
1079   pad1->SetRightMargin(0.10);
1080   pad1->SetTopMargin(0.10);
1081   pad1->SetFillColor(kWhite);
1082   std::string titl1[3] = {"Reconstructed", "Propagated", "Saved"};
1083   TLegend  *legend1 = new TLegend(0.55, 0.54, 0.90, 0.63);
1084   legend1->SetFillColor(kWhite);
1085   double ymax(0), xmax(0);
1086   for (int k=0; k<3; ++k) {
1087     int total(0), totaltk(0);
1088     for (int i=1; i<=h_tk[k]->GetNbinsX(); ++i) {
1089       if (ymax < h_tk[k]->GetBinContent(i)) ymax = h_tk[k]->GetBinContent(i);
1090       if (i > 1) total += (int)(h_tk[k]->GetBinContent(i));
1091       totaltk += (int)(h_tk[k]->GetBinContent(i))*(i-1);
1092       if (h_tk[k]->GetBinContent(i) > 0) {
1093     if (xmax < h_tk[k]->GetBinLowEdge(i)+h_tk[k]->GetBinWidth(i))
1094       xmax = h_tk[k]->GetBinLowEdge(i)+h_tk[k]->GetBinWidth(i);
1095       }
1096     }
1097     h_tk[k]->SetLineColor(color[k]);
1098     h_tk[k]->SetMarkerColor(color[k]);
1099     h_tk[k]->SetLineStyle(lines[k]);
1100     std::cout << h_tk[k]->GetTitle() << " Entries " << h_tk[k]->GetEntries()
1101           << " Events " << total << " Tracks " << totaltk << std::endl;
1102     legend1->AddEntry(h_tk[k],titl1[k].c_str(),"l");
1103   }
1104   int i1 = (int)(0.1*xmax) + 1;
1105   xmax   = 10.0*i1;
1106   int i2 = (int)(0.01*ymax) + 1;
1107
1108   ymax   = 100.0*i2;
1109   for (int k=0; k<3; ++k) {
1110     h_tk[k]->GetXaxis()->SetRangeUser(0,  xmax);
1111     h_tk[k]->GetYaxis()->SetRangeUser(0.1,ymax);
1112     h_tk[k]->GetXaxis()->SetTitle("# Tracks");
1113     h_tk[k]->GetYaxis()->SetTitle("Events");
1114     h_tk[k]->GetYaxis()->SetLabelOffset(0.005);
1115     h_tk[k]->GetYaxis()->SetTitleOffset(1.20);
1116     if (k == 0) h_tk[k]->Draw("hist");
1117     else        h_tk[k]->Draw("hist sames");
1118   }
1119   pad1->Update();
1120   pad1->SetLogy();
1121   ymax = 0.90;
1122   for (int k=0; k<3; ++k) {
1123     TPaveStats* st1 = (TPaveStats*)h_tk[k]->GetListOfFunctions()->FindObject("stats");
1124     if (st1 != NULL) {
1125       st1->SetLineColor(color[k]);
1126       st1->SetTextColor(color[k]);
1127       st1->SetY1NDC(ymax-0.09); st1->SetY2NDC(ymax);
1128       st1->SetX1NDC(0.55);      st1->SetX2NDC(0.90);
1129       ymax -= 0.09;
1130     }
1131   }
1132   pad1->Modified();
1133   pad1->Update();
1134   legend1->Draw("same");
1135   pad1->Update();
1136
1137   if (!ifOld_) {
1138     TCanvas *pad2 = new TCanvas("c_ieta", "c_ieta", 500, 500);
1139     pad2->SetRightMargin(0.10);
1140     pad2->SetTopMargin(0.10);
1141     pad2->SetFillColor(kWhite);
1142     pad2->SetLogy();
1143     std::string titl2[2] = {"All Tracks", "Selected Tracks"};
1144     TLegend  *legend2 = new TLegend(0.55, 0.28, 0.90, 0.35);
1145     legend2->SetFillColor(kWhite);
1146     i2    = (int)(0.001*ymaxk) + 1;
1147     ymax  = 1000.0*i2;
1148     for (int k=0; k<2; ++k) {
1149       h_eta[k]->GetYaxis()->SetRangeUser(1,ymax);
1150       h_eta[k]->SetLineColor(color[k]);
1151       h_eta[k]->SetMarkerColor(color[k]);
1152       h_eta[k]->SetLineStyle(lines[k]);
1153       h_eta[k]->GetXaxis()->SetTitle("i#eta");
1154       h_eta[k]->GetYaxis()->SetTitle("Tracks");
1155       h_eta[k]->GetYaxis()->SetLabelOffset(0.005);
1156       h_eta[k]->GetYaxis()->SetTitleOffset(1.20);
1157       legend2->AddEntry(h_eta[k],titl2[k].c_str(),"l");
1158       if (k == 0) h_eta[k]->Draw("hist");
1159       else        h_eta[k]->Draw("hist sames");
1160     }
1161     pad2->Update();
1162     double ymin = 0.10;
1163     for (int k=0; k<2; ++k) {
1164       TPaveStats* st1 = (TPaveStats*)h_eta[k]->GetListOfFunctions()->FindObject("stats");
1165       if (st1 != NULL) {
1166     st1->SetLineColor(color[k]);
1167     st1->SetTextColor(color[k]);
1168     st1->SetY1NDC(ymin); st1->SetY2NDC(ymin+0.09);
1169     st1->SetX1NDC(0.55); st1->SetX2NDC(0.90);
1170     ymin += 0.09;
1171       }
1172     }
1173     pad2->Modified();
1174     pad2->Update();
1175     legend2->Draw("same");
1176     pad2->Update();
1177
1178     TCanvas *pad3 = new TCanvas("c_effi", "c_effi", 500, 500);
1179     pad3->SetRightMargin(0.10);
1180     pad3->SetTopMargin(0.10);
1181     pad3->SetFillColor(kWhite);
1182     pad3->SetLogy();
1183     h_eff->SetStats(0);
1184     h_eff->SetMarkerStyle(20);
1185     h_eff->GetXaxis()->SetTitle("i#eta");
1186     h_eff->GetYaxis()->SetTitle("Efficiency");
1187     h_eff->GetYaxis()->SetLabelOffset(0.005);
1188     h_eff->GetYaxis()->SetTitleOffset(1.20);
1189     h_eff->Draw();
1190     pad3->Modified();
1191     pad3->Update();
1192   }
1193 }
1194
1195 template<class Hist> void
1196 DrawHist(std::vector<Hist*> hist, std::vector<int> comb,
1197      int m, bool typex, bool save) {
1198   std::string titles[27] = {"Data with Method 0", "Data with Method 2",
1199                 "Data with Method 3", "Data with Method 0",
1200                 "Data with Method 2", "Data with Method 3",
1201                 "Data with Method 0", "Data with Method 2",
1202                 "Data with Method 3", "Data with Method 0",
1203                 "Data with Method 2", "Data with Method 3",
1204                 "MC (#pi) Method 0",  "MC (#pi) Method 2",
1205                             "MC (#pi) Method 3",  "MC (#pi) No PU",
1206                 "MC (QCD) Method 0",  "MC (QCD) Method 0",
1207                 "MC (QCD) Method 0",  "MC (QCD) Method 0",
1208                 "p = 40:60 GeV Tracks (Method 0)",
1209                 "p = 40:60 GeV Tracks (Method 2)",
1210                 "p = 40:60 GeV Tracks (Method 3)",
1211                 "p = 40:60 GeV Tracks (Method 0)",
1212                 "p = 20:30 GeV Tracks (Method 0)",
1213                 "p = 40:60 GeV Tracks", "p = 20:30 GeV Tracks"};
1214   int colors[6]={1,2,4,7,6,9};
1215   int mtype[6]={20,21,22,23,24,33};
1216   std::string dtitl[17] = {"Data (25 ns) Method 0",  "Data (25 ns) Method 2",
1217                "Data (25 ns) Method 3",  "Data (50 ns) Method 0",
1218                "Data (50 ns) Method 2",  "Data (50 ns) Method 3",
1219                "Data (B+C) Method 0",    "Data (B+C) Method 2",
1220                "Data (B+C) Method 3",    "#pi MC (No PU) Method 0",
1221                "#pi MC (No PU) Method 2","#pi MC (No PU) Method 3",
1222                "#pi MC (25 ns) Method 0","#pi MC (25 ns) Method 2",
1223                "#pi MC (25 ns) Method 3","QCD MC (25 ns) Method 0",
1224                "QCD MC (50 ns) Method 0"};
1225   std::string stitl[30] = {"p 2:4 GeV all PV","p 4:7 GeV all PV",
1226                "p 7:10 GeV all PV","p 10:15 GeV all PV",
1227                "p 15:20 GeV all PV","p 20:30 GeV all PV",
1228                "p 30:40 GeV all PV","p 40:60 GeV all PV",
1229                "p 60:100 GeV all PV","all p all PV",
1230                "p 2:4 GeV PV 0:12","p 4:7 GeV PV 0:12",
1231                "p 7:10 GeV PV 0:12","p 10:15 GeV PV 0:12",
1232                "p 15:20 GeV PV 0:12","p 20:30 GeV PV 0:12",
1233                "p 30:40 GeV PV 0:12","p 40:60 GeV PV 0:12",
1234                "p 60:100 GeV PV 0:12","all p PV 0:12",
1235                "p 2:4 GeV PV > 12","p 4:7 GeV PV > 12",
1236                "p 7:10 GeV PV > 12","p 10:15 GeV PV > 12",
1237                "p 15:20 GeV PV > 12","p 20:30 GeV PV > 12",
1238                "p 30:40 GeV PV > 12","p 40:60 GeV PV > 12",
1239                "p 60:100 GeV PV > 12","all p PV > 12"};
1240
1241   char name[20];
1242   std::vector<double> resultv, resulte;
1243   for (unsigned int j=0; j<hist.size(); ++j) {
1244     hist[j]->SetTitle(titles[m].c_str());
1245     hist[j]->SetLineColor(colors[j]);
1246     hist[j]->SetMarkerColor(colors[j]);
1247     hist[j]->SetMarkerStyle(mtype[j]);
1248     hist[j]->GetYaxis()->SetRangeUser(0.4,1.6);
1249     hist[j]->GetYaxis()->SetTitle("<E_{HCAL}/(p-E_{ECAL})>");
1250     hist[j]->GetXaxis()->SetTitle("i#eta");
1251     /*
1252     if (!typex) {
1253       Hist* hist1 = (Hist*)hist[j]->Clone();
1254       TFitResultPtr Fit = hist1->Fit("pol0","+WWQRD","");
1255       resultv.push_back(Fit->Value(1));
1256       std::cout << "Fit " << Fit->Value(1) << std::endl;
1257     }
1258     */
1259   }
1260   if (typex) sprintf (name, "c_respX%d", m);
1261   else       sprintf (name, "c_respZ%d", m);
1262   TCanvas *pad = new TCanvas(name, name, 700, 500);
1263   pad->SetRightMargin(0.10); pad->SetTopMargin(0.10);
1264   double yl = 0.75 - 0.025*hist.size();
1265   TLegend *legend = new TLegend(0.40, yl, 0.90, 0.75);
1266   legend->SetFillColor(kWhite);
1267   for (unsigned int j=0; j<hist.size(); ++j) {
1268     if (j == 0) hist[j]->Draw("");
1269     else        hist[j]->Draw("sames");
1270     char title[100];
1271     int j1 = comb[j]/100 - 1;
1272     int j2 = comb[j]%100 - 1;
1273     sprintf (title, "%s (%s)", dtitl[j1].c_str(), stitl[j2].c_str());
1274     legend->AddEntry(hist[j],title,"lp");
1275   }
1276   pad->Update();
1277   double xmax = 0.90;
1278   for (unsigned int k=0; k<hist.size(); ++k) {
1279     TPaveStats* st1 = (TPaveStats*)hist[k]->GetListOfFunctions()->FindObject("stats");
1280     if (st1 != NULL) {
1281       st1->SetLineColor(colors[k]);
1282       st1->SetTextColor(colors[k]);
1283       st1->SetY1NDC(0.78); st1->SetY2NDC(0.90);
1284       st1->SetX1NDC(xmax-0.12); st1->SetX2NDC(xmax);
1285       xmax -= 0.12;
1286     }
1287   }
1288   if (!typex) {
1289     xmax = 0.90;
1290     double ymax = 0.895;
1291     int nbox(0);
1292     for (unsigned int k=0; k<hist.size(); ++k) {
1293       double mean(0), rms(0), ent(0);
1294       for (int i=1; i<=hist[k]->GetNbinsX(); ++i) {
1295     double error = hist[k]->GetBinError(i);
1296     double value = hist[k]->GetBinContent(i);
1297     mean += (value/(error*error));
1298     rms  += (value*value/(error*error));
1299     ent  += (1.0/(error*error));
1300       }
1301       mean /= ent;
1302       rms  /= ent;
1303       double err = std::sqrt((rms-mean*mean)/std::sqrt(ent));
1304       TPaveText *txt1 = new TPaveText(xmax-0.20,ymax-0.04,xmax,ymax,"blNDC");
1305       txt1->SetFillColor(0);
1306       txt1->SetLineColor(1);
1307       txt1->SetTextColor(colors[k]);
1308       char txt[80];
1309       sprintf (txt, "Mean = %5.3f #pm %5.3f", mean, err);
1310       txt1->AddText(txt);
1311       txt1->Draw("same");
1312       xmax -= 0.20; nbox++;
1313       if (nbox == 3) {
1314     ymax -= 0.04; xmax = 0.90; nbox = 0;
1315       }
1316       pad->Modified();
1317       pad->Update();
1318     }
1319   }
1320   pad->Modified();
1321   pad->Update();
1322   legend->Draw("same");
1323   pad->Update();
1324   if (save) {
1325     sprintf (name, "%s.gif", pad->GetName());
1326     pad->Print(name);
1327   }
1328 }
1329
1330 void PlotHists(std::string fname, int mode, bool save=false, bool typex=true) {
1331
1332   std::string dname[17] = {"D250","D252","D253","D500","D502","D503",
1333                "DAL0","DAL2","DAL3","PNP0","PNP2","PNP3",
1334                "P250","P252","P253","Q250","Q500"};
1335   std::string snamex[30] = {"X00","X01","X02","X03","X04","X05","X06","X07",
1336                 "X08","X09","X10","X11","X12","X13","X14","X15",
1337                 "X16","X17","X18","X19","X20","X21","X22","X23",
1338                 "X24","X25","X26","X27","X28","X29"};
1339   std::string snamez[30] = {"Z00","Z01","Z02","Z03","Z04","Z05","Z06","Z07",
1340                 "Z08","Z09","Z10","Z11","Z12","Z13","Z14","Z15",
1341                 "Z16","Z17","Z18","Z19","Z20","Z21","Z22","Z23",
1342                 "Z24","Z25","Z26","Z27","Z28","Z29"};
1343
1344   int comb[210] = {108,118,128,408,418,428,
1345            208,218,228,508,518,528,
1346            308,318,328,608,618,628,
1347            706,702,708,704,0,0,
1348            806,802,808,804,0,0,
1349            906,902,908,904,0,0,
1350            708,718,728,0,0,0,
1351            808,818,828,0,0,0,
1352            908,918,928,0,0,0,
1353            706,716,726,0,0,0,
1354            806,816,826,0,0,0,
1355            906,916,926,0,0,0,
1356            1308,2718,2728,0,0,0,
1357            1408,1418,1428,0,0,0,
1358            1508,1518,1528,0,0,0,
1359            1008,1108,1208,0,0,0,
1360            1608,1618,1628,0,0,0,
1361            1708,1718,1728,0,0,0,
1362            1606,1616,1626,0,0,0,
1363            1706,1716,1726,0,0,0,
1364            708,718,728,2708,2718,2728,
1365            808,818,828,1408,1418,1428,
1366            908,918,928,1508,1518,1528,
1367            708,1608,1708,1008,0,0,
1368            706,1606,1706,0,0,0,
1369                    1308,1408,1508,0,0,0,
1370            1306,1406,1506,0,0,0,
1371            701,711,721,0,0,0,
1372            702,712,722,0,0,0,
1373            703,713,723,0,0,0,
1374            704,714,724,0,0,0,
1375            705,715,725,0,0,0,
1376            707,717,727,0,0,0,
1377            709,719,729,0,0,0,
1378            710,720,730,0,0,0};
1379   int ncombs[35]={6,6,6,4,4,4,3,3,3,3,3,3,3,3,3,3,3,3,3,3,6,6,6,4,3,3,3,
1380           3,3,3,3,3,3,3,3};
1381
1382   gStyle->SetCanvasBorderMode(0); gStyle->SetCanvasColor(kWhite);
1383   gStyle->SetPadColor(kWhite);    gStyle->SetFillColor(kWhite);
1384   if (typex) {
1385     gStyle->SetOptStat(1100);
1386   } else {
1387     gStyle->SetOptStat(0);
1388     gStyle->SetOptFit(11);
1389   }
1390   TFile      *file = new TFile(fname.c_str());
1391   int modemin = (mode < 0 || mode > 34) ? 0 : mode;
1392   int modemax = (mode < 0 || mode > 34) ? 34 : mode;
1393   for (int m=modemin; m <= modemax; ++m) {
1394     int nh = ncombs[m];
1395     std::vector<TProfile*> histp;
1396     std::vector<TH1D*>     histh;
1397     std::vector<int>       icomb;
1398     char name[20];
1399     bool ok(true);
1400     for (int j=0; j<nh; ++j) {
1401       int j1 = comb[m*6+j]/100 - 1;
1402       int j2 = comb[m*6+j]%100 - 1;
1403       icomb.push_back(comb[m*6+j]);
1404       if (typex) {
1405     sprintf (name,"%seta%s",dname[j1].c_str(),snamex[j2].c_str());
1406     TProfile* hist1 = (TProfile*)file->FindObjectAny(name);
1407     std::cout << name << " read out at " << hist1 << std::endl;
1408     if (hist1 != 0) {
1409       TProfile* hist  = (TProfile*)hist1->Clone();
1410       histp.push_back(hist);
1411     } else {
1412       ok = false;
1413     }
1414       } else {
1415     sprintf (name,"%s%s",dname[j1].c_str(),snamez[j2].c_str());
1416     TH1D* hist1 = (TH1D*)file->FindObjectAny(name);
1417     std::cout << name << " read out at " << hist1 << std::endl;
1418     if (hist1 != 0) {
1419       TH1D* hist  = (TH1D*)hist1->Clone();
1420       histh.push_back(hist);
1421     } else {
1422       ok = false;
1423     }
1424       }
1425     }
1426     std::cout << "Mode " << m << " Flag " << ok << std::endl;
1427     if (ok) {
1428       if (typex) DrawHist(histp, icomb, m, true, save);
1429       else       DrawHist(histh, icomb, m, false, save);
1430     }
1431   }
1432 }
1433
1434 std::pair<double,double> GetMean(TH1D* hist, double xmin, double xmax) {
1435
1436   double mean(0), rms(0), err(0), wt(0);
1437   for (int i=1; i<=hist->GetNbinsX(); ++i) {
1438     if (((hist->GetBinLowEdge(i)) >= xmin) ||
1439     ((hist->GetBinLowEdge(i)+hist->GetBinWidth(i)) <= xmax)) {
1440       double cont = hist->GetBinContent(i);
1441       double valu = hist->GetBinLowEdge(i)+0.5*+hist->GetBinWidth(i);
1442       wt         += cont;
1443       mean       += (valu*cont);
1444       rms        += (valu*valu*cont);
1445     }
1446   }
1447   if (wt > 0) {
1448     mean /= wt;
1449     rms  /= wt;
1450     err   = std::sqrt((rms-mean*mean)/wt);
1451   }
1452   return std::pair<double,double>(mean,err);
1453 }
1454
1455 void FitHists(std::string infile, std::string outfile, std::string dname,
1456           int mode, bool append=true, bool saveAll=false) {
1457
1458   int iname[10]     = {0, 1, 2, 3, 4, 5, 6, 7, 8, 9};
1459   int checkmode[10] = {1000, 1000, 1000, 1000, 1000, 10, 1000, 100, 1000, 1};
1460   double xbins[9]   = {-21.0, -16.0, -12.0, -6.0, 0.0, 6.0, 12.0, 16.0, 21.0};
1461   double vbins[6]   = {0.0, 7.0, 10.0, 13.0, 16.0, 50.0};
1462   double dlbins[9]  = {0.0, 0.10, 0.20, 0.50, 1.0, 2.0, 2.5, 3.0, 10.0};
1463   std::string sname[4] = {"ratio","etaR", "dl1R","nvxR"};
1464   std::string lname[4] = {"Z", "E", "L", "V"};
1465   int         numb[4]  = {8, 8, 8, 5};
1466   TFile      *file = new TFile(infile.c_str());
1467   std::vector<TH1D*> hists;
1468   char name[100];
1469   if (file != 0) {
1470     for (int m1=0; m1<4; ++m1) {
1471       for (int m2=0; m2<10; ++m2) {
1472     sprintf (name, "%s%s%d0", dname.c_str(), sname[m1].c_str(), iname[m2]);
1473     TH1D* hist0 = (TH1D*)file->FindObjectAny(name);
1474     bool ok = ((hist0 != 0) && (hist0->GetEntries() > 25));
1475     if ((mode/checkmode[m2])%10 > 0 && ok) {
1476       TH1D* histo(0);
1477       for (int j=0; j<=numb[m1]; ++j) {
1478         sprintf (name, "%s%s%d%d", dname.c_str(), sname[m1].c_str(), iname[m2], j);
1479         TH1D* hist1 = (TH1D*)file->FindObjectAny(name);
1480         TH1D* hist  = (TH1D*)hist1->Clone();
1481         double value(0), error(0);
1482         if (hist->GetEntries() > 0) {
1483           value = hist->GetMean(); error = hist->GetRMS();
1484         }
1485         if (j == 0) {
1486           sprintf (name, "%s%s%d", dname.c_str(), lname[m1].c_str(), iname[m2]);
1487           if (m1 <= 1)      histo = new TH1D(name, hist->GetTitle(), numb[m1], xbins);
1488           else if (m1 == 2) histo = new TH1D(name, hist->GetTitle(), numb[m1], dlbins);
1489           else              histo = new TH1D(name, hist->GetTitle(), numb[m1], vbins);
1490         }
1491         if (hist->GetEntries() > 4) {
1492           double mean = hist->GetMean(), rms = hist->GetRMS();
1493           double LowEdge = mean - 1.5*rms;
1494           double HighEdge = mean + 2.0*rms;
1495           if (LowEdge < 0.15) LowEdge = 0.15;
1496           char option[20];
1497           if (hist0->GetEntries() > 100) sprintf (option, "+QRS");
1498           else                           sprintf (option, "+QRWLS");
1499           double minvalue(0.30);
1500           if (iname[m2] == 0) {
1501         sprintf (option, "+QRWLS");
1502         HighEdge = 0.9;
1503         minvalue = 0.10;
1504           }
1505           TFitResultPtr Fit = hist->Fit("gaus",option,"",LowEdge,HighEdge);
1506           value = Fit->Value(1);
1507           error = Fit->FitResult::Error(1);
1508           std::pair<double,double> meaner = GetMean(hist,0.2,2.0);
1509 //        std::cout << "Fit " << value << ":" << error << ":" << hist->GetMeanError() << " Mean " << meaner.first << ":" << meaner.second;
1510           if (value < minvalue || value > 2.0 || error > 0.5) {
1511         value = meaner.first; error = meaner.second;
1512           }
1513 //        std::cout << " Final " << value << ":" << error << std::endl;
1514         }
1515         if (j == 0) {
1516           hists.push_back(hist);
1517         } else {
1518           if (saveAll) hists.push_back(hist);
1519           histo->SetBinContent(j, value);
1520           histo->SetBinError(j, error);
1521           if (j == numb[m1]) {
1522         hists.push_back(histo);
1523           }
1524         }
1525       }
1526     }
1527       }
1528     }
1529     TFile* theFile(0);
1530     if (append) {
1531       theFile = new TFile(outfile.c_str(), "UPDATE");
1532     } else {
1533       theFile = new TFile(outfile.c_str(), "RECREATE");
1534     }
1535
1536     theFile->cd();
1537     for (unsigned int i=0; i<hists.size(); ++i) {
1538       TH1D* hnew = (TH1D*)hists[i]->Clone();
1539       hnew->Write();
1540     }
1541     theFile->Close();
1542     file->Close();
1543   }
1544 }
1545
1546 void FitCombineHist(std::string infile, std::string outfile, std::string dname1,
1547             std::string dname2, std::string dname, bool append=true) {
1548
1549   double xbins[43]= {-21.5,-20.5,-19.5,-18.5,-17.5,-16.5,-15.5,-14.5,-13.5,
1550              -12.5,-11.5,-10.5,-9.5,-8.5,-7.5,-6.5,-5.5,-4.5,-3.5,
1551              -2.5,-1.5,0.0,1.5,2.5,3.5,4.5,5.5,6.5,7.5,8.5,9.5,10.5,
1552              11.5,12.5,13.5,14.5,15.5,16.5,17.5,18.5,19.5,20.5,21.5};
1553   std::string sname("ratio"), lname("Z");
1554   int         iname(7), numb(42);
1555
1556   TFile      *file = new TFile(infile.c_str());
1557   std::vector<TH1D*> hists;
1558   char name[100];
1559   if (file != 0) {
1560     sprintf (name, "%s%s%d0", dname1.c_str(), sname.c_str(), iname);
1561     TH1D* hist1 = (TH1D*)file->FindObjectAny(name);
1562     bool ok1 = (hist1 != 0);
1563     sprintf (name, "%s%s%d0", dname2.c_str(), sname.c_str(), iname);
1564     TH1D* hist2 = (TH1D*)file->FindObjectAny(name);
1565     bool ok2 = (hist2 != 0);
1566     if (ok1 || ok2) {
1567       int    nbin;
1568       double xmin, xmax;
1569       if (ok1) {
1570     nbin = hist1->GetNbinsX();
1571     xmin = hist1->GetBinLowEdge(1);
1572     xmax = hist1->GetBinLowEdge(nbin)+hist1->GetBinWidth(nbin);
1573       } else {
1574     nbin = hist2->GetNbinsX();
1575     xmin = hist2->GetBinLowEdge(1);
1576     xmax = hist2->GetBinLowEdge(nbin)+hist1->GetBinWidth(nbin);
1577       }
1578       TH1D* hist0(0);
1579       TH1D* histo(0);
1580       for (int j=0; j<=numb; ++j) {
1581     sprintf (name, "%s%s%d%d", dname1.c_str(), sname.c_str(), iname, j);
1582     hist1 = (TH1D*)file->FindObjectAny(name);
1583     sprintf (name, "%s%s%d%d", dname2.c_str(), sname.c_str(), iname, j);
1584     hist2 = (TH1D*)file->FindObjectAny(name);
1585     if (hist1 != 0) {
1586       nbin = hist1->GetNbinsX();
1587       xmin = hist1->GetBinLowEdge(1);
1588       xmax = hist1->GetBinLowEdge(nbin)+hist1->GetBinWidth(nbin);
1589     } else {
1590       nbin = hist2->GetNbinsX();
1591       xmin = hist2->GetBinLowEdge(1);
1592       xmax = hist2->GetBinLowEdge(nbin)+hist1->GetBinWidth(nbin);
1593     }
1594     TH1D* hist;
1595     sprintf (name, "%s%s%d%d", dname.c_str(), sname.c_str(), iname, j);
1596     if (hist1 != 0) hist = new TH1D(name, hist1->GetTitle(),nbin,xmin,xmax);
1597     else            hist = new TH1D(name, hist2->GetTitle(),nbin,xmin,xmax);
1598     double total(0);
1599     for (int k=1; k<=nbin; ++k) {
1600       double value(0), error(0);
1601       if (ok1) {
1602         value += (hist1->GetBinContent(k));
1603         error += ((hist1->GetBinError(k))*(hist1->GetBinError(k)));
1604       }
1605       if (ok2) {
1606         value += (hist2->GetBinContent(k));
1607         error += ((hist2->GetBinError(k))*(hist2->GetBinError(k)));
1608       }
1609       hist->SetBinContent(k,value); hist->SetBinError(k,sqrt(error));
1610       total += value;
1611     }
1612     double value(0), error(0);
1613     if (hist->GetEntries() > 0) {
1614       value = hist->GetMean(); error = hist->GetRMS();
1615     }
1616     if (j == 0) {
1617       hist0 = hist;
1618       sprintf (name, "%s%s%d", dname.c_str(), lname.c_str(), iname);
1619       histo = new TH1D(name, hist->GetTitle(), numb, xbins);
1620     }
1621     if (total > 4) {
1622       double mean = hist->GetMean(), rms = hist->GetRMS();
1623       double LowEdge = mean - 1.5*rms;
1624       double HighEdge = mean + 2.0*rms;
1625       if (LowEdge < 0.15) LowEdge = 0.15;
1626       char option[20];
1627       if (total > 100) {
1628         sprintf (option, "+QRS");
1629       } else {
1630             sprintf (option, "+QRWLS");
1631         HighEdge= mean+1.5*rms;
1632       }
1633       double minvalue(0.30);
1634       TFitResultPtr Fit = hist->Fit("gaus",option,"",LowEdge,HighEdge);
1635       value = Fit->Value(1);
1636       error = Fit->FitResult::Error(1);
1637       std::pair<double,double> meaner = GetMean(hist,0.2,2.0);
1638 //    std::cout << "Fit " << value << ":" << error << ":" << hist->GetMeanError() << " Mean " << meaner.first << ":" << meaner.second;
1639       if (value < minvalue || value > 2.0 || error > 0.5) {
1640         value = meaner.first; error = meaner.second;
1641       }
1642 //    std::cout << " Final " << value << ":" << error << std::endl;
1643     }
1644     hists.push_back(hist);
1645     if (j != 0) {
1646       histo->SetBinContent(j, value);
1647       histo->SetBinError(j, error);
1648       if (j == numb) hists.push_back(histo);
1649     }
1650       }
1651     }
1652     TFile* theFile(0);
1653     if (append) {
1654       theFile = new TFile(outfile.c_str(), "UPDATE");
1655     } else {
1656       theFile = new TFile(outfile.c_str(), "RECREATE");
1657     }
1658
1659     theFile->cd();
1660     for (unsigned int i=0; i<hists.size(); ++i) {
1661       TH1D* hnew = (TH1D*)hists[i]->Clone();
1662       hnew->Write();
1663     }
1664     theFile->Close();
1665     file->Close();
1666   }
1667 }
1668
1669 void PlotFits(std::string fname, int mode, int mode2=0, bool save=false) {
1670
1671   std::string dname[48] = {"DJT0", "DJT2", "DEG0", "DEG2", "DV00", "DV02",
1672                "D250", "D252", "D500", "D502", "DAL0", "DAL2",
1673                "DAL3", "D120", "DV10", "DV12", "DV20", "DV22",
1674                "DV30", "DV32", "DNC0", "DNC2", "DHR0", "DHR2",
1675                "DPR0", "DPR2", "D750", "D752", "D5B0", "D5B2",
1676                "D5C0", "D5C2", "D5D0", "D5D2", "DED0", "DED2",
1677                "DMD0", "DMD2", "Q250", "Q500", "PNP0", "PNP2",
1678                "PNP3", "P250", "P252", "P253", "P750", "P752"};
1679   std::string snamex[10] = {"ratio00","ratio10","ratio20","ratio30",
1680                 "ratio40","ratio50","ratio60","ratio70",
1681                 "ratio80","ratio90"};
1682   std::string sname2[10] = {"Z0", "Z1", "Z2", "Z3", "Z4", "Z5", "Z6", "Z7",
1683                 "Z8", "Z9"};
1684   std::string dtitl[48] = {"Jet data (2015B,C,D) Method 0",
1685                "Jet data (2015B,C,D) Method 2",
1686                "e#gamma data (2015B,C,D) Method 0",
1687                "e#gamma data (2015B,C,D) Method 2",
1688                "(Jet + e#gamma) (2015B,C,D) Method 0",
1689                "(Jet + e#gamma) (2015B,C,D) Method 2",
1690                "Jet data (25 ns) Method 0",
1691                "Jet data (25 ns) Method 2",
1692                "Jet data (50 ns) Method 0",
1693                "Jet data (50 ns) Method 2",
1694                "Jet data (B+C) Method 0",
1695                "Jet data (B+C) Method 2",
1696                "Jet data (B+C) Method 3",
1697                "2012 Jet data",
1698                "Jet data (2015B+C) Version 1 Method 0",
1699                "Jet data (2015B+C) Version 1 Method 2",
1700                "Jet data (2015B+C) Version 2 Method 0",
1701                "Jet data (2015B+C) Version 2 Method 2",
1702                "Jet data (2015B+C) Version 3 Method 0",
1703                "Jet data (2015B+C) Version 3 Method 2",
1704                "Jet data (New Constants) Method 0",
1705                "Jet data (New Constants) Method 2",
1706                "Jet data (Reference) Method 0",
1707                "Jet data (Reference) Method 2",
1708                "Jet data (Old Reference) Method 0",
1709                "Jet data (Old Reference) Method 2",
1710                "Jet data (75X) Method 0",
1711                "Jet data (75X) Method 2",
1712                "Jet data (2015B) Method 0",
1713                "Jet data (2015B) Method 2",
1714                "Jet data (2015C) Method 0",
1715                "Jet data (2015C) Method 2",
1716                "Jet data (2015D) Method 0",
1717                "Jet data (2015D) Method 2",
1718                "e#gamma data (2015D) Method 0",
1719                "e#gamma data (2015D) Method 2",
1720                "Single #mu data Method 0",
1721                "Single #mu data Method 2",
1722                "QCD MC (25 ns) Method 0",
1723                "QCD MC (50 ns) Method 0",
1724                "#pi MC (No PU) Method 0",
1725                "#pi MC (No PU) Method 2",
1726                "#pi MC (No PU) Method 3",
1727                "#pi MC (25 ns) Method 0",
1728                "#pi MC (25 ns) Method 2",
1729                "#pi MC (25 ns) Method 3",
1730                "#pi MC (25 ns) Method 0",
1731                "#pi MC (25 ns) Method 2"};
1732   std::string stitl[10] = {"p 2:4 GeV","p 4:7 GeV",
1733                "p 7:10 GeV","p 10:15 GeV",
1734                "p 15:20 GeV","p 20:30 GeV",
1735                "p 30:40 GeV","p 40:60 GeV",
1736                "p 60:100 GeVV","all p"};
1737   std::string titles[40] = {"p 20:30 GeV Method 0", "p 40:60 GeV Method 0",
1738                 "All p Method 0",       "p 40:60 GeV Method 0",
1739                 "p 40:60 GeV Method 2", "p 20:30 GeV Method 0",
1740                 "p 20:30 GeV Method 2", "p 20:30 GeV Method 0",
1741                 "p 40:60 GeV Method 0", "All p Method 0",
1742                 "p 20:30 GeV Method 2", "p 40:60 GeV Method 2",
1743                 "All p Method 2",       "p 20:30 GeV Method 2",
1744                 "p 40:60 GeV Method 2", "All p Method 2",
1745                 "Method 0 (p 40:60 GeV)","Method 2 (p 40:60 GeV)",
1746                 "Method 2 (p 40:60 GeV)",
1747                 "#pi MC (No PU) p 20:30 GeV",
1748                 "#pi MC (No PU) p 40:60 GeV",
1749                 "#pi MC (No PU) all p",
1750                 "#pi MC (25 ns) p 20:30 GeV",
1751                 "#pi MC (25 ns) p 40:60 GeV",
1752                 "#pi MC (25 ns) all p",
1753                 "Data (p 40:60 GeV)", "Data (p 40:60 GeV)",
1754                 "Data (p 40:60 GeV)", "Data (p 40:60 GeV)",
1755                 "MC (p 40:60 GeV)",   "Data (p 2:4 GeV)",
1756                 "Data (p 4:7 GeV)",   "Data (p 7:10 GeV)",
1757                 "Data (p 10:15 GeV)", "Data (p 15:20 GeV)",
1758                 "Data (p 30:40 GeV)", "Data (p 60:100 GeV)",
1759                 "Data (All momenta)", "Method 0 (p 40:60 GeV)"
1760                 "Method 2 (p 40:60 GeV)"};
1761   int comb[240] = {1106,706,906,4406,3906,4006,
1762            1108,708,908,4408,3908,4008,
1763            1110,710,910,4410,3910,4010,
1764            1108,2708,2108,2308,2508,0,
1765            1208,2808,2208,2408,2608,0,
1766            1106,2706,2106,2306,2506,0,
1767            1206,2806,2206,2406,2606,0,
1768            2906,3106,3306,1406,0,0,
1769            2908,3108,3308,1408,0,0,
1770            2910,3110,3310,1410,0,0,
1771            3006,3206,3406,0,0,0,
1772            3008,3208,3408,0,0,0,
1773            3010,3210,3410,0,0,0,
1774            1206,806,1006,0,0,0,
1775            1208,808,1008,0,0,0,
1776            1210,810,1010,0,0,0,
1777            508,4708,0,0,0,0,
1778            608,4808,0,0,0,0,
1779            608,4508,0,0,0,0,
1780            4106,4206,4306,0,0,0,
1781            4108,4208,4308,0,0,0,
1782            4110,4210,4310,0,0,0,
1783            4406,4506,4606,4806,0,0,
1784            4408,4508,4608,4808,0,0,
1785            4410,4510,4610,4810,0,0,
1786            1108,1208,0,0,0,0,
1787            1508,1708,1908,0,0,0,
1788            1608,1808,2008,0,0,0,
1789            1508,1608,1708,1808,1908,2008,
1790            3908,4008,0,0,0,0,
1791            1101,1201,0,0,0,0,
1792            1102,1202,0,0,0,0,
1793            1103,1203,0,0,0,0,
1794            1104,1204,0,0,0,0,
1795            1105,1205,0,0,0,0,
1796            1107,1207,0,0,0,0,
1797            1109,1209,0,0,0,0,
1798            1110,1210,0,0,0,0,
1799            108,308,4708,0,0,0,
1800            208,408,4808,0,0,0};
1801   int ncombs[40]={6,6,6,5,5,5,5,3,3,3,3,3,3,3,3,3,2,2,2,3,3,3,4,4,4,
1802           2,3,3,6,2,2,2,2,2,2,2,2,2,3,3};
1803   int colors[6]={1,2,4,7,6,9};
1804   int mtype[6]={20,21,22,23,24,33};
1805   gStyle->SetCanvasBorderMode(0); gStyle->SetCanvasColor(kWhite);
1806   gStyle->SetPadColor(kWhite);    gStyle->SetFillColor(kWhite);
1807   if (mode2 == 0) gStyle->SetOptStat(1110);
1808   else            gStyle->SetOptStat(10);
1809   gStyle->SetOptFit(1);
1810
1811   TFile      *file = new TFile(fname.c_str());
1812   int modemin = (mode < 0 || mode > 39) ? 0 : mode;
1813   int modemax = (mode < 0 || mode > 39) ? 39 : mode;
1814   for (int m=modemin; m <= modemax; ++m) {
1815     int nh = ncombs[m];
1816     TH1D* hist[6];
1817     char name[40], namen[50];
1818     bool ok(true);
1819     double total[6];
1820     double ylow (0), yhigh(0), totalf(0);
1821     int ngood(0);
1822     for (int j=0; j<nh; ++j) {
1823       int j1 = comb[m*6+j]/100 - 1;
1824       int j2 = comb[m*6+j]%100 - 1;
1825       if (mode2 == 0)sprintf (name,"%s%s",dname[j1].c_str(),snamex[j2].c_str());
1826       else           sprintf (name,"%s%s",dname[j1].c_str(),sname2[j2].c_str());
1827       TH1D* h = (TH1D*)file->FindObjectAny(name);
1828       total[j] = 0;
1829       if (h != 0) {
1830     ngood++;
1831     sprintf (namen, "%sX", name);
1832     int    nbin = h->GetNbinsX();
1833     double xlow = h->GetBinLowEdge(1);
1834     double xhigh=  h->GetBinLowEdge(nbin)+h->GetBinWidth(nbin);
1835     hist[j] = new TH1D(namen,h->GetTitle(),nbin,xlow,xhigh);
1836     for (int i=1; i<=h->GetNbinsX(); ++i) {
1837       double value = h->GetBinContent(i);
1838       hist[j]->SetBinContent(i,value);
1839       hist[j]->SetBinError(i,h->GetBinError(i));
1840       if (mode2 == 0) {
1841         if (h->GetBinLowEdge(i) >= 0.25 &&
1842         h->GetBinLowEdge(i) < 2.0) total[j] += value;
1843       } else {
1844         total[j] += value;
1845       }
1846     }
1847     if (ngood == 1) totalf = total[j];
1848     if (mode2 == 0) {
1849       double scale = (total[j] > 0) ? totalf/total[j] : 1.0;
1850       for (int i=1; i<=hist[j]->GetNbinsX(); ++i) {
1851         hist[j]->SetBinContent(i,scale*hist[j]->GetBinContent(i));
1852         hist[j]->SetBinError(i,scale*hist[j]->GetBinError(i));
1853         if (hist[j]->GetBinLowEdge(i) >= 0.25 &&
1854         hist[j]->GetBinLowEdge(i) < 2.0) {
1855           if ((hist[j]->GetBinContent(i)) > yhigh)
1856         yhigh = (hist[j]->GetBinContent(i));
1857         }
1858       }
1859     } else {
1860       yhigh = 1.6;
1861     }
1862     if (mode2 == 0) {
1863       hist[j]->GetXaxis()->SetRangeUser(0.25,2.0);
1864       hist[j]->GetXaxis()->SetTitle("<E_{HCAL}/(p-E_{ECAL})>");
1865       hist[j]->GetYaxis()->SetTitle("Tracks");
1866       double mean = hist[j]->GetMean(), rms = hist[j]->GetRMS();
1867       double LowEdge = mean - 1.5*rms;
1868       double HighEdge = mean + 2.0*rms;
1869       if (LowEdge < 0.15) LowEdge = 0.15;
1870       char option[20];
1871       if (total[j] > 100) sprintf (option, "+QRS");
1872       else                sprintf (option, "+QRWLS");
1873       if (j2 <= 1) {
1874         sprintf (option, "+QRWLS");
1875         HighEdge = 0.9;
1876       }
1877       hist[j]->Fit("gaus",option,"",LowEdge,HighEdge);
1878       hist[j]->GetFunction("gaus")->SetLineColor(colors[j]);
1879     } else {
1880       hist[j]->GetYaxis()->SetTitle("<E_{HCAL}/(p-E_{ECAL})>");
1881       hist[j]->GetXaxis()->SetTitle("i#eta");
1882       hist[j]->Fit("pol0","q");
1883       hist[j]->GetFunction("pol0")->SetLineColor(colors[j]);
1884     }
1885     hist[j]->SetTitle(titles[m].c_str());
1886     hist[j]->SetLineColor(colors[j]);
1887     hist[j]->SetMarkerColor(colors[j]);
1888     hist[j]->SetMarkerStyle(mtype[j]);
1889       } else {
1890     ok = false;
1891       }
1892     }
1893     if (mode2 == 0) {
1894       if (yhigh < 150.) {
1895     int iy = int(0.1*yhigh)+2;
1896     yhigh  = double(iy*10);
1897       } else {
1898     int iy = int(0.01*yhigh)+2;
1899     yhigh  = double(iy*100);
1900       }
1901     } else {
1902       ylow = 0.4;
1903     }
1904     if (ngood > 0) {
1905       sprintf (name, "c_fitres%d", m);
1906       TCanvas *pad = new TCanvas(name, name, 700, 500);
1907       pad->SetRightMargin(0.10); pad->SetTopMargin(0.10);
1908       double ymin = (mode2==0) ? 0.76 : 0.84;
1909       double yminl= ymin - ngood*0.025 - 0.01;
1910       TLegend *legend = new TLegend(0.50, yminl, 0.90, ymin-0.01);
1911       legend->SetFillColor(kWhite);
1912       bool first(true);
1913       for (int j=0; j<nh; ++j) {
1914     if (hist[j] != 0) {
1915       hist[j]->GetYaxis()->SetRangeUser(ylow,yhigh);
1916       if (first) hist[j]->Draw("");
1917       else       hist[j]->Draw("sames");
1918       first = false;
1919       char title[100];
1920       int j1 = comb[m*6+j]/100 - 1;
1921       int j2 = comb[m*6+j]%100 - 1;
1922       sprintf (title, "%s (%s)", dtitl[j1].c_str(), stitl[j2].c_str());
1923       legend->AddEntry(hist[j],title,"lp");
1924     }
1925       }
1926       pad->Update();
1927       double xmax = 0.90;
1928       for (int k=0; k<nh; ++k) {
1929     if (hist[k] != 0) {
1930       TPaveStats* st1 = (TPaveStats*)hist[k]->GetListOfFunctions()->FindObject("stats");
1931       if (st1 != NULL) {
1932         st1->SetLineColor(colors[k]);
1933         st1->SetTextColor(colors[k]);
1934         st1->SetY1NDC(ymin);      st1->SetY2NDC(0.90);
1935         st1->SetX1NDC(xmax-0.12); st1->SetX2NDC(xmax);
1936         xmax -= 0.12;
1937       }
1938     }
1939       }
1940       pad->Modified();
1941       pad->Update();
1942       legend->Draw("same");
1943       pad->Update();
1944       if (save) {
1945     sprintf (name, "%s.gif", pad->GetName());
1946     pad->Print(name);
1947       }
1948     }
1949   }
1950 }
1951
1952 void PlotHist(std::string fname, int num, int mode, bool save) {
1953
1954   std::string name1[48] = {"DJT0", "DJT2", "DEG0", "DEG2", "DV00", "DV02",
1955                "D250", "D252", "D500", "D502", "DAL0", "DAL2",
1956                "DAL3", "D120", "DV10", "DV12", "DV20", "DV22",
1957                "DV30", "DV32", "DNC0", "DNC2", "DHR0", "DHR2",
1958                "DPR0", "DPR2", "D750", "D752", "D5B0", "D5B2",
1959                "D5C0", "D5C2", "D5D0", "D5D2", "DED0", "DED2",
1960                "DMD0", "DMD2", "Q250", "Q500", "PNP0", "PNP2",
1961                "PNP3", "P250", "P252", "P253", "P750", "P752"};
1962   std::string name2[40] = {"Z7", "E7", "L7", "V7", "Z9", "E9", "L9", "V9",
1963                "Z5", "E5", "L5", "V5", "Z6", "E6", "L6", "V6",
1964                "Z8", "E8", "L8", "V8", "Z4", "E4", "L4", "V4",
1965                "Z3", "E3", "L3", "V3", "Z2", "E2", "L2", "V2",
1966                "Z1", "E1", "L1", "V1", "Z0", "E0", "L0", "V0"};
1967   std::string name3[40] = {"ratio70", "etaR70", "dl1R70", "nvxR70",
1968                "ratio90", "etaR90", "dl1R90", "nvxR90",
1969                "ratio50", "etaR50", "dl1R50", "nvxR50",
1970                "ratio60", "etaR60", "dl1R60", "nvxR60",
1971                "ratio80", "etaR80", "dl1R80", "nvxR80",
1972                "ratio40", "etaR40", "dl1R40", "nvxR40",
1973                "ratio30", "etaR30", "dl1R30", "nvxR30",
1974                "ratio20", "etaR20", "dl1R20", "nvxR20",
1975                "ratio10", "etaR10", "dl1R10", "nvxR10",
1976                "ratio00", "etaR00", "dl1R00", "nvxR00"};
1977   std::string name4[29] = {"ratio71", "ratio72", "ratio73", "ratio74",
1978                "ratio75", "ratio76", "ratio77", "ratio78",
1979                "etaR71",  "etaR72",  "etaR73",  "etaR74",
1980                "etaR75",  "etaR76",  "etaR77",  "etaR78",
1981                "dl1R71",  "dl1R72",  "dl1R73",  "dl1R74",
1982                "dl1R75",  "dl1R76",  "dl1R77",  "dl1R78",
1983                "nvxR71",  "nvxR72",  "nvxR73",  "nvxR74",
1984                "nvxR75"};
1985   std::string name5[42] = {"ratio71",  "ratio72",  "ratio73",  "ratio74",
1986                "ratio75",  "ratio76",  "ratio77",  "ratio78",
1987                "ratio79",  "ratio710", "ratio711", "ratio712",
1988                "ratio713", "ratio714", "ratio715", "ratio716",
1989                "ratio717", "ratio718", "ratio719", "ratio720",
1990                "ratio721", "ratio722", "ratio723", "ratio724",
1991                "ratio725", "ratio726", "ratio727", "ratio728",
1992                "ratio729", "ratio730", "ratio731", "ratio732",
1993                "ratio733", "ratio734", "ratio735", "ratio736",
1994                "ratio737", "ratio738", "ratio739", "ratio740",
1995                "ratio741", "ratio742"};
1996   std::string titl1[48] = {"Jet data (2015B,C,D) Method 0",
1997                "Jet data (2015B,C,D) Method 2",
1998                "e#gamma data (2015B,C,D) Method 0",
1999                "e#gamma data (2015B,C,D) Method 2",
2000                "(Jet + e#gamma) (2015B,C,D) Method 0",
2001                "(Jet + e#gamma) (2015B,C,D) Method 2",
2002                "Jet data (25 ns) Method 0",
2003                "Jet data (25 ns) Method 2",
2004                "Jet data (50 ns) Method 0",
2005                "Jet data (50 ns) Method 2",
2006                "Jet data (B+C) Method 0",
2007                "Jet data (B+C) Method 2",
2008                "Jet data (B+C) Method 3",
2009                "2012 Jet data",
2010                "Jet data (2015B+C) Version 1 Method 0",
2011                "Jet data (2015B+C) Version 1 Method 2",
2012                "Jet data (2015B+C) Version 2 Method 0",
2013                "Jet data (2015B+C) Version 2 Method 2",
2014                "Jet data (2015B+C) Version 3 Method 0",
2015                "Jet data (2015B+C) Version 3 Method 2",
2016                "Jet data (New Constants) Method 0",
2017                "Jet data (New Constants) Method 2",
2018                "Jet data (Reference) Method 0",
2019                "Jet data (Reference) Method 2",
2020                "Jet data (Old Reference) Method 0",
2021                "Jet data (Old Reference) Method 2",
2022                "Jet data (75X) Method 0",
2023                "Jet data (75X) Method 2",
2024                "Jet data (2015B) Method 0",
2025                "Jet data (2015B) Method 2",
2026                "Jet data (2015C) Method 0",
2027                "Jet data (2015C) Method 2",
2028                "Jet data (2015D) Method 0",
2029                "Jet data (2015D) Method 2",
2030                "e#gamma data (2015D) Method 0",
2031                "e#gamma data (2015D) Method 2",
2032                "Single #mu data Method 0",
2033                "Single #mu data Method 2",
2034                "QCD MC (25 ns) Method 0",
2035                "QCD MC (50 ns) Method 0",
2036                "#pi MC (No PU) Method 0",
2037                "#pi MC (No PU) Method 2",
2038                "#pi MC (No PU) Method 3",
2039                "#pi MC (25 ns) Method 0",
2040                "#pi MC (25 ns) Method 2",
2041                "#pi MC (25 ns) Method 3",
2042                "#pi MC (25 ns) Method 0",
2043                "#pi MC (25 ns) Method 2"};
2044   std::string titl2[40] = {"i#eta", "i#eta", "d_{L1}", "# Vertex",
2045                "i#eta", "i#eta", "d_{L1}", "# Vertex",
2046                "i#eta", "i#eta", "d_{L1}", "# Vertex",
2047                "i#eta", "i#eta", "d_{L1}", "# Vertex",
2048                "i#eta", "i#eta", "d_{L1}", "# Vertex",
2049                "i#eta", "i#eta", "d_{L1}", "# Vertex",
2050                "i#eta", "i#eta", "d_{L1}", "# Vertex",
2051                "i#eta", "i#eta", "d_{L1}", "# Vertex",
2052                "i#eta", "i#eta", "d_{L1}", "# Vertex",
2053                "i#eta", "i#eta", "d_{L1}", "# Vertex"};
2054   static const int nmax1(154), nmax2(58), nmax3(84);
2055   int ncomb1[nmax1]     = {501,  502, 503, 504, 601, 602, 603, 604, 101, 102,
2056                103,  104, 201, 202, 203, 204, 301, 302, 303, 304,
2057                401,  402, 403, 404,4401,4402,4501,4502,4701,4702,
2058                4801,4802, 509, 510, 511, 512, 609, 610, 611, 612,
2059                109,  110, 111, 112, 209, 210, 211, 212, 309, 310,
2060                311,  312, 409, 410, 411, 412, 701, 702, 801, 802,
2061                901,  902,1001,1002,1101,1102,1201,1202,1401,1402,
2062                1501,1502,1601,1602,1701,1702,1801,1802,1901,1902,
2063                2001,2002,2101,2102,2201,2202,2301,2302,2401,2402,
2064                2501,2502,2601,2602,2701,2702,2801,2802,2901,2902,
2065                3001,3002,3101,3102,3201,3202,3301,3302,3401,3402,
2066                3501,3502,3601,3602,3701,3702,3801,3802,3901,3902,
2067                4001,4002,4101,4102,4201,4202,4301,4302,4601,4602,
2068                505,  506, 507, 508, 605, 606, 607, 608, 105, 106,
2069                107,  108, 205, 206, 207, 208, 305, 306, 307, 308,
2070                405,  406, 407, 408};
2071   int ncomb2[nmax2]     = {201,202,203,204,205,206,207,208,209,210,
2072                211,212,213,214,215,216,217,218,219,220,
2073                221,222,223,224,225,226,227,228,229,401,
2074                            402,403,404,405,406,407,408,409,410,411,
2075                            412,413,414,415,416,417,418,419,420,421,
2076                            422,423,424,425,426,427,428,429};
2077   int ncomb3[nmax3]     = {601,  602, 603, 604, 605, 606, 607, 608, 609, 610,
2078                611,  612, 613, 614, 615, 616, 617, 618, 619, 620,
2079                621,  622, 623, 624, 625, 626, 627, 628, 629, 630,
2080                            631,  632, 633, 634, 635, 636, 637, 638, 639, 640,
2081                            641,  642,
2082                4501,4502,4503,4504,4505,4506,4507,4508,4509,4510,
2083                4511,4512,4513,4514,4515,4516,4517,4518,4519,4520,
2084                4521,4522,4523,4524,4525,4526,4527,4528,4529,4530,
2085                            4531,4532,4533,4534,4535,4536,4537,4538,4539,4540,
2086                            4541,4542};
2087
2088   gStyle->SetCanvasBorderMode(0); gStyle->SetCanvasColor(kWhite);
2089   gStyle->SetPadColor(kWhite);    gStyle->SetFillColor(kWhite);
2090   gStyle->SetOptTitle(1);
2091   int iopt(1110), nmax(nmax1-1);
2092   if (mode == 2) {
2093     iopt = 1100; nmax = (nmax2-1);
2094   } else if (mode == 3) {
2095     iopt = 1100; nmax = (nmax3-1);  gStyle->SetOptTitle(0);
2096   } else if (mode == 0) {
2097     iopt = 10;
2098   }
2099   gStyle->SetOptStat(iopt);  gStyle->SetOptFit(1);
2100   TFile      *file = new TFile(fname.c_str());
2101   char name[100], namep[100], title[100];
2102   int kmin = (num >= 0 && num <= nmax) ? num : 0;
2103   int kmax = (num >= 0 && num <= nmax) ? num : nmax;
2104   for (int k=kmin; k<=kmax; ++k) {
2105     int i1(0), i2(0);
2106     if (mode == 3) {
2107       i1 = ((ncomb3[k]/100)%100-1); i2 = ((ncomb3[k]%100)-1);
2108     } else if (mode == 2) {
2109       i1 = ((ncomb2[k]/100)%100-1); i2 = ((ncomb2[k]%100)-1);
2110     } else {
2111       i1 = ((ncomb1[k]/100)%100-1); i2 = ((ncomb1[k]%100)-1);
2112     }
2113     if (mode == 3) {
2114       int eta = ((k%42) < 21) ? ((k%42)-21) : ((k%42)-20);
2115       sprintf (title, "%s (i#eta = %d)", titl1[i1].c_str(), eta);
2116     } else {
2117       sprintf (title, "%s", titl1[i1].c_str());
2118     }
2119     if (mode == 0) {
2120       sprintf (name,  "%s%s",  name1[i1].c_str(), name2[i2].c_str());
2121       sprintf (namep, "%s%s%d",name1[i1].c_str(), name2[i2].c_str(), mode);
2122     } else if (mode == 3) {
2123       sprintf (name,  "%s%s",  name1[i1].c_str(), name5[i2].c_str());
2124       sprintf (namep, "%s%s%d",name1[i1].c_str(), name5[i2].c_str(), mode);
2125     } else if (mode == 2) {
2126       sprintf (name,  "%s%s",  name1[i1].c_str(), name4[i2].c_str());
2127       sprintf (namep, "%s%s%d",name1[i1].c_str(), name4[i2].c_str(), mode);
2128     } else {
2129       sprintf (name,  "%s%s",  name1[i1].c_str(), name3[i2].c_str());
2130       sprintf (namep, "%s%s%d",name1[i1].c_str(), name3[i2].c_str(), mode);
2131     }
2132     TH1D* hist1 = (TH1D*)file->FindObjectAny(name);
2133     if (hist1 != 0) {
2134       TH1D* hist = (TH1D*)(hist1->Clone());
2135       TCanvas *pad = new TCanvas(namep, namep, 700, 500);
2136       pad->SetRightMargin(0.10);
2137       pad->SetTopMargin(0.10);
2138       if (mode == 0) {
2139     hist->GetXaxis()->SetTitle(titl2[i2].c_str());
2140     hist->GetYaxis()->SetTitle("<E_{HCAL}/(p-E_{ECAL})>");
2141     hist->GetYaxis()->SetRangeUser(0.4,1.6);
2142     hist->Fit("pol0","q");
2143       } else {
2144     hist->GetYaxis()->SetTitle("Tracks");
2145     hist->GetXaxis()->SetTitle("E_{HCAL}/(p-E_{ECAL})");
2146     hist->GetXaxis()->SetRangeUser(0.0,3.0);
2147       }
2148       hist->GetYaxis()->SetLabelOffset(0.005);
2149       hist->GetYaxis()->SetTitleOffset(1.20);
2150       hist->Draw();
2151       pad->Update();
2152       TPaveStats* st1 = (TPaveStats*)hist->GetListOfFunctions()->FindObject("stats");
2153       if (st1 != NULL) {
2154     double ymin = (mode == 0) ? 0.78 : 0.66;
2155     st1->SetY1NDC(ymin); st1->SetY2NDC(0.90);
2156     st1->SetX1NDC(0.65); st1->SetX2NDC(0.90);
2157       }
2158       TPaveText *txt1 = new TPaveText(0.50,0.60,0.90,0.65,"blNDC");
2159       txt1->SetFillColor(0);
2160       char txt[100];
2161       sprintf (txt, "%s", title);
2162       txt1->AddText(txt);
2163       txt1->Draw("same");
2164       pad->Modified();
2165       pad->Update();
2166       if (save) {
2167     sprintf (name, "%s.pdf", pad->GetName());
2168     pad->Print(name);
2169       }
2170     }
2171   }
2172 }
2173
2174 void doPlot(std::string outfile="histodata.root") {
2175   IsoTrkOfflineAnalyzer p1("/afs/cern.ch/work/s/sunanda/public/CMSSW_7_6_0_pre4/src/output_2015.root","HcalIsoTrkAnalyzerM0","DV10",100,true);
2176   p1.Loop();
2177   p1.SavePlot(outfile,false);
2178   IsoTrkOfflineAnalyzer p2("/afs/cern.ch/work/s/sunanda/public/CMSSW_7_6_0_pre4/src/output_2015.root","HcalIsoTrkAnalyzerM2","DV12",100,true);
2179   p2.Loop();
2180   p2.SavePlot(outfile,true);
2181   IsoTrkOfflineAnalyzer p3("/afs/cern.ch/work/s/sunanda/public/CMSSW_7_6_0_pre4/src/output_2015_All.root","HcalIsoTrkAnalyzerM0","DV20",100,true);
2182   p3.Loop();
2183   p3.SavePlot(outfile,true);
2184   IsoTrkOfflineAnalyzer p4("/afs/cern.ch/work/s/sunanda/public/CMSSW_7_6_0_pre4/src/output_2015_All.root","HcalIsoTrkAnalyzerM2","DV22",100,true);
2185   p4.Loop();
2186   p4.SavePlot(outfile,true);
2187   IsoTrkOfflineAnalyzer p5("/afs/cern.ch/work/s/sunanda/public/CMSSW_7_6_0_pre4/src/output_2015_JetHT.root","HcalIsoTrkAnalyzerM0","DV30",100,true);
2188   p5.Loop();
2189   p5.SavePlot(outfile,true);
2190   IsoTrkOfflineAnalyzer p6("/afs/cern.ch/work/s/sunanda/public/CMSSW_7_6_0_pre4/src/output_2015_JetHT.root","HcalIsoTrkAnalyzerM2","DV32",100,true);
2191   p6.Loop();
2192   p6.SavePlot(outfile,true);
2193   IsoTrkOfflineAnalyzer p7("/afs/cern.ch/work/s/sunanda/public/CMSSW_7_6_0_pre4/src/HLTnewcondition.root","HcalIsoTrkAnalyzerM0","DNC0",100,true);
2194   p7.Loop();
2195   p7.SavePlot(outfile,true);
2196   IsoTrkOfflineAnalyzer p8("/afs/cern.ch/work/s/sunanda/public/CMSSW_7_6_0_pre4/src/HLTnewcondition.root","HcalIsoTrkAnalyzerM2","DNC2",100,true);
2197   p8.Loop();
2198   p8.SavePlot(outfile,true);
2199   IsoTrkOfflineAnalyzer p9("/afs/cern.ch/work/s/sunanda/public/CMSSW_7_6_0_pre4/src/HLTreference.root","HcalIsoTrkAnalyzerM0","DHR0",100,true);
2200   p9.Loop();
2201   p9.SavePlot(outfile,true);
2202   IsoTrkOfflineAnalyzer p10("/afs/cern.ch/work/s/sunanda/public/CMSSW_7_6_0_pre4/src/HLTreference.root","HcalIsoTrkAnalyzerM2","DHR2",100,true);
2203   p10.Loop();
2204   p10.SavePlot(outfile,true);
2205   IsoTrkOfflineAnalyzer p11("/afs/cern.ch/work/s/sunanda/public/CMSSW_7_6_0_pre4/src/PRrefernce.root","HcalIsoTrkAnalyzerM0","DPR0",100,true);
2206   p11.Loop();
2207   p11.SavePlot(outfile,true);
2208   IsoTrkOfflineAnalyzer p12("/afs/cern.ch/work/s/sunanda/public/CMSSW_7_6_0_pre4/src/PRrefernce.root","HcalIsoTrkAnalyzerM2","DPR2",100,true);
2209   p12.Loop();
2210   p12.SavePlot(outfile,true);
2211   IsoTrkOfflineAnalyzer m1("/afs/cern.ch/work/s/sunanda/public/CMSSW_7_6_0_pre4/src/output_pi50M2NoPU.root", "method0", "PNP0",200,false);
2212   m1.Loop();
2213   m1.SavePlot(outfile,true);
2214   IsoTrkOfflineAnalyzer m2("/afs/cern.ch/work/s/sunanda/public/CMSSW_7_6_0_pre4/src/output_pi50M2NoPU.root", "method2", "PNP2",200,false);
2215   m2.Loop();
2216   m2.SavePlot(outfile,true);
2217   IsoTrkOfflineAnalyzer m3("/afs/cern.ch/work/s/sunanda/public/CMSSW_7_6_0_pre4/src/output_pi50M3NoPU.root", "method3", "PNP3",200,false);
2218   m3.Loop();
2219   m3.SavePlot(outfile,true);
2220   IsoTrkOfflineAnalyzer m4("/afs/cern.ch/work/s/sunanda/public/CMSSW_7_6_0_pre4/src/output_pi50M2Bx25ns.root", "method0", "P250",200,false);
2221   m4.Loop();
2222   m4.SavePlot(outfile,true);
2223   IsoTrkOfflineAnalyzer m5("/afs/cern.ch/work/s/sunanda/public/CMSSW_7_6_0_pre4/src/output_pi50M2Bx25ns.root", "method2", "P252",200,false);
2224   m5.Loop();
2225   m5.SavePlot(outfile,true);
2226   IsoTrkOfflineAnalyzer m6("/afs/cern.ch/work/s/sunanda/public/CMSSW_7_6_0_pre4/src/output_pi50M3Bx25ns.root", "method3", "P253",200,false);
2227   m6.Loop();
2228   m6.SavePlot(outfile,true);
2229   IsoTrkOfflineAnalyzer m7("/afs/cern.ch/work/g/gwalia/public/output_pi50_25nsBX_method2_753_2711.root", "method0", "P750",200,false);
2230   m7.Loop();
2231   m7.SavePlot(outfile,true);
2232   IsoTrkOfflineAnalyzer m8("/afs/cern.ch/work/g/gwalia/public/output_pi50_25nsBX_method2_753_2711.root", "method2", "P752",200,false);
2233   m8.Loop();
2234   m8.SavePlot(outfile,true);
2235   IsoTrkOfflineAnalyzer q1("/afs/cern.ch/work/s/sunanda/public/CMSSW_7_6_0_pre4/src/output_QCD25nsBX.root","method2","Q250",200,false);
2236   q1.Loop();
2237   q1.SavePlot(outfile,true);
2238   IsoTrkOfflineAnalyzer q2("/afs/cern.ch/work/s/sunanda/public/CMSSW_7_6_0_pre4/src/output_QCD50nsBX.root","method2","Q500",200,false);
2239   q2.Loop();
2240   q2.SavePlot(outfile,true);
2241   IsoTrkOfflineAnalyzer c0("/afs/cern.ch/work/g/gwalia/public/JetHT_2015B_75X.root","HcalIsoTrkAnalyzerM0","D750",100,true);
2242   c0.Loop();
2243   c0.SavePlot(outfile,true);
2244   IsoTrkOfflineAnalyzer c1("/afs/cern.ch/work/g/gwalia/public/JetHT_2015B_75X.root","HcalIsoTrkAnalyzerM2","D752",100,true);
2245   c1.Loop();
2246   c1.SavePlot(outfile,true);
2247   IsoTrkOfflineAnalyzer c2("/afs/cern.ch/work/g/gwalia/public/JetHT_76X_2015B.root","HcalIsoTrkAnalyzerM0","D5B0",100,true);
2248   c2.Loop();
2249   c2.SavePlot(outfile,true);
2250   IsoTrkOfflineAnalyzer c3("/afs/cern.ch/work/g/gwalia/public/JetHT_76X_2015B.root","HcalIsoTrkAnalyzerM2","D5B2",100,true);
2251   c3.Loop();
2252   c3.SavePlot(outfile,true);
2253   IsoTrkOfflineAnalyzer c4("/afs/cern.ch/work/g/gwalia/public/JetHT_76X_2015C.root","HcalIsoTrkAnalyzerM0","D5C0",100,true);
2254   c4.Loop();
2255   c4.SavePlot(outfile,true);
2256   IsoTrkOfflineAnalyzer c5("/afs/cern.ch/work/g/gwalia/public/JetHT_76X_2015C.root","HcalIsoTrkAnalyzerM2","D5C2",100,true);
2257   c5.Loop();
2258   c5.SavePlot(outfile,true);
2259   IsoTrkOfflineAnalyzer c6("/afs/cern.ch/work/s/shghosh/public/MYRESULTS/2015DJSON1.root","HcalIsoTrkAnalyzerM0","D5D0",100,true);
2260   c6.Loop();
2261   c6.SavePlot(outfile,true);
2262   IsoTrkOfflineAnalyzer c7("/afs/cern.ch/work/s/shghosh/public/MYRESULTS/2015DJSON1.root","HcalIsoTrkAnalyzerM2","D5D2",100,true);
2263   c7.Loop();
2264   c7.SavePlot(outfile,true);
2265   IsoTrkOfflineAnalyzer d1("/afs/cern.ch/work/s/sunanda/public/CMSSW_7_6_0_pre4/src/output_2015_25nsAll.root","HcalIsoTrkAnalyzerM0","D250",100,true);
2266   d1.Loop();
2267   d1.SavePlot(outfile,true);
2268   IsoTrkOfflineAnalyzer d2("/afs/cern.ch/work/s/sunanda/public/CMSSW_7_6_0_pre4/src/output_2015_25nsAll.root","HcalIsoTrkAnalyzerM2","D252",100,true);
2269   d2.Loop();
2270   d2.SavePlot(outfile,true);
2271   IsoTrkOfflineAnalyzer d3("/afs/cern.ch/work/s/sunanda/public/CMSSW_7_6_0_pre4/src/output_2015_50nsAll.root","HcalIsoTrkAnalyzerM0","D500",100,true);
2272   d3.Loop();
2273   d3.SavePlot(outfile,true);
2274   IsoTrkOfflineAnalyzer d4("/afs/cern.ch/work/s/sunanda/public/CMSSW_7_6_0_pre4/src/output_2015_50nsAll.root","HcalIsoTrkAnalyzerM2","D502",100,true);
2275   d4.Loop();
2276   d4.SavePlot(outfile,true);
2277   IsoTrkOfflineAnalyzer d5("/afs/cern.ch/work/s/sunanda/public/CMSSW_7_6_0_pre4/src/output_2015_All.root","HcalIsoTrkAnalyzerM0","DAL0",100,true);
2278   d5.Loop();
2279   d5.SavePlot(outfile,true);
2280   IsoTrkOfflineAnalyzer d6("/afs/cern.ch/work/s/sunanda/public/CMSSW_7_6_0_pre4/src/output_2015_All.root","HcalIsoTrkAnalyzerM2","DAL2",100,true);
2281   d6.Loop();
2282   d6.SavePlot(outfile,true);
2283   IsoTrkOfflineAnalyzer d7("/afs/cern.ch/work/s/sunanda/public/CMSSW_7_6_0_pre4/src/2012Total.root","HcalIsoTrkAnalyzer","D120",100,true);
2284   d7.Loop();
2285   d7.SavePlot(outfile,true);
2286   IsoTrkOfflineAnalyzer d8("/afs/cern.ch/work/s/sunanda/public/CMSSW_7_6_0_pre4/src/2015BCD.root","HcalIsoTrkAnalyzerM0","DV00",100,true);
2287   d8.Loop();
2288   d8.SavePlot(outfile,true);
2289   IsoTrkOfflineAnalyzer d9("/afs/cern.ch/work/s/sunanda/public/CMSSW_7_6_0_pre4/src/2015BCD.root","HcalIsoTrkAnalyzerM2","DV02",100,true);
2290   d9.Loop();
2291   d9.SavePlot(outfile,true);
2292   IsoTrkOfflineAnalyzer n1("/afs/cern.ch/work/s/shghosh/public/MYRESULTS/SingleMuon2015D.root","HcalIsoTrkAnalyzerM0","DMD0",100,true);
2293   n1.Loop();
2294   n1.SavePlot(outfile,true);
2295   IsoTrkOfflineAnalyzer n2("/afs/cern.ch/work/s/shghosh/public/MYRESULTS/SingleMuon2015D.root","HcalIsoTrkAnalyzerM2","DMD2",100,true);
2296   n2.Loop();
2297   n2.SavePlot(outfile,true);
2298   IsoTrkOfflineAnalyzer n3("/afs/cern.ch/work/s/shghosh/public/MYRESULTS/DoubleEG2015D.root","HcalIsoTrkAnalyzerM0","DED0",100,true);
2299   n3.Loop();
2300   n3.SavePlot(outfile,true);
2301   IsoTrkOfflineAnalyzer n4("/afs/cern.ch/work/s/shghosh/public/MYRESULTS/DoubleEG2015D.root","HcalIsoTrkAnalyzerM2","DED2",100,true);
2302   n4.Loop();
2303   n4.SavePlot(outfile,true);
2304   IsoTrkOfflineAnalyzer n5("/afs/cern.ch/work/s/shghosh/public/MYRESULTS/DoubleEGCOMB.root","HcalIsoTrkAnalyzerM0","DEG0",100,true);
2305   n5.Loop();
2306   n5.SavePlot(outfile,true);
2307   IsoTrkOfflineAnalyzer n6("/afs/cern.ch/work/s/shghosh/public/MYRESULTS/DoubleEGCOMB.root","HcalIsoTrkAnalyzerM2","DEG2",100,true);
2308   n6.Loop();
2309   n6.SavePlot(outfile,true);
2310   IsoTrkOfflineAnalyzer n7("/afs/cern.ch/work/g/gwalia/public/JetHT_BCD.root","HcalIsoTrkAnalyzerM0","DJT0",100,true);
2311   n7.Loop();
2312   n7.SavePlot(outfile,true);
2313   IsoTrkOfflineAnalyzer n8("/afs/cern.ch/work/g/gwalia/public/JetHT_BCD.root","HcalIsoTrkAnalyzerM2","DJT2",100,true);
2314   n8.Loop();
2315   n8.SavePlot(outfile,true);
2316   IsoTrkOfflineAnalyzer n9("/afs/cern.ch/work/g/gwalia/public/JetHT_Egamma_BCD.root","HcalIsoTrkAnalyzerM0","DXG0",100,true);
2317   n9.Loop();
2318   n9.SavePlot(outfile,true);
2319   IsoTrkOfflineAnalyzer n0("/afs/cern.ch/work/g/gwalia/public/JetHT_Egamma_BCD.root","HcalIsoTrkAnalyzerM2","DXG2",100,true);
2320   n0.Loop();
2321   n0.SavePlot(outfile,true);
2322 }
2323
2324 void doPlotN(std::string outfile="histodatan.root") {
2325   IsoTrkOfflineAnalyzer n5("/afs/cern.ch/work/s/shghosh/public/MYRESULTS/DoubleEGCOMB.root","HcalIsoTrkAnalyzerM0","DEG0",20100,true);
2326   n5.Loop();
2327   n5.SavePlot(outfile,false);
2328   IsoTrkOfflineAnalyzer n6("/afs/cern.ch/work/s/shghosh/public/MYRESULTS/DoubleEGCOMB.root","HcalIsoTrkAnalyzerM2","DEG2",20100,true);
2329   n6.Loop();
2330   n6.SavePlot(outfile,true);
2331   IsoTrkOfflineAnalyzer n7("/afs/cern.ch/work/g/gwalia/public/JetHT_BCD.root","HcalIsoTrkAnalyzerM0","DJT0",10100,true);
2332   n7.Loop();
2333   n7.SavePlot(outfile,true);
2334   IsoTrkOfflineAnalyzer n8("/afs/cern.ch/work/g/gwalia/public/JetHT_BCD.root","HcalIsoTrkAnalyzerM2","DJT2",10100,true);
2335   n8.Loop();
2336   n8.SavePlot(outfile,true);
2337   IsoTrkOfflineAnalyzer n9("/afs/cern.ch/work/g/gwalia/public/JetHT_Egamma_BCD.root","HcalIsoTrkAnalyzerM0","DXG0",10100,true);
2338   n9.Loop();
2339   n9.SavePlot(outfile,true);
2340   IsoTrkOfflineAnalyzer n0("/afs/cern.ch/work/g/gwalia/public/JetHT_Egamma_BCD.root","HcalIsoTrkAnalyzerM2","DXG2",10100,true);
2341   n0.Loop();
2342   n0.SavePlot(outfile,true);
2343   IsoTrkOfflineAnalyzer d8("/afs/cern.ch/work/s/sunanda/public/CMSSW_7_6_0_pre4/src/2015BCD.root","HcalIsoTrkAnalyzerM0","DV00",10100,true);
2344   d8.Loop();
2345   d8.SavePlot(outfile,true);
2346   IsoTrkOfflineAnalyzer d9("/afs/cern.ch/work/s/sunanda/public/CMSSW_7_6_0_pre4/src/2015BCD.root","HcalIsoTrkAnalyzerM2","DV02",10100,true);
2347   d9.Loop();
2348   d9.SavePlot(outfile,true);
2349   IsoTrkOfflineAnalyzer m7("/afs/cern.ch/work/g/gwalia/public/output_pi50_25nsBX_method2_753_2711.root", "method0", "P750",20200,false);
2350   m7.Loop();
2351   m7.SavePlot(outfile,true);
2352   IsoTrkOfflineAnalyzer m8("/afs/cern.ch/work/g/gwalia/public/output_pi50_25nsBX_method2_753_2711.root", "method2", "P752",20200,false);
2353   m8.Loop();
2354   m8.SavePlot(outfile,true);
2355 }
2356
2357 void doFit(std::string infile="histodata.root",
2358        std::string outfile="histfitnew3.root") {
2359   FitHists(infile,outfile,"D250",1111,false);
2360   FitHists(infile,outfile,"D252",1111,true);
2361   FitHists(infile,outfile,"D500",1111,true);
2362   FitHists(infile,outfile,"D502",1111,true);
2363   FitHists(infile,outfile,"DAL0",1111,true);
2364   FitHists(infile,outfile,"DAL2",1111,true);
2365   FitHists(infile,outfile,"D120",1111,true);
2366   FitHists(infile,outfile,"DV00",1111,true,true);
2367   FitHists(infile,outfile,"DV02",1111,true,true);
2368   FitHists(infile,outfile,"DJT0",1111,true,true);
2369   FitHists(infile,outfile,"DJT2",1111,true,true);
2370   FitHists(infile,outfile,"DEG0",1111,true,true);
2371   FitHists(infile,outfile,"DEG2",1111,true,true);
2372   FitHists(infile,outfile,"DV10",1111,true);
2373   FitHists(infile,outfile,"DV12",1111,true);
2374   FitHists(infile,outfile,"DV20",1111,true);
2375   FitHists(infile,outfile,"DV22",1111,true);
2376   FitHists(infile,outfile,"DV30",1111,true);
2377   FitHists(infile,outfile,"DV32",1111,true);
2378   FitHists(infile,outfile,"DNC0",1111,true);
2379   FitHists(infile,outfile,"DNC2",1111,true);
2380   FitHists(infile,outfile,"DHR0",1111,true);
2381   FitHists(infile,outfile,"DHR2",1111,true);
2382   FitHists(infile,outfile,"DPR0",1111,true);
2383   FitHists(infile,outfile,"DPR2",1111,true);
2384   FitHists(infile,outfile,"D750",1111,true);
2385   FitHists(infile,outfile,"D752",1111,true);
2386   FitHists(infile,outfile,"D5B0",1111,true);
2387   FitHists(infile,outfile,"D5B2",1111,true);
2388   FitHists(infile,outfile,"D5C0",1111,true);
2389   FitHists(infile,outfile,"D5C2",1111,true);
2390   FitHists(infile,outfile,"D5D0",1111,true);
2391   FitHists(infile,outfile,"D5D2",1111,true);
2392   FitHists(infile,outfile,"DED0",1111,true);
2393   FitHists(infile,outfile,"DED2",1111,true);
2394   FitHists(infile,outfile,"DMD0",1111,true);
2395   FitHists(infile,outfile,"DMD2",1111,true);
2396   FitHists(infile,outfile,"Q250",1111,true);
2397   FitHists(infile,outfile,"Q500",1111,true);
2398   FitHists(infile,outfile,"PNP0",1111,true);
2399   FitHists(infile,outfile,"PNP2",1111,true);
2400   FitHists(infile,outfile,"PNP3",1111,true);
2401   FitHists(infile,outfile,"P250",1111,true);
2402   FitHists(infile,outfile,"P252",1111,true);
2403   FitHists(infile,outfile,"P253",1111,true);
2404   FitHists(infile,outfile,"P750",1111,true,true);
2405   FitHists(infile,outfile,"P752",1111,true,true);
2406 }
2407
2408 void doFitN(std::string infile="histodatan.root",
2409         std::string outfile="histfitnew4.root") {
2410   FitCombineHist(infile,outfile,"DEG0","DJT0","DV00",false);
2411   FitCombineHist(infile,outfile,"DEG2","DJT2","DV02",true);
2412   FitCombineHist(infile,outfile,"P750","XXX0","P250",true);
2413   FitCombineHist(infile,outfile,"P752","XXX2","P252",true);
2414 }
2415
2416 void doGetEntry() {
2417   GetEntries m1("/afs/cern.ch/work/g/gwalia/public/output_pi50_25nsBX_method2_753_2711.root", "method2", false);
2418   m1.Loop();
2419 }