File indexing completed on 2021-02-14 14:14:05
0001
0002
0003
0004
0005
0006
0007
0008
0009
0010
0011
0012
0013
0014
0015
0016
0017
0018
0019
0020
0021
0022
0023
0024
0025
0026
0027
0028
0029
0030
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;
0049 Int_t fCurrent;
0050
0051
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
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,
0101 std::string dirnam="StudyHLT", std::string treeNam="testTree");
0102 virtual ~StudyHLT();
0103 virtual Int_t Cut(Long64_t entry);
0104 virtual Int_t GetEntry(Long64_t entry);
0105 virtual Long64_t LoadTree(Long64_t entry);
0106 virtual void Init(TTree *tree);
0107 virtual void Loop();
0108 virtual Bool_t Notify();
0109 virtual void Show(Long64_t entry = -1);
0110
0111 std::string outFile_;
0112 };
0113
0114 StudyHLT::StudyHLT(std::string inFile, std::string outFile, std::string dirnam,
0115 std::string treeNam) : outFile_(outFile) {
0116
0117 TFile *file = new TFile(inFile.c_str());
0118 TDirectory *dir = (TDirectory*)file->FindObjectAny(dirnam.c_str());
0119 std::cout << inFile << " file " << file << " " << dirnam << " " << dir
0120 << std::endl;
0121 TTree *tree = (TTree*)dir->Get(treeNam.c_str());
0122 std::cout << "Tree:" << treeNam << " at " << tree << std::endl;
0123 Init(tree);
0124 }
0125
0126 StudyHLT::~StudyHLT() {
0127 if (!fChain) return;
0128 delete fChain->GetCurrentFile();
0129 }
0130
0131 Int_t StudyHLT::GetEntry(Long64_t entry) {
0132
0133 if (!fChain) return 0;
0134 return fChain->GetEntry(entry);
0135 }
0136
0137 Long64_t StudyHLT::LoadTree(Long64_t entry) {
0138
0139 if (!fChain) return -5;
0140 Long64_t centry = fChain->LoadTree(entry);
0141 if (centry < 0) return centry;
0142 if (fChain->GetTreeNumber() != fCurrent) {
0143 fCurrent = fChain->GetTreeNumber();
0144 Notify();
0145 }
0146 return centry;
0147 }
0148
0149 void StudyHLT::Init(TTree *tree) {
0150
0151
0152
0153
0154
0155
0156
0157
0158
0159 tr_TrigName = 0;
0160 tr_TrkPt = 0;
0161 tr_TrkP = 0;
0162 tr_TrkEta = 0;
0163 tr_TrkPhi = 0;
0164 tr_TrkID = 0;
0165 tr_MaxNearP31X31 = 0;
0166 tr_MaxNearHcalP7x7 = 0;
0167 tr_FE7x7P = 0;
0168 tr_FE11x11P = 0;
0169 tr_FE15x15P = 0;
0170 tr_SE7x7P = 0;
0171 tr_SE11x11P = 0;
0172 tr_SE15x15P = 0;
0173 tr_H3x3 = 0;
0174 tr_H5x5 = 0;
0175 tr_H7x7 = 0;
0176
0177
0178 if (!tree) return;
0179 fChain = tree;
0180 fCurrent = -1;
0181 fChain->SetMakeClass(1);
0182
0183 fChain->SetBranchAddress("tr_goodRun", &tr_goodRun, &b_tr_goodRun);
0184 fChain->SetBranchAddress("tr_goodPV", &tr_goodPV, &b_tr_goodPV);
0185 fChain->SetBranchAddress("tr_eventWeight", &tr_eventWeight, &b_tr_eventWeight);
0186 fChain->SetBranchAddress("tr_TrigName", &tr_TrigName, &b_tr_TrigName);
0187 fChain->SetBranchAddress("tr_TrkPt", &tr_TrkPt, &b_tr_TrkPt);
0188 fChain->SetBranchAddress("tr_TrkP", &tr_TrkP, &b_tr_TrkP);
0189 fChain->SetBranchAddress("tr_TrkEta", &tr_TrkEta, &b_tr_TrkEta);
0190 fChain->SetBranchAddress("tr_TrkPhi", &tr_TrkPhi, &b_tr_TrkPhi);
0191 fChain->SetBranchAddress("tr_TrkID", &tr_TrkID, &b_tr_TrkID);
0192 fChain->SetBranchAddress("tr_MaxNearP31X31", &tr_MaxNearP31X31, &b_tr_MaxNearP31X31);
0193 fChain->SetBranchAddress("tr_MaxNearHcalP7x7", &tr_MaxNearHcalP7x7, &b_tr_MaxNearHcalP7x7);
0194 fChain->SetBranchAddress("tr_FE7x7P", &tr_FE7x7P, &b_tr_FE7x7P);
0195 fChain->SetBranchAddress("tr_FE11x11P", &tr_FE11x11P, &b_tr_FE11x11P);
0196 fChain->SetBranchAddress("tr_FE15x15P", &tr_FE15x15P, &b_tr_FE15x15P);
0197 fChain->SetBranchAddress("tr_SE7x7P", &tr_SE7x7P, &b_tr_SE7x7P);
0198 fChain->SetBranchAddress("tr_SE11x11P", &tr_SE11x11P, &b_tr_SE11x11P);
0199 fChain->SetBranchAddress("tr_SE15x15P", &tr_SE15x15P, &b_tr_SE15x15P);
0200 fChain->SetBranchAddress("tr_H3x3", &tr_H3x3, &b_tr_H3x3);
0201 fChain->SetBranchAddress("tr_H5x5", &tr_H5x5, &b_tr_H5x5);
0202 fChain->SetBranchAddress("tr_H7x7", &tr_H7x7, &b_tr_H7x7);
0203 fChain->SetBranchAddress("tr_iEta", &tr_iEta, &b_tr_iEta);
0204 Notify();
0205 }
0206
0207 Bool_t StudyHLT::Notify() {
0208
0209
0210
0211
0212
0213
0214 return kTRUE;
0215 }
0216
0217 void StudyHLT::Show(Long64_t entry) {
0218
0219
0220 if (!fChain) return;
0221 fChain->Show(entry);
0222 }
0223
0224 Int_t StudyHLT::Cut(Long64_t ) {
0225
0226
0227
0228 return 1;
0229 }
0230
0231 void StudyHLT::Loop() {
0232
0233
0234 TFile *f = new TFile(outFile_.c_str(), "RECREATE");
0235 if (fChain == 0) return;
0236
0237
0238 std::string titls[6] = {"NoIso","okEcal","EcalCharIso","HcalCharIso",
0239 "EcalNeutIso","HcalNeutIso"};
0240 char name[40], htit[400];
0241 TH1D *h_pt[6], *h_p[6], *h_eta[6], *h_phi[6];
0242 for (int i=0; i<6; ++i) {
0243 sprintf (name, "h_pt_%s", titls[i].c_str());
0244 h_pt[i] = new TH1D(name,"", 400, 0, 200);
0245 sprintf (name, "h_p_%s", titls[i].c_str());
0246 h_p[i] = new TH1D(name,"", 400, 0, 200);
0247 sprintf (name, "h_eta_%s", titls[i].c_str());
0248 h_eta[i] = new TH1D(name,"", 60, -3.0, 3.0);
0249 sprintf (name, "h_phi_%s", titls[i].c_str());
0250 h_phi[i] = new TH1D(name,"", 100, -3.15, 3.15);
0251 }
0252 static const int nPBin=10, nEtaBin=4, nPVBin=4;
0253 TH1D *h_energy[nPVBin+1][nPBin][nEtaBin][6];
0254 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};
0255 int etaBin[nEtaBin+1] = {1, 7, 13, 17, 23};
0256 int pvBin[nPVBin+1] = {1, 2, 3, 5, 100};
0257 std::string energyNames[6]={"E_{7x7}", "H_{3x3}", "(E_{7x7}+H_{3x3})",
0258 "E_{11x11}", "H_{5x5}", "{E_{11x11}+H_{5x5})"};
0259 for (int i=0; i<nPVBin; ++i) {
0260 for (int ip=0; ip<nPBin; ++ip) {
0261 for (int ie=0; ie<nEtaBin; ++ie) {
0262 for (int j=0; j<6; ++j) {
0263 sprintf(name, "h_energy_%d_%d_%d_%d", i, ip, ie, j);
0264 if (i < nPVBin) {
0265 sprintf(htit,"%s/p (p=%4.1f:%4.1f; i#eta=%d:%d, PV=%d:%d)",
0266 energyNames[j].c_str(),pBin[ip],pBin[ip+1],etaBin[ie],
0267 (etaBin[ie+1]-1), pvBin[i], pvBin[i+1]);
0268 } else {
0269 sprintf(htit,"%s/p (p=%4.1f:%4.1f; i#eta=%d:%d, All PV)",
0270 energyNames[j].c_str(),pBin[ip],pBin[ip+1],etaBin[ie],
0271 (etaBin[ie+1]-1));
0272 }
0273 h_energy[i][ip][ie][j] = new TH1D(name, htit, 500, -0.1, 4.9);
0274 h_energy[i][ip][ie][j]->Sumw2();
0275 }
0276 }
0277 }
0278 }
0279 Long64_t nentries = fChain->GetEntries();
0280
0281 Long64_t nbytes = 0, nb = 0;
0282 for (Long64_t jentry=0; jentry<nentries;jentry++) {
0283 Long64_t ientry = LoadTree(jentry);
0284 if (ientry < 0) break;
0285 nb = fChain->GetEntry(jentry); nbytes += nb;
0286 if ((tr_TrkPt->size() != tr_TrkP->size()) ||
0287 (tr_TrkPt->size() != tr_TrkEta->size()) ||
0288 (tr_TrkPt->size() != tr_TrkPhi->size())) {
0289 std::cout << "Error processing samples " << std::endl;
0290 break;
0291 }
0292
0293
0294
0295 for (unsigned int itk=0; itk !=tr_TrkPt->size(); ++itk){
0296
0297 h_pt[0]->Fill( ((*tr_TrkPt)[itk]), tr_eventWeight);
0298 h_p[0]->Fill( ((*tr_TrkP)[itk]), tr_eventWeight);
0299 h_eta[0]->Fill( ((*tr_TrkEta)[itk]), tr_eventWeight);
0300 h_phi[0]->Fill( ((*tr_TrkPhi)[itk]), tr_eventWeight);
0301 if ( (*tr_TrkPt)[itk] > 1.0 && abs((*tr_TrkEta)[itk]) < 2.5) {
0302 h_pt[1]->Fill( ((*tr_TrkPt)[itk]), tr_eventWeight);
0303 h_p[1]->Fill( ((*tr_TrkP)[itk]), tr_eventWeight);
0304 h_eta[1]->Fill( ((*tr_TrkEta)[itk]), tr_eventWeight);
0305 h_phi[1]->Fill( ((*tr_TrkPhi)[itk]), tr_eventWeight);
0306
0307 if ( (*tr_MaxNearP31X31)[itk] < 0) {
0308 h_pt[2]->Fill( ((*tr_TrkPt)[itk]), tr_eventWeight);
0309 h_p[2]->Fill( ((*tr_TrkP)[itk]), tr_eventWeight);
0310 h_eta[2]->Fill( ((*tr_TrkEta)[itk]), tr_eventWeight);
0311 h_phi[2]->Fill( ((*tr_TrkPhi)[itk]), tr_eventWeight);
0312
0313 if ( (*tr_MaxNearHcalP7x7)[itk] < 0) {
0314 h_pt[3]->Fill( ((*tr_TrkPt)[itk]), tr_eventWeight);
0315 h_p[3]->Fill( ((*tr_TrkP)[itk]), tr_eventWeight);
0316 h_eta[3]->Fill( ((*tr_TrkEta)[itk]), tr_eventWeight);
0317 h_phi[3]->Fill( ((*tr_TrkPhi)[itk]), tr_eventWeight);
0318
0319 if ( (*tr_SE11x11P)[itk] && (*tr_SE15x15P)[itk] &&
0320 ( (*tr_FE15x15P)[itk] - (*tr_FE11x11P)[itk]) < 2.0 ) {
0321 h_pt[4]->Fill( ((*tr_TrkPt)[itk]), tr_eventWeight);
0322 h_p[4]->Fill( ((*tr_TrkP)[itk]), tr_eventWeight);
0323 h_eta[4]->Fill( ((*tr_TrkEta)[itk]), tr_eventWeight);
0324 h_phi[4]->Fill( ((*tr_TrkPhi)[itk]), tr_eventWeight);
0325 if ( ((*tr_H7x7)[itk] - (*tr_H5x5)[itk] ) < 2.0) {
0326 h_pt[5]->Fill( ((*tr_TrkPt)[itk]), tr_eventWeight);
0327 h_p[5]->Fill( ((*tr_TrkP)[itk]), tr_eventWeight);
0328 h_eta[5]->Fill( ((*tr_TrkEta)[itk]), tr_eventWeight);
0329 h_phi[5]->Fill( ((*tr_TrkPhi)[itk]), tr_eventWeight);
0330 int ip(-1), ie(-1), nPV(-1);
0331 for (int i=0; i<nPBin; ++i) {
0332 if (((*tr_TrkP)[itk] >= pBin[i]) &&
0333 ((*tr_TrkP)[itk] < pBin[i+1])) { ip = i; break; }
0334 }
0335 for (int i=0; i<nEtaBin; ++i) {
0336 if (((*tr_iEta)[itk] >= etaBin[i]) &&
0337 ((*tr_iEta)[itk] < etaBin[i+1])) { ie = i; break; }
0338 }
0339 for (int i=0; i<nPVBin; ++i) {
0340 if (((tr_goodPV) >= pvBin[i]) &&
0341 ((tr_goodPV) < pvBin[i+1])) { nPV = i; break; }
0342 }
0343 double den = 1.0/((*tr_TrkP)[itk]);
0344 if ((ip >= 0) && (ie >= 0) && ((*tr_FE7x7P)[itk] > 0.02)
0345 && ((*tr_H3x3)[itk] > 0.1) && (nPV >= 0)) {
0346 h_energy[nPV][ip][ie][0]->Fill(den*(*tr_FE7x7P)[itk], tr_eventWeight);
0347 h_energy[nPV][ip][ie][1]->Fill(den*(*tr_H3x3)[itk], tr_eventWeight);
0348 h_energy[nPV][ip][ie][2]->Fill(den*((*tr_FE7x7P)[itk]+(*tr_H3x3)[itk]), tr_eventWeight);
0349 h_energy[nPV][ip][ie][3]->Fill(den*(*tr_FE11x11P)[itk], tr_eventWeight);
0350 h_energy[nPV][ip][ie][4]->Fill(den*(*tr_H5x5)[itk], tr_eventWeight);
0351 h_energy[nPV][ip][ie][5]->Fill(den*((*tr_FE11x11P)[itk]+(*tr_H5x5)[itk]), tr_eventWeight);
0352 h_energy[nPVBin][ip][ie][0]->Fill(den*(*tr_FE7x7P)[itk], tr_eventWeight);
0353 h_energy[nPVBin][ip][ie][1]->Fill(den*(*tr_H3x3)[itk], tr_eventWeight);
0354 h_energy[nPVBin][ip][ie][2]->Fill(den*((*tr_FE7x7P)[itk]+(*tr_H3x3)[itk]), tr_eventWeight);
0355 h_energy[nPVBin][ip][ie][3]->Fill(den*(*tr_FE11x11P)[itk], tr_eventWeight);
0356 h_energy[nPVBin][ip][ie][4]->Fill(den*(*tr_H5x5)[itk], tr_eventWeight);
0357 h_energy[nPVBin][ip][ie][5]->Fill(den*((*tr_FE11x11P)[itk]+(*tr_H5x5)[itk]), tr_eventWeight);
0358 }
0359 }
0360 }
0361 }
0362 }
0363 }
0364 }
0365 }
0366 f->Write();
0367 f->Close();
0368 }
0369
0370 void GetPUWeight(std::string mcFile, std::string dataFile, int type=0,
0371 std::string dirName="StudyHLT") {
0372
0373 std::string hName = (type == 0) ? "h_numberPV" : "h_goodPV";
0374 TFile *file1 = new TFile(mcFile.c_str());
0375 TDirectory *dir1 = (TDirectory*)file1->FindObjectAny(dirName.c_str());
0376 TH1D *histM = (TH1D*)dir1->Get(hName.c_str());
0377 TFile *file2 = new TFile(dataFile.c_str());
0378 TDirectory *dir2 = (TDirectory*)file2->FindObjectAny(dirName.c_str());
0379 TH1D *histD = (TH1D*)dir2->Get(hName.c_str());
0380 double scale = histM->Integral()/histD->Integral();
0381
0382 std::vector<double> weight;
0383 for (int k=1; k<=histM->GetNbinsX(); ++k) {
0384 double num = histD->GetBinContent(k);
0385 double den = histM->GetBinContent(k);
0386 double rat = (den > 0) ? (scale*num/den) : 0;
0387 weight.push_back(rat);
0388
0389 }
0390 char buff[100];
0391 for (int k=0; k<histM->GetNbinsX(); k+=10) {
0392 sprintf(buff,"%6.4f,%6.4f,%6.4f,%6.4f,%6.4f,%6.4f,%6.4f,%6.4f,%6.4f,%6.4f,",
0393 weight[k+0],weight[k+1],weight[k+2],weight[k+3],weight[k+4],
0394 weight[k+5],weight[k+6],weight[k+7],weight[k+8],weight[k+9]);
0395 std::cout << buff << std::endl;
0396 }
0397 }
0398
0399 TCanvas* PlotHist(std::string fileName, int type) {
0400 std::string names[5] = {"h_goodPV", "h_numberPV", "h_maxNearP_Ecal",
0401 "h_maxNearP_Hcal","h_p_HcalNeutIso"};
0402 std::string xtitl[5] = {"# Good Primary Vertex", "# Primary Vertex",
0403 "Highest Track Momentum in ECAL Isolation Zone (GeV)",
0404 "Highest Track Momentum in HCAL Isolation Zone (GeV)",
0405 "Track Momentum (GeV)"};
0406 double xmin[5] = {0,0,-2,-2,0};
0407 double xmax[5] = {10,10,10,10,60};
0408
0409 TCanvas* pad(0);
0410 gStyle->SetCanvasBorderMode(0); gStyle->SetCanvasColor(kWhite);
0411 gStyle->SetPadColor(kWhite); gStyle->SetFillColor(kWhite);
0412 gStyle->SetOptTitle(0);
0413 gStyle->SetOptStat(1110); gStyle->SetOptFit(0);
0414 TFile* file = new TFile(fileName.c_str());
0415 if (file) {
0416 TH1D *hist = (TH1D*)file->FindObjectAny(names[type].c_str());
0417 if (hist) {
0418 char name[100];
0419 sprintf(name,"%s",names[type].c_str());
0420 pad = new TCanvas(name,name,700,500);
0421 pad->SetRightMargin(0.10); pad->SetTopMargin(0.10);
0422 pad->SetLogy();
0423 hist->SetLineColor(1);
0424 hist->SetMarkerColor(1);
0425 hist->SetMarkerStyle(20);
0426 hist->GetXaxis()->SetTitle(xtitl[type].c_str());
0427 hist->GetYaxis()->SetTitle("Events");
0428 hist->GetYaxis()->SetLabelOffset(0.005);
0429 hist->GetYaxis()->SetTitleOffset(1.20);
0430 hist->GetXaxis()->SetRangeUser(xmin[type],xmax[type]);
0431 hist->Draw();
0432 pad->Modified(); pad->Update();
0433 TPaveStats* st1 = (TPaveStats*)hist->GetListOfFunctions()->FindObject("stats");
0434 std::cout << "Pad " << pad << " st " << st1 << std::endl;
0435 if (st1 != NULL) {
0436 st1->SetFillColor(kWhite);
0437 st1->SetLineColor(1);
0438 st1->SetTextColor(1);
0439 st1->SetY1NDC(0.78); st1->SetY2NDC(0.90);
0440 st1->SetX1NDC(0.60); st1->SetX2NDC(0.90);
0441 }
0442 pad->Modified(); pad->Update();
0443 }
0444 }
0445 return pad;
0446 }