Hcal/test/makeHFPhase1PMTParamsPlots.cc

0001 // Various standard headers
0002 #include <iostream>
0003 #include <string>
0004 #include <fstream>
0005 #include <cctype>
0006 #include <cassert>
0007 #include <algorithm>
0008 #include <set>
0009
0010 // ROOT headers
0011 #include "TROOT.h"
0012 #include "TFile.h"
0013
0014 // Serializaton headers
0015 #include "CondFormats/Serialization/interface/eos/portable_iarchive.hpp"
0016 #include "CondFormats/Serialization/interface/eos/portable_oarchive.hpp"
0017
0018 // Local headers
0019 #include "CondTools/Hcal/interface/CmdLine.h"
0020 #include "CondTools/Hcal/interface/visualizeHFPhase1PMTParams.h"
0021 #include "CondTools/Hcal/interface/parseHcalDetId.h"
0022
0023 using namespace cmdline;
0024
0025 static void print_usage(const char* progname, const char* options, const char* description) {
0026   std::cout << "\nUsage: " << progname << ' ' << options << " [-r referenceFile] idFile cutFile outputFile\n\n"
0027             << "Program arguments are:\n\n"
0028             << "  idfile       is the name of the text file containing the list of PMTs\n"
0029             << "               for which the cuts should be visualized. This file should\n"
0030             << "               contain the list of detector ids, one per line, in the format\n\n"
0031             << "               ieta  iphi  depth  subdetector\n\n"
0032             << "               Entries in this file which do not correspond to HF will be\n"
0033             << "               ignored, and depth 3 and 4 will be converted into 1 and 2.\n"
0034             << "               Duplicate entries will be ignored as well. Character '#'\n"
0035             << "               at the beginning of a line can be used to add comments.\n\n"
0036             << "  cutfile      is the name of the binary file containing the cuts\n"
0037             << "               to visualize. This file should be created by the\n"
0038             << "               \"write_HFPhase1PMTParams\" program.\n\n"
0039             << "  outputFile   is the name of the root file to write\n\n"
0040             << "Program options are:\n\n"
0041             << description << '\n'
0042             << "  -r           (default is not to plot reference cuts) reference binary file\n"
0043             << std::endl;
0044 }
0045
0046 static bool is_comment_or_blank(const std::string& s) {
0047   if (s.empty())
0048     return true;
0049   if (s[0] == '#' || s[0] == '\0')
0050     return true;
0051   for (const char* pc = s.c_str(); *pc; ++pc)
0052     if (!isspace(*pc))
0053       return false;
0054   return true;
0055 }
0056
0057 static std::vector<HcalDetId> exctractHFPmtIDs(std::vector<HcalDetId>& input) {
0058   std::set<HcalDetId> s;
0059   const unsigned n = input.size();
0060   for (unsigned i = 0; i < n; ++i)
0061     if (input[i].subdet() == HcalForward)
0062       s.insert(input[i].baseDetId());
0063   return std::vector<HcalDetId>(s.begin(), s.end());
0064 }
0065
0066 int main(int argc, char* argv[]) {
0067   using namespace std;
0068
0069   // Parse the input arguments
0070   CmdLine cmdline(argc, argv);
0071
0072   VisualizationOptions options;
0073   if (argc == 1) {
0074     print_usage(cmdline.progname(), options.options(), options.description());
0075     return 0;
0076   }
0077
0078   // Arguments which may be specified on the command line
0079   std::string idFile, cutFile, outputFile, referenceFile;
0080
0081   try {
0082     options.load(cmdline);
0083     cmdline.option("-r") >> referenceFile;
0084
0085     cmdline.optend();
0086     if (cmdline.argc() != 3)
0087       throw CmdLineError("wrong number of command line arguments");
0088     cmdline >> idFile >> cutFile >> outputFile;
0089   } catch (const CmdLineError& e) {
0090     cerr << "Error in " << cmdline.progname() << ": " << e.str() << endl;
0091     return 1;
0092   }
0093
0094   // Load the vector of ids
0095   vector<HcalDetId> idVec;
0096   {
0097     ifstream is(idFile.c_str());
0098     if (!is.is_open()) {
0099       cerr << "Failed to open file \"" << idFile << "\". Exiting." << endl;
0100       return 1;
0101     }
0102     string line;
0103     for (unsigned lineNumber = 1; getline(is, line); ++lineNumber)
0104       if (!is_comment_or_blank(line)) {
0105         HcalDetId id = parseHcalDetId(line);
0106         if (id.rawId())
0107           idVec.push_back(id);
0108         else
0109           cerr << "Invalid detector id entry on line " << lineNumber << " : \"" << line << "\". Ignored." << endl;
0110       }
0111   }
0112
0113   // We just need HF PMT ids, not all possible ids
0114   const vector<HcalDetId>& pmtIds = exctractHFPmtIDs(idVec);
0115   if (pmtIds.empty()) {
0116     cerr << "No valid HF PMT ids found in file " << idFile << ". Exiting." << endl;
0117     return 1;
0118   }
0119
0120   // Load the table of cuts
0121   HFPhase1PMTParams cuts;
0122   {
0123     std::ifstream is(cutFile.c_str(), std::ios::binary);
0124     if (!is.is_open()) {
0125       cerr << "Failed to open file \"" << cutFile << "\". Exiting." << endl;
0126       return 1;
0127     }
0128     eos::portable_iarchive ar(is);
0129     ar & cuts;
0130   }
0131
0132   // Load the reference table of cuts
0133   HFPhase1PMTParams refcuts;
0134   HFPhase1PMTParams* refptr = nullptr;
0135   if (!referenceFile.empty()) {
0136     std::ifstream is(referenceFile.c_str(), std::ios::binary);
0137     if (!is.is_open()) {
0138       cerr << "Failed to open file \"" << referenceFile << "\". Exiting." << endl;
0139       return 1;
0140     }
0141     eos::portable_iarchive ar(is);
0142     ar & refcuts;
0143     refptr = &refcuts;
0144   }
0145
0146   // Initialize root
0147   TROOT root("showCuts", "Cut Illustrator");
0148   root.SetBatch(kTRUE);
0149
0150   // Open the output root file. Note that, due to idiosyncratic root
0151   // object ownership policies, root objects (histograms, etc) allocated
0152   // on the heap will be owned by this file.
0153   TFile of(outputFile.c_str(), "RECREATE");
0154   if (!of.IsOpen()) {
0155     cerr << "Failed to open file \"" << outputFile << "\". Exiting." << endl;
0156     return 1;
0157   }
0158
0159   visualizeHFPhase1PMTParams(pmtIds, cuts, options, refptr);
0160
0161   // Write the output file
0162   of.Write();
0163   of.Close();
0164
0165   return 0;
0166 }