RecoMuon/test/makePlot.cpp

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 //#include "Riosteram.h"
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   //  PlotManager(std::vector<std::string> inFileNames, std::string outFileName);
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 // super-mkdir : recursively run mkdir in root file
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     // Extract command from buffer
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   // Set default values
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   // Check validity of objects
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   // Check bin size
0228   if (numeHist->GetNbinsX() != denoHist->GetNbinsX()) {
0229     cerr << "Bin size of two histograms are not same" << endl;
0230     return false;
0231   }
0232
0233   // Push to base directory
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   // Create new histogram
0247   TH1F* effHist = dynamic_cast<TH1F*>(numeHist->Clone(newHistName.c_str()));
0248
0249   // effHist->Divide(denoHist);
0250   // Set the error to binomial statistics
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   // Cosmetics
0265   //effHist->SetName(newHistName.c_str());
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   // Save histogram
0273   effHist->Write();
0274
0275   // Pop directory
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   // Check validity of objects
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   // Check bin size
0303   if (numeHist->GetNbinsX() != denoHist->GetNbinsX()) {
0304     cerr << "Bin size of two histograms are not same" << endl;
0305     return false;
0306   }
0307
0308   // Push to base directory
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   // Create new TGraphAsymmErrors
0322   TGraphAsymmErrors* effGraph = new TGraphAsymmErrors(numeHist, denoHist);
0323
0324   // Cosmetics
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   // Save histogram
0333   effGraph->Write();
0334
0335   // Pop directory
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   // Check validity of objects
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   // Check bin size
0363   if (numeHist->GetNbinsX() != denoHist->GetNbinsX()) {
0364     cerr << "Bin size of two histograms are not same" << endl;
0365     return false;
0366   }
0367
0368   // Push to base directory
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   // Create new histogram
0382   TH1F* fakeHist = dynamic_cast<TH1F*>(numeHist->Clone());
0383
0384   // effHist->Divide(denoHist);
0385   // Set the error to binomial statistics
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   // Cosmetics
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   // Save histogram
0408   fakeHist->Write();
0409
0410   // Pop directory
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   // Check validity of objects
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   // Push to base directory
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   // Create a function for resolution model
0448   TF1 gaus("gaus", "gaus");
0449   gaus.SetParameters(1.0, 0.0, 0.1);
0450   //gaus.SetRange(yMin, yMax);
0451
0452   // Do FitSlices.
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   // Cosmetics
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   // Save histograms
0476   meanHist->Write();
0477   widthHist->Write();
0478   chi2Hist->Write();
0479
0480   // Pop directory
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   // Push to base directory
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   // Save histogram
0514   saveObj->Write();
0515
0516   // Pop directory
0517   gDirectory->cd(pwd.c_str());
0518
0519   return true;
0520 }
0521
0522 TDirectory* mkdirs(TDirectory* dir, string path) {
0523   // Push to directory passed into argument
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 }