IsolatedParticles/macros/StudyHLT.C

0001 ///////////////////////////////////////////////////////////////////////////////
0002 //
0003 // Analysis script to prepare useful histograms from the Tree produced by
0004 // the StudyHLT EDAnalyzer from data and MC files
0005 //
0006 // StudyHLT t(inFile, outFile, dirname, treeName)
0007 // t.Loop()
0008 //
0009 // where
0010 //   infile    string    Name of the input ROOT tree file
0011 //   outfile   string    Name of the output ROOT histogram file
0012 //   dirname   string    Name of the directory ("StudyHLT")
0013 //   treeName  string    Name of the tree ("testTree")
0014 //
0015 // In addition there are useful methods:
0016 //
0017 // void GetPUWeight(mcFile, dataFile, type, dirName)
0018 //      Calculates PilUp weights using distribution of number of vertex
0019 // where
0020 //   mcFile    string    Name of the input ROOT tree MC file
0021 //   dataFile  string    Name of the input ROOT tree data file
0022 //   type      int       Variable to use all PV (0) or Good vertex (1)
0023 //   dirName   string    Name of the directory ("StudyHLT")
0024 //
0025 // TCanvas* PlotHist(fileName, type)
0026 //      Makes a plot from the histogram file
0027 // where
0028 //   fileName  string    Name of the input file
0029 //   type      int       Variable to be plotted: goodPV (0), numberPV (1),
0030 //                       maxNearP_Ecal(2), maxNearP_Hcal (3), Track Momentum(4)
0031 ///////////////////////////////////////////////////////////////////////////////
0032
0033 #include "TROOT.h"
0034 #include "TChain.h"
0035 #include "TFile.h"
0036 #include "TH1.h"
0037 #include "TH2.h"
0038 #include "TPaveStats.h"
0039 #include "TPaveText.h"
0040 #include "TStyle.h"
0041 #include "TCanvas.h"
0042
0043 #include <vector>
0044 #include <iostream>
0045
0046 class StudyHLT {
0047 public:
0048   TTree *fChain;   //!pointer to the analyzed TTree or TChain
0049   Int_t fCurrent;  //!current Tree number in a TChain
0050
0051   // Declaration of leaf types
0052   Int_t tr_goodRun;
0053   Int_t tr_goodPV;
0054   Double_t tr_eventWeight;
0055   std::vector<std::string> *tr_TrigName;
0056   std::vector<double> *tr_TrkPt;
0057   std::vector<double> *tr_TrkP;
0058   std::vector<double> *tr_TrkEta;
0059   std::vector<double> *tr_TrkPhi;
0060   std::vector<int> *tr_TrkID;
0061   std::vector<double> *tr_MaxNearP31X31;
0062   std::vector<double> *tr_MaxNearHcalP7x7;
0063   std::vector<double> *tr_FE7x7P;
0064   std::vector<double> *tr_FE11x11P;
0065   std::vector<double> *tr_FE15x15P;
0066   std::vector<bool> *tr_SE7x7P;
0067   std::vector<bool> *tr_SE11x11P;
0068   std::vector<bool> *tr_SE15x15P;
0069   std::vector<double> *tr_H7x7;
0070   std::vector<double> *tr_H5x5;
0071   std::vector<double> *tr_H3x3;
0072   std::vector<int> *tr_iEta;
0073
0074   // List of branches
0075   TBranch *b_tr_goodRun;          //!
0076   TBranch *b_tr_goodPV;           //!
0077   TBranch *b_tr_eventWeight;      //!
0078   TBranch *b_tr_TrigName;         //!
0079   TBranch *b_tr_TrkPt;            //!
0080   TBranch *b_tr_TrkP;             //!
0081   TBranch *b_tr_TrkEta;           //!
0082   TBranch *b_tr_TrkPhi;           //!
0083   TBranch *b_tr_TrkID;            //!
0084   TBranch *b_tr_MaxNearP31X31;    //!
0085   TBranch *b_tr_MaxNearHcalP7x7;  //!
0086   TBranch *b_tr_TrkQuality;       //!
0087   TBranch *b_tr_TrkokECAL;        //!
0088   TBranch *b_tr_TrkokHCAL;        //!
0089   TBranch *b_tr_FE7x7P;           //!
0090   TBranch *b_tr_FE11x11P;         //!
0091   TBranch *b_tr_FE15x15P;         //!
0092   TBranch *b_tr_SE7x7P;           //!
0093   TBranch *b_tr_SE11x11P;         //!
0094   TBranch *b_tr_SE15x15P;         //!
0095   TBranch *b_tr_H7x7;             //!
0096   TBranch *b_tr_H5x5;             //!
0097   TBranch *b_tr_H3x3;             //!
0098   TBranch *b_tr_iEta;             //!
0099
0100   StudyHLT(std::string inFile, std::string outFile, std::string dirnam = "StudyHLT", std::string treeNam = "testTree");
0101   virtual ~StudyHLT();
0102   virtual Int_t Cut(Long64_t entry);
0103   virtual Int_t GetEntry(Long64_t entry);
0104   virtual Long64_t LoadTree(Long64_t entry);
0105   virtual void Init(TTree *tree);
0106   virtual void Loop();
0107   virtual Bool_t Notify();
0108   virtual void Show(Long64_t entry = -1);
0109
0110   std::string outFile_;
0111 };
0112
0113 StudyHLT::StudyHLT(std::string inFile, std::string outFile, std::string dirnam, std::string treeNam)
0114     : outFile_(outFile) {
0115   TFile *file = new TFile(inFile.c_str());
0116   TDirectory *dir = (TDirectory *)file->FindObjectAny(dirnam.c_str());
0117   std::cout << inFile << " file " << file << " " << dirnam << " " << dir << std::endl;
0118   TTree *tree = (TTree *)dir->Get(treeNam.c_str());
0119   std::cout << "Tree:" << treeNam << " at " << tree << std::endl;
0120   Init(tree);
0121 }
0122
0123 StudyHLT::~StudyHLT() {
0124   if (!fChain)
0125     return;
0126   delete fChain->GetCurrentFile();
0127 }
0128
0129 Int_t StudyHLT::GetEntry(Long64_t entry) {
0130   // Read contents of entry.
0131   if (!fChain)
0132     return 0;
0133   return fChain->GetEntry(entry);
0134 }
0135
0136 Long64_t StudyHLT::LoadTree(Long64_t entry) {
0137   // Set the environment to read one entry
0138   if (!fChain)
0139     return -5;
0140   Long64_t centry = fChain->LoadTree(entry);
0141   if (centry < 0)
0142     return centry;
0143   if (fChain->GetTreeNumber() != fCurrent) {
0144     fCurrent = fChain->GetTreeNumber();
0145     Notify();
0146   }
0147   return centry;
0148 }
0149
0150 void StudyHLT::Init(TTree *tree) {
0151   // The Init() function is called when the selector needs to initialize
0152   // a new tree or chain. Typically here the branch addresses and branch
0153   // pointers of the tree will be set.
0154   // It is normally not necessary to make changes to the generated
0155   // code, but the routine can be extended by the user if needed.
0156   // Init() will be called many times when running on PROOF
0157   // (once per file to be processed).
0158
0159   // Set object pointer
0160   tr_TrigName = 0;
0161   tr_TrkPt = 0;
0162   tr_TrkP = 0;
0163   tr_TrkEta = 0;
0164   tr_TrkPhi = 0;
0165   tr_TrkID = 0;
0166   tr_MaxNearP31X31 = 0;
0167   tr_MaxNearHcalP7x7 = 0;
0168   tr_FE7x7P = 0;
0169   tr_FE11x11P = 0;
0170   tr_FE15x15P = 0;
0171   tr_SE7x7P = 0;
0172   tr_SE11x11P = 0;
0173   tr_SE15x15P = 0;
0174   tr_H3x3 = 0;
0175   tr_H5x5 = 0;
0176   tr_H7x7 = 0;
0177
0178   // Set branch addresses and branch pointers
0179   if (!tree)
0180     return;
0181   fChain = tree;
0182   fCurrent = -1;
0183   fChain->SetMakeClass(1);
0184
0185   fChain->SetBranchAddress("tr_goodRun", &tr_goodRun, &b_tr_goodRun);
0186   fChain->SetBranchAddress("tr_goodPV", &tr_goodPV, &b_tr_goodPV);
0187   fChain->SetBranchAddress("tr_eventWeight", &tr_eventWeight, &b_tr_eventWeight);
0188   fChain->SetBranchAddress("tr_TrigName", &tr_TrigName, &b_tr_TrigName);
0189   fChain->SetBranchAddress("tr_TrkPt", &tr_TrkPt, &b_tr_TrkPt);
0190   fChain->SetBranchAddress("tr_TrkP", &tr_TrkP, &b_tr_TrkP);
0191   fChain->SetBranchAddress("tr_TrkEta", &tr_TrkEta, &b_tr_TrkEta);
0192   fChain->SetBranchAddress("tr_TrkPhi", &tr_TrkPhi, &b_tr_TrkPhi);
0193   fChain->SetBranchAddress("tr_TrkID", &tr_TrkID, &b_tr_TrkID);
0194   fChain->SetBranchAddress("tr_MaxNearP31X31", &tr_MaxNearP31X31, &b_tr_MaxNearP31X31);
0195   fChain->SetBranchAddress("tr_MaxNearHcalP7x7", &tr_MaxNearHcalP7x7, &b_tr_MaxNearHcalP7x7);
0196   fChain->SetBranchAddress("tr_FE7x7P", &tr_FE7x7P, &b_tr_FE7x7P);
0197   fChain->SetBranchAddress("tr_FE11x11P", &tr_FE11x11P, &b_tr_FE11x11P);
0198   fChain->SetBranchAddress("tr_FE15x15P", &tr_FE15x15P, &b_tr_FE15x15P);
0199   fChain->SetBranchAddress("tr_SE7x7P", &tr_SE7x7P, &b_tr_SE7x7P);
0200   fChain->SetBranchAddress("tr_SE11x11P", &tr_SE11x11P, &b_tr_SE11x11P);
0201   fChain->SetBranchAddress("tr_SE15x15P", &tr_SE15x15P, &b_tr_SE15x15P);
0202   fChain->SetBranchAddress("tr_H3x3", &tr_H3x3, &b_tr_H3x3);
0203   fChain->SetBranchAddress("tr_H5x5", &tr_H5x5, &b_tr_H5x5);
0204   fChain->SetBranchAddress("tr_H7x7", &tr_H7x7, &b_tr_H7x7);
0205   fChain->SetBranchAddress("tr_iEta", &tr_iEta, &b_tr_iEta);
0206   Notify();
0207 }
0208
0209 Bool_t StudyHLT::Notify() {
0210   // The Notify() function is called when a new file is opened. This
0211   // can be either for a new TTree in a TChain or when when a new TTree
0212   // is started when using PROOF. It is normally not necessary to make changes
0213   // to the generated code, but the routine can be extended by the
0214   // user if needed. The return value is currently not used.
0215
0216   return kTRUE;
0217 }
0218
0219 void StudyHLT::Show(Long64_t entry) {
0220   // Print contents of entry.
0221   // If entry is not specified, print current entry
0222   if (!fChain)
0223     return;
0224   fChain->Show(entry);
0225 }
0226
0227 Int_t StudyHLT::Cut(Long64_t) {
0228   // This function may be called from Loop.
0229   // returns  1 if entry is accepted.
0230   // returns -1 otherwise.
0231   return 1;
0232 }
0233
0234 void StudyHLT::Loop() {
0235   //create a new root file to store the output.
0236   TFile *f = new TFile(outFile_.c_str(), "RECREATE");
0237   if (fChain == 0)
0238     return;
0239
0240   //declare the histograms in the same way as is done in studyHLT
0241   std::string titls[6] = {"NoIso", "okEcal", "EcalCharIso", "HcalCharIso", "EcalNeutIso", "HcalNeutIso"};
0242   char name[40], htit[400];
0243   TH1D *h_pt[6], *h_p[6], *h_eta[6], *h_phi[6];
0244   for (int i = 0; i < 6; ++i) {
0245     sprintf(name, "h_pt_%s", titls[i].c_str());
0246     h_pt[i] = new TH1D(name, "", 400, 0, 200);
0247     sprintf(name, "h_p_%s", titls[i].c_str());
0248     h_p[i] = new TH1D(name, "", 400, 0, 200);
0249     sprintf(name, "h_eta_%s", titls[i].c_str());
0250     h_eta[i] = new TH1D(name, "", 60, -3.0, 3.0);
0251     sprintf(name, "h_phi_%s", titls[i].c_str());
0252     h_phi[i] = new TH1D(name, "", 100, -3.15, 3.15);
0253   }
0254   static const int nPBin = 10, nEtaBin = 4, nPVBin = 4;
0255   TH1D *h_energy[nPVBin + 1][nPBin][nEtaBin][6];
0256   double pBin[nPBin + 1] = {1.0, 2.0, 3.0, 4.0, 5.0, 6.0, 7.0, 9.0, 11.0, 15.0, 20.0};
0257   int etaBin[nEtaBin + 1] = {1, 7, 13, 17, 23};
0258   int pvBin[nPVBin + 1] = {1, 2, 3, 5, 100};
0259   std::string energyNames[6] = {
0260       "E_{7x7}", "H_{3x3}", "(E_{7x7}+H_{3x3})", "E_{11x11}", "H_{5x5}", "{E_{11x11}+H_{5x5})"};
0261   for (int i = 0; i < nPVBin; ++i) {
0262     for (int ip = 0; ip < nPBin; ++ip) {
0263       for (int ie = 0; ie < nEtaBin; ++ie) {
0264         for (int j = 0; j < 6; ++j) {
0265           sprintf(name, "h_energy_%d_%d_%d_%d", i, ip, ie, j);
0266           if (i < nPVBin) {
0267             sprintf(htit,
0268                     "%s/p (p=%4.1f:%4.1f; i#eta=%d:%d, PV=%d:%d)",
0269                     energyNames[j].c_str(),
0270                     pBin[ip],
0271                     pBin[ip + 1],
0272                     etaBin[ie],
0273                     (etaBin[ie + 1] - 1),
0274                     pvBin[i],
0275                     pvBin[i + 1]);
0276           } else {
0277             sprintf(htit,
0278                     "%s/p (p=%4.1f:%4.1f; i#eta=%d:%d, All PV)",
0279                     energyNames[j].c_str(),
0280                     pBin[ip],
0281                     pBin[ip + 1],
0282                     etaBin[ie],
0283                     (etaBin[ie + 1] - 1));
0284           }
0285           h_energy[i][ip][ie][j] = new TH1D(name, htit, 500, -0.1, 4.9);
0286           h_energy[i][ip][ie][j]->Sumw2();
0287         }
0288       }
0289     }
0290   }
0291   Long64_t nentries = fChain->GetEntries();
0292
0293   Long64_t nbytes = 0, nb = 0;
0294   for (Long64_t jentry = 0; jentry < nentries; jentry++) {
0295     Long64_t ientry = LoadTree(jentry);
0296     if (ientry < 0)
0297       break;
0298     nb = fChain->GetEntry(jentry);
0299     nbytes += nb;
0300     if ((tr_TrkPt->size() != tr_TrkP->size()) || (tr_TrkPt->size() != tr_TrkEta->size()) ||
0301         (tr_TrkPt->size() != tr_TrkPhi->size())) {
0302       std::cout << "Error processing samples " << std::endl;
0303       break;
0304     }  //matches if
0305
0306     //fill the distributions
0307     //loop over all the reco tracks
0308     for (unsigned int itk = 0; itk != tr_TrkPt->size(); ++itk) {
0309       h_pt[0]->Fill(((*tr_TrkPt)[itk]), tr_eventWeight);
0310       h_p[0]->Fill(((*tr_TrkP)[itk]), tr_eventWeight);
0311       h_eta[0]->Fill(((*tr_TrkEta)[itk]), tr_eventWeight);
0312       h_phi[0]->Fill(((*tr_TrkPhi)[itk]), tr_eventWeight);
0313       if ((*tr_TrkPt)[itk] > 1.0 && abs((*tr_TrkEta)[itk]) < 2.5) {
0314         h_pt[1]->Fill(((*tr_TrkPt)[itk]), tr_eventWeight);
0315         h_p[1]->Fill(((*tr_TrkP)[itk]), tr_eventWeight);
0316         h_eta[1]->Fill(((*tr_TrkEta)[itk]), tr_eventWeight);
0317         h_phi[1]->Fill(((*tr_TrkPhi)[itk]), tr_eventWeight);
0318         //condition of charged Isolation
0319         if ((*tr_MaxNearP31X31)[itk] < 0) {
0320           h_pt[2]->Fill(((*tr_TrkPt)[itk]), tr_eventWeight);
0321           h_p[2]->Fill(((*tr_TrkP)[itk]), tr_eventWeight);
0322           h_eta[2]->Fill(((*tr_TrkEta)[itk]), tr_eventWeight);
0323           h_phi[2]->Fill(((*tr_TrkPhi)[itk]), tr_eventWeight);
0324           //condition for HCal Charged Isolation
0325           if ((*tr_MaxNearHcalP7x7)[itk] < 0) {
0326             h_pt[3]->Fill(((*tr_TrkPt)[itk]), tr_eventWeight);
0327             h_p[3]->Fill(((*tr_TrkP)[itk]), tr_eventWeight);
0328             h_eta[3]->Fill(((*tr_TrkEta)[itk]), tr_eventWeight);
0329             h_phi[3]->Fill(((*tr_TrkPhi)[itk]), tr_eventWeight);
0330             //condition of Neutral Isolation
0331             if ((*tr_SE11x11P)[itk] && (*tr_SE15x15P)[itk] && ((*tr_FE15x15P)[itk] - (*tr_FE11x11P)[itk]) < 2.0) {
0332               h_pt[4]->Fill(((*tr_TrkPt)[itk]), tr_eventWeight);
0333               h_p[4]->Fill(((*tr_TrkP)[itk]), tr_eventWeight);
0334               h_eta[4]->Fill(((*tr_TrkEta)[itk]), tr_eventWeight);
0335               h_phi[4]->Fill(((*tr_TrkPhi)[itk]), tr_eventWeight);
0336               if (((*tr_H7x7)[itk] - (*tr_H5x5)[itk]) < 2.0) {
0337                 h_pt[5]->Fill(((*tr_TrkPt)[itk]), tr_eventWeight);
0338                 h_p[5]->Fill(((*tr_TrkP)[itk]), tr_eventWeight);
0339                 h_eta[5]->Fill(((*tr_TrkEta)[itk]), tr_eventWeight);
0340                 h_phi[5]->Fill(((*tr_TrkPhi)[itk]), tr_eventWeight);
0341                 int ip(-1), ie(-1), nPV(-1);
0342                 for (int i = 0; i < nPBin; ++i) {
0343                   if (((*tr_TrkP)[itk] >= pBin[i]) && ((*tr_TrkP)[itk] < pBin[i + 1])) {
0344                     ip = i;
0345                     break;
0346                   }
0347                 }
0348                 for (int i = 0; i < nEtaBin; ++i) {
0349                   if (((*tr_iEta)[itk] >= etaBin[i]) && ((*tr_iEta)[itk] < etaBin[i + 1])) {
0350                     ie = i;
0351                     break;
0352                   }
0353                 }
0354                 for (int i = 0; i < nPVBin; ++i) {
0355                   if (((tr_goodPV) >= pvBin[i]) && ((tr_goodPV) < pvBin[i + 1])) {
0356                     nPV = i;
0357                     break;
0358                   }
0359                 }
0360                 double den = 1.0 / ((*tr_TrkP)[itk]);
0361                 if ((ip >= 0) && (ie >= 0) && ((*tr_FE7x7P)[itk] > 0.02) && ((*tr_H3x3)[itk] > 0.1) && (nPV >= 0)) {
0362                   h_energy[nPV][ip][ie][0]->Fill(den * (*tr_FE7x7P)[itk], tr_eventWeight);
0363                   h_energy[nPV][ip][ie][1]->Fill(den * (*tr_H3x3)[itk], tr_eventWeight);
0364                   h_energy[nPV][ip][ie][2]->Fill(den * ((*tr_FE7x7P)[itk] + (*tr_H3x3)[itk]), tr_eventWeight);
0365                   h_energy[nPV][ip][ie][3]->Fill(den * (*tr_FE11x11P)[itk], tr_eventWeight);
0366                   h_energy[nPV][ip][ie][4]->Fill(den * (*tr_H5x5)[itk], tr_eventWeight);
0367                   h_energy[nPV][ip][ie][5]->Fill(den * ((*tr_FE11x11P)[itk] + (*tr_H5x5)[itk]), tr_eventWeight);
0368                   h_energy[nPVBin][ip][ie][0]->Fill(den * (*tr_FE7x7P)[itk], tr_eventWeight);
0369                   h_energy[nPVBin][ip][ie][1]->Fill(den * (*tr_H3x3)[itk], tr_eventWeight);
0370                   h_energy[nPVBin][ip][ie][2]->Fill(den * ((*tr_FE7x7P)[itk] + (*tr_H3x3)[itk]), tr_eventWeight);
0371                   h_energy[nPVBin][ip][ie][3]->Fill(den * (*tr_FE11x11P)[itk], tr_eventWeight);
0372                   h_energy[nPVBin][ip][ie][4]->Fill(den * (*tr_H5x5)[itk], tr_eventWeight);
0373                   h_energy[nPVBin][ip][ie][5]->Fill(den * ((*tr_FE11x11P)[itk] + (*tr_H5x5)[itk]), tr_eventWeight);
0374                 }
0375               }  //HCal Neutral Iso
0376             }    //neutral isolation
0377           }      //HCal Charged Iso
0378         }        //charged Iso
0379       }
0380     }  //end for loop
0381   }
0382   f->Write();
0383   f->Close();
0384 }
0385
0386 void GetPUWeight(std::string mcFile, std::string dataFile, int type = 0, std::string dirName = "StudyHLT") {
0387   std::string hName = (type == 0) ? "h_numberPV" : "h_goodPV";
0388   TFile *file1 = new TFile(mcFile.c_str());
0389   TDirectory *dir1 = (TDirectory *)file1->FindObjectAny(dirName.c_str());
0390   TH1D *histM = (TH1D *)dir1->Get(hName.c_str());
0391   TFile *file2 = new TFile(dataFile.c_str());
0392   TDirectory *dir2 = (TDirectory *)file2->FindObjectAny(dirName.c_str());
0393   TH1D *histD = (TH1D *)dir2->Get(hName.c_str());
0394   double scale = histM->Integral() / histD->Integral();
0395   std::vector<double> weight;
0396   for (int k = 1; k <= histM->GetNbinsX(); ++k) {
0397     double num = histD->GetBinContent(k);
0398     double den = histM->GetBinContent(k);
0399     double rat = (den > 0) ? (scale * num / den) : 0;
0400     weight.push_back(rat);
0401   }
0402   char buff[100];
0403   for (int k = 0; k < histM->GetNbinsX(); k += 10) {
0404     sprintf(buff,
0405             "%6.4f,%6.4f,%6.4f,%6.4f,%6.4f,%6.4f,%6.4f,%6.4f,%6.4f,%6.4f,",
0406             weight[k + 0],
0407             weight[k + 1],
0408             weight[k + 2],
0409             weight[k + 3],
0410             weight[k + 4],
0411             weight[k + 5],
0412             weight[k + 6],
0413             weight[k + 7],
0414             weight[k + 8],
0415             weight[k + 9]);
0416     std::cout << buff << std::endl;
0417   }
0418 }
0419
0420 TCanvas *PlotHist(std::string fileName, int type) {
0421   std::string names[5] = {"h_goodPV", "h_numberPV", "h_maxNearP_Ecal", "h_maxNearP_Hcal", "h_p_HcalNeutIso"};
0422   std::string xtitl[5] = {"# Good Primary Vertex",
0423                           "# Primary Vertex",
0424                           "Highest Track Momentum in ECAL Isolation Zone (GeV)",
0425                           "Highest Track Momentum in HCAL Isolation Zone (GeV)",
0426                           "Track Momentum (GeV)"};
0427   double xmin[5] = {0, 0, -2, -2, 0};
0428   double xmax[5] = {10, 10, 10, 10, 60};
0429
0430   TCanvas *pad(0);
0431   gStyle->SetCanvasBorderMode(0);
0432   gStyle->SetCanvasColor(kWhite);
0433   gStyle->SetPadColor(kWhite);
0434   gStyle->SetFillColor(kWhite);
0435   gStyle->SetOptTitle(0);
0436   gStyle->SetOptStat(1110);
0437   gStyle->SetOptFit(0);
0438   TFile *file = new TFile(fileName.c_str());
0439   if (file) {
0440     TH1D *hist = (TH1D *)file->FindObjectAny(names[type].c_str());
0441     if (hist) {
0442       char name[100];
0443       sprintf(name, "%s", names[type].c_str());
0444       pad = new TCanvas(name, name, 700, 500);
0445       pad->SetRightMargin(0.10);
0446       pad->SetTopMargin(0.10);
0447       pad->SetLogy();
0448       hist->SetLineColor(1);
0449       hist->SetMarkerColor(1);
0450       hist->SetMarkerStyle(20);
0451       hist->GetXaxis()->SetTitle(xtitl[type].c_str());
0452       hist->GetYaxis()->SetTitle("Events");
0453       hist->GetYaxis()->SetLabelOffset(0.005);
0454       hist->GetYaxis()->SetTitleOffset(1.20);
0455       hist->GetXaxis()->SetRangeUser(xmin[type], xmax[type]);
0456       hist->Draw();
0457       pad->Modified();
0458       pad->Update();
0459       TPaveStats *st1 = (TPaveStats *)hist->GetListOfFunctions()->FindObject("stats");
0460       std::cout << "Pad " << pad << " st " << st1 << std::endl;
0461       if (st1 != NULL) {
0462         st1->SetFillColor(kWhite);
0463         st1->SetLineColor(1);
0464         st1->SetTextColor(1);
0465         st1->SetY1NDC(0.78);
0466         st1->SetY2NDC(0.90);
0467         st1->SetX1NDC(0.60);
0468         st1->SetX2NDC(0.90);
0469       }
0470       pad->Modified();
0471       pad->Update();
0472     }
0473   }
0474   return pad;
0475 }