Back to home page

Project CMSSW displayed by LXR

 
 

    


File indexing completed on 2021-02-14 12:51:15

0001 #ifndef STANDALONE
0002 #include <CondFormats/JetMETObjects/interface/JetResolutionObject.h>
0003 #include <CondFormats/JetMETObjects/interface/Utilities.h>
0004 #include <FWCore/Utilities/interface/EDMException.h>
0005 
0006 #else
0007 #include "JetResolutionObject.h"
0008 #include "Utilities.h"
0009 #include <exception>
0010 
0011 namespace edm {
0012   namespace errors {
0013     enum ErrorCode { NotFound = 8026, ConfigFileReadError = 7002, UnimplementedFeature = 8011, FileReadError = 8021 };
0014   };
0015 };  // namespace edm
0016 #endif
0017 
0018 #include <cmath>
0019 #include <iostream>
0020 #include <fstream>
0021 #include <iomanip>
0022 #include <algorithm>
0023 
0024 namespace JME {
0025 
0026   std::string getDefinitionLine(const std::string& line) {
0027     size_t first = line.find('{');
0028     size_t last = line.find('}');
0029 
0030     if (first != std::string::npos && last != std::string::npos && first < last)
0031       return std::string(line, first + 1, last - first - 1);
0032 
0033     return "";
0034   }
0035 
0036   void throwException(uint32_t code, const std::string& message) {
0037 #ifndef STANDALONE
0038     throw edm::Exception(static_cast<edm::errors::ErrorCodes>(code), message);
0039 #else
0040     std::stringstream error;
0041     error << message << " Error code: " << code;
0042     throw std::runtime_error(error.str());
0043 
0044 #endif
0045   }
0046 
0047   const bimap<Binning, std::string> JetParameters::binning_to_string = {{Binning::JetPt, "JetPt"},
0048                                                                         {Binning::JetEta, "JetEta"},
0049                                                                         {Binning::JetAbsEta, "JetAbsEta"},
0050                                                                         {Binning::JetE, "JetE"},
0051                                                                         {Binning::JetArea, "JetArea"},
0052                                                                         {Binning::Mu, "Mu"},
0053                                                                         {Binning::Rho, "Rho"},
0054                                                                         {Binning::NPV, "NPV"}};
0055 
0056   JetParameters::JetParameters(JetParameters&& rhs) { m_values = std::move(rhs.m_values); }
0057 
0058   JetParameters::JetParameters(std::initializer_list<typename value_type::value_type> init) {
0059     for (auto& i : init) {
0060       set(i.first, i.second);
0061     }
0062   }
0063 
0064   JetParameters& JetParameters::setJetPt(float pt) {
0065     m_values[Binning::JetPt] = pt;
0066     return *this;
0067   }
0068 
0069   JetParameters& JetParameters::setJetEta(float eta) {
0070     m_values[Binning::JetEta] = eta;
0071     m_values[Binning::JetAbsEta] = fabs(eta);
0072     return *this;
0073   }
0074 
0075   JetParameters& JetParameters::setJetE(float e) {
0076     m_values[Binning::JetE] = e;
0077     return *this;
0078   }
0079 
0080   JetParameters& JetParameters::setJetArea(float area) {
0081     m_values[Binning::JetArea] = area;
0082     return *this;
0083   }
0084 
0085   JetParameters& JetParameters::setMu(float mu) {
0086     m_values[Binning::Mu] = mu;
0087     return *this;
0088   }
0089 
0090   JetParameters& JetParameters::setNPV(float npv) {
0091     m_values[Binning::NPV] = npv;
0092     return *this;
0093   }
0094 
0095   JetParameters& JetParameters::setRho(float rho) {
0096     m_values[Binning::Rho] = rho;
0097     return *this;
0098   }
0099 
0100   JetParameters& JetParameters::set(const Binning& bin, float value) {
0101     m_values.emplace(bin, value);
0102 
0103     // Special case for eta
0104     if (bin == Binning::JetEta) {
0105       m_values.emplace(Binning::JetAbsEta, fabs(value));
0106     }
0107 
0108     return *this;
0109   }
0110 
0111   JetParameters& JetParameters::set(const typename value_type::value_type& value) {
0112     set(value.first, value.second);
0113     return *this;
0114   }
0115 
0116   std::vector<float> JetParameters::createVector(const std::vector<Binning>& binning) const {
0117     std::vector<float> values;
0118     for (const auto& bin : binning) {
0119       const auto& it = m_values.find(bin);
0120       if (it == m_values.cend()) {
0121         throwException(edm::errors::NotFound,
0122                        "JER parametrisation depends on '" + JetParameters::binning_to_string.left.at(bin) +
0123                            "' but no value for this parameter has been specified. Please call the appropriate 'set' "
0124                            "function of the JME::JetParameters object");
0125       }
0126 
0127       values.push_back(it->second);
0128     }
0129 
0130     return values;
0131   }
0132 
0133   JetResolutionObject::Definition::Definition(const std::string& definition) {
0134     std::vector<std::string> tokens = getTokens(definition);
0135 
0136     // We need at least 3 tokens
0137     if (tokens.size() < 3) {
0138       throwException(edm::errors::ConfigFileReadError,
0139                      "Definition line needs at least three tokens. Please check file format.");
0140     }
0141 
0142     size_t n_bins = std::stoul(tokens[0]);
0143 
0144     if (tokens.size() < (n_bins + 2)) {
0145       throwException(edm::errors::ConfigFileReadError, "Invalid file format. Please check.");
0146     }
0147 
0148     for (size_t i = 0; i < n_bins; i++) {
0149       m_bins_name.push_back(tokens[i + 1]);
0150     }
0151 
0152     size_t n_variables = std::stoul(tokens[n_bins + 1]);
0153 
0154     if (tokens.size() < (1 + n_bins + 1 + n_variables + 1)) {
0155       throwException(edm::errors::ConfigFileReadError, "Invalid file format. Please check.");
0156     }
0157 
0158     for (size_t i = 0; i < n_variables; i++) {
0159       m_variables_name.push_back(tokens[n_bins + 2 + i]);
0160     }
0161 
0162     m_formula_str = tokens[n_bins + n_variables + 2];
0163 
0164     std::string formula_str_lower = m_formula_str;
0165     std::transform(formula_str_lower.begin(), formula_str_lower.end(), formula_str_lower.begin(), ::tolower);
0166 
0167     if (formula_str_lower == "none") {
0168       m_formula_str = "";
0169 
0170       if ((tokens.size() > n_bins + n_variables + 3) && (std::atoi(tokens[n_bins + n_variables + 3].c_str()))) {
0171         size_t n_parameters = std::stoul(tokens[n_bins + n_variables + 3]);
0172 
0173         if (tokens.size() < (1 + n_bins + 1 + n_variables + 1 + 1 + n_parameters)) {
0174           throwException(edm::errors::ConfigFileReadError, "Invalid file format. Please check.");
0175         }
0176 
0177         for (size_t i = 0; i < n_parameters; i++) {
0178           m_formula_str += tokens[n_bins + n_variables + 4 + i] + " ";
0179         }
0180       }
0181     }
0182 
0183     init();
0184   }
0185 
0186   void JetResolutionObject::Definition::init() {
0187     if (!m_formula_str.empty()) {
0188       if (m_formula_str.find(' ') == std::string::npos)
0189 #ifndef STANDALONE
0190         m_formula = std::make_shared<reco::FormulaEvaluator>(m_formula_str);
0191 #else
0192         m_formula = std::make_shared<TFormula>("jet_resolution_formula", m_formula_str.c_str());
0193 #endif
0194       else
0195         m_parameters_name = getTokens(m_formula_str);
0196     }
0197     for (const auto& bin : m_bins_name) {
0198       const auto& b = JetParameters::binning_to_string.right.find(bin);
0199       if (b == JetParameters::binning_to_string.right.cend()) {
0200         throwException(edm::errors::UnimplementedFeature, "Bin name not supported: '" + bin + "'");
0201       }
0202       m_bins.push_back(b->second);
0203     }
0204 
0205     for (const auto& v : m_variables_name) {
0206       const auto& var = JetParameters::binning_to_string.right.find(v);
0207       if (var == JetParameters::binning_to_string.right.cend()) {
0208         throwException(edm::errors::UnimplementedFeature, "Variable name not supported: '" + v + "'");
0209       }
0210       m_variables.push_back(var->second);
0211     }
0212   }
0213 
0214   JetResolutionObject::Record::Record(const std::string& line, const Definition& def) {
0215     std::vector<std::string> tokens = getTokens(line);
0216 
0217     if (tokens.size() < (def.nBins() * 2 + def.nVariables() * 2 + 1)) {
0218       throwException(edm::errors::ConfigFileReadError, "Invalid record. Please check file format. Record: " + line);
0219     }
0220 
0221     size_t pos = 0;
0222 
0223     for (size_t i = 0; i < def.nBins(); i++) {
0224       Range r(std::stof(tokens[pos]), std::stof(tokens[pos + 1]));
0225       pos += 2;
0226       m_bins_range.push_back(r);
0227     }
0228 
0229     size_t n_parameters = std::stoul(tokens[pos++]);
0230 
0231     if (tokens.size() < (def.nBins() * 2 + def.nVariables() * 2 + 1 + (n_parameters - def.nVariables() * 2))) {
0232       throwException(edm::errors::ConfigFileReadError, "Invalid record. Please check file format. Record: " + line);
0233     }
0234 
0235     for (size_t i = 0; i < def.nVariables(); i++) {
0236       Range r(std::stof(tokens[pos]), std::stof(tokens[pos + 1]));
0237       pos += 2;
0238       m_variables_range.push_back(r);
0239       n_parameters -= 2;
0240     }
0241 
0242     for (size_t i = 0; i < n_parameters; i++) {
0243       m_parameters_values.push_back(std::stof(tokens[pos++]));
0244     }
0245   }
0246 
0247   JetResolutionObject::JetResolutionObject(const std::string& filename) {
0248     // Parse file
0249     std::ifstream f(filename);
0250 
0251     if (!f.good()) {
0252       throwException(edm::errors::FileReadError, "Can't read input file '" + filename + "'");
0253     }
0254 
0255     for (std::string line; std::getline(f, line);) {
0256       if ((line.empty()) || (line[0] == '#'))
0257         continue;
0258 
0259       std::string definition = getDefinitionLine(line);
0260 
0261       if (!definition.empty()) {
0262         m_definition = Definition(definition);
0263       } else {
0264         m_records.push_back(Record(line, m_definition));
0265       }
0266     }
0267 
0268     m_valid = true;
0269   }
0270 
0271   JetResolutionObject::JetResolutionObject(const JetResolutionObject& object) {
0272     m_definition = object.m_definition;
0273     m_records = object.m_records;
0274     m_valid = object.m_valid;
0275 
0276     m_definition.init();
0277   }
0278 
0279   JetResolutionObject::JetResolutionObject() {
0280     // Empty
0281   }
0282 
0283   void JetResolutionObject::dump() const {
0284     std::cout << "Definition: " << std::endl;
0285     std::cout << "    Number of binning variables: " << m_definition.nBins() << std::endl;
0286     std::cout << "        ";
0287     for (const auto& bin : m_definition.getBinsName()) {
0288       std::cout << bin << ", ";
0289     }
0290     std::cout << std::endl;
0291     std::cout << "    Number of variables: " << m_definition.nVariables() << std::endl;
0292     std::cout << "        ";
0293     for (const auto& bin : m_definition.getVariablesName()) {
0294       std::cout << bin << ", ";
0295     }
0296     std::cout << std::endl;
0297     std::cout << "    Formula: " << m_definition.getFormulaString() << std::endl;
0298 
0299     std::cout << std::endl << "Bin contents" << std::endl;
0300 
0301     for (const auto& record : m_records) {
0302       std::cout << "    Bins" << std::endl;
0303       size_t index = 0;
0304       for (const auto& bin : record.getBinsRange()) {
0305         std::cout << "        " << m_definition.getBinName(index) << " [" << bin.min << " - " << bin.max << "]"
0306                   << std::endl;
0307         index++;
0308       }
0309 
0310       std::cout << "    Variables" << std::endl;
0311       index = 0;
0312       for (const auto& r : record.getVariablesRange()) {
0313         std::cout << "        " << m_definition.getVariableName(index) << " [" << r.min << " - " << r.max << "] "
0314                   << std::endl;
0315         index++;
0316       }
0317 
0318       std::cout << "    Parameters" << std::endl;
0319       index = 0;
0320       for (const auto& par : record.getParametersValues()) {
0321         std::cout << "        Parameter #" << index << " = " << par << std::endl;
0322         index++;
0323       }
0324     }
0325   }
0326 
0327   void JetResolutionObject::saveToFile(const std::string& file) const {
0328     std::ofstream fout(file);
0329     fout.setf(std::ios::right);
0330 
0331     // Definition
0332     fout << "{" << m_definition.nBins();
0333 
0334     for (auto& bin : m_definition.getBinsName())
0335       fout << "    " << bin;
0336 
0337     fout << "    " << m_definition.nVariables();
0338 
0339     for (auto& var : m_definition.getVariablesName())
0340       fout << "    " << var;
0341 
0342     fout << "    " << (m_definition.getFormulaString().empty() ? "None" : m_definition.getFormulaString())
0343          << "    Resolution}" << std::endl;
0344 
0345     // Records
0346     for (auto& record : m_records) {
0347       for (auto& r : record.getBinsRange()) {
0348         fout << std::left << std::setw(15) << r.min << std::setw(15) << r.max << std::setw(15);
0349       }
0350       fout << (record.nVariables() * 2 + record.nParameters()) << std::setw(15);
0351 
0352       for (auto& r : record.getVariablesRange()) {
0353         fout << r.min << std::setw(15) << r.max << std::setw(15);
0354       }
0355 
0356       for (auto& p : record.getParametersValues()) {
0357         fout << p << std::setw(15);
0358       }
0359 
0360       fout << std::endl << std::setw(0);
0361     }
0362   }
0363 
0364   const JetResolutionObject::Record* JetResolutionObject::getRecord(const JetParameters& bins_parameters) const {
0365     // Find record for bins
0366     if (!m_valid)
0367       return nullptr;
0368 
0369     // Create vector of bins value. Throw if some values are missing
0370     std::vector<float> bins = bins_parameters.createVector(m_definition.getBins());
0371 
0372     // Iterate over all records, and find the one for which all bins are valid
0373     const Record* good_record = nullptr;
0374     for (const auto& record : m_records) {
0375       // Iterate over bins
0376       size_t valid_bins = 0;
0377       size_t current_bin = 0;
0378       for (const auto& bin : record.getBinsRange()) {
0379         if (bin.is_inside(bins[current_bin]))
0380           valid_bins++;
0381 
0382         current_bin++;
0383       }
0384 
0385       if (valid_bins == m_definition.nBins()) {
0386         good_record = &record;
0387         break;
0388       }
0389     }
0390 
0391     return good_record;
0392   }
0393 
0394   float JetResolutionObject::evaluateFormula(const JetResolutionObject::Record& record,
0395                                              const JetParameters& variables_parameters) const {
0396     if (!m_valid)
0397       return 1;
0398 
0399 #ifndef STANDALONE
0400     const auto* formula = m_definition.getFormula();
0401 #else
0402     // Set parameters
0403     auto const* pFormula = m_definition.getFormula();
0404     if (!pFormula)
0405       return 1;
0406     auto formula = *pFormula;
0407 #endif
0408     // Create vector of variables value. Throw if some values are missing
0409     std::vector<float> variables = variables_parameters.createVector(m_definition.getVariables());
0410 
0411     double variables_[4] = {0};
0412     for (size_t index = 0; index < m_definition.nVariables(); index++) {
0413       variables_[index] =
0414           clip(variables[index], record.getVariablesRange()[index].min, record.getVariablesRange()[index].max);
0415     }
0416     const std::vector<float>& parameters = record.getParametersValues();
0417 
0418 #ifndef STANDALONE
0419     //ArrayAdaptor only takes doubles
0420     std::vector<double> parametersD(parameters.begin(), parameters.end());
0421     return formula->evaluate(reco::formula::ArrayAdaptor(variables_, m_definition.nVariables()),
0422                              reco::formula::ArrayAdaptor(parametersD.data(), parametersD.size()));
0423 #else
0424     for (size_t index = 0; index < parameters.size(); index++) {
0425       formula.SetParameter(index, parameters[index]);
0426     }
0427 
0428     return formula.EvalPar(variables_);
0429 #endif
0430   }
0431 }  // namespace JME
0432 
0433 #ifndef STANDALONE
0434 #include "FWCore/Utilities/interface/typelookup.h"
0435 TYPELOOKUP_DATA_REG(JME::JetResolutionObject);
0436 #endif