RivetInterface/plugins/DQMRivetClient.cc

0001 /*
0002  *  Class:DQMRivetClient
0003  *
0004  *
0005  *
0006  *  \author Junghwan Goh - SungKyunKwan University
0007  */
0008
0009 #include "GeneratorInterface/RivetInterface/interface/DQMRivetClient.h"
0010
0011 #include "DQMServices/Core/interface/DQMStore.h"
0012 #include "FWCore/ParameterSet/interface/ParameterSet.h"
0013 #include "FWCore/ServiceRegistry/interface/Service.h"
0014 #include "FWCore/MessageLogger/interface/MessageLogger.h"
0015
0016 #include <TH1F.h>
0017 #include <TClass.h>
0018 #include <TString.h>
0019 #include <TPRegexp.h>
0020
0021 #include <cmath>
0022 #include <boost/tokenizer.hpp>
0023
0024 using namespace std;
0025 using namespace edm;
0026
0027 typedef DQMRivetClient::MonitorElement ME;
0028
0029 DQMRivetClient::DQMRivetClient(const ParameterSet& pset) {
0030   typedef std::vector<edm::ParameterSet> VPSet;
0031   typedef std::vector<std::string> vstring;
0032   typedef boost::escaped_list_separator<char> elsc;
0033
0034   usesResource("DQMStore");
0035
0036   elsc commonEscapes("\\", " \t", "\'");
0037
0038   // Parse Normalization commands
0039   vstring normCmds = pset.getUntrackedParameter<vstring>("normalizationToIntegral", vstring());
0040   for (vstring::const_iterator normCmd = normCmds.begin(); normCmd != normCmds.end(); ++normCmd) {
0041     if (normCmd->empty())
0042       continue;
0043     boost::tokenizer<elsc> tokens(*normCmd, commonEscapes);
0044
0045     vector<string> args;
0046     for (boost::tokenizer<elsc>::const_iterator iToken = tokens.begin(); iToken != tokens.end(); ++iToken) {
0047       if (iToken->empty())
0048         continue;
0049       args.push_back(*iToken);
0050     }
0051
0052     if (args.empty() or args.size() > 2) {
0053       LogInfo("DQMRivetClient") << "Wrong input to normCmds\n";
0054       continue;
0055     }
0056
0057     NormOption opt;
0058     opt.name = args[0];
0059     opt.normHistName = args.size() == 2 ? args[1] : args[0];
0060
0061     normOptions_.push_back(opt);
0062   }
0063
0064   VPSet normSets = pset.getUntrackedParameter<VPSet>("normalizationToIntegralSets", VPSet());
0065   for (VPSet::const_iterator normSet = normSets.begin(); normSet != normSets.end(); ++normSet) {
0066     NormOption opt;
0067     opt.name = normSet->getUntrackedParameter<string>("name");
0068     opt.normHistName = normSet->getUntrackedParameter<string>("normalizedTo", opt.name);
0069
0070     normOptions_.push_back(opt);
0071   }
0072
0073   //normalize to lumi
0074   vstring lumiCmds = pset.getUntrackedParameter<vstring>("normalizationToLumi", vstring());
0075   for (vstring::const_iterator lumiCmd = lumiCmds.begin(); lumiCmd != lumiCmds.end(); ++lumiCmd) {
0076     if (lumiCmd->empty())
0077       continue;
0078     boost::tokenizer<elsc> tokens(*lumiCmd, commonEscapes);
0079
0080     vector<string> args;
0081     for (boost::tokenizer<elsc>::const_iterator iToken = tokens.begin(); iToken != tokens.end(); ++iToken) {
0082       if (iToken->empty())
0083         continue;
0084       args.push_back(*iToken);
0085     }
0086
0087     if (args.size() != 2) {
0088       LogInfo("DQMRivetClient") << "Wrong input to lumiCmds\n";
0089       continue;
0090     }
0091
0092     DQMRivetClient::LumiOption opt;
0093     opt.name = args[0];
0094     opt.normHistName = args[1];
0095     opt.xsection = pset.getUntrackedParameter<double>("xsection", -1.);
0096     //opt.xsection = atof(args[2].c_str());
0097
0098     //std::cout << opt.name << " " << opt.normHistName << " " << opt.xsection << std::endl;
0099     lumiOptions_.push_back(opt);
0100   }
0101
0102   //multiply by a number
0103   vstring scaleCmds = pset.getUntrackedParameter<vstring>("scaleBy", vstring());
0104   for (vstring::const_iterator scaleCmd = scaleCmds.begin(); scaleCmd != scaleCmds.end(); ++scaleCmd) {
0105     if (scaleCmd->empty())
0106       continue;
0107     boost::tokenizer<elsc> tokens(*scaleCmd, commonEscapes);
0108
0109     vector<string> args;
0110     for (boost::tokenizer<elsc>::const_iterator iToken = tokens.begin(); iToken != tokens.end(); ++iToken) {
0111       if (iToken->empty())
0112         continue;
0113       args.push_back(*iToken);
0114     }
0115
0116     if (args.empty() or args.size() > 2) {
0117       LogInfo("DQMRivetClient") << "Wrong input to normCmds\n";
0118       continue;
0119     }
0120
0121     ScaleFactorOption opt;
0122     opt.name = args[0];
0123     opt.scale = atof(args[1].c_str());
0124     scaleOptions_.push_back(opt);
0125   }
0126
0127   outputFileName_ = pset.getUntrackedParameter<string>("outputFileName", "");
0128   subDirs_ = pset.getUntrackedParameter<vstring>("subDirs");
0129 }
0130
0131 void DQMRivetClient::endRun(const edm::Run& r, const edm::EventSetup& c) {
0132   typedef vector<string> vstring;
0133
0134   // Update 2009-09-23
0135   // Migrated all code from endJob to this function
0136   // endJob is not necessarily called in the proper sequence
0137   // and does not necessarily book histograms produced in
0138   // that step.
0139   // It more robust to do the histogram manipulation in
0140   // this endRun function
0141
0142   theDQM = nullptr;
0143   theDQM = Service<DQMStore>().operator->();
0144
0145   if (!theDQM) {
0146     LogInfo("DQMRivetClient") << "Cannot create DQMStore instance\n";
0147     return;
0148   }
0149
0150   // Process wildcard in the sub-directory
0151   set<string> subDirSet;
0152
0153   for (vstring::const_iterator iSubDir = subDirs_.begin(); iSubDir != subDirs_.end(); ++iSubDir) {
0154     string subDir = *iSubDir;
0155
0156     if (subDir[subDir.size() - 1] == '/')
0157       subDir.erase(subDir.size() - 1);
0158
0159     subDirSet.insert(subDir);
0160   }
0161
0162   for (set<string>::const_iterator iSubDir = subDirSet.begin(); iSubDir != subDirSet.end(); ++iSubDir) {
0163     const string& dirName = *iSubDir;
0164     for (vector<NormOption>::const_iterator normOption = normOptions_.begin(); normOption != normOptions_.end();
0165          ++normOption) {
0166       normalizeToIntegral(dirName, normOption->name, normOption->normHistName);
0167     }
0168   }
0169
0170   for (set<string>::const_iterator iSubDir = subDirSet.begin(); iSubDir != subDirSet.end(); ++iSubDir) {
0171     const string& dirName = *iSubDir;
0172     for (vector<LumiOption>::const_iterator lumiOption = lumiOptions_.begin(); lumiOption != lumiOptions_.end();
0173          ++lumiOption) {
0174       normalizeToLumi(dirName, lumiOption->name, lumiOption->normHistName, lumiOption->xsection);
0175     }
0176   }
0177
0178   for (set<string>::const_iterator iSubDir = subDirSet.begin(); iSubDir != subDirSet.end(); ++iSubDir) {
0179     const string& dirName = *iSubDir;
0180     for (vector<ScaleFactorOption>::const_iterator scaleOption = scaleOptions_.begin();
0181          scaleOption != scaleOptions_.end();
0182          ++scaleOption) {
0183       scaleByFactor(dirName, scaleOption->name, scaleOption->scale);
0184     }
0185   }
0186
0187   if (!outputFileName_.empty())
0188     theDQM->save(outputFileName_);
0189 }
0190
0191 void DQMRivetClient::endJob() {
0192   // Update 2009-09-23
0193   // Migrated all code from here to endRun
0194
0195   LogTrace("DQMRivetClient") << "inside of ::endJob()" << endl;
0196 }
0197
0198 void DQMRivetClient::normalizeToIntegral(const std::string& startDir,
0199                                          const std::string& histName,
0200                                          const std::string& normHistName) {
0201   if (!theDQM->dirExists(startDir)) {
0202     LogInfo("DQMRivetClient") << "normalizeToEntries() : "
0203                               << "Cannot find sub-directory " << startDir << endl;
0204     return;
0205   }
0206
0207   theDQM->cd();
0208
0209   ME* element = theDQM->get(startDir + "/" + histName);
0210   ME* normME = theDQM->get(startDir + "/" + normHistName);
0211
0212   if (!element) {
0213     LogInfo("DQMRivetClient") << "normalizeToEntries() : "
0214                               << "No such element '" << histName << "' found\n";
0215     return;
0216   }
0217
0218   if (!normME) {
0219     LogInfo("DQMRivetClient") << "normalizeToEntries() : "
0220                               << "No such element '" << normHistName << "' found\n";
0221     return;
0222   }
0223
0224   TH1F* hist = element->getTH1F();
0225   if (!hist) {
0226     LogInfo("DQMRivetClient") << "normalizeToEntries() : "
0227                               << "Cannot create TH1F from ME\n";
0228     return;
0229   }
0230
0231   TH1F* normHist = normME->getTH1F();
0232   if (!normHist) {
0233     LogInfo("DQMRivetClient") << "normalizeToEntries() : "
0234                               << "Cannot create TH1F from ME\n";
0235     return;
0236   }
0237
0238   const double entries = normHist->Integral();
0239   if (entries != 0) {
0240     hist->Scale(1. / entries, "width");
0241   } else {
0242     LogInfo("DQMRivetClient") << "normalizeToEntries() : "
0243                               << "Zero entries in histogram\n";
0244   }
0245
0246   return;
0247 }
0248
0249 void DQMRivetClient::normalizeToLumi(const std::string& startDir,
0250                                      const std::string& histName,
0251                                      const std::string& normHistName,
0252                                      double xsection) {
0253   normalizeToIntegral(startDir, histName, normHistName);
0254   theDQM->cd();
0255   ME* element = theDQM->get(startDir + "/" + histName);
0256   TH1F* hist = element->getTH1F();
0257   if (!hist) {
0258     LogInfo("DQMRivetClient") << "normalizeToEntries() : "
0259                               << "Cannot create TH1F from ME\n";
0260     return;
0261   }
0262   hist->Scale(xsection);
0263   return;
0264 }
0265
0266 void DQMRivetClient::scaleByFactor(const std::string& startDir, const std::string& histName, double factor) {
0267   if (!theDQM->dirExists(startDir)) {
0268     LogInfo("DQMRivetClient") << "normalizeToEntries() : "
0269                               << "Cannot find sub-directory " << startDir << endl;
0270     return;
0271   }
0272
0273   theDQM->cd();
0274
0275   ME* element = theDQM->get(startDir + "/" + histName);
0276
0277   if (!element) {
0278     LogInfo("DQMRivetClient") << "normalizeToEntries() : "
0279                               << "No such element '" << histName << "' found\n";
0280     return;
0281   }
0282
0283   TH1F* hist = element->getTH1F();
0284   if (!hist) {
0285     LogInfo("DQMRivetClient") << "normalizeToEntries() : "
0286                               << "Cannot create TH1F from ME\n";
0287     return;
0288   }
0289   hist->Scale(factor);
0290 }
0291
0292 #include "FWCore/Framework/interface/MakerMacros.h"
0293
0294 DEFINE_FWK_MODULE(DQMRivetClient);