File indexing completed on 2024-04-06 12:33:12
0001 #include "TROOT.h"
0002 #include "TSystem.h"
0003 #include "TStyle.h"
0004 #include "TFile.h"
0005 #include "TH1F.h"
0006 #include "TH2F.h"
0007 #include "TTree.h"
0008
0009 #include "TPDF.h"
0010 #include "TF1.h"
0011 #include "TGraphAsymmErrors.h"
0012
0013 #include <iostream>
0014 #include <fstream>
0015 #include <string>
0016 #include <cmath>
0017 #include <sstream>
0018 #include <algorithm>
0019 #include <map>
0020 #include <cctype>
0021 #include <boost/tokenizer.hpp>
0022
0023 using namespace std;
0024 using namespace boost;
0025
0026 class PlotManager {
0027 public:
0028 PlotManager(const string& inFileName, const string& outFileName);
0029
0030 ~PlotManager();
0031
0032 void processCommand(istream& in);
0033
0034 protected:
0035 bool saveEfficiency(const string& histName,
0036 const string& histTitle,
0037 const string& numeHistName,
0038 const string& denoHistName);
0039 bool saveBayesEfficiency(const string& graphName,
0040 const string& graphTitle,
0041 const string& numeHistName,
0042 const string& denoHistName);
0043 bool saveFakeRate(const string& histName,
0044 const string& histTitle,
0045 const string& numeHistName,
0046 const string& denoHistName);
0047 bool saveResolution(const string& histName,
0048 const string& histTitle,
0049 const string& srcHistName,
0050 const char sliceDirection = 'Y');
0051 bool dumpObject(const string& objName, const string& objTitle, const string& srcObjName);
0052
0053 protected:
0054 bool setSrcFile(const string& fileName);
0055 bool setOutFile(const string& fileName);
0056
0057 protected:
0058 bool isSetup_;
0059
0060 TFile* theSrcFile_;
0061 TFile* theOutFile_;
0062 };
0063
0064
0065 TDirectory* mkdirs(TDirectory* dir, string path);
0066 string dirname(const string& path);
0067 string basename(const string& path);
0068
0069 int main(int argc, const char* argv[]) {
0070 string srcFileName, outFileName;
0071 if (argc > 2) {
0072 srcFileName = argv[1];
0073 outFileName = argv[2];
0074 } else {
0075 cout << "Usage : " << argv[0] << " sourceFile.root outputFile.root < commands.txt" << endl;
0076 cout << "Usage : " << argv[0] << " sourceFile.root outputFile.root" << endl;
0077 return 0;
0078 }
0079
0080 PlotManager plotMan(srcFileName, outFileName);
0081 plotMan.processCommand(cin);
0082 }
0083
0084 void PlotManager::processCommand(istream& in) {
0085 string buffer, cmd;
0086 while (!in.eof()) {
0087 buffer.clear();
0088 cmd.clear();
0089 getline(in, buffer);
0090
0091
0092 stringstream ss(buffer);
0093 ss >> cmd;
0094 if (cmd.empty() || cmd[0] == '#')
0095 continue;
0096 int (*pf)(int) = std::toupper;
0097 transform(cmd.begin(), cmd.end(), cmd.begin(), pf);
0098
0099 buffer.erase(0, cmd.size() + 1);
0100
0101 typedef escaped_list_separator<char> elsc;
0102 tokenizer<elsc> tokens(buffer, elsc("\\", " \t", "\""));
0103
0104 vector<tokenizer<elsc>::value_type> args;
0105
0106 for (tokenizer<elsc>::const_iterator tok_iter = tokens.begin(); tok_iter != tokens.end(); ++tok_iter) {
0107 args.push_back(*tok_iter);
0108 }
0109
0110 if (cmd == "EFFICIENCY" && args.size() >= 4) {
0111 if (!saveEfficiency(args[0], args[1], args[2], args[3])) {
0112 cerr << "Error : cannot make efficiency plot" << endl;
0113 }
0114 continue;
0115 }
0116
0117 if (cmd == "BAYESIANEFFICIENCY" && args.size() >= 4) {
0118 if (!saveBayesEfficiency(args[0], args[1], args[2], args[3])) {
0119 cerr << "Error : cannot make bayesian efficiency plot" << endl;
0120 }
0121 continue;
0122 }
0123
0124 if (cmd == "FAKERATE" && args.size() == 4) {
0125 if (!saveFakeRate(args[0], args[1], args[2], args[3])) {
0126 cerr << "Error : cannot make fakerate plot" << endl;
0127 }
0128 continue;
0129 }
0130
0131 if (cmd == "RESOLUTIONX" && args.size() == 3) {
0132 if (!saveResolution(args[0], args[1], args[2], 'X')) {
0133 cerr << "Error : cannot make resolution-X plot" << endl;
0134 }
0135 continue;
0136 }
0137
0138 if (cmd == "RESOLUTIONY" && args.size() == 3) {
0139 if (!saveResolution(args[0], args[1], args[2], 'Y')) {
0140 cerr << "Error : cannot make resolution-Y plot" << endl;
0141 }
0142 continue;
0143 }
0144
0145 if (cmd == "DUMP" && args.size() == 3) {
0146 if (!dumpObject(args[0], args[1], args[2])) {
0147 cerr << "Error : cannot copy histogram" << endl;
0148 }
0149 continue;
0150 }
0151
0152 cerr << "Unknown command <" << cmd << ">" << endl;
0153 }
0154 }
0155
0156 PlotManager::PlotManager(const string& srcFileName, const string& outFileName) {
0157
0158 theSrcFile_ = nullptr;
0159 theOutFile_ = nullptr;
0160
0161 isSetup_ = (setSrcFile(srcFileName) && setOutFile(outFileName));
0162 }
0163
0164 PlotManager::~PlotManager() {
0165 if (theSrcFile_) {
0166 theSrcFile_->Close();
0167 }
0168 if (theOutFile_) {
0169 theOutFile_->Write();
0170 theOutFile_->Close();
0171 }
0172 }
0173
0174 bool PlotManager::setSrcFile(const string& fileName) {
0175 if (theSrcFile_)
0176 theSrcFile_->Close();
0177
0178 string pwd(gDirectory->GetPath());
0179
0180 theSrcFile_ = new TFile(fileName.c_str());
0181 if (!theSrcFile_ || theSrcFile_->IsZombie())
0182 return false;
0183
0184 gDirectory->cd(pwd.c_str());
0185 return true;
0186 }
0187
0188 bool PlotManager::setOutFile(const string& fileName) {
0189 if (theOutFile_)
0190 theOutFile_->Close();
0191
0192 string pwd(gDirectory->GetPath());
0193
0194 theOutFile_ = new TFile(fileName.c_str(), "RECREATE");
0195 if (!theOutFile_ || theOutFile_->IsZombie())
0196 return false;
0197
0198 gDirectory->cd(pwd.c_str());
0199 return true;
0200 }
0201
0202 bool PlotManager::saveEfficiency(const string& histName,
0203 const string& histTitle,
0204 const string& numeHistName,
0205 const string& denoHistName) {
0206 if (!isSetup_)
0207 return false;
0208
0209 TH1F* numeHist = dynamic_cast<TH1F*>(theSrcFile_->Get(numeHistName.c_str()));
0210 TH1F* denoHist = dynamic_cast<TH1F*>(theSrcFile_->Get(denoHistName.c_str()));
0211
0212
0213 if (numeHist == nullptr) {
0214 cerr << "Cannot get object : " << numeHistName << endl;
0215 return false;
0216 }
0217 if (denoHist == nullptr) {
0218 cerr << "Cannot get object : " << denoHistName << endl;
0219 return false;
0220 }
0221
0222 if (!numeHist->IsA()->InheritsFrom("TH1") || !denoHist->IsA()->InheritsFrom("TH1") ||
0223 numeHist->IsA()->InheritsFrom("TH2") || denoHist->IsA()->InheritsFrom("TH2")) {
0224 return false;
0225 }
0226
0227
0228 if (numeHist->GetNbinsX() != denoHist->GetNbinsX()) {
0229 cerr << "Bin size of two histograms are not same" << endl;
0230 return false;
0231 }
0232
0233
0234 string pwd(gDirectory->GetPath());
0235
0236 string newHistPath = dirname(histName);
0237 string newHistName = basename(histName);
0238
0239 if (newHistPath.empty()) {
0240 theOutFile_->cd();
0241 } else if (theOutFile_->cd(newHistPath.c_str()) == kFALSE) {
0242 cout << "Cannot find directory, do mkdirs : " << newHistPath << endl;
0243 mkdirs(theOutFile_, newHistPath)->cd();
0244 }
0245
0246
0247 TH1F* effHist = dynamic_cast<TH1F*>(numeHist->Clone(newHistName.c_str()));
0248
0249
0250
0251 int nBinsX = effHist->GetNbinsX();
0252 for (int bin = 1; bin <= nBinsX; bin++) {
0253 float nNume = numeHist->GetBinContent(bin);
0254 float nDeno = denoHist->GetBinContent(bin);
0255 float eff = nDeno == 0 ? 0 : nNume / nDeno;
0256 float err = 0;
0257 if (nDeno != 0 && eff <= 1) {
0258 err = sqrt(eff * (1 - eff) / nDeno);
0259 }
0260 effHist->SetBinContent(bin, eff);
0261 effHist->SetBinError(bin, err);
0262 }
0263
0264
0265
0266 effHist->SetTitle(histTitle.c_str());
0267 effHist->SetMinimum(0.0);
0268 effHist->SetMaximum(1.0);
0269 effHist->GetXaxis()->SetTitle(numeHist->GetXaxis()->GetTitle());
0270 effHist->GetYaxis()->SetTitle("Efficiency");
0271
0272
0273 effHist->Write();
0274
0275
0276 gDirectory->cd(pwd.c_str());
0277
0278 return true;
0279 }
0280
0281 bool PlotManager::saveBayesEfficiency(const string& graphName,
0282 const string& graphTitle,
0283 const string& numeHistName,
0284 const string& denoHistName) {
0285 if (!isSetup_)
0286 return false;
0287
0288 TH1F* numeHist = (TH1F*)(theSrcFile_->Get(numeHistName.c_str()));
0289 TH1F* denoHist = (TH1F*)(theSrcFile_->Get(denoHistName.c_str()));
0290
0291
0292 if (numeHist == nullptr || denoHist == nullptr) {
0293 cerr << "Cannot get object : " << graphName << endl;
0294 return false;
0295 }
0296
0297 if (!numeHist->IsA()->InheritsFrom("TH1") || !denoHist->IsA()->InheritsFrom("TH1") ||
0298 numeHist->IsA()->InheritsFrom("TH2") || denoHist->IsA()->InheritsFrom("TH2")) {
0299 return false;
0300 }
0301
0302
0303 if (numeHist->GetNbinsX() != denoHist->GetNbinsX()) {
0304 cerr << "Bin size of two histograms are not same" << endl;
0305 return false;
0306 }
0307
0308
0309 string pwd(gDirectory->GetPath());
0310
0311 string newGraphPath = dirname(graphName);
0312 string newGraphName = basename(graphName);
0313
0314 if (newGraphPath.empty()) {
0315 theOutFile_->cd();
0316 } else if (theOutFile_->cd(newGraphPath.c_str()) == kFALSE) {
0317 cout << "Cannot find directory, do mkdirs" << endl;
0318 mkdirs(theOutFile_, newGraphPath)->cd();
0319 }
0320
0321
0322 TGraphAsymmErrors* effGraph = new TGraphAsymmErrors(numeHist, denoHist);
0323
0324
0325 effGraph->SetName(newGraphName.c_str());
0326 effGraph->SetTitle(graphTitle.c_str());
0327 effGraph->SetMinimum(0.8);
0328 effGraph->SetMaximum(1.0);
0329 effGraph->GetXaxis()->SetTitle(numeHist->GetXaxis()->GetTitle());
0330 effGraph->GetYaxis()->SetTitle("Efficiency");
0331
0332
0333 effGraph->Write();
0334
0335
0336 gDirectory->cd(pwd.c_str());
0337
0338 return true;
0339 }
0340
0341 bool PlotManager::saveFakeRate(const string& histName,
0342 const string& histTitle,
0343 const string& numeHistName,
0344 const string& denoHistName) {
0345 if (!isSetup_)
0346 return false;
0347
0348 TH1F* numeHist = (TH1F*)(theSrcFile_->Get(numeHistName.c_str()));
0349 TH1F* denoHist = (TH1F*)(theSrcFile_->Get(denoHistName.c_str()));
0350
0351
0352 if (numeHist == nullptr || denoHist == nullptr) {
0353 cerr << "Cannot get object : " << histName << endl;
0354 return false;
0355 }
0356
0357 if (!numeHist->IsA()->InheritsFrom("TH1") || !denoHist->IsA()->InheritsFrom("TH1") ||
0358 numeHist->IsA()->InheritsFrom("TH2") || denoHist->IsA()->InheritsFrom("TH2")) {
0359 return false;
0360 }
0361
0362
0363 if (numeHist->GetNbinsX() != denoHist->GetNbinsX()) {
0364 cerr << "Bin size of two histograms are not same" << endl;
0365 return false;
0366 }
0367
0368
0369 string pwd(gDirectory->GetPath());
0370
0371 string newHistPath = dirname(histName);
0372 string newHistName = basename(histName);
0373
0374 if (newHistPath.empty()) {
0375 theOutFile_->cd();
0376 } else if (theOutFile_->cd(newHistPath.c_str()) == kFALSE) {
0377 cout << "Cannot find directory, do mkdirs" << endl;
0378 mkdirs(theOutFile_, newHistPath)->cd();
0379 }
0380
0381
0382 TH1F* fakeHist = dynamic_cast<TH1F*>(numeHist->Clone());
0383
0384
0385
0386 int nBinsX = fakeHist->GetNbinsX();
0387 for (int bin = 1; bin <= nBinsX; bin++) {
0388 float nNume = numeHist->GetBinContent(bin);
0389 float nDeno = denoHist->GetBinContent(bin);
0390 float fakeRate = nDeno == 0 ? 0 : 1.0 - nNume / nDeno;
0391 float err = 0;
0392 if (nDeno != 0 && fakeRate <= 1) {
0393 err = sqrt(fakeRate * (1 - fakeRate) / nDeno);
0394 }
0395 fakeHist->SetBinContent(bin, fakeRate);
0396 fakeHist->SetBinError(bin, err);
0397 }
0398
0399
0400 fakeHist->SetName(newHistName.c_str());
0401 fakeHist->SetTitle(histTitle.c_str());
0402 fakeHist->SetMinimum(0.8);
0403 fakeHist->SetMaximum(1.0);
0404 fakeHist->GetXaxis()->SetTitle(numeHist->GetXaxis()->GetTitle());
0405 fakeHist->GetYaxis()->SetTitle("Efficiency");
0406
0407
0408 fakeHist->Write();
0409
0410
0411 gDirectory->cd(pwd.c_str());
0412
0413 return true;
0414 }
0415
0416 bool PlotManager::saveResolution(const string& histName,
0417 const string& histTitle,
0418 const string& srcHistName,
0419 const char sliceDirection) {
0420 if (!isSetup_)
0421 return false;
0422
0423 TH2F* srcHist = dynamic_cast<TH2F*>(theSrcFile_->Get(srcHistName.c_str()));
0424
0425
0426 if (srcHist == nullptr) {
0427 cerr << "Cannot get object : " << histName << endl;
0428 return false;
0429 }
0430
0431 if (srcHist->IsA()->InheritsFrom("TH2"))
0432 return false;
0433
0434
0435 string pwd(gDirectory->GetPath());
0436
0437 string newHistPath = dirname(histName);
0438 string newHistName = basename(histName);
0439
0440 if (newHistPath.empty()) {
0441 theOutFile_->cd();
0442 } else if (theOutFile_->cd(newHistPath.c_str()) == kFALSE) {
0443 cout << "Cannot find directory, do mkdirs" << endl;
0444 mkdirs(theOutFile_, newHistPath)->cd();
0445 }
0446
0447
0448 TF1 gaus("gaus", "gaus");
0449 gaus.SetParameters(1.0, 0.0, 0.1);
0450
0451
0452
0453 if (sliceDirection == 'X')
0454 srcHist->FitSlicesX(&gaus);
0455 else
0456 srcHist->FitSlicesY(&gaus);
0457
0458 TH1F* meanHist = dynamic_cast<TH1F*>(theOutFile_->Get((srcHistName + "_1").c_str()));
0459 TH1F* widthHist = dynamic_cast<TH1F*>(theOutFile_->Get((srcHistName + "_2").c_str()));
0460 TH1F* chi2Hist = dynamic_cast<TH1F*>(theOutFile_->Get((srcHistName + "_chi2").c_str()));
0461
0462
0463 meanHist->SetName((newHistName + "_Mean").c_str());
0464 widthHist->SetName((newHistName + "_Width").c_str());
0465 chi2Hist->SetName((newHistName + "_Chi2").c_str());
0466
0467 meanHist->SetTitle((histTitle + " Mean").c_str());
0468 widthHist->SetTitle((histTitle + " Width").c_str());
0469 chi2Hist->SetTitle((histTitle + " Chi2").c_str());
0470
0471 meanHist->GetYaxis()->SetTitle("Gaussian mean");
0472 widthHist->GetYaxis()->SetTitle("Gaussian width");
0473 chi2Hist->GetYaxis()->SetTitle("Gaussian fit #Chi^{2}");
0474
0475
0476 meanHist->Write();
0477 widthHist->Write();
0478 chi2Hist->Write();
0479
0480
0481 gDirectory->cd(pwd.c_str());
0482
0483 return true;
0484 }
0485
0486 bool PlotManager::dumpObject(const string& objName, const string& objTitle, const string& srcObjName) {
0487 if (!isSetup_)
0488 return false;
0489
0490
0491 string pwd(gDirectory->GetPath());
0492
0493 string newObjPath = dirname(objName);
0494 string newObjName = basename(objName);
0495
0496 if (newObjPath.empty()) {
0497 theOutFile_->cd();
0498 } else if (theOutFile_->cd(newObjPath.c_str()) == kFALSE) {
0499 cout << "Cannot find dir, do mkdirs : " << newObjPath << endl;
0500 mkdirs(theOutFile_, newObjPath)->cd();
0501 }
0502
0503 TNamed* srcObj = dynamic_cast<TNamed*>(theSrcFile_->Get(srcObjName.c_str()));
0504
0505 if (srcObj == nullptr) {
0506 cerr << "Cannot get object : " << srcObjName << endl;
0507 return false;
0508 }
0509
0510 TNamed* saveObj = dynamic_cast<TNamed*>(srcObj->Clone(newObjName.c_str()));
0511 saveObj->SetTitle(objTitle.c_str());
0512
0513
0514 saveObj->Write();
0515
0516
0517 gDirectory->cd(pwd.c_str());
0518
0519 return true;
0520 }
0521
0522 TDirectory* mkdirs(TDirectory* dir, string path) {
0523
0524 string pwd(gDirectory->GetPath());
0525
0526 while (true) {
0527 if (path.empty())
0528 break;
0529
0530 string::size_type slashPos = path.find_first_of('/');
0531 if (slashPos != string::npos) {
0532 string newDirName = path.substr(0, slashPos);
0533 TDirectory* tmpDir = dir->GetDirectory(newDirName.c_str());
0534 dir = (tmpDir == nullptr) ? dir->mkdir(newDirName.c_str()) : tmpDir;
0535
0536 path.erase(0, slashPos + 1);
0537 } else {
0538 TDirectory* tmpDir = dir->GetDirectory(path.c_str());
0539 dir = (tmpDir == nullptr) ? dir->mkdir(path.c_str()) : tmpDir;
0540
0541 break;
0542 }
0543 }
0544
0545 return dir;
0546 }
0547
0548 string dirname(const string& path) {
0549 string::size_type slashPos = path.find_last_of('/');
0550 if (slashPos == string::npos) {
0551 return "";
0552 }
0553
0554 return path.substr(0, slashPos);
0555 }
0556
0557 string basename(const string& path) {
0558 string::size_type slashPos = path.find_last_of('/');
0559 if (slashPos == string::npos) {
0560 return path;
0561 }
0562
0563 return path.substr(slashPos + 1, path.size());
0564 }